#!/bin/bash


### First we have already created a folder named proteins, which contains the file with the 
queries (humanselenoprot_prot.fa) and the program extractu.pl.

### We also have a folder named fastafetch, with the programs forfastafetch.pl and forfastasubseq.pl.

### BLAST. We start in our "home" folder. We create a folder named blast. 
We copy the file containing the queries (humanselenoprot_prot.fa) and run tblastn in it. We set 
the-m9 option and filter out the hits with an e-value smaller than 0.001.

mkdir blast

cd blast

cp ../proteins/humanselenoprot_prot.fa .

export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH #In order to run the different programs of the script, we need to export their paths.
blastall -p tblastn -i humanselenoprot_prot.fa -d /cursos/BI/genomes/project_2013/Saimiri_boliviensis/genome.fa -m9 -e 0.001 -o blast

### We make a file without the lines starting with "#".

egrep -v "#" blast > ../inputfromblast.txt

### Next, Us in protein sequences are replaced by Xs.  For this, the programs extractu.pl
is required in the proteins folder.Then, individual protein sequences are extracted into indindividual files.
We do this with the programs fastaindex and fastafetch, as we will do for the genomic sequences.

cd ../proteins
./extractu.pl < humanselenoprot_prot.fa > protnou.fa

export PATH=/cursos/BI/soft/exonerate/i386/bin:$PATH

fastaindex protnou.fa protnou.index

for protein in `cut -f 1 ../inputfromblast.txt | uniq` ; do {
fastafetch protnou.fa protnou.index $protein > "${protein}.fa"
}
done



### In order to run fastafetch, first we use the program forfastafetch.pl to get a file with the queries
and the names of the scaffolds we want to retrieve.

cd ../fastafetch
./forfastafetch.pl < ../inputfromblast.txt > protein_genomic.txt

cut -f 2 protein_genomic.txt > inputforfastafetch.txt #This makes a file with the names of the 
scaffolds.

### Then we run FASTAFETCH (first fastaindex).

export PATH=/cursos/BI/soft/exonerate/i386/bin:$PATH
fastaindex /cursos/BI/genomes/project_2013/Saimiri_boliviensis/genome.fa genome.index

for name in `cat inputforfastafetch.txt`; do {

fastafetch /cursos/BI/genomes/project_2013/Saimiri_boliviensis/genome.fa genome.index $name > "${name}.fa"

} done  

### We change the names of the scaffolds to make the handling of the files easier. Each scaffold
is named with the name of the protein and the hit. When two or more proteins have hits on the 
same scaffold, a copy is generated for each protein.

nom=a
while read protein genome
do
  
if [ "$nom" != "$protein" ]; then
		num=1
		nom=$protein
		cp "${genome}.fa" "${protein}_${num}_genomic.fa"
      else
		num=$[$num+1]
		cp "${genome}.fa" "${protein}_${num}_genomic.fa"
fi
done < protein_genomic.txt

### Then we generate a file with the information we need for running fastasubseq with the program
forfastasubseq.pl. It is an extension of the program forfastafetch.pl and in this case the name 
of the scaffold is the new name.

./forfastasubseq.pl < ../inputfromblast.txt > inputforfastasubseq.txt

### Then we run fastasubseq. First we make a folder named subsequences. 

mkdir ../subsequences

while read query genomic start length orientation
do

fastasubseq $genomic $start $length > ../subsequences/"${genomic}_exonerate.fa" 2> error.temp
error=`cat error.temp` #this is to ammend those cases in which the length goes beyond the end of the scaffold, which gives an error.
if [ "$error" != "" ]; then
   final=`cat error.temp | egrep "before end" | cut -d '(' -f 2 | sed 's/)//'`
   length=$(( $final - $start))
   fastasubseq $genomic $start $length > ../subsequences/"${genomic}_exonerate.fa" 
fi
done < inputforfastasubseq.txt


### Then we run exonerate (from home).

cd ..
mkdir outputexoneratealign
mkdir outputexonerate
mkdir multiplepredictions

export PATH=/cursos/BI/soft/exonerate/i386/bin:$PATH

export PATH=/cursos/BI/bin:$PATH

while read query genomic start length orientation
do

genomicexonerate="${genomic}_exonerate.fa"

exonerate -m p2g --showtargetgff -q ./proteins/$query -t ./subsequences/$genomicexonerate > outputexoneratealign/"${genomicexonerate}.gff"

done  < ./fastafetch/inputforfastasubseq.txt

###Then we transform the exonerate output into fasta and translate the sequence. For that, 
we need to take into account that exonerate might give more than one prediction. In this
case, we only consider those predictions in the same orientation as the blast hit.

ls outputexoneratealign > listexonerate.txt

for name in `cat listexonerate.txt`; do {
subsequence=`echo $name | sed 's/\.gff//'` 
   
numprediction=`egrep "C4 Alignment" ./outputexoneratealign/$name | uniq -c | cut -d "C" -f 1`

if [ $numprediction == 1 ]; then 
	cat ./outputexoneratealign/$name | egrep -w exon > ./outputexonerate/"${name}.gff"

	fastaseqfromGFF.pl ./subsequences/$subsequence ./outputexonerate/"${name}.gff" > ./outputexonerate/"${name}fromgff.fa"  

else
	i=0
	region1=no
	region2=no
	while read file #read the file line by line
	do
  	if [ $region1 = "no" ]; then
      		field=`echo $file | cut -d ":" -f 1` 

      		if [ "$field" = "Target range" ]; then #here is where exonerate gives the start and ending position of the prediction
		    start="`echo $file | cut -d ":" -f 2 | cut -d "-" -f 1 | cut -d " " -f 2`" 			
		    end="`echo $file | cut -d ">" -f 2 | cut -d " " -f 2`" 
	    
	  		if [ "$start" -lt "$end" ]; then
	  		    orientation="+"

	  		else
	  		    orientation="-"

	  		fi

			hit=`basename $name _exonerate.fa.gff`; 
			orientationblast=`cat ./fastafetch/inputforfastasubseq.txt | egrep "$hit" | cut -f 5`
			echo orientationblast $orientationblast
	  		if [ "$orientationblast" = "$orientation" ]; then
             		     #then do the other steps
	     			region1=yes

	  		fi
		fi
 	fi
	if [ "$region1" = "yes" ]; then
		    
	     	if [ "$region2" = "no" ]; then
			 compare=`echo $file | cut -c 7-11`
			 
			if [ "$compare" = "START" ]; then
		     		i=$[$i + 1]
		     		region2=yes
		 	fi
	      
	     	fi

	     	if [ "$region2" = "yes" ]; then

			 echo "$file" >> ./multiplepredictions/"${name}${i}"
 
			 compare=`echo $file | cut -c 7-9`
     
	  		 if [ "$compare" = "END" ]; then
		    		 region2=no
		     		region1=no
   
			 fi  

	     	 fi 
     
	fi       
	done < ./outputexoneratealign/$name

fi
    
  
}

done

ls multiplepredictions > listgff.txt 

for file in `cat listgff.txt` ; do {

subsequence=`echo $file | sed 's/\.gff[0-9]//'` 


cat ./multiplepredictions/$file | egrep -w exon > ./outputexonerate/"${file}.gff"

fastaseqfromGFF.pl ./subsequences/$subsequence ./outputexonerate/"${file}.gff" > ./outputexonerate/"${file}fromgff.fa"  

} done  

# We make a file with the names of the fasta files generated in the previous step. These files 
contain the exonerate predictions that will be translated next.

ls outputexonerate | egrep "fromgff" > ./outputexonerate/totranslate.txt

mkdir translated

cd outputexonerate

export PATH=/cursos/BI/soft/exonerate/i386/bin:$PATH

for sequence in `cat totranslate.txt` ; do {
i=1
while [ $i -lt 4 ] #although we expect the good translation will be for the frame 1, we also consider the frames 2 and 3
do
fastatranslate  -F $i $sequence > ../translated/"${sequence}_${i}_translated.fa"
i=$[i+1]
done
}
done 


### T-COFFEE. Exonerate predictions and queries are aligned using T-COFFEE. 


cd ..
mkdir tcoffeeexonerate
cd tcoffeeexonerate
ls ../translated > namestcoffeeexonerate.txt

export PATH=/cursos/BI/bin:$PATH 

for sequence in `cat namestcoffeeexonerate.txt` ; do {

protein=`echo $sequence | sed 's/_[0-9][0-9]*_genomic\.fa_exonerate\.fa\.gff[0-9]*fromgff\.fa_[0-9]_translated\.fa//'
`

t_coffee ../translated/$sequence ../proteins/"${protein}.fa"
} done

cd ..

### SECISEARCH

mkdir secis
cd secis
ls ../subsequences > subsequences.txt

export PATH=/cursos/BI/bin:$PATH
while read subseq
do
SECISearch.pl ../subsequences/$subseq
done < subsequences.txt

cd ..

### GENEWISE preparation.

ls subsequences > subsequencesgenewise.txt

mkdir genewise

cd genewise

mkdir outputgenewise
mkdir outputgenewiserev

###First we run genewise for the forward strand.

for subsequence in `cat ../subsequencesgenewise.txt`; do {
query=`echo $subsequence | sed 's/_[0-9][0-9]*_genomic\.fa_exonerate\.fa//'`
genewise -pep -pretty -cdna -gff ../proteins/${query}.fa ../subsequences/$subsequence > ./outputgenewise/"${subsequence}_genewise.gff"
} done

###Then we run genewise for the reverse strand.

for subsequence in `cat ../subsequencesgenewise.txt`; do {
query=`echo $subsequence | sed 's/_[0-9][0-9]*_genomic\.fa_exonerate\.fa//'`
genewise -pep -pretty -cdna -trev -gff ../proteins/${query}.fa ../subsequences/$subsequence > ./outputgenewiserev/"${subsequence}_genewiserev.gff"
} done


###Then we need to get the protein prediction from the genewise output in order to align it with
the query using T-Coffee.

ls outputgenewise | egrep "genewise.gff" > namesgenewise.txt
ls outputgenewiserev | egrep "genewiserev.gff" > namesgenewiserev.txt

cd outputgenewise
cp ../namesgenewise.txt .

for name in `cat namesgenewise.txt` ; do {
fet=false
region=no
while read file 
do
if [ "$fet" = "false" ]; then 

    if [ $region = "no" ]; then
	compare="`echo $file | cut -c 1-3`"
 
	if [ "$compare" = ">gi" ]; then
	    region=yes
 	    i=0
	fi

    else

	compare="`echo $file | cut -c 1-2`"
	if [ $compare != "//" ]; then
	    i=$[i + 1]
	else
 
	  fet=true
	fi

    fi  
fi
done < $name

cat $name | egrep "pep" -A $i > ${name}_fortcoffee.fa

} done

cd ../outputgenewiserev
cp ../namesgenewiserev.txt .

for name in `cat namesgenewiserev.txt` ; do {
fet=false
region=no
while read file
do
if [ "$fet" = "false" ]; then

    if [ $region = "no" ]; then
        compare="`echo $file | cut -c 1-3`"

        if [ "$compare" = ">gi" ]; then
            region=yes
            i=0
        fi

    else

        comparar="`echo $file | cut -c 1-2`"
        if [ $comparar != "//" ]; then
            i=$[i + 1]
        else

          fet=true
        fi

    fi
fi
done < $name

cat $name | egrep "pep" -A $i > ${name}_fortcoffee.fa

} done

cd ..


mkdir tcoffeegenewise
mkdir tcoffeegenewiserev

ls outputgenewise | egrep "fortcoffee.fa" > ./tcoffeegenewise/namestcoffee.txt
ls outputgenewiserev | egrep "fortcoffee.fa" > ./tcoffeegenewiserev/namestcoffeerev.txt

cd tcoffeegenewise

export PATH=/cursos/BI/bin:$PATH 

for sequence in `cat namestcoffee.txt` ; do {

protein=`echo $sequence | sed 's/_[0-9][0-9]*_genomic\.fa_exonerate\.fa_genewise\.gff_fortcoffee\.fa//'`

t_coffee ../outputgenewise/$sequence ../../proteins/"${protein}.fa"
} done

cd ../tcoffeegenewiserev

for sequence in `cat namestcoffeerev.txt` ; do {

protein=`echo $sequence | sed 's/_[0-9][0-9]*_genomic\.fa_exonerate\.fa_genewiserev\.gff_fortcoffee\.fa//'`

t_coffee ../outputgenewiserev/$sequence ../../proteins/"${protein}.fa"
} done

#!/usr/bin/perl -w

use strict;

#This program is used to substitute the seleneocysteines (U) in the protein queries by X.

my @line;
my $i;
my $newline;


while (<STDIN>){
chomp($_);
$i = 0;
	@line = split ("", $_);
		if ($line[$i] eq ">"){
		$newline = join ("", @line);
		}
		else {
		 while ($i < scalar(@line)){
		 if ($line[$i] eq "U"){

			$line[$i] = "X";
		 }
		 $i = $i + 1;
		 }
		$newline = join ("", @line);	
}
print "$newline\n";
}

#!/usr/bin/perl -w

use strict;

#this program generates a file with the names of the queries and the scaffolds for which we got 
blast hits. 


     
my @linia;
my $a;
my $human;#protein name
my $nom;#scaffold name
my @output;
my $humancomparacio;

$human = "a";
$nom = "a";
$inici = 1;
$fi = 1;
$humancomparacio = "a";

while (<STDIN>){ 
chomp($_);
$a = $_;
@linia = split ("\t", $a);
if ($linia[0] ne $humancomparacio) {#first line or new protein
					if ($humancomparacio ne "a") {#if not first line
					    
					    @output = join("\t",$human, $nom);
					    print "@output\n";
					}
		      $humancomparacio=$linia[0];
		      $human = $linia[0];
		      $nom = $linia[1];
		
		
		      
				    }
else{#it's still the same protein
	 if ($linia[1] ne $nom){#new scaffold
				    @output = join("\t", $human, $nom);
				    print "@output\n";

				    $human = $linia[0];
				    $nom = $linia [1];

				}
	 

	}
}

@output = join("\t", $human, $nom);
print "@output\n";

#!/usr/bin/perl -w

use strict;

#This program uses the blast tabular output (without the lines starting by #) and produces a file
with five columns:query name, scaffold name, start for fastasubseq, length for fastasubseq,
gene orientation. 


my @linia;
my $a;
my $human;#protein name
my $nom;#scaffold name
my $inici;#start
my $fi;#end
my @output;
my $length;
my $numero;#scaffold counter
my $humancomparacio;#variable with the name of the protein query. It is used to know whether 
every line refers to the same protein as the previous line or a new one.
my $orientacio;#orientation
$human = "a";
$nom = "a";
$inici = 1;
$fi = 1;
$humancomparacio = "a";
while (<STDIN>){ 
chomp($_);
$a = $_;
@linia = split ("\t", $a);
if ($linia[0] ne $humancomparacio) {#the first file line or new protein
					if ($humancomparacio ne "a"){ #unless it is the first line
					    $length = ($fi - $inici + 300000);
					    $inici = ($inici - 150000);
					    if ($inici < 0) { $inici = 1;}
					    $nom = $human."_".$numero."_"."genomic.fa";
					    $human = $human.".fa";
					    
					    @output = join("\t", $human, $nom, $inici, $length, $orientacio);
					    print "@output\n";
					}
		      $numero = 1;
		      $human = $linia[0];
		      $humancomparacio = $linia[0];
		      $nom = $linia[1];
		      if ($linia[8] < $linia[9]){
			  $inici = $linia[8];
			  $fi = $linia[9];
			  $orientacio = "+";
		      }
		      else {
			  $inici = $linia[9];
			  $fi = $linia[8];
			  $orientacio = "-";
		      }
		      
		  }
else{#the protein is not new

	 if ($linia[1] ne $nom) {#new scaffold of the same protein
					#then we print what we have from the previous scaffolds
					$length = ($fi - $inici + 300000);
					$inici = ($inici - 150000);
					if ($inici < 0) { $inici = 1;}
					$nom = $human."_".$numero."_"."genomic.fa";
					$human = $human.".fa";

					@output = join("\t", $human, $nom, $inici, $length, $orientacio);
					print "@output\n";
					$numero = $numero + 1;#we need to add 1 because it is a new scaffold
					$human = $linia[0];
					$nom = $linia[1];
					if ($linia[8] < $linia[9]){
					    $inici = $linia[8];
					    $fi = $linia[9];
					    $orientacio = "+";
					}
					else {
					    $inici = $linia[9];
					    $fi = $linia[8];
					    $orientacio = "-";
					}
	}
	 if ($linia[1] eq $nom){ #the protein and the scaffold are the same, there is no need to add anything else
    if ($linia[8] < $linia[9]) {#forward
	if ($linia[8] < $inici) {#substitute the start value if it is smaller
	    $inici = $linia[8];
	}
	if ($linia[9] > $fi){#substitute the end value if it is greater
	    $fi = $linia[9]; 
	}
    }
    else #reverse
    {
	if ($linia[9] < $inici){
	    $inici = $linia[9];
	}
	if ($linia[8] > $fi){
	    $fi = $linia[8];

	}
    }
 }
   
}
}
$length = ($fi - $inici + 300000);
$inici = ($inici - 150000); if ($inici < 0) {$inici = 1;}
$nom = $human."_".$numero."_"."genomic.fa";
$human = $human.".fa";

@output = join("\t", $human, $nom, $inici, $length, $orientacio);
print "@output\n";

#print "$human","$nom","$inici","$fi","\n";