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