Skip to main content

Functional Enrichment Workshop

Making sense of gene and proteins lists with functional enrichment analysis

16 clusterProfiler and enrichplot

clusterProfiler is a comprehensive suite of enrichment tools. It has functions to run ORA or GSEA over commonly used databases (GO, KEGG, KEGG Modules, DAVID, Pathway Commons, WikiPathways) as well as universal enrichment functions to perform ORA or GSEA with novel species or custom gene sets. We will use these universal tools in the final activity of this workshop, focusing on the supported organisms and databases for the present activity.

One of the key advantages of using R over web tools is flexibility with visualisations.

The same authors have released a plotting package enrichplot dedicated to plotting enrichment results. In this activity, we will perform GSEA with clusterProfiler then explore many different visualisation options. At the end of the activity, we will have a poll to see which of the many plot types are the favourites! 😁

 

16.1 Supported databases, species and namespaces

One of the challenges when working with clusterProfiler for FEA is that each enrichment function has different supported organisms and different namespace requirements, so you can not necessarily use all of the functions over the same gene list. In this activity, we will review the FEA functions and investigate their requirements, before performing a gene ID conversion with the bitr function to enable compatibility with our Pezzini et al 2017 dataset.

Let’s start by reviewing the dedicated database enrichment functions and what inbuilt support they have.

Database GO KEGG KEGG Modules Pathway Commons WikiPathways DAVID
GSEA function gseGO gseKEGG gseMKEGG gsePC gseWP NA
ORA function enrichGO enrichKEGG enrichMKEGG enrichPC enrichWP enrichDAVID
Supported species Those with Bioconductor annotation package (20) KEGG organisms (thousands) KEGG organisms (thousands) Human, or convert to UniProt IDs Entrez Those supported by DAVID
Supported namespaces Differs depending on species kegg (compatible with Entrez), ncbi-geneid, ncbi-proteinid, or uniprot kegg (compatible with Entrez), ncbi-geneid, ncbi-proteinid, or uniprot hgnc or uniprot Entrez Those supported by DAVID

To obtain this information, 3 sources were required:

  1. The clusterProfiler user guide
  2. External websites (Bioconductor and KEGG Organisms)
  3. clusterProfiler functions (eg get_wp_organisms(), keytypes(org.Hs.eg.db))

This highlights the need to carefully review the tool you are using and explore the user guide and functionality.

 

16.2 Activity overview

Since we have covered ORA with gprofiler, we will perform a GSEA with clusterProfiler using gseKEGG.

  1. Explore the functions of clusterProfiler including which FEA functions support which organisms and which namespaces
  2. Load input dataset (a gene matrix with adjusted P values and log2 fold change values)
  3. Extract the gene IDs and sort by log2 fold change to create the GSEA ranked gene list R object
  4. Use bitr to convert gene IDs from ENSEMBL to ENTREZ for compatibility with gseKEGG
  5. Perform GSEA with gseKEGG
  6. Visualise results with many different plot types from enrichplot

 

16.3 Analysis workflow

The complete executable analysis is included below. Run the code chunks interactively while working through the chapter, or render the full Bookdown project to execute the analysis and retain the code, parameters, results, figures, and software-version information in the book output.

library(clusterProfiler)
#> 
#> clusterProfiler v4.20.0 Learn more at https://yulab-smu.top/contribution-knowledge-mining/
#> 
#> Please cite:
#> 
#> Guangchuang Yu, Li-Gen Wang, Yanyan Han and Qing-Yu He. clusterProfiler: an R package for
#> comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology.
#> 2012, 16(5):284-287
#> 
#> Attaching package: 'clusterProfiler'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(enrichplot)
#> enrichplot v1.32.0 Learn more at https://yulab-smu.top/contribution-knowledge-mining/
#> 
#> Please cite:
#> 
#> G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722
library(readr)
library(biomaRt)
library(org.Hs.eg.db)
#> Loading required package: AnnotationDbi
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> Loading required package: generics
#> 
#> Attaching package: 'generics'
#> The following object is masked from 'package:lubridate':
#> 
#>     as.difftime
#> The following object is masked from 'package:dplyr':
#> 
#>     explain
#> The following objects are masked from 'package:base':
#> 
#>     as.difftime, as.factor, as.ordered, intersect, is.element, setdiff, setequal, union
#> 
#> Attaching package: 'BiocGenerics'
#> The following object is masked from 'package:dplyr':
#> 
#>     combine
#> The following object is masked from 'package:fs':
#> 
#>     path
#> The following objects are masked from 'package:stats':
#> 
#>     IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#> 
#>     anyDuplicated, aperm, append, as.data.frame, basename, cbind, colnames, dirname,
#>     do.call, duplicated, eval, evalq, Filter, Find, get, grep, grepl, is.unsorted,
#>     lapply, Map, mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
#>     Position, rank, rbind, Reduce, rownames, sapply, saveRDS, table, tapply, unique,
#>     unsplit, which.max, which.min
#> Loading required package: Biobase
#> Welcome to Bioconductor
#> 
#>     Vignettes contain introductory material; view with 'browseVignettes()'. To cite
#>     Bioconductor, see 'citation("Biobase")', and for packages 'citation("pkgname")'.
#> Loading required package: IRanges
#> Loading required package: S4Vectors
#> 
#> Attaching package: 'S4Vectors'
#> The following object is masked from 'package:clusterProfiler':
#> 
#>     rename
#> The following object is masked from 'package:tidyr':
#> 
#>     expand
#> The following objects are masked from 'package:lubridate':
#> 
#>     second, second<-
#> The following objects are masked from 'package:dplyr':
#> 
#>     first, rename
#> The following object is masked from 'package:utils':
#> 
#>     findMatches
#> The following objects are masked from 'package:base':
#> 
#>     expand.grid, I, unname
#> 
#> Attaching package: 'IRanges'
#> The following object is masked from 'package:clusterProfiler':
#> 
#>     slice
#> The following object is masked from 'package:lubridate':
#> 
#>     %within%
#> The following objects are masked from 'package:dplyr':
#> 
#>     collapse, desc, slice
#> 
#> Attaching package: 'AnnotationDbi'
#> The following object is masked from 'package:clusterProfiler':
#> 
#>     select
#> The following object is masked from 'package:dplyr':
#> 
#>     select
#> 
library(dplyr)
library(ggupset)
library(ggplot2)

16.4 Load input data and extract ranked gene list

We will once again load the Pezzini RNAseq dataset.

data <- read_tsv("Pezzini_DE.txt", col_names = TRUE, show_col_types = FALSE)
head(data)

This time, instead of filtering for DEGs, we will extract all genes, and sort them by fold change (largest to smallest) to make our ranked gene list for GSEA.

What class of R object does our ranked gene list need to be in?

?gseKEGG

We can see geneList = order ranked geneList. Not informative!

Unfortunately, the required type of object is detailed under the enrichKEGG function, not the gseKEGG function! For enrichKEGG, the parameter gene is described as requiring “a vector of entrez gene id”, yet for gseKEGG, the description for geneList is “order ranked geneList”. There is a little bit of sleuthing required at times!

# Ranked gene list vector for GSEA
ranked <- setNames(data$Log2FC[order(-data$Log2FC)], data$Gene.ID[order(-data$Log2FC)])

# Inspect the vector
head(ranked)
tail(ranked)

16.5 Check gseKEGG function arguments and requirements

Review the parameters:

?gseKEGG

Most of those defaults look suitable to start.

We have human so the default organism = "hsa" argument is correct. If you were working with a species other than human, you first need to obtain your organism code. You can derive this from KEGG Organisms or using the clusterProfiler function search_kegg_organism.

Pick your favourite species and search for the KEGG organism code by editing the variable ‘fave’ then executing the code chunk:

fave <- "horse"

# search by common_name or scientific_name or kegg.code
search_kegg_organism(fave, by = "common_name")

Now back to the parameters - the defaults for P value correction and filtering, and gene set size limits are acceptable. We don’t want to use internal data (i.e. we want to search against the latest KEGG online) so FALSE is apt here, and we do want to use the fgsea algorithm for analysis.

For the seed parameter, we should provide a value, to ensure results are the same each time the command is run. By setting a seed, you fix the sequence of random numbers generated within the GSEA algorithm.

We also need to check the keyType (gene namespace) against the input data that we have. From the help page, we can see the supported namespaces for the KEGG database are one of ‘kegg’, ‘ncbi-geneid’, ‘ncib-proteinid’ or ‘uniprot’.

Going back to where we loaded the input data and ran head to view the first few rows, we can see our input has ENSEMBL gene IDs as well as official gene symbols. ENSEMBL gene IDs are generally preferable for bioinformatics analyses because they are more unique and stable compared to gene symbol.

Since our input data does not match any of the valid namespaces, we need to convert gene IDs! clusterProfiler has the bitr function to do this. BiomaRt is also a popular R package for this task.

16.6 Convert gene IDs

Check the usage for the bitr function:

?bitr

We need to understand what the valid fromType and toType values are, and it turns out we need an Org.db to use bitr! This is a Bioconductor annotation package, of which there are currently only 20. So while the gseKEGG function supports all organisms in KEGG, performing gene ID conversion within clusterProfiler may not be possible for non-model species and you would need to seek a different method.

We have already loaded the org.Hs.eg.db annotation library. We can use this to search keytypes:

keytypes(org.Hs.eg.db)

Converting to ENTREZ will be compatible with kegg. So our fromType is ENSEMBL and toType is ENTREZID. Note that these are case sensitive!

If we were to run enrichGO or gseaGO that require a Bioconductor Org.db package, we would not need to do a conversion, as both of the gene ID types within our input data (ENSEMBL and SYMBOL) are natively supported.

Gene ID conversion often results in duplicates. The below code performs the reformatting and handles duplicates by keeping only the first occurrence of duplicate Entrez IDs within the input. Note that this is not ideal and for a real experiment you should print out a list of duplicates, carefully review these and choose which and how to retain based on your biological context.

First, convert the ENSEMBL gene IDs from our ‘ranked’ list to ENTREZ gene IDs:

# Convert the Ensembl IDs to Entrez IDs, dropping NAs
converted_ids <- bitr(names(ranked),
  fromType = "ENSEMBL",
  toType = "ENTREZID",
  OrgDb = org.Hs.eg.db,
  drop = TRUE
)

<1% failing to map is pretty good.

The 1:many mappings warning means that some of our gene IDs matched more than 1 ENTREZ ID. We need to ensure that the final list we provide to GSEA does not contain duplicates. This needs to happen at two stages: first of all, ensuring that each ENSEMBL ID is mapped to only one ENTREZ ID, and then once the final converted vector has been created, check it for duplicated ENTREZ IDs, which could occur when two different ENSEMBL IDs from our input map to the the same ENTREZ ID.

The below code does this by selecting the first of each set of duplicates. This is not ideal. In a real experiment, you should print out the duplicates and directly manage how to handle them by reviewing the gene IDs involved and deciding whether it is valid to select one ID over another, or at times you may choose to merge the values for duplicate genes.

# Handle any `1:many mapping` (ENSEMBL ID matches >1 ENTREZID) by keeping only the first Entrez ID for each Ensembl ID (not ideal, see note above)
converted_ids <- converted_ids[!duplicated(converted_ids$ENSEMBL), ]

# Keep only the rows with valid Entrez IDs
converted_ids <- converted_ids[!is.na(converted_ids$ENTREZID), ]

# Filter the ranked gene list to include only the successfully mapped Ensembl IDs
ranked_entrez <- ranked[names(ranked) %in% converted_ids$ENSEMBL]

# Replace the Ensembl IDs with the corresponding Entrez IDs in the ranked vector
names(ranked_entrez) <- converted_ids$ENTREZID[match(names(ranked_entrez), converted_ids$ENSEMBL)]

# Remove duplicates from the vector by keeping only the first occurrence of each Entrez IDs (not ideal, see note above)
ranked_entrez <- ranked_entrez[!duplicated(names(ranked_entrez))]

# Inspect the vector
head(ranked_entrez)
tail(ranked_entrez)

Looks good! Now we are ready to enrich!

16.7 Run GSEA over KEGG database

Now run GSEA over the KEGG database with gseKEGG function, setting a seed value and kegg as our gene ID type. this may take a minute or two to run.

As in the previous gprofiler2 chapter, we’ll spell out every parameter explicitly below, even where the default would do - this keeps your parameter choices visible and documented, and protects you if a package’s defaults ever change between versions.

gsea_kegg <- gseKEGG(
  ranked_entrez, # ranked gene list vector object
  organism = "hsa", # species
  keyType = "kegg", # gene ID type/namespace
  exponent = 1, # weight of each step
  minGSSize = 10, # minimum gene set size
  maxGSSize = 500, # maximum gene set size
  eps = 1e-10, # sets the boundary for calculating the p value
  pvalueCutoff = 0.05, # adjusted P value cutoff for enriched terms
  pAdjustMethod = "BH", # multiple testing correction with BH method
  verbose = TRUE, # print output as it runs
  use_internal_data = FALSE, # use latest online KEGG data not local data
  seed = 123, # set seed for random number generation, for reproducibility
  by = "fgsea" # GSEA algorithm
)

We have 3 warnings, one about ties in the ranked list, which we could resolve manually by reviewing the raw data, or just ignore it as it is only 0.01% of the list.

The second warning may be resolved by following the suggestion to set nPermSimple = 10000. How frustrating that this parameter is not described within the gseKEGG help menu, nor the clusterProfiler PDF at all!

Let’s rerun following both suggestions, to use permutations (nPermSimple = 1000) and set eps = 0. We will keep the rest of the arguments the same as before.

gsea_kegg <- gseKEGG(
  ranked_entrez,
  organism = "hsa",
  keyType = "kegg",
  exponent = 1,
  minGSSize = 10,
  maxGSSize = 500,
  eps = 0, # changed to zero
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH",
  verbose = TRUE,
  use_internal_data = FALSE,
  seed = 123,
  by = "fgsea",
  nPermSimple = 10000 # added permutations
)

Great, those last 2 warnings have resolved and we only have the expected one about tied rankings.

16.8 Visualise GSEA KEGG results

16.8.1 Tabular results

First, let’s preview the results.

print(gsea_kegg)

At the line ## 'data.frame': 43 obs. of 11 variables, 43 is the number of significant enriched terms, and 11 is the number of columns in the output dataframe.

Extract results to a TSV file. This will print the significant enrichments, sorted by adjusted P value.

write.table(gsea_kegg, file = "clusterProfiler_gseKEGG_results.tsv", sep = "\t", quote = FALSE, row.names = FALSE)

You can view the table in the editor pane by clicking on it from the Files pane.

16.8.2 Plots

One of the key advantages of using R over web tools is flexibility with visualisations. Next we will produce a range of plot types from the enrichplot package, developed by the same authors as clusterProfiler.

In this package, some of the plots can be used to plot ORA or GSEA results, and others only for one type or the other.

To determine which plot functions are compatible with which FEA results type, you can check the help menu. If you see ## S4 method for signature 'enrichResult' this plot is compatible with ORA results object. If you see ## S4 method for signature 'gseaResult' this plot is compatible with GSEA results object.

Review the help menu for the dotplot function, which can plot both ORA and GSEA:

?enrichplot::dotplot

And review the help menu for ridgeplot which can only plot GSEA results:

?enrichplot::ridgeplot

Unfortunately, enrichplot is only compatible with the enrichment results from the packages produced by this development team, although the desire to use enrichplot with the output of other tools is widespread. The R package multienrichjam has a function enrichDF2enrichResult that converts ORA dataframe results from other FEA tools to the format required for enrichplot.

multienrichjam has a lot of dependencies and has not been installed on these VMs so we will not be performing this today. However, this functionality and flexibility is pretty cool, so if you wanted to install this on your own computer outside the workshop, below is the code for installing :-)

# github remotes package required:
# install.packages("remotes")
# library(remotes)

# install multienrichjam:
# remotes::install_github("jmw86069/multienrichjam", dependencies = TRUE)
# library(multienrichjam)

# check function help:
# ?enrichDF2enrichResult

16.8.2.1 dotplot

A dotplot is an easy to interpret plot. The size of the dot is proportional to the number of genes (count) in that term. A higher GeneRatio indicates that a larger portion of the gene set is significantly altered or associated with the observed changes in gene expression.

We can change the number of top terms that are plotted with the showCategory parameter.

enrichplot::dotplot(
  gsea_kegg,
  x = "geneRatio",
  color = "p.adjust",
  orderBy = "x",
  showCategory = 15,
  font.size = 8
) +
  ggtitle("GSEA KEGG")

16.8.2.2 upsetplot

The upset plot shows the pattern of shared genes among the enriched terms (since a gene can belong to multiple terms or pathways).

For GSEA, the plot includes fold change distributions on the Y axis.

n controls the number of terms plotted - here we have restricted to 10 for ease of viewing. Plot view is quite small within a notebook, but if you wanted to plot more terms, you could print to an A4 size output file for enhanced clarity.

# n = the number of top enrichments to plot
enrichplot::upsetplot(gsea_kegg, n = 15)

16.8.2.3 emapplot

The enrichment map plot organises enriched terms into a network with edges connecting overlapping gene sets. Overlapping gene sets cluster together.

It is required to run the function pairwise_termsim before producing the plot.

When running pairwise_termsim, we don’t need to save the output to a new object, because the information is added to the object and does not change any other attributes.

Each node is a term, and the number of genes associated with the term is shown by the dot size, with P values by dot colour.

# calculate pairwise similarities between the enriched terms
gsea_kegg <- enrichplot::pairwise_termsim(gsea_kegg)

# plot
enrichplot::emapplot(gsea_kegg, showCategory = 25)

16.8.2.4 treeplot

The treeplot provides the same information as the emapplot but presented in a different way.

Terms that share more genes or biological functions will be closer together in the tree structure. Clades are colour coded and ‘cluster tags’ assigned. You can control the number of words in the tag (default is 4). The user guide describes the argument nWords however running that will throw an error (it says ‘warning’ but it is fatal so to me that’s an error!):

“Warning: Use ‘cluster.params = list(label_words_n = your_value)’ instead of ‘nWords’. The nWords parameter will be removed in the next version.”

This plot also requires the pairwise similarity matrix calculation that emapplot does. Since we have already run it, it is hashed out in the code chunk below.

Plot view within the notebook affects the clarity of the cluster tags, however printing this plot to a file would resolve that.

Note also that enrichplot’s treeplot API has itself changed between versions - the cluster.params argument used above has since been replaced by a simpler nCluster argument:

# calculate pairwise similarities between the enriched terms
# gsea_kegg <- enrichplot::pairwise_termsim(gsea_kegg)

# earlier enrichplot versions: enrichplot::treeplot(gsea_kegg, showCategory = 15, color = "p.adjust", cluster.params = list(label_words_n = 5))
enrichplot::treeplot(gsea_kegg, showCategory = 15, color = "p.adjust", nCluster = 5)

16.8.2.5 cnetplot

The cnetplot depicts the linkages of genes and terms as a network. This is helpful to understanding which genes are involved in the enriched terms, adding a level of detail not offered by the plots we have generated so far.

For GSEA, where all genes (not just DEGs) are used, only the ‘core’ enriched genes are used to create the network plot. These are the ‘leading edge genes’, those genes up to the point where the Enrichment Score (ES) gets maximised from the base zero. In other words, the subset of genes that are most strongly associated with a specific term.

There are a few parameters to play around with here to get a readable plot.

Try changing the number of terms plotted (showCategory).

Try changing the node_label which controls whether labels are put on terms (category), genes (item), or both (all).

Since the information that is attempted to be plotted is complex, having a large number of terms and attempting to label everything won’t look very informative! If you want to plot both gene IDs and term names, you will need to plot a small number of categories.

Like treeplot, the cnetplot API has also changed between enrichplot versions - the cex_label_gene, cex_label_category and colorEdge arguments used previously have been removed, and label sizing/edge colouring is now handled by the companion ggtangle package:

# cex_label_gene to reduce the font size for gene ID
# node_label = one of ’all’, ’none’, ’category’ and ’item’
# showCategory = number of terms to plot

# earlier enrichplot versions: enrichplot::cnetplot(gsea_kegg, showCategory = 8, cex_label_gene = 0.5, cex_label_category = 0.8, colorEdge = TRUE, node_label = "category")

p <- enrichplot::cnetplot(
  gsea_kegg,
  showCategory = 8,
  color_edge = "category",
  node_label = "none"
) +
  ggtangle::geom_cnet_label(
    node_label = "category",
    size = 3.2
  ) +
  ggtangle::geom_cnet_label(
    node_label = "item",
    size = 2
  )

p

With 8 terms and no gene IDs, we don’t get any more detail than from the treeplot or emapplot.

This plot is really useful for showing a detailed look at a small number of terms. Just plotting the top 3 terms may not look very helpful (try plotting the top 3!)

A useful application is plotting the interaction between specific terms. This is helpful to obtain gene ID level resolution for term interactions of interest.

The 3 terms listed below are for ‘IL-17 signaling pathway’, ‘Viral protein interaction with cytokine and cytokine receptor’ and ‘Chemokine signaling pathway’ which are among the top 10 enrichments with a relationship of shared genes, evident from the plot above.

Run the code below or select a handful of terms of your choosing from the results table we printed earlier. We need the KEGG ID (column 1).

# Select terms of interest
select_terms <- c("hsa04657", "hsa04061", "hsa04062")

# Filter the gseaResult object based on the term IDs
# must keep it as a gsea result type and not as a dataframe when filtering
gsea_kegg_select <- gsea_kegg
gsea_kegg_select@result <- gsea_kegg@result[gsea_kegg@result$ID %in% select_terms, ]

# Plot with cnetplot
# earlier enrichplot versions: enrichplot::cnetplot(gsea_kegg_select, cex_label_gene = 0.5, cex_label_category = 0.8, colorEdge = TRUE, node_label = "all")

p <- enrichplot::cnetplot(
  gsea_kegg_select,
  color_edge = "category",
  node_label = "none"
) +
  ggtangle::geom_cnet_label(
    node_label = "item",
    size = 2
  ) +
  ggtangle::geom_cnet_label(
    node_label = "category",
    size = 3.2
  )

p

Now we can clearly see the gene IDs for this network.

16.8.2.6 heatplot

Like the upset plot, the heatplot shows shared genes across enriched terms. When there are a lot of genes, it can be poorly readable, especially within the notebook. This function does not include a parameter to restrict to ‘core’ genes like the cnetplot or ridgeplot so even when plotting a very small number of enriched terms, the X-axis is too cramped to be readable.

enrichplot::heatplot(gsea_kegg, showCategory = 3)

The gseKEGG result contains a data column containing the core/leading edge genes:

colnames(gsea_kegg@result)

So we can extract those and provide these to the heatmap to try to improve clarity. The downside of course is that we then lose some of the information about which genes are shared across terms.

# create an object that contains the core genes grouped by term and their fold change values
ranked_core_genes <- ranked[names(ranked) %in% gsea_kegg$core_enrichment]

# Plot with just core genes
enrichplot::heatplot(gsea_kegg, foldChange = ranked_core_genes, showCategory = 3)

Saving to a file only improves things slightly.

png("clusterprofiler_gseKEGG_heatplot.png", width = 11.7, height = 8.3, units = "in", res = 300)
# enrichplot::heatplot(gsea_kegg, showCategory = 3)
enrichplot::heatplot(gsea_kegg, foldChange = ranked_core_genes, showCategory = 3)
dev.off()

16.8.2.7 ridgeplot

This function is specific for GSEA. It does not apply to ORA, since it plots NES.

The height of the ridge in the ridge plot represents the density of genes from the gene set at different points in the ranked list.

The area under the curve represents the distribution of these genes across the ranked list.

Higher peaks and more concentrated areas under the curve at the top of the list indicate stronger enrichment.

enrichplot::ridgeplot(gsea_kegg,
  showCategory = 15,
  fill = "p.adjust",
  core_enrichment = TRUE,
  label_format = 30,
  orderBy = "NES",
  decreasing = FALSE
)

16.8.2.8 gseaplot

Another GSEA-specific plot, which shows the contribution to the normalised enrichment score (NES) of genes across the gene list.

The gseaplot can plot one term at a time. You would run GSEA plots for key terms of interest in your results.

The term ‘hsa04020’ (Calcium signaling pathway) has a positive NES, an adjusted P value of 0.00046 and a set size of 279.

gene_set_id <- "hsa04020"

enrichplot::gseaplot(gsea_kegg,
  gene_set_id,
  by = "all",
  title = paste("GSEA KEGG result:", gene_set_id),
  color = "black",
  color.line = "green",
  color.vline = "#FA5860"
)

The term ‘hsa04062’ (Chemokine signaling pathway) has a negative NES, a highly significant adjusted P value of 1.11E-06 and a set size of 129.

gene_set_id <- "hsa04062"

enrichplot::gseaplot(gsea_kegg,
  gene_set_id,
  by = "all",
  title = paste("GSEA KEGG result:", gene_set_id),
  color = "black",
  color.line = "green",
  color.vline = "#FA5860"
)

The term ‘hsa05163’ has a slightly negative NES and a borderline adjusted P value of 0.046.

gene_set_id <- "hsa05163"

enrichplot::gseaplot(gsea_kegg,
  gene_set_id,
  by = "all",
  title = paste("GSEA KEGG result:", gene_set_id),
  color = "black",
  color.line = "green",
  color.vline = "#FA5860"
)

16.8.2.9 volcano plot

The volplot function within enrichplot does not support GSEA result objects, but we can use ggplot2 for this.

ggplot2 is a highly flexible R package for visualisations, and unlike enrichplot is not constrained to a specific plot type.

The volcano plot has the advantage of showing whether the genes in the enriched terms were predominantly from the top (upregulated) or bottom (downregulated) end of the list.

ggplot2::ggplot(gsea_kegg@result, aes(x = enrichmentScore, y = -log10(p.adjust), color = p.adjust)) +
  geom_point(alpha = 0.7, size = 2) + # Adjust point size
  scale_color_gradient(low = "blue", high = "red") + # Color by p.adjust values
  theme_minimal() +
  labs(
    title = "GSEA KEGG",
    x = "Enrichment Score (NES)",
    y = "-log10(Adjusted P-value)",
    color = "Adjusted P-value"
  ) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + # Rotate x-axis labels for readability
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") + # Add vertical line at x=0
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black") + # Add horizontal line at p=0.05 cutoff
  geom_text(aes(label = Description),
    hjust = 0.5,
    vjust = -0.5, # Move labels higher off the points
    size = 3,
    check_overlap = TRUE,
    alpha = 0.7
  ) # Add labels for each pathway term

16.8.3 R version and R package versions

We can also capture the version of R and other session details including all loaded packages and versions with the sessionInfo() function:

16.9 End of activity summary

  • We have explored the supported organisms, namespaces and databases of the clusterProfiler enrichment functions
  • We have extracted a ranked gene list for GSEA and converted the gene IDs for compatibility with gseKEGG
  • We have performed GSEA on the KEGG database with gseKEGG and visualised the results with multiple plot types
  • We have captured all version details relevant to the session within the R notebook and knit the file to HTML for record keeping