misc/samtools.pl

   1 #!/usr/bin/perl -w
   2
   3 # Author: lh3
   4
   5 use strict;
   6 use warnings;
   7 use Getopt::Std;
   8
   9 my $version = '0.1.0';
  10 &usage if (@ARGV < 1);
  11
  12 my $command = shift(@ARGV);
  13 my %func = (indelFilter=>\&indelFilter, showALEN=>\&showALEN);
  14
  15 die("Unknown command \"$command\".\n") if (!defined($func{$command}));
  16 &{$func{$command}};
  17 exit(0);
  18
  19 #
  20 # showALEN
  21 #
  22
  23 sub showALEN {
  24   die(qq/Usage: samtools.pl showALEN <in.sam>\n/) if (@ARGV == 0 && -t STDIN);
  25   while (<>) {
  26         my @t = split;
  27         my $l = 0;
  28         $_ = $t[5];
  29         s/(\d+)[SMI]/$l+=$1/eg;
  30         print join("\t", @t[0..5]), "\t$l\t", join("\t", @t[6..$#t]), "\n";
  31   }
  32 }
  33
  34 #
  35 # indelFilter
  36 #
  37
  38 sub indelFilter {
  39   my %opts = (D=>100, m=>10, r=>undef);
  40   getopts('D:m:r', \%opts);
  41
  42   die(qq/
  43 Usage:   samtools.pl indelFilter [options] <in.indel>\n
  44 Options: -D INT    maximum read depth [$opts{D}]
  45          -m INT    minimum distance between two adjacent indels [$opts{m}]
  46 \n/) if (@ARGV == 0 && -t STDIN);
  47
  48   my (@arr1, @arr2);
  49   my ($curr, $last) = (\@arr1, \@arr2);
  50   my $is_ref = defined($opts{r})? 1 : 0;
  51   while (<>) {
  52         my @t = split;
  53         next if ($t[2] ne '*');
  54         if (!$is_ref) {
  55           next if ($t[3] eq '*/*');
  56           next if ($t[5] == 0);
  57         }
  58         next if ($t[7] > $opts{D});
  59         @$curr = ($t[0], $t[1], $t[5], $_);
  60         my $do_swap = 1;
  61         if (defined $last->[0]) {
  62           if ($curr->[0] eq $last->[0] && $last->[1] + $opts{m} > $curr->[1]) {
  63                 $do_swap = 0 if ($last->[2] > $curr->[2]);
  64           } else { # then print
  65                 print $last->[3];
  66           }
  67         }
  68         if ($do_swap) {
  69           my $tmp = $curr; $curr = $last; $last = $tmp;
  70         }
  71   }
  72   print $last->[3] if (defined $last->[0]);
  73 }
  74
  75 sub usage {
  76   die(qq/
  77 Usage:   samtools.pl <command> [<arguments>]\n
  78 Command: indelFilter   filter indels generated by `pileup -c'
  79          showALEN      print alignment length (ALEN) following CIGAR
  80 \n/);
  81 }