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   my $version = '0.1.0';
  14   &usage if (@ARGV < 1);
  15   my $command = shift(@ARGV);
  16   my %func = (subsam=>\&subsam, listsam=>\&listsam, fillac=>\&fillac, qstats=>\&qstats, varFilter=>\&varFilter,
  17                           hapmap2vcf=>\&hapmap2vcf, ucscsnp2vcf=>\&ucscsnp2vcf, filter4vcf=>\&filter4vcf);
  18   die("Unknown command \"$command\".\n") if (!defined($func{$command}));
  19   &{$func{$command}};
  20 }
  21
  22 sub subsam {
  23   die(qq/Usage: vcfutils.pl subsam <in.vcf> [samples]\n/) if (@ARGV == 0);
  24   my ($fh, %h);
  25   my $fn = shift(@ARGV);
  26   my @col;
  27   open($fh, ($fn =~ /\.gz$/)? "gzip -dc $fn |" : $fn) || die;
  28   $h{$_} = 1 for (@ARGV);
  29   while (<$fh>) {
  30         if (/^##/) {
  31           print;
  32         } elsif (/^#/) {
  33           my @t = split;
  34           my @s = @t[0..8]; # all fixed fields + FORMAT
  35           for (9 .. $#t) {
  36                 if ($h{$t[$_]}) {
  37                   push(@s, $t[$_]);
  38                   push(@col, $_);
  39                 }
  40           }
  41           pop(@s) if (@s == 9); # no sample selected; remove the FORMAT field
  42           print join("\t", @s), "\n";
  43         } else {
  44           my @t = split;
  45           if (@col == 0) {
  46                 print join("\t", @t[0..7]), "\n";
  47           } else {
  48                 print join("\t", @t[0..8], map {$t[$_]} @col), "\n";
  49           }
  50         }
  51   }
  52   close($fh);
  53 }
  54
  55 sub listsam {
  56   die(qq/Usage: vcfutils.pl listsam <in.vcf>\n/) if (@ARGV == 0 && -t STDIN);
  57   while (<>) {
  58         if (/^#/ && !/^##/) {
  59           my @t = split;
  60           print join("\n", @t[9..$#t]), "\n";
  61           exit;
  62         }
  63   }
  64 }
  65
  66 sub fillac {
  67   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);
  68   while (<>) {
  69         if (/^#/) {
  70           print;
  71         } else {
  72           my @t = split;
  73           my @c = (0);
  74           my $n = 0;
  75           my $s = -1;
  76           @_ = split(":", $t[8]);
  77           for (0 .. $#_) {
  78                 if ($_[$_] eq 'GT') { $s = $_; last; }
  79           }
  80           if ($s < 0) {
  81                 print join("\t", @t), "\n";
  82                 next;
  83           }
  84           for (9 .. $#t) {
  85                 if ($t[$_] =~ /^0,0,0/) {
  86                 } elsif ($t[$_] =~ /^([^\s:]+:){$s}(\d+).(\d+)/) {
  87                   ++$c[$2]; ++$c[$3];
  88                   $n += 2;
  89                 }
  90           }
  91           my $AC = "AC=" . join("\t", @c[1..$#c]) . ";AN=$n";
  92           my $info = $t[7];
  93           $info =~ s/(;?)AC=(\d+)//;
  94           $info =~ s/(;?)AN=(\d+)//;
  95           if ($info eq '.') {
  96                 $info = $AC;
  97           } else {
  98                 $info .= ";$AC";
  99           }
 100           $t[7] = $info;
 101           print join("\t", @t), "\n";
 102         }
 103   }
 104 }
 105
 106 sub qstats {
 107   my %opts = (r=>'', s=>0.02, v=>undef);
 108   getopts('r:s:v', \%opts);
 109   die("Usage: vcfutils.pl qstats [-r ref.vcf] <in.vcf>\n
 110 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);
 111   my %ts = (AG=>1, GA=>1, CT=>1, TC=>1);
 112   my %h = ();
 113   my $is_vcf = defined($opts{v})? 1 : 0;
 114   if ($opts{r}) { # read the reference positions
 115         my $fh;
 116         open($fh, $opts{r}) || die;
 117         while (<$fh>) {
 118           next if (/^#/);
 119           if ($is_vcf) {
 120                 my @t = split;
 121                 $h{$t[0],$t[1]} = $t[4];
 122           } else {
 123                 $h{$1,$2} = 1 if (/^(\S+)\s+(\d+)/);
 124           }
 125         }
 126         close($fh);
 127   }
 128   my $hsize = scalar(keys %h);
 129   my @a;
 130   while (<>) {
 131         next if (/^#/);
 132         my @t = split;
 133         next if (length($t[3]) != 1 || uc($t[3]) eq 'N');
 134         $t[3] = uc($t[3]); $t[4] = uc($t[4]);
 135         my @s = split(',', $t[4]);
 136         $t[5] = 3 if ($t[5] < 0);
 137         next if (length($s[0]) != 1);
 138         my $hit;
 139         if ($is_vcf) {
 140           $hit = 0;
 141           my $aa = $h{$t[0],$t[1]};
 142           if (defined($aa)) {
 143                 my @aaa = split(",", $aa);
 144                 for (@aaa) {
 145                   $hit = 1 if ($_ eq $s[0]);
 146                 }
 147           }
 148         } else {
 149           $hit = defined($h{$t[0],$t[1]})? 1 : 0;
 150         }
 151         push(@a, [$t[5], ($t[4] eq '.' || $t[4] eq $t[3])? 0 : 1, $ts{$t[3].$s[0]}? 1 : 0, $hit]);
 152   }
 153   push(@a, [-1, 0, 0, 0]); # end marker
 154   die("[qstats] No SNP data!\n") if (@a == 0);
 155   @a = sort {$b->[0]<=>$a->[0]} @a;
 156   my $next = $opts{s};
 157   my $last = $a[0];
 158   my @c = (0, 0, 0, 0);
 159   my @lc;
 160   $lc[1] = $lc[2] = 0;
 161   for my $p (@a) {
 162         if ($p->[0] == -1 || ($p->[0] != $last && $c[0]/@a > $next)) {
 163           my @x;
 164           $x[0] = sprintf("%.4f", $c[1]-$c[2]? $c[2] / ($c[1] - $c[2]) : 100);
 165           $x[1] = sprintf("%.4f", $hsize? $c[3] / $hsize : 0);
 166           $x[2] = sprintf("%.4f", $c[3] / $c[1]);
 167           my $a = $c[1] - $lc[1];
 168           my $b = $c[2] - $lc[2];
 169           $x[3] = sprintf("%.4f", $a-$b? $b / ($a-$b) : 100);
 170           print join("\t", $last, @c, @x), "\n";
 171           $next = $c[0]/@a + $opts{s};
 172           $lc[1] = $c[1]; $lc[2] = $c[2];
 173         }
 174         ++$c[0]; $c[1] += $p->[1]; $c[2] += $p->[2]; $c[3] += $p->[3];
 175         $last = $p->[0];
 176   }
 177 }
 178
 179 sub varFilter {
 180   my %opts = (d=>1, D=>10000, l=>30, Q=>25, q=>10, G=>25, s=>100, w=>10, W=>10, N=>2, p=>undef, F=>.001);
 181   getopts('pq:d:D:l:Q:w:W:N:G:F:', \%opts);
 182   die(qq/
 183 Usage:   vcfutils.pl varFilter [options] <in.vcf>
 184
 185 Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
 186          -q INT    minimum RMS mapping quality for gaps [$opts{q}]
 187          -d INT    minimum read depth [$opts{d}]
 188          -D INT    maximum read depth [$opts{D}]
 189
 190          -G INT    min indel score for nearby SNP filtering [$opts{G}]
 191          -w INT    SNP within INT bp around a gap to be filtered [$opts{w}]
 192
 193          -W INT    window size for filtering dense SNPs [$opts{W}]
 194          -N INT    max number of SNPs in a window [$opts{N}]
 195
 196          -l INT    window size for filtering adjacent gaps [$opts{l}]
 197
 198          -p        print filtered variants
 199 \n/) if (@ARGV == 0 && -t STDIN);
 200
 201   # calculate the window size
 202   my ($ol, $ow, $oW) = ($opts{l}, $opts{w}, $opts{W});
 203   my $max_dist = $ol > $ow? $ol : $ow;
 204   $max_dist = $oW if ($max_dist < $oW);
 205   # the core loop
 206   my @staging; # (indel_filtering_score, flt_tag)
 207   while (<>) {
 208         my @t = split;
 209         next if (/^#/);
 210         next if ($t[4] eq '.'); # skip non-var sites
 211         my $is_snp = 1;
 212         if (length($t[3]) > 1) {
 213           $is_snp = 0;
 214         } else {
 215           my @s = split(',', $t[4]);
 216           for (@s) {
 217                 $is_snp = 0 if (length > 1);
 218           }
 219         }
 220         # clear the out-of-range elements
 221         while (@staging) {
 222       # Still on the same chromosome and the first element's window still affects this position?
 223           last if ($staging[0][3] eq $t[0] && $staging[0][4] + $staging[0][2] + $max_dist >= $t[1]);
 224           varFilter_aux(shift(@staging), $opts{p}); # calling a function is a bit slower, not much
 225         }
 226         my ($flt, $score) = (0, -1);
 227
 228         # collect key annotations
 229         my ($dp, $mq, $af) = (-1, -1, 1);
 230         if ($t[7] =~ /DP=(\d+)/i) {
 231           $dp = $1;
 232         } elsif ($t[7] =~ /DP4=(\d+),(\d+),(\d+),(\d+)/i) {
 233           $dp = $1 + $2 + $3 + $4;
 234         }
 235         if ($t[7] =~ /MQ=(\d+)/i) {
 236           $mq = $1;
 237         }
 238         if ($t[7] =~ /AF=([^\s;=]+)/i) {
 239           $af = $1;
 240         } elsif ($t[7] =~ /AF1=([^\s;=]+)/i) {
 241           $af = $1;
 242         }
 243         # the depth filter
 244         if ($dp >= 0) {
 245           if ($dp < $opts{d}) {
 246                 $flt = 2;
 247           } elsif ($dp > $opts{D}) {
 248                 $flt = 3;
 249           }
 250         }
 251
 252         # site dependent filters
 253         my $dlen = 0;
 254         if ($flt == 0) {
 255           if (!$is_snp) { # an indel
 256         # If deletion, remember the length of the deletion
 257                 $dlen = length($t[3]) - 1;
 258                 $flt = 1 if ($mq < $opts{q});
 259                 # filtering SNPs
 260                 if ($t[5] >= $opts{G}) {
 261                   for my $x (@staging) {
 262             # Is it a SNP and is it outside the SNP filter window?
 263                         next if ($x->[0] >= 0 || $x->[4] + $x->[2] + $ow < $t[1]);
 264                         $x->[1] = 5 if ($x->[1] == 0);
 265                   }
 266                 }
 267                 # the indel filtering score
 268                 $score = $t[5];
 269                 # check the staging list for indel filtering
 270                 for my $x (@staging) {
 271           # Is it a SNP and is it outside the gap filter window
 272                   next if ($x->[0] < 0 || $x->[4] + $x->[2] + $ol < $t[1]);
 273                   if ($x->[0] < $score) {
 274                         $x->[1] = 6;
 275                   } else {
 276                         $flt = 6; last;
 277                   }
 278                 }
 279           } else { # a SNP
 280                 $flt = 1 if ($mq < $opts{Q});
 281                 # check adjacent SNPs
 282                 my $k = 1;
 283                 for my $x (@staging) {
 284                   ++$k if ($x->[0] < 0 && -($x->[0] + 1) > $opts{F} && $x->[4] + $x->[2] + $oW >= $t[1] && ($x->[1] == 0 || $x->[1] == 4 || $x->[1] == 5));
 285                 }
 286                 # filtering is necessary
 287                 if ($k > $opts{N}) {
 288                   $flt = 4;
 289                   for my $x (@staging) {
 290                          $x->[1] = 4 if ($x->[0] < 0 && $x->[4] + $x->[2] + $oW >= $t[1] && $x->[1] == 0);
 291                   }
 292                 } else { # then check gap filter
 293                   for my $x (@staging) {
 294                         next if ($x->[0] < 0 || $x->[4] + $x->[2] + $ow < $t[1]);
 295                         if ($x->[0] >= $opts{G}) {
 296                           $flt = 5; last;
 297                         }
 298                   }
 299                 }
 300           }
 301         }
 302         push(@staging, [$score < 0? -$af-1 : $score, $flt, $dlen, @t]);
 303   }
 304   # output the last few elements in the staging list
 305   while (@staging) {
 306         varFilter_aux(shift @staging, $opts{p});
 307   }
 308 }
 309
 310 sub varFilter_aux {
 311   my ($first, $is_print) = @_;
 312   if ($first->[1] == 0) {
 313         print join("\t", @$first[3 .. @$first-1]), "\n";
 314   } elsif ($is_print) {
 315         print STDERR join("\t", substr("UQdDWGgsiX", $first->[1], 1), @$first[3 .. @$first-1]), "\n";
 316   }
 317 }
 318
 319 sub filter4vcf {
 320   my %opts = (d=>3, D=>2000, 1=>1e-4, 2=>1e-100, 3=>0, 4=>1e-4, Q=>10, q=>3);
 321   getopts('d:D:1:2:3:4:Q:q:', \%opts);
 322   die(qq/
 323 Usage:   vcfutils.pl filter4vcf [options] <in.vcf>
 324
 325 Options: -d INT     min total depth (given DP or DP4) [$opts{d}]
 326          -D INT     max total depth [$opts{D}]
 327          -q INT     min SNP quality [$opts{q}]
 328          -Q INT     min RMS mapQ (given MQ) [$opts{Q}]
 329          -1 FLOAT   min P-value for strand bias (given PV4) [$opts{1}]
 330          -2 FLOAT   min P-value for baseQ bias [$opts{2}]
 331          -3 FLOAT   min P-value for mapQ bias [$opts{3}]
 332          -4 FLOAT   min P-value for end distance bias [$opts{4}]\n
 333 /) if (@ARGV == 0 && -t STDIN);
 334
 335   my %ts = (AG=>1, GA=>1, CT=>1, TC=>1);
 336
 337   my @n = (0, 0);
 338   while (<>) {
 339         next if (/^#/);
 340         next if (/PV4=([^,]+),([^,]+),([^,]+),([^,;\t]+)/ && ($1<$opts{1} || $2<$opts{2} || $3<$opts{3} || $4<$opts{4}));
 341         my $depth = -1;
 342         $depth = $1 if (/DP=(\d+)/);
 343         $depth = $1+$2+$3+$4 if (/DP4=(\d+),(\d+),(\d+),(\d+)/);
 344         next if ($depth > 0 && ($depth < $opts{d} || $depth > $opts{D}));
 345         next if (/MQ=(\d+)/ && $1 < $opts{Q});
 346         my @t = split;
 347         next if ($t[5] >= 0 && $t[5] < $opts{q});
 348         ++$n[0];
 349         my @s = split(',', $t[4]);
 350         ++$n[1] if ($ts{$t[3].$s[0]});
 351         print;
 352   }
 353 }
 354
 355 sub ucscsnp2vcf {
 356   die("Usage: vcfutils.pl <in.ucsc.snp>\n") if (@ARGV == 0 && -t STDIN);
 357   print "##fileformat=VCFv4.0\n";
 358   print join("\t", "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO"), "\n";
 359   while (<>) {
 360         my @t = split("\t");
 361         my $indel = ($t[9] =~ /^[ACGT](\/[ACGT])+$/)? 0 : 1;
 362         my $pos = $t[2] + 1;
 363         my @alt;
 364         push(@alt, $t[7]);
 365         if ($t[6] eq '-') {
 366           $t[9] = reverse($t[9]);
 367           $t[9] =~ tr/ACGTRYMKWSNacgtrymkwsn/TGCAYRKMWSNtgcayrkmwsn/;
 368         }
 369         my @a = split("/", $t[9]);
 370         for (@a) {
 371           push(@alt, $_) if ($_ ne $alt[0]);
 372         }
 373         if ($indel) {
 374           --$pos;
 375           for (0 .. $#alt) {
 376                 $alt[$_] =~ tr/-//d;
 377                 $alt[$_] = "N$alt[$_]";
 378           }
 379         }
 380         my $ref = shift(@alt);
 381         my $af = $t[13] > 0? ";AF=$t[13]" : '';
 382         my $valid = ($t[12] eq 'unknown')? '' : ";valid=$t[12]";
 383         my $info = "molType=$t[10];class=$t[11]$valid$af";
 384         print join("\t", $t[1], $pos, $t[4], $ref, join(",", @alt), 0, '.', $info), "\n";
 385   }
 386 }
 387
 388 sub hapmap2vcf {
 389   die("Usage: vcfutils.pl <in.ucsc.snp> <in.hapmap>\n") if (@ARGV == 0);
 390   my $fn = shift(@ARGV);
 391   # parse UCSC SNP
 392   warn("Parsing UCSC SNPs...\n");
 393   my ($fh, %map);
 394   open($fh, ($fn =~ /\.gz$/)? "gzip -dc $fn |" : $fn) || die;
 395   while (<$fh>) {
 396         my @t = split;
 397         next if ($t[3] - $t[2] != 1); # not SNP
 398         @{$map{$t[4]}} = @t[1,3,7];
 399   }
 400   close($fh);
 401   # write VCF
 402   warn("Writing VCF...\n");
 403   print "##fileformat=VCFv4.0\n";
 404   while (<>) {
 405         my @t = split;
 406         if ($t[0] eq 'rs#') { # the first line
 407           print join("\t", "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT", @t[11..$#t]), "\n";
 408         } else {
 409           next unless ($map{$t[0]});
 410           next if (length($t[1]) != 3); # skip non-SNPs
 411           my $a = \@{$map{$t[0]}};
 412           my $ref = $a->[2];
 413           my @u = split('/', $t[1]);
 414           if ($u[1] eq $ref) {
 415                 $u[1] = $u[0]; $u[0] = $ref;
 416           } elsif ($u[0] ne $ref) { next; }
 417           my $alt = $u[1];
 418           my %w;
 419           $w{$u[0]} = 0; $w{$u[1]} = 1;
 420           my @s = (@$a[0,1], $t[0], $ref, $alt, 0, '.', '.', 'GT');
 421           my $is_tri = 0;
 422           for (@t[11..$#t]) {
 423                 if ($_ eq 'NN') {
 424                   push(@s, './.');
 425                 } else {
 426                   my @a = ($w{substr($_,0,1)}, $w{substr($_,1,1)});
 427                   if (!defined($a[0]) || !defined($a[1])) {
 428                         $is_tri = 1;
 429                         last;
 430                   }
 431                   push(@s, "$a[0]/$a[1]");
 432                 }
 433           }
 434           next if ($is_tri);
 435           print join("\t", @s), "\n";
 436         }
 437   }
 438 }
 439
 440 sub usage {
 441   die(qq/
 442 Usage:   vcfutils.pl <command> [<arguments>]\n
 443 Command: subsam       get a subset of samples
 444          listsam      list the samples
 445          fillac       fill the allele count field
 446          qstats       SNP stats stratified by QUAL
 447          varFilter    filtering short variants
 448          filter4vcf   filtering VCFs produced by samtools+bcftools
 449          hapmap2vcf   convert the hapmap format to VCF
 450          ucscsnp2vcf  convert UCSC SNP SQL dump to VCF
 451 \n/);
 452 }