Release samtools-0.1.15 (r949:203)

[samtools.git] / bcftools / vcfutils.pl

diff --git a/bcftools/vcfutils.pl b/bcftools/vcfutils.pl
index 499a4e4ee5c8fd495898a0312ffbdc4765cab47b..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";
   }
 }
 
@@ -455,7 +460,7 @@ sub hapmap2vcf {
 }
 
 sub vcf2fq {
-  my %opts = (d=>3, D=>100000, Q=>10, l=>10);
+  my %opts = (d=>3, D=>100000, Q=>10, l=>5);
   getopts('d:D:Q:l:', \%opts);
   die(qq/
 Usage:   vcfutils.pl vcf2fq [options] <all-site.vcf>
@@ -484,11 +489,12 @@ Options: -d INT    minimum depth          [$opts{d}]
          $seq = $qual = '';
          @gaps = ();
        }
-       if ($t[1] - $last_pos != 1) {
+       die("[vcf2fq] unsorted input\n") if ($t[1] - $last_pos < 0);
+       if ($t[1] - $last_pos > 1) {
          $seq .= 'n' x ($t[1] - $last_pos - 1);
          $qual .= '!' x ($t[1] - $last_pos - 1);
        }
-       if (length($t[3]) == 1 && $t[4] =~ /^([A-Za-z.])(,[A-Za-z])*$/) { # a SNP or reference
+       if (length($t[3]) == 1 && $t[7] !~ /INDEL/ && $t[4] =~ /^([A-Za-z.])(,[A-Za-z])*$/) { # a SNP or reference
          my ($ref, $alt) = ($t[3], $1);
          my ($b, $q);
          $q = $1 if ($t[7] =~ /FQ=(-?[\d\.]+)/);
@@ -506,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];
@@ -546,7 +552,7 @@ Command: subsam       get a subset of samples
          ucscsnp2vcf  convert UCSC SNP SQL dump to VCF
 
          varFilter    filtering short variants (*)
-         vcf2fq       VCF->fastq (*)
+         vcf2fq       VCF->fastq (**)
 
 Notes: Commands with description endting with (*) may need bcftools
        specific annotations.