#!/usr/bin/perl -w

use strict;
use lib "/usr2/people/khuang/bin/Genomics/Perl/modules";
use Genome;
use Scaffold;
use GenomicsUtils;
use Getopt::Long;
my $help = qq{

    findGap.pl <outDir> <jobId> <unfinished contig fna file> <finished genome taxId> <gapSize> (options)

	<gapSize>: the minimum size of the nonsyntenic region
	(options):
	       -v: verbose

      Description: This script uses Nucmer to align draft genomes
      to a closely related finished genome. The script then parses
      the output from Nucmer to find regions in the unfinished genome
      that have no synteny in the finished genome

};

unless (@ARGV>=5) { die "$help";}

connect();
my ($outDir, $jobId, $contigF, $taxId, $gapSize) = @ARGV;
chdir("$outDir") or die "Can't cd to $outDir";
my $verbose = 0;
GetOptions( 'v' => \$verbose
);
my $bin = "/usr2/people/gtl/bin";
my $blastall = "$bin/blast/blastall";
my $nucmer = "$bin/MUMmer3.08/nucmer";
my $formatdb = "$bin/blast/formatdb";

#### parse the contig file into a hash table and formatdb the file 
my %contigs = ();
my $fastaHeader = "";
open CONTIG, "$contigF" or die "can't open $contigF:$!\n";
while (<CONTIG>){
    if (/>/){
	$fastaHeader = $_;
    }
    else{
	s/\n//;
	$contigs{$fastaHeader} .= $_;
    }
}
close CONTIG;


## Fetch the finished genome by taxonomy ID
my $genome = Genome::new(taxonomyId=>$taxId);
my @scaffolds = $genome->scaffolds(); 


## @scaffolds, a list of scaffold objects (genomic only)
## A finished genome only has one scaffold.
## However this script does handle multiple scaffolds,
## if you really want to compare an unfinished genome to another unfinished genome

open OUT, ">$outDir/$jobId.$taxId.gaps" or die "Can't write to $outDir/$jobId.$taxId.gaps: $!\n";
print OUT "## The reference genome is: ",$genome->name(), "\n";
print OUT "## Find non-syntenic region >= $gapSize bp \n";
print OUT join "\t", "## taxID","contig", "length", "start", "end", "\n";

my @gapFiles = ();
foreach my $scaffold (@scaffolds){
    my $sId = $scaffold->scaffoldId();
    my $seq = $scaffold->sequence();
    my $fnaF = "$jobId.$taxId.$sId.fna";
    open SCAFFOLD, ">$fnaF" or die "Can't write to $fnaF:$!\n";
    print SCAFFOLD ">", $sId, "\n", $seq, "\n";
    close SCAFFOLD;
    
## formatdb
    system("nice $formatdb -i $fnaF -p F -o T");

## nucmer
    system("nice $nucmer -o -p $jobId.$taxId.$sId $fnaF $contigF");

## parse nucmer results, just the coords file is needed
    open COORDS, "$jobId.$taxId.$sId.coords" or die "Can't open file $jobId.$taxId.$sId.coords:$!\n";
    my @coords = <COORDS>;
    close COORDS;
    
    
    foreach (keys %contigs){
	my $fna = $contigs{$_};
	my ($contig, $fnaL) = split / /;
	$contig =~ s/>//;

	#my $gapFile = "$jobId.$taxId.$sId.$contig.$gapSize.gaps";

	#open GAP, ">$gapFile" or die "Can't write to $gapFile:$!\n";

	my @lines = grep /$contig/, @coords;
	my %coords = ();
	foreach (@lines){
	    my @temp = split /\|/;
	    next if (/S1.*E1/)||(!/\|+/);
	    #print;
	    #die;
	    $temp[1] =~ s/^\s+//;
	    print $_."======\n" unless $temp[1];
	    my ($sS, $sE) = split /\s+/, $temp[1];
	    if ($sE>$sS){
		$coords{"$sS"} = $sE;
	    }
	    else{
		$coords{"$sE"} = $sS;
	    }
	    
	}
	
        my ($pStart, $pEnd); ## previous start and end positions
	print "========================\n" if $verbose;
	
	## now coordinates can be sorted by the begin position of the alignment regardless of strands
	foreach my $start (sort{$a<=>$b} keys %coords){

	    my $end = $coords{$start};
	    print join "\t", $contig,$start, $end, "\n" if $verbose;
	    if (defined($pStart)){
		next if ($end<$pEnd);
		if ($start-$pEnd-1 >=$gapSize){
		## the gap starts from $pEnd+1 to $start-1, so the gap length is ($start-1)-($pEnd+1)+1 => $start-$pEnd-1
		    my $gapL = $start-$pEnd-1;
		    my $gapSeq = substr($fna, $pEnd+1, $gapL);
		    #print join " ", $contig, ($start-$pEnd-1), $pEnd+1, $start-1,$gapSeq, "\n";
		    #print join " ", $contig, ($start-$pEnd-1), $pEnd+1, $start-1,"\n";
		
		    open FASTA, ">$jobId.temp.fna";
		    print FASTA ">".$contig." ".($pEnd+1)." ".($start-1)."\n";
		    print FASTA $gapSeq."\n";
		    close FASTA;

		    ## blast against the original fna file, 
		    ## because Nucmer might sometimes be to stingent and will not align regions with over the threshold mismatches
		    my @blast = `nice $blastall -F F -p blastn -i $jobId.temp.fna -d $fnaF -m 8 -v 1 -b 1`;
		    #print $blast[0] if $verbose;
		    ## Top BLASTn hit:
		    my ($q, $s, $id, $aL, $mis, $gap, $qS, $qE, $sS, $sE, $evalue, $score) = split /\t/, $blast[0];
		    #print ($aL/$gapL)."===\n" if $verbose;
		    if (($aL/$gapL>=0.5)&&($aL/$gapL<=0.8)){
			
			#print join " ", $contig, $gapL,  $pEnd+1, $start-1, $aL, $qS, $qE, "\n" if $verbose;

			## query start and stop, scale them back to the contig coords
			if ($qS>$qE){
			    my $swap = $qS;
			    $qS = $qE;
			    $qE = $swap;
			}
			$qS = ($pEnd+1)-1+$qS; 
			$qE = ($pEnd+1)-1+$qE;
		    
			## find un-aligned regions: here we get the end fragments not the gaps within the alignment
			my $fAL = ($qS-1)-($pEnd+1)+1;
			my $fBL = ($start-1)-($qE+1)+1;
			my $fA = substr($fna, $pEnd+1, $fAL);
			my $fB = substr($fna, $qE+1, $fBL);
			#print GAP ">",join " ", $contig."_".$fAL, $pEnd+1, $qS-1, "\n", $fA, "\n" if $fAL >= $gapSize;
			#print GAP ">", join " ", $contig."_".$fBL, $qE+1, $start-1, "\n", $fB, "\n" if $fBL >= $gapSize;
			print OUT join "\t", $taxId, $contig, $fAL, $pEnd+1, $qS-1, "\n" if $fAL >= $gapSize;
			print OUT join "\t", $taxId, $contig, $fBL, $qE+1, $start-1, "\n" if $fBL>=$gapSize;
		    }
		    elsif ($aL/$gapL<0.5){
		        #print GAP ">", join " ", $contig."_".$gapL,  $pEnd+1, $start-1, "\n", $gapSeq, "\n"; 
			print OUT join "\t", $taxId, $contig, $gapL, $pEnd+1, $start-1, "\n";
		    }
		}
	    }
	    $pStart = $start;
	    $pEnd = $end;
		    
	}
	#close GAP; 
	#push @gapFiles, $gapFile if (-s "$gapFile");
## BLASTx the gaps against all proteomes in the database

    }		
}
  
#open BLASTX, ">$taxId.blastx.cmd" or die "Can't write to $taxId.blastx.cmd:$!\n";

#foreach (@gapFiles){
    #print BLASTX "nice $blastall -F F -m 7 -p blastx -e 0.001 -o $_.blastx -i $_ -d all.genomes.faa -v 10 -b 10\n";
    #system("nice $blastall -F F -m 7 -p blastx -e 0.001 -o $_.blastx -i $_ -d all.genomes.faa -v 10 -b 10");   
#}
#close BLASTX;

disconnect();

exit;
