* samtools.pl-0.2.0

 * snpFilter

diff --git a/misc/samtools.pl b/misc/samtools.pl
index 8f109e0e57374f96a56984070ec71ce5fb43f623..6fba900efce35bee3907a2725cee4c5127789d32 100755 (executable)
--- a/misc/samtools.pl
+++ b/misc/samtools.pl
@@ -6,11 +6,11 @@ use strict;
 use warnings;
 use Getopt::Std;
 
-my $version = '0.1.1';
+my $version = '0.2.0';
 &usage if (@ARGV < 1);
 
 my $command = shift(@ARGV);
-my %func = (indelFilter=>\&indelFilter, showALEN=>\&showALEN);
+my %func = (snpFilter=>\&snpFilter, indelFilter=>\&indelFilter, showALEN=>\&showALEN);
 
 die("Unknown command \"$command\".\n") if (!defined($func{$command}));
 &{$func{$command}};
@@ -76,10 +76,74 @@ Options: -D INT    maximum read depth [$opts{D}]
   print $last->[3] if (defined $last->[0]);
 }
 
+#
+# snpFilter
+#
+
+sub snpFilter {
+  my %opts = (f=>'', Q=>40, d=>3, w=>10, D=>0, N=>2, W=>10, q=>20, s=>50);
+  getopts('fs:w:q:Q:d:D:W:N:', \%opts);
+  die(qq{
+Usage:   samtools.pl snpFilter [options] <cns2snp.snp>
+
+Options: -d INT        minimum depth to call a SNP [$opts{d}]
+         -D INT        maximum depth, 0 to ignore [$opts{D}]
+         -Q INT        required max mapping quality of the reads covering the SNP [$opts{Q}]
+         -q INT        minimum SNP quality [$opts{q}]
+
+         -f FILE       filtered samtools indels [null]
+         -s INT        minimum samtols indel score [$opts{s}]
+         -w INT        SNP within INT bp around an indel to be filtered [$opts{w}]
+
+         -W INT        window size for filtering dense SNPs [$opts{W}]
+         -N INT        maximum number of SNPs in a window [$opts{N}]
+\n}) unless (@ARGV);
+  my (%hash, $fh);
+  my $skip = $opts{w};
+  $opts{D} = 100000000 if ($opts{D} == 0);
+  if ($opts{f}) { # filtered samtools indel
+       my $n = 0;
+       open($fh, $opts{f}) || die;
+       while (<$fh>) {
+         my @t = split;
+         next if ($t[2] ne '*' || $t[3] eq '*/*' || $t[5] < $opts{s});
+         for (my $x = $t[1] - $skip + 1; $x < $t[1] + $skip; ++$x) {
+               $hash{$t[0],$x} = 1;
+         }
+       }
+       close($fh);
+  }
+  my (@last, $last_chr);
+  $last_chr = '';
+  while (<>) {
+       my @t = split;
+       next if ($t[2] eq '*' || $hash{$t[0],$t[1]});
+       my $is_good = ($t[7] >= $opts{d} && $t[7] <= $opts{D} && $t[6] >= $opts{Q} && $t[5] >= $opts{q})? 1 : 0;
+       next unless ($is_good); # drop
+       if ($t[0] ne $last_chr) { # a different chr, print
+         map { print $_->{L} if ($_->{F}) } @last;
+         @last = ();
+         $last_chr = $t[0];
+       }
+       # The following block implemented by Nathans Weeks.
+       push(@last, {L => $_, X => $t[1], F => 1}); # Enqueue current SNP
+       if ($#last == $opts{N}) {                   # number of SNPs in queue is N+1
+         if ($last[$#last]{X} - $last[0]{X} < $opts{W}) { # if all within window W
+               map {$_->{F} = 0} @last; # all SNPs in the window of size W are "bad"
+         }
+         print STDOUT $last[0]{L} if ($last[0]{F}); # print first SNP if good
+         shift @last # dequeue first SNP
+       }
+  }
+  # print the last few lines if applicable
+  map { print $_->{L} if ($_->{F}) } @last;
+}
+
 sub usage {
   die(qq/
 Usage:   samtools.pl <command> [<arguments>]\n
 Command: indelFilter   filter indels generated by `pileup -c'
+         snpFilter     filter SNPs generated by `pileup -c'
          showALEN      print alignment length (ALEN) following CIGAR
 \n/);
 }