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