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);
 219   getopts('pd:D:W:Q:w:a:1:2:3:4:', \%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
 293         my $score = $t[5] * 100 + $dp_alt;
 294         my $rlen = length($t[3]) - 1; # $indel_score<0 for SNPs
 295         if ($flt == 0) {
 296           if ($type == 3) { # an indel
 297                 # filtering SNPs and MNPs
 298                 for my $x (@staging) {
 299                   next if (($x->[0]&3) == 3 || $x->[1] || $x->[4] + $x->[2] + $ow < $t[1]);
 300                   $x->[1] = 5;
 301                 }
 302                 # check the staging list for indel filtering
 303                 for my $x (@staging) {
 304                   next if (($x->[0]&3) != 3 || $x->[1] || $x->[4] + $x->[2] + $ol < $t[1]);
 305                   if ($x->[0]>>2 < $score) {
 306                         $x->[1] = 6;
 307                   } else {
 308                         $flt = 6; last;
 309                   }
 310                 }
 311           } else { # SNP or MNP
 312                 for my $x (@staging) {
 313                   next if (($x->[0]&3) != 3 || $x->[4] + $x->[2] + $ow < $t[1]);
 314                   $flt = 5;
 315                   last;
 316                 }
 317                 # check MNP
 318                 for my $x (@staging) {
 319                   next if (($x->[0]&3) == 3 || $x->[4] + $x->[2] < $t[1]);
 320                   if ($x->[0]>>2 < $score) {
 321                         $x->[1] = 8;
 322                   } else {
 323                         $flt = 8; last;
 324                   }
 325                 }
 326           }
 327         }
 328         push(@staging, [$score<<2|$type, $flt, $rlen, @t]);
 329   }
 330   # output the last few elements in the staging list
 331   while (@staging) {
 332         varFilter_aux(shift @staging, $opts{p});
 333   }
 334 }
 335
 336 sub varFilter_aux {
 337   my ($first, $is_print) = @_;
 338   if ($first->[1] == 0) {
 339         print join("\t", @$first[3 .. @$first-1]), "\n";
 340   } elsif ($is_print) {
 341         print STDERR join("\t", substr("UQdDaGgPM", $first->[1], 1), @$first[3 .. @$first-1]), "\n";
 342   }
 343 }
 344
 345 sub gapstats {
 346   my (@c0, @c1);
 347   $c0[$_] = $c1[$_] = 0 for (0 .. 10000);
 348   while (<>) {
 349         next if (/^#/);
 350         my @t = split;
 351         next if (length($t[3]) == 1 && $t[4] =~ /^[A-Za-z](,[A-Za-z])*$/); # not an indel
 352         my @s = split(',', $t[4]);
 353         for my $x (@s) {
 354           my $l = length($x) - length($t[3]) + 5000;
 355           if ($x =~ /^-/) {
 356                 $l = -(length($x) - 1) + 5000;
 357           } elsif ($x =~ /^\+/) {
 358                 $l = length($x) - 1 + 5000;
 359           }
 360           $c0[$l] += 1 / @s;
 361         }
 362   }
 363   for (my $i = 0; $i < 10000; ++$i) {
 364         next if ($c0[$i] == 0);
 365         $c1[0] += $c0[$i];
 366         $c1[1] += $c0[$i] if (($i-5000)%3 == 0);
 367         printf("C\t%d\t%.2f\n", ($i-5000), $c0[$i]);
 368   }
 369   printf("3\t%d\t%d\t%.3f\n", $c1[0], $c1[1], $c1[1]/$c1[0]);
 370 }
 371
 372 sub ucscsnp2vcf {
 373   die("Usage: vcfutils.pl <in.ucsc.snp>\n") if (@ARGV == 0 && -t STDIN);
 374   print "##fileformat=VCFv4.0\n";
 375   print join("\t", "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO"), "\n";
 376   while (<>) {
 377         my @t = split("\t");
 378         my $indel = ($t[9] =~ /^[ACGT](\/[ACGT])+$/)? 0 : 1;
 379         my $pos = $t[2] + 1;
 380         my @alt;
 381         push(@alt, $t[7]);
 382         if ($t[6] eq '-') {
 383           $t[9] = reverse($t[9]);
 384           $t[9] =~ tr/ACGTRYMKWSNacgtrymkwsn/TGCAYRKMWSNtgcayrkmwsn/;
 385         }
 386         my @a = split("/", $t[9]);
 387         for (@a) {
 388           push(@alt, $_) if ($_ ne $alt[0]);
 389         }
 390         if ($indel) {
 391           --$pos;
 392           for (0 .. $#alt) {
 393                 $alt[$_] =~ tr/-//d;
 394                 $alt[$_] = "N$alt[$_]";
 395           }
 396         }
 397         my $ref = shift(@alt);
 398         my $af = $t[13] > 0? ";AF=$t[13]" : '';
 399         my $valid = ($t[12] eq 'unknown')? '' : ";valid=$t[12]";
 400         my $info = "molType=$t[10];class=$t[11]$valid$af";
 401         print join("\t", $t[1], $pos, $t[4], $ref, join(",", @alt), 0, '.', $info), "\n";
 402   }
 403 }
 404
 405 sub hapmap2vcf {
 406   die("Usage: vcfutils.pl <in.ucsc.snp> <in.hapmap>\n") if (@ARGV == 0);
 407   my $fn = shift(@ARGV);
 408   # parse UCSC SNP
 409   warn("Parsing UCSC SNPs...\n");
 410   my ($fh, %map);
 411   open($fh, ($fn =~ /\.gz$/)? "gzip -dc $fn |" : $fn) || die;
 412   while (<$fh>) {
 413         my @t = split;
 414         next if ($t[3] - $t[2] != 1); # not SNP
 415         @{$map{$t[4]}} = @t[1,3,7];
 416   }
 417   close($fh);
 418   # write VCF
 419   warn("Writing VCF...\n");
 420   print "##fileformat=VCFv4.0\n";
 421   while (<>) {
 422         my @t = split;
 423         if ($t[0] eq 'rs#') { # the first line
 424           print join("\t", "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT", @t[11..$#t]), "\n";
 425         } else {
 426           next unless ($map{$t[0]});
 427           next if (length($t[1]) != 3); # skip non-SNPs
 428           my $a = \@{$map{$t[0]}};
 429           my $ref = $a->[2];
 430           my @u = split('/', $t[1]);
 431           if ($u[1] eq $ref) {
 432                 $u[1] = $u[0]; $u[0] = $ref;
 433           } elsif ($u[0] ne $ref) { next; }
 434           my $alt = $u[1];
 435           my %w;
 436           $w{$u[0]} = 0; $w{$u[1]} = 1;
 437           my @s = (@$a[0,1], $t[0], $ref, $alt, 0, '.', '.', 'GT');
 438           my $is_tri = 0;
 439           for (@t[11..$#t]) {
 440                 if ($_ eq 'NN') {
 441                   push(@s, './.');
 442                 } else {
 443                   my @a = ($w{substr($_,0,1)}, $w{substr($_,1,1)});
 444                   if (!defined($a[0]) || !defined($a[1])) {
 445                         $is_tri = 1;
 446                         last;
 447                   }
 448                   push(@s, "$a[0]/$a[1]");
 449                 }
 450           }
 451           next if ($is_tri);
 452           print join("\t", @s), "\n";
 453         }
 454   }
 455 }
 456
 457 sub vcf2fq {
 458   my %opts = (d=>3, D=>100000, Q=>10, l=>5);
 459   getopts('d:D:Q:l:', \%opts);
 460   die(qq/
 461 Usage:   vcfutils.pl vcf2fq [options] <all-site.vcf>
 462
 463 Options: -d INT    minimum depth          [$opts{d}]
 464          -D INT    maximum depth          [$opts{D}]
 465          -Q INT    min RMS mapQ           [$opts{Q}]
 466          -l INT    INDEL filtering window [$opts{l}]
 467 \n/) if (@ARGV == 0 && -t STDIN);
 468
 469   my ($last_chr, $seq, $qual, $last_pos, @gaps);
 470   my $_Q = $opts{Q};
 471   my $_d = $opts{d};
 472   my $_D = $opts{D};
 473
 474   my %het = (AC=>'M', AG=>'R', AT=>'W', CA=>'M', CG=>'S', CT=>'Y',
 475                          GA=>'R', GC=>'S', GT=>'K', TA=>'W', TC=>'Y', TG=>'K');
 476
 477   $last_chr = '';
 478   while (<>) {
 479         next if (/^#/);
 480         my @t = split;
 481         if ($last_chr ne $t[0]) {
 482           &v2q_post_process($last_chr, \$seq, \$qual, \@gaps, $opts{l}) if ($last_chr);
 483           ($last_chr, $last_pos) = ($t[0], 0);
 484           $seq = $qual = '';
 485           @gaps = ();
 486         }
 487         if ($t[1] - $last_pos != 1) {
 488           $seq .= 'n' x ($t[1] - $last_pos - 1);
 489           $qual .= '!' x ($t[1] - $last_pos - 1);
 490         }
 491         if (length($t[3]) == 1 && $t[4] =~ /^([A-Za-z.])(,[A-Za-z])*$/) { # a SNP or reference
 492           my ($ref, $alt) = ($t[3], $1);
 493           my ($b, $q);
 494           $q = $1 if ($t[7] =~ /FQ=(-?[\d\.]+)/);
 495           if ($q < 0) {
 496                 $_ = $1 if ($t[7] =~ /AF1=([\d\.]+)/);
 497                 $b = ($_ < .5 || $alt eq '.')? $ref : $alt;
 498                 $q = -$q;
 499           } else {
 500                 $b = $het{"$ref$alt"};
 501                 $b ||= 'N';
 502           }
 503           $b = lc($b);
 504           $b = uc($b) if (($t[7] =~ /MQ=(\d+)/ && $1 >= $_Q) && ($t[7] =~ /DP=(\d+)/ && $1 >= $_d && $1 <= $_D));
 505           $q = int($q + 33 + .499);
 506           $q = chr($q <= 126? $q : 126);
 507           $seq .= $b;
 508           $qual .= $q;
 509         } else { # an INDEL
 510           push(@gaps, [$t[1], length($t[3])]);
 511         }
 512         $last_pos = $t[1];
 513   }
 514   &v2q_post_process($last_chr, \$seq, \$qual, \@gaps, $opts{l});
 515 }
 516
 517 sub v2q_post_process {
 518   my ($chr, $seq, $qual, $gaps, $l) = @_;
 519   for my $g (@$gaps) {
 520         my $beg = $g->[0] > $l? $g->[0] - $l : 0;
 521         my $end = $g->[0] + $g->[1] + $l;
 522         $end = length($$seq) if ($end > length($$seq));
 523         substr($$seq, $beg, $end - $beg) = lc(substr($$seq, $beg, $end - $beg));
 524   }
 525   print "\@$chr\n"; &v2q_print_str($seq);
 526   print "+\n"; &v2q_print_str($qual);
 527 }
 528
 529 sub v2q_print_str {
 530   my ($s) = @_;
 531   my $l = length($$s);
 532   for (my $i = 0; $i < $l; $i += 60) {
 533         print substr($$s, $i, 60), "\n";
 534   }
 535 }
 536
 537 sub usage {
 538   die(qq/
 539 Usage:   vcfutils.pl <command> [<arguments>]\n
 540 Command: subsam       get a subset of samples
 541          listsam      list the samples
 542          fillac       fill the allele count field
 543          qstats       SNP stats stratified by QUAL
 544
 545          hapmap2vcf   convert the hapmap format to VCF
 546          ucscsnp2vcf  convert UCSC SNP SQL dump to VCF
 547
 548          varFilter    filtering short variants (*)
 549          vcf2fq       VCF->fastq (*)
 550
 551 Notes: Commands with description endting with (*) may need bcftools
 552        specific annotations.
 553 \n/);
 554 }