Análisis de Enriquecimiento Funcional y mapas de vías de DE-IGs en Telencéfalo de Ratón
Dra. Katia Aviña Padilla
11-Nov-2024
Introducción
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.
# Cargar los paquetes necesarios
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("clusterProfiler", "org.Mm.eg.db", "enrichplot", "fgsea", "ggplot2"))## Bioconductor version 3.19 (BiocManager 1.30.25), R 4.4.1 (2024-06-14)## Warning: package(s) not installed when version(s) same as or greater than current; use
## `force = TRUE` to re-install: 'clusterProfiler' 'org.Mm.eg.db' 'enrichplot'
## 'fgsea' 'ggplot2'## Old packages: 'bigD', 'bio3d', 'curl', 'fs', 'future.apply', 'ggfun', 'gtable',
## 'Hmisc', 'httr2', 'igraph', 'knitr', 'Matrix', 'mvtnorm', 'parallelly',
## 'pkgbuild', 'progressr', 'ps', 'quantreg', 'R.oo', 'Rcpp', 'rmarkdown',
## 'rstudioapi', 'segmented', 'sf', 'spatstat.explore', 'spatstat.univar',
## 'spatstat.utils', 'terra', 'tinytex', 'waldo', 'withr', 'xfun', 'yulab.utils'# Cargar las librerías
library(clusterProfiler) #análisis de enriquecimiento funcional y el análisis de conjuntos de genes.## ## clusterProfiler v4.12.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':
##
## filterlibrary(org.Mm.eg.db) #base de datos de anotaciones para el organismo 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, 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## 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) # visualización de los resultados de enriquecimiento funcional generados con clusterProfiler.
library(ggplot2) #permite crear gráficos altamente personalizables
library(fgsea) #analizar datos de expresión génica donde se quiere evaluar el enriquecimiento de rutas o conjuntos de genes Paso 1: Cargar los Datos de IGs y los Resultados de Expresión Diferencial
###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"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# 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 6 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, title = "GSEA - GO BP")
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 columnsprint(head(IGs_entrez)) # Verificar los ENTREZ IDs generados## [1] "17202" "243923" "93897" "12570" "18014" "12741"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")
Paso 4: 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")
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")
Paso 5: Análisis de Enriquecimiento ORA
El Análisis de Enriquecimiento ORA (Over-Representation Analysis) es un método estadístico utilizado en la bioinformática para identificar términos biológicos, como procesos biológicos, funciones moleculares o vías metabólicas, que están sobrerrepresentados en un conjunto de genes o proteínas de interés. Este enfoque permite explorar las funciones celulares y procesos en los que está implicado un subconjunto específico de genes, proporcionando información valiosa sobre los mecanismos moleculares subyacentes a una condición biológica o un fenómeno estudiado.
gene_list <- setNames(igs_deg_data$log2FoldChange, igs_deg_data$ID)
gene_list <- sort(gene_list, decreasing = TRUE) # Orden descendente
head(gene_list)## ENSMUSG00000056043 ENSMUSG00000081683 ENSMUSG00000051726 ENSMUSG00000047259
## 1.257809 -1.139067 -1.651324 -1.742849
## ENSMUSG00000090071 ENSMUSG00000041378
## -1.966939 -2.556145# Crear un vector de log2FoldChange con los nombres de los genes (ENTREZ ID o ENSEMBL ID)
gene_list <- setNames(igs_deg_data$log2FoldChange, igs_deg_data$ID)
gene_list <- sort(gene_list, decreasing = TRUE) # Orden descendente para GSEA
head(gene_list)## ENSMUSG00000056043 ENSMUSG00000081683 ENSMUSG00000051726 ENSMUSG00000047259
## 1.257809 -1.139067 -1.651324 -1.742849
## ENSMUSG00000090071 ENSMUSG00000041378
## -1.966939 -2.556145# Análisis ORA para cada ontología utilizando IDs de Ensembl
ora_bp <- enrichGO(
gene = names(gene_list),
OrgDb = org.Mm.eg.db,
keyType = "ENSEMBL", # Especifica que los IDs son de Ensembl
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
readable = TRUE
)
ora_mf <- enrichGO(
gene = names(gene_list),
OrgDb = org.Mm.eg.db,
keyType = "ENSEMBL", # Especifica que los IDs son de Ensembl
ont = "MF",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
readable = TRUE
)
ora_cc <- enrichGO(
gene = names(gene_list),
OrgDb = org.Mm.eg.db,
keyType = "ENSEMBL", # Especifica que los IDs son de Ensembl
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
readable = TRUE
)
head(ora_mf)## ID Description GeneRatio BgRatio RichFactor
## GO:0070888 GO:0070888 E-box binding 2/8 59/22548 0.03389831
## FoldEnrichment zScore pvalue p.adjust qvalue
## GO:0070888 95.54237 13.69872 0.0001865725 0.005783747 0.002945882
## geneID Count
## GO:0070888 Neurog1/Bhlhe22 2head(ora_cc)## [1] ID Description GeneRatio BgRatio RichFactor
## [6] FoldEnrichment zScore pvalue p.adjust qvalue
## [11] geneID Count
## <0 rows> (or 0-length row.names)head(ora_bp)## ID Description GeneRatio BgRatio
## GO:0030900 GO:0030900 forebrain development 3/9 426/23040
## GO:0021987 GO:0021987 cerebral cortex development 2/9 121/23040
## GO:0021543 GO:0021543 pallium development 2/9 183/23040
## GO:0021562 GO:0021562 vestibulocochlear nerve development 1/9 10/23040
## GO:0021603 GO:0021603 cranial nerve formation 1/9 10/23040
## GO:0021957 GO:0021957 corticospinal tract morphogenesis 1/9 10/23040
## RichFactor FoldEnrichment zScore pvalue p.adjust
## GO:0030900 0.007042254 18.02817 7.012627 0.0004852890 0.04914891
## GO:0021987 0.016528926 42.31405 9.007189 0.0009612777 0.04914891
## GO:0021543 0.010928962 27.97814 7.243103 0.0021774297 0.04914891
## GO:0021562 0.100000000 256.00000 15.943728 0.0039001512 0.04914891
## GO:0021603 0.100000000 256.00000 15.943728 0.0039001512 0.04914891
## GO:0021957 0.100000000 256.00000 15.943728 0.0039001512 0.04914891
## qvalue geneID Count
## GO:0030900 0.02433418 Cdk5r2/Neurog1/Bhlhe22 3
## GO:0021987 0.02433418 Cdk5r2/Bhlhe22 2
## GO:0021543 0.02433418 Cdk5r2/Bhlhe22 2
## GO:0021562 0.02433418 Neurog1 1
## GO:0021603 0.02433418 Neurog1 1
## GO:0021957 0.02433418 Bhlhe22 1# Visualización de ORA
dotplot(ora_bp, showCategory = 10, title = "ORA - GO:BP") 
dotplot(ora_mf, showCategory = 10, title = "ORA - GO:MF") 
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(ora_bp, showCategory = 10, circular = TRUE, colorEdge = TRUE)## Warning in cnetplot.enrichResult(x, ...): Use 'color.params = list(edge = your_value)' instead of 'colorEdge'.
## The colorEdge parameter will be removed in the next version.