X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=samtools.git;a=blobdiff_plain;f=bcftools%2Fvcfutils.pl;fp=bcftools%2Fvcfutils.pl;h=2b7ba0b1d00932be1f79f08e4b3bd8e8db951295;hp=6ced66a2893e9b98884774e88607e9f2060af3d1;hb=b301e959d73eee0955c57004f344f17af00703f4;hpb=8a8212c76299ff302cdd41be704d8e0308bd2b1b

diff --git a/bcftools/vcfutils.pl b/bcftools/vcfutils.pl
index 6ced66a..2b7ba0b 100755
--- a/bcftools/vcfutils.pl
+++ b/bcftools/vcfutils.pl
@@ -86,7 +86,7 @@ sub fillac {
 	  print;
 	} else {
 	  my @t = split;
-	  my @c = (0);
+	  my @c = (0, 0);
 	  my $n = 0;
 	  my $s = -1;
 	  @_ = split(":", $t[8]);
@@ -215,8 +215,8 @@ Note: This command discards indels. Output: QUAL #non-indel #SNPs #transitions #
 }
 
 sub varFilter {
-  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);
+  my %opts = (d=>2, D=>10000000, a=>2, W=>10, Q=>10, w=>3, p=>undef, 1=>1e-4, 2=>1e-100, 3=>0, 4=>1e-4, G=>0, S=>1000, e=>1e-4);
+  getopts('pd:D:W:Q:w:a:1:2:3:4:G:S:e:', \%opts);
   die(qq/
 Usage:   vcfutils.pl varFilter [options] <in.vcf>
 
@@ -230,6 +230,7 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
          -2 FLOAT  min P-value for baseQ bias [$opts{2}]
          -3 FLOAT  min P-value for mapQ bias [$opts{3}]
          -4 FLOAT  min P-value for end distance bias [$opts{4}]
+		 -e FLOAT  min P-value for HWE (plus F<0) [$opts{e}]
          -p        print filtered variants
 
 Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools.
@@ -291,6 +292,12 @@ 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}));
 	$flt = 8 if ($flt == 0 && ((/MXGQ=(\d+)/ && $1 < $opts{G}) || (/MXSP=(\d+)/ && $1 >= $opts{S})));
+	# HWE filter
+	if ($t[7] =~ /G3=([^;,]+),([^;,]+),([^;,]+).*HWE=([^;,]+)/ && $4 < $opts{e}) {
+		my $p = 2*$1 + $2;
+		my $f = ($p > 0 && $p < 1)? 1 - $2 / ($p * (1-$p)) : 0;
+		$flt = 9 if ($f < 0);
+	}
 
 	my $score = $t[5] * 100 + $dp_alt;
 	my $rlen = length($t[3]) - 1; # $indel_score<0 for SNPs
@@ -499,7 +506,7 @@ Options: -d INT    minimum depth          [$opts{d}]
 	  my ($b, $q);
 	  $q = $1 if ($t[7] =~ /FQ=(-?[\d\.]+)/);
 	  if ($q < 0) {
-		$_ = $1 if ($t[7] =~ /AF1=([\d\.]+)/);
+		$_ = ($t[7] =~ /AF1=([\d\.]+)/)? $1 : 0;
 		$b = ($_ < .5 || $alt eq '.')? $ref : $alt;
 		$q = -$q;
 	  } else {