library(DOSE)
data(geneList)
de <- names(geneList)[abs(geneList) > 2]
edo <- enrichDO(de)14 Visualization of functional enrichment result
The enrichplot package implements several visualization methods to help interpret enrichment results. It supports visualizing enrichment results obtained from DOSE (Yu et al. 2015), clusterProfiler (Yu et al. 2012; Wu et al. 2021), ReactomePA (Yu and He 2016) and meshes (Yu 2018). Both over representation analysis (ORA) and gene set enrichment analysis (GSEA) are supported.
Note: Several visualization methods were first implemented in DOSE and rewrote from scratch using ggplot2. If you want to use the old methods, you can use the doseplot package.
14.1 Bar Plot
Bar plot is the most widely used method to visualize enriched terms. It depicts the enrichment scores (e.g. p values) and gene count or ratio as bar height and color Figure 14.1. Users can specify the number of terms (most significant) or selected terms (see also the FAQ) to display via the showCategory parameter.
library(enrichplot)
barplot(edo, showCategory=20) Other variables that derived using mutate can also be used as bar height or color as demonstrated in Section 18.4 and Figure 14.1.
library(clusterProfiler)
mutate(edo, qscore = -log(p.adjust, base=10)) |>
barplot(x="qscore")library(clusterProfiler) is not optional here. mutate() on an enrichResult works through an S3 method that clusterProfiler registers on load, so if only DOSE (or another downstream package) has been loaded, mutate() falls back to dplyr::mutate and fails with no applicable method for 'mutate' applied to an object of class "enrichResult".
14.1.1 Split top categories by ontology
When visualizing GO enrichment across all ontologies (ont = "ALL"), it is often more informative to show the top terms within each ontology (BP/CC/MF) rather than selecting the top terms globally. barplot() forwards additional arguments via ..., and one useful option is split. Setting split = "ONTOLOGY" will perform the “top N per ontology” selection and makes it straightforward to color (and later facet) the results by the ontology.
library(enrichplot)
library(ggplot2)
p_nofacet <- barplot(
ego_all,
x = "Count",
showCategory = 15,
split = "ONTOLOGY"
) +
aes(fill = ONTOLOGY) +
scale_fill_brewer(palette = "Set2")
p_nofacet +
geom_text(aes(label = Count), nudge_x = 2) +
scale_y_discrete()
split = "ONTOLOGY" (no faceting).
Note that barplot() wraps long term descriptions onto multiple lines by default (controlled by label_format, default is 30). If you prefer to keep the y-axis labels on a single line, you can override the default scale by adding + scale_y_discrete() as shown above.
p_facet <- barplot(
ego_all,
x = "Count",
showCategory = 15,
split = "ONTOLOGY"
) +
aes(fill = ONTOLOGY) +
scale_fill_brewer(palette = "Set2") +
enrichplot::autofacet(by = "row", scales = "free")
p_facet + scale_y_discrete()
autofacet() splitting by ontology.
If the fortified data contain an ONTOLOGY column (as in the output generated by split = "ONTOLOGY"), you can facet the plot directly with + enrichplot::autofacet(), and it will automatically split panels by ontology.
14.2 Dot plot
Dot plot is similar to bar plot with the capability to encode another score as dot size.
library(ggplot2)
edo2 <- gseDO(geneList)
dotplot(edo, showCategory=30, label_format=NULL) + ggtitle("dotplot for ORA")
dotplot(edo2, showCategory=30, label_format=NULL) + ggtitle("dotplot for GSEA")
Note: The dotplot() function also works with compareCluster() output.
14.2.1 Highlighting specific pathways
The label_format parameter in dotplot() accepts a function to format y-axis labels, which can be used to highlight specific pathways. This is useful for emphasizing pathways of interest in publication figures.
First, perform enrichment analysis:
library(DOSE)
data(geneList)
de <- names(geneList)[1:100]
x <- enrichDO(de)Define a function to highlight pathways containing “cancer”:
f <- function(ids) {
i <- grepl('cancer', ids)
ids[i] <- paste0("<i style='color:#009E73'>", ids[i], "</i>")
return(ids)
}To render the HTML formatting, use the ggtext package with element_markdown():
library(ggplot2)
library(ggtext)
library(enrichplot)
dotplot(x, label_format=f) + theme(axis.text.y = element_markdown())
You can also highlight specific pathways by index:
f2 <- function(ids) {
ids[1] <- paste0("<i style='color:#009E73'>", ids[1], "</i>")
ids[2] <- paste0("<i style='color:#0072B2'>**", ids[2], "**</i>")
return(ids)
}
dotplot(x, label_format=f2) + theme(axis.text.y = element_markdown())
Or apply different colors to all pathways:
f3 <- function(ids) {
# cols <- rcartocolor::carto_pal(length(ids), "Vivid")
# cols <- colorspace::rainbow_hcl(length(ids))
cols <- rainbow(length(ids))
ids <- paste0("<i style='color:", cols, "'>**", ids, "**</i>")
return(ids)
}
dotplot(x, label_format=f3) + theme(axis.text.y = element_markdown())
14.2.2 Formula interface of dotplot
The x variable of dotplot() supports a formula interface, allowing users to use derived variables. For example, we can calculate GeneRatio/BgRatio or -log(p.adjust) and use them as the x-axis variable.
p1 <- dotplot(edo, x = ~GeneRatio/BgRatio)
p2 <- dotplot(edo, x = ~ -log(p.adjust))
plot_list(p1, p2, ncol=2, tag_levels='A')
x = ~GeneRatio/BgRatio. (B) x = ~ -log(p.adjust).
This feature provides great flexibility in visualization. For instance, the Rich Factor is a common metric defined as the ratio of the number of differentially expressed genes annotated in a pathway to the total number of genes annotated in that pathway.
\[Rich Factor = \frac{Count}{BgRatio \times N}\]
where \(N\) is the total number of genes in the background distribution.
We can easily visualize the Rich Factor using the formula interface.
N <- 8007 # total number of genes in the background
dotplot(edo, x = ~Count/(BgRatio * N))Alternatively, you can precompute derived variables with mutate() (e.g., add neglog10p = -log10(p.adjust) or richFactor = Count / (BgRatio * N)) and then pass the new column name to x. This approach is useful when you want to reuse the computed variables across multiple plots. See the dedicated section on using dplyr verbs with enrichment results in Section 18.4.
14.2.3 Adjusting dot sizes
To increase the size of the dots in the dot plot, we can use the scale_size function from the ggplot2 package. This allows us to specify the range of the dot sizes.
p1 <- dotplot(edo, showCategory=10) + ggtitle("Default")
p2 <- dotplot(edo, showCategory=10) + scale_size(range=c(2, 20)) + ggtitle("Adjusted size")
plot_list(p1, p2, ncol=2, tag_levels='A')
scale_size(range=c(2, 20)).
14.3 Manhattan plot
The manhattanplot() function provides a landscape-style overview of enrichment results. It arranges enriched terms along the x-axis by ontology (or by category in compareCluster() results), uses -log10() transformed significance values on the y-axis, and maps term size to a variable such as Count. The most significant terms can be labeled automatically with showCategory, making it useful for summarizing a large number of enriched terms in a single plot.
When GO enrichment is performed with ont = "ALL", manhattanplot() can display the top terms across BP, CC, and MF simultaneously by combining it with split = "ONTOLOGY".
manhattanplot(
ego_all,
color = "p.adjust",
showCategory = 12,
size = "Count",
split = "ONTOLOGY",
title = "GO enrichment landscape"
)
-log10(p.adjust).
The manhattanplot() function also supports gseaResult, enrichResultList, gseaResultList, and compareClusterResult objects, which makes it convenient for comparing the distribution of significant terms across multiple ontologies or experimental groups.
14.4 Dotplot2: Comparing two clusters
The dotplot2() function is a specialized version of dotplot designed for comparing two selected clusters from compareCluster results. This function provides focused visualization when users want to directly contrast specific biological conditions or experimental groups.
# Example using compareCluster result
library(clusterProfiler)
data(gcSample)
xx <- compareCluster(gcSample, fun="enrichKEGG",
organism="hsa", pvalueCutoff=0.05)
# Compare cluster A and cluster B
dotplot2(xx, vars = c("A", "B"))The dotplot2() function accepts the same parameters as dotplot() but requires cluster1 and cluster2 arguments to specify which clusters to compare. This focused visualization helps in identifying differential enrichment patterns between specific experimental conditions.
14.5 Gene-Concept Network
Both the barplot() and dotplot() only displayed most significant or selected enriched terms, while users may want to know which genes are involved in these significant terms. To consider the potential biological complexities in which a gene may belong to multiple annotation categories and provide information about numeric changes when available, we developed the cnetplot() function to extract the complex associations. The cnetplot() depicts the linkages of genes and biological concepts (e.g. GO terms or KEGG pathways) as a network. GSEA result is also supported with only core enriched genes displayed. For GSEA results, cnetplot specifically visualizes the core enriched genes identified through leading edge analysis, which represent the subset of genes most responsible for driving the enrichment signal.
Users can use the fc_threshold parameter in the cnetplot function to filter genes to only include those with |foldChange| > fc_threshold. This allows users to plot only high fold-change genes (e.g., fc_threshold = 1). In addition, users can use the node_label parameter to label genes that are shared between groups of terms (i.e. node_label = "share").
Too many labels overlapping? cnetplot() draws gene and category labels with ggrepel, which hides labels that would collide. If you want every label shown, raise the limit — options(ggrepel.max.overlaps = Inf). If the network is simply too dense, it is usually better to reduce what gets drawn rather than to force all labels out: label only the categories with node_label = "category", or keep only the genes you care about with fc_threshold (e.g. fc_threshold = 1 keeps |foldChange| > 1) or by passing a foldChange vector that contains just those genes.
## convert gene ID to Symbol
edox <- setReadable(edo, 'org.Hs.eg.db', 'ENTREZID')
p1 <- cnetplot(edox, foldChange=geneList)
p2 <- cnetplot(edox, categorySizeBy=~ -log10(pvalue), foldChange=geneList)
## only plot high fold-change genes
p3 <- cnetplot(edox, foldChange=geneList, fc_threshold = 2)
p4 <- cnetplot(edox, foldChange=geneList, node_label = "share")
plot_list(p1, p2, p3, p4, ncol=2, tag_levels = 'A')
If you would like label subset of the nodes, you can use the node_label parameter, which supports 4 possible selections (i.e. “category”, “gene”, “all” and “none”), as demonstrated in Figure 14.14.
The node_label parameter also supports enhanced functionality:
- Vector selection: Specify specific genes to label using a character vector, e.g.,
node_label = c("gene1", "gene2") - ‘exclusive’ labeling: Label only genes that belong exclusively to one category using
node_label = "exclusive" - ‘share’ labeling: Label genes that are shared between multiple categories using
node_label = "share"(as shown in the example above) - Conditional filtering: Use comparison operators to filter genes based on fold change, e.g.,
node_label = "> 1"ornode_label = "< -1"to label genes with absolute fold change greater than 1
p1 <- cnetplot(edox, node_label="category")
p2 <- cnetplot(edox, node_label="gene")
p3 <- cnetplot(edox, node_label="all")
p4 <- cnetplot(edox, node_label="none",
color_category='firebrick',
color_item='steelblue')
plot_list(p1, p2, p3, p4, ncol=2, tag_levels = 'A')
The cnetplot function can be used as a general method to visualize data relationships in a network diagram. Please refer to the vignette of ggtangle.
14.5.1 Customizing gene colors in cnetplot
Users can color specific genes by providing a named vector of fold changes containing only those genes.
foldChange <- c(rep(1, 6), rep(-1, 4))
names(foldChange) <- c("MARCO", "GZMB", "CXCL11", "CXCL10",
"LAG3", "CCL8", "PDK1", "GABRP", "MELK", "CENPE")
p <- cnetplot(edox, foldChange=foldChange)
p
If you want to remove the color legend:
p <- p + guides(color='none')
p
To create a custom legend, you can use a dummy data frame and geom_point:
d <- data.frame(type=c('upregulated', 'downregulated'), x=0, y=0)
g <- p + geom_point(aes(alpha=type, x=x, y=y), data=d, shape=16, size=0)
g + guides(alpha=guide_legend(
override.aes=list(color=c("blue", "red"),
alpha=1,
size=3),
title = "VIP genes",
reverse = TRUE
))
14.5.2 Tuning cnetplot labels and edge colors
Besides the size of the category / item nodes (size_category / size_item) and the fold-change coloring above, two dedicated knobs are commonly needed:
color_edge = "category"colors each edge by the category (term) it belongs to, instead of mapping every edge to a single color.color_edgealso accepts a single color, e.g.color_edge = "grey".geom_cnet_label()is the label layer behindcnetplot()(fromggtangle). Becausecnetplot()returns aggplotobject, you can add anothergeom_cnet_label()layer on top withsize,colorandfontfaceto fine-tune the label text of a specific node type (e.g. only the category labels), which is useful when the last few labels are clipped or you want a distinct label style.
p1 <- cnetplot(edox, node_label = "none", showCategory = 4)
p2 <- cnetplot(edox, node_label = "none", showCategory = 4, color_edge = "category") +
geom_cnet_label(node_label = "category", size = 5, color = "firebrick", fontface = "bold")
plot_list(p1, p2, ncol = 2, tag_levels = "A")
geom_cnet_label() (B).
14.5.2.1 Node and edge colors
The default edge color (a red/green pair on gene-ID plots in some releases) can be replaced either with a single color (color_edge = "grey55") or, when you want the edges colored separately for each term and shown in the legend, with color_edge = "category". Category node colors are set with color_category, and — unlike some older versions — setting foldChange for the gene nodes no longer overrides color_category: the term nodes keep their color_category fill while the gene nodes get the fold-change palette.
p1 <- cnetplot(edox, node_label = "none", color_edge = "grey55",
color_category = "steelblue")
p2 <- cnetplot(edox, node_label = "none", color_edge = "category")
p3 <- cnetplot(edox, node_label = "category", foldChange = geneList) +
scale_color_gradient(low = "blue", high = "red")
plot_list(p1, p2, p3, ncol = 3, tag_levels = "A")
scale_color_gradient() (C).
Because these plots are ordinary ggplot objects, any ggplot2 scale can be plugged on top to re-map the aesthetics — scale_color_*() for gene / label colors and scale_fill_*() for node fills. Match the scale to the aesthetic: foldChange puts a continuous value (the fold change) on the colour aesthetic, so re-colour it with a continuous scale such as scale_color_gradient() / scale_color_gradient2(); a discrete scale_color_manual() would fail with “Continuous value supplied to a discrete scale”. For emapplot(), node colors are mapped to the fill aesthetic, so use scale_fill_manual(values = ...) (or scale_fill_gradient*()) there to customize the node colors.
Note: The cnetplot() function also works with compareCluster() output.
14.6 Heatmap-like functional classification
The heatplot is similar to cnetplot, but displays the relationships as a heatmap. The gene-concept network may become too complicated if users want to show a large number of significant terms. The heatplot can simplify the result and make it easier to identify expression patterns.
p1 <- heatplot(edox, showCategory=5)
p2 <- heatplot(edox, foldChange=geneList, showCategory=5)
plot_list(p1, p2, ncol=1, tag_levels = 'A')
foldChange=geneList (B)
The showTop parameter can be used to limit the number of genes displayed in the heatmap. This is particularly useful when dealing with large gene sets where only the top genes based on fold change or significance need to be visualized. For example, showTop = 20 will display only the top 20 genes in the heatmap.
14.7 Tree plot
The treeplot() function performs hierarchical clustering of enriched terms. It relies on the pairwise similarities of the enriched terms calculated by the pairwise_termsim() function, which by default uses Jaccard’s similarity index (JC). Users can also use semantic similarity values when supported (e.g., GO, DO, and MeSH).
The default agglomeration method in treeplot() is ward.D, and users can specify other methods via the cluster_method parameter (e.g., ‘average’, ‘complete’, ‘median’, ‘centroid’, etc.; see also the documentation of the hclust() function). The treeplot() function will cut the tree into several subtrees (specified by the nCluster parameter, default is 5) and label subtrees using high-frequency words. This reduces the complexity of the enriched result and improves user interpretation ability.
For fine-grained control over text appearance, the fontsize_tiplab and fontsize_cladelab parameters allow users to adjust the font size of tip (leaf) labels and clade labels respectively. For example, fontsize_tiplab = 8, fontsize_cladelab = 10 will set tip labels to 8pt and clade labels to 10pt.
14.7.1 Colouring the clades and wrapping labels
Two things that are easy to get wrong here:
Colouring the clades. Pass group_color a vector of colours, and it is used for both the clade labels and the highlight bars, so the two stay in sync:
treeplot(edox2, showCategory = 20,
group_color = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442"))Wrapping the tip labels. label_format controls the wrap length of the clade labels (the word-cloud summaries of each cluster), not the tip labels — passing it a number or a function will not change the terms along the tips:
## wraps the clade labels only
treeplot(edox2, label_format = 15)The tip labels are taken straight from the Description of the enrichment result, so to wrap them you have to shorten the descriptions before building the tree. pairwise_termsim() names the similarity matrix from Description, and treeplot() labels the tips from that matrix, so wrapping at this stage propagates all the way through:
edox_wrapped <- edox
edox_wrapped@result$Description <- stringr::str_wrap(edox_wrapped@result$Description, width = 30)
treeplot(pairwise_termsim(edox_wrapped), showCategory = 20)str_wrap() inserts newlines at word boundaries, which the labels render as line breaks. Note that this rewrites the term descriptions of the object, so do it on a copy if you still need the original wording elsewhere.
library(ggtree)
edox2 <- pairwise_termsim(edox)
p1 <- treeplot(edox2, cladelab_offset=8, tiplab_offset=.3, fontsize_cladelab =5) +
hexpand(.2)
p2 <- treeplot(edox2, cluster_method = "average",
cladelab_offset=14, tiplab_offset=.3, fontsize_cladelab =5) +
hexpand(.3)
aplot::plot_list(p1, p2, tag_levels='A', ncol=2)
hclust_method = "average" (B)
14.8 Semantic Space Plot
While treeplot() visualizes semantic similarity using a hierarchical structure, the ssplot() (Semantic Space Plot) projects enriched terms into a low-dimensional space (e.g., using Multidimensional Scaling, MDS). This provides a complementary spatial view where terms with high semantic similarity are clustered together.
Like treeplot(), ssplot() requires pairwise_termsim() to be run first. It automatically groups terms into clusters (default nCluster is determined automatically) and labels them with representative words.
ssplot(edox2, nCluster=5) + ggtitle("ssplot (nCluster=5)")
ssplot() is particularly useful for visualizing the overall semantic landscape of enrichment results and identifying distinct functional modules in a continuous space.
14.9 Enrichment Map
Enrichment map organizes enriched terms into a network with edges connecting overlapping gene sets. In this way, mutually overlapping gene sets are tend to cluster together, making it easy to identify functional module.
14.9.1 Handling GO term redundancy
GO annotations often contain redundant terms that can dominate enrichment results, potentially obscuring other biological stories. The simplify() function from the clusterProfiler package uses semantic similarity (via the GOSemSim package) to remove redundant GO terms, providing a clearer view of distinct functional modules.
# Load required packages
library(clusterProfiler)
library(enrichplot)
library(DOSE)
# Prepare example data
data(geneList, package="DOSE")
de <- names(geneList)[abs(geneList) > 2]
# Perform GO enrichment analysis
ego <- enrichGO(de, OrgDb = "org.Hs.eg.db", ont="BP", readable=TRUE)
# Remove redundant GO terms using simplify()
ego_simplified <- simplify(ego, cutoff=0.7, by="p.adjust", select_fun=min)
# Visualize both original and simplified results
ego <- pairwise_termsim(ego)
ego_simplified <- pairwise_termsim(ego_simplified)
p1 <- emapplot(ego, node_label_size=.8, size_edge=.5) +
scale_fill_continuous(low = "#e06663", high = "#327eba", name = "p.adjust",
guide = guide_colorbar(reverse = TRUE, order=1), trans='log10') +
ggtitle("Original GO terms")
p2 <- emapplot(ego_simplified, node_label_size=.8, size_edge=.5) +
scale_fill_continuous(low = "#e06663", high = "#327eba", name = "p.adjust",
guide = guide_colorbar(reverse = TRUE, order=1), trans='log10') +
ggtitle("After removing redundant terms")
# Combine plots
library(patchwork)
p1 + p2 + plot_layout(ncol = 2)
The simplify() function removes redundant GO terms based on semantic similarity (default cutoff = 0.7). This reveals distinct functional modules that might be obscured by redundant terms in the original enrichment results. The pairwise_termsim() function calculates pairwise similarities between terms, which is required for emapplot() visualization.
The emapplot function supports results obtained from hypergeometric test and gene set enrichment analysis. The size_category parameter can be used to resize nodes and the layout parameter can adjust the layout, as demonstrated in Figure 14.24.
edo <- pairwise_termsim(edo)
p1 <- emapplot(edo) # node_label = "category" (default)
p2 <- emapplot(edo, node_label = "none")
p3 <- emapplot(edo, node_label = "none", size_category=1.5)
p4 <- emapplot(edo, node_label = "none", layout="with_fr")
plot_list(p1, p2, p3, p4,
ncol=2, tag_levels = 'A',
design="AAAAAA\nBBCCDD",
heights = c(1, .3))
emapplot. default (A), node_label="none" (B), size_category=1.5 (C), and layout="with_fr" (D)
The node_label parameter controls how the labels were displayed. The enriched terms will be displayed by default with node_label="category" and it can be disabled by setting node_label="none".
The node_label_size parameter allows users to adjust the font size of node labels in the enrichment map. This is particularly useful when dealing with many overlapping terms or when labels need to be more readable. For example, node_label_size = 3 will increase the label size compared to the default.
If node_label="group", the emapplot function will cluster the enriched terms into different groups and only group names (determined by wordcloud) will be displayed. If node_label="all", then the enriched terms and group names will be displayed simultaneously.
p5 <- emapplot(edo, node_label = "group")
p6 <- emapplot(edo, node_label = "all")
plot_list(p5, p6,
ncol=1, tag_levels = 'A')
emapplot with enriched terms clustering. node_label="group" (A) and node_label="all" (B).
14.10 Biological theme comparison
The emapplot function also supports results obtained from compareCluster function of clusterProfiler package. In addition to size_category and layout parameters, the number of circles in the bottom left corner can be adjusted using the legend_n parameteras, and proportion of clusters in the pie chart can be adjusted using the pie parameter, when pie="count", the proportion of clusters in the pie chart is determined by the number of genes, as demonstrated in Figure 14.26.
library(clusterProfiler)
data(gcSample)
xx <- compareCluster(gcSample, fun="enrichKEGG",
organism="hsa", pvalueCutoff=0.05)
xx <- pairwise_termsim(xx)
p1 <- emapplot(xx)
p2 <- emapplot(xx)
p3 <- emapplot(xx, pie="count")
p4 <- emapplot(xx, pie="count", size_category=1.5, layout="kk")
plot_list(p1, p2, p3, p4, ncol=2, tag_levels = 'A')
compareCluster function of clusterProfiler package. default (A), legend_n=2 (B), pie="count" (C) and pie="count", size_category=1.5, layout="kk" (D).
14.11 UpSet Plot
The upsetplot is an alternative to cnetplot for visualizing the complex association between genes and gene sets. It emphasizes the gene overlapping among different gene sets.
upsetplot(edo)
For over-representation analysis, upsetplot will calculate the overlaps among different gene sets as demonstrated in Figure 14.27. For GSEA result, it will plot the fold change distributions of different categories (e.g. unique to pathway, overlaps among different pathways).
kk2 <- gseKEGG(geneList = geneList,
organism = 'hsa',
minGSSize = 120,
pvalueCutoff = 0.05,
verbose = FALSE)
upsetplot(kk2) 
14.12 ridgeline plot for expression distribution of GSEA result
The ridgeplot will visualize expression distributions of core enriched genes for GSEA enriched categories. It helps users to interpret up/down-regulated pathways.
The stat parameter controls the geometry used to summarize the distributions. By default, ridgeplot() uses stat = "density_ridges", and users can also set stat = "binline" to display binned profiles instead.
ridgeplot(edo2)
ridgeplot(edo2, stat = "binline")14.13 running score and preranked list of GSEA result
Running score and preranked list are traditional methods for visualizing GSEA result. The enrichplot package supports both of them to visualize the distribution of the gene set and the enrichment score.
p1 <- gseaplot(edo2, geneSetID = 1, by = "runningScore", title = edo2$Description[1])
p2 <- gseaplot(edo2, geneSetID = 1, by = "preranked", title = edo2$Description[1])
p3 <- gseaplot(edo2, geneSetID = 1, title = edo2$Description[1])
plot_list(p1, p2, p3, ncol=1, tag_levels='A')
by = "runningScore"). by = "runningScore" (A), by = "preranked" (B), default (C)
The gseaplot function also allows users to customize the colors of the running score line and the rank line.
p_default <- gseaplot(edo2, geneSetID = 1, title = edo2$Description[1])
p_custom <- gseaplot(edo2, geneSetID = 1, color="#DAB546", color.line='firebrick',
color.vline="steelblue", title = edo2$Description[1])
plot_list(p_default, p_custom, ncol=2, tag_levels='A')
color, color.line, and color.vline.
14.13.1 Unified color setting with set_enrichplot_color()
For consistent color settings across different plot types in enrichplot, the set_enrichplot_color() function provides a unified approach to customize color mappings. This avoids repetitive color code specifications in different plotting functions.
14.13.1.1 Function introduction
The set_enrichplot_color() function can be added to any enrichplot visualization using the + operator (ggplot2 style). It offers flexible control over color transformations and palettes.
14.13.1.2 Parameters
type: Type of color aesthetic to modify (default:"fill", can also be"color"for line colors)transform: Transformation to apply to the values used for color mapping"log10": Apply log10 transformation (default for recent versions)"identity": Use original values without transformation
colors: Custom color palette as a vector of colors (e.g.,c("red", "blue"))
14.13.1.3 Usage examples
# Example data
library(DOSE)
data(geneList)
de <- names(geneList)[abs(geneList) > 2]
edo <- enrichDO(de)
# A: Default log10 transformation
p1 <- dotplot(edo, showCategory=15) +
set_enrichplot_color(type='fill', transform='log10') +
ggtitle("Default log10 transform")
# B: Custom colors with log10 transformation
p2 <- dotplot(edo, showCategory=15) +
set_enrichplot_color(type='fill', transform='log10', colors=c("red", "blue")) +
ggtitle("Red-blue palette")
# C: Identity transformation (no transformation)
p3 <- dotplot(edo, showCategory=15) +
set_enrichplot_color(type='fill', transform='identity') +
ggtitle("Identity transform")
library(aplot)
plot_list(p1, p2, p3, ncol=3, tag_levels='A')
set_enrichplot_color(). (A) Default log10 transformation. (B) Custom red-blue palette with log10 transformation. (C) Identity transformation (no transformation).
Note: Starting from enrichplot v1.29.2, log10 transformation of p-values is applied by default in functions like dotplot(). The set_enrichplot_color() function allows users to override this default behavior or customize color palettes consistently across all visualization types.
Another method to plot GSEA result is the gseaplot2 function:
Reading the three panels. gseaplot2() stacks three panels that share the same x-axis, which is the position in the ranked gene list (not a gene identifier):
- Running Enrichment Score — the running-sum statistic. It rises when a gene in the set is encountered and falls otherwise; the enrichment score is its maximum deviation from zero (marked by the dashed vertical line).
- The middle band — the black vertical ticks are the positions of the genes in the set. The coloured rectangles behind them divide the hits (not the ranked list) into equal-count bins, running from the top of the list (red) to the bottom (blue); their widths therefore encode local hit density, with narrow rectangles indicating a region where the set’s genes are packed tightly, and the final rectangle extending to the end of the list. A gene set enriched at the top shows narrow red rectangles crowded on the left.
- Ranked List Metric — the ranking statistic itself (e.g. the log fold change or the signed correlation), with the gene set members drawn as vertical segments.
The coloured rectangles used to be wrong — the bug is fixed in this build. Earlier versions computed the hit bins from a cumulative count taken from the bottom of the ranked list and never reversed it, so the gradient came out backwards: for a strongly enriched gene set nearly all hits collapsed into a single rectangle and the rest were empty. That is the bug behind enrichplot#20 and enrichplot#221. It is fixed in the development version of enrichplot, which is the version this book builds against: measured on a gene set with 40 hits in a 1000-gene list, all 9 rectangles are populated and the largest holds 15% of the hits, instead of one rectangle holding all of them. On an older release the master branch still carries the old order, so there the red/blue split is not a faithful picture of where the hits are; if the rectangles look collapsed, that is why. The coloured rectangles are in any case a property of gseaplot2(), not of the GSEA method — the original Broad GSEA plots colour the ranked-list-metric panel by the sign of the metric instead, so the two conventions do not agree either.
gseaplot2(edo2, geneSetID = 1, title = edo2$Description[1])
The gseaplot2 also supports multile gene sets to be displayed on the same figure:
gseaplot2(edo2, geneSetID = 1:3)
User can also displaying the pvalue table on the plot via pvalue_table parameter:
gseaplot2(edo2, geneSetID = 1:3, pvalue_table = TRUE,
color = c("#E495A5", "#86B875", "#7DB0DD"), ES_geom = "dot")
The pvalue_table can be customized using the following parameters:
pvalue_table_rownamesto specify row names of the table (if NULL, no row names will be displayed)pvalue_table_columnsto specify column names of the table
gseaplot2(edo2, geneSetID = 1, pvalue_table = TRUE,
pvalue_table_rownames = NULL,
pvalue_table_columns = c("ID", "NES", "p.adjust"))
User can specify subplots to only display a subset of plots:
p1 <- gseaplot2(edo2, geneSetID = 1:3, subplots = 1)
p2 <- gseaplot2(edo2, geneSetID = 1:3, subplots = 1:2)
plot_list(p1, p2, ncol=1, tag_levels = 'A')
subplots = 1 (A),subplots = 1:2 (B)
14.13.2 Labeling genes in GSEA plot
Users can use geom_gsea_gene() to label specific genes in the GSEA plot.
library(ggplot2)
library(ggrepel)
# Get gene set ID
id <- edo2$ID[1]
# Randomly select genes to label
set.seed(123)
genes <- sample(edo2[[id]], 5)
# Label genes on gseaplot2
p <- gseaplot2(edo2, geneSetID = 1, title = edo2$Description[1])
# Add geom_gsea_gene layer to the first subplot (running score)
p[[1]] <- p[[1]] + geom_gsea_gene(genes, geom=geom_label)
p
If users prefer to label with gene symbols, they can convert the gene IDs to symbols (e.g., using setReadable) before plotting.
library(clusterProfiler)
# Assuming org.Hs.eg.db is available
if (require("org.Hs.eg.db")) {
edo2_symbol <- setReadable(edo2, 'org.Hs.eg.db', 'ENTREZID')
id <- edo2_symbol$ID[1]
genes_symbol <- sample(edo2_symbol[[id]], 5)
p_symbol <- gseaplot2(edo2_symbol, geneSetID = 1, title = edo2_symbol$Description[1])
p_symbol[[1]] <- p_symbol[[1]] + geom_gsea_gene(genes_symbol, geom=geom_text_repel)
p_symbol
}The gsearank function plot the ranked list of genes belong to the specific gene set.
gsearank(edo2, 1, title = edo2[1, "Description"])
Multiple gene sets can be aligned using cowplot:
library(ggplot2)
pp <- lapply(1:3, function(i) {
anno <- edo2[i, c("NES", "pvalue", "p.adjust")]
lab <- paste0(names(anno), "=", round(anno, 3), collapse="\n")
gsearank(edo2, i, edo2[i, 2]) + xlab(NULL) +ylab(NULL) +
annotate("text", 10000, edo2[i, "enrichmentScore"] * .75, label = lab, hjust=0, vjust=0)
})
plot_list(gglist=pp, ncol=1)
14.13.3 Extracting data from gsearank plot
The gsearank function can also output the data used for plotting by setting output = "table". This allows users to inspect the running score and other metrics for genes in the gene set.
gsearank(edo2, 1, output = "table") |> head() gene rank in geneList running ES core enrichment
1 9837 60 0.06243348 YES
2 1503 194 0.09767718 YES
3 7037 235 0.13600863 YES
4 3932 276 0.17107219 YES
5 3559 298 0.20599919 YES
6 51311 316 0.24005196 YESUsers can also merge this table with gene information (e.g. Symbol) using bitr or other methods.
rank_table <- gsearank(edo2, 1, output = "table")
# Assuming 'gene' column contains Entrez IDs
if (require("org.Hs.eg.db") && require("clusterProfiler")) {
gene_info <- bitr(rank_table$gene, fromType="ENTREZID", toType="SYMBOL", OrgDb="org.Hs.eg.db")
rank_table_symbol <- merge(rank_table, gene_info, by.x="gene", by.y="ENTREZID")
head(rank_table_symbol)
}14.14 pubmed trend of enriched terms
One of the problem of enrichment analysis is to find pathways for further investigation. Here, we provide pmcplot function to plot the number/proportion of publications trend based on the query result from PubMed Central. Of course, users can use pmcplot in other scenarios. All text that can be queried on PMC is valid as input of pmcplot.
## Europe PMC is an external service and does intermittently answer 503.
## retry() (defined in _common.R) tries a few times and returns NULL if it is
## still unavailable, so a third-party outage skips this figure instead of
## failing the entire book build.
terms <- edo$Description[1:5]
p <- retry(pmcplot(terms, 2010:2020))
p2 <- retry(pmcplot(terms, 2010:2020, proportion=FALSE))
if (!is.null(p) && !is.null(p2)) {
plot_list(p, p2, ncol=2)
}
14.15 Volcano plot for ORA results
The volplot() function provides volcano plot visualization for over-representation analysis (ORA) results, allowing users to visualize the relationship between fold change (or other effect size metrics) and statistical significance.
# Example using enrichment result with fold change information
# Note: volplot requires results with fold change values
volplot(edo)
The volplot() function accepts standard ggplot2 aesthetics and can be customized with additional parameters for point size, color, and labeling thresholds.
14.16 Horizontal GSEA plot
The hplot() function creates horizontal versions of GSEA plots, which can be useful for comparing multiple gene sets side-by-side or for presentations where horizontal layout is preferred.
# Example using GSEA result
edo2 <- gseDO(geneList)
hplot(edo2, geneSetID = 1:3)
The hplot() function supports the same parameters as gseaplot2() but arranges the output horizontally. This is particularly useful for comparing multiple pathways or when vertical space is limited.
14.17 GO Graph
The goplot() function can accept the output of enrichGO and visualize the enriched GO induced graph. See Figure 14.44 for an example.
library(clusterProfiler)
library(org.Hs.eg.db)
data(geneList, package = "DOSE")
de <- names(geneList)[abs(geneList) > 2]
ego <- enrichGO(de, OrgDb = org.Hs.eg.db, ont = "BP", readable = TRUE)goplot(ego)
The plotGOgraph() function is another option to visualize the GO topology.
plotGOgraph(ego)14.18 Importing results from other tools
The visualization methods described above work with any enrichResult or gseaResult object, and results produced by other enrichment tools can be brought into this framework with two routes:
- Tool-specific importers:
import_enrichr(),import_gprofiler2(),import_webgestalt()andimport_fgsea()map the output tables of popular tools to the canonical schema (see Section 14.18.5) and return enrichment objects directly. - Generic constructors:
as_enrichResult()andas_gseaResult(), provided by theenrichitpackage and re-exported byenrichplot, convert any result table that follows the canonical column schema.
Both routes return standard enrichResult / gseaResult objects, so every visualization introduced in this chapter (barplot(), dotplot(), cnetplot(), gseaplot(), …) works on them.
The import functions ship with enrichplot (development version 2.0.0), and the as_enrichResult() / as_gseaResult() constructors are provided by the enrichit package, which enrichplot re-exports. Every example in this section therefore runs as written.
14.18.1 enrichr
The enrichR::enrichr() function queries the enrichr web service and returns a list of result tables (one per database) with columns Term, Overlap, P.value, Adjusted.P.value and Genes:
library(enrichR)
res <- enrichr(de_genes, databases = c("KEGG_2021_Human", "GO_Biological_Process_2023"))import_enrichr() accepts either that list (use the db parameter to select a database) or a single result table. Here we use a small table that mimics the enrichr output, and pass the query genes (gene) and the background gene list (universe) so that GeneRatio, BgRatio and related statistics can be derived:
enrichr_table <- data.frame(
Term = c("GO:0006915;;apoptotic process", "COVID-19"),
Overlap = c("5/150", "3/60"),
P.value = c(1.2e-4, 3.4e-3),
Adjusted.P.value = c(0.011, 0.079),
Genes = c("CASP3; BAX; BCL2; TP53; FAS", "ACE2; TMPRSS2; IL6"),
stringsAsFactors = FALSE
)
de_genes <- c("CASP3", "BAX", "BCL2", "TP53", "FAS",
"ACE2", "TMPRSS2", "IL6", "TNF", "MYC")
edo_er <- import_enrichr(enrichr_table, gene = de_genes,
universe = paste0("Gene", 1:2000))
edo_er#
# over-representation test
#
#...@organism UNKNOWN
#...@ontology UNKNOWN
#...@gene chr [1:10] "CASP3" "BAX" "BCL2" "TP53" "FAS" "ACE2" "TMPRSS2" "IL6" "TNF" ...
#...pvalues adjusted by 'BH' with cutoff < 1
#...2 enriched terms found
'data.frame': 2 obs. of 12 variables:
$ ID : chr "GO:0006915" "COVID-19"
$ Description : chr "apoptotic process" "COVID-19"
$ GeneRatio : chr "5/10" "3/10"
$ BgRatio : chr "150/2000" "60/2000"
$ RichFactor : num 0.0333 0.05
$ FoldEnrichment: num 6.67 10
$ zScore : num 5.11 5.02
$ pvalue : num 0.00012 0.0034
$ p.adjust : num 0.011 0.079
$ qvalue : num 0.011 0.079
$ geneID : chr "CASP3/BAX/BCL2/TP53/FAS" "ACE2/TMPRSS2/IL6"
$ Count : num 5 3
#...Citation
S Xu, E Hu, Y Cai, Z Xie, X Luo, L Zhan, W Tang, Q Wang, B Liu, R Wang, W Xie, T Wu, L Xie, G Yu. Using clusterProfiler to characterize multiomics data. Nature Protocols. 2024, 19(11):3292-3320 The result is a standard enrichResult object, ready for any visualization:
aplot::plot_list(
dotplot(edo_er, showCategory = 10) + ggtitle("enrichr result"),
cnetplot(edo_er, showCategory = 2),
ncol = 2, tag_levels = "A"
)
14.18.2 g:Profiler
gprofiler2::gost() returns a list whose result element is a data.frame with columns term_id, term_name, p_value, query_size, intersection_size, term_size, effective_domain_size and intersections:
library(gprofiler2)
gostres <- gost(de_genes, organism = "hsapiens")gost_table <- data.frame(
term_id = c("GO:0006915", "GO:0002376"),
term_name = c("apoptotic process", "immune system process"),
p_value = c(1.2e-4, 3.4e-3),
adjusted_p_value = c(0.011, 0.079),
query_size = c(10L, 10L),
intersection_size = c(5L, 3L),
term_size = c(150L, 60L),
effective_domain_size = c(2000L, 2000L),
intersections = I(list(c("CASP3", "BAX", "BCL2", "TP53", "FAS"),
c("CASP3", "BAX", "TNF"))),
stringsAsFactors = FALSE
)
x_gp <- import_gprofiler2(list(result = gost_table))
as.data.frame(x_gp)[, c("ID", "Description", "GeneRatio", "BgRatio", "pvalue", "Count")] ID Description GeneRatio BgRatio pvalue Count
GO:0006915 GO:0006915 apoptotic process 5/10 150/2000 0.00012 5
GO:0002376 GO:0002376 immune system process 3/10 60/2000 0.00340 3If the table has no query_size column, pass the query genes with the gene parameter for exact GeneRatio values.
dotplot(x_gp, showCategory = 10) + ggtitle("g:Profiler result")
14.18.3 WebGestalt
WebGestaltR() returns a summary table with columns geneSet, description, size, overlap, rawPValue, adjPValue and userIds:
library(WebGestaltR)
wg <- WebGestaltR(enrichMethod = "ORA", organism = "hsapiens",
enrichDatabase = "pathway_KEGG", interestGene = de_genes,
interestGeneType = "genesymbol", referenceGene = background,
referenceGeneType = "genesymbol")wg_table <- data.frame(
geneSet = c("hsa04210", "hsa04110"),
description = c("Apoptosis", "Cell cycle"),
size = c(60L, 120L),
overlap = c(5L, 3L),
rawPValue = c(1.2e-4, 3.4e-3),
adjPValue = c(0.011, 0.079),
userIds = c("CASP3;BAX;BCL2;TP53;FAS", "CASP3;BAX;TNF"),
stringsAsFactors = FALSE
)
x_wg <- import_webgestalt(wg_table, gene = de_genes,
universe = paste0("Gene", 1:2000))dotplot(x_wg, showCategory = 10) + ggtitle("WebGestalt result")
14.18.4 fgsea
The fgsea() function (Korotkevich et al. 2019) takes a named statistics vector and a list of gene sets, and returns a data.frame with columns pathway, ES, NES, pval, padj, size and leadingEdge:
library(fgsea)
data(geneList, package = "DOSE")
pathways <- split(gene_ids, term_ids) # named list of gene sets
fgres <- fgsea(pathways, stats = geneList)import_fgsea() additionally requires the ranked statistics vector (stats), as most GSEA visualizations consume the ranked list. Passing the pathways list used as input (via the geneSets parameter) enables exact running-score plots. Here we convert a small fgsea-style result table, with ranked statistics taken from the DOSE example data:
data(geneList, package = "DOSE")
stats <- sort(geneList, decreasing = TRUE)
set1 <- names(stats)[1:40]
set2 <- names(stats)[200:240]
fgres <- data.frame(
pathway = c("Set1", "Set2"),
ES = c(0.62, -0.48),
NES = c(2.1, -1.7),
pval = c(1e-4, 0.01),
padj = c(2e-4, 0.01),
size = c(40L, 41L),
stringsAsFactors = FALSE
)
fgres$leadingEdge <- list(
names(stats)[1:10],
rev(names(stats))[1:8]
)
gx <- import_fgsea(fgres, stats = stats, geneSets = list(Set1 = set1, Set2 = set2))
gx#
# Gene Set Enrichment Analysis
#
#...@organism UNKNOWN
#...@setType UNKNOWN
#...@geneList Named num [1:12495] 4.57 4.51 4.42 4.14 3.88 ...
- attr(*, "names")= chr [1:12495] "4312" "8318" "10874" "55143" ...
#...nPerm
#...2 enriched terms found
'data.frame': 2 obs. of 11 variables:
$ ID : chr "Set1" "Set2"
$ enrichmentScore: num 0.62 -0.48
$ NES : num 2.1 -1.7
$ pvalue : num 1e-04 1e-02
$ p.adjust : num 2e-04 1e-02
$ setSize : int 40 41
$ Description : chr "Set1" "Set2"
$ qvalue : num 2e-04 1e-02
$ core_enrichment: chr "4312/8318/10874/55143/55388/991/6280/2305/9493/1062" "4969/57758/79901/79838/10974/10551/5241/4239"
$ rank : int 40 240
$ leading_edge : chr "tags=100%, list=0%, signal=100%" "tags=100%, list=2%, signal=98%"
#...Citation
S Xu, E Hu, Y Cai, Z Xie, X Luo, L Zhan, W Tang, Q Wang, B Liu, R Wang, W Xie, T Wu, L Xie, G Yu. Using clusterProfiler to characterize multiomics data. Nature Protocols. 2024, 19(11):3292-3320 gseaplot(gx, geneSetID = "Set1")
14.18.5 Arbitrary result tables
Tools that do not have a dedicated importer (e.g. DAVID, or in-house scripts) can be converted with the generic constructors. as_enrichResult() requires a data.frame with ID, pvalue and at least one of geneID, Count or GeneRatio; the remaining canonical columns are derived when the query genes (gene), the universe (universe) and optionally the gene sets (geneSets) are supplied. Column aliases (e.g. PValue, term_id, padj, FDR) are recognized automatically:
| column | required | content |
|---|---|---|
ID |
yes | term identifier |
Description |
no | term name (defaults to ID) |
pvalue |
yes | raw p-value |
p.adjust |
no | adjusted p-value (computed with pAdjustMethod if absent) |
qvalue |
no | q-value (estimated if absent) |
geneID |
no | overlap genes, separated by / |
Count |
no | number of overlap genes |
GeneRatio |
no | k/n: overlap size / query size |
BgRatio |
no | M/N: gene set size / universe size |
custom <- data.frame(
term_id = c("custom:001", "custom:002"),
name = c("my in-house gene set", "another custom set"),
PValue = c(0.001, 0.02),
FDR = c(0.01, 0.05),
Genes = c("CASP3, BAX, BCL2", "ACE2, TMPRSS2, IL6"),
stringsAsFactors = FALSE
)
x_ct <- as_enrichResult(custom, gene = de_genes)
as.data.frame(x_ct)[, c("ID", "Description", "geneID", "Count", "GeneRatio",
"pvalue", "p.adjust")] ID Description geneID Count GeneRatio
custom:001 custom:001 my in-house gene set CASP3/BAX/BCL2 3 3/10
custom:002 custom:002 another custom set ACE2/TMPRSS2/IL6 3 3/10
pvalue p.adjust
custom:001 0.001 0.01
custom:002 0.020 0.05dotplot(x_ct, showCategory = 10) + ggtitle("custom table")
For GSEA-type results (with ranked statistics), use as_gseaResult() instead, passing the ranked vector via the geneList parameter:
gx2 <- as_gseaResult(fgres, geneList = stats,
geneSets = list(Set1 = set1, Set2 = set2))References
Korotkevich, Gennady, Vladimir Sukhov, and Alexey Sergushichev. 2019. “Fast Gene Set Enrichment Analysis.” bioRxiv, October 22, 060012. https://doi.org/10.1101/060012.
Wu, Tianzhi, Erqiang Hu, Shuangbin Xu, et al. 2021. “clusterProfiler 4.0: A Universal Enrichment Tool for Interpreting Omics Data.” The Innovation 2 (3): 100141. https://doi.org/10.1016/j.xinn.2021.100141.
Yu, Guangchuang. 2018. “Using Meshes for MeSH Term Enrichment and Semantic Analyses.” Bioinformatics 34 (21): 3766–67. https://doi.org/10.1093/bioinformatics/bty410.
Yu, Guangchuang, and Qing-Yu He. 2016. “ReactomePA: An r/Bioconductor Package for Reactome Pathway Analysis and Visualization.” Molecular BioSystems 12 (2): 477–79. https://doi.org/10.1039/C5MB00663E.
Yu, Guangchuang, Le-Gen Wang, Yanyan Han, and Qing-Yu He. 2012. “clusterProfiler: An r Package for Comparing Biological Themes Among Gene Clusters.” OMICS: A Journal of Integrative Biology 16 (5): 284–87. https://doi.org/10.1089/omi.2011.0118.
Yu, Guangchuang, Li-Gen Wang, Guang-Rong Yan, and Qing-Yu He. 2015. “DOSE: An r/Bioconductor Package for Disease Ontology Semantic and Enrichment Analysis.” Bioinformatics 31 (4): 608–9. https://doi.org/10.1093/bioinformatics/btu684.