clusterProfiler

Author: GuangchuangYu

02-Enrichment.qmd

@@ -54,4 +54,6 @@ GO and KEGG are the most frequently used for functional analysis. They are typic
 
 Other gene sets include but are not limited to Disease Ontology ([DO](http://disease-ontology.org/)), Disease Gene Network ([DisGeNET](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4397996/)), [wikiPathways](https://www.wikipathways.org), Molecular Signatures Database ([MSigDb](http://software.broadinstitute.org/gsea/msigdb)). 
 
+In the current `clusterProfiler` ecosystem, classical ORA and GSEA are now complemented by enrichit-backed topology-aware and multi-omics workflows. The algorithm engine lives in `enrichit`, while `clusterProfiler` provides high-level interfaces such as `nseGO()`, `mnseGO()`, `nseKEGG()`, and `mnseKEGG()` for biological interpretation.
+
 

020-enrichment-algorithms.qmd

@@ -2,6 +2,7 @@
 
 ```{r include=FALSE}
 knitr::opts_chunk$set(message=FALSE, warning=FALSE, eval=FALSE)
+source("_common.R")
 ```
 
 Functional enrichment analysis is a staple in bioinformatics for interpreting lists of genes identified from omics experiments. Behind the scenes of our ecosystem (`clusterProfiler`, `DOSE`, etc.), the heavy lifting is powered by the `enrichit` package. 
@@ -10,6 +11,8 @@ We built `enrichit` with a clear philosophy: **minimal dependencies, maximum per
 
 In this chapter, we will walk through the algorithmic backbone of enrichment analysis, moving from the standard ORA and GSEA to our newly introduced Network-based, multi-omics, and Weighted enrichment frameworks.
 
+At the package architecture level, `enrichit` serves as the algorithm engine, while `clusterProfiler` provides the knowledge-base-aware high-level entry points. In practice, this means that topology-aware workflows are now available through wrappers such as `nseGO()`, `mnseGO()`, `nseKEGG()`, and `mnseKEGG()`, while multi-omics helpers such as `aggregate_omics()` and `aggregate_enrichment()` are also re-exported by `clusterProfiler` for downstream biological interpretation.
+
 ## Over-Representation Analysis (ORA) {#ora-algorithm}
 
 Over Representation Analysis (ORA) [@boyle2004] is a widely used approach to determine whether known biological functions or processes are over-represented (= enriched) in an experimentally-derived gene list, *e.g.* a list of differentially expressed genes (DEGs).
@@ -96,6 +99,7 @@ Our core philosophy is to leverage existing strengths seamlessly:
 One RWR diffusion + one GSEA test. You get not only network-aware scores but also rigorous p-values and FDRs, with blazing speed.
 
 ```{r}
+#| eval: false
 library(clusterProfiler)
 
 # 1. Fetch the full background PPI network for Human (taxID 9606)
@@ -118,6 +122,7 @@ Historically, RWR mathematically only accepts non-negative energy. This means yo
 To solve this, `enrichit` introduces a dual-track propagation architecture (`mode = "signed"`). When you pass signed statistics (like `-log10(p) * sign(logFC)`), the algorithm automatically splits the input into up-regulated and down-regulated vectors. It runs two independent RWRs on the same network, and then subtracts the two diffused energy vectors to reconstruct a fully signed, network-aware score vector. This seamlessly bridges bidirectional network diffusion with dual-tail GSEA.
 
 ```{r}
+#| eval: false
 # Using signed statistics (e.g., -log10(p) * sign(logFC))
 res_signed <- nsea(
     geneList = signed_stats, 
@@ -182,6 +187,45 @@ mnsea_res <- mnsea(
 )
 ```
 
+The same topology-aware engine can now be accessed from `clusterProfiler` using knowledge-base-aware wrappers. For example, the following high-level workflow reuses a single shared example object and performs GO- and KEGG-aware network enrichment without manually constructing the annotation backend.
+
+```{r}
+library(clusterProfiler)
+library(org.Hs.eg.db)
+
+demo <- clusterprofiler_enrichit_demo()
+
+go_nse <- nseGO(
+  geneList = demo$geneList_evidence,
+  network = demo$network,
+  OrgDb = org.Hs.eg.db,
+  keyType = "ENTREZID",
+  ont = "BP",
+  minGSSize = 5,
+  maxGSSize = 500,
+  threshold = 1e-6,
+  maxIter = 50,
+  verbose = FALSE,
+  pvalueCutoff = 1,
+  method = "sample",
+  nPerm = 30
+)
+
+kegg_nse <- nseKEGG(
+  geneList = demo$geneList_evidence,
+  network = demo$network,
+  organism = demo$kegg_gson,
+  minGSSize = 5,
+  maxGSSize = 500,
+  threshold = 1e-6,
+  maxIter = 50,
+  verbose = FALSE,
+  pvalueCutoff = 1,
+  method = "sample",
+  nPerm = 30
+)
+```
+
 ### Explanation-ready Outputs for Multi-layer Results
 
 Once you start doing topology-aware enrichment, the next question is obvious: *which layer actually drove the hit?* Instead of baking plotting logic into the engine, `enrichit` prepares explanation-ready tables that downstream packages such as `enrichplot` can consume directly.
@@ -230,6 +274,7 @@ By passing a `weight` vector (e.g., network centrality), the model calculates an
 For GSEA, we seamlessly fuse your ranked statistics (like $\log FC$) with your external weights before passing them to the C++ engine. This amplifies the signals of core nodes while attenuating unreliable ones, without reinventing the wheel.
 
 ```{r}
+#| eval: false
 library(igraph)
 
 # Fetch network as an igraph object to calculate topological weights
@@ -257,6 +302,7 @@ The beauty of weighted enrichment lies in its **infinite possibilities**. Your `
 Finally, for ORA results, `enrichit` provides a Bayesian framework (`bayes_enrich()`) to estimate the posterior probability of explanatory terms. This helps to select the most likely active pathways while automatically penalizing redundancy among highly overlapping GO terms.
 
 ```{r}
+#| eval: false
 # After obtaining an ORA result object
 bayes_res <- bayes_enrich(ora_result)
 summary_res <- bayes_summary(bayes_res)
@@ -422,6 +468,8 @@ early_ora_res <- ora(
 
 By decoupling ID mapping, statistical integration, and pathway enrichment, `enrichit` provides a clean, robust, and highly extensible pipeline for multi-omics functional analysis.
 
+These workflow helpers are also available from `clusterProfiler`, which means the aggregated result can be sent directly to high-level functions such as `enrichGO()`, `gseGO()`, or `nseGO()` without changing the underlying multi-omics logic.
+
 ### 5. Late Fusion at the Pathway Level
 
 Early fusion is elegant when your features align cleanly. But sometimes that assumption collapses immediately: transcriptomics might have a background of 20,000 genes, proteomics might only detect 4,000 proteins, and region-based omics may not even live in gene space before annotation. In these cases, forcing everything into one feature matrix is a great way to manufacture missing values and mapping noise.
@@ -465,6 +513,7 @@ After running a multi-omics integration, you may want to know whether a signific
 First, you can classify the pattern of each pathway by comparing the merged results with single-omics results:
 
 ```{r}
+#| eval: false
 # Run ORA for individual omics
 rna_ora_res <- ora(
     gene = omics_gene_ids[rna_pvals < 0.05],
@@ -492,6 +541,7 @@ The `Omics_Pattern` column will label the pathway as `Shared`, `RNA-Specific`, `
 Then, you can extract the gene-level contribution for a specific pathway to see the original statistics of the core enrichment genes:
 
 ```{r}
+#| eval: false
 # Get original multi-omics statistics for genes in the pathway
 contribution_df <- get_omics_contribution(
     res = classified_ora_res, 

021-go.qmd

@@ -145,9 +145,66 @@ ego3 <- gseGO(geneList     = geneList,
 
 The format of input data, `geneList`, is documented in the @sec-id-convert. Beware that only gene sets with size in `[minGSSize, maxGSSize]` will be tested.
 
+## Topology-aware GO enrichment
+
+In classical GO enrichment, genes are treated as independent members of a gene set. When a biological network is available, `clusterProfiler` can now call the `enrichit` engine through the high-level wrappers `nseGO()` and `mnseGO()`. The former performs single-network enrichment, while the latter extends the same idea to a multi-layer setting.
+
+To keep the example reproducible and consistent with the KEGG chapter, we reuse one shared demonstration object.
+
+```{r}
+#| label: topology-demo-go
+demo <- clusterprofiler_enrichit_demo()
+```
+
+```{r}
+#| label: nseGO-example
+ego_nse <- nseGO(
+  geneList = demo$geneList_evidence,
+  network = demo$network,
+  OrgDb = org.Hs.eg.db,
+  keyType = "ENTREZID",
+  ont = "BP",
+  minGSSize = 5,
+  maxGSSize = 500,
+  threshold = 1e-6,
+  maxIter = 50,
+  verbose = FALSE,
+  pvalueCutoff = 1,
+  method = "sample",
+  nPerm = 30
+)
+head(as.data.frame(ego_nse)[, c("ID", "Description", "NES", "p.adjust")], 3)
+```
+
+The `mnseGO()` function keeps the same GO-aware interface, but accepts a list of layer-specific networks together with explicit inter-layer couplings.
+
+```{r}
+#| label: mnseGO-example
+ego_mnse <- mnseGO(
+  seed_list = demo$seed_list,
+  networks = list(RNA = demo$network, PROT = demo$network_2),
+  couplings = demo$couplings,
+  OrgDb = org.Hs.eg.db,
+  keyType = "ENTREZID",
+  ont = "BP",
+  collapse = "weighted_mean",
+  layer_weights = c(RNA = 1, PROT = 1.2),
+  minGSSize = 5,
+  maxGSSize = 500,
+  threshold = 1e-6,
+  maxIter = 50,
+  verbose = FALSE,
+  pvalueCutoff = 1,
+  method = "sample",
+  nPerm = 30
+)
+head(as.data.frame(ego_mnse)[, c("ID", "Description", "NES", "p.adjust")], 3)
+```
+
+These wrappers preserve the familiar GO semantics (`OrgDb`, `ont`, and `keyType`) while delegating the network propagation and enrichment engine to `enrichit`.
+
 
 
-## GO analysis for non-model organisms {#clusterprofiler-go-non-model}
 
 Both the `enrichGO()` and `gseGO()` functions require an `OrgDb` object as the background annotation. For organisms that don't have `OrgDb` provided by `Bioconductor`, users can query one (if available) online via `r Biocpkg("AnnotationHub")`. If there is no `OrgDb` available, users can obtain GO annotation from other sources, e.g. from `r Biocpkg("biomaRt")`, or annotate the genes using [Blast2GO](https://www.blast2go.com/) or the [Trinotate](https://rnabio.org/module-07-trinotate/0007/02/01/Trinotate/) pipeline. Then the `enricher()` or `GSEA()` functions can be used to perform GO analysis for these organisms, similar to the examples using WikiPathways and MSigDB. Another solution is to create an `OrgDb` on your own using the `r Biocpkg("AnnotationForge")` package.
 

022-kegg.qmd

@@ -23,6 +23,7 @@ The `r Biocpkg("clusterProfiler")` package provides the `search_kegg_organism()`
 
 ```{r}
 #| label: search-kegg-organism
+#| eval: false
 library(clusterProfiler)
 search_kegg_organism('ece', by='kegg_code')
 ecoli <- search_kegg_organism('Escherichia coli', by='scientific_name')
@@ -205,12 +206,15 @@ This ensures that the analysis runs offline and is fully reproducible using the
 
 The `r Biocpkg("clusterProfiler")` package provides the `bitr_kegg()` function to support ID conversion via the KEGG API. This is particularly useful for species that do not have an `OrgDb` object but are supported by KEGG.
 
+The examples in this section require live KEGG REST access. To keep the book build reproducible when the KEGG API is slow or unreachable, we show them without executing them during rendering.
+
 ### Basic Usage
 
 Here is an example of converting KEGG IDs to NCBI Protein IDs and then to UniProt IDs for human genes.
 
 ```{r}
 #| label: bitr-kegg
+#| eval: false
 library(clusterProfiler)
 data(gcSample)
 hg <- gcSample[[1]]
@@ -225,12 +229,13 @@ np2up <- bitr_kegg(eg2np[,2], fromType='ncbi-proteinid', toType='uniprot', organ
 head(np2up)
 ```
 
-The ID type (both `fromType` and `toType`) should be one of 'kegg', 'ncbi-geneid', 'ncbi-proteinid', or 'uniprot'. The 'kegg' ID is the primary ID used in the KEGG database. A rule of thumb is that 'kegg' ID corresponds to `entrezgene` ID for eukaryote species and `Locus` ID for prokaryotes.
+The ID type (both `fromType` and `toType`) should be one of 'kegg', 'ko', 'ncbi-geneid', 'ncbi-proteinid', or 'uniprot'. The 'kegg' ID is the primary gene ID used in the KEGG database, while 'ko' refers to KEGG Orthology identifiers. A rule of thumb is that 'kegg' gene ID corresponds to `entrezgene` ID for eukaryote species and `Locus` ID for prokaryotes.
 
 For many prokaryote species, `entrezgene` IDs are not available. For example, `ece:Z5100` (E. coli O157:H7 EDL933) has `NCBI-ProteinID` and `UniProt` links, but not `NCBI-GeneID`. Attempting to convert `Z5100` to `ncbi-geneid` will result in an error stating that `ncbi-geneid` is not supported. However, conversion to `ncbi-proteinid` or `uniprot` is possible.
 
 ```{r}
 #| label: bitr-kegg-prokaryote
+#| eval: false
 # Example for prokaryote ID conversion
 bitr_kegg("Z5100", fromType="kegg", toType='ncbi-proteinid', organism='ece')
 bitr_kegg("Z5100", fromType="kegg", toType='uniprot', organism='ece')
@@ -242,13 +247,15 @@ A common confusion exists between `K` numbers (KEGG Orthology) and `ko` numbers
 
 ```{r}
 #| label: bitr-kegg-path-ko
+#| eval: false
 bitr_kegg("K00844", "kegg", "Path", "ko")
 ```
 
 To retrieve the names of the pathways (e.g., "Glutathione metabolism" instead of "ko00480"), we can define a simple helper function `ko2name` to query the KEGG API.
 
 ```{r}
 #| label: ko2name-example
+#| eval: false
 x <- bitr_kegg("K00799", "kegg", "Path", "ko")
 y <- ko2name(x$Path)
 merge(x, y, by.x='Path', by.y='ko')
@@ -258,6 +265,7 @@ merge(x, y, by.x='Path', by.y='ko')
 
 
 ```{r}
+#| eval: false
 data(geneList, package="DOSE")
 gene <- names(geneList)[abs(geneList) > 2]
 
@@ -274,6 +282,7 @@ Input ID type can be `kegg`, `ncbi-geneid`, `ncbi-proteinid` or `uniprot` (see a
 For GO analysis, the `readable` parameter controls whether to translate the IDs to human-readable gene names. This parameter is not available for KEGG analysis. However, the `setReadable()` function can translate input gene IDs to gene names if the corresponding `OrgDb` object is available.
 
 ```{r}
+#| eval: false
 library(org.Hs.eg.db)
 # 'kk' is the enrichKEGG result from the previous section
 kk <- setReadable(kk, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
@@ -283,6 +292,7 @@ head(kk)
 ## KEGG pathway gene set enrichment analysis {#clusterprofiler-kegg-pathway-gsea}
 
 ```{r}
+#| eval: false
 kk2 <- gseKEGG(geneList     = geneList,
                organism     = 'hsa',
                minGSSize    = 120,
@@ -291,6 +301,58 @@ kk2 <- gseKEGG(geneList     = geneList,
 head(kk2)
 ```
 
+## Topology-aware KEGG enrichment
+
+The same high-level design is now available for topology-aware KEGG analysis. `nseKEGG()` performs single-network enrichment and `mnseKEGG()` extends the workflow to multiple network layers, while keeping the familiar `clusterProfiler` interface.
+
+For a lightweight and reproducible example, we reuse the same shared demo object as in the GO chapter and pass a local KEGG-style `GSON` object through the `organism` argument.
+
+```{r}
+#| label: topology-demo-kegg
+demo <- clusterprofiler_enrichit_demo()
+```
+
+```{r}
+#| label: nseKEGG-example
+kk_nse <- nseKEGG(
+  geneList = demo$geneList_evidence,
+  network = demo$network,
+  organism = demo$kegg_gson,
+  minGSSize = 5,
+  maxGSSize = 500,
+  threshold = 1e-6,
+  maxIter = 50,
+  verbose = FALSE,
+  pvalueCutoff = 1,
+  method = "sample",
+  nPerm = 30
+)
+head(as.data.frame(kk_nse)[, c("ID", "Description", "NES", "p.adjust")], 3)
+```
+
+```{r}
+#| label: mnseKEGG-example
+kk_mnse <- mnseKEGG(
+  seed_list = demo$seed_list,
+  networks = list(RNA = demo$network, PROT = demo$network_2),
+  couplings = demo$couplings,
+  organism = demo$kegg_gson,
+  collapse = "weighted_mean",
+  layer_weights = c(RNA = 1, PROT = 1.2),
+  minGSSize = 5,
+  maxGSSize = 500,
+  threshold = 1e-6,
+  maxIter = 50,
+  verbose = FALSE,
+  pvalueCutoff = 1,
+  method = "sample",
+  nPerm = 30
+)
+head(as.data.frame(kk_mnse)[, c("ID", "Description", "NES", "p.adjust")], 3)
+```
+
+In a real analysis, users can pass a KEGG organism code such as `"hsa"` directly. Here we use a local `GSON` object to keep the example self-contained and reproducible.
+
 ## KEGG Pathway Classification
 
 KEGG pathways are organized into a hierarchical structure with top-level categories such as "Metabolism", "Genetic Information Processing", "Environmental Information Processing", "Cellular Processes", "Organismal Systems", "Human Diseases", and "Drug Development". This classification helps users interpret enrichment results by providing higher-level biological context. For instance, users might be interested specifically in signaling pathways ("Environmental Information Processing") or metabolic pathways ("Metabolism") while excluding disease-related pathways.
@@ -310,6 +372,7 @@ To prevent this cached data from becoming outdated, the package utilizes **GitHu
 
 ```{r}
 #| label: kegg-module-example
+#| eval: false
 mkk <- enrichMKEGG(gene = gene,
                    organism = 'hsa',
                    pvalueCutoff = 1,
@@ -321,6 +384,7 @@ head(mkk)
 
 
 ```{r}
+#| eval: false
 mkk2 <- gseMKEGG(geneList = geneList,
                  organism = 'hsa',
                  pvalueCutoff = 1)