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);
  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 ldstats {
 107   my %opts = (t=>0.9);
 108   getopts('t:', \%opts);
 109   die("Usage: vcfutils.pl ldstats [-t $opts{t}] <in.vcf>\n") if (@ARGV == 0 && -t STDIN);
 110   my $cutoff = $opts{t};
 111   my ($last, $lastchr) = (0x7fffffff, '');
 112   my ($x, $y, $n) = (0, 0, 0);
 113   while (<>) {
 114         if (/^([^#\s]+)\s(\d+)/) {
 115           my ($chr, $pos) = ($1, $2);
 116           if (/NEIR=([\d\.]+)/) {
 117                 ++$n;
 118                 ++$y, $x += $pos - $last if ($lastchr eq $chr && $pos > $last && $1 > $cutoff);
 119           }
 120           $last = $pos; $lastchr = $chr;
 121         }
 122   }
 123   print "Number of SNP intervals in strong LD (r > $opts{t}): $y\n";
 124   print "Fraction: ", $y/$n, "\n";
 125   print "Length: $x\n";
 126 }
 127
 128 sub qstats {
 129   my %opts = (r=>'', s=>0.02, v=>undef);
 130   getopts('r:s:v', \%opts);
 131   die("Usage: vcfutils.pl qstats [-r ref.vcf] <in.vcf>\n
 132 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);
 133   my %ts = (AG=>1, GA=>1, CT=>1, TC=>1);
 134   my %h = ();
 135   my $is_vcf = defined($opts{v})? 1 : 0;
 136   if ($opts{r}) { # read the reference positions
 137         my $fh;
 138         open($fh, $opts{r}) || die;
 139         while (<$fh>) {
 140           next if (/^#/);
 141           if ($is_vcf) {
 142                 my @t = split;
 143                 $h{$t[0],$t[1]} = $t[4];
 144           } else {
 145                 $h{$1,$2} = 1 if (/^(\S+)\s+(\d+)/);
 146           }
 147         }
 148         close($fh);
 149   }
 150   my $hsize = scalar(keys %h);
 151   my @a;
 152   while (<>) {
 153         next if (/^#/);
 154         my @t = split;
 155         next if (length($t[3]) != 1 || uc($t[3]) eq 'N');
 156         $t[3] = uc($t[3]); $t[4] = uc($t[4]);
 157         my @s = split(',', $t[4]);
 158         $t[5] = 3 if ($t[5] eq '.' || $t[5] < 0);
 159         next if (length($s[0]) != 1);
 160         my $hit;
 161         if ($is_vcf) {
 162           $hit = 0;
 163           my $aa = $h{$t[0],$t[1]};
 164           if (defined($aa)) {
 165                 my @aaa = split(",", $aa);
 166                 for (@aaa) {
 167                   $hit = 1 if ($_ eq $s[0]);
 168                 }
 169           }
 170         } else {
 171           $hit = defined($h{$t[0],$t[1]})? 1 : 0;
 172         }
 173         push(@a, [$t[5], ($t[4] eq '.' || $t[4] eq $t[3])? 0 : 1, $ts{$t[3].$s[0]}? 1 : 0, $hit]);
 174   }
 175   push(@a, [-1, 0, 0, 0]); # end marker
 176   die("[qstats] No SNP data!\n") if (@a == 0);
 177   @a = sort {$b->[0]<=>$a->[0]} @a;
 178   my $next = $opts{s};
 179   my $last = $a[0];
 180   my @c = (0, 0, 0, 0);
 181   my @lc;
 182   $lc[1] = $lc[2] = 0;
 183   for my $p (@a) {
 184         if ($p->[0] == -1 || ($p->[0] != $last && $c[0]/@a > $next)) {
 185           my @x;
 186           $x[0] = sprintf("%.4f", $c[1]-$c[2]? $c[2] / ($c[1] - $c[2]) : 100);
 187           $x[1] = sprintf("%.4f", $hsize? $c[3] / $hsize : 0);
 188           $x[2] = sprintf("%.4f", $c[3] / $c[1]);
 189           my $a = $c[1] - $lc[1];
 190           my $b = $c[2] - $lc[2];
 191           $x[3] = sprintf("%.4f", $a-$b? $b / ($a-$b) : 100);
 192           print join("\t", $last, @c, @x), "\n";
 193           $next = $c[0]/@a + $opts{s};
 194           $lc[1] = $c[1]; $lc[2] = $c[2];
 195         }
 196         ++$c[0]; $c[1] += $p->[1]; $c[2] += $p->[2]; $c[3] += $p->[3];
 197         $last = $p->[0];
 198   }
 199 }
 200
 201 sub varFilter {
 202   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);
 203   getopts('pd:D:W:Q:w:a:1:2:3:4:', \%opts);
 204   die(qq/
 205 Usage:   vcfutils.pl varFilter [options] <in.vcf>
 206
 207 Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
 208          -d INT    minimum read depth [$opts{d}]
 209          -D INT    maximum read depth [$opts{D}]
 210          -a INT    minimum number of alternate bases [$opts{a}]
 211          -w INT    SNP within INT bp around a gap to be filtered [$opts{w}]
 212          -W INT    window size for filtering adjacent gaps [$opts{W}]
 213          -1 FLOAT  min P-value for strand bias (given PV4) [$opts{1}]
 214          -2 FLOAT  min P-value for baseQ bias [$opts{2}]
 215          -3 FLOAT  min P-value for mapQ bias [$opts{3}]
 216          -4 FLOAT  min P-value for end distance bias [$opts{4}]
 217          -p        print filtered variants
 218
 219 Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools.
 220 \n/) if (@ARGV == 0 && -t STDIN);
 221
 222   # calculate the window size
 223   my ($ol, $ow) = ($opts{W}, $opts{w});
 224   my $max_dist = $ol > $ow? $ol : $ow;
 225   # the core loop
 226   my @staging; # (indel_filtering_score, flt_tag, indel_span; chr, pos, ...)
 227   while (<>) {
 228         my @t = split;
 229     if (/^#/) {
 230           print; next;
 231         }
 232         next if ($t[4] eq '.'); # skip non-var sites
 233         # check if the site is a SNP
 234         my $is_snp = 1;
 235         if (length($t[3]) > 1) {
 236           $is_snp = 0;
 237         } else {
 238           my @s = split(',', $t[4]);
 239           for (@s) {
 240                 $is_snp = 0 if (length > 1);
 241           }
 242         }
 243         # clear the out-of-range elements
 244         while (@staging) {
 245       # Still on the same chromosome and the first element's window still affects this position?
 246           last if ($staging[0][3] eq $t[0] && $staging[0][4] + $staging[0][2] + $max_dist >= $t[1]);
 247           varFilter_aux(shift(@staging), $opts{p}); # calling a function is a bit slower, not much
 248         }
 249         my $flt = 0;
 250         # parse annotations
 251         my ($dp, $mq, $dp_alt) = (-1, -1, -1);
 252         if ($t[7] =~ /DP=(\d+)/i) {
 253           $dp = $1;
 254         } elsif ($t[7] =~ /DP4=(\d+),(\d+),(\d+),(\d+)/i) {
 255           $dp = $1 + $2 + $3 + $4;
 256           $dp_alt = $3 + $4;
 257         }
 258         $mq = $1 if ($t[7] =~ /MQ=(\d+)/i);
 259         # the depth and mapQ filter
 260         if ($dp >= 0) {
 261           if ($dp < $opts{d}) {
 262                 $flt = 2;
 263           } elsif ($dp > $opts{D}) {
 264                 $flt = 3;
 265           }
 266         }
 267         $flt = 4 if ($dp_alt >= 0 && $dp_alt < $opts{a});
 268         $flt = 1 if ($flt == 0 && $mq >= 0 && $mq < $opts{Q});
 269         $flt = 7 if ($flt == 0 && /PV4=([^,]+),([^,]+),([^,]+),([^,;\t]+)/
 270                                  && ($1<$opts{1} || $2<$opts{2} || $3<$opts{3} || $4<$opts{4}));
 271
 272         # site dependent filters
 273         my ($rlen, $indel_score) = (0, -1); # $indel_score<0 for SNPs
 274         if ($flt == 0) {
 275           if (!$is_snp) { # an indel
 276                 $rlen = length($t[3]) - 1;
 277                 $indel_score = $t[5] * 100 + $dp_alt;
 278                 # filtering SNPs
 279                 for my $x (@staging) {
 280                   next if ($x->[0] >= 0 || $x->[1] || $x->[4] + $x->[2] + $ow < $t[1]);
 281                   $x->[1] = 5;
 282                 }
 283                 # check the staging list for indel filtering
 284                 for my $x (@staging) {
 285                   next if ($x->[0] < 0 || $x->[1] || $x->[4] + $x->[2] + $ol < $t[1]);
 286                   if ($x->[0] < $indel_score) {
 287                         $x->[1] = 6;
 288                   } else {
 289                         $flt = 6; last;
 290                   }
 291                 }
 292           } else { # a SNP
 293                 for my $x (@staging) {
 294                   next if ($x->[0] < 0 || $x->[1] || $x->[4] + $x->[2] + $ow < $t[1]);
 295                   $flt = 5;
 296                   last;
 297                 }
 298           }
 299         }
 300         push(@staging, [$indel_score, $flt, $rlen, @t]);
 301   }
 302   # output the last few elements in the staging list
 303   while (@staging) {
 304         varFilter_aux(shift @staging, $opts{p});
 305   }
 306 }
 307
 308 sub varFilter_aux {
 309   my ($first, $is_print) = @_;
 310   if ($first->[1] == 0) {
 311         print join("\t", @$first[3 .. @$first-1]), "\n";
 312   } elsif ($is_print) {
 313         print STDERR join("\t", substr("UQdDaGgP", $first->[1], 1), @$first[3 .. @$first-1]), "\n";
 314   }
 315 }
 316
 317 sub gapstats {
 318   my (@c0, @c1);
 319   $c0[$_] = $c1[$_] = 0 for (0 .. 10000);
 320   while (<>) {
 321         next if (/^#/);
 322         my @t = split;
 323         next if (length($t[3]) == 1 && $t[4] =~ /^[A-Za-z](,[A-Za-z])*$/); # not an indel
 324         my @s = split(',', $t[4]);
 325         for my $x (@s) {
 326           my $l = length($x) - length($t[3]) + 5000;
 327           $c0[$l] += 1 / @s;
 328         }
 329   }
 330   for (my $i = 0; $i < 10000; ++$i) {
 331         next if ($c0[$i] == 0);
 332         printf("%d\t%.2f\n", ($i-5000), $c0[$i]);
 333   }
 334 }
 335
 336 sub ucscsnp2vcf {
 337   die("Usage: vcfutils.pl <in.ucsc.snp>\n") if (@ARGV == 0 && -t STDIN);
 338   print "##fileformat=VCFv4.0\n";
 339   print join("\t", "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO"), "\n";
 340   while (<>) {
 341         my @t = split("\t");
 342         my $indel = ($t[9] =~ /^[ACGT](\/[ACGT])+$/)? 0 : 1;
 343         my $pos = $t[2] + 1;
 344         my @alt;
 345         push(@alt, $t[7]);
 346         if ($t[6] eq '-') {
 347           $t[9] = reverse($t[9]);
 348           $t[9] =~ tr/ACGTRYMKWSNacgtrymkwsn/TGCAYRKMWSNtgcayrkmwsn/;
 349         }
 350         my @a = split("/", $t[9]);
 351         for (@a) {
 352           push(@alt, $_) if ($_ ne $alt[0]);
 353         }
 354         if ($indel) {
 355           --$pos;
 356           for (0 .. $#alt) {
 357                 $alt[$_] =~ tr/-//d;
 358                 $alt[$_] = "N$alt[$_]";
 359           }
 360         }
 361         my $ref = shift(@alt);
 362         my $af = $t[13] > 0? ";AF=$t[13]" : '';
 363         my $valid = ($t[12] eq 'unknown')? '' : ";valid=$t[12]";
 364         my $info = "molType=$t[10];class=$t[11]$valid$af";
 365         print join("\t", $t[1], $pos, $t[4], $ref, join(",", @alt), 0, '.', $info), "\n";
 366   }
 367 }
 368
 369 sub hapmap2vcf {
 370   die("Usage: vcfutils.pl <in.ucsc.snp> <in.hapmap>\n") if (@ARGV == 0);
 371   my $fn = shift(@ARGV);
 372   # parse UCSC SNP
 373   warn("Parsing UCSC SNPs...\n");
 374   my ($fh, %map);
 375   open($fh, ($fn =~ /\.gz$/)? "gzip -dc $fn |" : $fn) || die;
 376   while (<$fh>) {
 377         my @t = split;
 378         next if ($t[3] - $t[2] != 1); # not SNP
 379         @{$map{$t[4]}} = @t[1,3,7];
 380   }
 381   close($fh);
 382   # write VCF
 383   warn("Writing VCF...\n");
 384   print "##fileformat=VCFv4.0\n";
 385   while (<>) {
 386         my @t = split;
 387         if ($t[0] eq 'rs#') { # the first line
 388           print join("\t", "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT", @t[11..$#t]), "\n";
 389         } else {
 390           next unless ($map{$t[0]});
 391           next if (length($t[1]) != 3); # skip non-SNPs
 392           my $a = \@{$map{$t[0]}};
 393           my $ref = $a->[2];
 394           my @u = split('/', $t[1]);
 395           if ($u[1] eq $ref) {
 396                 $u[1] = $u[0]; $u[0] = $ref;
 397           } elsif ($u[0] ne $ref) { next; }
 398           my $alt = $u[1];
 399           my %w;
 400           $w{$u[0]} = 0; $w{$u[1]} = 1;
 401           my @s = (@$a[0,1], $t[0], $ref, $alt, 0, '.', '.', 'GT');
 402           my $is_tri = 0;
 403           for (@t[11..$#t]) {
 404                 if ($_ eq 'NN') {
 405                   push(@s, './.');
 406                 } else {
 407                   my @a = ($w{substr($_,0,1)}, $w{substr($_,1,1)});
 408                   if (!defined($a[0]) || !defined($a[1])) {
 409                         $is_tri = 1;
 410                         last;
 411                   }
 412                   push(@s, "$a[0]/$a[1]");
 413                 }
 414           }
 415           next if ($is_tri);
 416           print join("\t", @s), "\n";
 417         }
 418   }
 419 }
 420
 421 sub usage {
 422   die(qq/
 423 Usage:   vcfutils.pl <command> [<arguments>]\n
 424 Command: subsam       get a subset of samples
 425          listsam      list the samples
 426          fillac       fill the allele count field
 427          qstats       SNP stats stratified by QUAL
 428          varFilter    filtering short variants
 429          hapmap2vcf   convert the hapmap format to VCF
 430          ucscsnp2vcf  convert UCSC SNP SQL dump to VCF
 431 \n/);
 432 }