Statistical enrichment results are evidence to inspect, not biological conclusions by themselves. This chapter describes how to move from enriched terms, core genes, networks, expression trends, and prior knowledge to a transparent interpretation. AI assistance is presented as an optional evidence-synthesis layer after those checks, not as a replacement for the enrichment analysis or literature review.
Enrichment results identify terms and genes associated with the input, but they do not by themselves establish a mechanism. Interpretation therefore starts with the result object, core genes, background, database provenance, and any supporting expression or network evidence.
interpret functioninterpret() accepts enrichResult, gseaResult, or compareClusterResult objects and asks an LLM to organize the supplied evidence into a structured report. The report is a draft for evidence review, not a replacement for the enrichment test or literature validation.
The function requires an API key for a provider supported by aisdk. If model is not supplied, it uses the current model from aisdk::get_model() / aisdk::set_model(); an explicit model argument overrides that default for one call. The examples below use DeepSeek through DEEPSEEK_API_KEY, but other aisdk providers can be used.
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)
}
cat(render_interpretation_blocks(res), sep = "\n")Overview:
Enrichment analysis reveals a dual proliferative-inflammatory signature. Cell cycle and DNA replication are highlighted by Meier-Gorlin syndrome (DOID:0060306; ORC1, ORC6, CDT1, CDC6, MCM5, CDC45) and multiple carcinoma terms (stomach cancer DOID:10534, breast cancer DOID:1612, esophageal cancer DOID:5041, bladder cancer DOID:11054, in situ carcinoma DOID:8719, female reproductive organ cancer DOID:120) with top marker genes including CDC45, MCM10, CDC20, CDCA8, KIF23, CENPE, NDC80, TOP2A, NCAPH, CCNB2, MELK, FOXM1, and MYBL2. Inflammatory and fibrotic processes are represented by myocarditis (DOID:820), cardiomyopathy (DOID:0050700), scleroderma (DOID:419), rheumatic disease (DOID:1575), and bronchiolitis obliterans (DOID:2799), with S100A8 (6279), S100A9 (6280), and MMP1 (4312). Viral infectious disease (DOID:934) and type 2 diabetes (DOID:9352) further indicate immune-metabolic crosstalk.
Regulatory drivers:
Key mechanisms:
Crosstalk and interactions:
Shared genes across terms reveal a network connecting proliferation, inflammation, and tissue remodeling. MMP1 (4312) is present in stomach cancer (DOID:10534), breast cancer (DOID:1612), esophageal cancer (DOID:5041), cardiomyopathy (DOID:0050700), scleroderma (DOID:419), and type 2 diabetes (DOID:9352). S100A8 (6279) and S100A9 (6280) co-occur in myocarditis (DOID:820) and cardiomyopathy (DOID:0050700). TOP2A (7153) is shared between breast cancer (DOID:1612), bladder cancer (DOID:11054), and scleroderma (DOID:419). FOXM1 (2305) links stomach cancer (DOID:10534) and type 2 diabetes (DOID:9352) with MMP1 (4312) and NDC80 (10403). These overlaps suggest that inflammatory mediators (S100A8/9, MMP1) may cooperate with proliferative drivers (FOXM1, TOP2A) to exacerbate both neoplastic and fibrotic pathologies.
Hypothesis:
We hypothesize that the enriched signature reflects a dual proliferative-inflammatory program. DNA replication licensing and mitotic regulators (FOXM1, MYBL2, CDC20, TOP2A, etc.) drive carcinoma progression, while S100A8/S100A9 and MMP1 promote an inflammatory and matrix-degrading microenvironment. This combined program contributes to the observed enrichment in cancers, cardiomyopathies, scleroderma, and type 2 diabetes, suggesting that targeting both proliferation and inflammation may be therapeutically beneficial.
Narrative draft:
The enrichment results define a molecular landscape in which cell cycle and DNA replication are strongly activated, as evidenced by Meier-Gorlin syndrome (DOID:0060306) and multiple carcinoma terms (stomach DOID:10534, breast DOID:1612, esophageal DOID:5041, bladder DOID:11054, in situ carcinoma DOID:8719, female reproductive organ cancer DOID:120). Top marker genes include DNA replication factors (CDC45, MCM10) and mitotic regulators (CDC20, CDCA8, KIF23, CENPE, NDC80, TOP2A, NCAPH, CCNB2, MELK, FOXM1, MYBL2). In parallel, inflammatory and fibrotic diseases such as myocarditis (DOID:820), cardiomyopathy (DOID:0050700), scleroderma (DOID:419), rheumatic disease (DOID:1575), and bronchiolitis obliterans (DOID:2799) are enriched, with shared inflammatory mediators S100A8 (6279), S100A9 (6280), and MMP1 (4312). Viral infectious disease (DOID:934) and type 2 diabetes (DOID:9352) further connect immune and metabolic pathways. The co-occurrence of MMP1 in cancers and fibrotic diseases, and the link between FOXM1 and MMP1 in stomach cancer and type 2 diabetes, suggest a crosstalk between proliferative and inflammatory programs. We propose that a FOXM1/MYBL2-driven mitotic program and an S100A8/S100A9-MMP1 inflammatory axis cooperate to drive the observed disease enrichment, offering a dual-targeting opportunity.
Network evidence:
The refined network centers on two modules: (1) a DNA replication licensing complex (ORC1, ORC6, CDC6, CDT1, MCM5, CDC45) supported by Meier-Gorlin syndrome (DOID:0060306), and (2) an inflammatory/fibrotic module (S100A8, S100A9, MMP1, MMP9, TOP2A, FOXM1, NDC80) supported by myocarditis (DOID:820), cardiomyopathy (DOID:0050700), scleroderma (DOID:419), stomach cancer (DOID:10534), bladder cancer (DOID:11054), and type 2 diabetes (DOID:9352). FOXM1 connects the proliferative and inflammatory modules by co-occurring with MMP1 and NDC80 in stomach cancer (DOID:10534) and with MMP1 in type 2 diabetes (DOID:9352).
## Interpretation Result
35.1.1 Cluster: Default
35.1.2 1. Overview
Enrichment analysis reveals a dual proliferative-inflammatory signature. Cell cycle and DNA replication are highlighted by Meier-Gorlin syndrome (DOID:0060306; ORC1, ORC6, CDT1, CDC6, MCM5, CDC45) and multiple carcinoma terms (stomach cancer DOID:10534, breast cancer DOID:1612, esophageal cancer DOID:5041, bladder cancer DOID:11054, in situ carcinoma DOID:8719, female reproductive organ cancer DOID:120) with top marker genes including CDC45, MCM10, CDC20, CDCA8, KIF23, CENPE, NDC80, TOP2A, NCAPH, CCNB2, MELK, FOXM1, and MYBL2. Inflammatory and fibrotic processes are represented by myocarditis (DOID:820), cardiomyopathy (DOID:0050700), scleroderma (DOID:419), rheumatic disease (DOID:1575), and bronchiolitis obliterans (DOID:2799), with S100A8 (6279), S100A9 (6280), and MMP1 (4312). Viral infectious disease (DOID:934) and type 2 diabetes (DOID:9352) further indicate immune-metabolic crosstalk.
35.1.3 2. Regulatory Drivers
FOXM1 (2305) and MYBL2 (4605) emerge as master transcriptional regulators of the mitotic gene program, as they are top marker genes and FOXM1 is co-enriched in stomach cancer (DOID:10534) and type 2 diabetes (DOID:9352). DNA replication licensing is driven by the ORC1/ORC6-CDC6-CDT1-MCM5-CDC45 module, directly supported by their co-enrichment in Meier-Gorlin syndrome (DOID:0060306). S100A8 (6279) and S100A9 (6280) act as inflammatory amplifiers, co-enriched in myocarditis (DOID:820) and cardiomyopathy (DOID:0050700).
35.1.4 3. Key Mechanisms
- Cell cycle and DNA replication: Meier-Gorlin syndrome (DOID:0060306) and carcinoma terms (stomach DOID:10534, breast DOID:1612, esophageal DOID:5041, bladder DOID:11054) are enriched with replication licensing (CDC45 8318, MCM5 4174, CDT1 81620, CDC6 990, ORC1 4998, ORC6 23594) and mitotic markers (CDC20 991, CDCA8 55143, KIF23 9493, CENPE 1062, NDC80 10403, TOP2A 7153, NCAPH 23397, CCNB2 9133, MELK 9833). 2) Inflammation and fibrosis: myocarditis (DOID:820), cardiomyopathy (DOID:0050700), scleroderma (DOID:419), rheumatic disease (DOID:1575), and bronchiolitis obliterans (DOID:2799) involve S100A8 (6279), S100A9 (6280), MMP1 (4312), and MMP9 (4318). 3) Immune-metabolic crosstalk: viral infectious disease (DOID:934) and type 2 diabetes (DOID:9352) share MMP1 (4312) and FOXM1 (2305).
35.1.5 4. Crosstalk & Interactions
Shared genes across terms reveal a network connecting proliferation, inflammation, and tissue remodeling. MMP1 (4312) is present in stomach cancer (DOID:10534), breast cancer (DOID:1612), esophageal cancer (DOID:5041), cardiomyopathy (DOID:0050700), scleroderma (DOID:419), and type 2 diabetes (DOID:9352). S100A8 (6279) and S100A9 (6280) co-occur in myocarditis (DOID:820) and cardiomyopathy (DOID:0050700). TOP2A (7153) is shared between breast cancer (DOID:1612), bladder cancer (DOID:11054), and scleroderma (DOID:419). FOXM1 (2305) links stomach cancer (DOID:10534) and type 2 diabetes (DOID:9352) with MMP1 (4312) and NDC80 (10403). These overlaps suggest that inflammatory mediators (S100A8/9, MMP1) may cooperate with proliferative drivers (FOXM1, TOP2A) to exacerbate both neoplastic and fibrotic pathologies.
35.1.6 5. Hypothesis
We hypothesize that the enriched signature reflects a dual proliferative-inflammatory program. DNA replication licensing and mitotic regulators (FOXM1, MYBL2, CDC20, TOP2A, etc.) drive carcinoma progression, while S100A8/S100A9 and MMP1 promote an inflammatory and matrix-degrading microenvironment. This combined program contributes to the observed enrichment in cancers, cardiomyopathies, scleroderma, and type 2 diabetes, suggesting that targeting both proliferation and inflammation may be therapeutically beneficial.
35.1.7 6. Narrative Draft
The enrichment results define a molecular landscape in which cell cycle and DNA replication are strongly activated, as evidenced by Meier-Gorlin syndrome (DOID:0060306) and multiple carcinoma terms (stomach DOID:10534, breast DOID:1612, esophageal DOID:5041, bladder DOID:11054, in situ carcinoma DOID:8719, female reproductive organ cancer DOID:120). Top marker genes include DNA replication factors (CDC45, MCM10) and mitotic regulators (CDC20, CDCA8, KIF23, CENPE, NDC80, TOP2A, NCAPH, CCNB2, MELK, FOXM1, MYBL2). In parallel, inflammatory and fibrotic diseases such as myocarditis (DOID:820), cardiomyopathy (DOID:0050700), scleroderma (DOID:419), rheumatic disease (DOID:1575), and bronchiolitis obliterans (DOID:2799) are enriched, with shared inflammatory mediators S100A8 (6279), S100A9 (6280), and MMP1 (4312). Viral infectious disease (DOID:934) and type 2 diabetes (DOID:9352) further connect immune and metabolic pathways. The co-occurrence of MMP1 in cancers and fibrotic diseases, and the link between FOXM1 and MMP1 in stomach cancer and type 2 diabetes, suggest a crosstalk between proliferative and inflammatory programs. We propose that a FOXM1/MYBL2-driven mitotic program and an S100A8/S100A9-MMP1 inflammatory axis cooperate to drive the observed disease enrichment, offering a dual-targeting opportunity.
35.1.8 7. Refined Regulatory Network
Key Interactions: ORC1 – CDC6 (binding) - Both enriched in Meier-Gorlin syndrome (DOID:0060306); ORC1 recruits CDC6 to replication origins. CDC6 – CDT1 (binding) - Both enriched in Meier-Gorlin syndrome (DOID:0060306); CDC6 and CDT1 cooperate to load MCM. CDT1 – MCM5 (binding) - Both enriched in Meier-Gorlin syndrome (DOID:0060306); CDT1 loads MCM5 onto chromatin. MCM5 – CDC45 (binding) - Both enriched in Meier-Gorlin syndrome (DOID:0060306); form CMG helicase. ORC1 – ORC6 (binding) - Both enriched in Meier-Gorlin syndrome (DOID:0060306); ORC complex. S100A8 – S100A9 (binding) - Both enriched in myocarditis (DOID:820) and cardiomyopathy (DOID:0050700); form heterodimer. S100A9 – MMP1 (activation) - Both enriched in cardiomyopathy (DOID:0050700); S100A9 induces MMP1. MMP1 – MMP9 (activation) - Both enriched in scleroderma (DOID:419); co-enriched matrix metalloproteinases. MMP1 – FOXM1 (activation) - Both enriched in stomach cancer (DOID:10534) and type 2 diabetes (DOID:9352); FOXM1 regulates MMP1. FOXM1 – NDC80 (activation) - Both enriched in stomach cancer (DOID:10534); FOXM1 drives mitotic genes including NDC80. MMP1 – TOP2A (co-expression) - Both enriched in bladder cancer (DOID:11054) and scleroderma (DOID:419).
Network Evidence: The refined network centers on two modules: (1) a DNA replication licensing complex (ORC1, ORC6, CDC6, CDT1, MCM5, CDC45) supported by Meier-Gorlin syndrome (DOID:0060306), and (2) an inflammatory/fibrotic module (S100A8, S100A9, MMP1, MMP9, TOP2A, FOXM1, NDC80) supported by myocarditis (DOID:820), cardiomyopathy (DOID:0050700), scleroderma (DOID:419), stomach cancer (DOID:10534), bladder cancer (DOID:11054), and type 2 diabetes (DOID:9352). FOXM1 connects the proliferative and inflammatory modules by co-occurring with MMP1 and NDC80 in stomach cancer (DOID:10534) and with MMP1 in type 2 diabetes (DOID:9352).
interpret() supports three tasks:
task = "interpretation": organize enrichment evidence into a biological narrative;task = "annotation": combine marker genes and enrichment terms for cell-type annotation;task = "phenotyping": propose a state or phenotype label for a group.The input can be one result object, a list of compatible objects, or a compareCluster result processed by cluster. For example, list(go_res, kegg_res) combines compatible knowledge sources, while a compareCluster result returns one report per cluster. The returned fields are structured so that reports can be checked, compared, and rendered consistently.
For single-cell data, the marker table can come from Seurat::FindAllMarkers() or another marker-calling method. Keep the cluster and gene columns, use the same identifier namespace in the marker and CellMarker tables, and pass the grouped markers through compareCluster() before annotation. The executable example below uses the repository-local CellMarker table so that it does not depend on a separate count matrix or a pasted annotation report. A complete Seurat marker workflow is available in the pbmc3k tutorial.
The cell-type example reads the local datasets/Cell_marker_Human.xlsx marker table, builds four marker profiles, runs compareCluster() with enricher(), and passes the result to interpret(task = "annotation") when DEEPSEEK_API_KEY is available.
library(clusterProfiler)
library(readxl)
library(dplyr)
cell_marker_path <- file.path("datasets", "Cell_marker_Human.xlsx")
if (!file.exists(cell_marker_path)) {
stop("The cell-marker dataset is missing: ", cell_marker_path)
}
cell_marker <- readxl::read_excel(cell_marker_path)
TERM2GENE_cell <- cell_marker |>
transmute(
cell_name = .data[["Cell name"]],
marker = toupper(.data[["Marker"]])
) |>
filter(
!is.na(cell_name), !is.na(marker),
nzchar(cell_name), nzchar(marker)
) |>
distinct()
# A compact marker panel keeps the example transparent and the prompt small.
marker_sets <- list(
T_cell = c("CD3D", "CD3E", "TRBC1", "IL7R", "LTB", "MAL", "CCR7", "LST1"),
B_cell = c("CD79A", "MS4A1", "CD37", "CD79B", "HLA-DRA", "CD74", "BANK1"),
Monocyte = c("LYZ", "S100A8", "S100A9", "CTSS", "FCN1", "LGALS3", "TYROBP", "LST1"),
NK_cell = c("NKG7", "GNLY", "PRF1", "GZMB", "KLRD1", "FCGR3A", "XCL1", "XCL2")
)
marker_table <- do.call(
rbind,
lapply(names(marker_sets), function(cluster) {
data.frame(
cluster = cluster,
gene = marker_sets[[cluster]],
stringsAsFactors = FALSE
)
})
)
cell_marker_enrichment <- compareCluster(
gene ~ cluster,
data = marker_table,
fun = enricher,
TERM2GENE = TERM2GENE_cell,
pvalueCutoff = 1,
qvalueCutoff = 1,
minGSSize = 2,
maxGSSize = 500
)
annotation_reports <- interpret(
cell_marker_enrichment,
task = "annotation",
context = paste(
"Human single-cell marker panels. Each cluster is a compact set of",
"top marker genes and was tested against the local CellMarker database.",
"Use the enrichment terms as supporting evidence, but prefer markers",
"that distinguish closely related cell types. Return a concise label",
"and explain which supplied markers and cell-marker terms support it."
),
n_pathways = 12,
max_tokens = 4096,
temperature = 0.1
)
if (inherits(annotation_reports, "interpretation")) {
annotation_reports <- list(Default = annotation_reports)
}
if (!is.list(annotation_reports) || length(annotation_reports) == 0) {
stop("DeepSeek returned no cell-type annotation reports.")
}
for (report in annotation_reports) {
if (annotation_failed(report)) {
stop(
"DeepSeek cell-type annotation failed for cluster ",
if (!is.null(report$cluster)) report$cluster else "(unknown)",
": ",
if (!is.null(report$overview)) report$overview else "no diagnostic returned"
)
}
}The result below is generated from the live object returned by DeepSeek. It is rendered as an HTML block, not Markdown headings: the preview is readable in the chapter, while the raw print() output remains available in a collapsed section without entering the book’s table of contents.
cat(render_interpretation_blocks(annotation_reports), sep = "\n")Cell type:
Follicular B cell
Confidence:
Medium
Reasoning:
The cluster expresses core B-cell markers CD79A, CD79B, MS4A1, CD37, HLA-DRA, CD74, and BANK1. Among the top enriched terms, Follicular B cell has the highest gene coverage (5/7: CD74/CD79A/MS4A1/BANK1/CD79B) and the same adjusted p-value as Mature B cell. The presence of BANK1 and CD79B supports follicular B-cell identity, while the absence of plasma-cell markers (e.g., SDC1, XBP1, MZB1) and germinal-center markers (e.g., BCL6, AICDA) argues against those alternatives. Shared B-cell markers (CD79A, MS4A1, CD37, HLA-DRA, CD74) also match Mature B cell and Naive B cell terms, but the combination of BANK1 with MS4A1 and CD79B favors Follicular B cell.
Supporting markers/pathways:
## Cell Type Annotation
### Cluster: B_cell
**Cell Type:** Follicular B cell
**Confidence:** Medium
**Reasoning:**
The cluster expresses core B-cell markers CD79A, CD79B, MS4A1, CD37, HLA-DRA, CD74, and BANK1. Among the top enriched terms, Follicular B cell has the highest gene coverage (5/7: CD74/CD79A/MS4A1/BANK1/CD79B) and the same adjusted p-value as Mature B cell. The presence of BANK1 and CD79B supports follicular B-cell identity, while the absence of plasma-cell markers (e.g., SDC1, XBP1, MZB1) and germinal-center markers (e.g., BCL6, AICDA) argues against those alternatives. Shared B-cell markers (CD79A, MS4A1, CD37, HLA-DRA, CD74) also match Mature B cell and Naive B cell terms, but the combination of BANK1 with MS4A1 and CD79B favors Follicular B cell.
**Supporting Markers/Pathways:**
- CD79A
- MS4A1
- CD37
- CD79B
- HLA-DRA
- CD74
- BANK1
- Follicular B cell (CD74/CD79A/MS4A1/BANK1/CD79B)
- Mature B cell (CD37/HLA-DRA/CD74/CD79A)
- Naive B cell (CD37/CD74/CD79A/MS4A1)
Cell type:
Classical monocyte
Confidence:
High
Reasoning:
Enrichment terms include broad Myeloid cell (7/8) and Macrophage (7/8), but the most specific monocyte terms are Classical monocyte (4/8: FCN1/S100A8/S100A9/LYZ) and CD14 monocyte (3/8: FCN1/S100A8/LYZ). The top markers are canonical classical monocyte genes (FCN1, S100A8, S100A9, LYZ, CTSS) and lack macrophage-specific markers (CD68, CD163, MRC1) and DC-specific markers (CD1C, FCER1A, CLEC10A). Therefore, the cluster is Classical monocyte.
Supporting markers/pathways:
## Cell Type Annotation
### Cluster: Monocyte
**Cell Type:** Classical monocyte
**Confidence:** High
**Reasoning:**
Enrichment terms include broad Myeloid cell (7/8) and Macrophage (7/8), but the most specific monocyte terms are Classical monocyte (4/8: FCN1/S100A8/S100A9/LYZ) and CD14 monocyte (3/8: FCN1/S100A8/LYZ). The top markers are canonical classical monocyte genes (FCN1, S100A8, S100A9, LYZ, CTSS) and lack macrophage-specific markers (CD68, CD163, MRC1) and DC-specific markers (CD1C, FCER1A, CLEC10A). Therefore, the cluster is Classical monocyte.
**Supporting Markers/Pathways:**
- LYZ
- S100A8
- S100A9
- CTSS
- FCN1
- LGALS3
- TYROBP
- LST1
Cell type:
Natural killer cell
Confidence:
High
Reasoning:
The top enriched term is Natural killer cell with 8/8 supplied markers (FCGR3A, PRF1, XCL2, GZMB, NKG7, GNLY, XCL1, KLRD1) and the strongest p.adjust (1.03e-12). Cytotoxic cell (5/8), CD8+ T cell (6/8), and NKT cell (4/8) are lower-ranked alternatives. The supplied marker list contains NK-associated FCGR3A (CD16), KLRD1 (CD94), and XCL1/XCL2, plus cytotoxic effector genes NKG7, GNLY, PRF1, GZMB. It lacks canonical CD8+ T-cell markers such as CD3/CD8A/CD8B in the input, so the CD8+ T-cell and cytotoxic T-cell terms are explained by shared cytotoxic genes rather than by T-cell receptor identity. Therefore, the cluster is best assigned as Natural killer cell.
Supporting markers/pathways:
## Cell Type Annotation
### Cluster: NK_cell
**Cell Type:** Natural killer cell
**Confidence:** High
**Reasoning:**
The top enriched term is Natural killer cell with 8/8 supplied markers (FCGR3A, PRF1, XCL2, GZMB, NKG7, GNLY, XCL1, KLRD1) and the strongest p.adjust (1.03e-12). Cytotoxic cell (5/8), CD8+ T cell (6/8), and NKT cell (4/8) are lower-ranked alternatives. The supplied marker list contains NK-associated FCGR3A (CD16), KLRD1 (CD94), and XCL1/XCL2, plus cytotoxic effector genes NKG7, GNLY, PRF1, GZMB. It lacks canonical CD8+ T-cell markers such as CD3/CD8A/CD8B in the input, so the CD8+ T-cell and cytotoxic T-cell terms are explained by shared cytotoxic genes rather than by T-cell receptor identity. Therefore, the cluster is best assigned as Natural killer cell.
**Supporting Markers/Pathways:**
- FCGR3A
- PRF1
- XCL2
- GZMB
- NKG7
- GNLY
- XCL1
- KLRD1
- Natural killer cell
- Cytotoxic cell
Cell type:
Naive CD8+ T cell
Confidence:
Medium
Reasoning:
The cluster expresses pan-T cell markers CD3D, CD3E, and TRBC1, indicating T cell lineage. Enrichment top terms include Memory CD8+ T cell (5/8, p.adj=1.8e-10), Naive CD8+ T cell (6/8, p.adj=5.7e-10), and T cell (7/8, p.adj=5.7e-10). The presence of CCR7, IL7R, MAL, and LTB supports a naive T cell state, as MAL and CCR7 are naive T cell markers and MAL is absent from the Memory CD8+ T cell term. Among naive subsets, Naive CD8+ T cell includes the same six markers (CD3D, LTB, MAL, CD3E, CCR7, IL7R) and has a higher gene ratio than Naive CD4 T cell, which lacks LTB. However, CD8A/CD8B are not provided, so CD4 vs CD8 discrimination is uncertain; therefore confidence is Medium.
Supporting markers/pathways:
## Cell Type Annotation
### Cluster: T_cell
**Cell Type:** Naive CD8+ T cell
**Confidence:** Medium
**Reasoning:**
The cluster expresses pan-T cell markers CD3D, CD3E, and TRBC1, indicating T cell lineage. Enrichment top terms include Memory CD8+ T cell (5/8, p.adj=1.8e-10), Naive CD8+ T cell (6/8, p.adj=5.7e-10), and T cell (7/8, p.adj=5.7e-10). The presence of CCR7, IL7R, MAL, and LTB supports a naive T cell state, as MAL and CCR7 are naive T cell markers and MAL is absent from the Memory CD8+ T cell term. Among naive subsets, Naive CD8+ T cell includes the same six markers (CD3D, LTB, MAL, CD3E, CCR7, IL7R) and has a higher gene ratio than Naive CD4 T cell, which lacks LTB. However, CD8A/CD8B are not provided, so CD4 vs CD8 discrimination is uncertain; therefore confidence is Medium.
**Supporting Markers/Pathways:**
- CD3D
- CD3E
- TRBC1
- IL7R
- LTB
- MAL
- CCR7
When the book is rendered without DEEPSEEK_API_KEY, the two chunks above are skipped and the rest of the chapter still renders. When the key is present, a render performs a real DeepSeek call; keep the cache enabled for ordinary builds and clear interpretation_cache/ deliberately when a fresh report is wanted.
interpret_agent() applies three sequential stages to an interpretation request:
The implementation refers to these stages as Cleaner, Detective, and Storyteller. They are components of one shared-model pipeline, not separate statistical analyses. The stages share one evidence bundle, and their output remains a model-generated draft that should be checked against the original result object.
context <- "scRNA-seq analysis of CD8+ T cells in tumor microenvironment, comparing exhausted and naive states."
res <- interpret_agent(edo, context = context)Set add_ppi = TRUE to request protein interactions from STRING. Pass a named gene_fold_change vector to retain expression direction in the interpretation. These inputs add network and expression evidence; they do not establish pathway activity or causality.
# Prepare a named vector of fold changes
gene_list <- c("CD8A" = 2.5, "PDCD1" = 1.8, "GZMB" = 3.2)
res <- interpret(
edo,
task = "interpretation",
add_ppi = TRUE,
gene_fold_change = gene_list
)For cell-type annotation, prior supplies an existing label for comparison. The returned annotation can agree with, refine, or disagree with that label; inspect the markers and enrichment terms before accepting a change.
my_priors <- c("Cluster1" = "T cells")
res <- interpret(edo, prior = my_priors, task = "annotation")interpret_hierarchical() accepts major- and minor-cluster enrichment objects together with a mapping. It constrains each minor-cluster interpretation to the supplied major-cluster context; the mapping itself must be checked against the experimental design.
cluster_mapping <- c(
"SubCluster1_1" = "MajorCluster1",
"SubCluster1_2" = "MajorCluster1",
"SubCluster2_1" = "MajorCluster2"
)
res_hier <- interpret_hierarchical(
x_minor = enrich_minor,
x_major = enrich_major,
mapping = cluster_mapping
)When a result has no enriched terms, interpret_agent() can be given the input genes and optional PPI data. The implementation warns that no significant pathways were found, passes the raw genes to the Detective and Storyteller stages, and marks the returned object with data_source = "gene_list_only". The output should therefore be reported as a gene-based hypothesis; its confidence is taken from the returned object rather than assigned a fixed value.
res <- interpret_agent(tough_genes, add_ppi = TRUE)The plot() method displays relationships stored in an interpretation object. When the returned object contains them, the plot can distinguish key drivers, activation or inhibition labels, and physical interactions. These labels are model outputs and require review against the supplied expression and interaction evidence.
plot(res)The complete hand-off is:
enrichment result → optional fold-change/PPI/prior evidence
→ interpretation object → report or plotThe final report and plot are views of the interpretation object. They do not constitute a second enrichment test, and interaction edges remain external supporting evidence.
context should identify the experiment, biological setting, contrast direction, and any relevant cell-type or tissue resolution. interpret() cannot infer these details from pathway names alone.
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()Calls that omit model use the session-wide setting selected with set_model(). Pass model for a one-off provider or model.
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 providerModel choice affects latency, cost, context capacity, and output variability. Fix the model when comparing reports, because different models can produce different summaries from the same input.
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. For provider failover, 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 and reuse it.
Unsupported entities. The model is instructed to ground claims in the supplied table, but it can still produce incorrect output. Before quoting a report, check that every gene named in the narrative appears in the geneID (or, for GSEA, core_enrichment) column and every pathway named appears in Description. If add_ppi = TRUE, identify the cited edges as STRING evidence, not measurements from the experiment.
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 from individual 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.
Treat the output as an unverified draft. Check every gene and pathway named in the narrative against the input table; then assess whether each mechanistic statement has supporting evidence.
This example builds a GSEA result offline from the geneList shipped with DOSE, passes a defined context to interpret(), and renders the returned report. The interpret_agent() example above is disabled because it makes multiple model calls.
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)
}
cat(render_interpretation_blocks(report), sep = "\n")Overview:
Grade III versus grade I breast tumours show coordinated down-regulation of extracellular matrix and basement-membrane components (interstitial matrix GO:0005614; basement membrane GO:0005604; elastic fiber GO:0071953; collagen trimers GO:0098644/GO:0005581) and up-regulation of mitotic spindle/kinetochore/kinesin modules (spindle midzone GO:0051233; spindle microtubule GO:0005876; outer kinetochore GO:0000940; kinesin complex GO:0005871) plus proteasome and immunoglobulin complexes (GO:0000502; GO:0019814).
Regulatory drivers:
Key mechanisms:
Two major themes emerge. First, loss of ECM/basement-membrane architecture in grade III is supported by negative NES for interstitial matrix and its phases (GO:0005614, GO:0140151, GO:0140149, GO:0140150), elastic fiber (GO:0071953), microfibril (GO:0001527), basement membrane (GO:0005604), and collagen trimer/complex (GO:0005581, GO:0098644). Second, increased proliferation/mitotic machinery and protein turnover/immune signals are supported by positive NES for spindle midzone (GO:0051233), spindle microtubule (GO:0005876), outer kinetochore (GO:0000940), kinesin complex (GO:0005871), proteasome complex (GO:0000502), and immunoglobulin complex (GO:0019814).
Crosstalk and interactions:
The mitotic spindle/kinetochore/kinesin terms form a coherent module for chromosome segregation, and co-enrichment with proteasome complex suggests potential proteolytic control of mitotic progression, although direct regulatory edges are not provided. ECM/basement-membrane loss may remove structural constraints and correlate with the proliferative shift, but the enrichment data alone do not establish causal direction. Immunoglobulin complex enrichment may reflect an immune-microenvironment component co-occurring with high-grade proliferation, but no direct regulator links are present in the provided terms.
Hypothesis:
Progression from grade I to grade III breast cancer involves loss of normal epithelial/ECM architecture, evidenced by lower interstitial matrix, basement membrane, elastic fiber, microfibril, and collagen trimer components, together with increased mitotic spindle/kinetochore/kinesin-driven proliferation, proteasome activity, and immunoglobulin-complex signal. In this dataset, ECM remodelling manifests primarily as down-regulation/loss of normal matrix components rather than increased expression of matrix genes.
Narrative draft:
In the MAINZ cohort, grade III tumours are distinguished from grade I tumours by a coordinated shift in gene expression: extracellular-matrix and basement-membrane components are depleted (interstitial matrix, solid/gel/non-collagenous interstitial matrix phases, elastic fiber, microfibril, basement membrane, collagen trimer and complex of collagen trimers; all NES < 0), whereas mitotic spindle and chromosome-segregation modules are enriched (spindle midzone, spindle microtubule, outer kinetochore, kinesin complex; all NES > 0). Proteasome complex and immunoglobulin complex are also up-regulated in grade III. These changes are consistent with loss of differentiated epithelial architecture and basement-membrane integrity accompanying a proliferative, mitotically active tumour state. The ECM signature suggests that matrix remodelling is dominated by loss of normal interstitial and basement-membrane constituents rather than by increased matrix gene expression. The co-enrichment of proteasome and immunoglobulin complexes may indicate proteolytic and immune-microenvironment contributions, but upstream transcriptional regulators are not identifiable from the provided terms.
Network evidence:
No PPI data were provided; the refined network is therefore a module-level representation inferred from co-enriched GO cellular components. It contains an ECM/basement-membrane module (collagen trimer, interstitial matrix, basement membrane, elastic fiber, microfibril) with negative NES and a mitotic-spindle module (outer kinetochore, spindle microtubule, spindle midzone, kinesin complex) with positive NES; proteasome and immunoglobulin complexes are additional positively enriched modules but lack direct edge evidence in the provided terms.
## Interpretation Result
### Cluster: Default
### 1. Overview
Grade III versus grade I breast tumours show coordinated down-regulation of extracellular matrix and basement-membrane components (interstitial matrix GO:0005614; basement membrane GO:0005604; elastic fiber GO:0071953; collagen trimers GO:0098644/GO:0005581) and up-regulation of mitotic spindle/kinetochore/kinesin modules (spindle midzone GO:0051233; spindle microtubule GO:0005876; outer kinetochore GO:0000940; kinesin complex GO:0005871) plus proteasome and immunoglobulin complexes (GO:0000502; GO:0019814).
### 2. Regulatory Drivers
No transcription-factor or upstream-regulator term is present among the provided top enriched terms, so specific master regulators cannot be nominated from this enrichment result alone. The positive enrichment of mitotic execution modules (spindle midzone, spindle microtubule, outer kinetochore, kinesin complex) suggests that proliferative cell-cycle regulators drive the grade III phenotype, while negative enrichment of ECM/basement-membrane programmes (interstitial matrix, basement membrane, collagen trimers) implies altered differentiation-state regulators, but these regulators are not directly represented in the provided data.
### 3. Key Mechanisms
Two major themes emerge. First, loss of ECM/basement-membrane architecture in grade III is supported by negative NES for interstitial matrix and its phases (GO:0005614, GO:0140151, GO:0140149, GO:0140150), elastic fiber (GO:0071953), microfibril (GO:0001527), basement membrane (GO:0005604), and collagen trimer/complex (GO:0005581, GO:0098644). Second, increased proliferation/mitotic machinery and protein turnover/immune signals are supported by positive NES for spindle midzone (GO:0051233), spindle microtubule (GO:0005876), outer kinetochore (GO:0000940), kinesin complex (GO:0005871), proteasome complex (GO:0000502), and immunoglobulin complex (GO:0019814).
### 4. Crosstalk & Interactions
The mitotic spindle/kinetochore/kinesin terms form a coherent module for chromosome segregation, and co-enrichment with proteasome complex suggests potential proteolytic control of mitotic progression, although direct regulatory edges are not provided. ECM/basement-membrane loss may remove structural constraints and correlate with the proliferative shift, but the enrichment data alone do not establish causal direction. Immunoglobulin complex enrichment may reflect an immune-microenvironment component co-occurring with high-grade proliferation, but no direct regulator links are present in the provided terms.
### 5. Hypothesis
Progression from grade I to grade III breast cancer involves loss of normal epithelial/ECM architecture, evidenced by lower interstitial matrix, basement membrane, elastic fiber, microfibril, and collagen trimer components, together with increased mitotic spindle/kinetochore/kinesin-driven proliferation, proteasome activity, and immunoglobulin-complex signal. In this dataset, ECM remodelling manifests primarily as down-regulation/loss of normal matrix components rather than increased expression of matrix genes.
### 6. Narrative Draft
In the MAINZ cohort, grade III tumours are distinguished from grade I tumours by a coordinated shift in gene expression: extracellular-matrix and basement-membrane components are depleted (interstitial matrix, solid/gel/non-collagenous interstitial matrix phases, elastic fiber, microfibril, basement membrane, collagen trimer and complex of collagen trimers; all NES < 0), whereas mitotic spindle and chromosome-segregation modules are enriched (spindle midzone, spindle microtubule, outer kinetochore, kinesin complex; all NES > 0). Proteasome complex and immunoglobulin complex are also up-regulated in grade III. These changes are consistent with loss of differentiated epithelial architecture and basement-membrane integrity accompanying a proliferative, mitotically active tumour state. The ECM signature suggests that matrix remodelling is dominated by loss of normal interstitial and basement-membrane constituents rather than by increased matrix gene expression. The co-enrichment of proteasome and immunoglobulin complexes may indicate proteolytic and immune-microenvironment contributions, but upstream transcriptional regulators are not identifiable from the provided terms.
### 7. Refined Regulatory Network
Key Interactions:
outer kinetochore -- spindle microtubule (binding/attachment) - Co-enriched positive NES; the outer kinetochore mediates chromosome attachment to spindle microtubules (GO:0000940, GO:0005876).
kinesin complex -- spindle midzone (transport/localization) - Co-enriched positive NES; kinesin motors localize to and organize the spindle midzone (GO:0005871, GO:0051233).
spindle midzone -- spindle microtubule (organization) - Co-enriched positive NES; the spindle midzone is a microtubule-rich region of the mitotic spindle (GO:0051233, GO:0005876).
collagen trimer -- basement membrane (structural component) - Co-enriched negative NES; collagen trimers are structural constituents of basement membranes (GO:0005581, GO:0005604).
collagen trimer -- interstitial matrix (structural component) - Co-enriched negative NES; collagens are major components of interstitial matrix (GO:0005581, GO:0005614).
elastic fiber -- microfibril (structural component) - Co-enriched negative NES; microfibrils are structural constituents of elastic fibers (GO:0071953, GO:0001527).
**Network Evidence:**
No PPI data were provided; the refined network is therefore a module-level representation inferred from co-enriched GO cellular components. It contains an ECM/basement-membrane module (collagen trimer, interstitial matrix, basement membrane, elastic fiber, microfibril) with negative NES and a mitotic-spindle module (outer kinetochore, spindle microtubule, spindle midzone, kinesin complex) with positive NES; proteasome and immunoglobulin complexes are additional positively enriched modules but lack direct edge evidence in the provided terms.
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.