* samtools-0.1.9-16 (r822)

[samtools.git] / bcftools / vcfutils.pl

diff --git a/bcftools/vcfutils.pl b/bcftools/vcfutils.pl
index 1270c6c7ff5bf8cdc6d615c507b282a8680a89c9..bbc479bfffa47aa2b1089b5b6e56414b696be413 100755 (executable)
--- a/bcftools/vcfutils.pl
+++ b/bcftools/vcfutils.pl
@@ -249,12 +249,13 @@ Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools.
        my $flt = 0;
        # parse annotations
        my ($dp, $mq, $dp_alt) = (-1, -1, -1);
-       if ($t[7] =~ /DP=(\d+)/i) {
-         $dp = $1;
-       } elsif ($t[7] =~ /DP4=(\d+),(\d+),(\d+),(\d+)/i) {
+       if ($t[7] =~ /DP4=(\d+),(\d+),(\d+),(\d+)/i) {
          $dp = $1 + $2 + $3 + $4;
          $dp_alt = $3 + $4;
        }
+       if ($t[7] =~ /DP=(\d+)/i) {
+         $dp = $1;
+       }
        $mq = $1 if ($t[7] =~ /MQ=(\d+)/i);
        # the depth and mapQ filter
        if ($dp >= 0) {
@@ -324,13 +325,21 @@ sub gapstats {
        my @s = split(',', $t[4]);
        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);
-       printf("%d\t%.2f\n", ($i-5000), $c0[$i]);
+       $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 {