#!/usr/bin/perl -w

use strict;

my $helpStr = qq{
    
    parseGap.pl [gapSize] [contig.fna db] [file.gaps] [file2.gaps]....

	This script find the overlapping region of the input files.

	[contig.fna]: the formatted db of contigs
	[*.gaps]: generated from 'findGap.pl > *.gaps'

};

unless (@ARGV>=3){die "$helpStr";}
my ($gapSize, $contigF, @files) = @ARGV;
my %gaps = (); ## parse and store *.gaps in this hash
my $binDir = "/usr2/people/gtl/bin";
## parsing all *.gaps
foreach (@files){
    open GAP, "$_" or die "Can't open $_: $!\n";
    while (<GAP>){
	next if /^\#/;
	s/\n//;
	my ($taxId, $contig, $l, $start, $end) = split /\t/;
	push @{ $gaps{$contig}{$start."\t".$end} },$taxId;
    }
    close GAP;
}

## now find out which ones overlap

foreach my $contig (keys %gaps){
    &checkOverlap($gapSize, $contig, \%gaps);
}

#map{my $contig = $_; map{print join "\t", $contig, $_, "\n" if @{$gaps{$contig}{$_}} ==scalar(@ARGV);} keys %{$gaps{$contig}} } keys %gaps;

foreach my $contig (keys %gaps){
    #print $contig."\n";
    my @temp = `$binDir/blast/fastacmd -s "$contig" -d $contigF`;
    shift @temp;
    chomp @temp;
    my $seq = join "", @temp;
    foreach my $key (keys %{$gaps{$contig}}){
	if (@{$gaps{$contig}{$key}} ==scalar(@files)){
	    my ( $start, $end )= split /\t/, $key;
	    print ">".$contig." ".$start." ".$end."\n";
	    my $gSeq = substr($seq, $start, ($end-$start+1));
	    print $gSeq."\n";
	}
    }

}


exit;

########################################################


sub checkOverlap($$\%){
    my ($gapSize, $contig, $gapsR) = @_;
    my $count;
    do{
	#print $contig."\n";
	$count = 0;
	my %temp = %{ $gapsR->{$contig} };
	
	foreach my $key (keys %temp){
	   
	    my ($start, $end) = split /\t/, $key;
	    
	    my %uniqueTax = ();
	    
	    foreach my $key2 ( keys %temp ){
		
		my ($keyS, $keyE) = split /\t/, $key2;
		if ( ($keyS<$end) && ($start<$keyE) && ($keyS!=$start) && ($keyE!=$end) ){
		    ## if another gap overlaps with this one
		    #print $key, "\t", $key2, "\n";
		    map{$uniqueTax{$_}++;} @{$temp{$key}};
		    map{$uniqueTax{$_}++;} @{$temp{$key2}};
		    
		    if (keys %uniqueTax == scalar(@files)){
			## if all genomes have one of these two gaps
			## create a new gap of the overlapping region
			$count = 1;
			my ($newS, $newE);
			if ($start >= $keyS){
			    $newS = $start;
			}
			else{
			    $newS = $keyS;
			} 
			if ($end <= $keyE){
			    $newE = $end;
			}
			else{
			    $newE = $keyE;
			} 
			delete $gapsR->{$contig}{$key};
			delete $gapsR->{$contig}{$key2};
			if (($newE-$newS)>=$gapSize){
			#print $newS,"\t", $newE, "new fragment \n";
			    @{$gapsR->{$contig}{$newS."\t".$newE}} = keys %uniqueTax unless $gapsR->{$contig}{$newS."\t".$newE};
			}
		    }
		}
		
	    }
	}
	
    }while($count);
    return;
    
}
