La minería de datos genómicos es una herramienta poderosa para descubrir biomarcadores y comprender procesos biológicos. Este artículo detalla un flujo de trabajo para el análisis de expresión génica diferencial en conjuntos de datos GEO que contienen múltiples grupos experimentales, utilizando el entorno R.
Configuración del Entorno y Descarga de Datos
Para asegurar una descarga eficiente de los conjuntos de datos, es recomendable ajustar las opciones de tiempo de espera en R.
# Establece un tiempo de espera prolongado para descargas de archivos grandes
options(timeout = 100000)
A continuación, se cargan las librerías necesarias, se descarga el conjunto de datos GEO, y se realiza un preprocesamiento básico.
library(stringr)
library(tinyarray) # Para funciones de descarga y simplificación
library(ggplot2) # Para visualización
# Identificador del estudio GEO
id_geo = "GSE474"
# Descarga los datos de GEO
datos_geo = geo_download(id_geo)
# Extrae la matriz de expresión
datos_expresion = datos_geo$exp
# Verifica el rango de los valores para determinar si se necesita transformación logarítmica
print(range(datos_expresion))
# Aplica transformación log2(x+1) para normalizar los datos si no están en escala logarítmica
# y evitar valores infinitos con log2(0)
datos_expresion = log2(datos_expresion + 1)
# Genera un boxplot para inspeccionar la distribución de la expresión y detectar anomalías
boxplot(datos_expresion, las = 2, main = "Distribución de la Expresión Génica")
Definición de Grupos Experimentales
Es crucial definir correctamente los grupos de muestras para las comparaciones. En este caso, la información del grupo se extrae de los títulos de las muestras y se convierte en un factor, estableciendo el grupo de control como el primer nivel.
# Extrae la información del grupo a partir de la columna 'title' de los datos fenotípicos (pd)
# Asume que el nombre del grupo es la tercera parte después de dividir por '-'
vector_grupos = str_split(datos_geo$pd$title, "-", simplify = T)[,3]
# Convierte el vector en un factor, definiendo el orden de los niveles
# El primer nivel (NonObese) se considera el grupo de referencia para las comparaciones
vector_grupos = factor(vector_grupos, levels = c("NonObese","Obese","MObese"))
print(vector_grupos)
Anotación de Sondas (Probes)
Las matrices de expresión de microarrays suelen usar IDs de sonda. Para interpretar los resultados, estas sondas deben ser mapeadas a IDs de genes (símbolos o ENTREZID).
# Muestra el identificador de la plataforma (GPL) para encontrar el paquete de anotación adecuado
print(datos_geo$gpl)
# Generalmente, se utiliza la función find_anno para obtener sugerencias de paquetes,
# pero en este ejemplo se asume el paquete 'hgu133a.db'.
# Por ejemplo: find_anno(datos_geo$gpl) podría sugerir 'library(hgu133a.db)'
# Carga la librería de anotación específica para la plataforma (ej. Affymetrix Human Genome U133A Array)
library(hgu133a.db)
# Extrae la tabla de mapeo de sondas a símbolos de genes
tabla_anotacion = toTable(hgu133aSYMBOL)
head(tabla_anotacion)
Análisis de Expresión Diferencial con limma
El paquete limma es ampliamente utilizado para el análisis de expresión diferencial en datos de microarrays y RNA-seq, especialmente en diseños complejos como los de múltiples grupos.
library(limma)
# Diseño experimental: se crea una matriz de diseño para 'limma'
# La fórmula '~0+vector_grupos' genera columnas para cada nivel del factor
matriz_diseno = model.matrix(~0 + vector_grupos)
colnames(matriz_diseno) = levels(vector_grupos) # Asigna nombres de columna basados en los niveles de grupo
# Define los nombres de los grupos para construir los contrastes
nombres_grupos = levels(vector_grupos)
print(nombres_grupos)
# Construye la matriz de contrastes para realizar comparaciones pareadas
# Por ejemplo, Obese vs NonObese, MObese vs Obese, MObese vs NonObese
matriz_contraste <- makeContrasts(paste0(nombres_grupos[2],"-",nombres_grupos[1]),
paste0(nombres_grupos[3],"-",nombres_grupos[2]),
paste0(nombres_grupos[3],"-",nombres_grupos[1]),
levels = matriz_diseno)
# Ajusta el modelo lineal a los datos de expresión
ajuste_limma <- lmFit(datos_expresion, matriz_diseno)
# Aplica los contrastes definidos al modelo ajustado
ajuste_contraste <- contrasts.fit(ajuste_limma, matriz_contraste)
# Aplica eBayes para la contracción de varianzas, mejorando la estimación de la significancia
ajuste_final <- eBayes(ajuste_contraste)
# Extrae los resultados de expresión diferencial para cada contraste
# 'lapply' itera sobre las columnas de la matriz de diseño (que representan los contrastes)
# 'topTable' extrae los genes más diferencialmente expresados
resultados_deg = lapply(1:ncol(matriz_contraste), function(i){
topTable(ajuste_final, coef = i, number = Inf, sort.by = "P")
})
# Asigna los nombres de los contrastes a la lista de resultados
names(resultados_deg) = colnames(matriz_contraste)
print(names(resultados_deg))
# Muestra las primeras filas de los resultados para la primera comparación (Obese vs NonObese)
head(resultados_deg[[1]])
Flujo de Trabajo Simplificado con tinyarray
El paquete tinyarray puede simplificar el proceso de aálisis de expresión diferencial, realizando automáticamente comparaciones pareadas en conjuntos de datos con múltiples grupos.
# La función get_deg_all realiza un análisis diferencial completo para múltiples grupos.
# symmetry = T indica que se realicen todas las comparaciones pareadas posibles.
# logFC_cutoff = 0.585 es un umbral para el logFC (equivalente a un cambio de 1.5 veces).
# cluster_cols = F desactiva el clustering de columnas en los heatmaps resultantes.
# entriz = F indica que los IDs no son ENTREZID inicialmente (se usarán los IDs de sondas y luego la anotación).
resultados_simplificados = get_deg_all(datos_expresion, vector_grupos, tabla_anotacion,
symmetry = T, logFC_cutoff = 0.585,
cluster_cols = F, entriz = F)
# La función devuelve un objeto que contiene los resultados del análisis y las visualizaciones.
# Guardar la imagen del heatmap generado por tinyarray
ggplot2::ggsave("heatmap_deg_tinyarray.png", plot = resultados_simplificados$plots, width = 15, height = 10)
Análisis de Enriquecimiento Funcional
El análisis de enriquecimiento permite identificar funciones biológicas, vías o procesos celulares que están sobrerrepresentados en una lista de genes diferencialmente expresados.
# Primero, se obtiene una lista de genes diferencialmente expresados utilizando get_deg
# (similar a get_deg_all pero puede ser más flexible en la salida para enriquecimiento)
deg_para_enriquecimiento = get_deg(datos_expresion, vector_grupos, tabla_anotacion, logFC_cutoff = 0.585)
# Extraer los ENTREZID de los genes diferencialmente expresados para cada comparación
# Se filtran los elementos "stable" (genes no diferencialmente expresados o sin anotación)
genes_comparacion1 = deg_para_enriquecimiento$'Obese-NonObese'$ENTREZID[deg_para_enriquecimiento$'Obese-NonObese'$ENTREZID != "stable"]
genes_comparacion2 = deg_para_enriquecimiento$'MObese-Obese'$ENTREZID[deg_para_enriquecimiento$'MObese-Obese'$ENTREZID != "stable"]
genes_comparacion3 = deg_para_enriquecimiento$'MObese-NonObese'$ENTREZID[deg_para_enriquecimiento$'MObese-NonObese'$ENTREZID != "stable"]
# Decide si trabajar con la intersección o la unión de las listas de genes
# Un bloque 'if(F)' se usa aquí para alternar fácilmente entre las opciones.
# Para este ejemplo, se utiliza la unión de genes.
if(FALSE){ # Cambiar a TRUE para usar la intersección
genes_para_enriquecimiento = intersect_all(genes_comparacion1, genes_comparacion2, genes_comparacion3) # Intersección
} else {
genes_para_enriquecimiento = union_all(genes_comparacion1, genes_comparacion2, genes_comparacion3) # Unión
}
head(genes_para_enriquecimiento)
# Realiza el análisis de enriquecimiento rápido utilizando quick_enrich
# 'destdir = tempdir()' guarda los resultados temporales en un directorio temporal
# Es posible que este paso falle si hay problemas de conexión a internet para acceder a bases de datos de enriquecimiento.
resultados_enriquecimiento = quick_enrich(genes_para_enriquecimiento, destdir = tempdir())
# Muestra los nombres de los elementos en la lista de resultados de enriquecimiento
print(names(resultados_enriquecimiento))
# Muestra las primeras filas de la tabla de resultados del primer tipo de enriquecimiento
resultados_enriquecimiento[[1]][1:4,1:4]
# Combina gráficos de enriquecimiento si la salida de quick_enrich son objetos ggplot
# Se utiliza el paquete 'patchwork' para organizar múltiples gráficos
library(patchwork)
combinacion_plots = resultados_enriquecimiento[[3]] + resultados_enriquecimiento[[4]]
print(combinacion_plots)
# Guarda la combinación de gráficos de enriquecimiento
ggplot2::ggsave("enriquecimiento_combinado.png", plot = combinacion_plots, width = 12, height = 7)