Introducción a R - Ejercicios con datos genómicos
Ejercitándonos con datos genómicos
Como se empezaron a dar cuenta al final del ejercicio pasado, R permite manipular información que viene en tablas de manera muy flexible. Veamos entonces un caso más real.
Descarguen el archivo: ensembl_info.tab
La información en la tabla que acaban de descargar fue obtenida de la página de Ensembl, que reune una gran cantidad de información sobre genomas completamente secuenciados de vertebrados y algunos organismos modelo.
Carguemos la tabla dentro de R, para empezar a explorar su contenido.
ensTab = read.table("ensembl_info.tab")Ciertas caracterÃsticas de una tabla nos permiten hacer una exploración rápida. Por ejemplo sus dimensiones (renglones y columnas), o ver solo los primeros renglones:
dim(ensTab)## [1] 25 10head(ensTab)## coding_genes ncrna_genes pseudo_genes coding_gene_avg_length
## Anole_lizard 18595 7168 157 35332.78
## C_elegans 20362 922 1658 3088.45
## Cat 19493 1855 542 36566.75
## Chicken 18346 6490 43 24355.34
## Chimpanzee 18759 8681 572 53191.61
## Cow 19994 3825 797 39145.97
## ncrna_gene_avg_length pseudo_gene_avg_length
## Anole_lizard 10318.72 878.24
## C_elegans 334.04 1389.75
## Cat 103.16 792.30
## Chicken 12638.01 1012.88
## Chimpanzee 125.81 901.91
## Cow 111.60 859.42
## coding_trans_avg_length ncrna_trans_avg_length
## Anole_lizard 2572.26 548.16
## C_elegans 1399.78 255.94
## Cat 2169.50 103.16
## Chicken 2073.08 859.23
## Chimpanzee 2541.03 125.81
## Cow 2046.76 111.60
## pseudo_trans_avg_length genome_length
## Anole_lizard 868.57 1799143587
## C_elegans 855.84 100286401
## Cat 743.36 2455541136
## Chicken 993.67 1230258557
## Chimpanzee 883.45 3309577922
## Cow 815.04 2670422299Como arriba vemos que los nombres de renglones son informativos, al tener los nombres de las especies, podemos extraerlos todos:
rownames(ensTab)## [1] "Anole_lizard" "C_elegans" "Cat"
## [4] "Chicken" "Chimpanzee" "Cow"
## [7] "Dog" "Dolphin" "Duck"
## [10] "Elephant" "Fruitfly" "Fugu"
## [13] "Gorilla" "Horse" "Human"
## [16] "Microbat" "Mouse" "Panda"
## [19] "Pig" "Platypus" "Rat"
## [22] "Stickleback" "Tasmanian_devil" "Xenopus"
## [25] "Zebrafish"También existe una función colnames(), aunque en este caso no resulta tan interesante.
- Cómo obtenemos los datos correspondiente al genoma humano (el renglón)?
Ya sabemos que podemos especificar la posición numérica (o Ãndice) de los datos que queremos, pero si los elementos que queremos extraer tienen nombre es mucho más fácil usarlo:
ensTab["Human",]## coding_genes ncrna_genes pseudo_genes coding_gene_avg_length
## Human 22097 15502 16146 62651.06
## ncrna_gene_avg_length pseudo_gene_avg_length coding_trans_avg_length
## Human 14626.52 3763.29 3506.79
## ncrna_trans_avg_length pseudo_trans_avg_length genome_length
## Human 647.39 814.55 3096649726- Y si simplemente queremos el tamaño del genoma humano?
Podemos especificar tanto el nombre del renglón como el nombre de la columna que necesitan:
ensTab["Human","genome_length"]## [1] 3096649726Un número bastante grande. Quizá sea más fácil de visualizar si lo representamos como el tamaño pero en millones de bases (normalmente hablamos de megabases o Mb):
round(ensTab["Human","genome_length"]/1e6,1)## [1] 3096.6Es decir, son más de 3 mil Mb.
Pueden interpretar qué hace cada parte del código anterior? Experimenten si no.
Ejercicios:
- Les interesa explorar más de cerca el tamaño de los genomas. Generen un nuevo objeto que tenga solamente los valores de los tamaños de los genomas, expresados en Megabases, y manteniendo el nombre de cada genoma. Demuestren que es útil:
- Ordenen el resultado para ver fácilmente tanto los valores más pequeños como los más grandes, y sobreescriban el objeto original para preservarlo ordenado. Es el resultado que ustedes esperarÃan?
- Creen que dentro de las especies de la tabla, el humano es el organismo de mayor complejidad? DeberÃa esto estar reflejado en el tamaño su genoma?
- Cómo calculan cuántos genomas son más grandes que el genoma humano?
- Conviertan todos los tamaños de genomas a porcentajes del genoma humano. Es decir, el genoma humano ahora va a medir 100, uno que mida la mitad será 50 y uno que mida tan solo la décima parte medirá 10.
Ordenando un tabla completa
Al terminar los ejercicios anteriores se podrán dar cuenta que aunque el nuevo objeto tiene un orden más útil (ordenado del genoma más chico al más grande), la tabla con todos los datos sigue teniendo un orden quizá menos interesante (ordenado alfabéticamente). La función sort es muy útil para ordenar vectores (valores de un mismo tipo, en una sola dimensión). Sin embargo, verán que sort no funciona con la tabla completa.
Para estos casos existe una función diferente llamada order. Veamos un ejemplo muy sencillo:
x = c("hongo","bacteria","planta","virus")
y = c(100, 10, 1000, 1)
x[order(y)]## [1] "virus" "bacteria" "hongo" "planta"y[order(y)]## [1] 1 10 100 1000Ejercicio:
- Usen la función
orderpara reordenar la tablaensTabde acuerdo a la longitud de los genomas (y guárdenlo asÃ).
## coding_genes ncrna_genes pseudo_genes coding_gene_avg_length
## C_elegans 20362 922 1658 3088.45
## Fruitfly 13918 2832 257 6900.90
## Fugu 18523 703 162 7910.66
## Stickleback 20787 1617 52 9250.29
## Duck 15634 567 249 22641.10
## Chicken 18346 6490 43 24355.34
## ncrna_gene_avg_length pseudo_gene_avg_length
## C_elegans 334.04 1389.75
## Fruitfly 1650.93 1054.58
## Fugu 106.32 1173.02
## Stickleback 120.66 1038.38
## Duck 106.44 394.53
## Chicken 12638.01 1012.88
## coding_trans_avg_length ncrna_trans_avg_length
## C_elegans 1399.78 255.94
## Fruitfly 2350.62 820.22
## Fugu 1713.82 106.32
## Stickleback 1728.33 120.66
## Duck 1702.94 106.44
## Chicken 2073.08 859.23
## pseudo_trans_avg_length genome_length
## C_elegans 855.84 100286401
## Fruitfly 806.63 143725995
## Fugu 1109.69 393312790
## Stickleback 1034.13 461533448
## Duck 387.70 1105035747
## Chicken 993.67 1230258557Tamaño y composición de los genomas
Para continuar con nuestra discusión de tamaños de los genomas, observen las siguientes figuras. Qué pueden interpretar de ellas?
- Traten de generar figuras parecidas por su cuenta.
Si los genomas varÃan tanto de tamaño, pero el número de genes permanece más o menos constante, qué es lo que ocupa el resto del espacio en los genomas?
Para acercarnos a este dilema, trabajemos primero con solo un genoma. En este ejemplo voy a tomar el genoma humano, pero para seguir el ejercicio pueden usar cualquiera.
genome = ensTab["Human",]
genome## coding_genes ncrna_genes pseudo_genes coding_gene_avg_length
## Human 22097 15502 16146 62651.06
## ncrna_gene_avg_length pseudo_gene_avg_length coding_trans_avg_length
## Human 14626.52 3763.29 3506.79
## ncrna_trans_avg_length pseudo_trans_avg_length genome_length
## Human 647.39 814.55 3096649726- Cómo obtengo el espacio en el genoma ocupado por todos los genes codificantes?
## [1] 1384400473Nota: la columna coding_genes contiene el número de genes codificantes, mientras que coding_gene_avg_length contiene la longitud que en promedio ocupa esta clase de genes en el genoma
- Lo pueden expresar como porcentaje de la longitud total del genoma?
## [1] 44.7- Recuerden sus clases de biologÃa molecular. Todas las bases dentro de un gene, realmente forman parte del RNA mensajero? Qué proceso ocurre después de la transcripción para que el mRNA esté listo para ser traducido?
En la tabla que tienen, la columna coding_trans_avg_length se refiere precisamente a la longitud promedio del mRNA de dichos genes, no la del gen completo (que contiene también intrones).
- Si vuelven a calcular el espacio en el genoma ocupado por los genes codificantes, pero ahora usando correctamente las longitudes de los mRNA, qué obtienen?
## [1] 77489539- Lo pueden expresar como porcentaje de la longitud total del genoma?
## [1] 2.5- Se dan cuenta que ahora podrÃan calcular el espacio ocupado por todos los intrones en el genoma?
## [1] 1306910934- Lo pueden expresar como porcentaje de la longitud total del genoma?
## [1] 42.2- Finalmente, terminen este tipo de cálculos y generen una figura que represente el espacio ocupado en el genoma:
Ejercicio final
Al final de la sección anterior, probablemente tuvieron que escribir bastante código en R para calcular el espacio ocupado por distintos tipos de elementos en un genoma. Lo ideal serÃa poder aprovechar este código para poder obtener el resultado para cada uno de los genomas que están representados en la tabla.
Qué nos hace falta para poder hacer esto?
En principio solo falta saber cómo hacer un ciclo en R. Algunos sabrán hacer un ciclo for en algún lenguaje o en Linux. R también tiene este tipo de ciclos, aunque la sintaxis es un poquito diferente:
genomas = head(rownames(ensTab))
for (genoma in genomas) {
print(genoma)
}## [1] "C_elegans"
## [1] "Fruitfly"
## [1] "Fugu"
## [1] "Stickleback"
## [1] "Duck"
## [1] "Chicken"Ahora si, dentro de ese ciclo for tienen que acomodar el código que usaron anteriormente, para que genere una figura en cada ciclo. Qué podemos discutir sobre el resultado?
