This chapter answers the comparison question: which biological themes differ among gene lists, clusters, or experimental groups? compareCluster() reruns the same enrichment function for each group and returns one compareClusterResult, so the groups share an explicit analysis design and can be compared without treating independently generated p-values as if they were a common table.
| Aspect | Comparison analysis |
|---|---|
| Input | A named list of gene vectors, or one gene-per-row data frame with formula grouping columns. |
| Engine | One shared enrichment function (enrichGO(), enrichKEGG(), GSEA(), or another supported function) applied per group. |
| Output | A compareClusterResult retaining group labels, terms, gene mappings, and per-group statistics. |
| Key limitation | Groups must be compared under an explicit identifier namespace, universe/background, annotation source, and multiple-testing design. |
The clusterProfiler package was developed for biological theme comparison (Yu et al. 2012; Wu et al. 2021), and it provides a function, compareCluster, to automatically calculate enriched functional profiles of each gene clusters and aggregate the results into a single object. Comparing functional profiles can reveal functional consensus and differences among different experiments and helps in identifying differential functional modules in omics datasets.
The compareCluster() function applies selected function (via the fun parameter) to perform enrichment analysis for each gene list.
data(gcSample)
str(gcSample) List of 8
$ X1: chr [1:216] "4597" "7111" "5266" "2175" ...
$ X2: chr [1:805] "23450" "5160" "7126" "26118" ...
$ X3: chr [1:392] "894" "7057" "22906" "3339" ...
$ X4: chr [1:838] "5573" "7453" "5245" "23450" ...
$ X5: chr [1:929] "5982" "7318" "6352" "2101" ...
$ X6: chr [1:585] "5337" "9295" "4035" "811" ...
$ X7: chr [1:582] "2621" "2665" "5690" "3608" ...
$ X8: chr [1:237] "2665" "4735" "1327" "3192" ...Users can use a named list of gene IDs as the input that passed to the geneCluster parameter.
ck <- compareCluster(geneCluster = gcSample, fun = enrichGO,
OrgDb = org.Hs.eg.db, ont="BP")
ck <- setReadable(ck, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck) Cluster ID
1 X1 GO:0061844
2 X1 GO:0019730
3 X1 GO:0086005
4 X1 GO:1990266
5 X1 GO:0061337
6 X1 GO:0086091
Description
1 antimicrobial humoral immune response mediated by antimicrobial peptide
2 antimicrobial humoral response
3 ventricular cardiac muscle cell action potential
4 neutrophil migration
5 cardiac conduction
6 regulation of heart rate by cardiac conduction
GeneRatio BgRatio RichFactor FoldEnrichment oddsRatio zScore pvalue
1 9/198 142/18842 0.06338028 6.031370 6.627641 6.201898 1.953016e-05
2 10/198 189/18842 0.05291005 5.035006 5.487044 5.745341 3.230336e-05
3 5/198 37/18842 0.13513514 12.859678 15.068005 7.441356 4.040326e-05
4 7/198 93/18842 0.07526882 7.162702 7.908560 6.139580 5.620341e-05
5 7/198 95/18842 0.07368421 7.011909 7.727987 6.053734 6.440834e-05
6 5/198 42/18842 0.11904762 11.328764 13.028287 6.905722 7.556233e-05
p.adjust qvalue
1 0.03040124 0.004425159
2 0.03040124 0.004425159
3 0.03040124 0.004425159
4 0.03040124 0.004425159
5 0.03040124 0.004425159
6 0.03040124 0.004425159
geneID Count
1 H2BC11/CXCL8/BPI/MUC7/CXCL5/CXCL3/PLA2G1B/KNG1/CCL20 9
2 H2BC11/CXCL8/PI3/BPI/MUC7/CXCL5/CXCL3/PLA2G1B/KNG1/CCL20 10
3 KCND3/KCNE5/SCN3B/CTNNA3/KCNH2 5
4 CXCL8/S100A8/FUT7/CXCL5/GP2/CXCL3/PLA2G1B 7
5 KCND3/ATP2B2/KCNE5/ACE2/SCN3B/CTNNA3/KCNH2 7
6 KCND3/KCNE5/SCN3B/CTNNA3/KCNH2 5As an alternative to using named list, the compareCluster() function also supports passing a formula to describe more complicated experimental designs (e.g., \(Gene \sim time + treatment\)).
data(geneList, package='DOSE')
mydf <- data.frame(Entrez=names(geneList), FC=geneList)
mydf <- mydf[abs(mydf$FC) > 1,]
mydf$group <- "upregulated"
mydf$group[mydf$FC < 0] <- "downregulated"
mydf$othergroup <- "A"
mydf$othergroup[abs(mydf$FC) > 2] <- "B"
formula_res <- compareCluster(Entrez~group+othergroup, data=mydf,
fun=enrichGO, OrgDb = org.Hs.eg.db, ont="BP")
head(formula_res) Cluster group othergroup ID
1 downregulated.A downregulated A GO:0030198
2 downregulated.A downregulated A GO:0043062
3 downregulated.A downregulated A GO:0045229
4 downregulated.A downregulated A GO:0001501
5 downregulated.A downregulated A GO:0008015
6 downregulated.A downregulated A GO:0030199
Description GeneRatio BgRatio RichFactor
1 extracellular matrix organization 38/631 278/18842 0.1366906
2 extracellular structure organization 38/631 279/18842 0.1362007
3 external encapsulating structure organization 38/631 280/18842 0.1357143
4 skeletal system development 35/631 341/18842 0.1026393
5 blood circulation 42/631 480/18842 0.0875000
6 collagen fibril organization 14/631 65/18842 0.2153846
FoldEnrichment oddsRatio zScore pvalue p.adjust qvalue
1 4.081656 4.798328 9.635419 1.617110e-13 2.761734e-10 1.696848e-10
2 4.067027 4.778152 9.607167 1.814206e-13 2.761734e-10 1.696848e-10
3 4.052502 4.758143 9.579047 2.034177e-13 2.761734e-10 1.696848e-10
4 3.064865 3.436170 7.162604 4.392006e-09 4.472160e-06 2.747757e-06
5 2.612797 2.893481 6.662545 1.479521e-08 1.205218e-05 7.405025e-06
6 6.431501 8.079575 8.165131 2.435878e-08 1.653555e-05 1.015967e-05
geneID
1 4856/10631/2200/1294/1290/11096/3249/8292/2621/1289/1287/165/1513/1842/2191/4016/10218/4035/1471/2006/4313/7049/7177/2331/7373/50509/1281/10266/4319/3913/2192/1300/11117/81035/3908/1296/7075/1805
2 4856/10631/2200/1294/1290/11096/3249/8292/2621/1289/1287/165/1513/1842/2191/4016/10218/4035/1471/2006/4313/7049/7177/2331/7373/50509/1281/10266/4319/3913/2192/1300/11117/81035/3908/1296/7075/1805
3 4856/10631/2200/1294/1290/11096/3249/8292/2621/1289/1287/165/1513/1842/2191/4016/10218/4035/1471/2006/4313/7049/7177/2331/7373/50509/1281/10266/4319/3913/2192/1300/11117/81035/3908/1296/7075/1805
4 4856/2200/1462/4653/5125/55112/1289/3479/7227/63923/1009/3952/2487/6424/5364/2690/5744/116039/284266/658/50509/2697/1281/1746/7704/7481/2121/2202/1909/4488/6926/1300/4982/85477/81029
5 6833/2200/8864/5125/2946/5350/26509/477/7220/7349/4922/3952/857/3751/185/5348/6505/2006/6863/55800/4313/9590/5138/9365/2697/63895/18/1359/6444/10266/3357/150/3778/1909/11117/775/9607/6387/290/7122/23327/23171
6 4856/2200/1294/1290/1289/1287/165/4016/2331/7373/50509/1281/11117/1805
Count
1 38
2 38
3 38
4 35
5 42
6 14compareCluster() takes gene lists, not enrichment tables. This is the single most common misunderstanding about the formula interface. Both the named list and the formula form expect genes, and compareCluster() runs the enrichment test itself, once per group:
data must be a data frame with one row per gene — a gene identifier column, plus the grouping columns on the right-hand side of the formula.fun (enrichGO() here) from each group’s gene list. They are not carried over from any input table.So if you already have two tables of enrichment results (ID, Description, GeneRatio, pvalue, p.adjust, …), there is nothing to join: those numbers were produced by whatever run created them, on whatever background, and compareCluster() has no way to reuse them. To compare the two conditions, combine the gene lists that produced those tables and let compareCluster() re-run the test:
## one row per gene, with a column naming the group it belongs to
genes <- rbind(
data.frame(gene = condition_A_genes, group = "A"),
data.frame(gene = condition_B_genes, group = "B")
)
res <- compareCluster(gene ~ group, data = genes, fun = "enrichGO",
OrgDb = org.Hs.eg.db, ont = "BP")Re-running the test is not a limitation to work around — it is what makes the two columns comparable. Enrichment p-values depend on the gene list, the universe/background and the multiple-testing correction, so p-values computed separately for A and B are not on a common footing and plotting them side by side would be misleading.
If you want to carry a pre-computed gene-level statistic into the comparison instead (a ranked list rather than a cut-off list), use the GSEA route: compareCluster(..., fun = "gseGO").
The compareCluster() function is highly versatile and can be directly integrated into single-cell analysis workflows. For example, after identifying marker genes for each cluster using Seurat, the results can be directly used for functional enrichment analysis.
Typically, FindAllMarkers() from Seurat returns a data frame where rows are genes and columns include cluster information and statistical metrics.
This example is wrapped in retry() (defined in _common.R) because its fun reaches the remote KEGG service, so a transient failure is retried a few times before giving up; the enrichGO examples further on stay unwrapped because they run locally.
library(Seurat)
library(SeuratData)
library(dplyr)
library(clusterProfiler)
# Load example data
data("pbmc3k")
sce <- pbmc3k.final
# Identify marker genes
sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE,
min.pct = 0.25,
thresh.use = 0.25)
# Filter markers
markers <- sce.markers |>
group_by(cluster) |>
filter(p_val_adj < 0.001) |>
ungroup()
# ID conversion (if needed)
gid <- bitr(unique(markers$gene), 'SYMBOL', 'ENTREZID', OrgDb = 'org.Hs.eg.db')
markers <- full_join(markers, gid, by = c('gene' = 'SYMBOL'))
# Perform comparison using formula interface
x <- retry(compareCluster(ENTREZID ~ cluster, data = markers, fun = 'enrichKEGG'))
# Visualization
dotplot(x, label_format = 40) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))Enrichment plots need enrichment result objects, not a table of pathway statistics. A frequent question is how to hand enrichplot a data frame of per-pathway statistics (PValue, FDR, NES, logFC, …) obtained elsewhere. There is no way to do this directly: dotplot(), gseaplot2(), ridgeplot() and friends dispatch on enrichResult / gseaResult / compareClusterResult objects, and those carry the gene sets, the gene-to-set mapping and the background that the plots need. A bare table has none of that.
The way in is to re-run the enrichment inside clusterProfiler on the same gene list and gene sets, which gives you an object the plotting functions understand:
FindAllMarkers() gives you) → compareCluster(ENTREZID ~ cluster, data = markers, fun = "enrichGO"), then dotplot() / cnetplot() / emapplot(), as above.gseGO() / GSEA(), then gseaplot2() / ridgeplot() / gseaplot().If the statistics came from another tool (fgsea, GSEA desktop, a vendor report), re-running the equivalent analysis in clusterProfiler with the same gene list and gene sets is usually a couple of lines, and is the only way to get the object the visualisation layer requires.
Methods like COSG return marker genes as a list or data frame where columns represent clusters. Since a data frame is essentially a list of equal-length vectors, it can be passed directly to compareCluster().
library(COSG)
# Identify markers using COSG
marker_cosg <- cosg(
sce,
groups = 'all',
assay = 'RNA',
slot = 'data',
mu = 1,
n_genes_user = 100
)
# The first element is a data frame of gene symbols
# Columns correspond to clusters
head(marker_cosg[[1]])
# Directly use the data frame for enrichment
y <- compareCluster(marker_cosg[[1]],
fun = 'enrichGO',
OrgDb = 'org.Hs.eg.db',
keyType = 'SYMBOL',
ont = "MF")
# Visualization
dotplot(y, label_format = 60) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))This demonstrates that compareCluster() can seamlessly handle various data structures commonly produced by single-cell analysis tools, simplifying the downstream functional interpretation.
compareClusterResult with comparison-aware dot and network plots.compareCluster(), see compare functional profiles.