#!/usr/bin/perl -w
# Aquest programa ha estat millorat per Lorena Galera, Carlota Grifols, Meritxell Guardiola i Jaume Minoves
# Aquest programa analitza si una query donada es troba dins del genoma de Lethenteron_camtschaticum
# Per qualsevol consulta envieu un mail a carlota.grifols01@estudiant.upf.edu

use strict;

print "Benvingut al programa d'anotacio per Lethenteron_camtschaticum. Premi enter per continuar\n";
<>;




# Requisits del programa: La query ha d'estar ubicada a la mateixa carpeta que el programa.
# No és necessari canviar les U per X doncs el programa accesori escrit en bash prova.sh ja s'encarrega de fer-ho.




# INTRODUCCIÓ DE LA QUERY
my $query;
print "Introdueixi el nom de la query (sense extensió .fa):";
$query = <STDIN>;
chomp $query;




# AUTOMATITZACIÓ DE blastall
# Es genera un blast simplificat (indicat a la comanda blastall mitjançant -m 8), dins d'una carpeta amb el nom de la query, que informarà només dels scaffolds amb un e-value menor a 10-4, de la posició d'inici i de la posició d'acabament d'aquest.
# Si el blast ja ha estat generat en recerques anteriors no es tornarà a generar de novo.
if (-e "$query.fa") 
{
   print "Query trobada!\n";
   if (-e "./$query/${query}_blast.fa")
   {
	print "Blast ja ha estat generat en recerques anteriors!\n"
   }
   else 
   {
	system ("mkdir -p ./$query");
	print "Executant blastall per a $query, esperi si us plau.\n";
	system ("blastall -p tblastn -i $query.fa -d /cursos/BI/genomes/vertebrates/2014/Lethenteron_camtschaticum/genome.fa -o ./$query/${query}_blast.fa -e 0.0001 -m 8");
	print "L'arxiu ${query}_blast.fa ha estat creat i emmagatzemat a la carpeta $query\n";
   }
} 
else 
{   
   print "La query no ha estat trobada!\n";
   die;
   
}




# AUTOMATITZACIÓ DE fastaindex (Ordenar el genoma de Lethenteron camtschaticum pels seus fragments).
# Si l'index ja ha estat generat en recerques anteriors no es tornarà a generar de novo. 
if (-e "./genome.index")
   {
	print "La indexació del genoma ja ha estat generat en recerques anteriors!\n";
   }
   else 
   {
	print "Creant index...\n";	
	system ("fastaindex /cursos/BI/genomes/vertebrates/2014/Lethenteron_camtschaticum/genome.fa genome.index");
	print "La indexació del genoma ha estat generada!\n";
   }




# AUTOMATITZACIÓ DE fastafetch (Seleccionar l'scaffold on es troba el hit. Es crearan tants scaffolds com hits s'hagin trobat.)
# L'arxiu noms.txt contindrà els noms de cada hit. Es llegirà l'arxiu noms.txt i es farà servir cada nom per crear l'arxiu fetch l'scaffold corresponent.
# Si l'arxiu fetch ja existeix, no serà creat de novo.
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/vertebrates/2014/Lethenteron_camtschaticum/genome.fa ./genome.index '$_' > './$query/fetch/${_}.fa' ");
  	print "$_ creat!\n";
  	$i = $i + 1;
     }
}

print "$i arxius creats!\n";
print "Fent el subseq...\n";




# AUTOMATITZACIÓ DE fastasubseq (Acotar regió de l'sacaffold on es troba el hit)

# Es llegueix l'arxiu generat per blast i es capturen els valors de les columnes 8 i 9 corresponents a la 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>) 
{
   
  # Es captura una línia amb la comanda split.
  chomp;
  next if ($_ =~ m/^#/);
  my @blast = split("\t",$_);
  
  # Establiment de les posicions d'inici i final. Es té en compte 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]})) 
  {
    # Actualització de 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 
  {
    # Inicialització de les posicions d'inici i final.
    $hits{$blast[1]}{"START"} = $start_pos;
    $hits{$blast[1]}{"STOP"} = $stop_pos;
  }

# S'estableix automàticament un marge de +/- 50000 per cada hit. Si el scaffold es més petit es produiría error.
$start = $start_pos - 50000;
if ($start < 0) {
	$start = 0;
}
$length = $stop_pos - $start_pos + 100000;

# Correcció del fet que si l'scaffold és més petit es produeixi error. Es procedeix a fer que si $start + $length és més gran que la longitut de l'scaffold, aleshores farem que $length adquireixi el valor de (longitut de l'scaffold menys $stop). 

# S'emmagatzema al fitxer long0.fa el nombre de caràcters de l' scaffold sencer
system ("wc './$query/fetch/${blast[1]}.fa' -m > './numcarseq.fa' ");
system ("cut -d ' ' -f1 './numcarseq.fa' > './long0.fa' ");

# S'emmagatzema al fitxer long1.fa el nombre de caràcters del nom de l'scaffold
system ("grep '>' './$query/fetch/${blast[1]}.fa' > './nomseq.fa' ");
system ("wc './nomseq.fa' -m > './numcarnom.fa' ");
system ("cut -d ' ' -f1 './numcarnom.fa' > './long1.fa' ");

# Transport de la informació continguda als fitxers long0.fa i long1.fa cap al Perl

# Càlcul de la longitut de l' scaffold amb nom inclós i emmagatzematge a $longseq0[0]
my @longseq0;
open (SCAFFOLDSENCER,"<./long0.fa");
my @longseq0 = <SCAFFOLDSENCER> ;
close (SCAFFOLDSENCER);
$longseq0[0] = $longseq0[0] - 1;

# Càlcul de la longitut del nom de l' scaffold i emmagatzematge a $longseq1[0]
my @longseq1;
open (NOMSCAFFOLD,"<./long1.fa");
my @longseq1 = <NOMSCAFFOLD> ;
close (NOMSCAFFOLD);
$longseq1[0] = $longseq1[0] - 1;

# Càlcul de la longitut de l'scaffold sense el nom
my $longscaffoldsensenom;
$longscaffoldsensenom = $longseq0[0] - $longseq1[0];

# Es procedeix a esborrar tots els arxius intermitjos generats
system ("rm './numcarseq.fa' "); system ("rm './long0.fa' ");
system ("rm './nomseq.fa' "); system ("rm './numcarnom.fa' "); system ("rm './long1.fa' ");

# Es comparar $start + $length amb la longitut de l'scaffold sense el nom, i si és més gran farem que $length adquireixi el valor de (longitut de l'scaffold menys $stop).
if (($start + $length) > $longscaffoldsensenom) {
	$length = $longscaffoldsensenom - $stop_pos;
}


# Crida de la comanda Fastasubseq
system ("fastasubseq './$query/fetch/${blast[1]}.fa' $start $length > './$query/subseq/${blast[1]}_subseq.fa' ");

# Es generen els mateixos fitxers molts cops perquè a cada volta es varia l'start i el length.
}
close (BLASTALL);
print "Els arxius subseq han estat creats!\n";
print "Fent l'exonerate...\n";




# AUTOMATITZACIÓ DE exonerate

# Prèviament es procedeix a generar 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");

# A continuació es procedeix a executar totes les comandes associades a 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 generat! Premi ENTER per passar al T-coffee\n";
<>;




# AUTOMATITZACIÓ DE tcoffee (Alineament entre la proteïna identificada en el genoma de Lethenteron_camtschaticum i la query).
# Es fa servir només la pauta de lectura +1
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;
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 ");
	
$j = $j + 1;	
   print "tcoffee creat!\n";
   print "\n";
}
print "T-coffee creat! Premi ENTER per fer el genewise\n";
<>;




# AUTOMATITZACIÓ DE genewise (Genewise fa servir un algorisme diferent al que fa servir l'exonerate. Servirà per contrastar la informació obtinguda).
# A l'igual que en l'exonerate, es genera el genewise per cadascun dels 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!\n";
print "\n";




# FI DE PROGRAMA


