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