Uploaded samtools_0.1.16-1_amd64.changes.

[samtools.git] / bcftools / vcfutils.pl

diff --git a/bcftools/vcfutils.pl b/bcftools/vcfutils.pl
index bbc479bfffa47aa2b1089b5b6e56414b696be413..6ced66a2893e9b98884774e88607e9f2060af3d1 100755 (executable)
--- a/bcftools/vcfutils.pl
+++ b/bcftools/vcfutils.pl
@@ -14,11 +14,27 @@ sub main {
   my $command = shift(@ARGV);
   my %func = (subsam=>\&subsam, listsam=>\&listsam, fillac=>\&fillac, qstats=>\&qstats, varFilter=>\&varFilter,
                          hapmap2vcf=>\&hapmap2vcf, ucscsnp2vcf=>\&ucscsnp2vcf, filter4vcf=>\&varFilter, ldstats=>\&ldstats,
   my $command = shift(@ARGV);
   my %func = (subsam=>\&subsam, listsam=>\&listsam, fillac=>\&fillac, qstats=>\&qstats, varFilter=>\&varFilter,
                          hapmap2vcf=>\&hapmap2vcf, ucscsnp2vcf=>\&ucscsnp2vcf, filter4vcf=>\&varFilter, ldstats=>\&ldstats,

-                         gapstats=>\&gapstats);
+                         gapstats=>\&gapstats, splitchr=>\&splitchr, vcf2fq=>\&vcf2fq);

   die("Unknown command \"$command\".\n") if (!defined($func{$command}));
   &{$func{$command}};
 }
 
   die("Unknown command \"$command\".\n") if (!defined($func{$command}));
   &{$func{$command}};
 }
 

+sub splitchr {
+  my %opts = (l=>5000000);
+  getopts('l:', \%opts);
+  my $l = $opts{l};
+  die(qq/Usage: vcfutils.pl splitchr [-l $opts{l}] <in.fa.fai>\n/) if (@ARGV == 0 && -t STDIN);
+  while (<>) {
+       my @t = split;
+       my $last = 0;
+       for (my $i = 0; $i < $t[1];) {
+         my $e = ($t[1] - $i) / $l < 1.1? $t[1] : $i + $l;
+         print "$t[0]:".($i+1)."-$e\n";
+         $i = $e;
+       }
+  }
+}
+

 sub subsam {
   die(qq/Usage: vcfutils.pl subsam <in.vcf> [samples]\n/) if (@ARGV == 0);
   my ($fh, %h);
 sub subsam {
   die(qq/Usage: vcfutils.pl subsam <in.vcf> [samples]\n/) if (@ARGV == 0);
   my ($fh, %h);


@@ -199,8 +215,8 @@ Note: This command discards indels. Output: QUAL #non-indel #SNPs #transitions #
 }
 
 sub varFilter {
 }
 
 sub varFilter {

-  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);
-  getopts('pd:D:W:Q:w:a:1:2:3:4:', \%opts);
+  my %opts = (d=>2, D=>10000000, a=>2, W=>10, Q=>10, w=>10, p=>undef, 1=>1e-4, 2=>1e-100, 3=>0, 4=>1e-4, G=>0, S=>1000);
+  getopts('pd:D:W:Q:w:a:1:2:3:4:G:S:', \%opts);

   die(qq/
 Usage:   vcfutils.pl varFilter [options] <in.vcf>
 
   die(qq/
 Usage:   vcfutils.pl varFilter [options] <in.vcf>
 


@@ -230,14 +246,19 @@ Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools.
          print; next;
        }
        next if ($t[4] eq '.'); # skip non-var sites
          print; next;
        }
        next if ($t[4] eq '.'); # skip non-var sites

+    next if ($t[3] eq 'N'); # skip sites with unknown ref ('N')

        # check if the site is a SNP
        # check if the site is a SNP

-       my $is_snp = 1;
+       my $type = 1; # SNP

        if (length($t[3]) > 1) {
        if (length($t[3]) > 1) {

-         $is_snp = 0;
+         $type = 2; # MNP
+         my @s = split(',', $t[4]);
+         for (@s) {
+               $type = 3 if (length != length($t[3]));
+         }

        } else {
          my @s = split(',', $t[4]);
          for (@s) {
        } else {
          my @s = split(',', $t[4]);
          for (@s) {

-               $is_snp = 0 if (length > 1);
+               $type = 3 if (length > 1);

          }
        }
        # clear the out-of-range elements
          }
        }
        # clear the out-of-range elements


@@ -269,36 +290,47 @@ Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools.
        $flt = 1 if ($flt == 0 && $mq >= 0 && $mq < $opts{Q});
        $flt = 7 if ($flt == 0 && /PV4=([^,]+),([^,]+),([^,]+),([^,;\t]+)/
                                 && ($1<$opts{1} || $2<$opts{2} || $3<$opts{3} || $4<$opts{4}));
        $flt = 1 if ($flt == 0 && $mq >= 0 && $mq < $opts{Q});
        $flt = 7 if ($flt == 0 && /PV4=([^,]+),([^,]+),([^,]+),([^,;\t]+)/
                                 && ($1<$opts{1} || $2<$opts{2} || $3<$opts{3} || $4<$opts{4}));

+       $flt = 8 if ($flt == 0 && ((/MXGQ=(\d+)/ && $1 < $opts{G}) || (/MXSP=(\d+)/ && $1 >= $opts{S})));

 
 

-       # site dependent filters
-       my ($rlen, $indel_score) = (0, -1); # $indel_score<0 for SNPs
+       my $score = $t[5] * 100 + $dp_alt;
+       my $rlen = length($t[3]) - 1; # $indel_score<0 for SNPs

        if ($flt == 0) {
        if ($flt == 0) {

-         if (!$is_snp) { # an indel
-               $rlen = length($t[3]) - 1;
-               $indel_score = $t[5] * 100 + $dp_alt;
-               # filtering SNPs
+         if ($type == 3) { # an indel
+               # filtering SNPs and MNPs

                for my $x (@staging) {
                for my $x (@staging) {

-                 next if ($x->[0] >= 0 || $x->[1] || $x->[4] + $x->[2] + $ow < $t[1]);
+                 next if (($x->[0]&3) == 3 || $x->[1] || $x->[4] + $x->[2] + $ow < $t[1]);

                  $x->[1] = 5;
                }
                # check the staging list for indel filtering
                for my $x (@staging) {
                  $x->[1] = 5;
                }
                # check the staging list for indel filtering
                for my $x (@staging) {

-                 next if ($x->[0] < 0 || $x->[1] || $x->[4] + $x->[2] + $ol < $t[1]);
-                 if ($x->[0] < $indel_score) {
+                 next if (($x->[0]&3) != 3 || $x->[1] || $x->[4] + $x->[2] + $ol < $t[1]);
+                 if ($x->[0]>>2 < $score) {

                        $x->[1] = 6;
                  } else {
                        $flt = 6; last;
                  }
                }
                        $x->[1] = 6;
                  } else {
                        $flt = 6; last;
                  }
                }

-         } else { # a SNP
+         } else { # SNP or MNP

                for my $x (@staging) {
                for my $x (@staging) {

-                 next if ($x->[0] < 0 || $x->[1] || $x->[4] + $x->[2] + $ow < $t[1]);
-                 $flt = 5;
+                 next if (($x->[0]&3) != 3 || $x->[4] + $x->[2] + $ow < $t[1]);
+                 if ($x->[4] + length($x->[7]) - 1 == $t[1] && substr($x->[7], -1, 1) eq substr($t[4], 0, 1)
+                         && length($x->[7]) - length($x->[6]) == 1) {
+                       $x->[1] = 5;
+                 } else { $flt = 5; }

                  last;
                }
                  last;
                }

+               # check MNP
+               for my $x (@staging) {
+                 next if (($x->[0]&3) == 3 || $x->[4] + $x->[2] < $t[1]);
+                 if ($x->[0]>>2 < $score) {
+                       $x->[1] = 8;
+                 } else {
+                       $flt = 8; last;
+                 }
+               }

          }
        }
          }
        }

-       push(@staging, [$indel_score, $flt, $rlen, @t]);
+       push(@staging, [$score<<2|$type, $flt, $rlen, @t]);

   }
   # output the last few elements in the staging list
   while (@staging) {
   }
   # output the last few elements in the staging list
   while (@staging) {


@@ -311,7 +343,7 @@ sub varFilter_aux {
   if ($first->[1] == 0) {
        print join("\t", @$first[3 .. @$first-1]), "\n";
   } elsif ($is_print) {
   if ($first->[1] == 0) {
        print join("\t", @$first[3 .. @$first-1]), "\n";
   } elsif ($is_print) {

-       print STDERR join("\t", substr("UQdDaGgP", $first->[1], 1), @$first[3 .. @$first-1]), "\n";
+       print STDERR join("\t", substr("UQdDaGgPMS", $first->[1], 1), @$first[3 .. @$first-1]), "\n";

   }
 }
 
   }
 }
 


@@ -427,6 +459,87 @@ sub hapmap2vcf {
   }
 }
 
   }
 }
 

+sub vcf2fq {
+  my %opts = (d=>3, D=>100000, Q=>10, l=>5);
+  getopts('d:D:Q:l:', \%opts);
+  die(qq/
+Usage:   vcfutils.pl vcf2fq [options] <all-site.vcf>
+
+Options: -d INT    minimum depth          [$opts{d}]
+         -D INT    maximum depth          [$opts{D}]
+         -Q INT    min RMS mapQ           [$opts{Q}]
+         -l INT    INDEL filtering window [$opts{l}]
+\n/) if (@ARGV == 0 && -t STDIN);
+
+  my ($last_chr, $seq, $qual, $last_pos, @gaps);
+  my $_Q = $opts{Q};
+  my $_d = $opts{d};
+  my $_D = $opts{D};
+
+  my %het = (AC=>'M', AG=>'R', AT=>'W', CA=>'M', CG=>'S', CT=>'Y',
+                        GA=>'R', GC=>'S', GT=>'K', TA=>'W', TC=>'Y', TG=>'K');
+
+  $last_chr = '';
+  while (<>) {
+       next if (/^#/);
+       my @t = split;
+       if ($last_chr ne $t[0]) {
+         &v2q_post_process($last_chr, \$seq, \$qual, \@gaps, $opts{l}) if ($last_chr);
+         ($last_chr, $last_pos) = ($t[0], 0);
+         $seq = $qual = '';
+         @gaps = ();
+       }
+       die("[vcf2fq] unsorted input\n") if ($t[1] - $last_pos < 0);
+       if ($t[1] - $last_pos > 1) {
+         $seq .= 'n' x ($t[1] - $last_pos - 1);
+         $qual .= '!' x ($t[1] - $last_pos - 1);
+       }
+       if (length($t[3]) == 1 && $t[7] !~ /INDEL/ && $t[4] =~ /^([A-Za-z.])(,[A-Za-z])*$/) { # a SNP or reference
+         my ($ref, $alt) = ($t[3], $1);
+         my ($b, $q);
+         $q = $1 if ($t[7] =~ /FQ=(-?[\d\.]+)/);
+         if ($q < 0) {
+               $_ = $1 if ($t[7] =~ /AF1=([\d\.]+)/);
+               $b = ($_ < .5 || $alt eq '.')? $ref : $alt;
+               $q = -$q;
+         } else {
+               $b = $het{"$ref$alt"};
+               $b ||= 'N';
+         }
+         $b = lc($b);
+         $b = uc($b) if (($t[7] =~ /MQ=(\d+)/ && $1 >= $_Q) && ($t[7] =~ /DP=(\d+)/ && $1 >= $_d && $1 <= $_D));
+         $q = int($q + 33 + .499);
+         $q = chr($q <= 126? $q : 126);
+         $seq .= $b;
+         $qual .= $q;
+       } elsif ($t[4] ne '.') { # an INDEL
+         push(@gaps, [$t[1], length($t[3])]);
+       }
+       $last_pos = $t[1];
+  }
+  &v2q_post_process($last_chr, \$seq, \$qual, \@gaps, $opts{l});
+}
+
+sub v2q_post_process {
+  my ($chr, $seq, $qual, $gaps, $l) = @_;
+  for my $g (@$gaps) {
+       my $beg = $g->[0] > $l? $g->[0] - $l : 0;
+       my $end = $g->[0] + $g->[1] + $l;
+       $end = length($$seq) if ($end > length($$seq));
+       substr($$seq, $beg, $end - $beg) = lc(substr($$seq, $beg, $end - $beg));
+  }
+  print "\@$chr\n"; &v2q_print_str($seq);
+  print "+\n"; &v2q_print_str($qual);
+}
+
+sub v2q_print_str {
+  my ($s) = @_;
+  my $l = length($$s);
+  for (my $i = 0; $i < $l; $i += 60) {
+       print substr($$s, $i, 60), "\n";
+  }
+}
+

 sub usage {
   die(qq/
 Usage:   vcfutils.pl <command> [<arguments>]\n
 sub usage {
   die(qq/
 Usage:   vcfutils.pl <command> [<arguments>]\n


@@ -434,8 +547,14 @@ Command: subsam       get a subset of samples
          listsam      list the samples
          fillac       fill the allele count field
          qstats       SNP stats stratified by QUAL
          listsam      list the samples
          fillac       fill the allele count field
          qstats       SNP stats stratified by QUAL

-         varFilter    filtering short variants
+

          hapmap2vcf   convert the hapmap format to VCF
          ucscsnp2vcf  convert UCSC SNP SQL dump to VCF
          hapmap2vcf   convert the hapmap format to VCF
          ucscsnp2vcf  convert UCSC SNP SQL dump to VCF

+
+         varFilter    filtering short variants (*)
+         vcf2fq       VCF->fastq (**)
+
+Notes: Commands with description endting with (*) may need bcftools
+       specific annotations.

 \n/);
 }
 \n/);
 }