#!/usr/bin/env perl

use strict;
use warnings;
use Getopt::Long;
use autodie qw(:all);
use Data::Dumper;
use Bio::SeqIO;
use Bio::Perl;
use Bio::Seq;
use List::Util qw(any);
use Bio::AlignIO;
use File::Slurp qw(read_file);
use Cwd;

#use File::Compare::compare_text;
#------------------------------------------------------------------------
#Definitions
#
#my $pathtotargetgenomeIn = '/scratch/franziska/genomes/worms';	

my $ExceS_A_2 = "/scratch/franziska/ExceS_A/ExceS_A_2.pl";

my $ExceS_A_3 = "/scratch/franziska/ExceS_A/ExceS_A_3.pl";

my $space = "/scratch/franziska/genomes/";


my $genome_query;
my $genome_target;
my $mode;
my $hmmModels;
my $querysequenceIn;



my $distance_first = 1500000; 	#Can be defined by the user in dependence of the intron size expected for the first (few) exons.   
my $distance_last = 1500000;

my $core = 8;	

GetOptions("genomeQuery=s" => \$genome_query, "genomeTarget=s" => \$genome_target, "querysequence=s" => \$querysequenceIn, "mode=s" => \$mode, "hmmModels=s" => \$hmmModels);

#------------------------------------------------------------------------
#Main

my $ID =  `grep ">" $querysequenceIn`;
my @ID = split /\s/, $ID;
$ID = $ID[0];
$ID =~ s/>//;
$ID =~ s/_/./;

mkdir "$ID/";
chdir "$ID/";
my $dir = getcwd();
mkdir "SearchQuery/";
mkdir "SearchTarget/";
mkdir "TargetSequences/";
mkdir "HomologousExons/";

`cp $querysequenceIn SearchQuery/${ID}.fa`;


my @buildquery = "formatdb -i ${genome_query} -p F -o";
unless (-e "${genome_query}.nhr" && -e "${genome_query}.nin" && -e "${genome_query}.nsd" && -e "${genome_query}.nsi" && -e "${genome_query}.nsq"){
	system (@buildquery) == 0 or die "Formatdb with query genome failed. $!\n";
}

my @call = "blastall -p tblastn -d $genome_query -i $querysequenceIn -a $core -m 8 -G 12 -E 1 -o SearchQuery/${ID}.blastout";
system (@call) == 0 or die "blastall against query genome failed. $!\n";

my $locus = `sort -k11,11g SearchQuery/${ID}.blastout| awk '{if (\$3>90) print \$2}'| head -n1`;
chomp $locus;

print Dumper ($locus, $ID);

print Dumper ("perl $ExceS_A_3 --genome $genome_query --locus $locus --ID $ID --querysequence $querysequenceIn");
system("perl $ExceS_A_3 --genome $genome_query --locus $locus --ID $ID --querysequence $querysequenceIn");
`cp SearchQuery/${ID}.inf HomologousExons/paralog_exon_alignment.fa`;

my $contig;
my $strand;
my $first;
my $last;

if ($mode eq "alignment"){

    chdir $space or die "Couldn't access $space\n";
    system "awk '{print \$1}' $genome_target > tmp1";
    unless (-e "${space}${ID}_ORF.aas" && -e $genome_target) {
		my $ret = system "getorf -minsize 9 -reverse Y -stdout Y -sequence tmp1 -outseq stdout|sed 's/ - /-/g'| awk '{print \$1, \$3, \$2}' > tmp2";
		if ($ret != 0){
		    die;
		};
		open (TMP, "< tmp2") or die "Couldn't open tmp2, $!";
		open (DEST, ">> ${space}${ID}_ORF.aas") or die "Couldn't open ${space}${ID}_ORF.aas, $!";
		while (<TMP>){
		    chomp;
		    s/\(//;
		    s/_(\d+)\s\s\[/ $1\(FORWARD\)\[/;
		    s/\sREVERSE/\(REVERSE\)/;
		    s/_(\d+)\(REVERSE\) / $1\(REVERSE\)/;
		    print DEST "$_\n";
		}
		close TMP;
		close DEST;
    };
	chdir $dir;
    mkdir "hmmModels/";
 	`cp -r $hmmModels/* hmmModels/.`;

    my @files = glob "hmmModels/*";
    print "Progress of hmm search against the target genome:\n";
    foreach my $file(@files){
		$_ = $file;
		s/hmmModels\///;
		s/\.model//;
		$file = $_;
		print "$file\n";
		system "hmmsearch --cpu $core --domtblout SearchTarget/${file}_domtableout -E 0.01 -o SearchTarget/${file}_out hmmModels/${file}.model ${space}${ID}_ORF.aas >> SearchTarget/Target_vs_query.hmmsearchout";
	};

	my $hmm = "SearchTarget/hmmsearch_tableout";
    system "cat SearchTarget/*_domtableout| sed \'s/-\\([0-9]\\+\\)]/_\\1]/g\' > $hmm";

	if (-e $hmm){
	    open (IN, "< $hmm") or die "Couldn't open $hmm, $!";
	    open (HITS, "> SearchTarget/hmmsearch_tableout_hits");
	    while (<IN>){
			unless ($_ =~ /#/) {
			    chomp;
			    print HITS "$_\n";
			}
	    }
	    close HITS;
	    close IN;
	}
	$contig = `sort -k7,7g $hmm | awk '{if (\$8>90) print \$1}'| head -n1`;
	chomp $contig;	
	my $startfirst = `sort -k7,7g $hmm | awk '{if (\$8>90) print \$_}'|sort -k23,23g | awk '{print \$23}| head -n1`;
	my $startlast = `sort -k7,7g $hmm | awk '{if (\$8>90) print \$_}'|sort -k23,23g | awk '{print \$23}| tail -n1`;
	chomp $startfirst;
	chomp $startlast;
	if ($startfirst =~ m/FORWARD/){
		$strand = "+";
	}
	else {
		$strand = "-";
	}	

	if ($strand eq "+"){		
		my @last = split /\[/, $startlast;
		my @first = split /\[/, $startfirst;
		$first = (split /\_/, $first[1])[0];
		$last = (split /\_/, $last[1])[1];
		$last =~ s/\]//;
		$first -= $distance_first;
		$last += $distance_last;		
	}
	elsif($strand eq "-"){	
		my @last = split /\[/, $startfirst;
		my @first = split /\[/, $startlast;
		$first = (split /\_/, $first[1])[1];
		$first =~ s/\]//;
		$last = (split /\_/, $last[1])[0];

		$first -= $distance_first;
		$last += $distance_last;
	};
	if ($first < 0){
		$first = 0;
	}
}

elsif($mode eq "fasta"){
	my @buildtarget = "formatdb -i ${genome_target} -p F -o";
	unless (-e "${genome_target}.nhr" && -e "${genome_target}.nin" && -e "${genome_target}.nsd" && -e "${genome_target}.nsi" && -e "${genome_target}.nsq"){
		system (@buildtarget) == 0 or die "Formatdb with target genome failed. $!. You need a permission to write files in the target genome directory\n";
	}

	my $blast = "SearchTarget/Target_vs_query.blastout";
	my @call = "blastall -p tblastn -d $genome_target -i SearchQuery/${ID}.inf -e 0.01 -a $core -m8 -G 11 -E 1 -C F > $blast";
	system (@call) == 0 or die "blastall against target genome failed. $!\n";

	$contig = `sort -k11,11g SearchTarget/Target_vs_query.blastout| awk '{if (\$3>75) print \$2}'| head -n1`;
	chomp $contig;  

	my $startfirst = `sort -k11,11g SearchTarget/Target_vs_query.blastout| awk '{if (\$3>75) print \$_}' |uniq |sort -k9,9g | awk '{print \$9}' |head -n1`;
	my $startlast = `sort -k11,11g SearchTarget/Target_vs_query.blastout| awk '{if (\$3>75) print \$_}' |uniq |sort -k10,10g | awk '{print \$10}' |head -n1`;

	if ($startfirst < $startlast){
		$strand = "+";
	}
	else {
		$strand = "-";
	}
	
	if ($strand eq "+"){
		$first = `sort -k11,11g SearchTarget/Target_vs_query.blastout| awk '{if (\$3>75) print \$_}' |uniq |sort -k9,9g | awk '{print \$9}' |head -n1`;
		$last = `sort -k11,11g SearchTarget/Target_vs_query.blastout| awk '{if (\$3>75) print \$_}' |uniq |sort -k10,10g | awk '{print \$10}' |tail -n1`;
		$first -= $distance_first;
		$last += $distance_last;		
	}
	elsif($strand eq "-"){		
		$first = `sort -k11,11g SearchTarget/Target_vs_query.blastout| awk '{if (\$3>75) print \$_}' |uniq |sort -k10,10g | awk '{print \$10}' |head -n1`;
		$last = `sort -k11,11g SearchTarget/Target_vs_query.blastout| awk '{if (\$3>75) print \$_}' |uniq |sort -k9,9g | awk '{print \$9}' |tail -n1`;
		$first -= $distance_first;
		$last += $distance_last;
	};
	if ($first < 0){
		$first = 0;
	}

	unless (-s "TargetSequences/${contig}.fa"){  
		`fastacmd -d ${genome_target} -s ${contig} | sed \'s/>lcl|/>/g\' > TargetSequences/${contig}.fa`;
		`formatdb -i TargetSequences/${contig}.fa -p F -o`;	
	}

	my $genome_seqIO = Bio::SeqIO->new(-file   => "TargetSequences/${contig}.fa", -format => "Fasta");
	$a = $genome_seqIO->next_seq;
	$a = $a -> seq;

	if ($last > (length $a)){
		$last = (length $a);
	}
}

print Dumper ($contig);
print Dumper ("perl $ExceS_A_2 --contig $contig --query $ID --strand $strand --start $first --end $last --mode 2");
system("perl $ExceS_A_2 --contig $contig --query $ID --strand $strand --start $first --end $last --mode 2");
