perl/etc/samrstat.pl

   1 #!/bin/env perl
   2 use strict;
   3 use warnings;
   4
   5 die "Usage: $0 <filled sam> <output>\n" if @ARGV <1;
   6 my ($in,$out)=@ARGV;
   7 $out=$in.'o' unless $out;
   8 warn "From [$in] to [$out]\n";
   9
  10 sub openfile($) {
  11     my ($filename)=@_;
  12     my $infile;
  13     if ($filename=~/.bz2$/) {
  14             open( $infile,"-|","bzip2 -dc $filename") or die "Error opening $filename: $!\n";
  15     } elsif ($filename=~/.gz$/) {
  16         open( $infile,"-|","gzip -dc $filename") or die "Error opening $filename: $!\n";
  17     } else {open( $infile,"<",$filename) or die "Error opening $filename: $!\n";}
  18     return $infile;
  19 }
  20
  21 my (%Count,$CountRF,$Reads);
  22 ($CountRF,$Reads,$Count{0},$Count{1},$Count{-1})=(0,0,0,0,0);
  23 my $fh=openfile($in);
  24 while (<$fh>) {
  25         next if /^@\w\w\t\w\w:/;
  26         chomp;
  27         my @read=split /\t/, $_;
  28         my $Alternativehits='';
  29         for ( @read ) {
  30         if (/^XA:Z:([\w,+-;]+)$/) { #XA:Z:chrX,+1144912,100M,0;
  31                         $Alternativehits = $1;
  32                         last;
  33                 }
  34         }
  35         my @Alt=split ';',$Alternativehits;
  36         ++$Reads;
  37         my $flag=0;
  38         for (@Alt) {
  39                 my ($chr,$pos,$CIGAR,$NM)=split /,/;
  40                 if (substr($pos,-2) eq '01') {
  41                         if ($pos>0) {
  42                                 $flag=1;
  43                         } else {$flag=-1;}
  44                 }
  45         }
  46         ++$Count{$flag};
  47         ++$CountRF if $flag;
  48 }
  49 close $fh;
  50 open O,'>',"$out" or die "[x]Error opening $out: $!\n";
  51 print O $Count{1}+$Count{-1}," ,RF:$CountRF\t+$Count{1},-$Count{-1},z$Count{0}\t$Reads\t$in\n";
  52 close O;
  53 print $Count{1}+$Count{-1}," ,RF:$CountRF\t+$Count{1},-$Count{-1},z$Count{0}\t$Reads\n";
  54 __END__
  55 find . -name '*.r'|xargs -n1 ./samrstat.pl
  56 find . -name '*.ro'|xargs -n1 cat|sort -n