1 Introducción

El análisis de enriquecimiento funcional es una herramienta esencial en bioinformática, que permite interpretar listas de genes —como aquellos identificados como diferencialmente expresados (DEGs)— para descubrir funciones biológicas, procesos moleculares y rutas celulares relevantes asociadas en un experimento. Este tipo de análisis permite identificar si ciertos procesos biológicos, funciones moleculares o componentes celulares están sobrerrepresentados en el conjunto de genes de interés, en comparación con lo que se esperaría al azar. Mediante herramientas como clusterProfiler, es posible realizar análisis de enriquecimiento usando ontologías como Gene Ontology (GO) o bases de datos de rutas como KEGG, lo que facilita la identificación de mecanismos biológicos relevantes. Además, representaciones gráficas como cnetplot permiten visualizar redes de términos funcionales conectados a genes específicos, proporcionando una perspectiva integrada de la función biológica de los genes estudiados.

En esta práctica, realizaremos un análisis de enriquecimiento funcional de los genes sin intrones (IGs) que están diferencialmente expresados (DEGs) en el telencéfalo de ratón. Usaremos clusterProfiler para analizar el enriquecimiento en términos GO, KEGG, y Cnetplot para construir una red de ontologías.

1.1 Objetivos de Aprendizaje

Al finalizar esta práctica, serás capaz de:

  • Comprender el propósito del análisis de enriquecimiento funcional y su utilidad en estudios transcriptómicos.
  • Aplicar métodos GSEA y ORA para identificar funciones biológicas relevantes en genes diferencialmente expresados.
  • Interpretar visualizaciones de resultados de GO, KEGG y cnetplot.
  • Analizar el papel funcional de los “genes intronless (IGs)” en el contexto del desarrollo del telencéfalo en ratón.

1.2 Terminología clave

Conjuntos génicos y rutas (pathways) Un conjunto génico es una colección no ordenada de genes que están funcionalmente relacionados. Una ruta o pathway puede interpretarse como un conjunto génico si se ignoran las relaciones funcionales específicas entre los genes. Ontología Génica (Gene Ontology, GO) La Ontología Génica (GO) define conceptos o clases que se utilizan para describir la función de los genes, así como las relaciones entre estos conceptos. Esta ontología clasifica las funciones génicas en tres grandes categorías: MF: Función Molecular (Molecular Function) Actividades bioquímicas llevadas a cabo por los productos génicos. CC: Componente Celular (Cellular Component) Localización celular donde los productos génicos ejercen su función. BP: Proceso Biológico (Biological Process) Conjunto de eventos o rutas biológicas que involucran múltiples actividades moleculares. Los términos de GO están organizados en un grafo dirigido acíclico, en el cual las conexiones representan relaciones jerárquicas tipo padre-hijo.

1.3 Glosario de Términos Clave

  • Genes sin intrones (IGs): Genes compuestos por un solo exón, es decir, no contienen intrones que deban ser eliminados por el splicing.
  • DEGs (Differentially Expressed Genes): Genes cuya expresión difiere significativamente entre dos condiciones biológicas.
  • GO (Gene Ontology): Sistema de clasificación que describe funciones biológicas, moleculares y localización celular de genes/proteínas.
  • KEGG (Kyoto Encyclopedia of Genes and Genomes): Base de datos que contiene rutas metabólicas y de señalización.
  • GSEA (Gene Set Enrichment Analysis): Método para determinar si un conjunto de genes muestra diferencias estadísticamente significativas entre dos condiciones.
  • ORA (Over-Representation Analysis): Método que evalúa si ciertos términos funcionales están sobrerrepresentados en una lista de genes.
  • Cnetplot: Visualización de red funcional que relaciona genes con términos GO o KEGG enriquecidos.
# Cargar los paquetes necesarios
#if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
#BiocManager::install(c("clusterProfiler", "org.Mm.eg.db", "enrichplot", "fgsea", "ggplot2", "DOSE"))

#Cargar las librerías
library(clusterProfiler)   # análisis de enriquecimiento funcional
## Warning: package 'clusterProfiler' was built under R version 4.4.2
## 
## clusterProfiler v4.14.6 Learn more at https://yulab-smu.top/contribution-knowledge-mining/
## 
## Please cite:
## 
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan,
## X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal
## enrichment tool for interpreting omics data. The Innovation. 2021,
## 2(3):100141
## 
## Attaching package: 'clusterProfiler'
## The following object is masked from 'package:stats':
## 
##     filter
library(org.Mm.eg.db)      # anotaciones de Mus musculus
## Loading required package: AnnotationDbi
## Loading required package: stats4
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, saveRDS, setdiff,
##     table, tapply, union, unique, unsplit, which.max, which.min
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: IRanges
## Warning: package 'IRanges' was built under R version 4.4.2
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:clusterProfiler':
## 
##     rename
## The following object is masked from 'package:utils':
## 
##     findMatches
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## 
## Attaching package: 'IRanges'
## The following object is masked from 'package:clusterProfiler':
## 
##     slice
## 
## Attaching package: 'AnnotationDbi'
## The following object is masked from 'package:clusterProfiler':
## 
##     select
## 
library(enrichplot)        # visualizaciones de enriquecimiento
## Warning: package 'enrichplot' was built under R version 4.4.2
## enrichplot v1.26.6 Learn more at https://yulab-smu.top/contribution-knowledge-mining/
## 
## Please cite:
## 
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024,
## 5(6):100722
library(DOSE)              # soporte para cnetplot y otros
## Warning: package 'DOSE' was built under R version 4.4.3
## DOSE v4.0.1 Learn more at https://yulab-smu.top/contribution-knowledge-mining/
## 
## Please cite:
## 
## Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an
## R/Bioconductor package for Disease Ontology Semantic and Enrichment
## analysis. Bioinformatics. 2015, 31(4):608-609
library(ggplot2)        # visualización de los resultados de enriquecimiento funcional generados con clusterProfiler.    
library(fgsea)             #analizar datos de expresión génica donde se quiere evaluar el enriquecimiento de rutas o conjuntos de genes 
## Warning: package 'fgsea' was built under R version 4.4.3
library(cowplot)

2 Práctica de Enriquecimiento Funcional

2.1 Paso 1: Cargar los Datos de IGs y los Resultados de Expresión Diferencial

2.1.1 Cargar el archivo de genes sin intrones (IGs)

# Cargar el archivo de genes sin intrones (IGs)
igs_data <- read.table("/Users/katiaavinapadilla/Desktop/mus_musculus-intronless_genes.tsv", header = TRUE, sep = "\t", quote = "", fill = TRUE, comment.char = "")
IGs <- igs_data$id  # Extraer la columna con los IDs de los genes IGs
print(head(IGs))    # Verificar la extracción de IDs
## [1] "ENSMUSG00000094307" "ENSMUSG00000094860" "ENSMUSG00000094596"
## [4] "ENSMUSG00000095814" "ENSMUSG00000044597" "ENSMUSG00000094273"
# Cargar el archivo de resultados de expresión diferencial (DEGs)
deg_data <- read.csv("/Users/katiaavinapadilla/Desktop/DEGs_mouse.csv", header = TRUE, quote = "", fill = TRUE, comment.char = "")
degs <- deg_data$ID  # Nombre de la columna de IDs en `deg_data`
print(head(degs))    # Verificar la extracción de IDs de DEGs
## [1] "ENSMUSG00000024112" "ENSMUSG00000040447" "ENSMUSG00000039860"
## [4] "ENSMUSG00000030310" "ENSMUSG00000034227" "ENSMUSG00000034796"

2.2 Paso 2: Obtener los IGs que También son DEGs

# Obtener los IGs que también son DEGs
igs_degs <- intersect(IGs, degs)
print(head(igs_degs))  # Verificar el resultado de la intersección
## [1] "ENSMUSG00000047259" "ENSMUSG00000056043" "ENSMUSG00000081683"
## [4] "ENSMUSG00000090071" "ENSMUSG00000048904" "ENSMUSG00000041378"
# Filtrar el archivo deg_data para obtener solo los datos de los IGs que son DEGs
igs_deg_data <- deg_data[deg_data$ID %in% igs_degs, ]  # Ajustar ID si el nombre es diferente

# Verificar el resultado del filtrado
head(igs_deg_data)
##                     ID  baseMean log2FoldChange     lfcSE      stat
## 28  ENSMUSG00000048904 164.84163      -2.598370 0.3498218 -7.427697
## 85  ENSMUSG00000056947  48.35968      -3.119002 0.3908627 -7.979788
## 86  ENSMUSG00000025128  36.09606      -4.140594 0.5431203 -7.623714
## 138 ENSMUSG00000056043  83.37197       1.257809 0.2241705  5.610949
## 169 ENSMUSG00000051726  48.51698      -1.651324 0.2628171 -6.283168
## 219 ENSMUSG00000041378  53.46535      -2.556145 0.4564468 -5.600093
##           pvalue         padj
## 28  1.105050e-13 4.296740e-12
## 85  1.465850e-15 6.593221e-14
## 86  2.464789e-14 1.006297e-12
## 138 2.012198e-08 4.171773e-07
## 169 3.317412e-10 8.781629e-09
## 219 2.142364e-08 4.415766e-07

2.3 Análisis de Enriquecimiento de Conjuntos Génicos (GSEA)

Una estrategia común para analizar perfiles de expresión génica es identificar genes diferencialmente expresados que se consideran interesantes. El análisis de enriquecimiento por sobre-representación (ORA, por sus siglas en inglés) se basa en este conjunto de genes diferencialmente expresados. Sin embargo, este enfoque detecta únicamente los genes cuya diferencia es grande, y falla en detectar aquellos casos donde los cambios individuales son pequeños pero se manifiestan de forma coordinada dentro de un conjunto funcionalmente relacionado.

El Análisis de Enriquecimiento de Conjuntos Génicos (GSEA) (Subramanian et al., 2005) aborda directamente esta limitación. A diferencia de ORA, GSEA utiliza todos los genes del experimento. Agrega las estadísticas por gen dentro de un conjunto génico predefinido, lo que permite detectar patrones donde los genes cambian de manera sutil pero coordinada. Esto es importante porque muchas diferencias fenotípicas relevantes se manifiestan justamente mediante cambios pequeños pero consistentes en varios genes relacionados.

En GSEA, los genes se ordenan según su asociación con un fenotipo. Dado un conjunto de genes previamente definido S (por ejemplo, genes que comparten una categoría de Disease Ontology), el objetivo de GSEA es determinar si los miembros de S están distribuidos aleatoriamente a lo largo de la lista ordenada L o si se concentran principalmente en los extremos (inicio o final).

2.3.1 Elementos clave del método GSEA

2.3.1.1 1. Cálculo del Enrichment Score (ES)

El Enrichment Score (ES) representa el grado en que el conjunto S está sobrerrepresentado al inicio o al final de la lista ordenada L. El puntaje se calcula recorriendo la lista L: se incrementa una suma acumulada cada vez que se encuentra un gen perteneciente a S, y se decrementa cuando no. La magnitud del incremento depende del valor estadístico asociado al gen (por ejemplo, la correlación del gen con el fenotipo).
El ES se define como la desviación máxima respecto a cero durante este recorrido, y corresponde a una estadística similar al test de Kolmogorov-Smirnov ponderado.

2.3.1.2 2. Estimación del nivel de significancia del ES

El valor p del ES se calcula mediante una prueba de permutación. Se permutan las etiquetas fenotípicas de los genes en la lista L y se recalcula el ES para los datos permutados, generando así una distribución nula del ES. El valor p observado se estima comparándolo contra esta distribución nula.

2.3.1.3 3. Corrección por pruebas múltiples

Cuando se evalúan múltiples conjuntos génicos, el nivel de significancia estimado se ajusta para tener en cuenta la corrección por pruebas múltiples. Además, se calculan los q-values para controlar la tasa de falsos descubrimientos (FDR).

Referencia: Subramanian, A. et al. (2005). Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proceedings of the National Academy of Sciences, 102(43), 15545–15550.

# Cargar el paquete enrichplot
library(enrichplot)

# Preparar la lista de genes para GSEA
gene_list <- setNames(deg_data$log2FoldChange, deg_data$ID)
gene_list <- sort(gene_list, decreasing = TRUE)  # Orden descendente

# Realizar GSEA para GO - Procesos Biológicos (BP)
gsea_go_bp <- gseGO(geneList = gene_list, 
                    OrgDb = org.Mm.eg.db, 
                    ont = "BP", 
                    keyType = "ENSEMBL", 
                    pvalueCutoff = 0.05, 
                    pAdjustMethod = "BH")
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 7 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)
## leading edge analysis...
## done...
# Calcular la matriz de similitud de términos
gsea_go_bp <- pairwise_termsim(gsea_go_bp)

# Visualización de resultados de GSEA con dotplot
dotplot(gsea_go_bp, showCategory = 15, title = "GSEA - GO BP")

# Generar el emapplot para visualizar los resultados de GSEA
emapplot(gsea_go_bp, showCategory = 15) + ggtitle("GSEA - GO BP")

2.4 Paso 3: Convertir los IDs de IGs DEGs a ENTREZ ID

# Convertir los IGs DEGs a ENTREZ ID usando ENSEMBL o SYMBOL como tipo de ID de origen
IGs_entrez <- bitr(igs_degs, fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Mm.eg.db)$ENTREZID
## 'select()' returned 1:1 mapping between keys and columns
print(head(IGs_entrez))  # Verificar los ENTREZ IDs generados
## [1] "17202"  "243923" "93897"  "12570"  "18014"  "12741"

2.5 Paso 4: Análisis de Enriquecimiento KEGG

El Análisis de Enriquecimiento de Vías KEGG (Kyoto Encyclopedia of Genes and Genomes) es una herramienta fundamental en la bioinformática y biología de sistemas para comprender las funciones y rutas metabólicas en las que participan grupos específicos de genes. Este análisis permite identificar si un conjunto de genes o proteínas, comúnmente aquellos diferencialmente expresados o asociados a una condición biológica particular, está sobrerrepresentado en determinadas rutas metabólicas o de señalización descrita en la base de datos KEGG.

# Análisis de enriquecimiento KEGG

KEGG_IGs<- enrichKEGG(gene         = IGs_entrez,      # Vector de genes a analizar
                         organism     = "mmu",      # Código KEGG para ratón
                         pvalueCutoff = 0.05)       # Umbral de significancia para p-valor ajustado
## Reading KEGG annotation online: "https://rest.kegg.jp/link/mmu/pathway"...
## Reading KEGG annotation online: "https://rest.kegg.jp/list/pathway/mmu"...
# Visualización del enriquecimiento de GO con un gráfico de puntos
dotplot(KEGG_IGs) + ggtitle("KEGG Pathway Enrichment")

2.6 Paso 5: Análisis de Enriquecimiento GO

El análisis enrichGO compara la frecuencia de aparición de cada proceso en la lista de genes de interés con la frecuencia esperada en un conjunto de referencia, utilizando un fondo genómico amplio. Se aplican métodos estadísticos, como la corrección de p-valores mediante el método de Benjamini-Hochberg, para reducir el riesgo de obtener resultados falsos positivos, garantizando que los términos biológicos identificados sean realmente relevantes.

pvalueCutoff = 0.05: valor de corte para el p-valor, el cual es una medida estadística que indica la probabilidad de que el enriquecimiento observado ocurra por azar. Un p-valor ≤ 0.05 sugiere que hay menos de un 5% de probabilidad de que el enriquecimiento observado sea debido al azar, proporcionando una base para inferir que el término GO está significativamente asociado con el conjunto de genes analizado.

qvalueCutoff = 0.05: Este parámetro define el umbral para el q-valor, que es el p-valor ajustado para múltiples pruebas. Al realizar análisis de enriquecimiento, se prueban muchos términos GO al mismo tiempo, lo cual aumenta el riesgo de obtener resultados significativos solo por azar (falsos positivos). El q-valor representa el p-valor ajustado mediante métodos como el de Benjamini-Hochberg, que controlan la tasa de descubrimiento falso (FDR). En este caso, el q-valor también se ha establecido en 0.05, lo que significa que solo se considerarán significativos los términos GO que tengan un q-valor igual o menor a 0.05. Este ajuste es crucial en estudios genómicos para obtener resultados robustos y confiables, minimizando los falsos positivos.

# Análisis de enriquecimiento GO
egoBP <- enrichGO(gene = IGs_entrez, 
                OrgDb = org.Mm.eg.db, 
                ont = "BP", 
                pAdjustMethod = "BH", 
                pvalueCutoff = 0.05, 
                qvalueCutoff = 0.05, 
                readable = TRUE)



# Visualización del enriquecimiento de GO con un gráfico de puntos
dotplot(egoBP, showCategory = 10, title = "GO - Procesos Biológicos para DE-IGs")

2.7 Paso 5: Análisis de Enriquecimiento GO:MF

egoMF <- enrichGO(gene = IGs_entrez, 
                OrgDb = org.Mm.eg.db, 
                ont = "MF", 
                pAdjustMethod = "BH", 
                pvalueCutoff = 0.05, 
                qvalueCutoff = 0.05, 
                readable = TRUE)
# Visualización del enriquecimiento de GO con un gráfico de puntos
dotplot(egoMF, showCategory = 10, title = "GO -MF para DE-IGs")

2.8 Comparación entre Ontologías

Ontología Tipos de términos predominantes Implicación funcional
BP Desarrollo neuronal, proliferación celular Procesos complejos de desarrollo
MF E-box binding Regulación de transcripción

2.9 Paso 6: Análisis de Enriquecimiento Cnetplot

¿Cómo funciona Cnetplot?

Datos de Enriquecimiento: Primero, se realiza un análisis de enriquecimiento usando funciones como enrichGO o enrichKEGG en el paquete ClusterProfiler. Esto genera una lista de términos enriquecidos y los genes asociados a cada término. Generación de Cnetplot: A partir de los resultados del enriquecimiento, el cnetplot crea una red donde: Los nodos representan tanto genes como términos biológicos enriquecidos. Las conexiones (o edges) unen genes con los términos en los que están involucrados. Un gen puede estar conectado a varios términos si participa en múltiples procesos.

cnetplot(egoBP, showCategory = 10, circular = TRUE, colorEdge = TRUE)

gene_fc <- setNames(deg_data$log2FoldChange, deg_data$ID)
gene_fc <- gene_fc[names(gene_fc) %in% gsea_go_bp@geneList]  # asegurar consistencia

# Red funcional (circular) con color
cnetplot(gsea_go_bp,
         foldChange = gene_fc,
         showCategory = 10,
         circular = TRUE,
         colorEdge = TRUE,
         node_label = "all")

#1. Asegúrate de que tu objeto GSEA esté cargado
#gsea_go_bp

#Convierte ENSEMBL a símbolo
gsea_go_bp_sym <- setReadable(gsea_go_bp, OrgDb = org.Mm.eg.db, keyType = "ENSEMBL")

#Crea el vector foldChange con nombres como los usados internamente
gene_fc <- setNames(deg_data$log2FoldChange, deg_data$ID)
gene_fc <- gene_fc[names(gene_fc) %in% names(gsea_go_bp@geneList)]
p1 <- cnetplot(gsea_go_bp_sym, foldChange = gene_fc)

#Combina en una figura con etiquetas A, B, C
cowplot::plot_grid(p1, labels = LETTERS[1:3])

#Crear las visualizaciones
p1 <- cnetplot(gsea_go_bp_sym, foldChange = gene_fc, node_label = "category", cex_label_category = 1.2)
p2 <- cnetplot(gsea_go_bp_sym, foldChange = gene_fc, node_label = "gene", cex_label_gene = 0.8)
p3 <- cnetplot(gsea_go_bp_sym, foldChange = gene_fc, node_label = "all")
p4 <- cnetplot(gsea_go_bp_sym, foldChange = gene_fc, node_label = "none", 
               color_category = "firebrick", color_gene = "steelblue")

#Unir en panel 2x2
cowplot::plot_grid(p1, p2, p3, p4, ncol = 2, labels = LETTERS[1:4])

3 Preguntas de Reflexión

  1. ¿Qué ventaja podrían tener los genes sin intrones en procesos de desarrollo rápido como la formación del telencéfalo?
  2. ¿Qué diferencias observas entre los resultados de GO:BP y KEGG? ¿Alguno te parece más informativo?
  3. ¿Cómo interpretarías que un solo gen esté conectado con múltiples términos GO en el cnetplot?

4 Conclusión y Discusión

Este análisis permitió identificar rutas y procesos sobrerrepresentados entre los genes sin intrones diferencialmente expresados en el telencéfalo en desarrollo. La presencia de términos relacionados con neurogénesis, regulación transcripcional y procesos sinápticos sugiere que los IGs pueden desempeñar un papel relevante en el establecimiento de redes neuronales durante etapas críticas del desarrollo embrionario.

Discusión sugerida:
¿Podría la ausencia de intrones facilitar una respuesta más rápida o directa en ciertos contextos celulares como el desarrollo temprano? ¿Qué tipo de genes o familias encontraste más representadas entre los IGs?

⚠️ Errores comunes - No verificar que los IDs estén en el formato correcto (ENSEMBL vs ENTREZ) antes de usarlos en enrichGO o enrichKEGG. - Olvidar ordenar el geneList de forma descendente para GSEA. - No usar readable = TRUE, lo que dificulta la interpretación de los genes en los resultados.