htswanalysis/scripts/count_reads.pm

   1 #!/usr/bin/perl -w
   2
   3 use strict;
   4 use warnings;
   5
   6 my $align_file = shift;
   7 my $all_gz_file = shift;
   8
   9 if(!defined($all_gz_file)) { print STDERR "Usage: count_reads.pm [align file] [all file.gz]\n"; exit(1); }
  10
  11 my $pf_count = 0;
  12 my $adapter_count = 0;
  13 my $align_count = 0;
  14
  15 open(ALIGN,$align_file);
  16 while(<ALIGN>) {
  17   chomp;
  18   my(@a) = split(/ /,$_);
  19   $pf_count++;
  20
  21   if(scalar(@a) > 3) {
  22     if($_ =~ /contam/) { $adapter_count++; } else { $align_count++; }
  23   }
  24 }
  25 close(ALIGN);
  26
  27 my $all_count = `cat $all_gz_file | gzip -d | wc -l`;
  28 chomp $all_count; $all_count =~ s/\s//;
  29
  30 print "#N_reads N_pf N_adapter N_align\n";
  31 print "$all_count\t$pf_count\t$adapter_count\t$align_count\n";