#!/bin/bash




#aquest programa es va fer per analitzats un a un aquells scaffolds d'interès, quan els resultats obtinguts amb el hit de milor e-value no eren bons. 

### EXPORTS ###
export	PATH=/mnt/NFS_UPF/bioinfo/BI/bin:$PATH     #per 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 al genoma de Calorhinus 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="./proteinesespecials"
path_scaffolds="./scaffolds"
path_subseq="./subseqs"
path_exonerate="./exonerate"
path_GFF="./GFF"
path_secis="./secis"
path_tcoffee="./tcoffee"



#creem les carpetes que necessitarem (-p si ja està creada no es 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
mkdir -p $path_out

### FASTAFETCH ###



	echo "introduim la proteina que volem estudiar"
	read proteina
	echo "$proteina"

	echo "introdueix un scaffold"

	read scaffold

	echo "has introduit scaffold: $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 50000 a la posició final i restem 50000 a la posició d'inici
	echo "introdueix inici"

	read inici

	echo "introdueix final"

	read final

	echo "has introduit inici: $inici i final: $final"

	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 -n 1 -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"    #obtenim només les linies que contenten la paraula ego

    
### 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"
	echo "canviant els * per X..."
    	sed 's/*/X/g' $path_GFF/$proteina.pred.aa.fa > "$path_GFF/$proteina.pred.aa.be.fa"    #canviem els * per X, ja que sino t-coffee dona error
    
   
### T-COFFEE ###
    echo "s'està fent el t-coffee de $proteina..."
    t_coffee $path_GFF/$proteina.pred.aa.be.fa $path_proteins/$proteina > "$path_tcoffee/$proteina.tcoffee"