Enriquecimiento con ClusterProfiler

ClusterProfiler es un paquete de R que permite comparar tus datos con bases de datos externas y actualizadas como KEGG, GO, entre otros. Para más información, consulte: https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html

0. Instalación del paquete

En caso ya lo tenga instalado, puede obviar este paso

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

BiocManager::install("clusterProfiler")

1. Datos que serán usados

Usaremos datos de expresión reales entre dos grupos de tejidos.

Para más información sobre el código Entrez de los genes, accese: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3013746/

#Crearemos una lista con los nombres de los genes en Código Entrez.
genes <- c(4312, 3552, 4321, 2810, 4171, 7422, 4292, 4193, 387,6774, 6789,  595, 26524,7157, 5728, 5037,11035, 3205, 8829, 7994,60485, 4297, 6934,79923,596, 5308, 7099, 5243, 8434, 6935, 9314, 1958, 2247, 2353, 7291,10891,216)

#Ahora una para los valores de expresión (FOLD CHANGE).
fold_change <-c(28.000000,13.390000,9.970000,4.260000,2.140000,2.060000,-1.298701,-1.408451,-1.515152,-1.515152,-1.538462,-1.639344,-1.724138,-1.754386,-1.818182,-1.886792,-1.886792,-1.960784,-2.083333,-2.127660,-2.222222,-2.500000,-2.777778,-3.030303,-3.225806,-3.333333,-3.571429,-4.166667,-4.166667,-4.347826,-4.545455,-5.263158,-5.555556,-5.555556,-5.555556,-9.090909,-11.111111)

#Asignar los nombres de genes para cada resultado de expresión
names(fold_change)<-genes
#Resultado final
fold_change
##       4312       3552       4321       2810       4171       7422       4292 
##  28.000000  13.390000   9.970000   4.260000   2.140000   2.060000  -1.298701 
##       4193        387       6774       6789        595      26524       7157 
##  -1.408451  -1.515152  -1.515152  -1.538462  -1.639344  -1.724138  -1.754386 
##       5728       5037      11035       3205       8829       7994      60485 
##  -1.818182  -1.886792  -1.886792  -1.960784  -2.083333  -2.127660  -2.222222 
##       4297       6934      79923        596       5308       7099       5243 
##  -2.500000  -2.777778  -3.030303  -3.225806  -3.333333  -3.571429  -4.166667 
##       8434       6935       9314       1958       2247       2353       7291 
##  -4.166667  -4.347826  -4.545455  -5.263158  -5.555556  -5.555556  -5.555556 
##      10891        216 
##  -9.090909 -11.111111

2. Abrir el paquete y generar el análisis de enriquecimiento.

Para este caso, usaremos datos de la enciclopedia de Kyoto como colección base. Podemos ver más información de esta base de datos en el siguiente link: https://www.genome.jp/kegg/

Sin embargo, también podemos usar este paquete para comparar con datos de Gene Ontology y otros.

library(clusterProfiler)
## 
## clusterProfiler v3.16.1  For help: https://guangchuangyu.github.io/software/clusterProfiler
## 
## If you use clusterProfiler in published research, please cite:
## Guangchuang Yu, Li-Gen Wang, Yanyan Han, Qing-Yu He. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology. 2012, 16(5):284-287.
## 
## Attaching package: 'clusterProfiler'
## The following object is masked from 'package:stats':
## 
##     filter
#Generar un enriquecimiento con datos a partir de información del KEGG

KEGG_genes <- enrichKEGG(gene         = genes, organism     = "hsa", pvalueCutoff = 0.05)
## Reading KEGG annotation online:
## Reading KEGG annotation online:
head(KEGG_genes)
##                ID             Description GeneRatio  BgRatio       pvalue
## hsa05205 hsa05205 Proteoglycans in cancer     10/35 205/8093 9.232181e-09
## hsa05206 hsa05206     MicroRNAs in cancer     11/35 310/8093 4.000791e-08
## hsa05210 hsa05210       Colorectal cancer      7/35  86/8093 6.314302e-08
## hsa05215 hsa05215         Prostate cancer      7/35  97/8093 1.459816e-07
## hsa04115 hsa04115   p53 signaling pathway      6/35  73/8093 5.770007e-07
## hsa05219 hsa05219          Bladder cancer      5/35  41/8093 7.531107e-07
##              p.adjust       qvalue
## hsa05205 1.514078e-06 8.454734e-07
## hsa05206 3.280648e-06 1.831941e-06
## hsa05210 3.451819e-06 1.927524e-06
## hsa05215 5.985245e-06 3.342210e-06
## hsa04115 1.892562e-05 1.056822e-05
## hsa05219 2.058503e-05 1.149485e-05
##                                                       geneID Count
## hsa05205    7422/4193/387/6774/595/7157/79923/7099/2247/7291    10
## hsa05206 7422/4193/387/6774/595/7157/5728/596/5243/8434/6935    11
## hsa05210                     4292/387/595/7157/6934/596/2353     7
## hsa05215                    4193/595/7157/5728/6934/596/6935     7
## hsa04115                         2810/4193/595/7157/5728/596     6
## hsa05219                             4312/7422/4193/595/7157     5
#Ahora intentaremos el gráfico de puntos de estas asociaciones

dotplot(KEGG_genes)

#Podemos aumentar (o disminuir) el número de vias metabólicas a ser mostradas

#En este caso, escogeremos 20 lineas

dotplot(KEGG_genes, showCategory=20)

3. Mapas de vías

Una vez que las vías son determinadas (enriquecidas), podemos evaluar si existe alguna interacción entre ellas, para esto usaremos la función emapplot del paquete enrichplot:

#Instalar Enrichplot (en caso no lo tenga)
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("enrichplot")
#Ejecutar la función

enrichplot::emapplot(KEGG_genes)

#En caso hayan muchos elementos, también podemos seleccionar los más relevantes

enrichplot::emapplot(KEGG_genes, showCategory = 10)

4. GSEA (Gene Set Enrichment Analysis)

Para este análisis, también consideraremos los niveles de expresión de cada gen.

GSEA_genes <- gseKEGG(geneList     = fold_change, organism     = "hsa", nPerm= 1000,         minGSSize    = 10, pvalueCutoff = 1, verbose      = FALSE)
## Warning in .GSEA(geneList = geneList, exponent = exponent, minGSSize =
## minGSSize, : We do not recommend using nPerm parameter incurrent and future
## releases
## Warning in fgsea(pathways = geneSets, stats = geneList, nperm = nPerm, minSize
## = minGSSize, : You are trying to run fgseaSimple. It is recommended to use
## fgseaMultilevel. To run fgseaMultilevel, you need to remove the nperm argument
## in the fgsea function call.
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (13.51% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
head(GSEA_genes)
##                ID             Description setSize enrichmentScore        NES
## hsa05200 hsa05200      Pathways in cancer      14       0.4957503  1.1855711
## hsa05205 hsa05205 Proteoglycans in cancer      10      -0.2913877 -0.7315382
## hsa05206 hsa05206     MicroRNAs in cancer      11       0.2692308  0.6200320
##             pvalue  p.adjust   qvalues rank                    leading_edge
## hsa05200 0.2482270 0.7446809 0.7446809   15  tags=71%, list=41%, signal=68%
## hsa05205 0.8090129 0.8707865 0.8707865    6  tags=40%, list=16%, signal=46%
## hsa05206 0.8707865 0.8707865 0.8707865   30 tags=100%, list=81%, signal=27%
##                                              core_enrichment
## hsa05200     4312/7422/4292/4193/387/6774/6789/595/7157/5728
## hsa05205                                79923/7099/2247/7291
## hsa05206 7422/4193/387/6774/595/7157/5728/596/5243/8434/6935
#A partir de los datos observados, podemos generar un gráfico GSEA para cada una de las asociaciones

#Aquí lo haremos para la primera asociación
gseaplot(GSEA_genes,geneSetID=1,title=GSEA_genes$Description[1])

#Para la segunda asociación
gseaplot(GSEA_genes,geneSetID=2,title=GSEA_genes$Description[2])

#Y para la tercera asociación
gseaplot(GSEA_genes,geneSetID=2,title=GSEA_genes$Description[3])

Espero puedan replicar los resultados mostrados, en caso tengan algún problema al correr este código, envíe la captura de pantalla a