Mètodes
Per tal d’explicar el protocol emprat per realitzar el projecte, enumerarem punt per punt els passos que hem seguit:
- Exportació dels arxius necessaris per als diferents programes
- Cerca d’homologia mitjançant BLAST en els genomes problema
- Extracció de la regió que conté el millor hit
- Estudi de la estructura gènica del gen predit
- Conversió de l’hipotètic gen en cDNA i seqüència proteica
- Estudi de la proteïna predita i la query original
- Estudi d’un multialineament amb totes les proteïnes extretes
- Representació gràfica dels alineaments amb Jalview
- *. Elaboració de l’automatització del procés amb BASH
Cal comentar que per tal d’agilitzar tot el procés, davant de totes les vegades que s’ha de repetir, vam elaborar un programa BASH.
*. Elaboració de l’automatització del procés amb BASH
1. Exportació dels arxius necessaris per als diferents programes
Abans de començar hem d’exportar els arxius necessaris per tal de fer servir tota la programació que necessitarem. Farem servir les comandes següents:
Per exportar els arxius del BLAST:
$ export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH
$ cp /cursos/BI/bin/ncbiblast/.ncbirc.
Per tal d’exportar tot el paquet de l’Exonerate:
$ export PATH=/cursos/BI/bin/exonerate/i386/bin:$PATH
Pel FastaseqfromGFF, genewise i t_coffee:
$ export PATH=/cursos/BI/bin:$PATH $ export WISECONFIGDIR=/cursos/BI/soft/wise-2.2.0/wisecfg
2. Cerca d’homologia mitjançant BLAST en els genomes problema
Primer de tot, es farà una cerca d’homologia llençant les nostres querys sobre els diferents genomes. Abans de començar hem de fer una base de dades del nostre genoma, per fer-ho farem servir la següent comanda:
$ formatdb –i /cursos/BI/genomes/protists/2012/protista/genome.fa -p F -n nombbddBLAST.fa
En la comanda, protista fa referència al nom del organisme amb el que volem treballar i nombbddBLAST al nom que posarem a la base de dades del genoma. Quan ja hem creat la base de dades podem passar a fer el BLAST. En aquest cas, farem servir el tblastn, ja que volem llençar una query proteica sobre una seqüència nucleotídica, tenint en compte les 6 possibilitats de traducció d’aquesta. Per tal de portar-ho a terme farem servir la següent comanda:
$ blastall -p tblastn -i query.fa -d nombbddBLAST -o fitxerdesortida
En aquesta comanda tenim tota la informació de l’alineament. A més, si afegim la opció –m 9 a la comanda anterior podem obtenir una taula resum. Així, la comanda seria:
$ blastall -p tblastn -i query.fa -d nombbddBLAST -o fitxerdesortida –m 9
3. Extracció de la regió que conté el millor hit
Per tal d’estudiar el millor hit primer hem d’extreure de tot el fitxer multiFASTA la regió que el conté. Ja que el programa fastasubseq per tal d’extreure una regió concreta només funciona amb fitxers d’una seqüència, caldrà fer servir:
$ fastaindex /cursos/BI/genomes/2012/protista/genome.fa dm2.index
$ fastafetch /cursos/BI/genomes/2012/protista/genome.fa dm2.index nomseq > nomseq.fa
$ fastasubseq nomseq.fa nucleotidstart length > genomic.fa
On dm2.index fa referència a l’arxiu que crearem amb l’indexat. Amb aquest, podrem extreure mitjançant el fastafetch la regió que cont&eacutte; el millor hit, i aquesta que volem la indiquem amb nomseq. A més, del resultat del BLAST hem de mirar el nucleòtid on començarà la nostra regió i la llargada que volem (que ficarem en nucleotidstart i en length), pensant en que la regió que extraurem ha d’incloure el hit que volem estudiar. Es recomanable que per hits no gaire bons fiquem un marge més ampli, d’unes 10kb per cada extrem (en total, una length de 20kb).
4. Estudi de la estructura gènica del gen predit
Ara volem predir l’estructura gènica del gen de la nostra proteïna problema. Farem servir l’exonerate. Per fer-ho, el cridem de la següent manera:
$ exonerate -m p2g --showtargetgff -q query.fa -t genomic.fa
On query.fa fa referència a la proteïna inicial amb la qual fèiem el tblastn i genomic.fa a la regió genòmica que vam extreure amb el fastasubseq. En alguns casos el resultat és tan dolent que farem servir l’opció --exhaustive yes, per tal de forçar la predicció tot i que d’aquesta manera s’alenteix el procés.
5. Conversió de l’hipotètic gen en cDNA i seqüència proteica
Abans de poder obtenir la seqüència proteica hem d’aconseguir el cDNA del gen. Per fer-ho, inicialment hem d’obtenir els exons a partir de l’exonerate. Ho farem amb la següent comanda:
$ exonerate -m p2g --showtargetgff -q query.fa -t genomic.fa | egrep –w exon > exones.gff
Amb el pipe i l’egrep estem extraient només les coordenades dels exons. Amb aquestes dades, podrem fer servir el programa fastaseqfromGFF de la següent manera:
$ fastaseqfromGFF.pl genomic.fa exones.gff > cDNA.fa
En aquest pas, agafem les coordenades on es troben els exons i extraiem només aquesta regió de la zona genòmica que hem obtingut anteriorment. Així, descartem els introns.
Un cop tenim el cDNA volem obtenir la seqüència proteica derivada del nostre cDNA. Ho farem mitjançant la comanda:
$ fastatranslate –f cDNA.fa
Amb aquesta comanda el que farem serà obtenir les 6 pautes de lectura en les quals pot resultar la nostra seqüència de cDNA. MitjançCODIGOant l’estudi de les 6 pautes, s’escull la frameadeqüada i aquesta serà la nostra hipotèticaproteïna.
6. Estudi de la proteïna predita i la query original
Per comparar la proteïna que hem obtingut de tot el procés amb la original (la query), farem servir el programa t_coffee. Així, emprarem la següent comanda:
$ t_coffee fitxerFASTAsequencia1 fitxerFASTAsequencia2
On fitxerFASTA fa referència a ambdues seqüències que volem comparar.
7. Estudi d’un multialineament amb totes les proteïnes extretes
Per tal d’obtenir uns resultats més concloents, fem un multialineament amb totes les proteïnes que hem obtingut. Per tal de fer-ho, fem una concatenació de totes les proteïnes amb la query original i les guardem en un arxiu multiFASTA. Aquest serà analitzat amb t_coffee tal i com fèiem abans.
8. Representació gràfica dels alineaments amb Jalview
Finalment, emprarem el programa Jalview per obtenir una representació gràfica dels resultats. Com a input d’aquest programa hem de fer servir l’arxiu *.aln resultat de cada un dels anàlisis amb t_coffee. Amb aquest programa tindrem un alineament similar al de t_coffee, un gràfic amb el grau d’identitat de cada una de les posicions i a més, existeix la possibilitat de crear un arbre filogenètic de les seqüències obtingudes.
*. Elaboració de l’automatització del procés amb BASH
Per tal d’agilitzar el procés vam crear un programa amb BASH que automatitza el protocol abans explicat. El codi complet BASH del programa, completament desenvolupat pel grup 3B es detalla a continuació:
#!/bin/bash
#Exportacion de programas necesarios
export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH
cp /cursos/BI/bin/ncbiblast/.ncbirc ~/
export PATH=/cursos/BI/bin/exonerate/i386/bin:$PATH
export PATH=/cursos/BI/bin:$PATH
export WISECONFIGDIR=/cursos/BI/soft/wise-2.2.0/wisecfg
export PATH=/cursos/BI/bin/netblast/bin:$PATH
#Consideraciones
echo "."
echo "."
echo "AutoSEL v3.3-BETA-"
echo "BUG en traduccion, solo seleccionable las frames de la hebra guia"
echo "."
echo "Desarrollado por el grupo 3B - BioinformaticTEAM"
echo "."
echo "."
echo "INSTRUCCIONES"
echo "¡SOBRE TODO! Si tu selenoproteina contiene una U, cambiala por X"
echo "Por defecto, se cogen 20000 nucleotidos de margen. Cambios son notificados"
echo "Cuando aparezca el output del exonerate, debes mirar cuantos exones quieres extraer"
echo "En ocasiones, si el alineamiento es malo el exonerate no se realiza, se indicara que queremos 0 exones"
echo "Si al no encontrar buenos hits no se incluye el valor, quiere decir que no hay ninguno"
echo "Para seleccionar frame se muestran las 6 posibilidades, solo seleccionables las 3 ultimas. Para la reversa, se debera hacer a mano"
echo "En la carpeta t_coffees se incluyen los alineamientos individuales y un multialineamiento"
echo "."
echo "."
#Selección de query
echo "Introduce el nombre de tu query: (Solo el nombre, no pongas el .fa )"
read query
echo "."
echo "."
#Carpetas organizadoras
mkdir blast.coordenadas
mkdir blast.alineamiento
mkdir regiones.genomicas
mkdir output.exonerate
mkdir exones
mkdir cDNAs
mkdir proteins
mkdir t_coffees
#Proceso
for genome in F.cylindrus P.capsici A.laibachii_Nc14 G.niphandrodes I.multifiliis_strain_G5 S.arctica P.polycephalum D.fasciculatum D.discoideum_AX4 L.donovani_BPK282A1 L.tarentolae T.congolense C.fasciculata A.rara;
do {
echo "Procesando $genome"
echo "Blasteando..."
blastdb=`grep $genome /cursos/BI/genomes/protists/genomes_list_and_info.tab | cut -f 4`
blastall -F f -p tblastn -i $query.fa -d $blastdb -o $query.CONTRA$genome.tblastn -m 9
blastall -F f -p tblastn -i $query.fa -d $blastdb -o $query.CONTRA$genome.alineamiento.tblastn
sed '1,4d' $query.CONTRA$genome.tblastn > corte$genome
menor=`cut -f 11 corte$genome | head -1 | fold -1 | egrep -w - | wc -w`
if [ $menor -eq 0 ]
then
echo "No hay hits menores que 0. El mejor es de `cut -f 11 corte$genome | head -1`"
echo "."
echo "."
echo "Pulsa INTRO para continuar"
read continuar
echo "."
echo "."
mv $query.CONTRA$genome.tblastn ./blast.coordenadas
mv $query.CONTRA$genome.alineamiento.tblastn ./blast.alineamiento
rm corte$genome
else
evalue=`cut -f 11 corte$genome | head -1`
echo "E-value asignado como mejor de `cut -f 11 corte$genome | head -1`"
echo "Hit posicionado en `cut -f 2 corte$genome | head -1`"
echo "Extrayendo región..."
fastaindex $blastdb index$genome.index
fastafetch $blastdb index$genome.index "`cut -f 2 corte$genome | head -1`" > nomseq$genome.fa
rm index$genome.index
mv $query.CONTRA$genome.tblastn ./blast.coordenadas
mv $query.CONTRA$genome.alineamiento.tblastn ./blast.alineamiento
cp nomseq$genome.fa 1
longitud=`sed '1d' 1 | fold -1 | wc -w`
rm 1
hitstart=`cut -f 9 corte$genome | head -1`
hitend=`cut -f 10 corte$genome | head -1`
if [ $hitstart -gt $hitend ]
then
echo "REVERSE"
echo "El hit empieza en la posición $hitstart y acaba en $hitend"
hitstart=$hitend
hitend=`cut -f 9 corte$genome | head -1`
else
echo "El hit empieza en la posición $hitstart y acaba en $hitend"
fi
start=`expr $hitstart - 10000`
if [ $start -lt 0 ]
then
echo "Margen upstream de $hitstart nucleotidos"
start=1
fi
length=20000;
comprobar=`expr $start + 20000`
if [ $comprobar -gt $longitud ]
then
margendown=`expr $longitud - $hitend`
echo "Margen downstream de $margendown nucleotidos"
length=`expr $longitud - $start`
fi
fastasubseq nomseq$genome.fa $start $length > genomic$genome.fa
rm corte*
rm nomseq*
echo "Prediciendo estructura génica..."
echo "Evalue de $evalue, prefieres un exonerate exhaustivo? (si=1/no=0)"
read exhaustivo
if [ $exhaustivo -eq 1 ]
then
exonerate -m p2g --showtargetgff -q $query.fa -t genomic$genome.fa --exhaustive yes > exonerate.$query.$genome.gff
cat exonerate.$query.$genome.gff
echo "Extrayendo exones..."
egrep -w exon exonerate.$query.$genome.gff > exon.$query.$genome.gff
mv exonerate.$query.$genome.gff ./output.exonerate
else
exonerate -m p2g --showtargetgff -q $query.fa -t genomic$genome.fa > exonerate.$query.$genome.gff
cat exonerate.$query.$genome.gff
echo "Extrayendo exones..."
egrep -w exon exonerate.$query.$genome.gff > exon.$query.$genome.gff
mv exonerate.$query.$genome.gff ./output.exonerate
fi
echo "Escoge cuantos exones quieres extraer... (0 indica que no hay ninguno)"
read cuantos
if [ $cuantos -eq 0 ]
then
echo "No tienes exones, finalizando..."
mv genomic$genome.fa ./regiones.genomicas
mv exon.$query.$genome.gff ./exones
else
echo "Convirtiendo a cDNA..."
head -$cuantos exon.$query.$genome.gff > buenexon.$query.$genome.gff
fastaseqfromGFF.pl genomic$genome.fa buenexon.$query.$genome.gff > cDNA.$query.$genome.fa
rm buenexon.$query.$genome.gff
mv genomic$genome.fa ./regiones.genomicas
mv exon.$query.$genome.gff ./exones
echo "Traduciendo..."
fastatranslate -f cDNA.$query.$genome.fa > muestra.fa
cat muestra.fa
rm muestra.fa
echo "Selecciona la frame adecuada (OJO, solo puedes seleccionar las 3 directas) (1-3):"
read frame
fastatranslate -f cDNA.$query.$genome.fa -F $frame > protein.$query.$genome.fa
mv cDNA.$query.$genome.fa ./cDNAs
t_coffee protein.$query.$genome.fa $query.fa > t_coffee.$query.$genome
mv t_coffee.$query.$genome ./t_coffees
mv protein.$query.$genome.fa ./proteins
mv protein.$query.$genome.* ./t_coffees
fi
echo "Hecho"
echo "."
echo "."
echo "Pulsa INTRO para continuar"
read continuar
fi
} done
cat ./proteins/protein.$query.* > proteins.$query.MSA.fa
cat $query.fa >> proteins.$query.MSA.fa
t_coffee proteins.$query.MSA.fa > t_coffee.$query.MSA
mv t_coffee.$query.MSA ./t_coffees
mv proteins.$query.MSA.* ./t_coffees
mkdir resultados.$query
mv blast.* ./resultados.$query
mv regiones.genomicas ./resultados.$query
mv output.exonerate ./resultados.$query
mv exones ./resultados.$query
mv cDNAs ./resultados.$query
mv proteins ./resultados.$query
mv t_coffees ./resultados.$query
Els avantatges del nostre programa són que minimitza el nombre d’interaccions necessàries i per tant, també les probabilitats d’error al haver d’escriure alguna comanda o nom de seqüència.
A continuació es detalla pas per pas el procés portat a terme pel programa:
Primer de tot, s'haurà d'indicar el nom de la query que volem analitzar (sense el .fa). Seguidament, el programa iniciar´ la cerca a tots els genomes que s'hagin definit dins del for.
El primer pas que porta a terme autoSEL.pl serà realitzar un tblastn de la nostra query sobre els genomes, tant amb l'opció -m 9 com sense. Aquests dos arxius s'emmagatzemen dins de les carpetes blast.alineamiento i blast.coordenades, per tal de poder ser consultats més tard. El programa s'encarrega d'extreure i analitzar el millor hit obtingut. Primer de tot, mirarà si el primer hit es positiu o negatiu, i en el cas de ser positiu, detindrà la cerca en aquest genoma i pasarà al següent. Aquest pas ho farà mitjanÇant aquest fragment del programa:
sed '1,4d' $query.CONTRA$genome.tblastn > corte$genome menor=`cut -f 11 corte$genome | head -1 | fold -1 | egrep -w - | wc -w` if [ $menor -eq 0 ] then echo "No hay hits menores que 0. El mejor es de `cut -f 11 corte$genome | head -1`" echo "." echo "." echo "Pulsa INTRO para continuar" read continuar echo "." echo "." mv $query.CONTRA$genome.tblastn ./blast.coordenadas mv $query.CONTRA$genome.alineamiento.tblastn ./blast.alineamiento
Si no hi ha cap resultat al mostrar "el mejor es de ..." vol dir que no s'ha trobat cap hit, no es tracta de cap error
En el cas de trobar un hit positiu, el programa executarà el else corresponent al if del fragment anterior, que conté la resta del procés a portar a terme. Tot i que automàticament el programa agafarà el primer hit, mostrarà un output amb el valor agafat per tal de portar un control en tot moment del procés. Aquest pas es realitza de la següent manera:
else evalue=`cut -f 11 corte$genome | head -1` echo "E-value asignado como mejor de `cut -f 11 corte$genome | head -1`" echo "Hit posicionado en `cut -f 2 corte$genome | head -1`"
Un cop tenim el primer hit identificat, volem extreure la regió genònomica que conté el hit seleccionat. Primer identifiquem i extraiem el contig o regi´ on es troba el hit. Seguidament volem extreure només la regió del contig que conté específicament la zona alineada amb el hit, però deixant un marge (per defecte de 10kb) tant upstream com downstream. Abans de res, automàticament s'analitza si el hit está a l'hebra guia o si, en canvi, es tracta d'un hit en posició reversa. Per deixar el marge corresponent, el programa intenta donar les 10kb de marge, però alhora estudia si aquest marge es posible. En els casos en que el contig es massa petit o el hit es troba en algun dels extrems, el marge agafat es el màxim que ens permet el posicionament del nostre hit.
echo "Extrayendo región..."
fastaindex $blastdb index$genome.index
fastafetch $blastdb index$genome.index "`cut -f 2 corte$genome | head -1`" > nomseq$genome.fa
cp nomseq$genome.fa 1
longitud=`sed '1d' 1 | fold -1 | wc -w`
rm 1
hitstart=`cut -f 9 corte$genome | head -1`
hitend=`cut -f 10 corte$genome | head -1`
if [ $hitstart -gt $hitend ]
then
echo "REVERSE"
echo "El hit empieza en la posición $hitstart y acaba en $hitend"
hitstart=$hitend
hitend=`cut -f 9 corte$genome | head -1`
else
echo "El hit empieza en la posición $hitstart y acaba en $hitend"
fi
start=`expr $hitstart - 10000`
if [ $start -lt 0 ]
then
echo "Margen upstream de $hitstart nucleotidos"
start=1
fi
length=20000;
comprobar=`expr $start + 20000`
if [ $comprobar -gt $longitud ]
then
margendown=`expr $longitud - $hitend`
echo "Margen downstream de $margendown nucleotidos"
length=`expr $longitud - $start`
fi
fastasubseq nomseq$genome.fa $start $length > genomic$genome.fa
Per defecte, s'agafaran les 10kb, però en els casos en que s'hagi modificat apareixerà un missatge que ens notifica aquest canvi
Quan ja tenim preparada la regió que volem estudiar, pasem a portar a terme l'exonerate per tal de predir l'estructura génica. Com segons el valor del evalue, de vegades ens interesa que l'exonerate faci un estudi exhaustiu, abans d'iniciar la predició el programa ens torna a informar del valor de l'evalue i ens pregunta si volem o no fer l'exhaustiu. Hem de tenir en compte que el fet d'afegir aquest opció ens fa que es comprometi el temps de cerca, per la qual cosa, si el hit es un valor molt significatiu, no caldrà fer-ho exhaustiu i per tant, la predicció serà molt més ràpida.
echo "Prediciendo estructura génica..." echo "Evalue de $evalue, prefieres un exonerate exhaustivo? (si=1/no=0)" read exhaustivo if [ $exhaustivo -eq 1 ] then exonerate -m p2g --showtargetgff -q $query.fa -t genomic$genome.fa --exhaustive yes > exonerate.$query.$genome.gff cat exonerate.$query.$genome.gff echo "Extrayendo exones..." egrep -w exon exonerate.$query.$genome.gff > exon.$query.$genome.gff mv exonerate.$query.$genome.gff ./output.exonerate else exonerate -m p2g --showtargetgff -q $query.fa -t genomic$genome.fa > exonerate.$query.$genome.gff cat exonerate.$query.$genome.gff echo "Extrayendo exones..." egrep -w exon exonerate.$query.$genome.gff > exon.$query.$genome.gff mv exonerate.$query.$genome.gff ./output.exonerate fi
Amb l'output de l'exonerate veurem el nombre d'exons que ens ha predit, i ens demanarà quants volem extreure. En el cas de que l'exonerate no hagi predit res, podem indicar 0 exons per tal de detenir el procés. Després d'extreure l'exó el programa fastaseqfromGFF pasarà les coordenades a cDNA i fastatranslate pasarà aquest a proteïna. Abans de fer la traducció ens mostrarà tots els frames possibles i ens demanarà quin volem guardar. Amb la proteïna predita es farà un t_coffee per tal de poder estudiar el resultat i treure conclusions.
echo "Escoge cuantos exones quieres extraer... (0 indica que no hay ninguno)"
read cuantos
if [ $cuantos -eq 0 ]
then
echo "No tienes exones, finalizando..."
mv genomic$genome.fa ./regiones.genomicas
mv exon.$query.$genome.gff ./exones
else
echo "Convirtiendo a cDNA..."
head -$cuantos exon.$query.$genome.gff > buenexon.$query.$genome.gff
fastaseqfromGFF.pl genomic$genome.fa buenexon.$query.$genome.gff > cDNA.$query.$genome.fa
rm buenexon.$query.$genome.gff
mv genomic$genome.fa ./regiones.genomicas
mv exon.$query.$genome.gff ./exones
echo "Traduciendo..."
fastatranslate -f cDNA.$query.$genome.fa > muestra.fa
cat muestra.fa
rm muestra.fa
echo "Selecciona la frame adecuada (OJO, solo puedes seleccionar las 3 directas) (1-3):"
read frame
fastatranslate -f cDNA.$query.$genome.fa -F $frame > protein.$query.$genome.fa
mv cDNA.$query.$genome.fa ./cDNAs
t_coffee protein.$query.$genome.fa $query.fa > t_coffee.$query.$genome
mv t_coffee.$query.$genome ./t_coffees
mv protein.$query.$genome.fa ./proteins
mv protein.$query.$genome.* ./t_coffees
fi
} done
Nomès són seleccionables els frames directes, ja que no es coneix la comanda per extreure les reverses
Per acabar, es concatenen totes les proteïnes predites en un arxiu amb la query per tal de fer un multialineament. A més, per facilitar a l'usuari l'anàlisi dels seus resultats, s'organitzen tots els outputs en una jerarquia de carpetes molt entenedora.
cat ./proteins/protein.$query.* > proteins.$query.MSA.fa cat $query.fa >> proteins.$query.MSA.fa t_coffee proteins.$query.MSA.fa > t_coffee.$query.MSA mv t_coffee.$query.MSA ./t_coffees mv proteins.$query.MSA.* ./t_coffees mkdir resultados.$query mv blast.* ./resultados.$query mv regiones.genomicas ./resultados.$query mv output.exonerate ./resultados.$query mv exones ./resultados.$query mv cDNAs ./resultados.$query mv proteins ./resultados.$query mv t_coffees ./resultados.$query
Tot i això, el programa es troba a la seva versió 3.3 beta, i conté alguns comentaris que s’han de tenir en compte a l’hora de fer-ho servir. Totes les instruccions estan detallades al arxiu README.txt facilitat a la carpeta zip de descarrega del programa.
Clica el botó per tal de descarregar el programa autoSEL v3.3-BETA-
