#!/usr/bin/perl -w

# Author: lh3

use strict;
use warnings;
use Getopt::Std;

my $version = '0.1.1';
&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, s=>100);
  getopts('D:m:rs:', \%opts); # -s for scaling factor in score calculation

  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});
	# calculate indel score
	my $score = $t[5];
	$score += $opts{s} * $t[10] if ($t[8] ne '*');
	$score += $opts{s} * $t[11] if ($t[9] ne '*');
	@$curr = ($t[0], $t[1], $score, $_);
	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/);
}
