compute max SP and max GQ from sample genotypes

[samtools.git] / bcftools / vcfutils.pl

diff --git a/bcftools/vcfutils.pl b/bcftools/vcfutils.pl
index 63509e820fb2c56dc255e4ec6dbc6cea5fd74abb..5df8b0c05e1275e7edc4b8063efa5c585ee7a80a 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=>10000, 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>
 
@@ -289,6 +289,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
@@ -338,7 +339,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";
   }
 }