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.2.3';
  10 &usage if (@ARGV < 1);
  11
  12 my $command = shift(@ARGV);
  13 my %func = (snpFilter=>\&snpFilter, indelFilter=>\&indelFilter, showALEN=>\&showALEN, pileup2fq=>\&pileup2fq);
  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 #
  80 # snpFilter
  81 #
  82
  83 sub snpFilter {
  84   my %opts = (f=>'', Q=>40, d=>3, w=>10, D=>0, N=>2, W=>10, q=>20, s=>50);
  85   getopts('f:s:w:q:Q:d:D:W:N:', \%opts);
  86   die(qq{
  87 Usage:   samtools.pl snpFilter [options] <cns2snp.snp>
  88
  89 Options: -d INT        minimum depth to call a SNP [$opts{d}]
  90          -D INT        maximum depth, 0 to ignore [$opts{D}]
  91          -Q INT        required max mapping quality of the reads covering the SNP [$opts{Q}]
  92          -q INT        minimum SNP quality [$opts{q}]
  93
  94          -f FILE       filtered samtools indels [null]
  95          -s INT        minimum samtols indel score [$opts{s}]
  96          -w INT        SNP within INT bp around an indel to be filtered [$opts{w}]
  97
  98          -W INT        window size for filtering dense SNPs [$opts{W}]
  99          -N INT        maximum number of SNPs in a window [$opts{N}]
 100 \n}) unless (@ARGV);
 101   my (%hash, $fh);
 102   my $skip = $opts{w};
 103   $opts{D} = 100000000 if ($opts{D} == 0);
 104   if ($opts{f}) { # filtered samtools indel
 105         my $n = 0;
 106         open($fh, $opts{f}) || die;
 107         while (<$fh>) {
 108           my @t = split;
 109           next if ($t[2] ne '*' || $t[3] eq '*/*' || $t[5] < $opts{s});
 110           for (my $x = $t[1] - $skip + 1; $x < $t[1] + $skip; ++$x) {
 111                 $hash{$t[0],$x} = 1;
 112           }
 113         }
 114         close($fh);
 115   }
 116   my (@last, $last_chr);
 117   $last_chr = '';
 118   while (<>) {
 119         my @t = split;
 120         next if ($t[2] eq '*' || $hash{$t[0],$t[1]});
 121         next if ($t[2] eq $t[3]);
 122         my $is_good = ($t[7] >= $opts{d} && $t[7] <= $opts{D} && $t[6] >= $opts{Q} && $t[5] >= $opts{q})? 1 : 0;
 123         next unless ($is_good); # drop
 124         if ($t[0] ne $last_chr) { # a different chr, print
 125           map { print $_->{L} if ($_->{F}) } @last;
 126           @last = ();
 127           $last_chr = $t[0];
 128         }
 129         # The following block implemented by Nathans Weeks.
 130         push(@last, {L => $_, X => $t[1], F => 1}); # Enqueue current SNP
 131         if ($#last == $opts{N}) {                   # number of SNPs in queue is N+1
 132           if ($last[$#last]{X} - $last[0]{X} < $opts{W}) { # if all within window W
 133                 map {$_->{F} = 0} @last; # all SNPs in the window of size W are "bad"
 134           }
 135           print STDOUT $last[0]{L} if ($last[0]{F}); # print first SNP if good
 136           shift @last # dequeue first SNP
 137         }
 138   }
 139   # print the last few lines if applicable
 140   map { print $_->{L} if ($_->{F}) } @last;
 141 }
 142
 143 sub pileup2fq {
 144   my %opts = (d=>3, D=>255, Q=>25, G=>50, l=>10);
 145   getopts('d:D:Q:G:l:', \%opts);
 146   die(qq/
 147 Usage:   samtools.pl pileup2fq [options] <in.cns-pileup>
 148
 149 Options: -d INT    minimum depth        [$opts{d}]
 150          -D INT    maximum depth        [$opts{D}]
 151          -Q INT    min RMS mapQ         [$opts{Q}]
 152          -G INT    minimum indel score  [$opts{G}]
 153          -l INT    indel filter winsize [$opts{l}]
 154 /) if (@ARGV == 0 && -t STDIN);
 155
 156   my ($last_chr, $seq, $qual, @gaps, $last_pos);
 157   my $_Q = $opts{Q};
 158   my $_d = $opts{d};
 159   my $_D = $opts{D};
 160
 161   $last_chr = '';
 162   while (<>) {
 163         my @t = split;
 164         if ($last_chr ne $t[0]) {
 165           &p2q_post_process($last_chr, \$seq, \$qual, \@gaps, $opts{l}) if ($last_chr);
 166           $last_chr = $t[0];
 167           $last_pos = 0;
 168           $seq = ''; $qual = '';
 169           @gaps = ();
 170         }
 171         if ($t[1] - $last_pos != 1) {
 172           $seq .= 'n' x ($t[1] - $last_pos - 1);
 173           $qual .= '!' x ($t[1] - $last_pos - 1);
 174         }
 175         if ($t[2] eq '*') {
 176           push(@gaps, $t[1]) if ($t[5] >= $opts{G});
 177         } else {
 178           $seq .= ($t[6] >= $_Q && $t[7] >= $_d && $t[7] <= $_D)? uc($t[3]) : lc($t[3]);
 179           my $q = $t[4] + 33;
 180           $q = 126 if ($q > 126);
 181           $qual .= chr($q);
 182         }
 183         $last_pos = $t[1];
 184   }
 185   &p2q_post_process($last_chr, \$seq, \$qual, \@gaps, $opts{l});
 186 }
 187
 188 sub p2q_post_process {
 189   my ($chr, $seq, $qual, $gaps, $l) = @_;
 190   &p2q_filter_gaps($seq, $gaps, $l);
 191   print "\@$chr\n"; &p2q_print_str($seq);
 192   print "+\n"; &p2q_print_str($qual);
 193 }
 194
 195 sub p2q_filter_gaps {
 196   my ($seq, $gaps, $l) = @_;
 197   for my $g (@$gaps) {
 198         my $x = $g > $l? $g - $l : 0;
 199         substr($$seq, $x, $l + $l) = lc(substr($$seq, $x, $l + $l));
 200   }
 201 }
 202
 203 sub p2q_print_str {
 204   my ($s) = @_;
 205   my $l = length($$s);
 206   for (my $i = 0; $i < $l; $i += 60) {
 207         print substr($$s, $i, 60), "\n";
 208   }
 209 }
 210
 211 #
 212 # Usage
 213 #
 214
 215 sub usage {
 216   die(qq/
 217 Program: samtools.pl (helper script for SAMtools)
 218 Version: $version
 219 Contact: Heng Li <lh3\@sanger.ac.uk>\n
 220 Usage:   samtools.pl <command> [<arguments>]\n
 221 Command: indelFilter   filter indels generated by `pileup -c'
 222          snpFilter     filter SNPs generated by `pileup -c'
 223          pileup2fq     generate fastq from `pileup -c'
 224          showALEN      print alignment length (ALEN) following CIGAR
 225 \n/);
 226 }