Release samtools-0.1.15 (r949:203)

[samtools.git] / bcftools / vcfutils.pl

diff --git a/bcftools/vcfutils.pl b/bcftools/vcfutils.pl
index 63509e820fb2c56dc255e4ec6dbc6cea5fd74abb..6ced66a2893e9b98884774e88607e9f2060af3d1 100755 (executable)
--- a/bcftools/vcfutils.pl
+++ b/bcftools/vcfutils.pl
@@ -215,8 +215,8 @@ Note: This command discards indels. Output: QUAL #non-indel #SNPs #transitions #
 }
 
 sub varFilter {
-  my %opts = (d=>2, D=>10000, a=>2, W=>10, Q=>10, w=>10, p=>undef, 1=>1e-4, 2=>1e-100, 3=>0, 4=>1e-4);
-  getopts('pd:D:W:Q:w:a:1:2:3:4:', \%opts);
+  my %opts = (d=>2, D=>10000000, a=>2, W=>10, Q=>10, w=>10, p=>undef, 1=>1e-4, 2=>1e-100, 3=>0, 4=>1e-4, G=>0, S=>1000);
+  getopts('pd:D:W:Q:w:a:1:2:3:4:G:S:', \%opts);
   die(qq/
 Usage:   vcfutils.pl varFilter [options] <in.vcf>
 
@@ -246,6 +246,7 @@ Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools.
          print; next;
        }
        next if ($t[4] eq '.'); # skip non-var sites
+    next if ($t[3] eq 'N'); # skip sites with unknown ref ('N')
        # check if the site is a SNP
        my $type = 1; # SNP
        if (length($t[3]) > 1) {
@@ -289,6 +290,7 @@ Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools.
        $flt = 1 if ($flt == 0 && $mq >= 0 && $mq < $opts{Q});
        $flt = 7 if ($flt == 0 && /PV4=([^,]+),([^,]+),([^,]+),([^,;\t]+)/
                                 && ($1<$opts{1} || $2<$opts{2} || $3<$opts{3} || $4<$opts{4}));
+       $flt = 8 if ($flt == 0 && ((/MXGQ=(\d+)/ && $1 < $opts{G}) || (/MXSP=(\d+)/ && $1 >= $opts{S})));
 
        my $score = $t[5] * 100 + $dp_alt;
        my $rlen = length($t[3]) - 1; # $indel_score<0 for SNPs
@@ -311,7 +313,10 @@ Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools.
          } else { # SNP or MNP
                for my $x (@staging) {
                  next if (($x->[0]&3) != 3 || $x->[4] + $x->[2] + $ow < $t[1]);
-                 $flt = 5;
+                 if ($x->[4] + length($x->[7]) - 1 == $t[1] && substr($x->[7], -1, 1) eq substr($t[4], 0, 1)
+                         && length($x->[7]) - length($x->[6]) == 1) {
+                       $x->[1] = 5;
+                 } else { $flt = 5; }
                  last;
                }
                # check MNP
@@ -338,7 +343,7 @@ sub varFilter_aux {
   if ($first->[1] == 0) {
        print join("\t", @$first[3 .. @$first-1]), "\n";
   } elsif ($is_print) {
-       print STDERR join("\t", substr("UQdDaGgPM", $first->[1], 1), @$first[3 .. @$first-1]), "\n";
+       print STDERR join("\t", substr("UQdDaGgPMS", $first->[1], 1), @$first[3 .. @$first-1]), "\n";
   }
 }
 
@@ -507,7 +512,7 @@ Options: -d INT    minimum depth          [$opts{d}]
          $q = chr($q <= 126? $q : 126);
          $seq .= $b;
          $qual .= $q;
-       } else { # an INDEL
+       } elsif ($t[4] ne '.') { # an INDEL
          push(@gaps, [$t[1], length($t[3])]);
        }
        $last_pos = $t[1];