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, ldstats=>\&ldstats);
  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=>1, D=>10000, l=>10, Q=>10, s=>100, w=>10, p=>undef, F=>.001);
 203   getopts('pd:D:l:Q:w:G:F:', \%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          -w INT    SNP within INT bp around a gap to be filtered [$opts{w}]
 211          -l INT    window size for filtering adjacent gaps [$opts{l}]
 212          -p        print filtered variants
 213 \n/) if (@ARGV == 0 && -t STDIN);
 214
 215   # calculate the window size
 216   my ($ol, $ow) = ($opts{l}, $opts{w});
 217   my $max_dist = $ol > $ow? $ol : $ow;
 218   # the core loop
 219   my @staging; # (indel_filtering_score, flt_tag, indel_span; chr, pos, ...)
 220   while (<>) {
 221         my @t = split;
 222         next if (/^#/);
 223         next if ($t[4] eq '.'); # skip non-var sites
 224         my $is_snp = 1;
 225         if (length($t[3]) > 1) {
 226           $is_snp = 0;
 227         } else {
 228           my @s = split(',', $t[4]);
 229           for (@s) {
 230                 $is_snp = 0 if (length > 1);
 231           }
 232         }
 233         # clear the out-of-range elements
 234         while (@staging) {
 235       # Still on the same chromosome and the first element's window still affects this position?
 236           last if ($staging[0][3] eq $t[0] && $staging[0][4] + $staging[0][2] + $max_dist >= $t[1]);
 237           varFilter_aux(shift(@staging), $opts{p}); # calling a function is a bit slower, not much
 238         }
 239         my $flt = 0;
 240
 241         # parse annotations
 242         my ($dp, $mq, $dp_alt) = (-1, -1, -1);
 243         if ($t[7] =~ /DP=(\d+)/i) {
 244           $dp = $1;
 245         } elsif ($t[7] =~ /DP4=(\d+),(\d+),(\d+),(\d+)/i) {
 246           $dp = $1 + $2 + $3 + $4;
 247           $dp_alt = $3 + $4;
 248         }
 249         $mq = $1 if ($t[7] =~ /MQ=(\d+)/i);
 250         # the depth and mapQ filter
 251         if ($dp >= 0) {
 252           if ($dp < $opts{d}) {
 253                 $flt = 2;
 254           } elsif ($dp > $opts{D}) {
 255                 $flt = 3;
 256           }
 257         }
 258         $flt = 1 if ($flt == 0 && $mq >= 0 && $mq < $opts{Q});
 259
 260         # site dependent filters
 261         my ($rlen, $indel_score) = (0, -1); # $indel_score<0 for SNPs
 262         if ($flt == 0) {
 263           if (!$is_snp) { # an indel
 264                 $rlen = length($t[3]) - 1;
 265                 $indel_score = $t[5] * 100 + $dp_alt;
 266                 # filtering SNPs
 267                 for my $x (@staging) {
 268                   next if ($x->[0] >= 0 || $x->[1] || $x->[4] + $x->[2] + $ow < $t[1]);
 269                   $x->[1] = 5;
 270                 }
 271                 # check the staging list for indel filtering
 272                 for my $x (@staging) {
 273                   next if ($x->[0] < 0 || $x->[1] || $x->[4] + $x->[2] + $ol < $t[1]);
 274                   if ($x->[0] < $indel_score) {
 275                         $x->[1] = 6;
 276                   } else {
 277                         $flt = 6; last;
 278                   }
 279                 }
 280           } else { # a SNP
 281                 for my $x (@staging) {
 282                   next if ($x->[0] < 0 || $x->[1] || $x->[4] + $x->[2] + $ow < $t[1]);
 283                   $flt = 5;
 284                   last;
 285                 }
 286           }
 287         }
 288         push(@staging, [$indel_score, $flt, $rlen, @t]);
 289   }
 290   # output the last few elements in the staging list
 291   while (@staging) {
 292         varFilter_aux(shift @staging, $opts{p});
 293   }
 294 }
 295
 296 sub varFilter_aux {
 297   my ($first, $is_print) = @_;
 298   if ($first->[1] == 0) {
 299         print join("\t", @$first[3 .. @$first-1]), "\n";
 300   } elsif ($is_print) {
 301         print STDERR join("\t", substr("UQdDWGgsiX", $first->[1], 1), @$first[3 .. @$first-1]), "\n";
 302   }
 303 }
 304
 305 sub filter4vcf {
 306   my %opts = (d=>3, D=>2000, 1=>1e-4, 2=>1e-100, 3=>0, 4=>1e-4, Q=>10, q=>3);
 307   getopts('d:D:1:2:3:4:Q:q:', \%opts);
 308   die(qq/
 309 Usage:   vcfutils.pl filter4vcf [options] <in.vcf>
 310
 311 Options: -d INT     min total depth (given DP or DP4) [$opts{d}]
 312          -D INT     max total depth [$opts{D}]
 313          -q INT     min SNP quality [$opts{q}]
 314          -Q INT     min RMS mapQ (given MQ) [$opts{Q}]
 315          -1 FLOAT   min P-value for strand bias (given PV4) [$opts{1}]
 316          -2 FLOAT   min P-value for baseQ bias [$opts{2}]
 317          -3 FLOAT   min P-value for mapQ bias [$opts{3}]
 318          -4 FLOAT   min P-value for end distance bias [$opts{4}]\n
 319 /) if (@ARGV == 0 && -t STDIN);
 320
 321   my %ts = (AG=>1, GA=>1, CT=>1, TC=>1);
 322
 323   my @n = (0, 0);
 324   while (<>) {
 325         if (/^#/) {
 326           print;
 327           next;
 328         }
 329         next if (/PV4=([^,]+),([^,]+),([^,]+),([^,;\t]+)/ && ($1<$opts{1} || $2<$opts{2} || $3<$opts{3} || $4<$opts{4}));
 330         my $depth = -1;
 331         $depth = $1 if (/DP=(\d+)/);
 332         $depth = $1+$2+$3+$4 if (/DP4=(\d+),(\d+),(\d+),(\d+)/);
 333         next if ($depth > 0 && ($depth < $opts{d} || $depth > $opts{D}));
 334         next if (/MQ=(\d+)/ && $1 < $opts{Q});
 335         my @t = split;
 336         next if ($t[5] >= 0 && $t[5] < $opts{q});
 337         ++$n[0];
 338         my @s = split(',', $t[4]);
 339         ++$n[1] if ($ts{$t[3].$s[0]});
 340         print;
 341   }
 342 }
 343
 344 sub ucscsnp2vcf {
 345   die("Usage: vcfutils.pl <in.ucsc.snp>\n") if (@ARGV == 0 && -t STDIN);
 346   print "##fileformat=VCFv4.0\n";
 347   print join("\t", "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO"), "\n";
 348   while (<>) {
 349         my @t = split("\t");
 350         my $indel = ($t[9] =~ /^[ACGT](\/[ACGT])+$/)? 0 : 1;
 351         my $pos = $t[2] + 1;
 352         my @alt;
 353         push(@alt, $t[7]);
 354         if ($t[6] eq '-') {
 355           $t[9] = reverse($t[9]);
 356           $t[9] =~ tr/ACGTRYMKWSNacgtrymkwsn/TGCAYRKMWSNtgcayrkmwsn/;
 357         }
 358         my @a = split("/", $t[9]);
 359         for (@a) {
 360           push(@alt, $_) if ($_ ne $alt[0]);
 361         }
 362         if ($indel) {
 363           --$pos;
 364           for (0 .. $#alt) {
 365                 $alt[$_] =~ tr/-//d;
 366                 $alt[$_] = "N$alt[$_]";
 367           }
 368         }
 369         my $ref = shift(@alt);
 370         my $af = $t[13] > 0? ";AF=$t[13]" : '';
 371         my $valid = ($t[12] eq 'unknown')? '' : ";valid=$t[12]";
 372         my $info = "molType=$t[10];class=$t[11]$valid$af";
 373         print join("\t", $t[1], $pos, $t[4], $ref, join(",", @alt), 0, '.', $info), "\n";
 374   }
 375 }
 376
 377 sub hapmap2vcf {
 378   die("Usage: vcfutils.pl <in.ucsc.snp> <in.hapmap>\n") if (@ARGV == 0);
 379   my $fn = shift(@ARGV);
 380   # parse UCSC SNP
 381   warn("Parsing UCSC SNPs...\n");
 382   my ($fh, %map);
 383   open($fh, ($fn =~ /\.gz$/)? "gzip -dc $fn |" : $fn) || die;
 384   while (<$fh>) {
 385         my @t = split;
 386         next if ($t[3] - $t[2] != 1); # not SNP
 387         @{$map{$t[4]}} = @t[1,3,7];
 388   }
 389   close($fh);
 390   # write VCF
 391   warn("Writing VCF...\n");
 392   print "##fileformat=VCFv4.0\n";
 393   while (<>) {
 394         my @t = split;
 395         if ($t[0] eq 'rs#') { # the first line
 396           print join("\t", "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT", @t[11..$#t]), "\n";
 397         } else {
 398           next unless ($map{$t[0]});
 399           next if (length($t[1]) != 3); # skip non-SNPs
 400           my $a = \@{$map{$t[0]}};
 401           my $ref = $a->[2];
 402           my @u = split('/', $t[1]);
 403           if ($u[1] eq $ref) {
 404                 $u[1] = $u[0]; $u[0] = $ref;
 405           } elsif ($u[0] ne $ref) { next; }
 406           my $alt = $u[1];
 407           my %w;
 408           $w{$u[0]} = 0; $w{$u[1]} = 1;
 409           my @s = (@$a[0,1], $t[0], $ref, $alt, 0, '.', '.', 'GT');
 410           my $is_tri = 0;
 411           for (@t[11..$#t]) {
 412                 if ($_ eq 'NN') {
 413                   push(@s, './.');
 414                 } else {
 415                   my @a = ($w{substr($_,0,1)}, $w{substr($_,1,1)});
 416                   if (!defined($a[0]) || !defined($a[1])) {
 417                         $is_tri = 1;
 418                         last;
 419                   }
 420                   push(@s, "$a[0]/$a[1]");
 421                 }
 422           }
 423           next if ($is_tri);
 424           print join("\t", @s), "\n";
 425         }
 426   }
 427 }
 428
 429 sub usage {
 430   die(qq/
 431 Usage:   vcfutils.pl <command> [<arguments>]\n
 432 Command: subsam       get a subset of samples
 433          listsam      list the samples
 434          fillac       fill the allele count field
 435          qstats       SNP stats stratified by QUAL
 436          varFilter    filtering short variants
 437          filter4vcf   filtering VCFs produced by samtools+bcftools
 438          hapmap2vcf   convert the hapmap format to VCF
 439          ucscsnp2vcf  convert UCSC SNP SQL dump to VCF
 440 \n/);
 441 }