#!/usr/bin/perl -w
use strict;

my $fetch = "";
my $subseq = "";

#Paths:
my $genomepath = '/export/home/u114892/genomareal.fa';
my $genomeindex = '/export/home/u114892/genomareal.index';
my $fetchpath = '/export/home/u114892/fetch';

my $selenoproteina;
my $hit;
my $inputstart;
my $inputend;
my $especie;
my $proteinpath = '/export/home/u114892';
my $tipus;

#Hit joining:
open (HITS, "<$proteinpath/hitsseleccionats");
my $blasttxt = join '', <HITS>;
close HITS;

my $sentidofusion;
my @blastfilas = split ('\n', $blasttxt);
#my $comptadorblast = 1;

foreach my $blastfila (@blastfilas){
my @blastinputfusion = split ('\t', $blastfila);
my $selenoproteinafusion = $blastinputfusion[2];
    my $hitfusion = $blastinputfusion[3];
my $especiefusion = $blastinputfusion[1];
my $inputstartfusion = $blastinputfusion[10];
my $inputendfusion = $blastinputfusion[11];
if ($inputstartfusion < $inputendfusion){
    $sentidofusion = 'forward';
}
else {
    $sentidofusion = 'reverse';
}
 open (LECTURAFUSION,">> $proteinpath/fusion/$selenoproteinafusion.$hitfusion.$especiefusion.$sentidofusion");  
print LECTURAFUSION "$blastfila\n";
    close (LECTURAFUSION);
}


my $maxstart;
my $maxend;

opendir (my $fusiondir, "$proteinpath/fusion");
my @fusiondirs = readdir $fusiondir;
closedir $fusiondir;

my $tipusfusion;
my $selenoproteinafusion;
my $hitfusion;
my $especiefusion;
my @hitrepecosas;
my $hitarchivo; 

print @fusiondirs;
foreach my $hitarchivopath (@fusiondirs){

    open (FILEHITARCHIVO, "<$proteinpath/fusion/$hitarchivopath");
$hitarchivo = join '', <FILEHITARCHIVO>;
close FILEHITARCHIVO;
  
    my @hitsrepes = split ('\n', $hitarchivo);

    my @determinarmaxs = split ('\t', $hitsrepes[0]);
	my $inputstartfusion = $determinarmaxs[10];
	my $inputendfusion = $determinarmaxs[11];
	if ($inputstartfusion < $inputendfusion){
	    $sentidofusion = 'forward';
	}
	else {
	    $sentidofusion = 'reverse';
	}

	if ($sentidofusion eq 'forward'){
	    $maxstart = 1000000000000;
	    $maxend = 0;
	}
	else {
	    $maxstart = 0;
	    $maxend = 100000000000000;
	}

    
    foreach my $hitrepe (@hitsrepes){
	@hitrepecosas = split ('\t', $hitrepe);
	$tipusfusion = $hitrepecosas[0];
    $selenoproteinafusion = $hitrepecosas[2];
    $hitfusion = $hitrepecosas[3];
    $especiefusion = $hitrepecosas[1];
    $inputstartfusion = $hitrepecosas[10];
    $inputendfusion = $hitrepecosas[11];
    if ($inputstartfusion < $inputendfusion){
	$sentidofusion = 'forward';
    }
    else {
	$sentidofusion = 'reverse';
    }
    if ($sentidofusion eq 'forward'){
       if ($maxstart > $inputstartfusion){
	    $maxstart = $inputstartfusion;
	}
       if ($maxend < $inputendfusion){
	   $maxend = $inputendfusion;
       }
    }
    if ($sentidofusion eq 'reverse'){
	if ($maxstart < $inputstartfusion){
	    $maxstart = $inputstartfusion;
	}
	if ($maxend > $inputendfusion){
	    $maxend = $inputendfusion;
	}
    }
    }
    open (LECTURAFUSIONFINAL,">> $proteinpath/hitsfusio");
    print LECTURAFUSIONFINAL "$tipusfusion\t$selenoproteinafusion\t$hitfusion\t$especiefusion\t$maxstart\t$maxend\n";
    close (LECTURAFUSIONFINAL);
}

#Input with joined hits:
open (HITSFUSIO, "<$proteinpath/hitsfusio");
my $hitsfusiotxt = join '', <HITSFUSIO>;
close HITSFUSIO;

my @blastrows = split ('\n', $hitsfusiotxt);

foreach my $blastrow (@blastrows) {

    my @blastinput = split ('\t', $blastrow);

    $selenoproteina = $blastinput[1];
    $hit = $blastinput[2];
    $inputstart = $blastinput[4];
    $inputend = $blastinput[5];
    $especie = $blastinput[3];
    $tipus = $blastinput[0];

    print "$selenoproteina\n";
    print "$hit\n";
    print "$inputstart\n";
    print "$inputend\n";
    print "$especie\n";


#Position elongation for subseq step:
my $start;
my $lenght;
my $subseqpath = '/export/home/u114892/subseq';

if ($inputstart < $inputend) {
    $start = ( $inputstart - 30000);
    $lenght = ( 60000 + ($inputend - $inputstart));
}

else {
    $start = ( $inputend - 30000);
    $lenght = ( 60000 + ($inputstart - $inputend));
}

    if ($start < 0){
	$start = 0;
    }
    print $start;
    print $lenght;

#Fastafetch:
$fetch = "fastafetch $genomepath $genomeindex $hit > $fetchpath/$especie/$selenoproteina.$hit.$inputstart.$inputend.fetch.fa";

    system ("$fetch");

#Changes in elongation to that cases where fatafetch is too short for our elongation:
    open (FETCHFILE, "<$fetchpath/$especie/$selenoproteina.$hit.$inputstart.$inputend.fetch.fa");
    my $fetchfile = join '', <FETCHFILE>;
    close FETCHFILE;

    my $comptadorfetch;

    my @fetchnucleotids = split ('', $fetchfile);

    foreach my $fetchnucleotid (@fetchnucleotids) {

            $comptadorfetch = $comptadorfetch + 1;
        }

    if ($lenght > $comptadorfetch - 250){
	$lenght = $comptadorfetch - 250;
        $start = $inputstart - ($lenght/2);
    }

#Subseq:
$subseq = "fastasubseq $fetchpath/$especie/$selenoproteina.$hit.$inputstart.$inputend.fetch.fa $start $lenght > $subseqpath/$especie/$selenoproteina.$hit.$inputstart.$inputend.subseq.fa";


#Exonerate:
my $exonerate = "";
my $especieproteinpath = "SP$especie";   
my $exoneratepath = '/export/home/u114892/exonerate';

my $numeroespecie;
if ($especie eq "human"){
    $numeroespecie = "1";
} 
if ($especie eq "horse"){
    $numeroespecie = "2";
}

$exonerate = "exonerate -m p2g --showtargetgff -q $proteinpath/$especieproteinpath.u/$numeroespecie-$selenoproteina.fa.u.fa -t $subseqpath/$especie/$selenoproteina.$hit.$inputstart.$inputend.subseq.fa | egrep -w exon > $exoneratepath/$especie.$selenoproteina.$hit.$inputstart.$inputend.exonerate.gff";



#Fastaseq:
my $fastaseq = "";

$fastaseq = "fastaseqfromGFF.pl $subseqpath/$especie/$selenoproteina.$hit.$inputstart.$inputend.subseq.fa $exoneratepath/$especie.$selenoproteina.$hit.$inputstart.$inputend.exonerate.gff > $proteinpath/cDNAs/cDNA.$especie.$selenoproteina.$inputstart.$inputend.$hit.fa"; 

#Fastatranslate:
my $fastatranslate = "";

$fastatranslate = "fastatranslate $proteinpath/cDNAs/cDNA.$especie.$selenoproteina.$inputstart.$inputend.$hit.fa > $proteinpath/proteins/prot.$especie.$selenoproteina.$inputstart.$inputend.$hit.fa";

my $showfastatranslate = "";

$showfastatranslate = "cat $proteinpath/proteins/prot.$especie.$selenoproteina.$inputstart.$inputend.$hit.fa";

system ("$subseq");
system ("$exonerate");
system ("$fastaseq");
system ("$fastatranslate");
system ("$showfastatranslate");

#Division of the fastatranslate sequences to make the tcoffees individually for each obtained sequence:
my $tcoffee = "";

open (FILE, "<$proteinpath/proteins/prot.$especie.$selenoproteina.$inputstart.$inputend.$hit.fa");
my $prottranslate = join '', <FILE>;
close FILE;

my @tfiles = split ('>', $prottranslate);
my $comptador = 0;

foreach my $tfile (@tfiles) {
    open (LECTURA,"> $proteinpath/tfiles/$comptador.tfile.$especie.$selenoproteina.$inputstart.$inputend.$hit.fa");
print LECTURA ">";    
print LECTURA "$tfile";


#Changing the tfiles '*' for an 'U' to refer to selenocysteine:
    system ("sed -e 's/*/U/g' $proteinpath/tfiles/$comptador.tfile.$especie.$selenoproteina.$inputstart.$inputend.$hit.fa > $proteinpath/tfiles.u/$comptador.tfile.$especie.$selenoproteina.$inputstart.$inputend.$hit.fa.u\n");

#T-coffee:
$tcoffee = "t_coffee $proteinpath/$especieproteinpath/$numeroespecie-$selenoproteina.fa $proteinpath/tfiles.u/$comptador.tfile.$especie.$selenoproteina.$inputstart.$inputend.$hit.fa.u -outfile $proteinpath/tcoffees/$especie.$selenoproteina.$hit.$comptador.$inputstart.$inputend.tcoffee";

system ("$tcoffee");


#Removing of the tcoffees with an score smaller than 900:
open (TCOFFEE, "<$proteinpath/tcoffees/$especie.$selenoproteina.$hit.$comptador.$inputstart.$inputend.tcoffee");
my $tcoffile = join '', <TCOFFEE>;
close TCOFFEE;

my @score = split ('SCORE=', $tcoffile);
my @score2 = split (',',$score[1]);

    if ($score2[0] < 900){
 system ("rm -f $proteinpath/tcoffees/$especie.$selenoproteina.$hit.$comptador.$inputstart.$inputend.tcoffee\n");
system ("rm -f $proteinpath/tcoffees/$especie.$selenoproteina.$hit.$comptador.$inputstart.$inputend.tcoffee.html\n");  
}                                                                                          
#$posicioseq = <STDIN>;                                                                                                                                                                                     
#chomp ($posicioseq);                                                                                                                                          
    $comptador = $comptador + 1;
}
}
