clusterProfiler 4.0: A universal enrichment tool for interpreting omics data

Kyoto encyclopedia of genes and genomes

KEGG is an encyclopedia of genes and genomes.¹⁹ Molecular functions are represented by networks of interactions and reactions mainly in the form of KEGG pathways and modules. A KEGG module is a collection of manually defined function units. In some situations, KEGG modules have a more straightforward interpretation. Both KEGG pathways and KEGG modules are supported by clusterProfiler. Many software tools that support KEGG analysis have stopped updating since July 2011 when KEGG initiated an academic subscription model for FTP downloading. These tools use relatively old KEGG data, and the result might be inaccurate and misleading. Fortunately, the KEGG web resource is freely available. The clusterProfiler package does not pack any KEGG data. Instead, it queries the latest online KEGG database through web API to perform functional analysis. The advantage of this feature is obvious: it allows clusterProfiler to use up-to-date data and support all the species that have KEGG annotation (more than 6,000 species are listed in http://www.genome.jp/kegg/catalog/org_list.html). Moreover, clusterProfiler supports the KEGG Orthology database and can be used to perform functional characterization of the microbiomes.²⁰

In the following example, GSEA was performed with KEGG pathway. Figure 2A shows the plotting of GSEA enrichment results to visualize the top five perturbed pathways, i.e., the top five highest absolute values of the normalized enrichment score (NES).⁹ The NES indicates the shift of genes belonging to a certain pathway toward either end of the ranked list and represents pathway activation or suppression. To further explore the pathway crosstalk effects, we visualized gene expression distribution of core enrichment genes using an UpSet plot (Figure 2B). The result shows that the expression values of genes in the intersection of cell-cycle and DNA-replication pathways are higher than those uniquely belonging to either of the two pathways. These overlapping genes are mainly minichromosome maintenance (MCM) genes, which can potentially serve as biomarkers for tumor diagnosis.²¹ The intersection of the interleukin-17 (IL-17) signaling pathway and the proteasome pathway is only associated with one gene, interferon-γ (IFN-γ). The IL-17 signaling pathway induces an inflammatory response,²² while IFN-γ regulates proteasome formation.²³ These effects ultimately reshape the tumor microenvironment.

Open in a separate window

Figure 2

KEGG pathway enrichment analysis

In GSEA enrichment plot (A), the curves represent the running sum of the enrichment scores, the middle part of the graph shows the position of genes that are related to certain pathways, and the bottom part of the graph displays how the metric (e.g., fold change) is distributed along with the list. The UpSet plot (B) visualizes the metric distribution of core enrichment genes. It differentiates genes that uniquely belong to a pathway or are associated with two or more pathways.

Universal interface for biomedical gene sets

With the advancement of the sequencing technology, the investigation into functions for transcriptomes from non-model organisms is increasingly demanded. However, most tools in this field are designed for GO and KEGG analyses with support limited to one or several model organisms. Besides, there are increasingly more biological knowledge databases available for exploring functional characteristics from different perspectives, such as Disease Ontology,²⁴ Reactome Pathway,²⁵ Medical Subject Headings,²⁶, and WikiPathway.²⁷ There is an urgent need for integration and support of these databases. To address these issues, clusterProfiler provides two general functions, enricher and GSEA for ORA and GSEA, with user-provided gene annotations. These two functions allow the application of all ontologies or pathways curated in diverse databases as the background in customized analyses. Therefore, users could easily import external annotations (e.g., electronic annotations using Blast2GO²⁸ and KAAS²⁹ for GO and KEGG annotations, respectively) for newly sequenced species. Moreover, it is convenient to perform functional analysis using up-to-date annotations from all popular databases, such as InterPro, Clusters of Orthologous Groups, and Mouse Phenotype Ontology, to name a few, without waiting for the updates of other tools. It would be suitable for the timely analysis of gene sets with emerging interests, such as human cell markers³⁰ and COVID-19-related gene sets.

The gene set annotation required by enricher and GSEA is a two-column data frame with one column representing gene set names (ID or descriptive name) and the other showing the corresponding genes. The gene matrix transposed (GMT) format is widely used to distribute gene set annotations. There are many gene set libraries available online (e.g., https://maayanlab.cloud/Enrichr/#stats), including MSigDB (Molecular Signatures Database), Disease Signatures, and CCLE (Cancer Cell Line Encyclopedia). To enable the utilization of these gene sets in clusterProfiler as the background annotation to explore the underlying biological mechanisms, clusterProfiler provides a parser function, read.gmt, to import GMT files that can be directly passed to the enricher and GSEA functions. In the following example, we used the GSEA function to perform gene set enrichment analysis using WikiPathways (Figure 5B). The annotation data were parsed by using read.gmt.wp, which is a customized version of read.gmt for importing GMT files from WikiPathways.

Open in a separate window

Figure 5

Visualizing enrichment results using ggplot2

A lollipop chart to visualize the rich factors from ORA (A) and a bar chart to visualize normalized enrichment scores from GSEA (B).

Functional interpretation of genomic ROIs

With the increasing availability of genomic sequences, non-coding genomic regions (e.g., cis-regulatory elements, non-coding RNAs, and transposons) have posed a demanding challenge to exploration of their roles in various biological processes.¹ Unlike coding genes, non-coding genomic regions are typically not well functionally annotated. Analyzing biological functions of the proximal genes is a common strategy in research on the biological meaning of a set of non-coding genomic regions. Software tools, such as the Genomic Regions Enrichment of Annotations Tool (GREAT),³¹ are implemented to follow this strategy. However, these tools only support a limited number of species. For example, GREAT is designed for human and mouse only. In addition, many tools only take the host or nearest genes into consideration but ignore long-distance regulations. Our in-house developed package, ChIPseeker,¹⁰ is originally designed for chromatin immunoprecipitation (ChIP) peak annotation, comparison, and visualization and has been employed to analyze genome-wide ROIs, such as open chromatin regions obtained by DNase-seq³² and ATAC-seq.³³ To facilitate biological interpretation of genome-wide regions, we implemented a function, seq2gene, in ChIPseeker to associate genomic regions with coding genes through many-to-many mapping. It automatically maps genomic regions to host genes (either located in exon or intron), proximal genes (located in the promoter region), and flanking genes (located upstream and downstream within user-specified distance). The seq2gene function supports a wide variety of species if a genomic annotation, such as the TxDb (UCSC-based) or EnsDb (Ensembl-based) object, is available. After mapping genomic regions to coding genes, clusterProfiler can be employed to perform functional enrichment analysis of the coding genes to assign biological meanings to the set of genomic regions. The combination of ChIPseeker and clusterProfiler allows more biological ontology or pathway databases to be utilized to explore functions of genomic regions for a wide variety of species.

A dataset of ChIP-seq with antibody against CBX6 (GEO: GSM1295076) was used in the above example. The genomic binding regions were mapped to coding genes using the seq2gene function with UCSC genomic annotation. The Entrez gene IDs were converted into gene symbols using the bitr function implemented in clusterProfiler. To identify and characterize transcript cofactors, we performed functional enrichment analysis using the ENCODE and ChEA transcript factor gene sets. The result was visualized as a category-gene network (Figure 3), which showed that genes associated with CBX6 (obtained by the seq2gene function) significantly overlap with genes regulated by POU5F1, TRIM28, SUZ12, and EZH2. OCT4 (POU5F1)³⁴ and KAP1 (TRIM28)³⁵ have been reported to interact with polycomb repressive complex 1 (PRC1), and CBX6 is a known subunit of PRC1.³⁶ SUZ12 and EZH2 are core components of PRC2 and negatively regulate CBX6.³⁷ These pieces of evidence support the effectiveness of these analyses including the mapping of genomic ROIs to coding genes and functional enrichment, which suggest that this method can be used to identify unknown cofactors (Figure 3) and characterize functions of genomic regions.

Open in a separate window

Figure 3

Functional enrichment analysis of genomic regions of interest

Genomic regions are linked to coding genes, which are then used to identify transcript cofactors by testing significant overlap of target genes.

Comparison among different conditions

The clusterProfiler library is designed to allow the comparison of functional enrichment results from multiple experimental conditions or multiple time points. With an input of a collection of gene lists, the compareCluster function applies a function (e.g., enricher) with user settings to perform functional enrichment analysis for each of the gene lists and aggregates the results into a single object. Thus, enrichment results of multiple groups are easily explored and plotted together for comparison with a user-friendly interface. Comparing functional profiles can reveal functional consensus and differences among different experiments and helps in identifying differential functional modules in omics datasets. In the updated version, compareCluster provides a new interface supporting a formula that is widely used in R for specifying statistical models; this allows more complicated experimental designs to be supported (e.g., time-course experiment with different treatments). With the infrastructure of clusterProfiler to support a wide range of ontology and pathway annotations and multiple organisms, the comparison can be applied to many circumstances.

The dataset, DE_GSE8057, was derived from the GEO: GSE8057 dataset in the GEO database. The GSE8057 dataset contains expression data from ovarian cancer cells at multiple time points (0, 2, 6, and 24 h) and under two treatment conditions (cisplatin and oxaliplatin).³⁸ The DE_GSE8057 dataset contains DEGs obtained from different treatments and time points versus control samples. Eight groups of DEG lists (specified by the formula Gene ~ Time + Treatment) were analyzed simultaneously using compareCluster with WikiPathways. The result (Figure 4) indicates that the two drugs have distinct effects at the beginning but consistent effects in the later stages. Several pathways including DNA damage and cell-cycle progression were perturbed by either cisplatin or oxaliplatin drug exposure. The finding is consistent with the discovery obtained by data-driven modeling.³⁸

Open in a separate window

Figure 4

Comparing functional profiles among different levels of conditions

The compareCluster function performed enrichment analysis simultaneously for eight lists of DEGs. The results were visualized as a dot plot with an x axis representing one level of conditions (time course) and a facet panel indicating another level of conditions (drug treatments).

Data frame interface for accessing enriched results

The outputs of ORA and GSEA are enrichResult and gseaResult objects, respectively, while the output of compareCluster is a compareClusterResult object. These S4 objects contain input data, analysis settings, and enriched results, which allow more informative data to be available for downstream interpretation and visualization. To enable easy access to the enriched result, clusterProfiler implements as.data.frame methods to convert the S4 objects to data frames that can be easily exported as CSV files. In addition, clusterProfiler provides a data frame interface that mimics data frame operations to access rows, columns, and subsets of rows and columns from the S4 objects of the enriched result. Users can use head and tail to print part of the result. The nrow, ncol, and dim methods are also supported to access basic information such as how many pathways are enriched.

The [ and $ operators for subsetting are also supported. We redefined the [[ operator to help users access which genes are annotated by a selected ontology or pathway. For GSEA output, the [[ operator will return core enriched genes (i.e., genes in the leading edge) of the selected gene set.

Tidy interface for data operation

To facilitate data manipulation and exploration of the enrichment result, clusterProfiler extends the dplyr verbs to support enrichResult, gseaResult, and compareClusterResult objects. Following the concept of tidiness, these verbs provide robust and standardized operations for data transformation and can be assembled into a workflow using the pipe operator (%>%). This allows users to explore the results effectively and develop reproducible and human-readable pipelines. For example, it allows the filtering of enriched results using different criteria (e.g., adjusted p values less than 0.001, and the number of input genes annotated to the enriched term should be greater than 10).

For ORA results, clusterProfiler provides geneRatio (ratio of input genes that are annotated in a term) and BgRatio (ratio of all genes that are annotated in this term). However, other concepts are widely used to help in interpreting enrichment results, such as the rich factor and fold enrichment. A rich factor is defined as the ratio of input genes (e.g., DEGs) that are annotated in a term to all genes that are annotated in this term. The fold enrichment is defined as the ratio of the frequency of input genes annotated in a term to the frequency of all genes annotated to that term, and it is easy to calculate by dividing geneRatio by BgRatio. Here, as an example, we used the mutate verb to create a new column of richFactor based on information available in the clusterProfiler output.

The following example uses the GSEA enrichment result generated in the previous session. The result was sorted by absolute values of NESs using the arrange verb. NES is an indicator to interpret the degree of enrichment. A positive NES indicates that members of the gene set tend to appear at the top of the rank (pathway activation), and a negative NES indicates the opposite circumstance (pathway suppression). We used the group_by verb to group the result based on the sign of NES, and the slice verb was used to extract the first five enriched pathways for each group (i.e., five activated pathways that have the largest NES values and five suppressed pathways that have the smallest NES values). These verbs return the same object type as their input and do not affect downstream analysis and visualization.

Visualization using Ggplot2

The enrichplot package is originally derived from DOSE and clusterProfiler packages and serves as a de facto visualization tool for visualizing enrichment results for outputs from clusterProfiler as well as DOSE, ReactomePA, and meshes. These methods allow users without programming skills to generate effective visualization to explore and interpret results. All the visualization methods implemented are based on ggplot2, which allows customization using the grammar of graphics. Moreover, we also extend ggplot2 to support enrichment results so that users can use the ggplot2 syntax directly to visualize enrichment results. The following example demonstrates the application of ggplot2 grammar of graphics to visualize the GO enrichment result (ORA) as a lollipop chart using the rich factor that was generated in the previous session using the dplyr verbs (Figure 5A).

In Figure 5B, the most significant activated and suppressed pathways (GSEA) selected by a series of dplyr verb operations are visualized as a bar chart using the ggplot2 syntax. Although visualization methods used to generate Figure 5 are not provided in clusterProfiler, it is easy to generate such graphs using the tidy interface and ggplot2. The combination of the tidy interface for data wrangling and the support of ggplot2 for visualization creates many possibilities for users to explore and visualize enrichment results using consistent grammar. It allows novel exploratory approaches to reveal unanticipated mechanisms as well as prototyping new visualization methods.

This work was supported by a startup fund from Southern Medical University.