En este caso, veremos una de las maneras más sencillas de obtener datos de expresión desde una base de datos que reune miles de experimentos ya publicados: el Gene Expression Omnibus (GEO) del NCBI.
Para muchas revistas es un requisito subir los datos de expresión generados para una publicación antes de ser aceptada. Esto ha generado que recursos como GEO sean invaluables para obtener datos para generar o confirmar hipótesis, o para explorarlos nuevamente a la luz de nuevos datos o preguntas.
En algunos casos, aunque no siempre, los datos originales o crudos están disponibles para ser analizados. Sin embargo, siempre deben subirse los datos ya analizados, que típicamente consiste en una tabla o matriz de expresión. Esta tabla consiste de genes (o transcritos o sondas) como renglones, y condiciones o muestras como columnas. Estas tablas pueden descargarse navegando en el sitio web, pero también pueden obtenerse usando el paquete GEOquery de Bioconductor. Suponiendo que no lo tienen instalado, es buen momento para aprender a instalar paquetes:
source("http://bioconductor.org/biocLite.R")
biocLite("GEOquery")
Con un poco de suerte, no se requerirá mayor esfuerzo. El Vignette que acompaña a GEOquery puede ser algo confuso al principio.
openVignette(package="GEOquery")
Recomiendo que le den una revisada, pero trabajemos primero con un ejemplo simple. Supongamos que explorando el sitio de GEO o leyendo un paper, nos topamos con un trabajo donde usaron microarreglos para analizar: SIRT1 deficiency effect on the brain. La entrada de la base de datos de GEO es accesible con el siguiente identificador: GDS4895. En esa liga nos podemos enterar del resumen del trabajo, la referencia de la publicación, la especie, el nombre del microarreglo usado, el número de muestras analizadas, etc.
Teniendo un identificador de GEO, podemos descargar los datos asociados con una simple función:
gdsId = "GDS4895" gds = getGEO(gdsId)
El objeto resultante puede ser algo intimidante. El Vignette de GEOquery lo describe con mayor detalle, y muestra algunos ejemplos de cómo accesar información contenida dentro. Pero aprovechemos la situación para describir una manera de obtener ayuda dentro de R para trabajar con un objeto de cierta clase. Primero determinemos la clase:
class(gds)
[1] "GDS" attr(,"package") [1] "GEOquery"
La siguiente es la manera en R para pedir ayuda sobre una clase en particular de objetos:
?"GDS-class"
Ojo con las comillas, para evitar que R crea que queremos hacer una resta. En este caso la ayuda puede no ser tan informativa, pero sugiere que la clase GDS extiende la clase GEOData, tal que podemos buscar más ayuda:
?"GEOData-class"
Ahora si podemos ver una lista de Methods que son funciones que podemos aplicar sobre este tipo de objetos.
Accession signature(object = "GEOData"): returns the GEO acccession for the current object Columns signature(object = "GEOData"): returns the column descriptions for the current object Meta signature(object = "GEOData"): returns the metadata for the current object Table signature(object = "GEOData"): returns the "Table" for the current object
En particular dos de ellos son bastante útiles:
Columns(gds)
sample genotype/variation description
1 GSM712769 SIRT1 null Value for GSM712769: BSKO 1; src: Brain of BSKO mice
2 GSM712798 SIRT1 null Value for GSM712798: BSKO 2; src: Brain of BSKO mice
3 GSM712800 SIRT1 null Value for GSM712800: BSKO 3; src: Brain of BSKO mice
4 GSM712802 SIRT1 null Value for GSM712802: BSKO 4; src: Brain of BSKO mice
5 GSM712797 wild type Value for GSM712797: WT 1; src: Brain of WT mice
6 GSM712799 wild type Value for GSM712799: WT 2; src: Brain of WT mice
7 GSM712801 wild type Value for GSM712801: WT 3; src: Brain of WT mice
8 GSM712803 wild type Value for GSM712803: WT 4; src: Brain of WT mice
head(Table(gds))
ID_REF IDENTIFIER GSM712769 GSM712798 GSM712800 GSM712802 GSM712797 GSM712799
1 1415670_at Copg1 -0.000136375 0.806718 0.000136 -0.0481071 0.106815 -0.0114985
2 1415671_at Atp6v0d1 0.0509319 -0.240773 0.0573463 0.119117 -0.0727654 -0.0301743
3 1415672_at Golga7 -0.114847 0.134076 -0.148633 -0.113399 -0.0276966 0.0763321
4 1415673_at Psph 0.0617933 0.0967522 -0.192098 -0.00198555 -0.0124388 -0.0176058
5 1415674_a_at Trappc4 0.0637188 -1.28701 -0.135022 0.0238371 -0.134872 -0.0238371
6 1415675_at Dpm2 -0.00965786 -0.726266 0.0884523 0.0903587 -0.154342 -0.0743914
Desafortunadamente el objeto sigue siendo un poco complicado de manipular. La mejora manera es convertirlo a otra clase mucho mas usado dentro del mundo de Bioconductor, que nos va a dar mucha flexibilidad, llamado el ExpressionSet.
eset = GDS2eSet(gds) class(eset)
[1] "ExpressionSet" attr(,"package") [1] "Biobase"
Podemos navegar las entrañas del sistema de ayudas, pero resumo algunas de las funciones más relevantes:
pData(eset) head(fData(eset)) head(exprs(eset))
Este objeto es mucho mejor para varios tipos de análisis:
boxplot(exprs(eset))
