Revert one of my earlier changes - Heng was right, CIGAR P not sensible in a padded...

[samtools.git] / misc / samtools.pl

diff --git a/misc/samtools.pl b/misc/samtools.pl
index e6dd8e8eac12e659ea0ec12e5bb4ca0572252ed1..d03c1c7f148dcca835946e7ccd9ee52f3b6dc00b 100755 (executable)
--- a/misc/samtools.pl
+++ b/misc/samtools.pl
@@ -10,7 +10,7 @@ my $version = '0.3.3';
 &usage if (@ARGV < 1);
 
 my $command = shift(@ARGV);
-my %func = (showALEN=>\&showALEN, pileup2fq=>\&pileup2fq, varFilter=>\&varFilter,
+my %func = (showALEN=>\&showALEN, pileup2fq=>\&pileup2fq, varFilter=>\&varFilter, plp2vcf=>\&plp2vcf,
                        unique=>\&unique, uniqcmp=>\&uniqcmp, sra2hdr=>\&sra2hdr, sam2fq=>\&sam2fq);
 
 die("Unknown command \"$command\".\n") if (!defined($func{$command}));
@@ -46,10 +46,12 @@ sub showALEN {
 # G close to a high-quality indel (SNP only)
 # Q low RMS mapping quality (SNP only)
 # g close to another indel with higher quality (indel only)
+# s low SNP quality (SNP only)
+# i low indel quality (indel only)
 
 sub varFilter {
-  my %opts = (d=>3, D=>100, l=>30, Q=>25, q=>10, G=>25, s=>100, w=>10, W=>10, N=>2, p=>undef);
-  getopts('pq:d:D:l:Q:w:W:N:G:', \%opts);
+  my %opts = (d=>3, D=>100, l=>30, Q=>25, q=>10, G=>25, s=>100, w=>10, W=>10, N=>2, p=>undef, S=>'', i=>'');
+  getopts('pq:d:D:l:Q:w:W:N:G:S:i:', \%opts);
   die(qq/
 Usage:   samtools.pl varFilter [options] <in.cns-pileup>
 
@@ -57,6 +59,8 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
          -q INT    minimum RMS mapping quality for gaps [$opts{q}]
          -d INT    minimum read depth [$opts{d}]
          -D INT    maximum read depth [$opts{D}]
+         -S INT    minimum SNP quality [$opts{S}]
+         -i INT    minimum indel quality [$opts{i}]
 
          -G INT    min indel score for nearby SNP filtering [$opts{G}]
          -w INT    SNP within INT bp around a gap to be filtered [$opts{w}]
@@ -80,7 +84,8 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
        next if (uc($t[2]) eq uc($t[3]) || $t[3] eq '*/*'); # skip non-var sites
        # clear the out-of-range elements
        while (@staging) {
-         last if ($staging[0][2] eq $t[0] && $staging[0][3] + $max_dist >= $t[1]);
+      # Still on the same chromosome and the first element's window still affects this position?  
+         last if ($staging[0][3] eq $t[0] && $staging[0][4] + $staging[0][2] + $max_dist >= $t[1]);
          varFilter_aux(shift(@staging), $opts{p}); # calling a function is a bit slower, not much
        }
        my ($flt, $score) = (0, -1);
@@ -90,14 +95,31 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
        } elsif ($t[7] > $opts{D}) {
          $flt = 3;
        }
+    if ($t[2] eq '*') { # an indel
+        if ($opts{i} && $opts{i}>$t[5]) { $flt = 8; }
+    }
+    elsif ($opts{S} && $opts{S}>$t[5]) { $flt = 7; }    # SNP
+
        # site dependent filters
+    my $len=0;
        if ($flt == 0) {
          if ($t[2] eq '*') { # an indel
+        # If deletion, remember the length of the deletion
+        my ($a,$b) = split(m{/},$t[3]);
+        my $alen = length($a) - 1;
+        my $blen = length($b) - 1;
+        if ( $alen>$blen )
+        {
+            if ( substr($a,0,1) eq '-' ) { $len=$alen; }
+        }
+        elsif ( substr($b,0,1) eq '-' ) { $len=$blen; }
+
                $flt = 1 if ($t[6] < $opts{q});
                # filtering SNPs
                if ($t[5] >= $opts{G}) {
                  for my $x (@staging) {
-                       next if ($x->[0] >= 0 || $x->[3] + $ow < $t[1]);
+            # Is it a SNP and is it outside the SNP filter window?
+                       next if ($x->[0] >= 0 || $x->[4] + $x->[2] + $ow < $t[1]);
                        $x->[1] = 5 if ($x->[1] == 0);
                  }
                }
@@ -107,7 +129,8 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
                $score += $opts{s} * $t[11] if ($t[9] ne '*');
                # check the staging list for indel filtering
                for my $x (@staging) {
-                 next if ($x->[0] < 0 || $x->[3] + $ol < $t[1]);
+          # Is it a SNP and is it outside the gap filter window
+                 next if ($x->[0] < 0 || $x->[4] + $x->[2] + $ol < $t[1]);
                  if ($x->[0] < $score) {
                        $x->[1] = 6;
                  } else {
@@ -119,17 +142,17 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
                # check adjacent SNPs
                my $k = 1;
                for my $x (@staging) {
-                 ++$k if ($x->[0] < 0 && $x->[3] + $oW >= $t[1] && ($x->[1] == 0 || $x->[1] == 4 || $x->[1] == 5));
+                 ++$k if ($x->[0] < 0 && $x->[4] + $x->[2] + $oW >= $t[1] && ($x->[1] == 0 || $x->[1] == 4 || $x->[1] == 5));
                }
                # filtering is necessary
                if ($k > $opts{N}) {
                  $flt = 4;
                  for my $x (@staging) {
-                        $x->[1] = 4 if ($x->[0] < 0 && $x->[3] + $oW >= $t[1] && $x->[1] == 0);
+                        $x->[1] = 4 if ($x->[0] < 0 && $x->[4] + $x->[2] + $oW >= $t[1] && $x->[1] == 0);
                  }
                } else { # then check gap filter
                  for my $x (@staging) {
-                       next if ($x->[0] < 0 || $x->[3] + $ow < $t[1]);
+                       next if ($x->[0] < 0 || $x->[4] + $x->[2] + $ow < $t[1]);
                        if ($x->[0] >= $opts{G}) {
                          $flt = 5; last;
                        }
@@ -137,7 +160,7 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
                }
          }
        }
-       push(@staging, [$score, $flt, @t]);
+       push(@staging, [$score, $flt, $len, @t]);
   }
   # output the last few elements in the staging list
   while (@staging) {
@@ -148,9 +171,9 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
 sub varFilter_aux {
   my ($first, $is_print) = @_;
   if ($first->[1] == 0) {
-       print join("\t", @$first[2 .. @$first-1]), "\n";
+       print join("\t", @$first[3 .. @$first-1]), "\n";
   } elsif ($is_print) {
-       print STDERR join("\t", substr("UQdDWGgX", $first->[1], 1), @$first[2 .. @$first-1]), "\n";
+       print STDERR join("\t", substr("UQdDWGgsiX", $first->[1], 1), @$first[3 .. @$first-1]), "\n";
   }
 }
 
@@ -450,6 +473,44 @@ sub uniqcmp_aux {
   close($fh);
 }
 
+sub plp2vcf {
+  while (<>) {
+       my @t = split;
+       next if ($t[3] eq '*/*');
+       if ($t[2] eq '*') { # indel
+         my @s = split("/", $t[3]);
+         my (@a, @b);
+         my ($ref, $alt);
+         for (@s) {
+               next if ($_ eq '*');
+               if (/^-/) {
+                 push(@a, 'N'.substr($_, 1));
+                 push(@b, 'N');
+               } elsif (/^\+/) {
+                 push(@a, 'N');
+                 push(@b, 'N'.substr($_, 1));
+               }
+         }
+         if ($a[0] && $a[1]) {
+               if (length($a[0]) < length($a[1])) {
+                 $ref = $a[1];
+                 $alt = ($b[0] . ('N' x (length($a[1]) - length($a[0])))) . ",$b[1]";
+               } elsif (length($a[0]) > length($a[1])) {
+                 $ref = $a[0];
+                 $alt = ($b[1] . ('N' x (length($a[0]) - length($a[1])))) . ",$b[0]";
+               } else {
+                 $ref = $a[0];
+                 $alt = ($b[0] eq $b[1])? $b[0] : "$b[0],$b[1]";
+               }
+         } else {
+               $ref = $a[0]; $alt = $b[0];
+         }
+         print join("\t", @t[0,1], '.', $ref, $alt, $t[5], '.', '.'), "\n";
+       } else { # SNP
+       }
+  }
+}
+
 #
 # Usage
 #