#!/bin/bash

#
# script to automaticaly blast different querys against diferent database
# script adaptet/modified from original found at http://bioinformatica.upf.edu/2011/projectes11/Bx/web/materialsimetodes.html 
#


# The directory "query_protein" must exist and contain the query proteins in fasta-format
# The perl-scripts uni.pl llarg.pl substitucion.pl must be present


mkdir output_blast
mkdir output_blast_detailed
mkdir genomes
mkdir regions
mkdir subseq
mkdir exonerate_output
mkdir cDNA
mkdir translations
mkdir t_coffee_output_exonerate
mkdir t_coffee_output_genewise
mkdir output_genewise
mkdir output_genewise/pretty
mkdir output_genewise/fasta
mkdir output_genewise/fasta/whole
mkdir output_genewise/fasta/longer/
mkdir output_secissearch



#specify genome paths

genome_path="/cursos/BI/genomes/protists/2012"

#read in proteins, cut away pathname

for proteins_fa in ./query_proteins/*; do {
basename_fa=`basename $proteins_fa`
protein=${basename_fa%.fa}
proteins=$proteins" $protein"
} done

#Two loops to combine genomes and proteins, extract query header

rm results_blast.txt 2>> error.txt

       
for protein in $proteins ; do {
	for genomes in A.laibachii_Nc14 A.rara C.fasciculata D.discoideum_AX4 D.fasciculatum F.cylindrus G.niphandrodes I.multifiliis_strain_G5 L.donovani_BPK282A1 L.tarentolae P.capsici P.polycephalum S.arctica T.congolense ; do {
		
		echo $genomes
		pattern=`cat ./query_proteins/$protein.fa | egrep ^">"`
		pattern=${pattern#*>}
		echo $pattern > temppattern
		pattern=`cut -f 1 -d " " temppattern`
		rm temppattern
		echo fasta-header is $pattern
                       
	#blasts standard and -m9 option
    
		blastall -p tblastn -i ./query_proteins/$protein.fa -d $genome_path/$genomes/genome.fa -o ./output_blast/$protein.$genomes.txt -m 9 -F F
		blastall -p tblastn -i ./query_proteins/$protein.fa -d $genome_path/$genomes/genome.fa -o ./output_blast_detailed/$protein.$genomes.txt -F F
			             
	#Directory and hit-file structure, extract target header to 
                       
		echo ${protein}_$genomes: >> results_blast.txt
		egrep "$pattern" output_blast/$protein.$genomes.txt | egrep e- >> results_blast.txt
		egrep "$pattern" output_blast/$protein.$genomes.txt | egrep e- > temp
		x=`egrep "$pattern" output_blast/$protein.$genomes.txt | egrep e-`
		if [ "$x" != "" ]; then
			mkdir hits 2>> error.txt
			mkdir hits/$protein 2>> error.txt
			mkdir hits/$protein/$genomes 2>> error.txt
			egrep -v '^#' temp > temp_stripped
			./uni.pl $protein $genomes > temp_stripped
		fi
		echo ${protein}_$genomes ready
	} done
} done

echo "###BLAST READY" >> error.txt

echo "#################"
echo "###BLAST READY###"
echo "#################"


#automaticaly anotate various blast -m9 hits with exonerate

# Genome indexing with fastaindex
	for genome in /cursos/BI/genomes/protists/2012/* ; do {
		
		stripped=`basename $genome`
		mkdir ./genomes/$stripped 
		echo $genome
		fastaindex $genome_path/$stripped/genome.fa genomes/$stripped/$stripped.index
	} done


# select if exonerate should run in exhaustive mode

if [ "$1" = "-exhaustive" ]; then
exhaustive="--exhaustive yes"
fi

for proteins in hits/* ; do {
	for genome in $proteins/* ; do {
		for hit in $genome/*; do {

		#chop of pathname-part from filenames
		              
			proteins=`basename $proteins`
			genome=`basename $genome`
			hit=`basename $hit`
			hit=${hit%.txt}
			echo ${proteins}_${genome}_${hit} 

		# detection of fasta-header to extract multifasta fragment from -m9 blast

			header_pattern=`cat ./hits/$proteins/$genome/$hit.txt | egrep ^">"`
			header_pattern=${header_pattern#*>}
			echo $header_pattern > tempheader
			header_pattern=`cut -f 1 -d " " tempheader`
			rm tempheader 2>> error.txt       

		# extract header and perform fastafetch
										        
			region=`cat $header_pattern ./hits/$proteins/$genome/$hit.txt | cut -f 2`
			fastafetch $genome_path/$genome/genome.fa ./genomes/$genome/$genome.index $region > ./regions/$region.fa

		# extract beginning and end of hit from -m9 blast

			beg=`cat $header_pattern ./hits/$proteins/$genome/$hit.txt | cut -f 9`
			ending=`cat $header_pattern ./hits/$proteins/$genome/$hit.txt | cut -f 10`

			if [ "$beg" -gt "$ending" ]; then
				beg="$ending"
			fi
												    
			start=$(( $beg - 20000 ))
			length="50000"

		# perform fastasubseq on fetched fragment, various conditions
			
			fastasubseq ./regions/$region.fa $start $length > ./subseq/${proteins}_${genome}_${hit}_genomic.fa 2> error_msg.temp

		# variable scenarios in case sequence is not long enough/too short...
													     
			error=`cat error_msg.temp`
			if [ "$error" != "" ]; then
			error_long=`cat error_msg.temp | egrep "before end"`
			error_short=`cat error_msg.temp | egrep "Unknown flag"`
																         
		# if ending is too near chop of lenght to max lenght
				if [ "$error_long" != "" ]; then
					long=$error_long
					long=${long##*(}
					long=${long%)}
					length=$(( $long - $start))
					fastasubseq ./regions/$region.fa $start $length > ./subseq/${proteins}_${genome}_${hit}_genomic.fa
				fi
																	        
		# If beginning is too near set start to 0
				if [ "$error_short" != "" ]; then
					start="0"
					fastasubseq ./regions/$region.fa $start $length > ./subseq/${proteins}_${genome}_${hit}_genomic.fa 2> error_msg_2.temp
					error_2=`cat error_msg_2.temp`
			
					if [ "$error_2" != "" ]; then
						error_2_long=`cat error_msg_2.temp | egrep "before end"`
						if [ "$error_2_long" != "" ]; then
							long=$error_2_long
							long=${long##*(}
							long=${long%)}
							long=$(( $long - $start))
							fastasubseq ./regions/$region.fa $start $length > ./subseq/${proteins}_${genome}_${hit}_genomic.fa 2> error_msg_3.temp
						fi
					fi 
				fi
			fi

		
		#												   
		# perform genprediciton with exonerate
		#												   
			echo "performing: exonerate"
			
			exonerate -m p2g --showtargetgff -q ./query_proteins/$proteins.fa -t ./subseq/${proteins}_${genome}_${hit}_genomic.fa $exhaustive > ./exonerate_output/exonerate_${proteins}_${genome}_${hit}.gff 2> error_de_substitucio.temp
														         
		# run perl program "substitucion.pl" in case there are selenocysteins marked as U in the sequence

			error_substitucio=`cat error_de_substitucio.temp | egrep "Unknown amino acid"`
			if [ "$error_substitucio" = "** FATAL ERROR **: Unknown amino acid [U]" ]; then
				if [ "$int_UX" = "0" ]; then
					./substitucion.pl < ./query_proteins/$proteins.fa > ${protein}_temporal.txt
					int_UX="1"
				fi
				
				exonerate -m p2g --showtargetgff -q ${proteins}_temporal.txt -t ./subseq/${proteins}_${genome}_${hit}_genomic.fa $exhaustive > ./exonerate_output/exonerate_${proteins}_${genome}_${hit}.gff 2> error_de_substitucio2.temp
			fi
			
			gene=`grep gene ./exonerate_output/exonerate_${proteins}_${genome}_${hit}.gff`
			if [ "$gene" != "" ]; then
				echo ${proteins}_${genome}_${hit} : >> results_exonerate.txt
			fi

			echo "performing: GFF-extraction and translation"

		# extract and translate cDNA from predictet gene-structure with fastaseqfromGFF.pl and fastatranslate in frame 1
																        																	   
			exon=`grep -w exon ./exonerate_output/exonerate_${proteins}_${genome}_${hit}.gff`
			grep -w exon ./exonerate_output/exonerate_${proteins}_${genome}_${hit}.gff > exon.temp
			if [ "$exon" != "" ]; then
				fastaseqfromGFF.pl ./subseq/${proteins}_${genome}_${hit}_genomic.fa exon.temp > ./cDNA/cDNA_${proteins}_${genome}_${hit}.fa
				fastatranslate --frame 1 ./cDNA/cDNA_${proteins}_${genome}_${hit}.fa > ./translations/aa_${proteins}_${genome}_${hit}.fa
			fi
																		    
			echo exonerate ${proteins}_${genome}_${hit} ready
			
		# for each subseq of a hit, realize genewise annotations in "pretty" and fasta mode (check options)
	
		echo "*** performing : genewise "
	
		genewise -pep -both query_proteins/$proteins.fa subseq/${proteins}_${genome}_${hit}_genomic.fa > output_genewise/fasta/whole/genewise_fasta_${proteins}_${genome}_${hit}

		genewise -pep -both -pretty -gff -cdna query_proteins/$proteins.fa subseq/${proteins}_${genome}_${hit}_genomic.fa > output_genewise/pretty/genewise_pretty_${proteins}_${genome}_${hit}
				
																		       
		} done
	} done
} done

rm *.temp
echo "##################################"
echo "###EXONERATE AND GENEWISE READY###"
echo "##################################"

	#select longer genewise annotation

for genewise_out in ./output_genewise/fasta/whole/* ; do {
			
			
			temp_strp=`basename $genewise_out`	
			echo selecting longer hit from $temp_strp			
			./llarg.pl < $genewise_out > ./output_genewise/fasta/longer/$temp_strp.fa
			
} done

	# perform global alignment with t_coffee & secissearch

	#chop of pathname part from filesnames 

for proteins in hits/* ; do {
	for genome in $proteins/* ; do {
		for hit in $genome/*; do {
		
		proteins=`basename $proteins`
        genome=`basename $genome`
        hit=`basename $hit`
		hit=${hit%.txt}
		
		echo "performing: tcoffee"

	#perform global alignment with t_coffee

		if [ -f ./translations/aa_${proteins}_${genome}_${hit}.fa ]; then
			t_coffee query_proteins/$proteins.fa ./translations/aa_${proteins}_${genome}_${hit}.fa > t_coffee_output_exonerate/gloaln_${proteins}_${genome}_${hit}.txt
		fi
		
		if [ -f ./output_genewise/fasta/longer/genewise_fasta_${proteins}_${genome}_${hit}.fa ]; then
			t_coffee query_proteins/$proteins.fa ./output_genewise/fasta/longer/genewise_fasta_${proteins}_${genome}_${hit}.fa > t_coffee_output_genewise/gloaln_${proteins}_${genome}_${hit}.txt
		fi

		echo "*** performing SECISearch of subseq of $hit , from $proteins in $genome"
		cd output_secissearch
		SECISearch.pl -p s -T -H -x ../subseq/${proteins}_${genome}_${hit}_genomic.fa 
		cd ..
		
		} done 
	echo tcoffee of $genome ready
	
	} done
	echo tcoffee of $proteins ready
	
} done

echo "###SECISSEARCH READY###"

echo "##################################"
echo "###T-COFFEE FOR EXONERATE READY###"
echo "##################################"

echo NAME / EXONERATE / GENEWISE >>exonerate_genewise_compare.list

for proteins in hits/* ; do {
	for genome in $proteins/* ; do {
		for hit in $genome/*; do {
			
		proteins=`basename $proteins`
        genome=`basename $genome`
        hit=`basename $hit`
               
        if [[ "$hit" = "*" ]] ; then
        	echo shit ; else
        
        hit=${hit%.txt}
		
		if [ -f ./t_coffee_output_exonerate/gloaln_${proteins}_${genome}_${hit}.txt ]; then
			ali_exo=`head -n 1 ./t_coffee_output_exonerate/gloaln_${proteins}_${genome}_${hit}.txt | cut -d ',' -f 3,5` ; else
			ali_exo="NONE"
		fi
		
		if [ -f ./t_coffee_output_genewise/gloaln_${proteins}_${genome}_${hit}.txt ]; then
			ali_gen=`head -n 1 ./t_coffee_output_genewise/gloaln_${proteins}_${genome}_${hit}.txt | cut -d ',' -f 3,5` ; else
			ali_gen="NONE"
		fi
		
			
				
		echo  $proteins $genome $hit / $ali_exo / $ali_gen >> exonerate_genewise_compare.list
		fi
		
		} done
	} done
} done

echo "lei lei!"

#!/usr/bin/perl
use warnings;
use strict;
my $contador=0;

while ()
{
    open (OUT,">hits/$ARGV[0]/$ARGV[1]/hit$contador.txt");
	print OUT $_;
    close OUT;
    $contador++;
}

#!/usr/bin/perl

use strict;

my @v;
my $string;
my $counter = 0;

while ()
{
    if ($counter == 0)
    {
	print $_;
    }
    if ($counter > 0)
    {

	@v=split(//,$_);
	
	while (<@v>)
	{
	    if ($_ eq "U")
	    {
		$string.= "X";
	    }
	    else
	    {
		$string.=$_;
	    }
	}
    }
    $counter++;
    $string.="\n";
}

print $string;

#!/usr/bin/perl

#CAL QUE L'INPUT SIGUI EN FORMAT GENEWISE

use strict;

my @v;
my $inici=0;
my $longv=0;
my $c1=0;
my $par1=0;
my $par2=0;
my $temp1="";
my $temp2="";


while ()
{
    @v=split(//,$_);
    $longv=scalar(@v);
    $c1=0;
   
	while ($c1 < $longv){
		if ($v[$c1] eq ">"){
	    	$inici++;
		}
	
		if ($inici==2){
			$temp2.=$v[$c1];
			$par2++;
		}
		else{
	    $temp1.=$v[$c1];
	    $par1++;
		}
		
	$c1++;
	
	}

}

if ($par1 > $par2){
	print $temp1;
}
else{
	print $temp2;
}