#!/usr/bin/perl -w
### Aquest script s'empra per estudiar les selenoproteïnes de l'esèpcie Miniopterus natalensis
### Aquest programa ha estat dissenyat per Ariadna Anunciación, Laura Castarlenas, Gemma COmas i Mar FLores 

use strict;
use warnings; 
use Data::Dumper;

######Determinem les dades. L'index del nostre genoma ja esta realitzat mitjançant fastaindex(mnindex.index)

##### El nostra genoma es troba a: /cursos/20428/BI/genomes/2016/Miniopterus_natalensis/genome.fa
##### L'index del genoma es troba a:/cursos/20428/BI/genomes/2016/Miniopterus_natalensis/genome.index

### Tenim totes les proteines descarregades en la carpeta proteinsmb per a la comparació amb microbat
 
system ("ls ./proteinsmb > namesprot.txt"); 

my @filename;
my $i;

$i = 0;

open(INFILE,"namesprot.txt");
while($filename[$i] = <INFILE>){
  chop($filename[$i]);
	$filename[$i]=~ s/.fa//g; # Els arxius es guarden sense l'extensió .fa
  print "$filename[$i]\n";
# Ens serveix per passar a la seguent proteina
  $i++;
}
close (INFILE);
### A partir d'aquí s'agafarà cada proteïna, una a una, per fer la cerca de selenoproteïnes

my $a = 0;
my $filenameX;
my $filenameB;
my $currentdir; ##per indicar el path de la carpeta 
my @hits;
my $anterior = 0;
my $anteriorstart = 0;
my $anteriorend = 0;

#...BLAST...

unless (-e "mn.fa.nin") {
	system ("makeblastdb -in /cursos/20428/BI/genomes/2016/Miniopterus_natalensis/genome.fa -dbtype nucl -out mn.fa");
	}

system ("mkdir -p proteinsmbX");
system ("mkdir -p outputs"); 
### Creem la carpeta outputs per no barrejar els inputs (proteinsmb) amb els outputs (blast i prediccions).


while ($a < scalar(@filename)) {
	open(INFILE,"./proteinsmb/$filename[$a].fa");
       $filenameX = "./proteinsmbX/$filename[$a]x.fa"; 

### Els canvis de notació es guardaran tots en aquesta carpeta (per evitar problemes d'anotació incorrecta)

	open(OUTPUT,">","$filenameX");
	while(<INFILE>){ 
	
		if (/^>/) {
			print OUTPUT $_; 
		}

		elsif ($_ !~/^>/){	# Es mira si hi ha un símbol @,# o %, i s'elimina
			s/U/X/g;
			s/@//g;
			s/#//g;
			s/%//g;

			print OUTPUT $_;	
		}
	}
	close(OUTPUT);

#...tBlastn...
	$currentdir = "./outputs/$filename[$a]";
	system ("mkdir -p $currentdir");

	$filenameB = "$filename[$a].blast"; 

	system ("blastall -p tblastn -i $filenameX -d mn.fa -o ./$currentdir/$filenameB -e 0.05 -m 8"); 

#...Recerca de hits...
	open (BLAST, "./$currentdir/$filenameB");
	while (<BLAST>) {
		chomp ($_); # Agafa el valor d'una línia 
		@hits = split(/\t/,$_); # Separem cada columna de la línia
		my $evalue;
		$evalue = substr($hits[10], 3, 4);
		print "$hits[10], $evalue\n";
		if ($hits[1] eq $anterior) { # Mirem aquells que tinguin el mateix nom que el hit anterior
			if ($hits[8] < $hits[9]){
				if ($hits[8] <= $anteriorstart) {
					$anteriorstart = $hits[8];
				}
				if ($hits[9] >= $anteriorend) {
					$anteriorend = $hits[9];
				}	
			}
				
			else {
				if ($hits[8] >= $anteriorend) {
					$anteriorend = $hits[8];
				}
				if ($hits[9] <= $anteriorstart) {
					$anteriorstart = $hits[9];
				}	
			}			
			$anterior = $hits[1];
		}	
		elsif ($evalue >= 15) { # Si ens trobem en un scaffold nou, només seleccionem els que tenen un e-value alt
			if ($hits[8] < $hits[9]){
				$anteriorstart = $hits[8];
				$anteriorend = $hits[9];
			}
				
			else {
				$anteriorstart = $hits[9];
				$anteriorend = $hits[8];
			}			
			$anterior = $hits[1];
			
		}
	my %hash; # Creem un diccionari on es guardi el nom de l'scaffold juntament amb el número start i end	
	$hash{$anterior}{"start"}=$anteriorstart;
	$hash{$anterior}{"end"}=$anteriorend;
	print Dumper(\%hash);


	foreach my $key (keys %hash){
		print "$key\n"; 
#...FastaFecth...
		system ("mkdir -p $currentdir/hits");
		system ("mkdir -p $currentdir/hits/fastasubseq");
		system ("mkdir -p $currentdir/hits/predictions");

		unless (-e "$currentdir/hits/$key") {
			system ("fastafetch /cursos/20428/BI/genomes/2016/Miniopterus_natalensis/genome.fa /cursos/20428/BI/genomes/2016/Miniopterus_natalensis/genome.index $key > $currentdir/hits/$key.fa");
		}
	}
#...FastaSubseq...

	foreach my $key (keys %hash){
		my $start = $hash{$key}{"start"} - 5000;
				if ($start <= 0) {
				$start = 0;
				}		
		my $end = $hash{$key}{"end"} - $hash{$key}{"start"} + 10000;


		print $hash{$key}{"start"}, "\t", $hash{$key}{"end"}, "\n";

		system ("fastasubseq $currentdir/hits/$key.fa $start $end > $currentdir/hits/fastasubseq/$key.fastasubseq.fa");
	
#...Exonerate...

		system ("exonerate --exhaustive yes -m p2g --showtargetgff -q $filenameX -t $currentdir/hits/fastasubseq/$key.fastasubseq.fa | egrep -w exon > $currentdir/hits/predictions/$key.pred.gff");


#...Predicció de cDNA exonerate...

		system ("fastaseqfromGFF.pl $currentdir/hits/fastasubseq/$key.fastasubseq.fa $currentdir/hits/predictions/$key.pred.gff > $currentdir/hits/predictions/$key.pred.nuc.gff");	

#...Genewise...

##	system ("genewise -pep -pretty -cdna -gff $filenameX $currentdir/hits/fastasubseq/$key.fastasubseq.fa > $currentdir/hits/predictions/$key.genewise.gff");
#Al final només utilitzem la predicció exonerate

#...FastaTranslate...

		system ("fastatranslate -f $currentdir/hits/predictions/$key.pred.nuc.gff -F 1 > $currentdir/hits/predictions/$key.pred.aa.fa");  

#...t_coffee...

		system ("t_coffee $filenameX $currentdir/hits/predictions/$key.pred.aa.fa > $currentdir/hits/predictions/$key.t_coffee");

		}

	}
$a = $a+1; # Passarem a la següent proteïna 

}

