projects / samtools.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
Imported Upstream version 0.1.10
[samtools.git] / bcftools / vcfutils.pl
diff --git a/bcftools/vcfutils.pl b/bcftools/vcfutils.pl
index d0b797143cf56fe0a67383ec4f610096eaecf2c0..bbc479bfffa47aa2b1089b5b6e56414b696be413 100755 (executable)
--- a/bcftools/vcfutils.pl
+++ b/bcftools/vcfutils.pl
@@ -10,11 +10,11 @@ use Getopt::Std;
 exit;
 
 sub main {
 exit;
 
 sub main {
-  my $version = '0.1.0';
   &usage if (@ARGV < 1);
   my $command = shift(@ARGV);
   my %func = (subsam=>\&subsam, listsam=>\&listsam, fillac=>\&fillac, qstats=>\&qstats, varFilter=>\&varFilter,
   &usage if (@ARGV < 1);
   my $command = shift(@ARGV);
   my %func = (subsam=>\&subsam, listsam=>\&listsam, fillac=>\&fillac, qstats=>\&qstats, varFilter=>\&varFilter,
-                         hapmap2vcf=>\&hapmap2vcf, ucscsnp2vcf=>\&ucscsnp2vcf, filter4vcf=>\&filter4vcf, ldstats=>\&ldstats);
+                         hapmap2vcf=>\&hapmap2vcf, ucscsnp2vcf=>\&ucscsnp2vcf, filter4vcf=>\&varFilter, ldstats=>\&ldstats,
+                         gapstats=>\&gapstats);
   die("Unknown command \"$command\".\n") if (!defined($func{$command}));
   &{$func{$command}};
 }
   die("Unknown command \"$command\".\n") if (!defined($func{$command}));
   &{$func{$command}};
 }
@@ -155,7 +155,7 @@ Note: This command discards indels. Output: QUAL #non-indel #SNPs #transitions #
        next if (length($t[3]) != 1 || uc($t[3]) eq 'N');
        $t[3] = uc($t[3]); $t[4] = uc($t[4]);
        my @s = split(',', $t[4]);
        next if (length($t[3]) != 1 || uc($t[3]) eq 'N');
        $t[3] = uc($t[3]); $t[4] = uc($t[4]);
        my @s = split(',', $t[4]);
-       $t[5] = 3 if ($t[5] < 0);
+       $t[5] = 3 if ($t[5] eq '.' || $t[5] < 0);
        next if (length($s[0]) != 1);
        my $hit;
        if ($is_vcf) {
        next if (length($s[0]) != 1);
        my $hit;
        if ($is_vcf) {
@@ -199,37 +199,38 @@ Note: This command discards indels. Output: QUAL #non-indel #SNPs #transitions #
 }
 
 sub varFilter {
 }
 
 sub varFilter {
-  my %opts = (d=>1, D=>10000, l=>30, Q=>25, q=>10, G=>25, s=>100, w=>10, W=>10, N=>2, p=>undef, F=>.001);
-  getopts('pq:d:D:l:Q:w:W:N:G:F:', \%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);
+  getopts('pd:D:W:Q:w:a:1:2:3:4:', \%opts);
   die(qq/
 Usage:   vcfutils.pl varFilter [options] <in.vcf>
 
 Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
   die(qq/
 Usage:   vcfutils.pl varFilter [options] <in.vcf>
 
 Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
-         -q INT    minimum RMS mapping quality for gaps [$opts{q}]
          -d INT    minimum read depth [$opts{d}]
          -D INT    maximum read depth [$opts{D}]
          -d INT    minimum read depth [$opts{d}]
          -D INT    maximum read depth [$opts{D}]
-
-         -G INT    min indel score for nearby SNP filtering [$opts{G}]
+         -a INT    minimum number of alternate bases [$opts{a}]
          -w INT    SNP within INT bp around a gap to be filtered [$opts{w}]
          -w INT    SNP within INT bp around a gap to be filtered [$opts{w}]
-
-         -W INT    window size for filtering dense SNPs [$opts{W}]
-         -N INT    max number of SNPs in a window [$opts{N}]
-
-         -l INT    window size for filtering adjacent gaps [$opts{l}]
-
+         -W INT    window size for filtering adjacent gaps [$opts{W}]
+         -1 FLOAT  min P-value for strand bias (given PV4) [$opts{1}]
+         -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}]
          -p        print filtered variants
          -p        print filtered variants
+
+Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools.
 \n/) if (@ARGV == 0 && -t STDIN);
 
   # calculate the window size
 \n/) if (@ARGV == 0 && -t STDIN);
 
   # calculate the window size
-  my ($ol, $ow, $oW) = ($opts{l}, $opts{w}, $opts{W});
+  my ($ol, $ow) = ($opts{W}, $opts{w});
   my $max_dist = $ol > $ow? $ol : $ow;
   my $max_dist = $ol > $ow? $ol : $ow;
-  $max_dist = $oW if ($max_dist < $oW);
   # the core loop
   # the core loop
-  my @staging; # (indel_filtering_score, flt_tag)
+  my @staging; # (indel_filtering_score, flt_tag, indel_span; chr, pos, ...)
   while (<>) {
        my @t = split;
   while (<>) {
        my @t = split;
-       next if (/^#/);
+    if (/^#/) {
+         print; next;
+       }
        next if ($t[4] eq '.'); # skip non-var sites
        next if ($t[4] eq '.'); # skip non-var sites
+       # check if the site is a SNP
        my $is_snp = 1;
        if (length($t[3]) > 1) {
          $is_snp = 0;
        my $is_snp = 1;
        if (length($t[3]) > 1) {
          $is_snp = 0;
@@ -245,24 +246,18 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
          last if ($staging[0][3] eq $t[0] && $staging[0][4] + $staging[0][2] + $max_dist >= $t[1]);
          varFilter_aux(shift(@staging), $opts{p}); # calling a function is a bit slower, not much
        }
          last if ($staging[0][3] eq $t[0] && $staging[0][4] + $staging[0][2] + $max_dist >= $t[1]);
          varFilter_aux(shift(@staging), $opts{p}); # calling a function is a bit slower, not much
        }
-       my ($flt, $score) = (0, -1);
-
-       # collect key annotations
-       my ($dp, $mq, $af) = (-1, -1, 1);
-       if ($t[7] =~ /DP=(\d+)/i) {
-         $dp = $1;
-       } elsif ($t[7] =~ /DP4=(\d+),(\d+),(\d+),(\d+)/i) {
+       my $flt = 0;
+       # parse annotations
+       my ($dp, $mq, $dp_alt) = (-1, -1, -1);
+       if ($t[7] =~ /DP4=(\d+),(\d+),(\d+),(\d+)/i) {
          $dp = $1 + $2 + $3 + $4;
          $dp = $1 + $2 + $3 + $4;
+         $dp_alt = $3 + $4;
        }
        }
-       if ($t[7] =~ /MQ=(\d+)/i) {
-         $mq = $1;
-       }
-       if ($t[7] =~ /AF=([^\s;=]+)/i) {
-         $af = $1;
-       } elsif ($t[7] =~ /AF1=([^\s;=]+)/i) {
-         $af = $1;
+       if ($t[7] =~ /DP=(\d+)/i) {
+         $dp = $1;
        }
        }
-       # the depth filter
+       $mq = $1 if ($t[7] =~ /MQ=(\d+)/i);
+       # the depth and mapQ filter
        if ($dp >= 0) {
          if ($dp < $opts{d}) {
                $flt = 2;
        if ($dp >= 0) {
          if ($dp < $opts{d}) {
                $flt = 2;
@@ -270,58 +265,40 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
                $flt = 3;
          }
        }
                $flt = 3;
          }
        }
+       $flt = 4 if ($dp_alt >= 0 && $dp_alt < $opts{a});
+       $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}));
 
        # site dependent filters
 
        # site dependent filters
-       my $dlen = 0;
+       my ($rlen, $indel_score) = (0, -1); # $indel_score<0 for SNPs
        if ($flt == 0) {
          if (!$is_snp) { # an indel
        if ($flt == 0) {
          if (!$is_snp) { # an indel
-        # If deletion, remember the length of the deletion
-               $dlen = length($t[3]) - 1;
-               $flt = 1 if ($mq < $opts{q});
+               $rlen = length($t[3]) - 1;
+               $indel_score = $t[5] * 100 + $dp_alt;
                # filtering SNPs
                # filtering SNPs
-               if ($t[5] >= $opts{G}) {
-                 for my $x (@staging) {
-            # Is it a SNP and is it outside the SNP filter window?
-                       next if ($x->[0] >= 0 || $x->[4] + $x->[2] + $ow < $t[1]);
-                       $x->[1] = 5 if ($x->[1] == 0);
-                 }
+               for my $x (@staging) {
+                 next if ($x->[0] >= 0 || $x->[1] || $x->[4] + $x->[2] + $ow < $t[1]);
+                 $x->[1] = 5;
                }
                }
-               # the indel filtering score
-               $score = $t[5];
                # check the staging list for indel filtering
                for my $x (@staging) {
                # check the staging list for indel filtering
                for my $x (@staging) {
-          # Is it a SNP and is it outside the gap filter window
-                 next if ($x->[0] < 0 || $x->[4] + $x->[2] + $ol < $t[1]);
-                 if ($x->[0] < $score) {
+                 next if ($x->[0] < 0 || $x->[1] || $x->[4] + $x->[2] + $ol < $t[1]);
+                 if ($x->[0] < $indel_score) {
                        $x->[1] = 6;
                  } else {
                        $flt = 6; last;
                  }
                }
          } else { # a SNP
                        $x->[1] = 6;
                  } else {
                        $flt = 6; last;
                  }
                }
          } else { # a SNP
-               $flt = 1 if ($mq < $opts{Q});
-               # check adjacent SNPs
-               my $k = 1;
                for my $x (@staging) {
                for my $x (@staging) {
-                 ++$k if ($x->[0] < 0 && -($x->[0] + 1) > $opts{F} && $x->[4] + $x->[2] + $oW >= $t[1] && ($x->[1] == 0 || $x->[1] == 4 || $x->[1] == 5));
-               }
-               # filtering is necessary
-               if ($k > $opts{N}) {
-                 $flt = 4;
-                 for my $x (@staging) {
-                        $x->[1] = 4 if ($x->[0] < 0 && $x->[4] + $x->[2] + $oW >= $t[1] && $x->[1] == 0);
-                 }
-               } else { # then check gap filter
-                 for my $x (@staging) {
-                       next if ($x->[0] < 0 || $x->[4] + $x->[2] + $ow < $t[1]);
-                       if ($x->[0] >= $opts{G}) {
-                         $flt = 5; last;
-                       }
-                 }
+                 next if ($x->[0] < 0 || $x->[1] || $x->[4] + $x->[2] + $ow < $t[1]);
+                 $flt = 5;
+                 last;
                }
          }
        }
                }
          }
        }
-       push(@staging, [$score < 0? -$af-1 : $score, $flt, $dlen, @t]);
+       push(@staging, [$indel_score, $flt, $rlen, @t]);
   }
   # output the last few elements in the staging list
   while (@staging) {
   }
   # output the last few elements in the staging list
   while (@staging) {
@@ -334,47 +311,35 @@ sub varFilter_aux {
   if ($first->[1] == 0) {
        print join("\t", @$first[3 .. @$first-1]), "\n";
   } elsif ($is_print) {
   if ($first->[1] == 0) {
        print join("\t", @$first[3 .. @$first-1]), "\n";
   } elsif ($is_print) {
-       print STDERR join("\t", substr("UQdDWGgsiX", $first->[1], 1), @$first[3 .. @$first-1]), "\n";
+       print STDERR join("\t", substr("UQdDaGgP", $first->[1], 1), @$first[3 .. @$first-1]), "\n";
   }
 }
 
   }
 }
 
-sub filter4vcf {
-  my %opts = (d=>3, D=>2000, 1=>1e-4, 2=>1e-100, 3=>0, 4=>1e-4, Q=>10, q=>3);
-  getopts('d:D:1:2:3:4:Q:q:', \%opts);
-  die(qq/
-Usage:   vcfutils.pl filter4vcf [options] <in.vcf>
-
-Options: -d INT     min total depth (given DP or DP4) [$opts{d}]
-         -D INT     max total depth [$opts{D}]
-         -q INT     min SNP quality [$opts{q}]
-         -Q INT     min RMS mapQ (given MQ) [$opts{Q}]
-         -1 FLOAT   min P-value for strand bias (given PV4) [$opts{1}]
-         -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}]\n
-/) if (@ARGV == 0 && -t STDIN);
-
-  my %ts = (AG=>1, GA=>1, CT=>1, TC=>1);
-
-  my @n = (0, 0);
+sub gapstats {
+  my (@c0, @c1);
+  $c0[$_] = $c1[$_] = 0 for (0 .. 10000);
   while (<>) {
   while (<>) {
-       if (/^#/) {
-         print;
-         next;
-       }
-       next if (/PV4=([^,]+),([^,]+),([^,]+),([^,;\t]+)/ && ($1<$opts{1} || $2<$opts{2} || $3<$opts{3} || $4<$opts{4}));
-       my $depth = -1;
-       $depth = $1 if (/DP=(\d+)/);
-       $depth = $1+$2+$3+$4 if (/DP4=(\d+),(\d+),(\d+),(\d+)/);
-       next if ($depth > 0 && ($depth < $opts{d} || $depth > $opts{D}));
-       next if (/MQ=(\d+)/ && $1 < $opts{Q});
+       next if (/^#/);
        my @t = split;
        my @t = split;
-       next if ($t[5] >= 0 && $t[5] < $opts{q});
-       ++$n[0];
+       next if (length($t[3]) == 1 && $t[4] =~ /^[A-Za-z](,[A-Za-z])*$/); # not an indel
        my @s = split(',', $t[4]);
        my @s = split(',', $t[4]);
-       ++$n[1] if ($ts{$t[3].$s[0]});
-       print;
+       for my $x (@s) {
+         my $l = length($x) - length($t[3]) + 5000;
+         if ($x =~ /^-/) {
+               $l = -(length($x) - 1) + 5000;
+         } elsif ($x =~ /^\+/) {
+               $l = length($x) - 1 + 5000;
+         }
+         $c0[$l] += 1 / @s;
+       }
+  }
+  for (my $i = 0; $i < 10000; ++$i) {
+       next if ($c0[$i] == 0);
+       $c1[0] += $c0[$i];
+       $c1[1] += $c0[$i] if (($i-5000)%3 == 0);
+       printf("C\t%d\t%.2f\n", ($i-5000), $c0[$i]);
   }
   }
+  printf("3\t%d\t%d\t%.3f\n", $c1[0], $c1[1], $c1[1]/$c1[0]);
 }
 
 sub ucscsnp2vcf {
 }
 
 sub ucscsnp2vcf {
@@ -470,7 +435,6 @@ Command: subsam       get a subset of samples
          fillac       fill the allele count field
          qstats       SNP stats stratified by QUAL
          varFilter    filtering short variants
          fillac       fill the allele count field
          qstats       SNP stats stratified by QUAL
          varFilter    filtering short variants
-         filter4vcf   filtering VCFs produced by samtools+bcftools
          hapmap2vcf   convert the hapmap format to VCF
          ucscsnp2vcf  convert UCSC SNP SQL dump to VCF
 \n/);
          hapmap2vcf   convert the hapmap format to VCF
          ucscsnp2vcf  convert UCSC SNP SQL dump to VCF
 \n/);
woldlab customizations/backport of samtools package