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.1';
  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, s=>100);
  40   getopts('D:m:rs:', \%opts); # -s for scaling factor in score calculation
  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         # calculate indel score
  60         my $score = $t[5];
  61         $score += $opts{s} * $t[10] if ($t[8] ne '*');
  62         $score += $opts{s} * $t[11] if ($t[9] ne '*');
  63         @$curr = ($t[0], $t[1], $score, $_);
  64         my $do_swap = 1;
  65         if (defined $last->[0]) {
  66           if ($curr->[0] eq $last->[0] && $last->[1] + $opts{m} > $curr->[1]) {
  67                 $do_swap = 0 if ($last->[2] > $curr->[2]);
  68           } else { # then print
  69                 print $last->[3];
  70           }
  71         }
  72         if ($do_swap) {
  73           my $tmp = $curr; $curr = $last; $last = $tmp;
  74         }
  75   }
  76   print $last->[3] if (defined $last->[0]);
  77 }
  78
  79 sub usage {
  80   die(qq/
  81 Usage:   samtools.pl <command> [<arguments>]\n
  82 Command: indelFilter   filter indels generated by `pileup -c'
  83          showALEN      print alignment length (ALEN) following CIGAR
  84 \n/);
  85 }