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);
  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 usage {
 458   die(qq/
 459 Usage:   vcfutils.pl <command> [<arguments>]\n
 460 Command: subsam       get a subset of samples
 461          listsam      list the samples
 462          fillac       fill the allele count field
 463          qstats       SNP stats stratified by QUAL
 464          varFilter    filtering short variants
 465          hapmap2vcf   convert the hapmap format to VCF
 466          ucscsnp2vcf  convert UCSC SNP SQL dump to VCF
 467 \n/);
 468 }