Búsqueda de querys y genomas


Búsqueda de queries

Las selenoproteínas en las cuales se basa el trabajo han sido proporcionadas por los profesores de la asignatura de bioinformática. En este caso las selenoproteínas son SelR, SelN y SelT. Definimos como query aquella secuencia de aminoácidos conocida que creemos que tiene un posible relación con un genoma determinado (en nuestro caso un grupo de protistas) y, con los siguentes pasos, comprobamos esta relación. Para obtener sus secuencias proteicas accedemos a la base de datos SelenoDB, donde encontramos las distintas secuencias para cada uno de los organismos que presentan estas selenoproteínas. Estas secuencias las guardamos por separado en ficheros tipo FASTA que nombramos con esta nomenclatura:

selprot_query.fa

Por ejemplo,

selprotN_human.fa

Queries para SelR

Queries para SelT

Queries para SelN


Búsqueda de genomas

De la misma manera, el listado de los organismos también han sido proporcionados por los profesores de la asignatura. Los nombres e información relativa a cada uno de ellos se puede consultar en la sección de protistas de esta página. Sus genomas están contenidos en el siguiente archivo:

/cursos/BI/genomes/protists/2012/nombredelprotista/genome.fa

Donde se puede encontrar todo el genoma contenido de cada uno en archivos tipo multiFASTA.

tBLASTn

El BLAST (Basic Local Alignment Search Tool) es una herramienta informática que nos permite llevar a cabo alineamientos de tipo local de secuencias biológicas, tanto nucleotídicas como de proteínas. El programa compara una secuencia problema (query) con secuencias que se encuentran en bases de datos y calcula la significancia estadística de estas coincidencias. BLAST se puede usar para inferir las relaciones funcionales y evolutivas entre las secuencias, así como identificar miembros de familias de genes. El algoritmo que utiliza este programa permite encontrar regiones de similitud entre secuencias de forma muy rápida. No obstante, es un programa heurístico, de modo que no nos garantiza que la solución encontrada sea la correcta. Existen distintos tipos de BLAST según la naturaleza de la secuencia query y la secuencia de la base de datos.

En este proyecto, se ha utilizado el software tBLASTn, el cual permite comparar la secuencia proteica de la query con las secuencias de nucleótidos de los genomas de los protistas. Para ello, el programa traduce todas las secuencias de nucleótidos que hay en la base de datos de interés en los 6 posibles marcos de lectura y los compara con la secuencia proteica de la query.

En primer lugar, se extrae el software del NCBI BLAST mediante las siguientes comandas en la ventana del shell:

$ export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH

$ cp /cursos/BI/bin/ncbiblast/.ncbirc ~/

Seguidamente, se utiliza la herramienta informática formatdb que es necesaria para utilizar los archivos con secuencias en formato FASTA de los genomas de los protistas, de modo que se necesita formatear la base de datos para poder hacer el BLAST. Así pues, esta herramienta nos genera un archivo multiFASTA. Este proceso lo llevamos a cabo con la siguiente comanda:

$ formatdb -i /cursos/BI/genomes/protists/2012/nom_protista/genome.fa -p F -n genoma_protista.fa

Donde -i indica el directorio donde se encuentra el genoma del protista, -p F nos está formateando la base de datos y -n nos permite darle un nombre a la base de datos. Por lo tanto, mediante esta comanda, le estamos diciendo que coja el fichero FASTA donde se encuentra el genoma del protista y que nos haga una base de datos.

A continuación, ejecutamos una búsqueda tBLASTn de regiones de similitud entre la secuencia query contra una base de datos determinada mediante la comanda siguiente:

$ blastall -p tblastn -i query.fa -d genoma_protista.fa -o query.CONTRA.genoma_protista.fa

Donde -p indica el tipo de BLAST que queremos utilizar, -i es el archivo donde se ubica la secuencia query, -d hace referencia al archivo que se utiliza como base de datos y -o es el fichero donde se exportará la información de salida, es decir, donde el BLAST almacenará los resultados de la búsqueda.

En esta comanda del BLAST, también le hemos añadido el parámetro -m 8 y, de este modo, los resultados del BLAST aparecen en una tabla que contiene las características de cada uno de los hits que se obtienen en cada alineamiento, como las posiciones de inicio y fin de la query y del genoma, el E-value, etc.

$ blastall -p tblastn -i query.fa -d genoma_protista.fa -o query.CONTRA.genoma_protista.fa -m8

Así pues, el programa tBLASTn, para cada uno de los alineamientos, nos da como resultado una serie de hits, que son las regiones de los genomas de los protistas de la base de datos que presentan similitud con la secuencia de la selenoprotíena (query). Además, cada uno de estos hits posee un valor E-value que define la probabilidad de que un alineamiento determinado, con una secuencia perteneciente a una base de datos de un determinado tamaño, sea fruto del azar. Este valor es estadísticamente significativo cuando es menor o igual a 10-4. Cuando encontramos un hit significativo, miramos la presencia de selenocisteína en el alineamiento, de modo que la región de la secuencia query, alineada con una región del genoma protista, debe contener el aminoácido selenocisteína representado por una U. De este modo, podremos continuar el proceso con este hit para poder determinar si se encuentra en el genoma protista.

Con el objetivo de llevar a cabo este proceso anterior de forma más eficiente y obtener los resultados más rápidamente, hemos automatizado dicho procedimiento. Para ello, hemos elaborado el programa automaticblast.sh en bash.

Extracción de la región genómica

El objetivo de este análisis es obtener una región que sea de nuestro interés. Por ello, una vez conocemos los resultados significativos obtenidos mediante tBLASTn, que habrán sido guardados en archivos independientes, obtenemos un punto de inicio y extensión de la diferentes regiones alineadas. Usando esta información podemos extraer los fragmentos genómicos que son de nuestro interés, ya que es probable que contengan la proteína que estamos buscando.

Primero, debemos usar la comanda fastaindex. Este paso es necesario para poder realizar el siguiente paso (fastafetch) y se encarga de realizar un índice de los contigs que contiene la base de datos de cada organismo. Se requiere de la siguiente comanda:

$ fastaindex /directorio/genome.fa organismo.index

Donde organismo.index es un fichero donde están todos los índices de los contigs.

Como en el caso del BLAST, hemos realizado un programa para conseguir realizar esta tarea de forma automática: automaticindex.

Seguidamente, se puede realizar fastafetch que selecciona la región de interés (es decir, el contig deseado que conocemos gracias al análisis del tBLASTn) y lo almacena en un fichero nuevo de tipo FASTA, con el que sí que podremos trabajar para realizar los siguientes pasos. El comando necesario es:

$ fastafetch /directorio/genome.fa salida.index nombsec > nombsec.fa

Lo que estamos haciendo es indicarle una ubicación de genoma, seguidamente la ubicación del índice, nombsec es la región que queremos extraer, en nuestro caso un contig concreto. Finalmente, redireccionamos el resultado a nombsec.fa, siendo el fichero de salida.

Para este paso, también hemos diseñado un programa que nos permitirá automatizar este paso: automaticfastafetch

Una vez tenemos el fichero fasta con la región genómica de interés delimitamos más la zona de esta región, haciendo más fácil el trabajo con la secuencia. Para hacerlo, cogemos las posiciones obtenidas en el hit del BLAST y entonces expandir los márgenes tanto en upstream como en downstream para asegurarnos que nuestro gen se encontrará contenido en esta trozo. Usamos la comanda fastasubseq, que nos permitirá delimitar la región genómica un más concreta que con el contig total. La comanda es la siguiente:

$ fastasubseq /directorio/genome.fa start length > subseq_nombreorganismo.fa

Indicamos primero la ubicación del genoma en cuestión, start indica el nucleótido de la subsecuencia y length indica la largada. Finalmente, todo se redirige a un fichero de salida subsec_query_organismo.fa. Como ha pasado con los pasos anteriores, hemos podido automatizarlo mediante un programa.

Predicción de genes

De la extracción de la región genómica obtenemos una larga secuencia de 20000 nucleótidos. Comprobamos la presencia de exones dentro de los hits y así miramos si es un secuencia codificante o no.

Exonerate

Exonerate es uno de los programas que nos permite realizar este proceso, permitiéndonos obtener la secuencia del cDNA, ya que que nos permite saber las coordenadas de los exones y, así, la anotación completa del genoma de la selenoproteína.

Primero debemos llamar el programa, mediante el comando:

$ export PATH=/cursos/BI/bin/exonerate/i386/bin:$PATH

Para ejecutar el programa usamos el siguiente comando:

$ exonerate-m p2g -showtargetgff -q query.fa -t subseq_query_organismo.fa > query_organismo.gff

En este comando -m nos indica el tipo de alineamiento que queremos utilizar, seguido del argumento p2g que significa que se hará un alineamiento del tipo protein2genome que compara un secuencia proteica con una DNA. El argumento -showtargetgff nos permite obtener el resultado en formato GFF. Seguidamente, -q nos indica la ubicación de la secuencia query, mientras que -t nos indica la secuencia target, que será el resultado del fastasubseq. También podemos concatenar esta comanda con la opción grep para obtener solo las lineas de interés.

$ exonerate-m p2g -showtargetgff -q query.fa -t subseq_query_organismo.fa | grep -w "exon" > query_organismo.gff

De esta forma, sólo obtenemos aquellas líneas que contengan la palabra exon.

Hemos diseñado un programa que nos permitirá automatizar este paso: exonerate

Posteriormente debemos correr un programa en Perl dado por los profesores (fastaseqfromGFF.pl) que nos permitirá contruir el cDNA de la secuencia (únicamente los exones), por lo que hemos utilizado anteriormente la opción grep. Ejecutaremos el programa de la siguiente manera:

./fastaseqfromGFF.pl subseq_query_organismo.fa query_organismo.gff > query_organismo.cdna

Como podemos ver comparamos el resultado del fastasubseq con el resultado del Exonerate y lo redirigimos a un nuevo fichero que contendrá el cDNA de nuestra proteina predicha.

Como en los casos anteriores, hemos automatizado este proceso: automaticcdna.sh


GeneWise

GeneWise tiene la misma finalidad que el programa anterior, Exonerate. Nos permite obtener la secuencia de aminoácidos de la proteína predicha, el cDNA, el número de exones y el alineamiento de ésta con la query en la región delimitada.

A diferencia de Exonerate, GeneWise es sensible a la dirección del alineamiento. Así, hemos de ejecutar el programa teniendo en cuenta si es la forward strand o la reverse strand. Primero debemos llamar al programa:

$ export PATH=/disc8/bin:$PATH

$ export WISECONFIGDIR=/disc8/soft/wise-2.2.0/wisecfg

Las comandas para ejecutar el programa, para forward y reverse respectivamente, son:

$ genewise -pep -pretty -cdna -gff query.fa fastasubseqqueryVSprotista.fa > genewisequeryVSprotista.fa

$ genewise -pep -pretty -cdna -gff -trev query.fa fastasubseqqueryVSprotista.fa > genewisequeryVSprotistarev.fa

El programa GeneWise ha sido automatizado mediante el diseño en un programa en bash

De los datos obtenidos de GeneWise como en Exonerate obtenemos la secuencia nucleotídica que forma el cDNA de la secuencia predicha. Con lo que necesitamos ejecutar una comanda llamada fastatranslate que hará una traducción a aminoácidos en los seis marcos de lectura posibles.

$ fastatranslate query_organismo.cdna > query_organismo.mRNA

Usando esta comanda obtenemos los seis marcos de lectura que se redirigen a un fichero. Debemos escoger el mejor de ellos, descartando aquellos que presenten muchos * en su secuencia ya que significa que hay muchos codones STOP, y por lo tanto, no nos interesan. Para escoger solo uno usamos la opción -F y lo redirigimos a un fichero FASTA donde solo estará el marco de lectura elegido (X en la comanda):

$ fastatranslate query_organismo.cdna -F X > query_organismo.mRNA

T-COFFEE

T-COFFEE nos permite alinear las secuencias que hemos obtenido y nos da un resultado de su significancia. En nuestro caso, alineamos la query con la proteína predicha y el programa comprueba, mediante una aproximación progresiva, si estas secuencias alinean o no. La comanda utilizada es:

$ t_coffee query_organismo.mRNA > tcoffe_query_organismo.aln

Debido a que se ha obtenido el cDNA a partir tanto de Exonerate como de GeneWise, al realizar los programas para automatizar el proceso, se han diseñado dos: uno para Exonerate y para GeneWise.

BLASTp

Usamos BLASTp para comprobar que la proteína que hemos predicho mediante los pasos realizados anteriormente se corresponde con la query con la que hemos comparado. Usamos como query el resultado obtenido al realizar el fastatranslate. Este paso también nos permite ver si existen homólogos en las bases de datos del NCBI.

SECIS

Una vez hemos completado todos estos pasos y pensamos que podemos afirmar que nuestra proteína predicha es una selenoproteína, hemos de comprobar que realmente ese codón TGA codifique para una selenocisteína y que no son, en realidad, un codón STOP. Para ello, buscamos en nuestra query la presencia de elementos SECIS. Estos elementos son los encargados de recodificar el codón TGA a selenocisteína. Para ello, usamos el programa SECISsearch.

Para realizar la búsqueda de elementos SECIS hemos realizado el siguiente programa automático en bash.