#!/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';

#Inputnput variables:                                                                                                                                                                
my $selenoproteina;
my $hit;
my $inputstart;
my $inputend;

#Inputs insertion:                                                                                                                                                               

print "Insereix la selenoproteïna:";
$selenoproteina = <STDIN>;
chomp ($selenoproteina);

print "Insereix el hit:";
$hit = <STDIN>;
chomp ($hit);

print "Insereix la posició incial:";
$inputstart = <STDIN>;
chomp ($inputstart);

print "Insereix la posició final:";
$inputend = <STDIN>;
chomp ($inputend);

#Per tenir els hits separats per espècie i alinear les seqüències amb l'espècie a comparar:                                                                                       
my $especie;

print "Espècie (human o horse)?";
$especie = <STDIN>;
chomp ($especie);


$fetch = "fastafetch $genomepath $genomeindex $hit > $fetchpath/$especie/$hit.fetch.fa";


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

if ($inputstart < $inputend) {
    $start = ( $inputstart - 2000);
    $lenght = ( 4000 + ($inputend - $inputstart));
}
else {
    $start = ( $inputend - 2000);
    $lenght = ( 4000 + ($inputstart - $inputend));
}

#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/$hit.fetch.fa $start $lenght > $subseqpath/$especie/$selenoproteina.$hit.subseq.fa";


#Exonerate:                                                                                                                                                                       
my $exonerate = "";
my $especieproteinpath = "SP$especie";
my $proteinpath = '/export/home/u114892';
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.subseq.fa | egrep -w exon > $exoneratepath/$especie.$selenoproteina.$hit.exonerate.gff";



#Fastaseq:                                                                                                                                                                        
my $fastaseq = "";

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

#Fastatranslate:                                                                                                                                                                  
my $fastatranslate = "";

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

my $showfastatranslate = "";

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

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

#Division of the fastatranslate sequences to make the tcoffees individually for each obtained sequence:                                                                                                                                                                      

open (FILE, "<$proteinpath/proteins/prot.$especie.$selenoproteina.$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.$hit.fa");
print LECTURA ">";
print LECTURA "$tfile";
    close (LECTURA);

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

$comptador = $comptador + 1;
}

#Manual selection for each sequence for the T-coffee:
my $posicioseq;
print "De les anteriors, digues quina posició ocupa la seqüència que vols comparar:";
$posicioseq = <STDIN>;
chomp ($posicioseq);

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

system ("$tcoffee");

print "Vols comparar-ne alguna més (si o no)?";
my $sino = <STDIN>;
chomp ($sino);

while ($sino eq 'si' || $sino eq 'sí' || $sino eq 'Si' || $sino eq 'Sí') {
    system ("cat $proteinpath/proteins/prot.$especie.$selenoproteina.$hit.fa");
    print "De les anteriors, digues quina posició ocupa la seqüència que vols comparar:";
    $posicioseq = <STDIN>;
    chomp ($posicioseq);
system ("t_coffee $proteinpath/$especieproteinpath/$numeroespecie-$selenoproteina.fa $proteinpath/tfiles.u/$posicioseq.tfile.$especie.$selenoproteina.$hit.fa.u -outfile $proteinpath/tcoffeesaparte/$posicioseq.tfile.$especie.$selenoproteina.$hit.fa");
    print "Vols comparar-ne alguna més (si o no)?";
    $sino = <STDIN>;
    chomp ($sino);
}
