* samtools-0.1.19-11 (r817)

[samtools.git] / bcftools / vcfutils.pl

diff --git a/bcftools/vcfutils.pl b/bcftools/vcfutils.pl
index c1177eacb9603cc2cb9f93ab27d54b5add1e0779..1270c6c7ff5bf8cdc6d615c507b282a8680a89c9 100755 (executable)
--- a/bcftools/vcfutils.pl
+++ b/bcftools/vcfutils.pl
@@ -13,7 +13,8 @@ sub main {
   &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=>\&varFilter, ldstats=>\&ldstats);
+                         hapmap2vcf=>\&hapmap2vcf, ucscsnp2vcf=>\&ucscsnp2vcf, filter4vcf=>\&varFilter, ldstats=>\&ldstats,
+                         gapstats=>\&gapstats);
   die("Unknown command \"$command\".\n") if (!defined($func{$command}));
   &{$func{$command}};
 }
@@ -313,6 +314,25 @@ sub varFilter_aux {
   }
 }
 
+sub gapstats {
+  my (@c0, @c1);
+  $c0[$_] = $c1[$_] = 0 for (0 .. 10000);
+  while (<>) {
+       next if (/^#/);
+       my @t = split;
+       next if (length($t[3]) == 1 && $t[4] =~ /^[A-Za-z](,[A-Za-z])*$/); # not an indel
+       my @s = split(',', $t[4]);
+       for my $x (@s) {
+         my $l = length($x) - length($t[3]) + 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]);
+  }
+}
+
 sub ucscsnp2vcf {
   die("Usage: vcfutils.pl <in.ucsc.snp>\n") if (@ARGV == 0 && -t STDIN);
   print "##fileformat=VCFv4.0\n";