bcftools/vcfutils.pl

   1 #!/usr/bin/perl -w
   2
   3 # Author: lh3
   4
   5 use strict;
   6 use warnings;
   7 use Getopt::Std;
   8
   9 &main;
  10 exit;
  11
  12 sub main {
  13   &usage if (@ARGV < 1);
  14   my $command = shift(@ARGV);
  15   my %func = (subsam=>\&subsam, listsam=>\&listsam, fillac=>\&fillac, qstats=>\&qstats, varFilter=>\&varFilter,
  16                           hapmap2vcf=>\&hapmap2vcf, ucscsnp2vcf=>\&ucscsnp2vcf, filter4vcf=>\&varFilter, ldstats=>\&ldstats,
  17                           gapstats=>\&gapstats, splitchr=>\&splitchr, vcf2fq=>\&vcf2fq);
  18   die("Unknown command \"$command\".\n") if (!defined($func{$command}));
  19   &{$func{$command}};
  20 }
  21
  22 sub splitchr {
  23   my %opts = (l=>5000000);
  24   getopts('l:', \%opts);
  25   my $l = $opts{l};
  26   die(qq/Usage: vcfutils.pl splitchr [-l $opts{l}] <in.fa.fai>\n/) if (@ARGV == 0 && -t STDIN);
  27   while (<>) {
  28         my @t = split;
  29         my $last = 0;
  30         for (my $i = 0; $i < $t[1];) {
  31           my $e = ($t[1] - $i) / $l < 1.1? $t[1] : $i + $l;
  32           print "$t[0]:".($i+1)."-$e\n";
  33           $i = $e;
  34         }
  35   }
  36 }
  37
  38 sub subsam {
  39   die(qq/Usage: vcfutils.pl subsam <in.vcf> [samples]\n/) if (@ARGV == 0);
  40   my ($fh, %h);
  41   my $fn = shift(@ARGV);
  42   my @col;
  43   open($fh, ($fn =~ /\.gz$/)? "gzip -dc $fn |" : $fn) || die;
  44   $h{$_} = 1 for (@ARGV);
  45   while (<$fh>) {
  46         if (/^##/) {
  47           print;
  48         } elsif (/^#/) {
  49           my @t = split;
  50           my @s = @t[0..8]; # all fixed fields + FORMAT
  51           for (9 .. $#t) {
  52                 if ($h{$t[$_]}) {
  53                   push(@s, $t[$_]);
  54                   push(@col, $_);
  55                 }
  56           }
  57           pop(@s) if (@s == 9); # no sample selected; remove the FORMAT field
  58           print join("\t", @s), "\n";
  59         } else {
  60           my @t = split;
  61           if (@col == 0) {
  62                 print join("\t", @t[0..7]), "\n";
  63           } else {
  64                 print join("\t", @t[0..8], map {$t[$_]} @col), "\n";
  65           }
  66         }
  67   }
  68   close($fh);
  69 }
  70
  71 sub listsam {
  72   die(qq/Usage: vcfutils.pl listsam <in.vcf>\n/) if (@ARGV == 0 && -t STDIN);
  73   while (<>) {
  74         if (/^#/ && !/^##/) {
  75           my @t = split;
  76           print join("\n", @t[9..$#t]), "\n";
  77           exit;
  78         }
  79   }
  80 }
  81
  82 sub fillac {
  83   die(qq/Usage: vcfutils.pl fillac <in.vcf>\n\nNote: The GT field MUST BE present and always appear as the first field.\n/) if (@ARGV == 0 && -t STDIN);
  84   while (<>) {
  85         if (/^#/) {
  86           print;
  87         } else {
  88           my @t = split;
  89           my @c = (0);
  90           my $n = 0;
  91           my $s = -1;
  92           @_ = split(":", $t[8]);
  93           for (0 .. $#_) {
  94                 if ($_[$_] eq 'GT') { $s = $_; last; }
  95           }
  96           if ($s < 0) {
  97                 print join("\t", @t), "\n";
  98                 next;
  99           }
 100           for (9 .. $#t) {
 101                 if ($t[$_] =~ /^0,0,0/) {
 102                 } elsif ($t[$_] =~ /^([^\s:]+:){$s}(\d+).(\d+)/) {
 103                   ++$c[$2]; ++$c[$3];
 104                   $n += 2;
 105                 }
 106           }
 107           my $AC = "AC=" . join("\t", @c[1..$#c]) . ";AN=$n";
 108           my $info = $t[7];
 109           $info =~ s/(;?)AC=(\d+)//;
 110           $info =~ s/(;?)AN=(\d+)//;
 111           if ($info eq '.') {
 112                 $info = $AC;
 113           } else {
 114                 $info .= ";$AC";
 115           }
 116           $t[7] = $info;
 117           print join("\t", @t), "\n";
 118         }
 119   }
 120 }
 121
 122 sub ldstats {
 123   my %opts = (t=>0.9);
 124   getopts('t:', \%opts);
 125   die("Usage: vcfutils.pl ldstats [-t $opts{t}] <in.vcf>\n") if (@ARGV == 0 && -t STDIN);
 126   my $cutoff = $opts{t};
 127   my ($last, $lastchr) = (0x7fffffff, '');
 128   my ($x, $y, $n) = (0, 0, 0);
 129   while (<>) {
 130         if (/^([^#\s]+)\s(\d+)/) {
 131           my ($chr, $pos) = ($1, $2);
 132           if (/NEIR=([\d\.]+)/) {
 133                 ++$n;
 134                 ++$y, $x += $pos - $last if ($lastchr eq $chr && $pos > $last && $1 > $cutoff);
 135           }
 136           $last = $pos; $lastchr = $chr;
 137         }
 138   }
 139   print "Number of SNP intervals in strong LD (r > $opts{t}): $y\n";
 140   print "Fraction: ", $y/$n, "\n";
 141   print "Length: $x\n";
 142 }
 143
 144 sub qstats {
 145   my %opts = (r=>'', s=>0.02, v=>undef);
 146   getopts('r:s:v', \%opts);
 147   die("Usage: vcfutils.pl qstats [-r ref.vcf] <in.vcf>\n
 148 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);
 149   my %ts = (AG=>1, GA=>1, CT=>1, TC=>1);
 150   my %h = ();
 151   my $is_vcf = defined($opts{v})? 1 : 0;
 152   if ($opts{r}) { # read the reference positions
 153         my $fh;
 154         open($fh, $opts{r}) || die;
 155         while (<$fh>) {
 156           next if (/^#/);
 157           if ($is_vcf) {
 158                 my @t = split;
 159                 $h{$t[0],$t[1]} = $t[4];
 160           } else {
 161                 $h{$1,$2} = 1 if (/^(\S+)\s+(\d+)/);
 162           }
 163         }
 164         close($fh);
 165   }
 166   my $hsize = scalar(keys %h);
 167   my @a;
 168   while (<>) {
 169         next if (/^#/);
 170         my @t = split;
 171         next if (length($t[3]) != 1 || uc($t[3]) eq 'N');
 172         $t[3] = uc($t[3]); $t[4] = uc($t[4]);
 173         my @s = split(',', $t[4]);
 174         $t[5] = 3 if ($t[5] eq '.' || $t[5] < 0);
 175         next if (length($s[0]) != 1);
 176         my $hit;
 177         if ($is_vcf) {
 178           $hit = 0;
 179           my $aa = $h{$t[0],$t[1]};
 180           if (defined($aa)) {
 181                 my @aaa = split(",", $aa);
 182                 for (@aaa) {
 183                   $hit = 1 if ($_ eq $s[0]);
 184                 }
 185           }
 186         } else {
 187           $hit = defined($h{$t[0],$t[1]})? 1 : 0;
 188         }
 189         push(@a, [$t[5], ($t[4] eq '.' || $t[4] eq $t[3])? 0 : 1, $ts{$t[3].$s[0]}? 1 : 0, $hit]);
 190   }
 191   push(@a, [-1, 0, 0, 0]); # end marker
 192   die("[qstats] No SNP data!\n") if (@a == 0);
 193   @a = sort {$b->[0]<=>$a->[0]} @a;
 194   my $next = $opts{s};
 195   my $last = $a[0];
 196   my @c = (0, 0, 0, 0);
 197   my @lc;
 198   $lc[1] = $lc[2] = 0;
 199   for my $p (@a) {
 200         if ($p->[0] == -1 || ($p->[0] != $last && $c[0]/@a > $next)) {
 201           my @x;
 202           $x[0] = sprintf("%.4f", $c[1]-$c[2]? $c[2] / ($c[1] - $c[2]) : 100);
 203           $x[1] = sprintf("%.4f", $hsize? $c[3] / $hsize : 0);
 204           $x[2] = sprintf("%.4f", $c[3] / $c[1]);
 205           my $a = $c[1] - $lc[1];
 206           my $b = $c[2] - $lc[2];
 207           $x[3] = sprintf("%.4f", $a-$b? $b / ($a-$b) : 100);
 208           print join("\t", $last, @c, @x), "\n";
 209           $next = $c[0]/@a + $opts{s};
 210           $lc[1] = $c[1]; $lc[2] = $c[2];
 211         }
 212         ++$c[0]; $c[1] += $p->[1]; $c[2] += $p->[2]; $c[3] += $p->[3];
 213         $last = $p->[0];
 214   }
 215 }
 216
 217 sub varFilter {
 218   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, G=>0, S=>1000);
 219   getopts('pd:D:W:Q:w:a:1:2:3:4:G:S:', \%opts);
 220   die(qq/
 221 Usage:   vcfutils.pl varFilter [options] <in.vcf>
 222
 223 Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
 224          -d INT    minimum read depth [$opts{d}]
 225          -D INT    maximum read depth [$opts{D}]
 226          -a INT    minimum number of alternate bases [$opts{a}]
 227          -w INT    SNP within INT bp around a gap to be filtered [$opts{w}]
 228          -W INT    window size for filtering adjacent gaps [$opts{W}]
 229          -1 FLOAT  min P-value for strand bias (given PV4) [$opts{1}]
 230          -2 FLOAT  min P-value for baseQ bias [$opts{2}]
 231          -3 FLOAT  min P-value for mapQ bias [$opts{3}]
 232          -4 FLOAT  min P-value for end distance bias [$opts{4}]
 233          -p        print filtered variants
 234
 235 Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools.
 236 \n/) if (@ARGV == 0 && -t STDIN);
 237
 238   # calculate the window size
 239   my ($ol, $ow) = ($opts{W}, $opts{w});
 240   my $max_dist = $ol > $ow? $ol : $ow;
 241   # the core loop
 242   my @staging; # (indel_filtering_score, flt_tag, indel_span; chr, pos, ...)
 243   while (<>) {
 244         my @t = split;
 245     if (/^#/) {
 246           print; next;
 247         }
 248         next if ($t[4] eq '.'); # skip non-var sites
 249         # check if the site is a SNP
 250         my $type = 1; # SNP
 251         if (length($t[3]) > 1) {
 252           $type = 2; # MNP
 253           my @s = split(',', $t[4]);
 254           for (@s) {
 255                 $type = 3 if (length != length($t[3]));
 256           }
 257         } else {
 258           my @s = split(',', $t[4]);
 259           for (@s) {
 260                 $type = 3 if (length > 1);
 261           }
 262         }
 263         # clear the out-of-range elements
 264         while (@staging) {
 265       # Still on the same chromosome and the first element's window still affects this position?
 266           last if ($staging[0][3] eq $t[0] && $staging[0][4] + $staging[0][2] + $max_dist >= $t[1]);
 267           varFilter_aux(shift(@staging), $opts{p}); # calling a function is a bit slower, not much
 268         }
 269         my $flt = 0;
 270         # parse annotations
 271         my ($dp, $mq, $dp_alt) = (-1, -1, -1);
 272         if ($t[7] =~ /DP4=(\d+),(\d+),(\d+),(\d+)/i) {
 273           $dp = $1 + $2 + $3 + $4;
 274           $dp_alt = $3 + $4;
 275         }
 276         if ($t[7] =~ /DP=(\d+)/i) {
 277           $dp = $1;
 278         }
 279         $mq = $1 if ($t[7] =~ /MQ=(\d+)/i);
 280         # the depth and mapQ filter
 281         if ($dp >= 0) {
 282           if ($dp < $opts{d}) {
 283                 $flt = 2;
 284           } elsif ($dp > $opts{D}) {
 285                 $flt = 3;
 286           }
 287         }
 288         $flt = 4 if ($dp_alt >= 0 && $dp_alt < $opts{a});
 289         $flt = 1 if ($flt == 0 && $mq >= 0 && $mq < $opts{Q});
 290         $flt = 7 if ($flt == 0 && /PV4=([^,]+),([^,]+),([^,]+),([^,;\t]+)/
 291                                  && ($1<$opts{1} || $2<$opts{2} || $3<$opts{3} || $4<$opts{4}));
 292         $flt = 8 if ($flt == 0 && ((/MXGQ=(\d+)/ && $1 >= $opts{G}) || (/MXSP=(\d+)/ && $1 < $opts{S})));
 293
 294         my $score = $t[5] * 100 + $dp_alt;
 295         my $rlen = length($t[3]) - 1; # $indel_score<0 for SNPs
 296         if ($flt == 0) {
 297           if ($type == 3) { # an indel
 298                 # filtering SNPs and MNPs
 299                 for my $x (@staging) {
 300                   next if (($x->[0]&3) == 3 || $x->[1] || $x->[4] + $x->[2] + $ow < $t[1]);
 301                   $x->[1] = 5;
 302                 }
 303                 # check the staging list for indel filtering
 304                 for my $x (@staging) {
 305                   next if (($x->[0]&3) != 3 || $x->[1] || $x->[4] + $x->[2] + $ol < $t[1]);
 306                   if ($x->[0]>>2 < $score) {
 307                         $x->[1] = 6;
 308                   } else {
 309                         $flt = 6; last;
 310                   }
 311                 }
 312           } else { # SNP or MNP
 313                 for my $x (@staging) {
 314                   next if (($x->[0]&3) != 3 || $x->[4] + $x->[2] + $ow < $t[1]);
 315                   $flt = 5;
 316                   last;
 317                 }
 318                 # check MNP
 319                 for my $x (@staging) {
 320                   next if (($x->[0]&3) == 3 || $x->[4] + $x->[2] < $t[1]);
 321                   if ($x->[0]>>2 < $score) {
 322                         $x->[1] = 8;
 323                   } else {
 324                         $flt = 8; last;
 325                   }
 326                 }
 327           }
 328         }
 329         push(@staging, [$score<<2|$type, $flt, $rlen, @t]);
 330   }
 331   # output the last few elements in the staging list
 332   while (@staging) {
 333         varFilter_aux(shift @staging, $opts{p});
 334   }
 335 }
 336
 337 sub varFilter_aux {
 338   my ($first, $is_print) = @_;
 339   if ($first->[1] == 0) {
 340         print join("\t", @$first[3 .. @$first-1]), "\n";
 341   } elsif ($is_print) {
 342         print STDERR join("\t", substr("UQdDaGgPMS", $first->[1], 1), @$first[3 .. @$first-1]), "\n";
 343   }
 344 }
 345
 346 sub gapstats {
 347   my (@c0, @c1);
 348   $c0[$_] = $c1[$_] = 0 for (0 .. 10000);
 349   while (<>) {
 350         next if (/^#/);
 351         my @t = split;
 352         next if (length($t[3]) == 1 && $t[4] =~ /^[A-Za-z](,[A-Za-z])*$/); # not an indel
 353         my @s = split(',', $t[4]);
 354         for my $x (@s) {
 355           my $l = length($x) - length($t[3]) + 5000;
 356           if ($x =~ /^-/) {
 357                 $l = -(length($x) - 1) + 5000;
 358           } elsif ($x =~ /^\+/) {
 359                 $l = length($x) - 1 + 5000;
 360           }
 361           $c0[$l] += 1 / @s;
 362         }
 363   }
 364   for (my $i = 0; $i < 10000; ++$i) {
 365         next if ($c0[$i] == 0);
 366         $c1[0] += $c0[$i];
 367         $c1[1] += $c0[$i] if (($i-5000)%3 == 0);
 368         printf("C\t%d\t%.2f\n", ($i-5000), $c0[$i]);
 369   }
 370   printf("3\t%d\t%d\t%.3f\n", $c1[0], $c1[1], $c1[1]/$c1[0]);
 371 }
 372
 373 sub ucscsnp2vcf {
 374   die("Usage: vcfutils.pl <in.ucsc.snp>\n") if (@ARGV == 0 && -t STDIN);
 375   print "##fileformat=VCFv4.0\n";
 376   print join("\t", "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO"), "\n";
 377   while (<>) {
 378         my @t = split("\t");
 379         my $indel = ($t[9] =~ /^[ACGT](\/[ACGT])+$/)? 0 : 1;
 380         my $pos = $t[2] + 1;
 381         my @alt;
 382         push(@alt, $t[7]);
 383         if ($t[6] eq '-') {
 384           $t[9] = reverse($t[9]);
 385           $t[9] =~ tr/ACGTRYMKWSNacgtrymkwsn/TGCAYRKMWSNtgcayrkmwsn/;
 386         }
 387         my @a = split("/", $t[9]);
 388         for (@a) {
 389           push(@alt, $_) if ($_ ne $alt[0]);
 390         }
 391         if ($indel) {
 392           --$pos;
 393           for (0 .. $#alt) {
 394                 $alt[$_] =~ tr/-//d;
 395                 $alt[$_] = "N$alt[$_]";
 396           }
 397         }
 398         my $ref = shift(@alt);
 399         my $af = $t[13] > 0? ";AF=$t[13]" : '';
 400         my $valid = ($t[12] eq 'unknown')? '' : ";valid=$t[12]";
 401         my $info = "molType=$t[10];class=$t[11]$valid$af";
 402         print join("\t", $t[1], $pos, $t[4], $ref, join(",", @alt), 0, '.', $info), "\n";
 403   }
 404 }
 405
 406 sub hapmap2vcf {
 407   die("Usage: vcfutils.pl <in.ucsc.snp> <in.hapmap>\n") if (@ARGV == 0);
 408   my $fn = shift(@ARGV);
 409   # parse UCSC SNP
 410   warn("Parsing UCSC SNPs...\n");
 411   my ($fh, %map);
 412   open($fh, ($fn =~ /\.gz$/)? "gzip -dc $fn |" : $fn) || die;
 413   while (<$fh>) {
 414         my @t = split;
 415         next if ($t[3] - $t[2] != 1); # not SNP
 416         @{$map{$t[4]}} = @t[1,3,7];
 417   }
 418   close($fh);
 419   # write VCF
 420   warn("Writing VCF...\n");
 421   print "##fileformat=VCFv4.0\n";
 422   while (<>) {
 423         my @t = split;
 424         if ($t[0] eq 'rs#') { # the first line
 425           print join("\t", "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT", @t[11..$#t]), "\n";
 426         } else {
 427           next unless ($map{$t[0]});
 428           next if (length($t[1]) != 3); # skip non-SNPs
 429           my $a = \@{$map{$t[0]}};
 430           my $ref = $a->[2];
 431           my @u = split('/', $t[1]);
 432           if ($u[1] eq $ref) {
 433                 $u[1] = $u[0]; $u[0] = $ref;
 434           } elsif ($u[0] ne $ref) { next; }
 435           my $alt = $u[1];
 436           my %w;
 437           $w{$u[0]} = 0; $w{$u[1]} = 1;
 438           my @s = (@$a[0,1], $t[0], $ref, $alt, 0, '.', '.', 'GT');
 439           my $is_tri = 0;
 440           for (@t[11..$#t]) {
 441                 if ($_ eq 'NN') {
 442                   push(@s, './.');
 443                 } else {
 444                   my @a = ($w{substr($_,0,1)}, $w{substr($_,1,1)});
 445                   if (!defined($a[0]) || !defined($a[1])) {
 446                         $is_tri = 1;
 447                         last;
 448                   }
 449                   push(@s, "$a[0]/$a[1]");
 450                 }
 451           }
 452           next if ($is_tri);
 453           print join("\t", @s), "\n";
 454         }
 455   }
 456 }
 457
 458 sub vcf2fq {
 459   my %opts = (d=>3, D=>100000, Q=>10, l=>5);
 460   getopts('d:D:Q:l:', \%opts);
 461   die(qq/
 462 Usage:   vcfutils.pl vcf2fq [options] <all-site.vcf>
 463
 464 Options: -d INT    minimum depth          [$opts{d}]
 465          -D INT    maximum depth          [$opts{D}]
 466          -Q INT    min RMS mapQ           [$opts{Q}]
 467          -l INT    INDEL filtering window [$opts{l}]
 468 \n/) if (@ARGV == 0 && -t STDIN);
 469
 470   my ($last_chr, $seq, $qual, $last_pos, @gaps);
 471   my $_Q = $opts{Q};
 472   my $_d = $opts{d};
 473   my $_D = $opts{D};
 474
 475   my %het = (AC=>'M', AG=>'R', AT=>'W', CA=>'M', CG=>'S', CT=>'Y',
 476                          GA=>'R', GC=>'S', GT=>'K', TA=>'W', TC=>'Y', TG=>'K');
 477
 478   $last_chr = '';
 479   while (<>) {
 480         next if (/^#/);
 481         my @t = split;
 482         if ($last_chr ne $t[0]) {
 483           &v2q_post_process($last_chr, \$seq, \$qual, \@gaps, $opts{l}) if ($last_chr);
 484           ($last_chr, $last_pos) = ($t[0], 0);
 485           $seq = $qual = '';
 486           @gaps = ();
 487         }
 488         die("[vcf2fq] unsorted input\n") if ($t[1] - $last_pos < 0);
 489         if ($t[1] - $last_pos > 1) {
 490           $seq .= 'n' x ($t[1] - $last_pos - 1);
 491           $qual .= '!' x ($t[1] - $last_pos - 1);
 492         }
 493         if (length($t[3]) == 1 && $t[7] !~ /INDEL/ && $t[4] =~ /^([A-Za-z.])(,[A-Za-z])*$/) { # a SNP or reference
 494           my ($ref, $alt) = ($t[3], $1);
 495           my ($b, $q);
 496           $q = $1 if ($t[7] =~ /FQ=(-?[\d\.]+)/);
 497           if ($q < 0) {
 498                 $_ = $1 if ($t[7] =~ /AF1=([\d\.]+)/);
 499                 $b = ($_ < .5 || $alt eq '.')? $ref : $alt;
 500                 $q = -$q;
 501           } else {
 502                 $b = $het{"$ref$alt"};
 503                 $b ||= 'N';
 504           }
 505           $b = lc($b);
 506           $b = uc($b) if (($t[7] =~ /MQ=(\d+)/ && $1 >= $_Q) && ($t[7] =~ /DP=(\d+)/ && $1 >= $_d && $1 <= $_D));
 507           $q = int($q + 33 + .499);
 508           $q = chr($q <= 126? $q : 126);
 509           $seq .= $b;
 510           $qual .= $q;
 511         } else { # an INDEL
 512           push(@gaps, [$t[1], length($t[3])]);
 513         }
 514         $last_pos = $t[1];
 515   }
 516   &v2q_post_process($last_chr, \$seq, \$qual, \@gaps, $opts{l});
 517 }
 518
 519 sub v2q_post_process {
 520   my ($chr, $seq, $qual, $gaps, $l) = @_;
 521   for my $g (@$gaps) {
 522         my $beg = $g->[0] > $l? $g->[0] - $l : 0;
 523         my $end = $g->[0] + $g->[1] + $l;
 524         $end = length($$seq) if ($end > length($$seq));
 525         substr($$seq, $beg, $end - $beg) = lc(substr($$seq, $beg, $end - $beg));
 526   }
 527   print "\@$chr\n"; &v2q_print_str($seq);
 528   print "+\n"; &v2q_print_str($qual);
 529 }
 530
 531 sub v2q_print_str {
 532   my ($s) = @_;
 533   my $l = length($$s);
 534   for (my $i = 0; $i < $l; $i += 60) {
 535         print substr($$s, $i, 60), "\n";
 536   }
 537 }
 538
 539 sub usage {
 540   die(qq/
 541 Usage:   vcfutils.pl <command> [<arguments>]\n
 542 Command: subsam       get a subset of samples
 543          listsam      list the samples
 544          fillac       fill the allele count field
 545          qstats       SNP stats stratified by QUAL
 546
 547          hapmap2vcf   convert the hapmap format to VCF
 548          ucscsnp2vcf  convert UCSC SNP SQL dump to VCF
 549
 550          varFilter    filtering short variants (*)
 551          vcf2fq       VCF->fastq (**)
 552
 553 Notes: Commands with description endting with (*) may need bcftools
 554        specific annotations.
 555 \n/);
 556 }