#!/bin/bash

##############################################AUTOMATITZACIÓ DEL BLAST###################################################

export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH;
cp /cursos/BI/bin/ncbiblast/.ncbirc ~/;

echo Quines selenoproteïnes vols?;
read selenoprot;
#selenoprot="Sel15_1 Sel15_2 Sel15_3 EEFSec_1 EEFSec_2 GPx1_1 GPx1_2 GPx2_1 GPx3_1 GPx3_2 GPx3_3 GPx3_4 GPx4_1 GPx4_2 GPx5_1 GPx5_2 GPx6_1 GPx6_2 GPx7_1 GPx8_1 GPx8_2 GPx8_3 DI1_1 DI1_2 DI1_3 DI1_4 DI1_5 DI1_6 DI1_7 DI2_1 DI2_2 DI2_3 DI2_4 DI2_5 DI3_1 MsrA_1 MsrA_2 MsrA_3 MsrA_4 PSTK_1 SBP2_1 SBP2_2 SecS_1 SecS_2 SPS1_1 SPS1_2 SPS1_3 SPS1_4 SPS2_1 SelH_1 SelH_2 SelI SelK SelM_1 SelM_2 SelN_1 SelN_2 SelO SelP_1 SelP_2 SelP_3 SelP_4 SelP_5 SelR1_1 SelR1_2 SelR1_3 SelR1_4 SelR2 SelR3_1 SelR3_2 SelR3_3 SelR3_4 SelS_1 SelS_2 SelS_3 SelS_4 SelS_5 SelT_1 SelT_2 SelT_3 SelT_4 SelU1_1 SelU1_2 SelU1_3 SelU1_4 SelU2_1 SelU2_2 SelU3_1 SelU3_2 SelV SelW1_1 SelW1_2 SelW1_3 SelW2_1 SelW2_2 TR1_1 TR1_2 TR1_3 TR1_4 TR2_1 TR2_2 TR2_3 TR2_4 TR3_1 TR3_2 TR3_3 SECp43";

for query in $selenoprot;
do {
dquery=`grep $query /home/U73966/luke/querylist.tab | cut -f 2`;
echo $dquery;
blastall -p tblastn -i $dquery -d /cursos/BI/genomes/project_2014/Pteropus_alecto/genome.fa -o ./blast/$query.tblastn -m8;
} done


#############################AUTOMATITZACIÓ DE LA SELECCIÓ DE HITS, FASTASUBSEQ, EXONERATE I T_COFFE##E#####################

pathGen=/cursos/BI/genomes/project_2014/Pteropus_alecto
echo quina selenoproteïna vols?

read selenoproteines #"Sel15_1 Sel15_2 Sel15_3 EEFSec_1 EEFSec_2 GPx1_1 GPx1_2 GPx2_1 GPx3_1 GPx3_2 GPx3_3 GPx3_4 GPx4_1 GPx4_2 GPx5_1 GPx5_2 GPx6_1 GPx6_2 GPx7_1 GPx8_1 GPx8_2 GPx8_3 DI1_1 DI1_2 DI1_3 DI1_4 DI1_5 DI1_6 DI1_7 DI2_1 DI2_2 DI2_3 DI2_4 DI2_5 DI3_1 MsrA_1 MsrA_2 MsrA_3 MsrA_4 PSTK_1 SBP2_1 SBP2_2 SecS_1 SecS_2 SPS1_1 SPS1_2 SPS1_3 SPS1_4 SPS2_1 SelH_1 SelH_2 SelI SelK SelM_1 SelM_2 SelN_1 SelN_2 SelO SelP_1 SelP_2 SelP_3 SelP_4 SelP_5 SelR1_1 SelR1_2 SelR1_3 SelR1_4 SelR2 SelR3_1 SelR3_2 SelR3_3 SelR3_4 SelS_1 SelS_2 SelS_3 SelS_4 SelS_5 SelT_1 SelT_2 SelT_3 SelT_4 SelU1_1 SelU1_2 SelU1_3 SelU1_4 SelU2_1 SelU2_2 SelU3_1 SelU3_2 SelV SelW1_1 SelW1_2 SelW1_3 SelW2_1 SelW2_2 TR1_1 TR1_2 TR1_3 TR1_4 TR2_1 TR2_2 TR2_3 TR2_4 TR3_1 TR3_2 TR3_3 SECp43"#


for selenoprot in $selenoproteines
do

COUNTER=1

echo $selenoprot
while read id_query id_subject ident alig_lenght mm gap q_start q_end s_start s_end evalue bit_score
do
    echo $COUNTER
    conte=`echo $evalue | grep "e-"`
    if [ $conte ]; then
	echo $evalue
	initial_length=$(( $s_end - $s_start ))
	if [ "$initial_length" -lt '0' ]
	then length=$(expr $initial_length \* '-1'); fi
        start=$s_start-50000
	total_length=100000+$length
	export PATH=/cursos/BI/soft/exonerate/i386/bin:$PATH
	fastafetch $pathGen/genome.fa genome.index $id_subject > ./fastafetch/$selenoprot.$COUNTER.fa 
	fastasubseq ./fastafetch/$selenoprot.$COUNTER.fa $start $total_length > ./fastasubseq/$selenoprot.$COUNTER.genome.fa
	./canviUperX.pl < ./queries/$selenoprot.fa > ./queriescanviades/$selenoprot.fa
	exonerate -m p2g --showtargetgff -q ./queriescanviades/$selenoprot.fa -t ./fastasubseq/$selenoprot.$COUNTER.genome.fa -E > ./GFF/$selenoprot.$COUNTER.gff
	egrep -w exon ./GFF/$selenoprot.$COUNTER.gff > ./exonsGFF/$selenoprot.$COUNTER.exons.gff;
	echo ./exonsGFF/$selenoprot.$COUNTER.exons.gff;
	export PATH=/cursos/BI/bin:$PATH
	./Fastaseq.gff.pl ./fastasubseq/$selenoprot.$COUNTER.genome.fa ./exonsGFF/$selenoprot.$COUNTER.exons.gff > ./cdna/$selenoprot.$COUNTER.cdna;
	echo ./cdna/$selenoprot.$COUNTER.cdna;
	fastatranslate ./cdna/$selenoprot.$COUNTER.cdna -F 1 > ./translate/$selenoprot.$COUNTER.translate.fa;
	echo ./translate/$selenoprot.$COUNTER.translate.fa;
	t_coffee ./translate/$selenoprot.$COUNTER.translate.fa ./queries/$selenoprot.fa 
	COUNTER=$((COUNTER + 1))
    fi
done <./blast/$selenoprot.tblastn
done


############################################AUTOMATITZACIÓ PEL GENEWISE I T_COFFEE##########################################

echo digues les selenoproteïnes

read selenoproteines

echo digues el num

read COUNTER

for selenoprot in $selenoproteines; do


#Ara entrem a fer el Genewise. 1r posem permisos, després fem el geniwise per la reverse (-trev), guardem una versió completa i una versió amb només la proteïna (.pep). Com que apareixen frases sobre la seqüència s'eliminaran amb la comanda sed -i. Finalment ja podrem fer el t-coffee. És fa el mateix amb la forward

    export PATH=/cursos/BI/bin:$PATH
    export WISECONFIGDIR=/cursos/BI/soft/wise-2.2.0/wisecfg

    genewise -trev -pep -pretty -cdna -gff ./queriescanviades/$selenoprot.fa ./fastasubseq/$selenoprot.$COUNTER.genome.fa > ./genewisecomplet/$selenoprot.$COUNTER.R.gff
    genewise -trev -pep ./queriescanviades/$selenoprot.fa ./fastasubseq/$selenoprot.$COUNTER.genome.fa > ./genewise/$selenoprot.$COUNTER.R.pep
    sed -i '/Making/d' ./genewise/$selenoprot.$COUNTER.R.pep
    t_coffee ./genewise/$selenoprot.$COUNTER.R.pep ./queries/$selenoprot.fa > ./genewise/$selenoprot.$COUNTER.R.aln 

    genewise -pep -pretty -cdna -gff ./queriescanviades/$selenoprot.fa ./fastasubseq/$selenoprot.$COUNTER.genome.fa > ./genewisecomplet/$selenoprot.$COUNTER.F.gff
    genewise -pep ./queriescanviades/$selenoprot.fa ./fastasubseq/$selenoprot.$COUNTER.genome.fa > ./genewise/$selenoprot.$COUNTER.F.pep
    sed -i '/Making/d' ./genewise/$selenoprot.$COUNTER.F.pep
    t_coffee ./genewise/$selenoprot.$COUNTER.F.pep ./queries/$selenoprot.fa > ./genewise/$selenoprot.$COUNTER.F.aln
done