* modify the filtering rule: also filter SNPs around filtered indels

 * added MNP filter

diff --git a/bcftools/vcfutils.pl b/bcftools/vcfutils.pl
index bbc479bfffa47aa2b1089b5b6e56414b696be413..6cc168fb363f59541511429b6e6538e0263acbe1 100755 (executable)
--- a/bcftools/vcfutils.pl
+++ b/bcftools/vcfutils.pl
@@ -231,13 +231,17 @@ Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools.
        }
        next if ($t[4] eq '.'); # skip non-var sites
        # check if the site is a SNP
-       my $is_snp = 1;
+       my $type = 1; # SNP
        if (length($t[3]) > 1) {
-         $is_snp = 0;
+         $type = 2; # MNP
+         my @s = split(',', $t[4]);
+         for (@s) {
+               $type = 3 if (length != length($t[3]));
+         }
        } else {
          my @s = split(',', $t[4]);
          for (@s) {
-               $is_snp = 0 if (length > 1);
+               $type = 3 if (length > 1);
          }
        }
        # clear the out-of-range elements
@@ -270,35 +274,42 @@ Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools.
        $flt = 7 if ($flt == 0 && /PV4=([^,]+),([^,]+),([^,]+),([^,;\t]+)/
                                 && ($1<$opts{1} || $2<$opts{2} || $3<$opts{3} || $4<$opts{4}));
 
-       # site dependent filters
-       my ($rlen, $indel_score) = (0, -1); # $indel_score<0 for SNPs
+       my $score = $t[5] * 100 + $dp_alt;
+       my $rlen = length($t[3]) - 1; # $indel_score<0 for SNPs
        if ($flt == 0) {
-         if (!$is_snp) { # an indel
-               $rlen = length($t[3]) - 1;
-               $indel_score = $t[5] * 100 + $dp_alt;
-               # filtering SNPs
+         if ($type == 3) { # an indel
+               # filtering SNPs and MNPs
                for my $x (@staging) {
-                 next if ($x->[0] >= 0 || $x->[1] || $x->[4] + $x->[2] + $ow < $t[1]);
+                 next if (($x->[0]&3) == 3 || $x->[1] || $x->[4] + $x->[2] + $ow < $t[1]);
                  $x->[1] = 5;
                }
                # check the staging list for indel filtering
                for my $x (@staging) {
-                 next if ($x->[0] < 0 || $x->[1] || $x->[4] + $x->[2] + $ol < $t[1]);
-                 if ($x->[0] < $indel_score) {
+                 next if (($x->[0]&3) != 3 || $x->[1] || $x->[4] + $x->[2] + $ol < $t[1]);
+                 if ($x->[0]>>2 < $score) {
                        $x->[1] = 6;
                  } else {
                        $flt = 6; last;
                  }
                }
-         } else { # a SNP
+         } else { # SNP or MNP
                for my $x (@staging) {
-                 next if ($x->[0] < 0 || $x->[1] || $x->[4] + $x->[2] + $ow < $t[1]);
+                 next if (($x->[0]&3) != 3 || $x->[4] + $x->[2] + $ow < $t[1]);
                  $flt = 5;
                  last;
                }
+               # check MNP
+               for my $x (@staging) {
+                 next if (($x->[0]&3) == 3 || $x->[4] + $x->[2] < $t[1]);
+                 if ($x->[0]>>2 < $score) {
+                       $x->[1] = 8;
+                 } else {
+                       $flt = 8; last;
+                 }
+               }
          }
        }
-       push(@staging, [$indel_score, $flt, $rlen, @t]);
+       push(@staging, [$score<<2|$type, $flt, $rlen, @t]);
   }
   # output the last few elements in the staging list
   while (@staging) {
@@ -311,7 +322,7 @@ sub varFilter_aux {
   if ($first->[1] == 0) {
        print join("\t", @$first[3 .. @$first-1]), "\n";
   } elsif ($is_print) {
-       print STDERR join("\t", substr("UQdDaGgP", $first->[1], 1), @$first[3 .. @$first-1]), "\n";
+       print STDERR join("\t", substr("UQdDaGgPM", $first->[1], 1), @$first[3 .. @$first-1]), "\n";
   }
 }