Traditional enrichment analysis typically results in a list of significant pathways or GO terms. While statistically sound, these lists often leave researchers asking, “So what?” What is the underlying biological mechanism? Who are the key drivers? Is this a pro-survival or pro-death signal?
To bridge the gap between statistical results and biological insights, clusterProfiler introduces an AI-powered interpretation module. By leveraging Large Language Models (LLMs) and a multi-agent system, clusterProfiler can now act as a virtual bioinformatician, converting dry enrichment lists into coherent, evidence-based biological narratives.
interpret FunctionThe core function for this feature is interpret(). It accepts enrichment results (e.g., from enrichKEGG, enrichGO, or compareCluster) and uses an LLM to generate a structured report.
To use this feature, you need to configure an API key for a provider supported by aisdk. Model selection is now delegated to aisdk: if model is not specified, interpret() uses the current default model managed by aisdk::get_model() / aisdk::set_model(), while an explicit model = ... argument overrides that default for a single call. In practice, this usually means loading credentials from your .env file or setting the relevant environment variables before calling set_model(). For provider configuration details, see the aisdk documentation. The examples on this page use DeepSeek, whose key aisdk reads from the DEEPSEEK_API_KEY environment variable; any other provider supported by aisdk can be used instead.
library(aisdk)
# Credentials come from the environment; a local .env file works too
if (file.exists(".env") && requireNamespace("dotenv", quietly = TRUE)) {
dotenv::load_dot_env()
}
# Set a session-wide default model. Any provider supported by aisdk works here;
# in CI the key is supplied as the DEEPSEEK_API_KEY repository secret.
set_model("deepseek:deepseek-v4-flash")
# Check the current default model
get_model()[1] "deepseek:deepseek-v4-flash"library(clusterProfiler)
library(DOSE)
# A self-contained enrichment result: `edo` is built offline from the demo gene
# list shipped with DOSE, so the example depends neither on a previous chapter
# (`edo` is defined in the enrichplot chapter) nor on a network service.
data(geneList, package = "DOSE")
de <- names(geneList)[abs(geneList) > 1.5]
edo <- enrichDO(de)
# Basic usage
# 'edo' is your enrichment result object
res <- interpret(edo)
if (interpret_failed(res)) {
stop("interpret() returned no report -- check the API key and quota:\n",
res$overview)
}
print(res)## Interpretation Result
### Cluster: Default
### 1. Overview
The enrichment is dominated by two intersecting programs: (i) a cell-cycle and DNA-replication module highlighted by Meier-Gorlin syndrome and multiple epithelial cancers (stomach cancer, breast cancer, esophageal cancer, bladder cancer, female reproductive organ cancer, in situ carcinoma) and (ii) an inflammatory and extracellular matrix (ECM) remodeling/cardiovascular module highlighted by myocarditis, cardiomyopathy, systemic scleroderma, rheumatic disease, viral infectious disease, type 2 diabetes, and congestive heart failure. Top marker genes include FOXM1, MYBL2, CDC20, CDC45, MCM10, CDCA8, KIF23, CENPE, MELK, CCNB2, NDC80, NCAPH, TOP2A alongside S100A8, S100A9, MMP1, MMP9, BCL2A1, MARCO, and NMU.
### 2. Regulatory Drivers
FOXM1 (2305) and MYBL2 (4605) are the strongest transcriptional drivers inferred from marker status. Both are cell-cycle transcription factors and co-occur with target genes such as CDC20, CDC45, MCM10, CCNB2, TOP2A, KIF23, CENPE, NDC80, NCAPH, CDCA8, and MELK in proliferative disease terms (Meier-Gorlin syndrome, stomach cancer, breast cancer, esophageal cancer, bladder cancer, female reproductive organ cancer, in situ carcinoma). S100A8/S100A9 (6279/6280) are not transcription factors but act as inflammatory alarmins and master regulators of myeloid inflammation, and they are markers in myocarditis and viral infectious disease.
### 3. Key Mechanisms
Mechanism 1 is DNA replication licensing and mitotic proliferation. Meier-Gorlin syndrome enrichment includes CDC45, ORC1, CDT1, CDC6, ORC6, and MCM2, while markers include CDC45, MCM10, CDC20, CDCA8, KIF23, CENPE, MYBL2, MELK, CCNB2, NDC80, NCAPH, TOP2A, and FOXM1. This module is shared with stomach cancer, breast cancer, esophageal cancer, bladder cancer, female reproductive organ cancer, and in situ carcinoma. Mechanism 2 is inflammatory alarmin and ECM remodeling. S100A8, S100A9, MMP1, MMP9, BCL2A1, MARCO, and NMU are markers, and enriched terms include myocarditis, systemic scleroderma, rheumatic disease, scleroderma, viral infectious disease, and bronchiolitis obliterans. Mechanism 3 is cardiovascular-metabolic overlap. Cardiomyopathy, extrinsic cardiomyopathy, congestive heart failure, and type 2 diabetes share MMP1/MMP9, S100A8/S100A9, and FOXM1-related proliferative signals.
### 4. Crosstalk & Interactions
The proliferative and inflammatory modules are not independent. FOXM1/MYBL2-driven cell-cycle activation (Meier-Gorlin syndrome, stomach cancer, breast cancer, esophageal cancer, bladder cancer, female reproductive organ cancer, in situ carcinoma) can expand epithelial and immune compartments, while S100A8/S100A9 alarmins and MMP1/MMP9 proteases (myocarditis, systemic scleroderma, viral infectious disease, cardiomyopathy) create a remodeled microenvironment that can further promote cancer and cardiovascular failure. MMP1/MMP9 are shared between scleroderma, type 2 diabetes, and cancers, suggesting ECM turnover as a crosstalk node. The S100A8/S100A9 heterodimer links viral infectious disease and myocarditis. DNA replication licensing genes (CDC45, MCM10) connect Meier-Gorlin syndrome to cancer proliferation.
### 5. Hypothesis
The sample represents a dual-pathology state in which FOXM1/MYBL2 activation drives unscheduled DNA replication and mitosis (Meier-Gorlin syndrome, stomach cancer, breast cancer, esophageal cancer, bladder cancer, female reproductive organ cancer, in situ carcinoma), while S100A8/S100A9 and MMP1/MMP9-mediated inflammation and ECM remodeling drive cardiovascular, autoimmune/fibrotic, and metabolic outcomes (myocarditis, cardiomyopathy, systemic scleroderma, rheumatic disease, type 2 diabetes, congestive heart failure). The convergence of these programs may create a feed-forward loop: inflammatory alarmins and proteases remodel tissue and stimulate proliferation, and proliferating cells amplify inflammatory signaling, explaining the co-enrichment of cancer and cardiovascular/fibrotic diseases.
### 6. Narrative Draft
Enrichment analysis of the top marker genes reveals a coordinated cell-cycle and inflammatory/ECM program. Cell-cycle and DNA-replication genes, including FOXM1, MYBL2, CDC20, CDC45, MCM10, CDCA8, KIF23, CENPE, MELK, CCNB2, NDC80, NCAPH, and TOP2A, are prominent markers and are enriched in Meier-Gorlin syndrome and epithelial cancers (stomach cancer, breast cancer, esophageal cancer, bladder cancer, female reproductive organ cancer, in situ carcinoma). In parallel, S100A8/S100A9 alarmins, MMP1/MMP9 proteases, BCL2A1, MARCO, and NMU mark an inflammatory and ECM-remodeling module enriched in myocarditis, cardiomyopathy, systemic scleroderma, rheumatic disease, viral infectious disease, bronchiolitis obliterans, type 2 diabetes, and congestive heart failure. These modules converge on shared nodes: FOXM1/MYBL2 drive proliferation, MMP1/MMP9 mediate ECM turnover in both cancer and fibrotic/metabolic disease, and S100A8/S100A9 link infection to myocardial inflammation. This suggests a feed-forward interaction between proliferative and inflammatory programs in the experimental model.
### 7. Refined Regulatory Network
Key Interactions:
FOXM1 -- CDC20 (activation) - FOXM1 is a marker transcription factor and CDC20 is a marker mitotic gene; both are co-enriched in multiple cancers (stomach cancer, breast cancer, esophageal cancer, bladder cancer, female reproductive organ cancer).
FOXM1 -- CDC45 (activation) - FOXM1 drives G1/S and DNA replication genes, and CDC45 is a marker and Meier-Gorlin syndrome gene (Meier-Gorlin syndrome, stomach cancer, breast cancer).
FOXM1 -- MCM10 (activation) - FOXM1 regulates replication licensing, and MCM10 is a marker DNA replication gene co-enriched in cancers (stomach cancer, breast cancer, esophageal cancer, bladder cancer).
FOXM1 -- MYBL2 (activation) - FOXM1 and MYBL2 cooperate in cell-cycle transcription; both are top markers and co-enriched in proliferative disease terms (Meier-Gorlin syndrome, stomach cancer, breast cancer).
MYBL2 -- CCNB2 (activation) - MYBL2 is a marker transcription factor controlling G2/M genes, and CCNB2 is a marker co-enriched in cancers (breast cancer, esophageal cancer, bladder cancer, female reproductive organ cancer).
MYBL2 -- TOP2A (activation) - MYBL2 regulates TOP2A; both are markers and co-enriched in cancers and scleroderma (breast cancer, stomach cancer, systemic scleroderma).
FOXM1 -- TOP2A (activation) - FOXM1 targets TOP2A, and both markers are co-enriched in cancers and scleroderma (stomach cancer, breast cancer, systemic scleroderma).
CDC45 -- MCM10 (binding) - CDC45 and MCM10 function in the DNA replication machinery; CDC45 is a Meier-Gorlin syndrome gene and both are markers in cancers (Meier-Gorlin syndrome, stomach cancer, breast cancer).
NDC80 -- CENPE (binding) - NDC80 and CENPE are mitotic kinetochore proteins; both are markers and co-enriched in cancers (stomach cancer, breast cancer, esophageal cancer, bladder cancer).
CDCA8 -- KIF23 (functional interaction) - Both are marker mitotic/cytokinesis genes and are co-enriched in cancers (breast cancer, stomach cancer, esophageal cancer).
S100A8 -- S100A9 (binding) - S100A8 and S100A9 form the calprotectin heterodimer; both are markers and co-enriched in myocarditis and viral infectious disease (myocarditis, viral infectious disease).
S100A9 -- MMP9 (activation) - S100A9 alarmin can induce MMP9 in inflammatory microenvironments; both are markers and co-enriched in cardiomyopathy (cardiomyopathy, type 2 diabetes, systemic scleroderma).
**Network Evidence:**
The refined network captures a FOXM1/MYBL2 transcriptional module (CDC20, CDC45, MCM10, CCNB2, TOP2A, KIF23, CENPE, NDC80, NCAPH, CDCA8, MELK) supported by enrichment in Meier-Gorlin syndrome and multiple carcinomas (stomach cancer, breast cancer, esophageal cancer, bladder cancer, female reproductive organ cancer, in situ carcinoma); a DNA replication licensing module (CDC45-MCM10) supported by Meier-Gorlin syndrome and cancer; a mitotic/kinetochore module (NDC80-CENPE, CDCA8-KIF23) supported by cancer; and an inflammatory/ECM module (S100A8-S100A9 heterodimer, MMP1/MMP9) supported by myocarditis, cardiomyopathy, systemic scleroderma, and viral infectious disease. interpret() is not just for explaining enrichment results. It breaks down LLM capabilities into three distinct tasks:
task = "interpretation": (Default) Converts enrichment results into a mechanistic narrative suitable for publication (What -> So What).task = "annotation": Performs cell type annotation for single-cell clusters using both marker genes and enrichment terms as evidence.task = "phenotyping": Assigns a “state/phenotype label” to a group (e.g., “Pro-inflammatory” or “Senescent-like”).To strengthen the evidence, interpret() supports “evidence synthesis” from multiple sources:
enrichResult, gseaResult, or compareClusterResult.list() of results (e.g., list(kegg_res, go_res) or list(cellmarker_res, go_res)).compareCluster result, it automatically splits by cluster and generates a report for each.The key features of interpret() include:
For example, we can use Seurat to identify marker genes for each cluster in a single-cell RNA-seq dataset. Then we can use compareCluster to perform enrichment analysis for each cluster. Finally, we can use interpret to annotate cell types based on the enrichment results and marker genes.
## Not run in the book build: needs Seurat plus the pbmc3k 10x matrix
## (data/filtered_gene_bc_matrices/hg19) and the CellMarker xlsx file.
library(Seurat)
dir = "data/filtered_gene_bc_matrices/hg19"
pbmc.data <- Read10X(data.dir = dir)
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k",
min.cells=3, min.features=200)
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc,
subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5
)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize",
scale.factor = 10000)
pbmc <- ScaleData(pbmc)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst",
nfeatures = 2000)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
pbmc <- RunUMAP(pbmc, dims = 1:10)
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE)
library(dplyr)
topN_marker <- function(markers, n) {
markers %>%
group_by(cluster) %>%
dplyr::filter(avg_log2FC > 1) %>%
slice_head(n = n) %>%
ungroup()
}
top20 <- topN_marker(pbmc.markers, 20)
# downloaded from: http://www.bio-bigdata.center/CellMarker_download_files/file/Cell_marker_Human.xlsx
# note: `gene` from FindAllMarkers() is a gene symbol, and so is CellMarker's
# `marker` column -- both must be the same identifier type for the join to work.
cm <- rio::import("Cell_marker_Human.xlsx")
x <- compareCluster(gene~cluster, data=top20, fun=enricher, TERM2GENE=cm[,c("cell_name", "marker")])
y <- interpret(x, task="annotation")The output y is a list of interpretation results, one for each cluster. We can extract the inferred cell types.
> sapply(y, \(x) x$cell_type)
0
"Naive T Cell"
1
"Classical Monocyte"
2
"CD4+ T cell"
3
"Follicular B cell"
4
"CD8+ Cytotoxic T Cell"
5
"CD16+ monocyte (Non-classical monocyte)"
6
"Natural Killer (NK) cell"
7
"Plasmacytoid Dendritic Cell (pDC)"
8
"Megakaryocyte"This result is highly consistent with the manual annotation from the Seurat pbmc3k tutorial:
| Cluster ID | Markers | Cell Type |
|---|---|---|
| 0 | IL7R, CCR7 | Naive CD4+ T |
| 1 | CD14, LYZ | CD14+ Mono |
| 2 | IL7R, S100A4 | Memory CD4+ |
| 3 | MS4A1 | B |
| 4 | CD8A | CD8+ T |
| 5 | FCGR3A, MS4A7 | FCGR3A+ Mono |
| 6 | GNLY, NKG7 | NK |
| 7 | FCER1A, CST3 | DC |
| 8 | PPBP | Platelet |
The full report provides detailed reasoning, confidence levels, and supporting evidence (markers/pathways) for each cluster assignment, offering transparency and explainability that simple label transfer methods lack.
## Output of the annotation example above, which is not run in the book build.
print(y)Instead of relying on a single prompt, clusterProfiler employs a Multi-Agent System (MAS) to ensure accuracy and depth. This system consists of three specialized agents that work in a pipeline:
You can activate this deep mode using the interpret_agent() function.
## Kept off: the multi-agent pipeline runs three LLM calls per cluster (Cleaner,
## Detective, Storyteller) on top of the single-call examples above, which would
## push this chapter's cost per build well past its budget. To enable it once the
## build is known to be green, replace `eval: false` with `eval: !expr has_ai` and
## wrap the result in `interpret_failed()` as the examples above do.
# Provide biological context to help Agent Cleaner
context <- "scRNA-seq analysis of CD8+ T cells in Tumor Microenvironment, comparing Exhausted vs. Naive states."
res <- interpret_agent(edo, context = context)Enrichment analysis often treats genes as a “bag of words,” ignoring their interactions and expression changes. The Knowledge-Guided Interpretation mode injects external knowledge to empower the AI’s reasoning.
By setting add_ppi = TRUE, the system fetches protein-protein interaction data (from STRING). The AI can then identify functional modules (e.g., a TCR signaling complex) rather than just isolated genes.
By providing gene_fold_change, the AI can infer the direction of pathway activity. For example, if the “Apoptosis” pathway is enriched but pro-apoptotic genes are downregulated while anti-apoptotic genes are upregulated, the AI will correctly interpret this as “Apoptosis Resistance.”
# Prepare a named vector of fold changes
gene_list <- c("CD8A" = 2.5, "PDCD1" = 1.8, "GZMB" = 3.2)
res <- interpret(edo,
task = "cell_type",
add_ppi = TRUE, # Enable PPI network analysis
gene_fold_change = gene_list # Inject expression data
)For cell type annotation, LLMs can sometimes hallucinate. To prevent this, clusterProfiler supports Reference-Guided Interpretation.
You can provide “prior knowledge” (e.g., results from SingleR, scGPT, or manual rough annotation) to the AI. The AI acts as a validator and refiner: * Validation: Checks if pathway evidence supports the prior label. * Refinement: Refines a broad label (e.g., “T cell”) into a specific state (e.g., “Proliferating CD8+ T cell”) based on pathway activity. * Correction: Flags potential misannotations if the evidence contradicts the prior.
## Kept off: assigns to `res` and prints nothing, so running it would only add
## an API call per build without changing the rendered page.
# Prior knowledge from SingleR
my_priors <- c("Cluster1" = "T cells")
res <- interpret(edo, prior = my_priors, task = "cell_type")For complex datasets, interpret_hierarchical() mimics the human thought process of annotating major lineages first (e.g., Myeloid) and then subtypes (e.g., M1 Macrophage). It enforces lineage constraints to prevent impossible annotations (e.g., a T cell subtype appearing within a Myeloid cluster).
This approach is also highly applicable to Single-cell Trajectory Inference. Developmental processes inherently follow a hierarchical structure (e.g., Stem Cell -> Progenitor -> Terminally Differentiated Cell). By utilizing this hierarchical relationship, interpret_hierarchical() can provide context-aware interpretations that respect the biological differentiation path, ensuring that downstream states are interpreted within the context of their upstream progenitors.
# Mapping between minor and major clusters
cluster_mapping <- c(
"SubCluster1_1" = "MajorCluster1",
"SubCluster1_2" = "MajorCluster1",
"SubCluster2_1" = "MajorCluster2"
)
# Hierarchical interpretation
res_hier <- interpret_hierarchical(
x_minor = enrich_minor,
x_major = enrich_major,
mapping = cluster_mapping
)In real-world research, enrichment analysis sometimes fails to return significant pathways due to small gene sets or background noise.
clusterProfiler introduces a Gene-Based Fallback Mode. When no enriched pathways are found, the agent does not simply error out. Instead, it: 1. Directly analyzes the function of the input genes. 2. Retrieves PPI networks for these genes. 3. Infers biological function based on gene connectivity and function, providing a “Medium Confidence” report instead of an empty result.
# When enrichment fails, this still works
res <- interpret_agent(tough_genes, add_ppi = TRUE)Text reports are great, but graphical representations are often preferred for presentations. The plot() method for interpretation results uses ggtangle (a grammar of graphics for networks) to visualize the AI-inferred regulatory network.
The resulting plot highlights: * Key Drivers: Central nodes identified by the AI. * Activation (Green) / Inhibition (Red): Regulatory relationships inferred from data. * Interactions (Grey): Physical associations.
# Visualize the interpretation result
plot(res)This feature creates a closed loop: from Enrichment (Statistics) to Interpretation (Insight) and finally to Visualization (Communication).
The examples above are deliberately small. This section collects the guidance that matters once you point interpret() at a real analysis: how to write context, how to choose and switch models, what multi-model workflows are (and are not) supported, and how to keep the report honest.
context Be?context is the only place where the model learns what experiment produced the gene list, and it is the highest-leverage input you control. An enrichment table is ambiguous on its own: “extracellular matrix organization” is up-regulated in fibrosis, in wound healing and in tumour invasion alike. The model cannot tell those apart, so it averages over them and returns a bland report. A specific context removes the ambiguity.
A weak context names only the assay:
context <- "scRNA-seq of T cells"A strong context answers six questions in one or two sentences:
| Question | Example |
|---|---|
| Species and tissue | human breast tumour, bulk tissue |
| Disease / biological setting | breast carcinoma, tumour grade |
| Experimental design | 200 samples from the MAINZ cohort |
| Contrast tested, and its direction | grade III vs. grade I; positive values are higher in grade III |
| Cell type / composition | bulk tissue, no cell-type resolution |
| Hypothesis under test | progression involves loss of epithelial differentiation and ECM remodelling |
Written out, that is still only a few lines:
context <- paste(
"Human breast-cancer transcriptome (MAINZ cohort, bulk tumour tissue, 200",
"samples). Genes are ranked by the log2 ratio of geometric-mean expression in",
"grade III versus grade I tumours, so positive values are higher in the",
"higher-grade tumours. No cell-type resolution is available. Hypothesis:",
"progression from grade I to grade III involves loss of epithelial",
"differentiation and increased extracellular-matrix remodelling."
)Two rules of thumb follow. First, state the direction of the contrast – “A vs. B” without saying which side is up is a common way to get a report with the regulation backwards. Second, do not put the answer in the context: if you write “we expect apoptosis to be suppressed”, the model will often agree with you regardless of the data. Describe the design, not the conclusion you are hoping for. (If your data are from single cells, add the cell type or the compartment the cluster belongs to; for bulk tissue, saying that no cell-type resolution exists is more useful than omitting the point.)
Model selection is delegated to aisdk. interpret() itself has no model registry and no provider logic: when model is left at its default (NULL) it uses whatever aisdk currently considers the default, and model = "provider:model" overrides that default for a single call.
library(aisdk)
# Session-wide default, used by every interpret() call that omits `model`
set_model("deepseek:deepseek-v4-flash")
get_model()The setup chunk on this page calls set_model() once, so every example that omits model uses DeepSeek. To send one call elsewhere, pass model directly:
res <- interpret(edo, model = "openai:gpt-4o")Credentials follow the same rule. aisdk reads each provider’s key from an environment variable: OPENAI_API_KEY, ANTHROPIC_API_KEY and GEMINI_API_KEY for the three providers bundled with aisdk, and the corresponding variable for the adapters in aisdk.providers (for example DEEPSEEK_API_KEY, which is what this chapter uses). Put the key in a project .Renviron or export it in the shell; interpret() never takes a key argument. If you never call set_model(), aisdk falls back to its built-in default, openai:gpt-4o – which is why the examples you may have seen elsewhere appear to use GPT-4o by default.
aisdk ships a catalogue of the providers and models it knows about:
list_models() # every provider and model aisdk knows about
list_models("openai") # just one providerWhich one should you use? interpret() only needs a model that can follow the structured-output schema. In practice the useful split is between small, cheap models – fine for first-pass triage – and larger reasoning models, which tend to produce more coherent mechanistic narratives but cost more and consume more of the output budget (see max_tokens below). Whatever you choose, keep it fixed once you start comparing reports: different models write different reports, and that is the subject of the next subsection.
There is no built-in consensus or ensemble mode in interpret() or interpret_agent(). A call goes to exactly one model. interpret_agent() is a multi-agent pipeline, not a multi-model one: its three agents share a single model.
What you can do today is run the same input through several models and compare the reports yourself. Because model = is a per-call argument, this needs no extra machinery:
res_deepseek <- interpret(edo, model = "deepseek:deepseek-v4-flash")
res_gpt <- interpret(edo, model = "openai:gpt-4o")Each call is a separate, billable request, and the two reports will not be identical even for the same model (see below). A pragmatic middle ground is to use one model to write the report and a second, independent call only to criticise it: ask a different model whether every gene and pathway named in the narrative is actually present in the input table. That is a manual workflow, but it catches the most common failure mode – plausible-sounding genes that were never in the data. If what you need is provider failover rather than consensus, note that aisdk provides generate_text_fallback() for its own generation calls, but interpret() does not expose a fallback argument.
Non-determinism. interpret() samples at temperature = 0.3 by default, so two calls with identical input will differ in wording and occasionally in emphasis. This is expected, not a bug. It has a direct consequence for reproducible builds: Quarto’s chunk cache (cache: true) is what keeps the published page stable – the LLM call runs once and its output is reused on later renders. Clear the cache and you will get a different report. If you need the exact same text in two places, cache the object rather than calling twice.
Hallucination. The system prompt instructs the model to ground every claim in the supplied table and to cite supporting terms, and the structured output forces it to name genes in refined_network and network_evidence. That reduces hallucination; it does not eliminate it. Before you quote a report, check the two things that are cheap to check: every gene named in the narrative appears in the geneID (or, for GSEA, core_enrichment) column of the input table, and every pathway named appears in its Description column. If you set add_ppi = TRUE, remember that the edges the model cites come from a STRING query, not from your data.
n_pathways. Only the top n_pathways terms (default 20) are sent to the model, ordered by adjusted p-value (p.adjust, falling back to pvalue). Raising it gives the model more to work with but dilutes the signal and enlarges the prompt; lowering it forces a tighter focus. For a compareClusterResult, interpret() splits the result by cluster and returns one report per cluster, so the number of calls – and the cost – scales with the number of clusters, not with n_pathways.
max_tokens. The default is 8192. Reasoning models spend part of that budget deliberating before they emit the JSON, and if they run out, the reply is truncated and cannot be parsed; the symptom is a warning mentioning finish_reason: length and a placeholder report. If you hit it, raise max_tokens (the diagnostic suggests 16384 or more) or switch to a model that does not reason before answering.
Seeing the cost. interpret() returns only the parsed report; token usage is not part of the returned object. Run with verbose = TRUE to print the raw response and the prompt/completion/total token counts – the practical way to watch spend while you tune n_pathways and max_tokens.
add_ppi. Setting add_ppi = TRUE queries STRING for the genes in the top terms and injects the edge list into the prompt, which lets the model identify functional modules rather than isolated genes. It requires network access, adds latency, and is a second failure point; the query is wrapped so that a STRING outage degrades to “no network evidence” instead of an error. Leave it at its default (FALSE) unless module-level reasoning is what you need.
gene_fold_change. Passing a named numeric vector of log fold changes adds the direction of regulation to the prompt. Note the limitation: the helper intersects your vector with the genes in the result’s geneID column, and a gseaResult table has no geneID column (it carries core_enrichment instead). So gene_fold_change is used for over-representation results (enrichResult, compareClusterResult) and silently ignored for GSEA results – which is harmless, because a GSEA result already encodes direction in the sign of its enrichment score.
Ranking. For over-representation results the terms are ordered by p.adjust; for GSEA the direction the model sees is the sign of NES. If your GSEA ranked list was built so that positive means “up in the condition of interest”, say so in context – otherwise the model has to guess the orientation.
Validating a report before trusting it. Treat the output as a first draft from a competent but over-eager collaborator. The minimum checks are the two above (every gene and pathway exists in the input); beyond that, ask whether the narrative is falsifiable by your data – if it claims a mechanism, is there a term or gene in the table that could have contradicted it? A report that only restates the pathway names in longer sentences has added no information, and sharpening context or raising n_pathways is usually the fix.
This example runs end to end: it builds a GSEA result offline from the geneList shipped with DOSE, interprets it with a specific context, and prints the report. It is the interpret() counterpart of the interpret_agent() call that users often start from; interpret_agent() itself is kept disabled on this page because it makes several LLM calls per build.
library(clusterProfiler)
library(org.Hs.eg.db)
data(geneList, package = "DOSE")
# A real enrichment result, computed offline: GSEA over GO cellular components.
# geneList is the DOSE demo list -- log2 ratios of grade III vs. grade I tumours
# from the MAINZ breast-cancer cohort.
gse <- gseGO(
geneList = geneList,
OrgDb = org.Hs.eg.db,
keyType = "ENTREZID",
ont = "CC",
minGSSize = 15,
maxGSSize = 100,
eps = 0,
pvalueCutoff = 0.05,
pAdjustMethod = "BH"
)
# A context that names the species, tissue, design, contrast and hypothesis.
context <- paste(
"Human breast-cancer transcriptome (MAINZ cohort, bulk tumour tissue).",
"Genes are ranked by the log2 ratio of geometric-mean expression in grade",
"III versus grade I tumours, so positive values are higher in the",
"higher-grade tumours. Hypothesis: progression from grade I to grade III",
"involves loss of epithelial differentiation and increased",
"extracellular-matrix remodelling."
)
report <- interpret(gse, context = context, n_pathways = 15)
if (interpret_failed(report)) {
stop("interpret() returned no report -- check the API key and quota:\n",
report$overview)
}
report## Interpretation Result
### Cluster: Default
### 1. Overview
Grade III tumours show coordinated down-regulation of extracellular matrix and basement membrane structural modules (interstitial matrix NES -2.65; elastic fiber -2.46; basement membrane -2.16; collagen trimer -2.16) and up-regulation of mitotic and proteostatic modules (spindle midzone NES 2.47; spindle microtubule 2.21; outer kinetochore 2.25; kinesin complex 2.18; proteasome complex 2.17) plus immunoglobulin complex (NES 2.20). This supports a progression signature of ECM/epithelial differentiation loss with increased proliferation and immune/plasma-cell infiltration.
### 2. Regulatory Drivers
The enrichment results do not include upstream regulator or transcription-factor terms, so no specific TF can be directly assigned. The reciprocal enrichment patterns imply two regulatory axes: (i) repression of ECM/basement-membrane structural genes, inferred from down-regulation of interstitial matrix, basement membrane, collagen trimer, elastic fiber, and microfibril; and (ii) activation of mitotic gene expression, inferred from up-regulation of spindle midzone, spindle microtubule, outer kinetochore, and kinesin complex. Candidate master regulators of the mitotic axis include proliferation-associated TFs (e.g., FOXM1, MYC, E2F), and of the ECM-repression axis include EMT-inducing TFs (e.g., SNAI1/2, TWIST1, ZEB1/2), but these are hypotheses not directly measured in the provided enrichment. The immunoglobulin complex up-regulation may reflect immune-cell infiltration rather than tumour-cell-intrinsic TF activity.
### 3. Key Mechanisms
Major themes: (1) Loss of extracellular architecture: negative enrichment of interstitial matrix, solid phase, non-collagenous component, gel phase, elastic fiber, microfibril, basement membrane, collagen trimer, and complex of collagen trimers indicates reduced deposition/remodelling of collagenous and elastic matrices and basement membrane. (2) Increased mitotic machinery: positive enrichment of spindle midzone, spindle microtubule, outer kinetochore, and kinesin complex indicates enhanced mitotic spindle assembly, chromosome segregation, and kinetochore function. (3) Proteostasis: proteasome complex up-regulation suggests increased protein degradation, possibly supporting mitotic progression. (4) Immune component: immunoglobulin complex up-regulation suggests B-cell/plasma-cell infiltration or antibody production in the tumour microenvironment.
### 4. Crosstalk & Interactions
The mitotic module is internally connected: kinesin complex motors act on spindle microtubules and concentrate at the spindle midzone, while outer kinetochore provides microtubule attachment sites; proteasome complex may regulate these processes by degrading cell-cycle regulators. ECM loss may promote invasion and release of sequestered growth factors, potentially reinforcing proliferative signals. Immunoglobulin complex up-regulation could interact with ECM fragments and immune microenvironment, but direct evidence for ECM-immune crosstalk is not provided. The inverse relationship between ECM structural terms and mitotic/proteasome terms suggests a phenotypic switch from differentiated/ECM-rich to proliferative/ECM-poor state.
### 5. Hypothesis
Progression from grade I to grade III breast cancer involves a dual biological program: reduced expression of interstitial matrix, elastic fiber, microfibril, basement membrane, and collagen trimers, consistent with loss of epithelial/stromal differentiation and increased invasive capacity; and increased expression of spindle midzone, spindle microtubule, outer kinetochore, kinesin complex, and proteasome complex, consistent with enhanced mitotic proliferation and chromosome segregation. The concurrent rise in immunoglobulin complex suggests an immune/plasma-cell component. Thus, higher tumour grade is characterized by ECM breakdown and proliferative/immune activation, supporting the original hypothesis of epithelial differentiation loss and ECM remodelling.
### 6. Narrative Draft
In this human breast-cancer transcriptome cohort, comparison of grade III versus grade I tumours revealed a coordinated down-regulation of extracellular matrix and basement membrane structural components, including interstitial matrix, elastic fiber, microfibril, basement membrane, collagen trimer, and complex of collagen trimers (all NES < 0). In contrast, mitotic and proteostatic modules were up-regulated, including spindle midzone, spindle microtubule, outer kinetochore, kinesin complex, and proteasome complex (all NES > 0). An immunoglobulin complex signature was also increased. These findings suggest that progression to grade III involves loss of ECM/epithelial differentiation, with reduced collagenous and elastic matrices and basement membrane, alongside increased mitotic spindle and kinetochore activity that drives proliferation and genomic instability. The proteasome up-regulation may support mitotic progression through regulated protein degradation, while the immunoglobulin signature may reflect an altered immune microenvironment. Together, these modules define a high-grade phenotype characterized by ECM remodelling and proliferative activation, consistent with the hypothesis that grade progression entails loss of epithelial differentiation and increased extracellular-matrix remodelling.
### 7. Refined Regulatory Network
Key Interactions:
COL1A1/COL1A2/COL3A1 -- Interstitial matrix (structural component) - GO:0005581 collagen trimer and GO:0005614 interstitial matrix are both negatively enriched (NES<0), indicating collagen trimers are structural constituents of interstitial matrix lost in grade III.
ELN -- Elastic fiber (structural component) - GO:0071953 elastic fiber is negatively enriched (NES -2.46), consistent with reduced elastin-based elastic fiber assembly.
FBN1 -- Microfibril (structural component) - GO:0001527 microfibril is negatively enriched (NES -2.27), consistent with reduced fibrillin microfibril formation.
LAMA1/LAMB1/LAMC1 -- Basement membrane (structural component) - GO:0005604 basement membrane is negatively enriched (NES -2.16), indicating loss of laminin-rich basement membrane.
KIF11/KIF15 -- Spindle midzone (motor activity) - GO:0005871 kinesin complex and GO:0051233 spindle midzone are both positively enriched (NES>0), linking kinesin motors to midzone function.
KIF11/KIF15 -- Spindle microtubule (motor activity) - GO:0005871 kinesin complex and GO:0005876 spindle microtubule are positively enriched, indicating kinesin-driven microtubule organization.
NDC80 complex -- Outer kinetochore (binding) - GO:0000940 outer kinetochore is positively enriched (NES 2.25), and the NDC80 complex is a core outer kinetochore component.
Proteasome complex -- Mitotic regulators (degradation) - GO:0000502 proteasome complex is positively enriched (NES 2.17) and may degrade cell-cycle regulators to promote mitotic progression.
Immunoglobulin complex -- Tumour microenvironment (immune signaling) - GO:0019814 immunoglobulin complex is positively enriched (NES 2.20), suggesting B-cell/plasma-cell infiltration or antibody production.
Interstitial matrix -- Tissue integrity (structural support) - GO:0005614 interstitial matrix is negatively enriched (NES -2.65), and its loss may reduce tissue integrity and promote invasion.
**Network Evidence:**
The refined network is supported by four enriched modules: (i) an ECM/basement-membrane module comprising interstitial matrix (GO:0005614), solid phase of interstitial matrix (GO:0140151), non-collagenous component of interstitial matrix (GO:0140149), gel phase of interstitial matrix (GO:0140150), elastic fiber (GO:0071953), microfibril (GO:0001527), basement membrane (GO:0005604), collagen trimer (GO:0005581), and complex of collagen trimers (GO:0098644), all negatively enriched; (ii) a mitotic spindle/kinetochore module comprising spindle midzone (GO:0051233), spindle microtubule (GO:0005876), outer kinetochore (GO:0000940), and kinesin complex (GO:0005871), all positively enriched; (iii) a proteasome module (GO:0000502) positively enriched; and (iv) an immunoglobulin complex module (GO:0019814) positively enriched. These modules jointly support a high-grade phenotype of ECM loss and proliferative/immune activation. The printed report follows the schema described at the top of this chapter. Read it in this order:
context was probably too vague.plot(report) draws the same network with ggtangle when it is available.