37  Task recipes

These short recipes start with a common input and point to the next analysis decision. The main chapters provide the statistical details and knowledge-source constraints.

37.1 RNA-seq differential genes

Use a thresholded differential-expression table for ORA when the cutoff is part of the study design. Keep the full signed statistic for GSEA as well. The canonical local example is datasets/de_table.tsv, generated from the Bioconductor airway data; its fields and regeneration command are recorded in datasets/readme.md. It uses gene_id, log2FoldChange, stat, and padj columns.

library(clusterProfiler)
library(org.Hs.eg.db)

de_table <- read.delim("datasets/de_table.tsv", stringsAsFactors = FALSE)
annotation_keys <- AnnotationDbi::keys(org.Hs.eg.db, keytype = "ENSEMBL")
eligible <- is.finite(de_table$stat) & de_table$gene_id %in% annotation_keys

universe <- de_table$gene_id[eligible]
selected <- de_table$gene_id[
    eligible & !is.na(de_table$padj) &
        de_table$padj < 0.05 &
        abs(de_table$log2FoldChange) >= 1
]

ora <- enrichGO(
    gene = selected,
    universe = universe,
    OrgDb = org.Hs.eg.db,
    keyType = "ENSEMBL",
    ont = "BP"
)

ranking <- de_table$stat[eligible]
names(ranking) <- de_table$gene_id[eligible]
ranking <- sort(ranking, decreasing = TRUE)
gsea <- gseGO(ranking, OrgDb = org.Hs.eg.db, keyType = "ENSEMBL", ont = "BP")
head(as.data.frame(ora)[, c("ID", "Description", "pvalue", "p.adjust", "Count")])
                   ID                          Description       pvalue
GO:0003013 GO:0003013           circulatory system process 6.250970e-09
GO:0001525 GO:0001525                         angiogenesis 1.735764e-08
GO:0007411 GO:0007411                        axon guidance 2.553542e-08
GO:0030198 GO:0030198    extracellular matrix organization 2.778551e-08
GO:0097485 GO:0097485           neuron projection guidance 2.913354e-08
GO:0043062 GO:0043062 extracellular structure organization 3.084790e-08
               p.adjust Count
GO:0003013 2.090662e-05    56
GO:0001525 2.090662e-05    55
GO:0007411 2.090662e-05    28
GO:0030198 2.090662e-05    35
GO:0097485 2.090662e-05    28
GO:0043062 2.090662e-05    35
head(as.data.frame(gsea)[, c("ID", "Description", "setSize", "NES", "pvalue", "p.adjust")])
                   ID                                        Description
GO:0050951 GO:0050951         sensory perception of temperature stimulus
GO:0098704 GO:0098704         carbohydrate import across plasma membrane
GO:0098708 GO:0098708            D-glucose import across plasma membrane
GO:0140271 GO:0140271               hexose import across plasma membrane
GO:0071276 GO:0071276                   cellular response to cadmium ion
GO:0045648 GO:0045648 positive regulation of erythrocyte differentiation
           setSize       NES       pvalue    p.adjust
GO:0050951      22 -2.110286 5.470924e-06 0.005057870
GO:0098704      58  1.983135 1.507972e-05 0.007604293
GO:0098708      58  1.983135 1.507972e-05 0.007604293
GO:0140271      58  1.983135 1.507972e-05 0.007604293
GO:0071276      18  1.971514 2.516067e-05 0.009304416
GO:0045648      28  1.936275 6.328391e-05 0.011581469

Check the identifier namespace and construct universe from genes eligible for detection, not from the significant subset. The first workflow shows the same route with plots, mapping-loss reporting, and object export.

37.2 No significant differential genes

When a thresholded list is empty or unstable, use the complete ranking rather than relaxing the cutoff after looking at the enrichment output:

ids <- de_table$gene_id[is.finite(de_table$stat)][1:40]
term2gene <- data.frame(
    term = rep(c("PathA", "PathB"), each = 20),
    gene = ids
)
term2name <- data.frame(
    term = c("PathA", "PathB"),
    name = c("Path A (demo)", "Path B (demo)")
)

ranking <- setNames(de_table$stat, de_table$gene_id)
ranking <- sort(ranking[is.finite(ranking)], decreasing = TRUE)
gsea <- GSEA(ranking, TERM2GENE = term2gene, TERM2NAME = term2name,
             pvalueCutoff = 1, minGSSize = 5)
head(as.data.frame(gsea)[, c("ID", "Description", "setSize", "NES", "pvalue", "p.adjust")])
         ID   Description setSize       NES    pvalue  p.adjust
PathA PathA Path A (demo)      20 -1.211922 0.2337662 0.4537037
PathB PathB Path B (demo)      20  1.014915 0.4537037 0.4537037

Report how the statistic was defined and whether its sign has a biological direction. Continue to ORA and GSEA foundations for score types and leading-edge interpretation.

37.3 Single-cell markers

For marker tables with one row per gene and cluster, retain the grouping column and run one enrichment analysis per cluster:

markers <- data.frame(
    gene = c("G1", "G2", "G3"),
    cluster = c("T_cell", "T_cell", "B_cell")
)

marker_term2gene <- data.frame(
    term = c("T_path", "T_path", "B_path"),
    gene = c("G1", "G2", "G3")
)
marker_term2name <- data.frame(
    term = c("T_path", "B_path"),
    name = c("T-cell program (demo)", "B-cell program (demo)")
)

cc <- compareCluster(
    geneCluster = split(markers$gene, markers$cluster),
    fun = "enricher",
    TERM2GENE = marker_term2gene,
    TERM2NAME = marker_term2name,
    pvalueCutoff = 1, qvalueCutoff = 1, minGSSize = 1
)
head(as.data.frame(cc)[, c("Cluster", "ID", "Description", "pvalue", "p.adjust")])
  Cluster     ID           Description    pvalue  p.adjust
1  B_cell B_path B-cell program (demo) 0.3333333 0.3333333
2  T_cell T_path T-cell program (demo) 0.3333333 0.3333333

Use the formula interface when marker metadata contain several grouping variables. Compare clusters with the compareCluster analysis and visualization chapters. Cell-type annotation requires marker evidence in addition to enrichment terms; the optional LLM route is described in interpretation.

37.4 Protein and UniProt identifiers

First map the identifier namespace, then keep the mapped and unmapped counts:

proteins <- c("P04637", "P38398", "Q9Y6N7")
background_map <- bitr(
    AnnotationDbi::keys(org.Hs.eg.db, keytype = "UNIPROT")[1:500],
    fromType = "UNIPROT",
    toType = "ENTREZID",
    OrgDb = org.Hs.eg.db
)

mapped <- bitr(
    proteins,
    fromType = "UNIPROT",
    toType = "ENTREZID",
    OrgDb = org.Hs.eg.db
)

mapped_ids <- unique(mapped$ENTREZID)
ora <- enrichGO(
    mapped_ids,
    OrgDb = org.Hs.eg.db,
    keyType = "ENTREZID",
    universe = unique(c(background_map$ENTREZID, mapped_ids)),
    ont = "BP",
    pvalueCutoff = 1, qvalueCutoff = 1, minGSSize = 1
)
head(as.data.frame(ora)[, c("ID", "Description", "pvalue", "p.adjust", "Count")])
                   ID                                               Description
GO:0010558 GO:0010558 negative regulation of macromolecule biosynthetic process
GO:0010605 GO:0010605    negative regulation of macromolecule metabolic process
GO:0071345 GO:0071345                    cellular response to cytokine stimulus
GO:0009890 GO:0009890               negative regulation of biosynthetic process
GO:0034097 GO:0034097                                      response to cytokine
GO:1901652 GO:1901652                                       response to peptide
                 pvalue   p.adjust Count
GO:0010558 0.0001111420 0.05039202     3
GO:0010605 0.0002778550 0.05039202     3
GO:0071345 0.0005557099 0.05039202     3
GO:0009890 0.0015559878 0.05039202     3
GO:0034097 0.0015559878 0.05039202     3
GO:1901652 0.0015559878 0.05039202     3

Do not treat unmapped proteins as absent biology without checking the annotation release and species. Preserve the mapping table with the result. The universal enrichment chapter covers custom namespaces and annotation tables.

37.5 ChIP-seq and genomic regions

For peaks or other genomic regions, the first decision is how regions map to genes. Use nearest, host, flanking, promoter, or seq2gene() mappings according to the biological question, then pass the resulting gene sets to enrichment. See genomic coordination enrichment; the downstream ORA/GSEA choice remains the same as for other gene-level inputs.

37.6 Protein-plus-transcript or other multi-omics data

Keep each layer’s identifiers, direction, and uncertainty separate until the integration method is chosen. Use early feature-level aggregation when the layers measure compatible features; use late pathway-level fusion when the layers should retain separate evidence. The multi-omics chapter documents aggregate_omics(), ID harmonization, and contribution tracing.

37.7 Non-model organisms

For organisms without a standard OrgDb, construct explicit TERM2GENE and TERM2NAME tables from a documented annotation source. The non-model recipe gives an eggNOG-based route; package reusable annotations as GSON when they will be shared or reused.

37.8 Next steps