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] =~ /DP4=(\d+),(\d+),(\d+),(\d+)/i) {
 253           $dp = $1 + $2 + $3 + $4;
 254           $dp_alt = $3 + $4;
 255         }
 256         if ($t[7] =~ /DP=(\d+)/i) {
 257           $dp = $1;
 258         }
 259         $mq = $1 if ($t[7] =~ /MQ=(\d+)/i);
 260         # the depth and mapQ filter
 261         if ($dp >= 0) {
 262           if ($dp < $opts{d}) {
 263                 $flt = 2;
 264           } elsif ($dp > $opts{D}) {
 265                 $flt = 3;
 266           }
 267         }
 268         $flt = 4 if ($dp_alt >= 0 && $dp_alt < $opts{a});
 269         $flt = 1 if ($flt == 0 && $mq >= 0 && $mq < $opts{Q});
 270         $flt = 7 if ($flt == 0 && /PV4=([^,]+),([^,]+),([^,]+),([^,;\t]+)/
 271                                  && ($1<$opts{1} || $2<$opts{2} || $3<$opts{3} || $4<$opts{4}));
 272
 273         # site dependent filters
 274         my ($rlen, $indel_score) = (0, -1); # $indel_score<0 for SNPs
 275         if ($flt == 0) {
 276           if (!$is_snp) { # an indel
 277                 $rlen = length($t[3]) - 1;
 278                 $indel_score = $t[5] * 100 + $dp_alt;
 279                 # filtering SNPs
 280                 for my $x (@staging) {
 281                   next if ($x->[0] >= 0 || $x->[1] || $x->[4] + $x->[2] + $ow < $t[1]);
 282                   $x->[1] = 5;
 283                 }
 284                 # check the staging list for indel filtering
 285                 for my $x (@staging) {
 286                   next if ($x->[0] < 0 || $x->[1] || $x->[4] + $x->[2] + $ol < $t[1]);
 287                   if ($x->[0] < $indel_score) {
 288                         $x->[1] = 6;
 289                   } else {
 290                         $flt = 6; last;
 291                   }
 292                 }
 293           } else { # a SNP
 294                 for my $x (@staging) {
 295                   next if ($x->[0] < 0 || $x->[1] || $x->[4] + $x->[2] + $ow < $t[1]);
 296                   $flt = 5;
 297                   last;
 298                 }
 299           }
 300         }
 301         push(@staging, [$indel_score, $flt, $rlen, @t]);
 302   }
 303   # output the last few elements in the staging list
 304   while (@staging) {
 305         varFilter_aux(shift @staging, $opts{p});
 306   }
 307 }
 308
 309 sub varFilter_aux {
 310   my ($first, $is_print) = @_;
 311   if ($first->[1] == 0) {
 312         print join("\t", @$first[3 .. @$first-1]), "\n";
 313   } elsif ($is_print) {
 314         print STDERR join("\t", substr("UQdDaGgP", $first->[1], 1), @$first[3 .. @$first-1]), "\n";
 315   }
 316 }
 317
 318 sub gapstats {
 319   my (@c0, @c1);
 320   $c0[$_] = $c1[$_] = 0 for (0 .. 10000);
 321   while (<>) {
 322         next if (/^#/);
 323         my @t = split;
 324         next if (length($t[3]) == 1 && $t[4] =~ /^[A-Za-z](,[A-Za-z])*$/); # not an indel
 325         my @s = split(',', $t[4]);
 326         for my $x (@s) {
 327           my $l = length($x) - length($t[3]) + 5000;
 328           if ($x =~ /^-/) {
 329                 $l = -(length($x) - 1) + 5000;
 330           } elsif ($x =~ /^\+/) {
 331                 $l = length($x) - 1 + 5000;
 332           }
 333           $c0[$l] += 1 / @s;
 334         }
 335   }
 336   for (my $i = 0; $i < 10000; ++$i) {
 337         next if ($c0[$i] == 0);
 338         $c1[0] += $c0[$i];
 339         $c1[1] += $c0[$i] if (($i-5000)%3 == 0);
 340         printf("C\t%d\t%.2f\n", ($i-5000), $c0[$i]);
 341   }
 342   printf("3\t%d\t%d\t%.3f\n", $c1[0], $c1[1], $c1[1]/$c1[0]);
 343 }
 344
 345 sub ucscsnp2vcf {
 346   die("Usage: vcfutils.pl <in.ucsc.snp>\n") if (@ARGV == 0 && -t STDIN);
 347   print "##fileformat=VCFv4.0\n";
 348   print join("\t", "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO"), "\n";
 349   while (<>) {
 350         my @t = split("\t");
 351         my $indel = ($t[9] =~ /^[ACGT](\/[ACGT])+$/)? 0 : 1;
 352         my $pos = $t[2] + 1;
 353         my @alt;
 354         push(@alt, $t[7]);
 355         if ($t[6] eq '-') {
 356           $t[9] = reverse($t[9]);
 357           $t[9] =~ tr/ACGTRYMKWSNacgtrymkwsn/TGCAYRKMWSNtgcayrkmwsn/;
 358         }
 359         my @a = split("/", $t[9]);
 360         for (@a) {
 361           push(@alt, $_) if ($_ ne $alt[0]);
 362         }
 363         if ($indel) {
 364           --$pos;
 365           for (0 .. $#alt) {
 366                 $alt[$_] =~ tr/-//d;
 367                 $alt[$_] = "N$alt[$_]";
 368           }
 369         }
 370         my $ref = shift(@alt);
 371         my $af = $t[13] > 0? ";AF=$t[13]" : '';
 372         my $valid = ($t[12] eq 'unknown')? '' : ";valid=$t[12]";
 373         my $info = "molType=$t[10];class=$t[11]$valid$af";
 374         print join("\t", $t[1], $pos, $t[4], $ref, join(",", @alt), 0, '.', $info), "\n";
 375   }
 376 }
 377
 378 sub hapmap2vcf {
 379   die("Usage: vcfutils.pl <in.ucsc.snp> <in.hapmap>\n") if (@ARGV == 0);
 380   my $fn = shift(@ARGV);
 381   # parse UCSC SNP
 382   warn("Parsing UCSC SNPs...\n");
 383   my ($fh, %map);
 384   open($fh, ($fn =~ /\.gz$/)? "gzip -dc $fn |" : $fn) || die;
 385   while (<$fh>) {
 386         my @t = split;
 387         next if ($t[3] - $t[2] != 1); # not SNP
 388         @{$map{$t[4]}} = @t[1,3,7];
 389   }
 390   close($fh);
 391   # write VCF
 392   warn("Writing VCF...\n");
 393   print "##fileformat=VCFv4.0\n";
 394   while (<>) {
 395         my @t = split;
 396         if ($t[0] eq 'rs#') { # the first line
 397           print join("\t", "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT", @t[11..$#t]), "\n";
 398         } else {
 399           next unless ($map{$t[0]});
 400           next if (length($t[1]) != 3); # skip non-SNPs
 401           my $a = \@{$map{$t[0]}};
 402           my $ref = $a->[2];
 403           my @u = split('/', $t[1]);
 404           if ($u[1] eq $ref) {
 405                 $u[1] = $u[0]; $u[0] = $ref;
 406           } elsif ($u[0] ne $ref) { next; }
 407           my $alt = $u[1];
 408           my %w;
 409           $w{$u[0]} = 0; $w{$u[1]} = 1;
 410           my @s = (@$a[0,1], $t[0], $ref, $alt, 0, '.', '.', 'GT');
 411           my $is_tri = 0;
 412           for (@t[11..$#t]) {
 413                 if ($_ eq 'NN') {
 414                   push(@s, './.');
 415                 } else {
 416                   my @a = ($w{substr($_,0,1)}, $w{substr($_,1,1)});
 417                   if (!defined($a[0]) || !defined($a[1])) {
 418                         $is_tri = 1;
 419                         last;
 420                   }
 421                   push(@s, "$a[0]/$a[1]");
 422                 }
 423           }
 424           next if ($is_tri);
 425           print join("\t", @s), "\n";
 426         }
 427   }
 428 }
 429
 430 sub usage {
 431   die(qq/
 432 Usage:   vcfutils.pl <command> [<arguments>]\n
 433 Command: subsam       get a subset of samples
 434          listsam      list the samples
 435          fillac       fill the allele count field
 436          qstats       SNP stats stratified by QUAL
 437          varFilter    filtering short variants
 438          hapmap2vcf   convert the hapmap format to VCF
 439          ucscsnp2vcf  convert UCSC SNP SQL dump to VCF
 440 \n/);
 441 }