perl/etc/WoodyMiaoLin/PbePopGen/count_FLAG_CIGAR.pl

   1 #!/usr/bin/perl
   2 use strict;
   3 use warnings;
   4 use threads;
   5
   6 open I, "-|", "ls -1 ../c.sort_rmdup_bam/*.bam";
   7 open O, ">", "count_FLAG_CIGAR_result.txt";
   8
   9 my %file;
  10 while (<I>) {
  11         chomp;
  12         $file{$_} = openfile($_);
  13 }
  14 close I;
  15 print STDERR "Read file list complete!\n";
  16
  17 my %thread;
  18 foreach (sort keys %file) {
  19         $thread{$_} = threads->new(\&countsam, $file{$_});
  20 }
  21
  22 my %count;
  23 foreach (sort keys %thread) {
  24         $count{$_} = $thread{$_}->join;
  25         print STDERR "Count $_ complete!\n";
  26 }
  27
  28 print O "#File\t#FLAG\t#CIGAR\n";
  29 foreach (sort keys %count) {
  30         print O $_, "\t", (join "\t", @{$count{$_}}), "\n";
  31 }
  32 close O;
  33 print STDERR "Done!\n";
  34
  35 sub openfile {
  36         my $filename = shift;
  37         my $infile;
  38         open $infile, "-|", "samtools-0.1.7 view $filename" or die "Error opening $filename: $!\n";
  39         return $infile;
  40 }
  41
  42 sub countsam {
  43         my $in = shift;
  44         my %FLAG;
  45         my %CIGAR;
  46         while (<$in>) {
  47                 my @a = split /\t/;
  48                 ++$FLAG{$a[1]};
  49                 my @b = $a[5] =~ /(\d+\D)/g;
  50                 foreach (@b) {
  51                         /(\d+)(\D)/;
  52                         $CIGAR{$2} += $1;
  53                 }
  54         }
  55         close $in;
  56         my $flag;
  57         my $cigar;
  58         $flag .= "#$_=$FLAG{$_}," foreach sort keys %FLAG;
  59         $cigar .= "$CIGAR{$_}$_," foreach sort keys %CIGAR;
  60         return [$flag, $cigar];
  61 }