#!/usr/bin/bash

### EXPORTS ###
export	PATH=/mnt/NFS_UPF/bioinfo/BI/bin:$PATH     #para el fastasubseqformGFF.pl

### PATHS ### 

path_files="/mnt/NFS_UPF/bioinfo/BI/genomes/2018/Callorhinus_ursinus"      #per genome.index I genome.length
path_db="$path_files/genome.fa"    #per el genoma de Callorhinus ursinus

#treballem des de la carpeta callorhinus, l'hem creat previament, en aquesta carpeta trobarem la carpeta proteinashumano on es troben totes les selenoproteines obtingudes a partir del selenoDB i el programa
path_out="./blast"
path_proteins="./proteinashumano"
path_scaffolds="./scaffolds"
path_subseq="./subseqs"
path_exonerate="./exonerate"
path_GFF="./GFF"
path_secis="./secis"
path_tcoffee="./tcoffee"


#creem les carpetes que utilitzarem (-p farà que si la carpeta ja està creada no es torni a crear)
mkdir -p $path_scaffolds
mkdir -p $path_subseq
mkdir -p $path_exonerate
mkdir -p $path_GFF
mkdir -p $path_secis
mkdir -p $path_tcoffee

#### TBLASTN ####

mkdir -p $path_out

#Entrem en el loop cada proteina

for prot in $(ls $path_proteins) 
do
	echo "Analysing the protein: $prot"
	if [ -e $path_out/$prot.blast ]     #document $prot.blast a la carpeta blast
		then
			echo "$prot.blast ja s'ha creat"
		else
			echo "Fent el blast per $prot..."
			tblastn -query $path_proteins/$prot -db $path_db -evalue 0.001 -outfmt 6 -out $path_out/$prot.blast 
	fi
done

### FASTAFETCH ###

for proteina in $(ls $path_proteins) 
do

    #ens basem en el e-value per escollir el millor scaffold. Del blast ja obtenim els resultats de millor a pitjor e-value
	scaffold=$(head -n 1 $path_out/$proteina.blast | cut -f2)
	inici=$(head -n 1 $path_out/$proteina.blast | cut -f9)
	final=$(head -n 1 $path_out/$proteina.blast | cut -f10)
	lengthseq=$(cat $path_files/genome.lengths | egrep $scaffold | cut -d ' ' -f1)

  	echo "$scaffold $inici $final $lengthseq" > $path_scaffolds/$proteina.scaffold

	cut -f2 $path_out/$proteina.blast | sed 's/\.1//g' | sort | uniq -c > $path_scaffolds/$proteina.scaffolds_list      #sort ordena els scaffolds per ordre alfabetic o numèric, uniq elimina aquells 	scaffolds que estan repetits. Posem la columna amb els scaffolds a un altre document guardat dins de la proteina que esta dins de la carpeta scaffold.


	echo "Fent el fastafetch per el scaffold: $scaffold"
	fastafetch $path_db $path_files/genome.index $scaffold > "$path_scaffolds/$proteina.fastafetch"   #obtenim la seqüència

	
### FASTASUBSEQ ###

	#volem una seqüència que cobreixi tota la proteïna (sumem 500000 a la posició final i li restem 500000 a la posición d'inici)


	if (($inici>$final)); then
	    inici1=$final
	    final1=$inici
	    echo "la seqüència és reversa"
	else
	    inici1=$inici
	    final1=$final
	    echo "la seqüència és forward"
	fi
	
	inici2=$(($inici1-50000))
	final2=$(($final1+50000))
	length=$(($final2-$inici2))

	if (($inici2>0)); then
	    start=$inici2
	else
	    start=0
	fi

	if (($final2>$lengthseq)); then
	    lengthfinal=$(($lengthseq-$start))
	else
	    lengthfinal=$length
	fi

	echo "S'està fent el fastasubseq de $scaffold de la $proteina..."

	fastasubseq $path_scaffolds/$proteina.fastafetch $start $lengthfinal > "$path_subseq/$proteina.subseq"

### EXONERATE ###
	echo "s'està fent el exonerate de $scaffold de la $proteina..."
  	
	exonerate -m p2g --showtargetgff -q $path_proteins/$proteina -t $path_subseq/$proteina.subseq > "$path_exonerate/$proteina.exonerate.gff"
	
	egrep -w exon $path_exonerate/$proteina.exonerate.gff > "$path_exonerate/$proteina.exon.gff"   #agafem només aquelles files que tenen la paraula exon


### FASTASEQFROMGFF###
    
 	echo "s'està fent el fastaseqfromGFF de la $proteina..."

	fastaseqfromGFF.pl $path_subseq/$proteina.subseq $path_exonerate/$proteina.exon.gff  > "$path_GFF/$proteina.pred.nuc.fa


### FASTATRANSLATE ###
	echo "s'està fent el fastatranslate de $proteina..."
 	fastatranslate -F 1 $path_GFF/$proteina.pred.nuc.fa > "$path_GFF/$proteina.pred.aa.fa"
	sed 's/*/X/g' $path_GFF/$proteina.pred.aa.fa > $path_GFF/$proteina.pred.aa.be.fa.    #canviem els * per X ja que sino el t-coffee (pas segment) dona error

### T-COFFEE ###

    echo "s'està fent el t-coffee de $proteina..."
    t_coffee $path_GFF/$proteina.pred.aa.fa $path_proteins/$proteina > "$path_tcoffee/$proteina.tcoffee"
done