* removed a comment line in kaln.c

[samtools.git] / bcftools / vcfutils.pl

diff --git a/bcftools/vcfutils.pl b/bcftools/vcfutils.pl
index 6fb1952cc6e578524078d6dcef12b65cd760a4ed..cdb96cb53ee287fbbc3b541baa50c7778ec176bb 100755 (executable)
--- a/bcftools/vcfutils.pl
+++ b/bcftools/vcfutils.pl
@@ -70,12 +70,21 @@ sub fillac {
          print;
        } else {
          my @t = split;
-         my @c;
+         my @c = (0);
          my $n = 0;
-         $c[1] = 0;
+         my $s = -1;
+         @_ = split(":", $t[8]);
+         for (0 .. $#_) {
+               if ($_[$_] eq 'GT') { $s = $_; last; }
+         }
+         if ($s < 0) {
+               print join("\t", @t), "\n";
+               next;
+         }
          for (9 .. $#t) {
-               if ($t[$_] =~ /^(\d+).(\d+)/) {
-                 ++$c[$1]; ++$c[$2];
+               if ($t[$_] =~ /^0,0,0/) {
+               } elsif ($t[$_] =~ /^([^\s:]+:){$s}(\d+).(\d+)/) {
+                 ++$c[$2]; ++$c[$3];
                  $n += 2;
                }
          }
@@ -95,7 +104,7 @@ sub fillac {
 }
 
 sub qstats {
-  my %opts = (r=>'', s=>0.01, v=>undef);
+  my %opts = (r=>'', s=>0.02, v=>undef);
   getopts('r:s:v', \%opts);
   die("Usage: vcfutils.pl qstats [-r ref.vcf] <in.vcf>\n
 Note: This command discards indels. Output: QUAL #non-indel #SNPs #transitions #joint ts/tv #joint/#ref #joint/#non-indel \n") if (@ARGV == 0 && -t STDIN);
@@ -147,14 +156,20 @@ Note: This command discards indels. Output: QUAL #non-indel #SNPs #transitions #
   my $next = $opts{s};
   my $last = $a[0];
   my @c = (0, 0, 0, 0);
+  my @lc;
+  $lc[1] = $lc[2] = 0;
   for my $p (@a) {
        if ($p->[0] == -1 || ($p->[0] != $last && $c[0]/@a > $next)) {
          my @x;
          $x[0] = sprintf("%.4f", $c[1]-$c[2]? $c[2] / ($c[1] - $c[2]) : 100);
          $x[1] = sprintf("%.4f", $hsize? $c[3] / $hsize : 0);
          $x[2] = sprintf("%.4f", $c[3] / $c[1]);
+         my $a = $c[1] - $lc[1];
+         my $b = $c[2] - $lc[2];
+         $x[3] = sprintf("%.4f", $a-$b? $b / ($a-$b) : 100);
          print join("\t", $last, @c, @x), "\n";
          $next = $c[0]/@a + $opts{s};
+         $lc[1] = $c[1]; $lc[2] = $c[2];
        }
        ++$c[0]; $c[1] += $p->[1]; $c[2] += $p->[2]; $c[3] += $p->[3];
        $last = $p->[0];