From 8d5d72a53eb4c459bcecb7535b2d0594d9194c9c Mon Sep 17 00:00:00 2001
From: Heng Li <lh3@live.co.uk>
Date: Tue, 17 Feb 2009 10:49:55 +0000
Subject: [PATCH]  * move indel_filter.pl to samtools.pl

---
 misc/indel_filter.pl | 42 -----------------------
 misc/samtools.pl     | 81 ++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 81 insertions(+), 42 deletions(-)
 delete mode 100755 misc/indel_filter.pl
 create mode 100755 misc/samtools.pl

diff --git a/misc/indel_filter.pl b/misc/indel_filter.pl
deleted file mode 100755
index d357fcc..0000000
--- a/misc/indel_filter.pl
+++ /dev/null
@@ -1,42 +0,0 @@
-#!/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]);
diff --git a/misc/samtools.pl b/misc/samtools.pl
new file mode 100755
index 0000000..0b77782
--- /dev/null
+++ b/misc/samtools.pl
@@ -0,0 +1,81 @@
+#!/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/);
+}
-- 
2.39.2