1.17. Ejemplo práctico 1: Analizando el genoma humano
1.17.3. Análisis de los genes humanos
El análisis del catálogo de genes de una especie es un ejemplo perfecto de la utilidad de los comandos del terminal en este capítulo. Un gen es una secuencia de ADN que contiene información para crear una molécula biológica. Los genes de los organismos eucariotas se componen de exones. Muchos genes humanos tienen combinaciones plausibles de exones, lo que resulta en transcritos alternativos.
Para codificar la ubicación de los genes, se utilizan ficheros tabulados que contienen información como la localización y las características del transcrito del gen, como el número de exones y las coordenadas exactas.
En este ejercicio, se utiliza la anotación del consorcio RefSeq para obtener información sobre el genoma humano, en un primer momento, pero también se utilizará información del genoma de ratón (mouse). La información se encuentra en un fichero llamado refGene.txt y se puede obtener en la página web del navegador genómico de UCSC mediante el enlace Annotation. La dirección de descargas se muestra en la tabla 19.
Tabla 19. Páginas de descarga de los ficheros refGene.txt en función de la especie.
| Acceso | Dirección |
| Página descarga human | https://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/ |
| Página descarga mouse | https://hgdownload.soe.ucsc.edu/goldenPath/mm39/database/ |
Las anotaciones, en el servidor de UCSC, suelen representarse gráficamente en el navegador genómico en forma de pista. El navegador genómico está gestionado internamente por el administrador de bases de datos relacionales de MySQL. En consecuencia, en esta página web encontraremos dos tipos de fichero para cada pista de anotaciones de UCSC. El primer fichero, cuya extensión será del tipo sql, contiene una especificación genérica de los atributos de las anotaciones. Este archivo es necesario para crear una tabla vacía en una base de datos relacional. El segundo fichero, que estará comprimido y poseerá la extensión txt, contiene propiamente la información de cada anotación de forma tabulada. Al acceder a la página de descarga se visualiza la siguiente información:
This directory contains a dump of the UCSC genome annotation database for the Dec. 2013 (GRCh38/hg38) assembly of the human genome (hg38, GRCh38 Genome Reference Consortium Human Reference 38 (GCA_000001405.15)) affyGnf1h.sql 2015-05-11 01:50 2.1K affyGnf1h.txt.gz 2015-05-11 01:50 596K … refGene.sql 2020-08-17 18:56 1.9K refGene.txt.gz 2020-08-17 18:56 8.3M … xenoRefSeqAli.sql 2020-08-17 19:17 2.1K xenoRefSeqAli.txt.gz 2020-08-17 19:17 17M
Ahora procedemos a descargar el fichero txt asociado a la pista refGene, que contiene el catálogo de genes humanos anotados por el consorcio RefSeq. Para ello, se utiliza el comando wget nuevamente para transferir el fichero al terminal.
$ wget https://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/refGene.txt.gz
--2023-04-26 16:35:43-- https://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/refGene.txt.gz S'està resolent hgdownload.soe.ucsc.edu (hgdownload.soe.ucsc.edu)... 128.114.119.163 S'està connectant a hgdownload.soe.ucsc.edu (hgdownload.soe.ucsc.edu)|128.114.119.163|:443... conectat. HTTP: s'ha enviat la petició, s'està esperant una resposta... 200 OK Mida: 8668756 (8,3M) [application/x-gzip] S'està desant a: «refGene.txt.gz» refGene.txt.gz 100%[================>] 8,27M 2,22MB/s in 3,7s 2023-04-26 16:35:49 (2,22 MB/s) - s'ha desat «refGene.txt.gz» [8668756/8668756]
A continuación, te invitamos a que revises el contenido del primer fichero: refGene.sql. ya que es útil para conocer las características de los genes anotados en cada columna del fichero. Los atributos que se utilizan con mayor frecuencia son los siguientes: name (código del transcrito), chrom (cromosoma), strand (hebra), txStart y txEnd (coordenadas de inicio y final), exonCount (número de exones) y name2 (nombre del gen). Es importante no confundir los campos name y name2: un gen puede tener varios transcritos, pero un transcrito solo puede pertenecer a un gen. Estos parámetros se encuentran en la cabecera del fichero descargado.
Vamos a visualizar el contenido del fichero refGene.txt. Este archivo contiene información sobre el catálogo completo de genes anotados en el genoma humano. Cada línea representa un transcrito de un gen determinado. En el caso de que un gen tenga varios transcritos, cada uno se codifica en líneas separadas, cada una con su propio código y sus coordenadas correspondientes. En todas las líneas, los atributos de cada transcrito están separados por el carácter tabulador o «\t». Antes de poder visualizar el contenido del archivo, debemos descomprimirlo utilizando el comando gzip.
# Descompresión del fichero con gzip
$ gzip -d refGene.txt.gz
# Visualización de las primeras cinco líneas con el comando head
$ head -5 refGene.txt
585 NR_024540 chr1 - 14361 29370 29370 29370 11 14361,14969,15795,16606,16857,17232,17605,17914,18267,24737,29320, 14829,15038,15947,16765,17055,17368,17742,18061,18366,24891,29370, 0 WASH7P unk unk -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 585 NR_106918 chr1 - 17368 17436 17436 17436 1 17368, 17436, 0 MIR6859-1 unk unk -1, 585 NR_107062 chr1 - 17368 17436 17436 17436 1 17368, 17436, 0 MIR6859-2 unk unk -1, 585 NR_107063 chr1 - 17368 17436 17436 17436 1 17368, 17436, 0 MIR6859-3 unk unk -1, 585 NR_128720 chr1 - 17368 17436 17436 17436 1 17368, 17436, 0 MIR6859-4 unk unk -1,
Ten en cuenta que las anotaciones de un genoma suelen actualizarse frecuentemente. Por este motivo, los datos mostrados en este tutorial pueden variar ligeramente con el paso del tiempo si lo comparas con una nueva descarga del mismo fichero.
Debido a la gran cantidad de información que contiene este archivo son difíciles de manejar. A continuación, se decide trabajar con unos determinados campos y se utilizará el comando gawk para extraer únicamente los campos necesarios de este ejercicio. En particular, encontramos interesantes los siguientes atributos: nombre del gen, identificador del transcrito, cromosoma, hebra, coordenadas y número de exones. Para evitar tener que ejecutar el comando previo cada vez que se necesiten estos datos, se redirecciona la salida hacia un archivo refgene-select.txt. Una vez obtenido, se cuentan el número total de transcritos humanos.
# Analizamos el fichero refGene.sql para determinar la posición de cada columna
$ gawk '{print $13,$2,$3,$4,$5,$6,$9;}' refGene.txt > refgene-select.txt
# Se visualizan las últimas 5 líneas
$ tail -5 refgene-select.txt
KIAA0825 NM_001385728 chr5 - 94519215 94618604 6 KIAA0825 NM_001385729 chr5 - 94519215 94618604 6 KIAA0825 NR_169752 chr5 - 94519215 94618604 4 KIAA0825 NR_169753 chr5 - 94519215 94618604 6 KIAA0825 NR_169754 chr5 - 94519215 94618604 6
# El número de transcritos es el número de líneas del fichero
$ wc -l refgene-select.txt
88819 refgene-select.txt
Continuamos con diferentes análisis,
# Se solicita el número de transcritos genéticos distribuidos en la hebra positiva
$ gawk '{if ($4 == "+") print $0;}' refgene-select.txt | wc -l
45121
# Se solicita el número de transcritos del cromosoma 21
$ grep chr21 refgene-select.txt | wc -l
1121
# Se muestran los primeros siete transcritos después de ordenar por nombre del gen
$ sort refgene-select.txt | head -7
A1BG NM_130786 chr19 - 58345182 58353492 8 A1BG-AS1 NR_015380 chr19 + 58351969 58355183 4 A1CF NM_001198818 chr10 - 50799408 50885675 14 A1CF NM_001198819 chr10 - 50799408 50885675 15 A1CF NM_001198820 chr10 - 50799408 50885675 14 A1CF NM_001370130 chr10 - 50799408 50885627 12 A1CF NM_001370131 chr10 - 50799408 50885627 12
# Se muestran los primeros seis transcritos después de ordenar por el número de genes que contiene cada transcrito
$ sort -rnk 7 refgene-select.txt | head -6
TTN NM_001267550 chr2 - 178525990 178807423 363 TTN NM_001256850 chr2 - 178525990 178807423 313 TTN NM_133378 chr2 - 178525990 178807423 312 TTN NM_133437 chr2 - 178525990 178807423 192 TTN NM_133432 chr2 - 178525990 178807423 192 TTN NM_003319 chr2 - 178525990 178807423 191
El fichero refgene-select.txt contiene el listado de los transcritos humanos. Como un gen puede dar lugar a varios transcritos alternativos, podemos contar cuántos genes tiene el genoma humano y averiguar cuáles poseen un mayor número de transcritos. Para conseguirlo, ordenamos los nombres de los genes e introducimos diferentes variantes del comando uniq para agruparlos. Por otra parte, la mejor forma de estudiar el comportamiento de una línea de comandos es la deconstrucción: eliminando las últimas instrucciones se puede mostrar el resultado parcial de la ejecución en pantalla.
# Ordena la columna por el nombre de los genes
$ gawk '{print $1;}' refgene-select.txt | sort | more
A1BG A1BG-AS1 A1CF A1CF A1CF A1CF A1CF …
# Ordena la columna por el nombre de los genes y que sean valores únicos
$ gawk '{print $1;}' refgene-select.txt | sort | uniq | more
A1BG A1BG-AS1 A1CF A2M A2M-AS1 A2ML1 …
# Determina el número de genes únicos en el fichero
$ gawk '{print $1;}' refgene-select.txt | sort | uniq | wc -l
28307
# Determina cuántos transcritos hay para cada uno de los genes
$ gawk '{print $1;}' refgene-select.txt | sort | uniq -c | more
1 A1BG 1 A1BG-AS1 8 A1CF 4 A2M 3 A2M-AS1 2 A2ML1 1 A2MP1 …
# Determina cuántos transcritos hay para cada uno de los genes y ordenalos por el número de exones, de mayor a menor
$ gawk '{print $1;}' refgene-select.txt | sort | uniq -c | sort -rn
260 KIR2DS2 215 LOC101928804 177 KIR2DS4 144 MAP4 144 KIR3DS1 129 KIR2DL …
Otra operación muy frecuente consiste en el cálculo de valores promedio. En el siguiente ejemplo, el estudiante calculará el promedio de exones por transcrito y la longitud media de estos. El funcionamiento de gawk es similar en ambos casos, ya que trabaja con la variable t como contador que acumula la suma total. La división por el número de líneas visitadas (NR), una vez que se completa la lectura del archivo, genera el valor promedio en cada caso.
# Cálculo de promedio de exones por transcrito
$ gawk 'BEGIN{t=0;}{t=t+$7;}END{print t/NR;}' refgene-select.txt
10.0961
# Cálculo de la longitud media de los transcritos
$ gawk 'BEGIN{t=0;}{t=t+$6-$5+1;}END{print t/NR;}' refgene-select.txt
64883
Para ilustrar la utilidad de asociar dos ficheros de texto tabulado, se va a combinar el catálogo de genes humanos introducido hasta ahora (archivo refGene.txt para H. sapiens se renombra a refGene-human.txt) con el catálogo equivalente del genoma del ratón (archivo refGene.txt para M. musculus y renombrado como refGene-mouse.txt después de haberlo descargado del servidor genómico de UCSC).
$ gzip -d refGene.txt.gz
$ mv refGene.txt refGene-human.txt
$ wget https://hgdownload.soe.ucsc.edu/goldenPath/mm39/database/refGene.txt
# Se renombre el fichero del genoma M.musculus
$ mv refGene.txt refGene-mouse.txt
# Para cada uno de los ficheros, se seleccionan las columnas nombre del gen y cromosoma y posteriormente se visualizan los 5 primeros
$ gawk '{print $13, $3;}' refGene-human.txt | sort | uniq > refgene-select-human.txt
$ head -5 refgene-select-human.txt
A1BG chr19 A1BG-AS1 chr19 A1CF chr10 A2M chr12 A2M-AS1 chr12
# Para el fichero genoma de ratón
$ gawk '{print $13, $3;}' refGene-mouse.txt | sort | uniq > refgene-select-mouse.txt
$ head -5 refgene-select-mouse.txt
0610005C13Rik chr7 0610009B22Rik chr11 0610009E02Rik chr2 0610009L18Rik chr11 0610010F05Rik chr11
Con esta transformación previa del fichero, se han preparado los ficheros para utilizar el comando join: a continuación, se muestra cómo asociar los genes que tienen el mismo nombre en ambas especies. Al comando join se añade la opción -i, para que se ignoran las diferencias de mayúsculas/minúsculas. Cada línea del resultado final contiene el nombre del gen y su ubicación en ambos genomas.
$ join -i refgene-select-human.txt refgene-select-human.txt > refgene-comun.txt
$ head -5 refgene-comun.txt
A1BG chr19 chr19 A1BG-AS1 chr19 chr19 A1CF chr10 chr10 A2M chr12 chr12 A2M-AS1 chr12 chr12
A partir de este análisis preliminar del catálogo de genes humanos, se pueden poner en práctica múltiples variantes de los comandos aquí mostrados. Por ejemplo, es posible añadir más atributos, más genomas o más ficheros con otras propiedades para ampliar esta exploración. Por lo tanto, os animamos a experimentar con cada bloque de comandos utilizado durante este ejercicio. Con esta aplicación práctica, se demuestra la validez del terminal de Gnu/Linux para analizar eficientemente datos biológicos. En futuros apartados, se introducirán sistemas más complejos para gestionar grandes conjuntos de datos, de modo que el usuario podrá reproducir el mismo caso práctico utilizando dichas técnicas.