#!/usr/bin/perl -w
### Aquest programa ha estat creat per Jan Huertas, Èrika Iglesias, Pau Jané, Pau Jurado i Teresa Puig
### Aquest script analitza si una query donada es troba dins del genoma de Myotis Davidii
### Per contactar, envieu un mail a jan.huertas01@estudiant.upf.edu

use strict;

print "Benvingut al fantabuloso programa d'anotacio per Myotis davidii. Premi enter per continuar\n";
<>;

#condicions prèvies:
#Necessites tenir la query creada, i les U canviades per X per tal que el programa funcioni correctament. A més, la query ha d'estar a la mateixa carpeta que el programa.
my $siono;
print "Ja has creat una query i canviat les U per X (SI o NO)?\n";
$siono = <STDIN>;
chomp $siono;
if ("$siono" eq "NO"){ 
    print "Necessites la query per poder-ho analitzar\n";
    exit
}
my $query;
print "introdueix el nom de la query (sense terminació .fa):";
$query = <STDIN>;
chomp $query;

#blastall
#Aquest programa crearà un blast simplificat (m 8), dins d'una carpeta amb el nom de la query, que donarà només el nom dels scaffolds amb un e-value menor a 10-4 i la posició d'inici i acabament d'aquest.
#Per agilitzar el procés, si el blast ja existeix no es farà.
if (-e "$query.fa") 
{
   print "la query existeix!\n";
   if (-e "./$query/${query}_blast.fa")
   {
	print "el blast ja existeix!\n"
   }
   else 
   {
	system ("mkdir -p ./$query");
	print "blastejant $query, por favor, mantengase a la espera\n";
	system ("blastall -p tblastn -i $query.fa -d /cursos/BI/genomes/project_2014/Myotis_davidii/genome.fa -o ./$query/${query}_blast.fa -e 0.0001 -m 8");
	print "ara ja tens el fitxer ${query}_blast.fa a la carpeta $query\n";
   }
} 
else 
{   
   print "la query no existeix!\n";
   die;
   
}

#index
#L'index es crea només si encara no s'ha creat. 
if (-e "./sortida.index")
   {
	print "l'index ja existeix!\n";
   }
   else 
   {
	print "creant index...\n";	
	system ("fastaindex /cursos/BI/genomes/project_2014/Myotis_davidii/genome.fa myotis.index");
	print "index creat!\n";
   }

#fastafetch
#Per tal de llegir els noms de cada hit, es crea un document amb només aquests (noms.txt). Es llegirà aquest document i es farà servir cada nom per crear l'arxiu fetch amb només aquell scaffold.
#Si el fetch ja existeix, no es crearà.
system ("cut ./$query/${query}_blast.fa -s -f 2 > ./$query/noms.txt");
my $i = 0;
open (NOMSBLAST,"<./$query/noms.txt");
my @v = <NOMSBLAST> ;
close (NOMSBLAST);
for (@v)
{
  chomp;
  unless (-e "./$query/fetch/${_}.fa")
     { 
        system ("mkdir -p ./$query/fetch");
	system ("fastafetch /cursos/BI/genomes/project_2014/Myotis_davidii/genome.fa ./sortida.index '$_' > './$query/fetch/${_}.fa' ");
  	print "$_ creat!\n";
  	$i = $i + 1;
     }
}

print "$i arxius creats!\n";
print "Fent el subseq...\n";

#fastasubseq
#En aquest cas llegim l'arxiu generat per blast i agafem els valors de les columnes 8 i 9, corresponents a posició d'inici i final del hit.
my %hits;
my $a;
my $start_pos;
my $stop_pos;
my $start;
my $length;
system ("mkdir -p ./$query/subseq");
open (BLASTALL, "<./$query/${query}_blast.fa");
while(<BLASTALL>) 
{
   
  # llegeixo linia amb split
  chomp;
  next if ($_ =~ m/^#/);
  my @blast = split("\t",$_);
  # Aquest pas estableix les posicions d'inici i final, considerant si el sentit es positiu o negatiu.
  if ($blast[8] < $blast[9]) {
	$start_pos = $blast[8];
	$stop_pos = $blast[9];
    }
    else {
	$start_pos = $blast[9];
	$stop_pos = $blast[8];
    }

  if (exists($hits{$blast[1]})) 
  {
    #actualitzo les posicions start i stop
    if ($hits{$blast[1]}{"START"} > $start_pos)
	{
	$hits{$blast[1]}{"START"} = $start_pos;
	}
    if ($hits{$blast[1]}{"STOP"} < $stop_pos)
	{
	$hits{$blast[1]}{"STOP"} = $stop_pos;
	}
  }
  else 
  {
    # inicialitzes les posicions d'inici i final
    $hits{$blast[1]}{"START"} = $start_pos;
    $hits{$blast[1]}{"STOP"} = $stop_pos;
  }
#el programa estableix automàticament un marge de +/- 50000 per cada hit. Pot ser que doni un error si el scaffold és més petit. En aquest cas, cal parar el programa i fer l'anàlisi manual.
$start = $start_pos - 50000;
$length = $stop_pos - $start_pos + 100000;
system ("fastasubseq './$query/fetch/${blast[1]}.fa' $start $length > './$query/subseq/${blast[1]}_subseq.fa' ");
# crea els mateixos fitxers molts cops perquè a cada volta varia els start i length
}
close (BLASTALL);
print "Ja tens els arxius subseq creats!\n";
print "Fent l'exonerate...\n";

#exonerate
#primer creem totes les carpetes necessàries.
system ("mkdir -p ./$query/exonerate");
system ("mkdir -p ./$query/exonerate/gff");
system ("mkdir -p ./$query/exonerate/cDNAgff");
system ("mkdir -p ./$query/exonerate/cDNAfa");
system ("mkdir -p ./$query/exonerate/aa_gen");
system ("mkdir -p ./$query/exonerate/transcripts");
system ("mkdir -p ./$query/genewise");
system ("mkdir -p ./$query/tcoffee");
system ("mkdir -p ./$query/tcoffee/exonerate");

#seguidament executem totes les comandes de l'exonerate per cadascun dels arxius de la carpeta subseq
opendir (SUBSEQDIR, "./$query/subseq/");
while (my $file = readdir(SUBSEQDIR))
  {
	next if ($file =~ m/^\./);	
	chomp $file;	
	my $protein_file = $file;
	$protein_file =~ s/\_subseq.fa//g;	
	chomp $protein_file;
	system ("exonerate -m p2g --showtargetgff -q './$query.fa' -t './$query/subseq/$file' > './$query/exonerate/gff/$protein_file.gff' ");
	system ("egrep -w exon './$query/exonerate/gff/${protein_file}.gff' > './$query/exonerate/cDNAgff/${protein_file}_cDNA.gff' ");
	system ("fastaseqfromGFF.pl './$query/subseq/$file' './$query/exonerate/cDNAgff/${protein_file}_cDNA.gff' > './$query/exonerate/cDNAfa/${protein_file}_cDNA.fa' ");
	system ("fastatranslate './$query/exonerate/cDNAfa/${protein_file}_cDNA.fa' > './$query/exonerate/aa_gen/${protein_file}_aa_gen.mfa' ");	
  }
print "exonerate fet! Prem ENTER per passar al T-coffee\n";
<>;

#tcoffee
#donat que el fastatranslate ens dóna un arxiu amb els 6 possibles translates (segons el reading frame), per tal d'optimitzar el procés, hem decidit fer només el transcrit bo. El programa ensenya en pantalla els 6 translates i comparant manualment amb la query l'usuari ha de seleccionar el que és més semblant a la query
my $j = 1;
opendir (TRANSLATE, "./$query/exonerate/aa_gen/");
while (my $fail = readdir(TRANSLATE))
  {
	next if ($fail =~ m/^\./);	
	chomp $fail;
	system ("cat './$query/exonerate/aa_gen/$fail' ");
	my $protein_fail = $fail;
	$protein_fail =~ s/\_aa_gen.mfa//g;
	chomp $protein_fail;
	print "compara $query amb aquests 6 alignments i escriu quina s'assembla més (-3,-2,-1,1,2,3):";
	my $transcript = <STDIN>;	
	chomp $transcript;
	if ($transcript eq "-3")
	{
	system ("fastarevcomp './$query/exonerate/cDNAfa/${protein_fail}_cDNA.fa' > './$query/exonerate/transcripts/${query}_revcomp.fa' ");
	system ("fastatranslate -F 3 './$query/exonerate/transcripts/${query}_revcomp.fa'  > './$query/exonerate/transcripts/${query}_${j}.fa' ");
	system("t_coffee $query.fa './$query/exonerate/transcripts/${query}_${j}.fa' > ./$query/tcoffee/exonerate/${query}_${j}_tcoffee.fa ");
	}	
	if ($transcript eq "-2")
	{
	system ("fastarevcomp './$query/exonerate/cDNAfa/${protein_fail}_cDNA.fa' > './$query/exonerate/transcripts/${query}_revcomp.fa' ");
	system ("fastatranslate -F 2 './$query/exonerate/transcripts/${query}_revcomp.fa'  > './$query/exonerate/transcripts/${query}_${j}.fa' ");
	system("t_coffee $query.fa './$query/exonerate/transcripts/${query}_${j}.fa' > ./$query/tcoffee/exonerate/${query}_${j}_tcoffee.fa ");
	}
	if ($transcript eq "-1")
	{
	system ("fastarevcomp './$query/exonerate/cDNAfa/${protein_fail}_cDNA.fa' > './$query/exonerate/transcripts/${query}_revcomp.fa' ");
	system ("fastatranslate -F 1 './$query/exonerate/transcripts/${query}_revcomp.fa'  > './$query/exonerate/transcripts/${query}_${j}.fa' ");
	system("t_coffee $query.fa './$query/exonerate/transcripts/${query}_${j}.fa' > ./$query/tcoffee/exonerate/${query}_${j}_tcoffee.fa ");
	}
	if ($transcript eq "1")
	{
	system ("fastatranslate -F 1 './$query/exonerate/cDNAfa/${protein_fail}_cDNA.fa' > './$query/exonerate/transcripts/${query}_${j}.fa' ");
	system("t_coffee $query.fa './$query/exonerate/transcripts/${query}_${j}.fa' > ./$query/tcoffee/exonerate/${query}_${j}_tcoffee.fa ");
	}
	if ($transcript eq "2")
	{
	system ("fastatranslate -F 2 './$query/exonerate/cDNAfa/${protein_fail}_cDNA.fa' > './$query/exonerate/transcripts/${query}_${j}.fa' ");
	system("t_coffee $query.fa './$query/exonerate/transcripts/${query}_${j}.fa' > ./$query/tcoffee/exonerate/${query}_${j}_tcoffee.fa ");
	}
	if ($transcript eq "3")
	{
	system ("fastatranslate -F 3 './$query/exonerate/cDNAfa/${protein_fail}_cDNA.fa' > './$query/exonerate/transcripts/${query}_${j}.fa' ");
	system("t_coffee $query.fa './$query/exonerate/transcripts/${query}_${j}.fa' > ./$query/tcoffee/exonerate/${query}_${j}_tcoffee.fa ");
	
	}
   $j = $j + 1;	
   print "\n";
   print "tcoffee creat!\n";
   print "\n";
   print "passar al següent T-coffee? (premi ENTER)\n";
   <>;
}
print "Tots els T-coffees creats! gracias por participar! dale ENTER si quieres aun más diversion, y hacer el genewise\n";
<>;

#genewise
#igual que en l'exonerate, es genera el genewise per tots els arxius de la carpeta subseq
opendir (SUBSEQDIR, "./$query/subseq/");
while (my $arxiu = readdir(SUBSEQDIR))
  {
	next if ($arxiu =~ m/^\./);	
	chomp $arxiu;	
	my $protein_arxiu = $arxiu;
	$protein_arxiu =~ s/\_subseq.fa//g;	
	chomp $protein_arxiu;
	system ("genewise -cdna -pep -pretty -gff -both './$query.fa' './$query/subseq/$arxiu' > './$query/genewise/$protein_arxiu.genewise.fa' ");
  }	

print "Genewise creat! si l'exonerate no ha funcionat correctament, usa la sequencia pep del genewise per fer el tcoffee (manualment) un altre cop! \n";
print "\n";
