bamqc_dna.pl

#!/usr/bin/perl
#uploadqc.pl

open OUT, ">sequence.stats.txt" or die $!;
print OUT join("\t",'Sample','total.raw','total.trimmed','ontarget','maprate','propair','Percent.ontarget','Mean.MapQ',
	       'Percent.QualReads','Median.Mismatch','Indel.Rate','Error.Rate','Percent.Dups','Library.Size',
	       'Median.Insert','Mean.Insert','Std.insert',
	       'all.avg.depth','all.median.depth','all.perc.100x','all.perc.200x','all.perc.500x',
	       'dedup.avg.depth','dedup.median.depth','dedup.perc.100x','dedup.perc.200x','dedup.perc.500x'),"\n";

my @statfiles = @ARGV;

foreach $sfile (@statfiles) {
  $sfile =~ m/(\S+)\.genomecov.txt/;
  my $prefix = $1;
  my %hash;
  open FLAG, "<$prefix\.trimreport.txt" or die $!;
  while (my $line = <FLAG>) {
    chomp($line);
    my ($file,$raw,$trim) = split(/\t/,$line);
    $hash{rawct} += $raw;
    $hash{trimct} += $trim;
  }
  open FLAG, "<$prefix\.flagstat.txt" or die $!;
  while (my $line = <FLAG>) {
    chomp($line);
    if ($line =~ m/(\d+) \+ \d+ in total/) {
      $hash{total} = $1;
      $hash{total} = $hash{trimct} if ($hash{trimct} && $hash{trimct} > $hash{total});
    }elsif ($line =~ m/(\d+) \+ \d+ read1/) {
      $hash{pairs} = $1;
    }elsif ($line =~ m/(\d+) \+ \d+ mapped/) {
      $hash{maprate} = 100*sprintf("%.4f",$1/$hash{total});
    }elsif ($line =~ m/(\d+) \+ \d+ properly paired/) {
      $hash{propair} = 100*sprintf("%.4f",$1/$hash{total});
    }elsif ($line =~ m/(\d+) \+ \d+ properly paired/) {
      $hash{propair} = 100*sprintf("%.4f",$1/$hash{total});
	}
  }
  # open FLAG, "<$prefix\.ontarget.flagstat.txt" or die $!;
  # while (my $line = <FLAG>) {
  #   chomp($line);
  #   if ($line =~ m/(\d+) \+ \d+ in total/) {
  #     $hash{ontarget} = $1;
  #   }elsif ($line =~ m/(\d+) \+ \d+ read1/) {
  #     $hash{ontarget_pairs} = $1;
  #   }elsif ($line =~ m/(\d+) \+ \d+ properly paired/) {
  #     $hash{ontarget_propair} = 100*sprintf("%.4f",$1/$hash{total});
  #   }elsif ($line =~ m/(\d+) \+ \d+ properly paired/) {
  #     $hash{ontarget_propair} = 100*sprintf("%.4f",$1/$hash{total});
  #   }
  # }
  open ASM, "<$prefix\.alignmentsummarymetrics.txt" or die $!;
  while (my $line = <ASM>) {
    chomp($line);
    if ($line =~ m/## METRICS/) {
      $header = <ASM>;
      chomp($header);
      my @stats = split(/\t/,$header);
      while (my $nums = <ASM>) {
	chomp($nums);
	next unless ($nums =~ m/^PAIR/);
	my @row = split(/\t/,$nums);
	my %info;
	foreach my $i (0..$#stats) {
	  $info{$stats[$i]} = $row[$i];
	}
	$hash{ontarget} = $info{PF_READS_ALIGNED};
	$hash{qualreads} = $info{PF_HQ_ALIGNED_READS};
	$hash{median_mismatch} = $info{PF_MISMATCH_RATE};