Configuración del entorno y paquetes necesarios

La secuenciación de ARN (RNA-seq) es una de las técnicas principales para cuantificar la expresión génica. Este método permite identificar diferencias entre condiciones experimentales, descubrir nuevas variantes y caracterizar mutaciones. El flujo de trabajo descrito aquí abarca la preparación de los datos y la ejecución del análisis estadístico.

Instalación de bibliotecas en R

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("DESeq2")
library(DESeq2)

BiocManager::install("ggplot2")
library(ggplot2)

BiocManager::install("clusterProfiler")
library(clusterProfiler)

BiocManager::install("biomaRt")
library(biomaRt)

BiocManager::install("ReactomePA")
library(ReactomePA)

BiocManager::install("org.Hs.eg.db")
library(org.Hs.eg.db)

install.packages("pheatmap")
library(pheatmap)

install.packages("RColorBrewer")
library(RColorBrewer)

install.packages("dplyr")
library(dplyr)

Preparación de los datos de recuento

Se utiliza una tabla de recuentos generada por featureCounts como antrada para el análisis con DESeq2. El archivo debe contener los identificadores de genes como nombres de fila.

# Cargar la matriz de recuentos
matriz_conteos <- read.table("ruta/archivo_conteos.txt", header = TRUE, 
                             skip = 1, row.names = 1)

# Limpiar nombres de columna
colnames(matriz_conteos) <- gsub(".bam", "", colnames(matriz_conteos), fixed = TRUE)
colnames(matriz_conteos) <- gsub("\\.\\.", "", colnames(matriz_conteos))

# Eliminar columnas no necesarias
matriz_conteos <- matriz_conteos[, -c(1:5)]

# Vista preliminar
head(matriz_conteos)

Integración de metadatos experimentales

Los metadatos deben incluir la información de las muestras, con el identificador de muestra como primera columna.

# Cargar archivo de metadatos
expmat <- read.delim("ruta/metadatos.txt", row.names = 1)
expmat$identificador <- rownames(expmat)

# Ordenar para coincidir con la matriz de conteos
expmat <- expmat[match(colnames(matriz_conteos), expmat$identificador), ]

head(expmat)

Construcción del objeto DESeqDataSet

# Crear el objeto DESeq2
dds_experimento <- DESeqDataSetFromMatrix(
  countData = matriz_conteos,
  colData = expmat,
  design = ~Grupo
)

# Ejecutar el análisis diferencial
dds_experimento <- DESeq(dds_experimento)

Extracción de resultados estadísticos

# Obtener resultados con corrección por FDR
res_analisis <- results(dds_experimento, pAdjustMethod = "fdr", alpha = 0.05)
summary(res_analisis)

# Inspeccionar metadatos de los resultados
mcols(res_analisis, use.names = TRUE)

Anotación de identificadores génicos

Para enriquecer el análisis, los símbolos de genes se asocian con nombres completos y códigos de acceso mediante bases de datos biológicas.

# Utilizar base de datos de anotación humana
library(org.Hs.eg.db)

# Agregar anotaciones al dataframe de resultados
res_analisis$nombre_completo <- mapIds(
  x = org.Hs.eg.db,
  keys = rownames(res_analisis),
  column = "GENENAME",
  keytype = "SYMBOL",
  multiVals = "first"
)

res_analisis$simbolo <- rownames(res_analisis)

res_analisis$entrez_id <- mapIds(
  x = org.Hs.eg.db,
  keys = rownames(res_analisis),
  column = "ENTREZID",
  keytype = "SYMBOL",
  multiVals = "first"
)

res_analisis$ensembl_id <- mapIds(
  x = org.Hs.eg.db,
  keys = rownames(res_analisis),
  column = "ENSEMBL",
  keytype = "SYMBOL",
  multiVals = "first"
)

# Filtrar genes significativos (valor p ajustado < 0.05)
genes_significativos <- subset(res_analisis, padj < 0.05)

# Exportar resultados
write.table(as.data.frame(counts(dds_experimento, normalized = TRUE)), 
            file = "conteos_normalizados.txt", sep = "\t", quote = FALSE)

write.table(as.data.frame(genes_significativos), 
            file = "resultados_anotados_significativos.txt", sep = "\t", quote = FALSE)

Visualización de datos

Análisis de componentes principales (PCA)

# Transformación logarítmica regularizada
rld_experimento <- rlog(dds_experimento, blind = FALSE)

# Generar PCA
plotPCA(rld_experimento, intgroup = "Grupo", ntop = 500) +
  geom_point(size = 4) +
  ggtitle("PCA - 500 genes más variables") +
  theme_minimal()

Mapa de calor de expresión

# Seleccionar genes significativos
matriz_expresion <- assay(rld_experimento[rownames(genes_significativos)[1:40], ])

# Definir anotaciones para columnas
info_columnas <- data.frame(
  Grupo = factor(colData(rld_experimento)$Grupo),
  fila = rownames(colData(rld_experimento))
)

# Configurar colores
colores_anotacion <- list(
  Grupo = c(Control = "#4DAF4A", Tratamiento = "#984EA3")
)

# Crear mapa de calor
pheatmap(
  matriz_expresion,
  color = colorRampPalette(brewer.pal(9, "Blues"))(255),
  scale = "row",
  annotation_col = info_columnas,
  annotation_colors = colores_anotacion,
  show_colnames = FALSE,
  clustering_distance_rows = "correlation"
)

Gráfico de volcan

# Preparar datos para el volcan plot
datos_volcan <- data.frame(
  gen = rownames(res_analisis),
  log_pval = -log10(res_analisis$padj),
  cambio_log2 = res_analisis$log2FoldChange
)

datos_volcan <- na.omit(datos_volcan)

# Clasificar por dirección y significancia
datos_volcan <- mutate(datos_volcan, 
  categoria = case_when(
    cambio_log2 > 0 & log_pval > 1.3 ~ "Sobreexpresado",
    cambio_log2 < 0 & log_pval > 1.3 ~ "Subexpresado",
    log_pval <= 1.3 ~ "No significativo"
  )
)

# Crear gráfico
ggplot(datos_volcan, aes(x = cambio_log2, y = log_pval, color = categoria)) +
  geom_point(alpha = 0.7, size = 1.5) +
  scale_color_manual(values = c("#E41A1C", "#377EB8", "#999999")) +
  geom_hline(yintercept = 1.3, linetype = "dashed") +
  labs(x = "Cambio en log2 (Tratamiento/Control)", 
       y = "-log10(p-valor ajustado)") +
  theme_bw()

Otras visualizaciones útiles

# Gráfico MA
plotMA(res_analisis, ylim = c(-5, 5))

# Estimación de dispersiones
plotDispEsts(dds_experimento)

# Expresión de gen individual
gen_objetivo <- rownames(res_analisis)[which.min(res_analisis$padj)]
plotCounts(dds_experimento, gene = gen_objetivo, intgroup = "Grupo")

Análisis de enriquecimiento de rutas

Preparación de datos para el análisis

# Filtrar genes con identificador Entrez
genes_entrez <- subset(genes_significativos, !is.na(entrez_id))

# Crear vector de cambios de expresión
vector_cambio <- genes_entrez$log2FoldChange
names(vector_cambio) <- genes_entrez$entrez_id

Enriquecimiento en rutas KEGG

# Analizar enriquecimiento en KEGG
enriquecimiento_kegg <- enrichKEGG(
  gene = names(vector_cambio),
  organism = 'hsa',
  pvalueCutoff = 0.05
)

# Visualizar resultados principales
barplot(enriquecimiento_kegg, showCategory = 10, font.size = 8)

Enriquecimiento en Gene Ontology

# Analizar enriquecimiento en GO (procesos biológicos)
enriquecimiento_go <- enrichGO(
  gene = names(vector_cambio),
  OrgDb = org.Hs.eg.db,
  ont = "BP",
  pAdjustMethod = "BH",
  pvalueCutoff = 0.01
)

# Mostrar categorías principales
dotplot(enriquecimiento_go, showCategory = 15)

Visualización de rutas metabólicas

# Cargar paquete de visualización de rutas
library(pathview)

# Visualizar ruta específica con cambios de expresión
pathview(
  gene.data = vector_cambio,
  pathway.id = "hsa04110",
  species = "hsa",
  limit = list(gene = 2, cpd = 1)
)