10 &usage if (@ARGV < 1);
12 my $command = shift(@ARGV);
13 my %func = (snpFilter=>\&snpFilter, indelFilter=>\&indelFilter, showALEN=>\&showALEN);
15 die("Unknown command \"$command\".\n") if (!defined($func{$command}));
24 die(qq/Usage: samtools.pl showALEN <in.sam>\n/) if (@ARGV == 0 && -t STDIN);
29 s/(\d+)[SMI]/$l+=$1/eg;
30 print join("\t", @t[0..5]), "\t$l\t", join("\t", @t[6..$#t]), "\n";
39 my %opts = (D=>100, m=>10, r=>undef, s=>100);
40 getopts('D:m:rs:', \%opts); # -s for scaling factor in score calculation
43 Usage: samtools.pl indelFilter [options] <in.indel>\n
44 Options: -D INT maximum read depth [$opts{D}]
45 -m INT minimum distance between two adjacent indels [$opts{m}]
46 \n/) if (@ARGV == 0 && -t STDIN);
49 my ($curr, $last) = (\@arr1, \@arr2);
50 my $is_ref = defined($opts{r})? 1 : 0;
53 next if ($t[2] ne '*');
55 next if ($t[3] eq '*/*');
58 next if ($t[7] > $opts{D});
59 # calculate indel score
61 $score += $opts{s} * $t[10] if ($t[8] ne '*');
62 $score += $opts{s} * $t[11] if ($t[9] ne '*');
63 @$curr = ($t[0], $t[1], $score, $_);
65 if (defined $last->[0]) {
66 if ($curr->[0] eq $last->[0] && $last->[1] + $opts{m} > $curr->[1]) {
67 $do_swap = 0 if ($last->[2] > $curr->[2]);
73 my $tmp = $curr; $curr = $last; $last = $tmp;
76 print $last->[3] if (defined $last->[0]);
84 my %opts = (f=>'', Q=>40, d=>3, w=>10, D=>0, N=>2, W=>10, q=>20, s=>50);
85 getopts('fs:w:q:Q:d:D:W:N:', \%opts);
87 Usage: samtools.pl snpFilter [options] <cns2snp.snp>
89 Options: -d INT minimum depth to call a SNP [$opts{d}]
90 -D INT maximum depth, 0 to ignore [$opts{D}]
91 -Q INT required max mapping quality of the reads covering the SNP [$opts{Q}]
92 -q INT minimum SNP quality [$opts{q}]
94 -f FILE filtered samtools indels [null]
95 -s INT minimum samtols indel score [$opts{s}]
96 -w INT SNP within INT bp around an indel to be filtered [$opts{w}]
98 -W INT window size for filtering dense SNPs [$opts{W}]
99 -N INT maximum number of SNPs in a window [$opts{N}]
103 $opts{D} = 100000000 if ($opts{D} == 0);
104 if ($opts{f}) { # filtered samtools indel
106 open($fh, $opts{f}) || die;
109 next if ($t[2] ne '*' || $t[3] eq '*/*' || $t[5] < $opts{s});
110 for (my $x = $t[1] - $skip + 1; $x < $t[1] + $skip; ++$x) {
116 my (@last, $last_chr);
120 next if ($t[2] eq '*' || $hash{$t[0],$t[1]});
121 my $is_good = ($t[7] >= $opts{d} && $t[7] <= $opts{D} && $t[6] >= $opts{Q} && $t[5] >= $opts{q})? 1 : 0;
122 next unless ($is_good); # drop
123 if ($t[0] ne $last_chr) { # a different chr, print
124 map { print $_->{L} if ($_->{F}) } @last;
128 # The following block implemented by Nathans Weeks.
129 push(@last, {L => $_, X => $t[1], F => 1}); # Enqueue current SNP
130 if ($#last == $opts{N}) { # number of SNPs in queue is N+1
131 if ($last[$#last]{X} - $last[0]{X} < $opts{W}) { # if all within window W
132 map {$_->{F} = 0} @last; # all SNPs in the window of size W are "bad"
134 print STDOUT $last[0]{L} if ($last[0]{F}); # print first SNP if good
135 shift @last # dequeue first SNP
138 # print the last few lines if applicable
139 map { print $_->{L} if ($_->{F}) } @last;
144 Usage: samtools.pl <command> [<arguments>]\n
145 Command: indelFilter filter indels generated by `pileup -c'
146 snpFilter filter SNPs generated by `pileup -c'
147 showALEN print alignment length (ALEN) following CIGAR