#!/usr/bin/perl -w
use strict;

my $p; #contador proteines
my $q; #path de la query
my $g; #path del genoma
my $blast; 
my $content;
my $contentfinal;
my $blastcut9;
my $blastcut10;
my $inici;
my $iniciprov;
my $iniciprov1;
my $iniciprov2;
my $numscaffold;
my $sizescaffold;
my $scaffold;
my $resta;
my $finalprov;
my $length;


$p=1;

while ($p < 54){ #tenim 53 proteïnes en total
    $q= "/home/u137265/Query/$p.fa"; #seqüència d'aminoàcids de la proteïna en Zebrafish
    $g= "/mnt/NFS_UPF/soft/genomes/2019/Datnioides_undecimradiatus/genome.fa"; #genoma del Datnioides undecimradiatus
    my $gi= "/mnt/NFS_UPF/soft/genomes/2019/Datnioides_undecimradiatus/genome.index"; #variable del índex dels scaffolds del nostre genoma
    my  $length = 200000; #la longitud del subseq que anem a agafar, 100000 per baix de la posició inicial del primer hit i 100000 per dalt
    my $cadena = "/home/u137265/$p/tipuscadena.fa"; #variable que ens dirà si la cadena és forward o reverse

#amb system faré que l'ordinador em faci programes. 
    system ("sed s/@//g $q > temp; mv temp $q -f"); #canvia caracters com @ per res
    system ("sed s/\%//g $q > temp; mv temp $q -f"); #canvia caracters com % per res
    system ("sed s/\#//g $q > temp; mv temp $q -f"); #canvia caracters com # per res
    system ("sed s/U/X/g $q > temp; mv temp $q -f"); #canvia caràcters com U, que és la selenocisteïna, per una X
    
    system ("tblastn -query $q -db $g -outfmt 6 -evalue 0.0001 -out /home/u137265/$p/blastn$p"); #farà el tblastn de la seqüència de la proteïna amb el genoma del peix
    system ("cat /home/u137265/$p/blastn$p"); #ens mostrarà per pantalla el blast per a que apuntem si hi ha valors allunyats i altres scaffolds
	

#############fastafetch##################

    my $scaffold = `cat /home/u137265/$p/blastn$p | head -1 | awk {'print \$2 }'` ; #esta variable contindrà el nom de la scaffold elegida, la del e-value més menut, i la seqüència d'eixa scaffold
    chomp ($scaffold);
    system ("fastafetch $g $gi $scaffold  > /home/u137265/$p/scaffold$p");

  
    
#############FASTASUBSEQ#################

  #fer que no es passi al restar 

    $iniciprov1 = `cat /home/u137265/$p/blastn$p | head -1 | awk {'print \$(NF-3) }'`; #agafarem la primera posició del tblastn, independentment de si és forward o reverse, del scaffold elegit
    chomp ($iniciprov1);

    $iniciprov2 = `cat /home/u137265/$p/blastn$p | head -1 | awk {'print \$(NF-2) }'` ; #agafarem la segona posició del tblastn, independentment de si és forward o reverse, del scaffold elegit
    chomp ($iniciprov2);

    
    if ($iniciprov1 < $iniciprov2){ #comprovarem quina de les dues posicions és més menuda per saber el sentit de la cadena
	$iniciprov = $iniciprov1;
	print "la cadena és forward";
	my $str = "la cadena és forward" ;
         open(FH,'>', $cadena) or die $!;
	 print FH $str;
 	 close (FH);
    }
    else {$iniciprov = $iniciprov2;
	  print "la cadena és reverse \n";
	my $str2 = "la cadena és reverse";
	 open(FH,'>', $cadena) or die $!;
	 print FH $str2;
 	 close (FH);
}

    print "l'scaffold comença a: $iniciprov \n";
    
    #açò ho fem per assegurar-se que la mesura de 100000 nucleòtids per baix de la posició d'inici és possible, és a dir, que la posició d'inici no siga menor a 100000 i per tant ens done error

    $inici= $iniciprov - 100000;
    if ($inici >= 0){
	$inici = $inici}
    else {$inici = 0};
    print "l'inici del scaffold és $inici \n";

    #length 
    system ("grep $scaffold /mnt/NFS_UPF/soft/genomes/2019/Datnioides_undecimradiatus/genome.lengths > /home/u137265/$p/lengths");
    print "la length de l'scaffold $scaffold es: ";
    system ("more /home/u137265/$p/lengths");

    $sizescaffold= `more /home/u137265/$p/lengths | awk {'print \$1}'`;
    chomp ($sizescaffold);

    #per veure si entra, he de fer una variable inici +200000 per veure si em passo de la grandària de la scaffold
    
    if ($sizescaffold > ($inici + 200000)){
	$length = 200000;
	}
    else {$length = ($sizescaffold - $inici)} #hem de tenir en compte que si li dic al length = sizescaffold li hem de restar les posicions d'inici, ja que sino se sortira del scaffold
	
#portem a terme tots els programes posteriors:
    
  
  system ("fastasubseq /home/u137265/$p/scaffold$p $inici  $length > /home/u137265/$p/subseq$p");

  
  
#############EXONERATE###################

  
  system ("exonerate -m p2g --exhaustive yes --showtargetgff -q $q -t /home/u137265/$p/subseq$p | egrep -w exon > /home/u137265/$p/prot$p.exonerate.gff");
  

  
############FASTASEQFROMGFF##############

  
  system ("fastaseqfromGFF.pl /home/u137265/$p/subseq$p  /home/u137265/$p/prot$p.exonerate.gff > /home/u137265/$p/prot$p.exonerate.fa");

  

############FASTATRANSLATE##############

  
  system ("fastatranslate -f /home/u137265/$p/prot$p.exonerate.fa -F 1 | sed 's/\*/X/g' > /home/u137265/$p/translate$p.fa");

    
#################TCOFEE#################

  system ("t_coffee /home/u137265/$p/translate$p.fa $q");
  
$p= $p + 1; #passem a la següent proteïna

}

    
$blast = "blast\n";
system ("more blast")
