#!/bin/bash

#Carreguem els mòduls i les paths
module load modulepath/goolf-1.7.20
module load BLAST+/2.2.30-goolf-1.7.20
module load Exonerate/2.2.0-goolf-1.7.20
module load T-Coffee/11.00.8cbe486-goolf-1.7.20
export PATH=/cursos/20428/BI/bin:$PATH
export PATH=/cursos/20428/BI/soft/genewise/x86_64/bin:$PATH
export WISECONFIGDIR=/cursos/20428/BI/soft/genewise/x86_64/wise2.2.0/wisecfg/

#Definim les nostres paths
PATH_output=/home/u114888/BI/output #Directori que contindrà tots els directoris i fitxers de soritda.  
PATH_genomes=/cursos/20428/BI/genomes/2017/Anas_zonorhyncha #On tenim el genome.fa i el genome.index de l'ànec.  
PATH_protselenoDB=/home/u114888/BI/protselenoDB #Selenoproteïnes de Chicken del SelenoDB

#Canviem el directori i entrem al path de les proteïnes del Chicken del selenoDB que acabem de crear i definim la variable de les proteïnes.

cd $PATH_protselenoDB

output="dio1 dio2 dio3 eefsec gpx3 gpx7 gpx8 msra msrb1 msrb3 pstk sbp2_1 sbp2_2 sbp2_3 secp43 secs_1 secs_2 sel15 selenoh selenoi selenok_1 selenok_2 selenok_3 selenon selenoo_1 selenoo_2 selenoo_3 selenoo_4 selenop selenos selenot selenou1 sephs txnrd2 txnrd3 txnrd"


#COMENÇA EL TBLASTN

#Per entrar en un loop per cada selenoproteïna:

for element in $output

#Fem anar el software blastall per cada selenoproteïna del Chicken amb la base de dades del genoma problema:

do
    tblastn -query $PATH_protselenoDB/$element.fa -db $PATH_genomes/genome.fa -out $PATH_output/$element/tblastn_$element.fa
    tblastn -query $PATH_protselenoDB/$element.fa -db $PATH_genomes/genome.fa -out $PATH_output/$element/tblastn2_$element.fa -outfmt 7
 
done

#ACABA EL TBLASTN: l'output és tblastn_proteïna.fa i tblastn2_proteïna.fa


#COMENÇA EL PREFASTAFETCH

for element in $output
do
    echo 'Processant...'

    cd $PATH_output/$element

    grep -v '^#' $PATH_output/$element/tblastn2_$element.fa > $PATH_output/$element/temporary_$element.fa

    awk '$3>75 && $4>25 && $11<0.000001 {print $2"\t"$3"\t"$4"\t"$9"\t"$10"\t"$11"\t"$12}' $PATH_output/$element/temporary_$element.fa > $PATH_output/$element/hits_$element.fa
    
    rm $PATH_output/$element/temporary_$element.fa

    echo  'Finalitzat' $element

done

#ACABA EL PREFASTAFETCH: l'output és hits_proteïna.fa


#COMENÇA EL FASTAFETCH

for element in $output
do
    echo 'Processant fastafetch...' $element
    cd $PATH_output/$element
   
    awk '$1 {print $1}' $PATH_output/$element/hits_$element.fa > $PATH_output/$element/temp_contigs_$element.txt

    sort -b -d -u $PATH_output/$element/temp_contigs_$element.txt -o $PATH_output/$element/temp_contigs_$element.txt

    contigs=($(cat $PATH_output/$element/temp_contigs_$element.txt))
        
    for contig in "${contigs[@]}"
    do
fastafetch $PATH_genomes/genome.fa $PATH_genomes/genome.index "${contig}" > $PATH_output/$element/"${contig}".fa
    done
    
    rm $PATH_output/$element/temp_contigs_$element.txt

    echo  'Fastafetch finalitzat' $element
done

#ACABA EL FASTAFETCH: l'output és contig.fa


#COMENÇA EL POSTFASTAFETCH
for element in $output
do
    cd $PATH_output/$element/
    
    rm $PATH_output/$element/hits2_$element.fa
    touch $PATH_output/$element/hits2_$element.fa

    awk '{print $1}' $PATH_output/$element/hits_$element.fa > $PATH_output/$element/tempcontig.fa

    sort -u -b -d $PATH_output/$element/tempcontig.fa -o $PATH_output/$element/tempcontig.fa #S'HA DE BORRAR
    
    contig=($(cat $PATH_output/$element/tempcontig.fa))
        
    for name in "${contig[@]}"
    do
	echo $name
	
	grep -v "^>" $PATH_output/$element/"${name}".fa > $PATH_output/$element/postgrep_$element.fa #S'HA DE BORRAR

	length=$(tr -d '[:space:]' < $PATH_output/$element/postgrep_$element.fa | wc -c | awk '{print $1}')
	
	while read line
	do
	    awk=$(echo $line | awk '{print $1}')     
	        
	        if [ "${awk}"  == "${name}" ]
		        then
		    
		    echo $line | awk -v VAR="${length}" '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"VAR}' >> $PATH_output/$element/hits2_$element.fa

    fi
    
done < $PATH_output/$element/hits_$element.fa

    done
    
    rm $PATH_output/$element/tempcontig.fa
    rm $PATH_output/$element/postgrep_$element.fa

done

#ACABA EL POSTFASTAFETCH: l'output és hits2_proteïna.fa


#COMENÇA EL FASTASUBSEQ

for element in $output
do

    cd $PATH_output/$element/
    ~/BI/programes/numeros.pl $PATH_output/$element/hits2_$element.fa > ./numeros_$element.txt
   
    awk -v var=$PATH_output/$element '{print "fastasubseq "var"/"$2".fa "$4" "$5" > "var"/fastasubseq_"$2"_"$1".fa"}' numeros_$element.txt > fastasubseq_commands.sh
    sh fastasubseq_commands.sh
    
    rm $PATH_output/$element/fastasubseq_commands.sh
done

#ACABA EL FASTASUBSEQ: l'output és fastasubseq_contig_NumeroDeHit.fa


#COMENÇA EL EXONERATE

for element in $output
do
 
   sed 's/U/*/g' $PATH_protselenoDB/$element.fa > $PATH_output/$element/$element.fa
   
   cd $PATH_output/$element
   
   subseq=($(ls fastasubseq*.fa))
    
    for seq in "${subseq[@]}"
    do

    exonerate -m p2g --showtargetgff -q $PATH_output/$element/$element.fa -t $PATH_output/$element/$seq | egrep -w exon > exon_$seq.gff

    done

done

#ACABA EXONERATE: l'output és exon_contig_NumeroDeHit.gff


#COMENÇA EL FASTASEQFROMGFF

for element in $output
do

    cd $PATH_output/$element
    subseq=($(ls fastasubseq*.fa))
    gff=($(ls *.gff))
    
    for seq in "${subseq[@]}"
    do
    
	fastaseqfromGFF.pl $PATH_output/$element/$seq  $PATH_output/$element/exon_$seq.gff > preprot_$seq.pl
	
    done

done

#ACABA EL FASTASEQFROMGFF: l'output és preprot_contig_NumeroDeContig.pl


#COMENÇA EL FASTATRANSLATE

for element in $output
do
    cd $PATH_output/$element
    
    preprot=($(ls fastasubseq*.fa))
    
    for prep in "${preprot[@]}"
    do
	
	fastatranslate -f $PATH_output/$element/preprot_$prep.pl -F 1 > $PATH_output/$element/prot_$prep
    
    done
done

#ACABA EL FASTATRANSLATE: l'output és prot_fastasubseq_contig_NumeroDeContig.fa


#COMENÇA EL TCOFFEE

for element in $output
do

    cd $PATH_output/$element

      preprot=($(ls fastasubseq*.fa))
    
      for prep in "${preprot[@]}"
      do
  sed 's/*/U/g' $PATH_output/$element/prot_$prep > prot2_$prep
  t_coffee $PATH_protselenoDB/$element.fa $PATH_output/$element/prot2_$prep > tcoffee_$prep.txt  
  
      done
    

done

#ACABA EL TCOFFEE: l'output és tcoffee_fastasubseq_contig_NumeroDeContig.fa