En esta sesión veremos cómo usar Bioconductor para procesar y normalizar datos de microarreglos Affymetrix. Para ello necesitamos tener instalados Bioconductor y los siguientes paquetes:
Pueden revisar la sesión Instalando R y Bioconductor para saber como instalar paquetes adicionales de Bioconductor.
Vamos a analizar datos de expresión en microarreglos Affymetrix tomados del siguiente artículo:
Requirement of bic/microRNA-155 for normal immune function. Rodriguez A, Vigorito E, Clare S, Warren MV, Couttet P, Soond DR, van Dongen S, Grocock RJ, Das PP, Miska EA, Vetrie D, Okkenhaug K, Enright AJ, Dougan G, Turner M, Bradley A Science. 2007;316;608-11. PMID: 17463290 (pubmed)
En este experimento se usaron microarreglos Mouse430_2 de Affymetrix para comparar células T-helper (Th1) wildtype contra células en las que miR-155 había sido deletado. El propósito del experimento es encontrar los genes cuya expresión se incrementa al perder la regulación negativa de miR-155 (blancos directos).
Antes de empezar, necesitamos bajar el siguiente archivo: mir155.zip. Guarden el archivo en su Escritorio y descomprímanlo.
Estaremos utilizando los datos crudos, que son archivos .CEL de Affymetrix y representan a las imagenes tomadas de los 6 microarreglos:
Además de estos archivos, necesitamos una tabla que contenga la relación entre el nombre de cada archivo y las características (tipo celular, genotipo, replica, etc) de cada muestra. Para ello se genera un archivo con el siguiente formato:
| filename | cell | type | name |
| th1_ko1.cel | th1 | ko | ko1 |
| th1_ko2.cel | th1 | ko | ko2 |
| th1_ko3.cel | th1 | ko | ko3 |
| th1_wt1.cel | th1 | wt | wt1 |
| th1_wt2.cel | th1 | wt | wt2 |
| th1_wt3.cel | th1 | wt | wt3 |
Podemos agregar todas las columnas que queramos, si creemos que la información puede ser útil en algún momento. En la carpeta que bajaron, se incluye un archivo pdata.txt que contiene la tabla anterior.
Primeramente necesitamos cargar las librerías o paquetes necesarios para realizar todo el análisis. Cada paquete nos da acceso a funciones nuevas dentro de R.
library(affy) library(pvclust) library(vsn)
Ahora nos vamos a cambiar al directorio donde tenemos guardados los archivos CEL de Affymetrix. Para cambiar de directorio de trabajo, podemos hacerlo vía los menús (busquen la opción Change Working Directory... en el menú File o Misc). También podemos hacerlo directamente en R si conocemos la ruta completa:
setwd("/Users/cei/Desktop/mir155")
Para saber si funcionó, o en general si necesitamos saber nuestro directorio de trabajo actual:
getwd()
[1] "/Users/cei/Desktop/mir155"
Y finalmente, para asegurarnos de que estamos en la carpeta correcta, veamos los archivos que ahí se encuentran:
list.files()
[1] "pdata.txt" "th1_ko1.cel" "th1_ko2.cel" "th1_ko3.cel" "th1_wt1.cel" "th1_wt2.cel" "th1_wt3.cel"
Si podemos ver esos 7 archivos en el listado, quiere decir que estamos en la carpeta correcta y podemos continuar. El primer paso es leer el archivo pdata.txt:
pd = read.table("pdata.txt", header=TRUE, as.is=TRUE)
pd
filename cell type name
1 th1_ko1.cel th1 ko ko1
2 th1_ko2.cel th1 ko ko2
3 th1_ko3.cel th1 ko ko3
4 th1_wt1.cel th1 wt wt1
5 th1_wt2.cel th1 wt wt2
6 th1_wt3.cel th1 wt wt3
Nota: header=TRUE indica que el primer renglón corresponde a los encabezados y as.is=TRUE hace que las columnas con texto se lean correctamente (como character y no como factor).
Recuerden que hay varias maneras de accesar las columnas de una tabla:
pd[,3] pd$type
[1] "ko" "ko" "ko" "wt" "wt" "wt" [1] "ko" "ko" "ko" "wt" "wt" "wt"
Ahora ya podemos leer los datos contenidos en los archivos CEL, usando la función ReadAffy que es parte del paquete affy:
affyData = ReadAffy(filenames=pd$filename)
Podemos agregar a este objecto la información que guardamos en pd. Para ello hay una función que nos permite accesar el PhenoData del objeto que creamos:
pData(affyData) = pd
También podemos especificar un nombre corto que nos facilite identificar cada muestra:
sampleNames(affyData) = pd$name
Podemos revisar el contenido de affyData, para enterarnos que contiene información de 6 muestras, 45101 genes y que la anotación (automáticamente asignada por Bioconductor) es del chip mouse4302.
affyData
AffyBatch object size of arrays=1002x1002 features (19 kb) cdf=Mouse430_2 (45101 affyids) number of samples=6 number of genes=45101 annotation=mouse4302 notes=
Con esto sabemos que hemos cargado exitósamente toda la información del experimento.
Una fase muy importante del análisis de los datos de cualquier experimento de alto rendimiento es el análisis de la calidad. Existen varios tipos de gráficas que podemos generar para observar las tendencias de los valores, en este caso de expresión. Por ejemplo:
boxplot(affyData, col=rainbow(6)) hist(affyData, col=rainbow(6))
Nota: la función rainbow() genera colores de acuerdo al orden en un arcoiris. Aquí queremos 6 colores, pues tenemos 6 muestras.
Nota: podemos abrir mas ventanas para graficar, para poder comparar las figuras. En Windows se usa windows(), en Mac quartz(), en Linux x11().
Con estas figuras estamos tratando de ver si alguna de las muestras claramente se comporta diferente al resto. En este caso todo se ve bien.
Puede resultar conveniente visualizar las intensidades directamente a como aparecen sobre el microarreglo. A veces se pueden ver distorsiones técnicas que pueden explicar el comportamiento anómalo de uno de las muestras. Precaución: esta instrucción puede ser muy lenta, mejor utilizarla con una muestra a la vez.
image(affyData[,3]) # Calcula una imagen para la 3ra muestra
Si quieren ver algunos efectos muy extraños, revisen los siguientes ejemplos.
Ahora observaremos como correlacionan las intensidades entre las muestras. Esperamos que las muestras wt correlacionen mejor entre si, que con las muestras ko, y vice versa.
Podemos hacerlo con un heatmap:
heatmap(cor(exprs(affyData)), symm=T)
o con un dendrograma:
corClust = pvclust(exprs(affyData), nboot=1, method.dist="correlation") plot(corClust)
Aunque parece que una de las muestras (ko1) tiene un comportamiento extraño, las otras pruebas de calidad no mostraron que fuera problemático. En estos casos puede convenir quedarnos con todas las muestras y repetir algunos análisis de calidad despues de la normalización.
Otra herramienta popular para ver cómo se relacionan las muestras es un "Análisis de Componentes Principales". Podemos realizar una gráfica de PCA en R de la siguiente manera:
pca = princomp(exprs(affyData)) plot(pca$loadings, main="Principal Component Analysis", col=rainbow(6), pch=19, cex=2) text(pca$loadings, colnames(exprs(affyData)), pos=3, cex=0.8)
Con este análisis podemos ver que ko1 se agrupa correctamente con los demás ko.
Vamos a normalizar los microarreglos usando el algoritmo RMA (Robust Microarray Average). Esto va a crear un nuevo objeto eset que contiene los valores normalizados de expresión:
eset = rma(affyData)
Hay una gran variedad de métodos para normalizar disponibles en Bioconductor. Una lista de los métodos comunes para datos de Affymetrix se pueden obtener con la siguiente función:
normalize.AffyBatch.methods()
Hasta podemos "inventar" diferentes estrategias de normalizado usando la función expresso. Con ello podemos elegir cada una de las partes del proceso de normalizado. Sin embargo, es mucho más lento y consume mucha memoria, haciéndolo impractico para >10 chips.
eset = expresso(affyData, bgcorrect.method="rma", normalize.method="quantiles", pmcorrect.method="pmonly", summary.method="medianpolish")
Después de normalizar la distribución de los valores de expresión debe ser muy similar para todas las muestras (o idéntica para algunas variedades de normalización por cuantiles). Podemos comparar los boxplots antes y después de normalizar:
par(mfrow=c(1,2)) # Una misma ventana con multiples-figuras por "row", con 1 file y 2 columnas boxplot(affyData, col=rainbow(6)) boxplot(data.frame(exprs(eset)), col=rainbow(6)) par(mfrow=c(1,1)) # Regresar a una figura por ventana
También podemos ver si el agrupamiento de las muestras mejora después de normalizar:
par(mfrow=c(1,2)) corClust = pvclust(exprs(affyData), nboot=1, method.dist="correlation") plot(corClust, main="Agrupamiento de muestras antes de normalizar") corClustAfter = pvclust(exprs(eset), nboot=1, method.dist="correlation") plot(corClustAfter, main="Agrupamiento de muestras despues de normalizar") par(mfrow=c(1,1))
Una buena estrategia de normalizado debe controlar la variabilidad a lo largo del rango de valores de expresión. Por ejemplo, no queremos que los genes altamente expresados cambien en promedio más que los genes de expresión media. Podemos revisar esto con la siguiente gráfica:
meanSdPlot(exprs(eset))
Una muy buena idea es revisar algunos genes individuales para ver cómo se comportan, particularmente si tenemos controles. Dado que este experimento contempla la deleción de un gen, podemos observar la expresión del mismo, cuya sonda tiene el identificador 1428027_at:
boxplot(data.frame(exprs(eset)), col="grey"); lines(exprs(eset)["1428027_at",], lwd=2, type="b", col="red")
Podemos observar el comportamiento esperando para esta sonda, con niveles muy bajos en todas las muestras tipo ko y niveles mucho más altos en las muestras wt.
Finalmente, podemos guardar los valores completos de expresión, normalizados, para análisis subsecuentes:
write.exprs(eset, file="expr_normalizada.txt")
También podemos guardar el objeto de R, lo cual facilitará su uso en sesiones futuras de R.
save(eset,file="eset.Rdata")