#!/usr/bin/perl -w
use strict;
# Abello Laie, Angelats Laura, Armenteros Elena, Ayneto Alba, Balcells Reyes
#per executar ./programa.mexicana.pl ./parseblastmexicana ./multifasta_major.fa



	my ($genoma, $ref_genoma, $parseblast, $genP, $genN, $genPi, $genNi, $MF, $MFi);

	$ref_genoma="/disc8/genomes/L.mexicana/genome.fa";
	$parseblast=$ARGV[0];
	$MF=$ARGV[1];
	$MFi="$MF\.index";

	$genP="genomemexicana.fa";
	$genN="genomemexicana.fa.revcomp";#genoma revers
	$genPi="$genP\.index";#index del genoma positiu
	$genNi="$genN\.index";#index del genoma negatiu

    	system ("cp $ref_genoma $genP");
    	print "cp\n";
    	system "fastarevcomp $genP > $genN"; #obte el genoma girat
    	print"fastarevcomp\n";
    	system "fastaindex   $genP $genPi"; #fa falta fer un fasta index i un fasta fetch perque les coordenades corresponen als contigs no al genoma. 
    	print"fastaindex1\n";
    	system "fastaindex   $genN $genNi";
    	print"fastaindex2\n";
    	system "fastaindex   $MF $MFi";

# inici de la lectura del parseblast

	open (F1,$parseblast) || die;
	
	my @v;
    	my $longitudvector;
    	my $start;
    	my $end;
    	my $length;
    	my $contig;
    	my $len;
    	my $nomprot;
   	my $line;
    	my ($cont, $g, $gi, $genome_seg, $prot,);
    	my  $pret; 
    	my $p;

	while (<F1>){

    		$line=$_;
    		chomp ($line);

		@v = split(/\s/,$line);

    		$start=$v[2];
    		$end=$v[3];
    		$len=$v[9]; #longitud del genoma L.mexicana
    		$contig =$v[1];
    		$nomprot =$v[0];

 		$nomprot=~/([^#]*)_\#/; #registra el nom nomes fins la # per evitar problemes al posar noms als arxius

		if ($1){
			$nomprot=$1;
			}
		else { 
			 print "ERROR";
			 die;
		     }

  		$prot="$nomprot\.seq"; 				#allarga la sequencia extreta amb el blast 3000 nucleotids per davant i per darrere
    		$start=(($start-3000)<0)?0:($start-3000);	#si el contig s'acaba abans dels 3000 nucleotids el allarga fins el principi del contig
    		$end=(($end+3000)>$len)?$len:($end+3000);	#si el contig s'acaba abans dels 3000 nucleotids el allarga fins el final del contig
   		$length = $end-$start;
    		$genome_seg="$nomprot\_mex\_$start\_$end";
    
		if ($v[4]eq "-"){
        		 $cont="./fastafetchN_$nomprot\_mex\.fa";   #En funcio del strand agafa el genoma invers o el normal  
       			 $g=$genN;
        		 $gi=$genNi;
        			}
   		else{	
       			$cont="./fastafetchP_$nomprot\_mex\.fa";       
        		$g=$genP;
        		$gi=$genPi;
        	     }           
  
     		print "$cont\n";       

#Extraccio del genome
    
		system "fastafetch $g $gi $contig > $cont";
    		print "fastafetch\n";
   		system "fastasubseq $cont $start $length > $genome_seg";#s'ultilitza tambe per la predicció de secis
   		print "fastasubseq\n"; 

#Extraccio de la proteina

    		system "fastafetch $MF $MFi $nomprot > $prot";

#primer genewise
		
		system "genewise -pep -pretty -cdna -gff $prot $genome_seg > $genome_seg\.genewise.fa";
		print "genewise -pep -pretty -cdna -gff $prot $genome_seg\n";
		print "done_genewise1\n";

#canvi de TGA per UCA

		&canviTGA ("$genome_seg");


#segon genewise

		system "genewise -pep -pretty -cdna -gff $prot $genome_seg\.2 > $genome_seg\.genewise.2.fa.";
		print "genewise -pep -pretty -cdna -gff $prot $genome_seg\.2\n";
		print "done_genewise2\n";

#tamany de genewise

		my $t1;
		my $t2;

		$t1= "tamany1"; #primer genewise
		$t2= "tamany2"; #segon genewise

		system "ls -l *genewise.fa* > $t1"; #primer genewise
		system "ls -l *genewise.2.fa* > $t2"; #segon genewise

       		&tamany ("$t1","$t2",$genome_seg);
# Exonerate

		#print "exonerate\n";
	
		#system "exonerate -m p2g --showtargetgff -q $prot -t $genome_seg | egrep -w exon > $genome_seg\.exonerate.gff";
		#print "exonerate -m p2g --showtargetgff -q $prot -t $genome_seg | egrep -w exon > $genome_seg\.exonerate.gff";	
	
		#system "fastaseqfromGFF.pl $genome_seg $genome_seg\.exonerate.gff > $genome_seg\.exonerate.fa";
		#print "fastaseqfromGFF.pl $genome_seg $genome_seg\.exonerate.gff > $genome_seg\.exonerate.fa\n";

# fem arxius per Tcofee
	
		&genewise2pep("$genome_seg\.genewise.2.fa.","$genome_seg"); 
		&genewise2genome("$genome_seg\.genewise.2.fa.","$genome_seg"); 

		system "rm tamany1"; #borra els arxius que no necessitem
		system "rm tamany2";
		system "rm $cont";
		system "rm $genome_seg\.2";	
	}

#canvi de TGA per UCA

	sub canviTGA 
	{
	my $file;
	my $line;
	my $prot;

		$file=$_[0];
		open (F2,$file);
		while (<F2>){
		 	$line=$_;
    			chomp ($line);

			$line=~s/tga/uca/g; #sustitució de tots els tga per uca
			$line=~s/TGA/UCA/g;
			open (OUT, ">>$file\.2");
        		print OUT $line;
			print OUT "\n";
			close (OUT);
         		    } 
	}

#comparem el tamany del genewise

	sub tamany

	{
	my ($file1, $file2);
	my ($line1, $line2);
	my (@v1,@v2);
	my ($tamany1,$tamany2);
	my $nom;

	$file1=$_[0];#primer genewise
		
	open (F3,$file1);
	while (<F3>){
		
		$line1=$_;
    		chomp ($line1);
		@v1 = split(/\s/,$line1);
		
		$tamany1=$v1[4]; #registre el tamany del primer genewise
		   }
	$file2=$_[1];#segon genewise
	
	open (F4,$file2);
	while (<F4>){
	
		$line2=$_;
    		chomp ($line2);
		@v2 = split(/\s/,$line2);

		$tamany2=$v2[4]; #registra el tamany del segon genewise
		   }

	$nom= $_[2];#nom selenoproteina

	if ($tamany1 != $tamany2){

		open (OUT, ">>resultats.selenoproteines.mex"); #arxiu de prediccio de proteines
		print OUT $nom;
		print OUT " probablemente es una selenoproteina";
		print OUT "\n";
		close (OUT);
		}

	else {
		open (OUT, ">>resultats.selenoproteines.mex");
		print OUT $nom;
		print OUT " probablemente NO es una selenoproteina";
		print OUT "\n";
		close (OUT);
	     }
	print"$tamany1 $tamany2\n";
	}

# fem arxiu per tcofee

	sub genewise2pep{
	
	my $file=$_[0];
	my ($l, @pep_list, $p, $pep);
	my $np;
	my $nomprot;
	my $parseblast;
	my @v;
	
	$np=$_[1];
    
	open (F5,$file);
	while (<F5>)
		{	

    		my $line=$_;
    		chomp ($line);
		
                $l.=$line; #posa tot el arxiu en una linea
	    		}
	print "$np\n";

	@pep_list=($l=~/(\>[^\/]+\.pep[^\/]+)\/\//g); #posa tot el arxiu en una linea
	foreach my $p (@pep_list)
		{
    		$p=~s/pep/pep\n/;
	
	open (OUT,">$np\.genewise.pep"); #arxiu amb la proteina obtinguda del genewise
	print OUT $p;
	print OUT "\n";
	close (OUT);         
		}
        		     }	

# arxius per tcofee

 	sub genewise2genome {
	
	my $file=$_[0];
	my ($l, @gem_list, $p, $pep);
	my $np;

	$np=$_[1];   
	open (F6,$file);

	while (<F6>){
			
    		my $line=$_;
    		chomp ($line);
	        $l.=$line; #posa tot el arxiu en una linea
	    	     }

	@gem_list=($l=~/(\>[^\/]+\.sp[^\/]+)\/\//g); #agafa el fragment entre // que conte .sp
	foreach my $g (@gem_list)
		{
	  	$g=~s/sp/sp\n/;
	
	open (OUT,">$np.genwise.gem"); #arxiu amb el cDNA obtingut del genewise
	print OUT $g;
	print OUT "\n";
	close (OUT); 
	        }
        		      }	


