Mention ‘SAMtools’ in libbam-dev's description, to make it easier to find with apt...

[samtools.git] / misc / samtools.pl

diff --git a/misc/samtools.pl b/misc/samtools.pl
index 9f48b8f7586899a9b9456b3fdc09fb8a28faa3d3..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}));
@@ -104,7 +104,6 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
     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;
@@ -474,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
 #