#!/usr/bin/perl -w

#Aquest codi correspon al del FindingSeCys.pl, creat per Mar Costa, Anna González, Clara Mayayo, Mirna Muntal i Isshak Mrabet. Si teniu qualsevol pregunta, no dubteu a enviar-nos un correu a mar.costa01@estudiant.upf.edu o clara.mayayo01@estudiant.upf.edu.

#El codi permet fer alineaments entre una seqüència determinada i l'homòloga que es troba en el genoma que s'està estudiant, que en aquest cas és el de Balaenoptera acutorostrata. Per a fer-ho amb qualsevol genoma, s'ha de modificar el path que ens porta al directory on hi és aquest

#El programa s'executa a través d'un bash, en el qual es realitzen tots els exports necessaris i crida directament aquest programa.

use strict;

my $id; #nom proteïna

my $header; #capçalera del fitxer de la proteïna

my $seq; #seqüència de la proteïna

my @vseq; #vector on s'emmagatzema la seqüència

my $l; #llargada de la seqüència (vector)

my $i = 0; #posicions del vector i comptador llargada

#Aïllament de la seqüència de la proteïna i canvi d'U a X: separa la capçalera i la seqüència del fitxer de la proteïna. Aquesta última l'emmagatzema en un vector que a l'analitzar les diferents posicions, substituirà per una X tant les Us com %, # i @.

while (<STDIN>){

if (/>(\S+)/){

$id = $1;

print "La proteïna que s'està analitzant és $id\n";

}

if ($_=~/^>(.*)/){

$header = $1;

$seq = '';

}

else {

@vseq = split(//,$_);

$l = scalar(@vseq);

for ($i = 0; $i < $l; $i++) {

if ($vseq[$i] eq "U"){

$vseq[$i] = "X";

}

if ($vseq[$i] eq "%"){

$vseq[$i] = "X";

}

if ($vseq[$i] eq "#"){

$vseq[$i] = "X";

}

if ($vseq[$i] eq "@"){

$vseq[$i] = "X";

}

}

$seq= join("",@vseq);

}

}

print "Comprova que no hi hagi cap U a la seqüència: $seq\n";

#Emmagatzema la seqüència en un fitxer

my $seq_fitxer = 'seq_fitxer.fa';

open (FITXER, ">", $seq_fitxer) || die "No es pot crear el fitxer\n";

print FITXER "$seq";

close (FITXER);

my $proteina = 'proteina.fa';

open (PROT, ">", $proteina) || die "No es pot crear el fitxer\n";

print PROT ">$id\n";

print PROT "$seq";

close (PROT);

#Execució del BLAST: es posa el genoma en una variable per tal de poder treballar-hi més fàcilment. Es fa còrrer el programa BLAST simplificat en -m8 per tal d'obtenir els resultats en forma tabular. Agafarem els scaffolds amb un valor de Evalue més petit que 0.0001.

my $genome = "/cursos/BI/genomes/vertebrates/2014/Balaenoptera_acutorostrata/genome.fa";

my $query_file = $id;

my $blast = "";

print "Fent el BLAST...\n";

$blast = "blastall -p tblastn -i 'seq_fitxer.fa' -d '$genome' -o $id.blast -m8 -e 0.0001";

system ($blast);

print "BLAST finalitzat!!\n";

#Indexant genoma: només en el cas que aquest no hagi estat indexat anteriorment.

print "Indexant el genoma...\n";

if (-e "genome.index"){

print "L'índex ja existeix\n";

}

else {

system ("fastaindex '$genome' genome.index");

print "Index creat!\n";

}

#Del document creat en el BLAST, s'extreuen els paràmetres d'interès: nom del scaffold, posició d'inici i posició final. S'analitzen tots els scaffolds possibles, però en el cas que aquest ja s'hagi estudiat anteriorment, no es tornarà a repetir.

my $scaffold = 1; #variable on s'emmagatzema el nom del scaffold

my $oldscaf = 0; #s'emmagatzema el scaffold analitzat prèviament

my @blast; #vector amb els resultats del blast

my $start_pos; #posició d'inici del hit

my $end_pos; #posició final del hit

my $c = 1; #comptador dels diferents scaffolds

my $entrada = "$id.blast";

open (SCAF, '<', $entrada);

while (<SCAF>){

@blast = split("\t",$_);

$scaffold = $blast[1];

$start_pos = $blast[8];

$end_pos = $blast[9];

if ($oldscaf eq $scaffold){

next;

}

print "$scaffold\n";

#Fastafetch: s'agafa la regió corresponent al scaffold del genoma de l'espècie on s'ha trobat el hit.

print "Realitzant el fastafetch...\n";

system ("fastafetch '$genome' genome.index '$scaffold' > '$id.$c.scaffold.fa' ");

#Es modifiquen les posicions d'inici i final depenent de si el hit està en la cadena positiva o negativa. En cada extrem, s'allarga la regió d'interès 100.000 nucleòtids.

my $length;

my $start;

my $end;

if ($start_pos < $end_pos){

$start = $start_pos - 100000;

$end = $end_pos + 100000;

$length = $end - $start;

}

else {

$end = $start_pos + 100000;

$start = $end_pos - 100000;

$length = $end - $start;

}

print "$start_pos\n";

print "$end_pos\n";

print "$length\n";

#Fastasubseq: s'extreu la part d'interès del scaffold a partir del punt d'inici i la llargada prèviament establerts.

print "Fent fastasubseq..\n";

system ("fastasubseq '$id.$c.scaffold.fa' '$start' '$length' > '$id.$c.subseq.fa'");

#Exonerate: s'extreuen els exons de la seqüència anterior. S'utilitza el programa Exonerate en mode 'exhaustive', per tal d'obtenir millors prediccions.

print "Fent l'exonerate..\n";

system ("exonerate -m p2g --showtargetgff --exhaustive yes -q 'proteina.fa' -t '$id.$c.subseq.fa' > '$id.$c.exonerate.fa' ");

system ("cat '$id.$c.exonerate.fa' | egrep -w exon > '$id.$c.exons.exonerate.gff' ");

system ("fastaseqfromGFF.pl '$id.$c.subseq.fa' '$id.$c.exons.exonerate.gff' > $id.$c.dna.fa");

#Fastatranslate: es tradueix la regió del DNA d'interès a aminoàcids, resultant en una predicció de la proteïna en l'espècie Balaenoptera acutorostrata.

print "Traduint...\n";

system ("fastatranslate -f '$id.$c.dna.fa' -F 1 > '$id.$c.predicted.fa'");

print "$id";

#Abans de realitzar el T_coffee es necessita canviar l'* del fitxer per una X, sinó el programa no ho podrà llegir. Per tal de fer això, s'emmagatzema la seqüència d'aminoàcids que s'ha predit en un vector, el qual s'analitzen totes les seves posicions i es substitueix per una X en el cas de trobar un *. Un cop fet el canvi, es redirecciona a un nou fitxer.

my $input = "$id.$c.predicted.fa";

my $output = "$id.$c.predicted2.fa";

my @vseq1;

my $t = 0;

my $d = 0;

my $seq1;

open (INPUT, '<', $input) or die 'ERROR: No es pot obrir!';

open (OUTPUT, '>', $output) or die 'ERROR: No es pot crear!';

while (<INPUT>){

@vseq1 = split(//,$_);

$t = scalar(@vseq1);

for ($d = 0; $d < $t; $d++) {

if ($vseq1[$d] eq "*"){

$vseq1[$d] = "X";

}

}

$seq1= join("",@vseq1);

print OUTPUT $seq1;

}

close (INPUT);

close (OUTPUT);

#Tcoffee: compara la proteïna inicial amb la que s'ha predit.

system ("t_coffee 'proteina.fa' '$id.$c.predicted2.fa' > '$id.$c.tcoffee.fa'");

$oldscaf = $scaffold;

$c = $c + 1;

}

close (SCAF);

#Es crea una carpeta per cada proteïna on es copien tots els fitxers que s'han creat durant el seu anàlisi (es poden moure perquè totes tenen el mateix prefix: el nom de la proteïna).

system("mkdir -p ./$id");

system("mv '$id'* './$id'");

#Durant l'anàlisi es creen alguns fitxers intermediaris que s'han d'eliminar de cara al processament de la següent proteïna.

system ("rm seq_fitxer.fa");

system ("rm proteina*");