#!/usr/bin/perl -w

#Aquest programa ha estat desenvolupat per Julia Yemelyanenko, Marc Sans, Elena Rodríguez, Gaudi Vall i Mireia Kun durant la realització del projecte de l'assignatura de Bioinformàtica de quart curs de Biologia Humana de la Universitat Pompeu Fabra, partint de programes desenvolupats per al projecte per estudiants de cursos anteriors. 


use strict;
system ('export PATH=/cursos/BI/bin/ncbiblast/x64/bin:$PATH');
system ('cp /cursos/BI/bin/ncbiblast/.ncbirc ~/');

my $query;
$query = $ARGV[0];
chomp $query;

################################################################################

#Blastall
#Aquest programa crearà un blast simplificat amb tota la informació només dels hits amb un e-value menor a 1e-3 dins de la carpeta Blast.Per agilitzar el procés, si el blast ja existeix no es farà.


if (-e "../inputs/$query.fasta") 
{
   if (-e "../outputs/Blast/${query}.blast")
   {
       print "Ja s'havia fet el blast per aquesta query\n";
   }
   else 
   {
	print "Executant el blast per $query\n";
	system ("blastall -a 4 -e 1e-3 -m 9 -p tblastn -i ../inputs/$query.fasta -d /cursos/BI/genomes/2015/Aotus_nancymaae/genome.fa -o ../outputs/Blast/${query}.blast");
	print "${query}.blast s'ha desat a Documents/BI/outputs/Blast/${query}.blast\n";
   }
} 
else 
{   
   print "Aquesta query no existeix\n";
   die;   
}

################################################################################

#Fastafetch

#En primer lloc creem un fitxer anomenat nomscontig.txt amb el nom de tots els scaffolds del genoma d'Aotus on s'ha trobat un o mes hits amb e-value < 10-3 segons el blastall fet previament.  

system ("cut ../outputs/Blast/${query}.blast -s -f 2 | sort -u > ../outputs/Blast/${query}_nomscontig.txt");

# Es llegira aquest fitxer i es fara servir cada nom per crear un arxiu fetch amb nomes la sequencia de nucleotids corresponent a aquell scaffold.El fastafetch nomes es creara si encara no existeix.

system ("mkdir ../outputs/fastafetch/${query}");
print "Executant el Fastafetch per $query\n";

open (NOMSBLAST, "<../outputs/Blast/${query}_nomscontig.txt");
my @v = <NOMSBLAST> ;
close (NOMSBLAST);
my %contig_length;
for (@v)
{
  chomp;
  my $contig = $_;
  unless (-e "../outputs/fastafetch/${query}/${_}.fa")
     {
        system ("fastafetch /cursos/BI/genomes/2015/Aotus_nancymaae/genome.fa /cursos/BI/genomes/2015/Aotus_nancymaae/genome.index '$_' > '../outputs/fastafetch/${query}/${_}.fa' ");
     }

#A partir del fitxer contig.fastafetch que conte la sequencia de nucleotids de cada scaffold, extraurem la llargada de cadascun i l'emmagatzemarem en una variable $contig_length, que sera necessaria per fer el fastasubseq.


  unless (exists ($contig_length{$contig}))
  {
      my $cl = 0;
      open (FASTA,"<../outputs/fastafetch/${query}/${_}.fa");
      while(<FASTA>) 
      {
	  chomp;
	  next if ($_ =~ m/^>/);
	  my @nucl = split(//, $_);
	  $cl = $cl + scalar(@nucl);   
      }
      $contig_length{$contig} = $cl;
  }
}



###############################################################################


#Fastasubseq

#Llegim l'arxiu generat per blast i agafem els valors de les columnes 8 i 9, corresponents a posicio d'inici i final del hit.
my %hits;
my $start_pos;
my $stop_pos;
my $start;
my $length;
system ("mkdir -p ../outputs/subseq/$query");
open (BLASTALL, "<../outputs/Blast/${query}.blast");
my $i = 1;
while(<BLASTALL>) 
{
 chomp;
 next if ($_ =~ m/^#/);
 my @blast = split("\t",$_);


#Establim les posicions d'inici i final segons si la sequencia es troba en sentit directe o reverse.
 if ($blast[8] < $blast[9]) {
	$start_pos = $blast[8];
	$stop_pos = $blast[9];
 }
 else {
	$start_pos = $blast[9];
	$stop_pos = $blast[8];
 }


#Mirem si ja s'ha trobat un hit en aquest scaffold abans; si ja se n'ha trobat un pero la posicio d'inici era mes gran, ens quedem amb la posicio d'inici del nou hit. I el mateix per la posicio final. 


if (exists($hits{$blast[1]})) 
 {

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;
	}
 }
#Si es el primer hit que es llegeix per aquest scaffold, inicialitzem les posicions d'inici i final
 else 
  {
    $hits{$blast[1]}{"START"} = $start_pos;
    $hits{$blast[1]}{"STOP"} = $stop_pos;
  }

 my $controlstart = $hits{$blast[1]}{"START"};
 my $controlstop = $hits{$blast[1]}{"STOP"};


#Establim un marge de +-50000 per cada hit. Si l'inici del hit es troba a una distància de l'inici de l'scaffold menor de 50000, establim 0 com a posició inicial.

$start = $hits{$blast[1]}{"START"} - 50000;
if ($start < 0) { $start = 0; }
$length = $hits{$blast[1]}{"STOP"} - $hits{$blast[1]}{"START"} + 100000;

system ("fastasubseq '../outputs/fastafetch/$query/${blast[1]}.fa' $start $length > '../outputs/subseq/$query/${blast[1]}_subseq.fa' ");

$i++;

}
close (BLASTALL);
print "Ja s'han creat els arxius subseq per a cada scaffold\n";





################################################################################


# Exonerate
# Creem una carpeta per a cada tipus de fitxer d'output.
system ("mkdir -p ../outputs/exonerate/$query");
system ("mkdir -p ../outputs/exonerate/$query/gff");
system ("mkdir -p ../outputs/exonerate/$query/cDNAgff");
system ("mkdir -p ../outputs/exonerate/$query/cDNAfa");
system ("mkdir -p ../outputs/exonerate/$query/aa_gen");
system ("mkdir ../outputs/tcoffee/$query");

# Executem totes les comandes de l'exonerate per cadascun dels arxius de la carpeta subseq
opendir (SUBSEQDIR, "../outputs/subseq/$query");

	print "Executant l'exonerate per $query\n";

while (my $file = readdir(SUBSEQDIR))
 {
	next if ($file =~ m/^\./);	
	chomp $file;	
	my $protein_file = $file;
	$protein_file =~ s/\_subseq.fa//g;	
	chomp $protein_file;

#Amb exonerate obtenim la sequencia de cDNA incloent nomes els exons (NO hem utilitzat l'opció EXHAUSTIVE)
	system ("exonerate -m p2g --showtargetgff -q '../inputs/$query.fasta' -t '../outputs/subseq/$query/$file' > '../outputs/exonerate/$query/gff/${protein_file}_alignment.gff'");

system (" egrep -w exon '../outputs/exonerate/$query/gff/${protein_file}_alignment.gff' > '../outputs/exonerate/$query/cDNAgff/${protein_file}_cDNA.gff' ");
	
#FastaseqfromGFF
#obtenim la sequencia nucleotidica de cDNA en format fasta a partir de la sequencia nucleotidica del scaffold i  el fitxer cDNA.gff que acabem d'obtenir
	system ("fastaseqfromGFF.pl '../outputs/subseq/$query/${protein_file}_subseq.fa' '../outputs/exonerate/$query/cDNAgff/${protein_file}_cDNA.gff' > '../outputs/exonerate/$query/cDNAfa/${protein_file}_cDNA.fa' ");
	
#Fastatranslate

	print "Executant el Fastatranslate per $query\n";

	system ("fastatranslate -F 1 '../outputs/exonerate/$query/cDNAfa/${protein_file}_cDNA.fa' > '../outputs/exonerate/$query/aa_gen/${protein_file}_aa_gen.fa' ");	

#Tcoffee

	print "Executant el Tcoffee per $query\n";
	system ("t_coffee '../inputs/$query.fasta' '../outputs/exonerate/$query/aa_gen/${protein_file}_aa_gen.fa' > '../outputs/tcoffee/$query/${protein_file}_tcoffee.fa' ");
}
