perl/etc/radseq/samkeepr1.pl

   1 #!/usr/bin/env perl
   2 =pod
   3 Author: Hu Xuesong @ BIOPIC <galaxy001@gmail.com>
   4 Version: 1.0.0 @ 20120330
   5 =cut
   6 use strict;
   7 use warnings;
   8 #use Data::Dump qw(ddx);
   9
  10 die "Usage: $0 <sam in> <out prefix>\n" if @ARGV<2;
  11 my $in=shift;
  12 my $outp=shift;
  13
  14 sub openfile($) {
  15     my ($filename)=@_;
  16     my $infile;
  17     if ($filename=~/.bz2$/) {
  18             open( $infile,"-|","bzip2 -dc $filename") or die "Error opening $filename: $!\n";
  19     } elsif ($filename=~/.gz$/) {
  20         open( $infile,"-|","gzip -dc $filename") or die "Error opening $filename: $!\n";
  21     } else {open( $infile,"<",$filename) or die "Error opening $filename: $!\n";}
  22     return $infile;
  23 }
  24
  25 my $samin = openfile($in);
  26 open O,'|-',"gzip -9c >$outp.sam.gz" or die "Error opening $outp.sam.gz with gzip: $!\n";
  27 open L,'>',$outp.'.log' or die "Error opening $outp.log: $!\n";
  28 select(L);
  29 $|=1;
  30 print L "From [$in] to [$outp.sam.gz]\n";
  31 my ($Total,$Out,$notOut)=(0,0,0);
  32
  33 while (<$samin>) {
  34         if (/^@\w\w\t\w\w:/) {
  35                 print O $_;
  36                 next;
  37         }
  38         my @read1=split /\t/;
  39         if (($read1[1] & 64) == 64) {
  40                 print O $_;
  41                 ++$Out;
  42         } elsif (($read1[1] & 128) == 128) {
  43                 ++$notOut;
  44         }
  45         ++$Total;
  46 }
  47
  48 close $samin;
  49 close O;
  50 print L "Read_1: $Out , ",$Out/$Total,"\nRead_2: $notOut , ",$notOut/$Total,"\nTotal: $Total\nRemain: ",$Total-$Out-$notOut,"\n";
  51 close L;