6  Network-aware and weighted enrichment

This chapter extends classical enrichment when interaction topology, propagation, or gene-level detection weights are part of the evidence model. It is separate from ordinary PPI visualization: the network enters the statistical ranking or weighting step before enrichment is tested.

6.1 Network-based Set Enrichment Analysis (NSEA)

Traditional enrichment methods treat pathways as independent “bags of genes,” completely ignoring the biological interactions between them. If you only look at the overlap, you are wasting a massive amount of network topological information.

Some previous methods (like EnrichNet) attempted to solve this by calculating network distances. However, they suffered from a fatal flaw: they only provided a relative score without rigorous p-values, and calculating empirical p-values via network permutation was excruciatingly slow. Life is too short to wait for permutations.

To break this bottleneck, we introduced Network-based Set Enrichment Analysis (NSEA), a true Network-based GSEA.

Our core philosophy is to leverage existing strengths seamlessly:

  1. Propagation: We take your input gene scores and diffuse them across a biological network using a Random Walk with Restart (RWR). To squeeze out every drop of performance, the RWR is powered by RcppEigen sparse matrix multiplication in C++. Even networks with hundreds of thousands of edges converge in a flash.
  2. Enrichment: After diffusion, every node in the network receives a network-aware score. The default NSEA significance route (significance = "whole_pipeline") builds its null by permuting the input labels, re-running the RWR, and recomputing the enrichment score for each permutation. This keeps the network diffusion inside the null model instead of treating the smoothed scores as exchangeable ordinary GSEA statistics.

The practical consequence is that NSEA is not simply “one diffusion followed by one ordinary GSEA test”: the observed score is produced by one diffusion, while the default p-values and NES are calibrated from repeated whole-pipeline null runs. significance = "internal" is retained as a legacy alternative that applies the internal GSEA permutation step to already diffused scores; it should not be interpreted as the same null model.

library(clusterProfiler)

# 1. Fetch the full background PPI network for Human (taxID 9606)
# Default output is "data.frame" (edge list), which nsea() natively supports
net <- get_ppi_network(9606)

# 2. Run NSEA
# (Assuming 'stats' is a named numeric vector of non-negative evidence scores)
nsea_result <- nsea(
    geneList = stats,
    network = net,
    gene_sets = pathways
)

6.1.1 Direction-aware Propagation (Signed Mode)

Historically, RWR mathematically only accepts non-negative energy. This means you could only pass absolute statistics (like -log10(p)), and the result would only tell you if a pathway was “hit”, but not whether it was up-regulated or down-regulated.

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.

# Use the offline demo: signed statistics, network edges and gene sets share
# one Entrez namespace. The full-STRING version above is the network-dependent
# variant; this one runs in the book build.
demo <- clusterprofiler_enrichit_demo()
gene_sets <- split(demo$kegg_gson@gsid2gene$gene,
                   demo$kegg_gson@gsid2gene$gsid)
names(gene_sets) <- demo$kegg_gson@gsid2name$name[
    match(names(gene_sets), demo$kegg_gson@gsid2name$gsid)]

res_signed <- nsea(
    geneList = demo$geneList_signed,
    network = demo$network,
    gene_sets = gene_sets,
    mode = "signed"
)
head(as.data.frame(res_signed)[, c("ID", "enrichmentScore", "NES", "pvalue", "p.adjust")])
                                                               ID
Cell cycle (demo)                               Cell cycle (demo)
Insulin signaling pathway (demo) Insulin signaling pathway (demo)
MAPK signaling pathway (demo)       MAPK signaling pathway (demo)
                                 enrichmentScore       NES    pvalue  p.adjust
Cell cycle (demo)                      0.3894050 1.4091496 0.1558442 0.3716284
Insulin signaling pathway (demo)       0.2797944 1.2059385 0.2477522 0.3716284
MAPK signaling pathway (demo)          0.2251367 0.9165147 0.5074925 0.5074925

6.1.2 Multi-layer Network-based Enrichment (MNSEA)

Single-layer propagation is already useful, but real multi-omics systems are rarely confined to one graph. Transcriptomics, proteomics, and phosphoproteomics often live on related yet non-identical network layers. Forcing everything into a single collapsed graph throws away layer-specific topology before the analysis even begins.

To address this, enrichit introduces multi-layer NSEA through mnsea(). The algorithm builds a supra-graph from layer-specific networks plus explicit inter-layer couplings, performs Random Walk with Restart on the combined topology, and then collapses layer-wise diffusion scores into a unified ranked vector for GSEA. In other words, it is still “propagate first, enrich later”, but now the propagation itself respects the multi-layer structure.

# Small multi-layer example with RNA and Protein layers
networks <- list(
  RNA = data.frame(
    from = c("Gene1", "Gene2", "Gene3"),
    to = c("Gene2", "Gene3", "Gene4"),
    weight = c(1, 1, 1)
  ),
  PROT = data.frame(
    from = c("Gene1", "Gene3", "Gene4"),
    to = c("Gene3", "Gene4", "Gene1"),
    weight = c(1, 1, 1)
  )
)

couplings <- data.frame(
  from_layer = c("RNA", "RNA", "RNA"),
  from_id = c("Gene1", "Gene3", "Gene4"),
  to_layer = c("PROT", "PROT", "PROT"),
  to_id = c("Gene1", "Gene3", "Gene4"),
  weight = c(1, 1, 1)
)

seed_list <- list(
  RNA = c(Gene1 = 1.5, Gene2 = 0.8, Gene4 = -0.9),
  PROT = c(Gene1 = 0.5, Gene3 = 1.2, Gene4 = -0.4)
)

ml_gene_sets <- list(
  SignalPath = c("Gene1", "Gene2", "Gene3"),
  DownPath = c("Gene4")
)

mnsea_res <- mnsea(
  seed_list = seed_list,
  networks = networks,
  couplings = couplings,
  gene_sets = ml_gene_sets,
  mode = "signed",
  collapse = "weighted_mean",
  layer_weights = c(RNA = 1, PROT = 1.5),
  minGSSize = 1,
  maxGSSize = 10,
  method = "sample",
  nPerm = 30,
  verbose = FALSE
)
head(as.data.frame(mnsea_res)[, c("ID", "NES", "pvalue", "setSize")])
                   ID       NES    pvalue setSize
DownPath     DownPath -1.159091 0.6111111       1
SignalPath SignalPath  1.000000 1.0000000       3

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.

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
)
head(as.data.frame(go_nse)[, c("ID", "Description", "setSize", "NES", "pvalue", "p.adjust")])
                   ID                                       Description setSize
GO:0032496 GO:0032496                    response to lipopolysaccharide      10
GO:0051656 GO:0051656           establishment of organelle localization      10
GO:0071216 GO:0071216              cellular response to biotic stimulus       8
GO:0071219 GO:0071219 cellular response to molecule of bacterial origin       8
GO:0071222 GO:0071222           cellular response to lipopolysaccharide       8
GO:0033993 GO:0033993                                 response to lipid      22
                NES     pvalue  p.adjust
GO:0032496 3.536267 0.03225806 0.4833589
GO:0051656 3.213823 0.03225806 0.4833589
GO:0071216 3.147129 0.03225806 0.4833589
GO:0071219 3.147129 0.03225806 0.4833589
GO:0071222 3.147129 0.03225806 0.4833589
GO:0033993 3.145848 0.03225806 0.4833589
head(as.data.frame(kegg_nse)[, c("ID", "Description", "setSize", "NES", "pvalue", "p.adjust")])
                           ID                      Description setSize      NES
hsa_demo_04010 hsa_demo_04010    MAPK signaling pathway (demo)      50 8.043335
hsa_demo_04110 hsa_demo_04110                Cell cycle (demo)      35 6.512758
hsa_demo_04910 hsa_demo_04910 Insulin signaling pathway (demo)      60 5.285313
                   pvalue   p.adjust
hsa_demo_04010 0.03225806 0.03225806
hsa_demo_04110 0.03225806 0.03225806
hsa_demo_04910 0.03225806 0.03225806

6.1.3 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.

The mnseaResult object stores precomputed pathway-level and feature-level contribution caches, and also supports extracting a pathway-specific explanation subnetwork.

# Pathway-level contribution table
pathway_tbl <- get_mnsea_contribution(
  mnsea_res,
  level = "pathway"
)

# Feature-level contribution table for one pathway
feature_tbl <- get_mnsea_contribution(
  mnsea_res,
  pathway_id = pathway_tbl$ID[1],
  level = "feature"
)

# Extract an explanation-ready subnetwork
subnet <- extract_mnsea_subnetwork(
  mnsea_res,
  pathway_id = pathway_tbl$ID[1]
)
pathway_tbl
          ID Description layer contribution     share n_feature
1   DownPath    DownPath   RNA    0.7720049 0.5306919         1
2   DownPath    DownPath  PROT    0.6827092 0.4693081         1
3 SignalPath  SignalPath   RNA    0.2075870 0.5174437         3
4 SignalPath  SignalPath  PROT    0.1935909 0.4825563         3
head(feature_tbl)
        ID Description Feature layer      score abs_score is_core
1 DownPath    DownPath   Gene4   RNA -0.7720049 0.7720049    TRUE
2 DownPath    DownPath   Gene4  PROT -0.6827092 0.6827092    TRUE

6.2 Weighted Enrichment Analysis

Have you ever noticed a weird phenomenon? No matter what odd phenotype you study, the enriched pathways always seem to feature the same “usual suspects”—like the ribosome, spliceosome, or broad immune responses.

This is driven by a systematic bias. In real biological networks, genes are not created equal. Hub genes participate in numerous functions and are easily perturbed. If you don’t account for this, traditional ORA will be easily fooled by these highly active Hub genes, yielding false-positive significance for broad pathways.

To correct this, enrichit introduces Weighted Enrichment Analysis, supporting both ORA and GSEA.

6.2.1 Weighted ORA: De-biasing with Wallenius Distribution

How do we correct the bias in ORA? We bring in the Wallenius’ noncentral hypergeometric distribution.

By passing a weight vector (e.g., network centrality), the model calculates an Odds ratio for each pathway. If a pathway is full of high-weight core genes, the model applies a penalty to the p-value. The false significance driven by Hub genes is brought back to reality, allowing truly specific pathways to stand out.

Note: In the early days of RNA-seq, goseq used gene length as a weight to correct ORA. However, in the era of 3’ poly-A sequencing and 10x single-cell, length bias is mostly obsolete. Our implementation is universal—your weights are dictated by your scientific question, not hardcoded to gene length.

6.2.2 Weighted GSEA: Fusing Weights with Statistics

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.

library(igraph)
library(enrichit)

# The offline demo keeps the gene list, the network nodes and the gene sets in
# the same Entrez namespace — the property the original example was missing.
demo <- clusterprofiler_enrichit_demo()

# ora()/gsea() take a named list of gene sets; the demo GSON uses the same IDs.
gene_sets <- split(demo$kegg_gson@gsid2gene$gene,
                   demo$kegg_gson@gsid2gene$gsid)
names(gene_sets) <- demo$kegg_gson@gsid2name$name[
    match(names(gene_sets), demo$kegg_gson@gsid2name$gsid)]

# Any gene-level weight works. Here, PageRank centrality on the demo network:
net <- graph_from_data_frame(demo$network[, c("from", "to")], directed = FALSE)
w <- page_rank(net)$vector      # names = Entrez IDs, same as the gene list

de_genes <- names(demo$geneList_evidence)[1:100]
res_wora <- ora(gene = de_genes,
                gene_sets = gene_sets,
                universe = names(demo$geneList_evidence),
                weight = w)

res_wgsea <- gsea(geneList = demo$geneList_signed,
                  gene_sets = gene_sets,
                  weight = w)
head(as.data.frame(res_wora)[, c("ID", "SetSize", "Count", "pvalue")])
                                ID SetSize Count       pvalue
1    MAPK signaling pathway (demo)      50    50 0.000000e+00
2                Cell cycle (demo)      35    35 0.000000e+00
3 Insulin signaling pathway (demo)      60    40 2.907958e-09
head(as.data.frame(res_wgsea)[, c("ID", "setSize", "enrichmentScore", "NES", "pvalue")])
                                ID setSize enrichmentScore      NES
1    MAPK signaling pathway (demo)      50       0.4811736 3.322632
2                Cell cycle (demo)      35       0.5476000 3.280942
3 Insulin signaling pathway (demo)      60       0.4085579 2.953280
        pvalue
1 1.936519e-08
2 6.586029e-08
3 1.783112e-06
NoteWeighted ORA needs BiasedUrn installed

BiasedUrn is a suggested package of enrichit, so a default install does not bring it in. Weighted ORA calls BiasedUrn::pWNCHypergeo() and stops when it is missing — it does not fall back to ordinary ORA:

Error: The package "BiasedUrn" is required for weighted ORA.

Install it with install.packages("BiasedUrn") before running this section. Weighted GSEA is unaffected; the weights are fused into the ranked statistics before the C++ engine, so no extra package is needed.

Warningweight must share the namespace of the genes

weight is matched by name, not by position. Its names must be the same identifiers as gene / geneList / universe — here, Entrez IDs. If they do not overlap:

  • ora() only warns (“Some background genes do not have weights. Assigning median weight.”) and continues with a fallback weighting, so the result still looks plausible;
  • gsea() stops with “No common genes found between geneList and weight”.

That asymmetry is easy to miss, especially in a chunk or session with warning = FALSE. Check the overlap explicitly before trusting a weighted result:

stopifnot(any(names(w) %in% names(demo$geneList_evidence)))

The demo network is deliberately small, so its centrality is nearly uniform; the numbers here demonstrate the calling contract rather than a biological effect. Use your own network and gene-level scores for a real analysis.

The beauty of weighted enrichment lies in its infinite possibilities. Your weight doesn’t have to be network centrality. It can be a tissue specificity score, an evolutionary conservation score, or a multi-omics confidence score. If you have prior knowledge that “these genes are more important”, simply pass them to the weight parameter and let enrichit handle the rest.

6.3 Bayesian term selection

Bayesian term selection is documented as a dedicated evidence-integration capability in Bayesian term selection. The algorithm backend can be used after ORA, but term selection is not itself a network-enrichment test.

6.4 Next steps

References