Skip to content
Snippets Groups Projects
BatchCoverage.pl 1.38 KiB
#!/usr/bin/perl

use warnings;
use strict;
use File::Basename;

sub parseConfig($$);

MAIN:
{
    my $file_config = shift;
    my $path_bam = shift;
    my $path_gbed = shift;
    my $tag = shift;

    my @bams = ();
    parseConfig(\@bams, $file_config);
    $path_bam =~ s/\/$//g;
    $path_gbed =~ s/\/$//g;

    foreach my $qry (@bams)
    {
        my $file_name = fileparse($qry->[0], $tag . ".bam");
        my $file_bam = $path_bam . "/" . $qry->[0];
        my $file_gbed = $path_gbed . "/" . $file_name . ".gbed.gz";
        my $command = "samtools view -buh -q 255 " . $file_bam . " | " .
        "bedtools genomecov" .
        " -ibam -" .
        " -scale " . $qry->[2] . 
        " -bg -split |" .
        " sort -k1,1 -k2,2n -k3,3n |" .
        " bgzip > " .
        $file_gbed;

        print $file_name,"\n";
        system($command);
        #print $command,"\n";

        $command = "tabix --zero-based --sequence 1 --begin 2 --end 3 " .
        $file_gbed;
        #print $command,"\n";
        system($command);

    }

}

sub parseConfig($$)
{
    my $bams_ref = $_[0];
    my $file_config = $_[1];

    open(my $fh, "<", $file_config) or die $!;
    while (<$fh>)
    {
        chomp($_);
        next if ($_ =~ m/^#/);
        my ($file_bam, $group, $rpm) = split("\t", $_, 3);
        $rpm = 1000000 / $rpm;

        push(@{$bams_ref}, [$file_bam, $group, $rpm]);
    }
    close($fh);
}