###### tags: `UAX` `Estudiantes`
# Laboratorio IV: Análisis de expresión génica con microarrays
<font color = 'gray'>
<p style="text-align:right;">Prof.: Laura J. Marcos-Zambrano </font>
### Preparar la sesión de RStudio
#### Instalar paquetes:
```r
#Instalar paquetes de BioConductor:
BiocManager::install("affy")
BiocManager::install("oligo")
BiocManager::install("annotate")
BiocManager::install("genefilter")
BiocManager::install("hgu133plus2.db")
BiocManager::install("limma")
#Instalar paquetes de CRAN
install.packages("gplots")
```
>Si aparece el mensaje: "Update all/some/none? [a/s/n]:"
Teclear *a*
#### Cargar librerias
```r
library("tidyverse")
library("GEOquery")
library("affy")
library("oligo")
library("annotate")
library("genefilter")
library("hgu133plus2.db")
library("gplots")
library("limma")
```
### Datos del Seminario V
Usaremos el dataset **GSE33146** que exploramos en el [Seminario V](https://hackmd.io/@laurichi13/SkJs3fsXo).
Necesitamos tener el `ExpresionSet` ya que contiene todos los metadatos del experimento dentro del objeto `phenoData`.
```r
gse33146 <- getGEO('GSE33146')
length(gse33146)
class(gse33146[[1]])
#Guardar el experimento en un objeto
gse33146 <- gse33146[[1]]
```
Explorando los metadatos.
>Los metadatos, literalmente «sobre datos», son datos que describen otros datos.
`varLabels` permite ver la información de los datos incluidos en el objeto `phenoData`
```
varLabels(gse33146)
```
En nuestro ejemplo, la información de mayor interés es "culture medium:ch1" ya que sabemos que el experimento [GSE33146](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE33146) prueba diferentes medios de cultivo para ver la diferencia en la expresión de genes.
### Descargar datos de GEO en formato CEL
Los archivos `.CEL` tienen los datos originales sin procesar. :eye-in-speech-bubble: Como queremos procesarlo nosotros mismos, eso es lo que necesitaremos.
```r=1
#Descargar data en formato CEL
getGEOSuppFiles(GEO = "GSE33146", makeDirectory = TRUE)
```
Este comando creará un nuevo directorio llamado "GSE33146", que contiene un archivo comprimido de los datos sin procesar de GSE33146.
```r=3
untar("./GSE33146/GSE33146_RAW.tar", exdir = "./GSE33146")
```
El comando `untar` descomprimirá el archivo comprimido (.tar) en el directorio principal.
### Procesando un grupo de archivos CEL de Affymetrix.
```r=4
celfilelist <- list.files(path = "./GSE33146/", pattern = ".CEL", full.names = TRUE)
```
El comando `list.files` enumera todos los archivos en GSE33146 que tienen ".CEL" en el nombre de archivo. Es decir los datos crudos de todas las muestras de ese experimento (En nuestro caso de ejemplo son 6)
```r=5
cf <- read.celfiles(filenames = celfilelist,phenoData=phenoData(gse33146))
print(class(cf))
```
El comando `read.celfiles` lee la lista de archivos CEL (objeto `celfilelist`) y los almacenará en el objeto nuevo llamado cf del tipo `ExpressionFeatureSet`.
>:memo: Los archivos `CEL` continen los datos crudos, pero no siempre tinen todos los metadatos, por eso incluimos los metadatos de `gse33146`(Nuestro `ExpresionSet` completo) mediante el comando `phenoData=phenoData(gse33146)`
¿Qué hay en el objeto `ExpressionFeatureSet`?
```r=7
print(cf)
```
Este objeto almacena todas las muestras con datos CRUDOS y contiene información sobre la anotación apropiada para este tipo de chip ("pd.hg.u133.plus.2")
### Exploración de datos "crudos"
Histograma de intensidad antes de la normalización
```r=8
hist(cf)
```
Boxplot de intensidad antes de la normalización
```r=9
boxplot(cf)
```
### Normalización mediante RMA (Robust Multiarray Analysis)
```r=10
cf.norm <- rma(object = cf)
```
Es importante normalizar los arrays juntos si se tiene la intención de comparar valores entre arrays. Existen diferentes métodos de normalización (los vimos en clase). RMA se utiliza generalmente para anlizar arrays de Affimetrix.
:eyes: Volvemos a ver los gráficos después de la normalización.
```r=11
hist(cf.norm)
boxplot(cf.norm)
```
### Filtrado de valores
Filtrar las sondas de baja expresión y baja varianza puede mejorar las estadísticas posteriores. Con la función `genefilter()`, podemos eliminar todas las sondas que tienen una expresión inferior a log2(100) en más del 20 % de las muestras y tienen un rango intercuartílico inferior a 0.25. (Estos son los valores comúnmente usados en este tipo de análisis)
```r=13
#Creamos el dataset con las características deseadas
filter <- genefilter(cf.norm, filterfun(pOverA(p = 0.2, A = log2(100)), function(x) (IQR(x) > 0.25)))
#Aplicamos el filtro en nuetros datos
cf.norm <- cf.norm[filter,]
```
Aquí podemos ver el número de sondas que se eliminan y las que quedan después de aplicar el filtro (TRUE).
```r=17
table(filter)
```
### Comparaciones entre grupos
#### Un nuevo tipo de objeto en R: *"La clase fórmula"*
La clase fórmula es el objeto más importante para el modelado estadístico en R. Se utiliza `~` en la especificación del modelo, donde `y ~ x` significa que la respuesta *y* está modelada por una variable lineal *x*. (Es decir "y" depende de "x")
Los modelos lineales expresados de esta manera incluyen un término de "intersección" implícito (El intercepto), que podemos eliminar del modelo indicando + 0 en la fórmula de nuestro modelo.
>Por una variedad de razones, a menudo es conveniente eliminar el término intercepto cuando se realiza un análisis de expresión génica diferencial.
### Identificar características expresadas diferencialmente
En un experimento de transcriptómica, ya sea usando microarrays o RNA-seq, a menudo ajustamos los datos de cada transcripción como respuesta frente a un modelo común de variables independientes. Las variables describen el experimento y especificamos solo el lado derecho de una fórmula (el lado izquierdo se usa para ajustar los datos). Por ejemplo, si tenemos un tratamiento `treat` que puede tomar varios valores, y el valor real de `treat`varía entre las muestras, podríamos especificar un modelo de la siguiente forma:
`~ treat`
o `0 + treat` eliminando el intercepto.
(Querría decir que los datos dependen del tratamiento... )
### Especificando nuestro modelo para el análisis de expresión génica diferencial
Para identificar genes expresados diferencialmente usando modelos lineales, necesitamos hacer dos cosas:
1. Generar una **matriz modelo** especificando el diseño,
2. Ajustar el modelo al diseño.
3. En todos los casos, excepto en la comparación simple de dos clases, también debemos especificar una **matriz de contraste**, a la que llegaremos en breve.
En nuestro conjunto de datos (GSE33146), solo tenemos un tratamiento: **el cambio de condiciones de cultivo.**
Generamos nuestra matriz modelo:
```r=18
design <- model.matrix( ~ cf.norm$'culture medium:ch1')
colnames(design)[2] <- "SCGM"
```
>"cf.norm$'culture medium:ch1'" es el medio de cultivo usado en cada muestra que sería la variable a comparar en nuestro experimento.
Si ves el objeto `design` Verás que la matriz tiene dos columnas: una que especifica el intercepto y otra que especifica el cambio de las condiciones de cultivo.
>El intercepto es una columna que se le añade un valor arbitrario de 1.
>:nerd_face: En estadística, cuando tenemos una variable de factor con *k* niveles, necesitamos convertirla en *k−1* variables indicadoras. Elegimos un nivel como línea de base y luego tenemos una variable indicadora para cada uno de los niveles restantes.
La función `lmFit()` del paquete `limma` ajusta un modelo lineal para cada fila de nuestra matriz de expresión. En el caso de una comparación de dos clases, esto es equivalente a una prueba t simple.
```r=20
fit <- lmFit(cf.norm,design)
```
### Corrección empírica de Bayes en limma
Empirical Bayes (eBayes) es un método que toma prestada información sobre la distribución entre genes para calcular una estadística de prueba robusta. En `limma`, esto se puede realizar usando la función `eBayes()`. La función requiere que proporcionemos un objeto devuelto al ajustar un modelo lineal (o matriz de contraste) a los datos.
```r=21
fitted.ebayes <- eBayes(fit)
```
### Extracción de genes expresados diferencialmente
Ahora tenemos ajustes de modelo para cada característica en la matriz, y podemos organizarlos en una tabla usando `topTable()`. De forma predeterminada, topTable proporcionará las 10 características principales ordenadas por la estadística "B", que es el log de probabilidades de expresión diferencial.
```r=22
topTable(fitted.ebayes)
```
Otra buena función es `decideTests()`. Esta función devuelve qué características probadas de una matriz pasan los criterios de prueba. De forma predeterminada, este es un valor *p* (Benjamini-Hochberg) de 0.05, y sin punto de corte de log fold change. En el ejemplo vamos a poner un punto de corte de 1 en el log fold change (lo que significa un cambio de 2 veces (2-fold) cualquier dirección).
```r=23
summary(decideTests(fitted.ebayes[,"SCGM"],lfc=1))
```
### Usando contrastes
El intercepto es básicamente inútil al hacer este tipo de pruebas, y lo que realmente queremos es comparar grupos experimentales, ¿no? Entonces, ¿podemos crear un diseño que refleje los cambios en los medios del grupo? No tenemos que hacer esto solo para dos condiciones, pero se vuelve esencial para diseños experimentales más complejos.
Para lograr esto tenemos que especificar contrastes dentro de nuestro diseño experimental. Los contrastes son las comparaciones que deseamos hacer. Si `treat` representa dos tratamientos alternativos, simplemente queremos saber si la expresión génica es diferente para un tratamiento frente al otro, calculamos las medias de grupo para cada tratamiento y el contraste (la diferencia) entre ellos.
```r=24
design <- model.matrix( ~ 0 + cf.norm[['culture medium:ch1']])
colnames(design) <- levels(as.factor(cf.norm[['culture medium:ch1']]))
```
>Al añadir `~ 0 +` al diseño del modelo estamos eliminando el intercepto
Puedes volver a ver la matriz `design` y ya no aparece el intercepto. Nuestros coeficientes ahora corresponden a las condiciones experimentales. Esto es más natural, pero crea un paso extra: decidir los contrastes entre grupos (Las comparaciones a realizar por el análisis estadístico). En este caso es fácil porque solo hay dos grupos. En otros casos, tendremos más opciones.
```r=26
contrast_matrix <- makeContrasts(SCGM - MEGM, levels=design)
contrast_matrix
```
En nuestra matriz de contraste, estamos interesados en encontrar la diferencia entre el grupo cultivado en SCGM (el fenotipo EMT) y el control (cultivado en MEGM). Por ello utilizamos SCGM-MEGM, siendo este último el de referencia. **Esta formulación significa que cuando el log-fold change es positivo, la expresión es mayor en SCGM que en MEGM, y viceversa.** Naturalmente, puede realizar análisis más complejos con múltiples comparaciones según la pregunta de interés. Puedes leer más sobre el uso de limma para análisis más complejos en la guía del usuario de limma [aqui](https://bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf) y, en particular, páginas 35-64), que demuestra el uso de limma para una amplia gama de preguntas e incluso plataformas de dos colores.
Ahora podemos pasar al ajuste del modelo otra vez :slightly_smiling_face:
```r=28
fit <- lmFit(cf.norm,design)
fit2 <- contrasts.fit(fit,contrasts=contrast_matrix)
fit2 <- eBayes(fit2)
summary(decideTests(fit2,lfc=1))
```
# Anotando los datos
Podemos añadir los nombres (símbolos) de los genes que se corresponden a las sondas usadas en el microarray. En este caso vamos usar las anotacione scorrespondientes al microarray "hgu133plus2.db" que es el usado en este experimento.
```r=32
#Primero seleccionamos los genes con una p significativa en nuetsro modelo
tab <- topTable(fit2, number = 3e6, adjust.method = "BH", sort.by = "p")
pr <- row.names(tab)
#Obtenemos los nombres correspondientes de los genes usados en el mciroarray
sy <- getSYMBOL(pr, data = "hgu133plus2.db")
#Creamos un dataframe juntando los dos datos
df <- data.frame(Symbol = sy, tab)
View(df)
```
### Visualizando datos
#### Heatmaps (Mapas de calor)
Vamos a cosntruir un heatmap seleccionando las 50 sondas más importantes.
```r=40
df.subset <- df[1:50,]
ex <- exprs(cf.norm)[row.names(df.subset),]
sy.subset <- df.subset$Symbol
```
Graficamos el heatmap
```r=43
heatmap.2(ex, trace = "none", scale="row", margins = c(10, 7), col = colorRampPalette(c("red", "black", "green")), labRow = sy.subset)
```
### Volcano plot
```r=44
#Transformar los datos al formato necesario del gráfico
dat <- data.frame(df, n_log10_adj_pval = -c(log10(df$adj.P.Val)), col=ifelse(abs(df$logFC)>2 & df$adj.P.Val<0.01, 'A', ifelse(abs(df$logFC)>2, 'B', 'C')))
dat$Symbol<-c(as.character(dat$Symbol[1]), rep(NA, 6), as.character(dat$Symbol[8]), rep(NA, 36), as.character(dat$Symbol[45]), rep(NA, length(dat$Symbol)-45))
#Construcción del gráfico
a <- ggplot(dat, aes(x = logFC, y = n_log10_adj_pval, col=col, label=Symbol))
a <- a+ylab("-log10(adjusted P value)\n")
a <- a+xlab("logFC")
a <- a+theme_classic(base_size = 12)
a <- a+theme(legend.position="none")
a <- a+geom_point()
a <- a+geom_vline(xintercept=c(2, -2))
a <- a+geom_hline(yintercept = -log10(0.01))
a <- a+geom_label()
plot(a)
```
---
<font color = 'gray'> :memo: Ejercicios basados en The GTK-teaching/lesson-theme, derived from The Carpentries & Rpubs from lesson: Microarray data, limma, and GEO Anton Wellstein, Garrett Grahamstyle </font>
Sign in with Wallet
Connect another wallet