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