#! /usr/bin/perl 
# delete the INDELs that convered by another one.

# author: X.Gan

use strict;
use Getopt::Std;
use List::Util qw[min max];

my %opts =( S=>undef, t=>10, p=>"");
getopts('St:p:', \%opts);
my $isSam =$opts{S} ;
my $thrsh= $opts{t};
my $pref=$opts{p};

die "Usage: get_unalign.pl bamfile > out.file\n" if (@ARGV <1); 

my $bamfile =shift;
my $oldchid;
my $seqstr;

my $k;
if ( $isSam) {
	open(BAM_F, $bamfile) or die "Can't find the file";
} else {
	my $cmd="samtools view $bamfile | ";
	open(BAM_F, $cmd) or die "Can't find the file";
}

my @sticky;
my %ct;
while ( <BAM_F> ) 
{
    my $secline = $_;
    my ($contig_id2, $oritNum, $chid2, $startI2, $cig2, $mystrq2);	
    my $rub;
    my $cig;
    ($contig_id2, $oritNum, $chid2, $startI2, $rub, $cig2,$rub, $rub, $rub, $mystrq2) = split(/\t/, $secline);
    my $orit2= ($oritNum ==0 ) ? "+": "-";

    if ( exists  $ct{$contig_id2} ) {
	$ct{$contig_id2} = [ @{$ct{$contig_id2}}, $orit2, $cig2 ];
    } else {
	$ct{$contig_id2} = [ $orit2, $mystrq2, $orit2, $cig2 ];
    }
}

foreach my $ccid (keys %ct) {
    my $fi=$ct{$ccid};

    my ($orit, $mystrq);
    $orit=shift @$fi;
    $mystrq=shift @$fi;
    my @rang_taken;
    while (@$fi >= 1 ) {
	my ($orit2,$cig2);
	$orit2=shift @$fi;	
	$cig2=shift @$fi;
	my @tk=&allel_taken($cig2, $thrsh, $orit2 eq $orit);
	
	if (@rang_taken) {
	    for (@tk) {
		&push_sortedrange(\@rang_taken, $_);
	    }
	} else {
		push(@rang_taken,@tk);
	}
    }

    my $end=length($mystrq);
    my $start=-1;  #last occupied index
    my $k=0;
    for my $rg (@rang_taken) {
	if ($rg ->[0]- $start -1 >= $thrsh ) {
	    my $ssq= substr($mystrq, $start+1, $rg->[0]-$start-1);
	    print ">$pref"."$ccid"."_".$k++, "\n";
	    print $ssq, "\n";
	}    
	$start=$rg->[1]-1;
    }
    if ($end -$start-1 > $thrsh ) {
	my $ssq= substr($mystrq, $start+1, $end-$start-1);
	print ">$pref"."$ccid"."_".$k++, "\n";
	print $ssq, "\n";
    }
}

print STDERR "$k entries were produced \n" if ($k>0);
exit 1;

sub push_sortedrange {
    my ($targ, $newg)=@_;
    my $i;
    for ($i=0; $i<@$targ; $i++) {
	next if  ( $targ->[$i]->[1] < $newg->[0] );
	
	if ($targ->[$i]->[0] > $newg->[1] ) {
	    &arrayInsertAfterPosition($targ, $i-1, $newg);
	    return;
	} else {
	    my $k=$i;
	    my $TSIZE=int @$targ;
	    while ($i < $TSIZE && ($targ->[$i]->[1] >= $newg->[0] || $newg->[1] >= $targ->[$i]->[0]) ) {
		@$newg = (min($newg->[0], $targ->[$i]->[0]), max($newg->[1],$targ->[$i]->[1]));
		$i++;
	    }
	    @{$targ->[$i-1]}=@$newg;

	    for (my $tt=$k; $tt<$i-1; $tt++) {
		delete ( $targ->[$tt] );  
	    }
	    return;
	}

	return;
    }
    push(@$targ, $newg);
}

sub arrayInsertAfterPosition
{
    my ($inArray, $inPosition, $inElement) = @_;
    my @res         = ();
    my @after       = ();
    my $arrayLength = int @{$inArray};

    if ($inPosition < 0) { @after = @{$inArray}; }
    else {
	if ($inPosition >= $arrayLength)    { $inPosition = $arrayLength - 1; }
	if ($inPosition < $arrayLength - 1) { @after = @{$inArray}[($inPosition+1)..($arrayLength-1)]; }
    }

    push (@res, @{$inArray}[0..$inPosition],
              $inElement,
	  @after);

    @$inArray= @res;
}



# or_id : 0 processing from the end; 1 processing from the start
sub allel_taken {   #for a $ol long loci, how many allels has been aligned to
    my ($ci, $thr, $or_id) =@_;
    my @range;
    my $start=0;
    my $cur=0;
    if ($or_id == 0 ) {    #check the end;
	$ci =~ s/(\d+)[HS]$//;
	$start=$cur=$1 if ($1);
	while ( $ci ) {
	    return @range if ($ci !~ /\d+[MNDIPHS]$/);
	    if  ( $ci =~ s/(\d+)[M]$// ) {
		$cur += $1;
	    } elsif  ( $ci =~ s/(\d+)[ND]$// ) {
		;
	    } elsif ( $ci =~ s/(\d+)[HSIP]$// ) {
	        my $nch=$1;
	        if ($nch >= $thr ) {
		    if ($cur >$start) { 	
			push(@range, [$start, $cur-1]);	
 		    }	
		    $cur += $1;
		    $start=$cur;
		} else {
		    $cur += $1;
		}
	   }   
	}
    } elsif ($or_id == 1) {
	$ci =~ s/^(\d+)[HS]//;
	$start=$cur=$1 if ($1);
	while ( $ci ) {
	    return @range if ($ci !~ /^\d+[MNDIPHS]/);
	    if  ( $ci =~ s/^(\d+)[M]// ) {
		$cur += $1;
	    } elsif  ( $ci =~ s/^(\d+)[ND]// ) {
		;
	    } elsif ( $ci =~ s/^(\d+)[HSIP]// ) {
	        my $nch=$1;
	        if ($nch >= $thr ) {
		    if ($cur >$start) { 	
			push(@range, [$start, $cur-1]);	
 		    }	
		    $cur += $1;
		    $start=$cur;
		} else {
		    $cur += $1;
		}
	   }   
	}
    }
    if ($cur >$start) { 	
	push(@range, [$start, $cur-1]);	
    }	

    return @range;
}

sub get_range {
	my ($a) = @_;
	my $mat=0;
	my $ext=0;
	my $sho=0;
	my $clip=0;

	$clip += $1 while ( $a =~ /(\d+)[HS]/g ) ;  #### I didn't take into account \d+S, be careful.
	$mat += $1 while ( $a =~ /(\d+)M/g ) ;  
	$ext += $1 while ( $a =~ /(\d+)[ND]/g ) ;
        $sho += $1 while ( $a =~ /(\d+)[IP]/g ) ;

	return ($mat +$sho, $mat + $ext);  ##first, # of aligned seq; second, # of loci taken in reference
}

