+++ /dev/null
-#!/usr/bin/perl -w
-
-# Contact: lh3
-
-use strict;
-use warnings;
-use Getopt::Std;
-
-my %opts = (D=>100, m=>10, r=>undef);
-getopts('D:m:r', \%opts);
-
-die(qq/
-Usage: indel_filter.pl [options] <in.indel>\n
-Options: -D INT maximum read depth [$opts{D}]
- -m INT minimum distance between two adjacent indels [$opts{m}]
-\n/) if (@ARGV == 0 && -t STDIN);
-
-my (@arr1, @arr2);
-my ($curr, $last) = (\@arr1, \@arr2);
-my $is_ref = defined($opts{r})? 1 : 0;
-while (<>) {
- my @t = split;
- next if ($t[2] ne '*');
- if (!$is_ref) {
- next if ($t[3] eq '*/*');
- next if ($t[5] == 0);
- }
- next if ($t[7] > $opts{D});
- @$curr = ($t[0], $t[1], $t[5], $_);
- my $do_swap = 1;
- if (defined $last->[0]) {
- if ($curr->[0] eq $last->[0] && $last->[1] + $opts{m} > $curr->[1]) {
- $do_swap = 0 if ($last->[2] > $curr->[2]);
- } else { # then print
- print $last->[3];
- }
- }
- if ($do_swap) {
- my $tmp = $curr; $curr = $last; $last = $tmp;
- }
-}
-print $last->[3] if (defined $last->[0]);
--- /dev/null
+#!/usr/bin/perl -w
+
+# Author: lh3
+
+use strict;
+use warnings;
+use Getopt::Std;
+
+my $version = '0.1.0';
+&usage if (@ARGV < 1);
+
+my $command = shift(@ARGV);
+my %func = (indelFilter=>\&indelFilter, showALEN=>\&showALEN);
+
+die("Unknown command \"$command\".\n") if (!defined($func{$command}));
+&{$func{$command}};
+exit(0);
+
+#
+# showALEN
+#
+
+sub showALEN {
+ die(qq/Usage: samtools.pl showALEN <in.sam>\n/) if (@ARGV == 0 && -t STDIN);
+ while (<>) {
+ my @t = split;
+ my $l = 0;
+ $_ = $t[5];
+ s/(\d+)[SMI]/$l+=$1/eg;
+ print join("\t", @t[0..5]), "\t$l\t", join("\t", @t[6..$#t]), "\n";
+ }
+}
+
+#
+# indelFilter
+#
+
+sub indelFilter {
+ my %opts = (D=>100, m=>10, r=>undef);
+ getopts('D:m:r', \%opts);
+
+ die(qq/
+Usage: samtools.pl indelFilter [options] <in.indel>\n
+Options: -D INT maximum read depth [$opts{D}]
+ -m INT minimum distance between two adjacent indels [$opts{m}]
+\n/) if (@ARGV == 0 && -t STDIN);
+
+ my (@arr1, @arr2);
+ my ($curr, $last) = (\@arr1, \@arr2);
+ my $is_ref = defined($opts{r})? 1 : 0;
+ while (<>) {
+ my @t = split;
+ next if ($t[2] ne '*');
+ if (!$is_ref) {
+ next if ($t[3] eq '*/*');
+ next if ($t[5] == 0);
+ }
+ next if ($t[7] > $opts{D});
+ @$curr = ($t[0], $t[1], $t[5], $_);
+ my $do_swap = 1;
+ if (defined $last->[0]) {
+ if ($curr->[0] eq $last->[0] && $last->[1] + $opts{m} > $curr->[1]) {
+ $do_swap = 0 if ($last->[2] > $curr->[2]);
+ } else { # then print
+ print $last->[3];
+ }
+ }
+ if ($do_swap) {
+ my $tmp = $curr; $curr = $last; $last = $tmp;
+ }
+ }
+ print $last->[3] if (defined $last->[0]);
+}
+
+sub usage {
+ die(qq/
+Usage: samtools.pl <command> [<arguments>]\n
+Command: indelFilter filter indels generated by `pileup -c'
+ showALEN print alignment length (ALEN) following CIGAR
+\n/);
+}