* samtools-0.1.4-9 (r340)

[samtools.git] / misc / samtools.pl

diff --git a/misc/samtools.pl b/misc/samtools.pl
index 497b128d7fc3fd55edfe4f9c8335c94d83e5dd29..c014c52bcf0880d39ba23f567a31396017220142 100755 (executable)
--- a/misc/samtools.pl
+++ b/misc/samtools.pl
@@ -6,7 +6,7 @@ use strict;
 use warnings;
 use Getopt::Std;
 
-my $version = '0.3.1';
+my $version = '0.3.2 (r321)';
 &usage if (@ARGV < 1);
 
 my $command = shift(@ARGV);
@@ -36,12 +36,13 @@ sub showALEN {
 #
 
 sub varFilter {
-  my %opts = (d=>3, D=>100, l=>30, Q=>25, G=>25, s=>100, w=>10, W=>10, N=>2, p=>undef);
+  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('pd:D:l:Q:w:W:N:G:', \%opts);
   die(qq/
 Usage:   samtools.pl varFilter [options] <in.cns-pileup>
 
-Options: -Q INT    minimum RMS mapping quality [$opts{Q}]
+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}]
 
@@ -49,12 +50,12 @@ Options: -Q INT    minimum RMS mapping quality [$opts{Q}]
          -w INT    SNP within INT bp around a gap to be filtered [$opts{w}]
 
          -W INT    window size for filtering dense SNPs [$opts{W}]
-         -N INT    max number of SNPs in a window [$opts{W}]
+         -N INT    max number of SNPs in a window [$opts{N}]
 
          -l INT    window size for filtering adjacent gaps [$opts{l}]
 
          -p        print filtered variants
-/) if (@ARGV == 0 && -t STDIN);
+\n/) if (@ARGV == 0 && -t STDIN);
 
   # calculate the window size
   my ($ol, $ow, $oW) = ($opts{l}, $opts{w}, $opts{W});
@@ -72,9 +73,7 @@ Options: -Q INT    minimum RMS mapping quality [$opts{Q}]
        }
        my ($flt, $score) = (0, -1);
        # first a simple filter
-       if ($t[6] < $opts{Q}) {
-         $flt = 1;
-       } elsif ($t[7] < $opts{d}) {
+       if ($t[7] < $opts{d}) {
          $flt = 2;
        } elsif ($t[7] > $opts{D}) {
          $flt = 3;
@@ -82,6 +81,7 @@ Options: -Q INT    minimum RMS mapping quality [$opts{Q}]
        # site dependent filters
        if ($flt == 0) {
          if ($t[2] eq '*') { # an indel
+               $flt = 1 if ($t[6] < $opts{q});
                # filtering SNPs
                if ($t[5] >= $opts{G}) {
                  for my $x (@staging) {
@@ -103,6 +103,7 @@ Options: -Q INT    minimum RMS mapping quality [$opts{Q}]
                  }
                }
          } else { # a SNP
+               $flt = 1 if ($t[6] < $opts{Q});
                # check adjacent SNPs
                my $k = 1;
                for my $x (@staging) {
@@ -213,6 +214,30 @@ sub p2q_print_str {
   }
 }
 
+#
+# varStats
+#
+
+sub varStats {
+  my %opts = (d=>'', c=>5);
+  getopts('d:c:', \%opts);
+  die("Usage: samtools.pl varStats [-d dbSNP.snp] [-c $opts{c}] <in.plp.snp>\n") if (@ARGV == 0 && -t STDIN);
+  my (@cnt, %hash);
+  my $col = $opts{c} - 1;
+  while (<>) {
+       my @t = split;
+       if ($t[2] eq '*') {
+       } else {
+         my $q = $t[$col];
+         $q = 99 if ($q > 99);
+         $q = int($q/10);
+         my $is_het = ($t[3] =~ /^[ACGT]$/)? 0 : 1;
+         ++$cnt[$q][$is_het];
+         $hash{$t[0],$t[1]} = $q;
+       }
+  }
+}
+
 #
 # Usage
 #