#!/usr/bin/perl -w

use strict; 
use warnings;

#Aquest programa ha estat creat per l'Elisabet Mae Henderson, la Julia Martinez, la Judit Puigferrer, la Maria Florencia Sacristan, la Maria Touceda. 
#El programa ha estat creat per automatitzar els processos d'indexat del genoma de Propithecus coquereli, per l'extraccio de les regions genomiques que els resultats del BLAST indiquen que hi ha una selenoproteïna i per l'alineament de les selenoproteines de H. sapiens amb les possibles selenoproteines de P. coquereli.
#Universitat Pompeu Fabra, 4t curs, Bioinformatica.

print "Programa d'anotacio de selenoproteines en Propithecus coquereli\n";


###### Introduccio de la variable $genomepath, que marca el pathway fins la carpeta on esta el genoma de Propithecus coquereli ######

my $genomepath = "/cursos/BI/genomes/2015/Propithecus_coquereli";
chomp $genomepath;

########################################################### Fastaindex ##############################################################
#Indexat del genoma de P.coquereli 
#L'index es crea nomes si encara no s'ha creat. 

if(-e "./P_coquereli.genome.index"){
  print "L'index ja existeix\n";
}
else{
  system("fastaindex $genomepath/genome.fa P_coquereli.genome.index");
  print "L'index s'ha creat\n";
}

######################################################### Fastafetch ################################################################
#TABLA = taula amb coordenades de les proteines, longitud de sequencia i scaffold
######################################################### Fastasubseq ###############################################################

open(TABLA, "<", "./P_coquereli/Pcoquereli_Coordtable.txt");

while(<TABLA>) {

    {
      local $/="\r\n";
      chomp;
    }

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

    my $start = $vectorsubseq[2]; 
    my $length = $vectorsubseq[3];
    my $scaffold = $vectorsubseq[4];
    my $protname = $vectorsubseq[0];

    $scaffold=~ s/,/\./;

    print "\n|$scaffold|\n";
  
    print "RUN:\nfastafetch -i ./P_coquereli.genome.index -f ./Propithecus_co quereli/genome.fa -q $scaffold > ./P_coquereli/P_coquereli.FF/$protname.scaffold.fa\n";    

    system("fastafetch -i ./P_coquereli.genome.index -f $genomepath/genome.fa -q $scaffold > ./P_coquereli/P_coquereli.FF/$protname.scaffold.fa");

    print "RUN:\nfastasubseq -f ./P_coquereli/P_coquereli.FF/$protname.scaffold.fa -s $start -l $length > ./P_coquereli/P_coquereli.FSS/U/$protname.U.fa\n";

    system("fastasubseq -f ./P_coquereli/P_coquereli.FF/$protname.scaffold.fa -s $start -l $length > ./P_coquereli/P_coquereli.FSS/U/$protname.U.fa"); 


########################################################## Us per Xs ################################################################

my $ent;
my $sort;
#shift coge el primer elemento del vector
my $path = shift;


print "OPEN TO READ:\n./P_coquereli/P_coquereli.QUERY/$protname.h.aa.fa\n";
open ($ent, "<./P_coquereli/P_coquereli.QUERY/$protname.h.aa.fa") || die "No es pot obrir el fitxer de entrada $protname.h.aa.fa\n";


print "OPEN TO WRITE:\n./P_coquereli/P_coquereli.QUERY/X/$protname.X.aa.fa";
open ($sort, ">./P_coquereli/P_coquereli.QUERY/X/$protname.X.aa.fa");

my $line_count = 0;

while (<$ent>){

  my $line=$_;
  if ($line_count!=0) {
    $line=~s/U/X/;
    $line=~s/@//;
    $line=~s/#//;
    $line=~s/%//;
  }

  print $sort $line;

  $line_count=$line_count+1;
}

print "\nUs canviades per Xs\n";

close ($ent);
close ($sort);


######################################################### Exonerate #################################################################
######################################################### Fastaseqfromgff ###########################################################

  print "RUN:\nexonerate -m p2g --showtargetgff -q P_coquereli/P_coquereli.QUERY/X/$protname.X.aa.fa -t ./P_coquereli/P_coquereli.FSS/U/$protname.U.fa | egrep -w exon > ./P_coquereli/P_coquereli.EX/GFF/$protname.gff\n";

  system ("exonerate -m p2g --showtargetgff -q P_coquereli/P_coquereli.QUERY/X/$protname.X.aa.fa -t ./P_coquereli/P_coquereli.FSS/U/$protname.U.fa | egrep -w exon > ./P_coquereli/P_coquereli.EX/GFF/$protname.gff");

  print "RUNNING:\nperl ./fastaseqfromGFF.pl ./P_coquereli/P_coquereli.FSS/U/$protname.U.fa ./P_coquereli/P_coquereli.EX/GFF/$protname.gff  > ./P_coquereli/P_coquereli.EX/FA/$protname.p.dna.fa\n";
  system ("perl ./fastaseqfromGFF.pl ./P_coquereli/P_coquereli.FSS/U/$protname.U.fa ./P_coquereli/P_coquereli.EX/GFF/$protname.gff  > ./P_coquereli/P_coquereli.EX/FA/$protname.p.dna.fa"); 


  print "RUNNING:\nfastatranslate -F1 ./P_coquereli/P_coquereli.EX/FA/$protname.p.dna.fa > ./P_coquereli/P_coquereli.EX/FA/AA/$protname.p.aa.fa";

  system ("fastatranslate -F 1 ./P_coquereli/P_coquereli.EX/FA/$protname.p.dna.fa > ./P_coquereli/P_coquereli.EX/FA/AA/$protname.p.aa.fa");

######################################################### TCoffee ###################################################################

  system ("t_coffee ./P_coquereli/P_coquereli.QUERY/X/$protname.X.aa.fa ./P_coquereli/P_coquereli.EX/FA/AA/$protname.p.aa.fa > ./P_coquereli/P_coquereli.TCOFFE/$protname.tcoffee.fa");

  unless (-e "$protname.X.aa.aln" ) {
  
    print "There was a problem with $protname\n";

    exit();

  }
}

print "El programa ha estat executat exitosament\n";