From: Petr Danecek <pd3@sanger.ac.uk>
Date: Wed, 17 Feb 2010 13:32:53 +0000 (+0000)
Subject: Fixed filtering of SNPs near indels. Added min indel and SNP quality filter.
X-Git-Url: https://git.donarmstrong.com/?a=commitdiff_plain;h=88abd40bb30422736fc8f1aa1f40ea6ece050bac;p=samtools.git

Fixed filtering of SNPs near indels. Added min indel and SNP quality filter.
---

diff --git a/misc/samtools.pl b/misc/samtools.pl
index e6dd8e8..9f48b8f 100755
--- a/misc/samtools.pl
+++ b/misc/samtools.pl
@@ -46,10 +46,12 @@ sub showALEN {
 # G close to a high-quality indel (SNP only)
 # Q low RMS mapping quality (SNP only)
 # g close to another indel with higher quality (indel only)
+# s low SNP quality (SNP only)
+# i low indel quality (indel only)
 
 sub varFilter {
-  my %opts = (d=>3, D=>100, l=>30, Q=>25, q=>10, G=>25, s=>100, w=>10, W=>10, N=>2, p=>undef);
-  getopts('pq:d:D:l:Q:w:W:N:G:', \%opts);
+  my %opts = (d=>3, D=>100, l=>30, Q=>25, q=>10, G=>25, s=>100, w=>10, W=>10, N=>2, p=>undef, S=>'', i=>'');
+  getopts('pq:d:D:l:Q:w:W:N:G:S:i:', \%opts);
   die(qq/
 Usage:   samtools.pl varFilter [options] <in.cns-pileup>
 
@@ -57,6 +59,8 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
          -q INT    minimum RMS mapping quality for gaps [$opts{q}]
          -d INT    minimum read depth [$opts{d}]
          -D INT    maximum read depth [$opts{D}]
+         -S INT    minimum SNP quality [$opts{S}]
+         -i INT    minimum indel quality [$opts{i}]
 
          -G INT    min indel score for nearby SNP filtering [$opts{G}]
          -w INT    SNP within INT bp around a gap to be filtered [$opts{w}]
@@ -80,7 +84,8 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
 	next if (uc($t[2]) eq uc($t[3]) || $t[3] eq '*/*'); # skip non-var sites
 	# clear the out-of-range elements
 	while (@staging) {
-	  last if ($staging[0][2] eq $t[0] && $staging[0][3] + $max_dist >= $t[1]);
+      # Still on the same chromosome and the first element's window still affects this position?  
+	  last if ($staging[0][3] eq $t[0] && $staging[0][4] + $staging[0][2] + $max_dist >= $t[1]);
 	  varFilter_aux(shift(@staging), $opts{p}); # calling a function is a bit slower, not much
 	}
 	my ($flt, $score) = (0, -1);
@@ -90,14 +95,32 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
 	} elsif ($t[7] > $opts{D}) {
 	  $flt = 3;
 	}
+    if ($t[2] eq '*') { # an indel
+        if ($opts{i} && $opts{i}>$t[5]) { $flt = 8; }
+    }
+    elsif ($opts{S} && $opts{S}>$t[5]) { $flt = 7; }    # SNP
+
 	# site dependent filters
+    my $len=0;
 	if ($flt == 0) {
 	  if ($t[2] eq '*') { # an indel
+        
+        # If deletion, remember the length of the deletion
+        my ($a,$b) = split(m{/},$t[3]);
+        my $alen = length($a) - 1;
+        my $blen = length($b) - 1;
+        if ( $alen>$blen )
+        {
+            if ( substr($a,0,1) eq '-' ) { $len=$alen; }
+        }
+        elsif ( substr($b,0,1) eq '-' ) { $len=$blen; }
+
 		$flt = 1 if ($t[6] < $opts{q});
 		# filtering SNPs
 		if ($t[5] >= $opts{G}) {
 		  for my $x (@staging) {
-			next if ($x->[0] >= 0 || $x->[3] + $ow < $t[1]);
+            # Is it a SNP and is it outside the SNP filter window?
+			next if ($x->[0] >= 0 || $x->[4] + $x->[2] + $ow < $t[1]);
 			$x->[1] = 5 if ($x->[1] == 0);
 		  }
 		}
@@ -107,7 +130,8 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
 		$score += $opts{s} * $t[11] if ($t[9] ne '*');
 		# check the staging list for indel filtering
 		for my $x (@staging) {
-		  next if ($x->[0] < 0 || $x->[3] + $ol < $t[1]);
+          # Is it a SNP and is it outside the gap filter window
+		  next if ($x->[0] < 0 || $x->[4] + $x->[2] + $ol < $t[1]);
 		  if ($x->[0] < $score) {
 			$x->[1] = 6;
 		  } else {
@@ -119,17 +143,17 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
 		# check adjacent SNPs
 		my $k = 1;
 		for my $x (@staging) {
-		  ++$k if ($x->[0] < 0 && $x->[3] + $oW >= $t[1] && ($x->[1] == 0 || $x->[1] == 4 || $x->[1] == 5));
+		  ++$k if ($x->[0] < 0 && $x->[4] + $x->[2] + $oW >= $t[1] && ($x->[1] == 0 || $x->[1] == 4 || $x->[1] == 5));
 		}
 		# filtering is necessary
 		if ($k > $opts{N}) {
 		  $flt = 4;
 		  for my $x (@staging) {
-			 $x->[1] = 4 if ($x->[0] < 0 && $x->[3] + $oW >= $t[1] && $x->[1] == 0);
+			 $x->[1] = 4 if ($x->[0] < 0 && $x->[4] + $x->[2] + $oW >= $t[1] && $x->[1] == 0);
 		  }
 		} else { # then check gap filter
 		  for my $x (@staging) {
-			next if ($x->[0] < 0 || $x->[3] + $ow < $t[1]);
+			next if ($x->[0] < 0 || $x->[4] + $x->[2] + $ow < $t[1]);
 			if ($x->[0] >= $opts{G}) {
 			  $flt = 5; last;
 			}
@@ -137,7 +161,7 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
 		}
 	  }
 	}
-	push(@staging, [$score, $flt, @t]);
+	push(@staging, [$score, $flt, $len, @t]);
   }
   # output the last few elements in the staging list
   while (@staging) {
@@ -148,9 +172,9 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
 sub varFilter_aux {
   my ($first, $is_print) = @_;
   if ($first->[1] == 0) {
-	print join("\t", @$first[2 .. @$first-1]), "\n";
+	print join("\t", @$first[3 .. @$first-1]), "\n";
   } elsif ($is_print) {
-	print STDERR join("\t", substr("UQdDWGgX", $first->[1], 1), @$first[2 .. @$first-1]), "\n";
+	print STDERR join("\t", substr("UQdDWGgsiX", $first->[1], 1), @$first[3 .. @$first-1]), "\n";
   }
 }