irGSEA: Integrated Single-Cell Rank-Based Gene Set Enrichment Analysis

To install the core irGSEA package, you must first install its dependencies from CRAN and Bioconductor, then install irGSEA from GitHub.

1. Install CRAN dependencies

cran.packages <- c("aplot", "BiocManager", "circlize", "cowplot","data.table", 
                   "devtools", "doParallel", "doRNG", "dplyr", "ggfun", "gghalves", 
                   "ggplot2", "ggplotify", "ggridges", "ggsci", "irlba",
                   "magrittr", "Matrix", "msigdbr", "pagoda2", "plyr", "pointr", 
                   "purrr", "RcppML", "readr", "reshape2", "reticulate", 
                   "rlang", "RMTstat", "RobustRankAggreg", "roxygen2", 
                   "Seurat", "SeuratObject", "stringr", "tibble", "tidyr", 
                   "tidyselect", "tidytree", "VAM")

for (i in cran.packages) {
  if (!requireNamespace(i, quietly = TRUE)) {
    install.packages(i, ask = F, update = F)
  }
}

2. Install Bioconductor dependencies

bioconductor.packages <- c("AUCell", "BiocParallel", "ComplexHeatmap", 
                           "decoupleR", "fgsea", "ggtree", "GSEABase", 
                           "GSVA", "Nebulosa", "scde", "singscore",
                           "SummarizedExperiment", "UCell",
                           "viper","sparseMatrixStats")

for (i in bioconductor.packages) {
  if (!requireNamespace(i, quietly = TRUE)) {
    BiocManager::install(i, ask = F, update = F)
  }
}

3. Install irGSEA from GitHub

if (!requireNamespace("irGSEA", quietly = TRUE)) { 
    devtools::install_github("chuiqin/irGSEA", force =T)
}

# install packages from CRAN
cran.packages <- c("aplot", "BiocManager", "circlize", "cowplot","data.table", 
                   "devtools", "doParallel", "doRNG", "dplyr", "ggfun", "gghalves", 
                   "ggplot2", "ggplotify", "ggridges", "ggsci", "irlba",
                   "magrittr", "Matrix", "msigdbr", "pagoda2", "plyr", "pointr", 
                   "purrr", "RcppML", "readr", "reshape2", "reticulate", 
                   "rlang", "RMTstat", "RobustRankAggreg", "roxygen2", 
                   "Seurat", "SeuratObject", "stringr", "tibble", "tidyr", 
                   "tidyselect", "tidytree", "VAM")

for (i in cran.packages) {
  if (!requireNamespace(i, quietly = TRUE)) {
    install.packages(i, ask = F, update = F)
  }
}

# install packages from Bioconductor
bioconductor.packages <- c("AUCell", "BiocParallel", "ComplexHeatmap", 
                           "decoupleR", "fgsea", "ggtree", "GSEABase", 
                           "GSVA", "Nebulosa", "scde", "singscore",
                           "SummarizedExperiment", "UCell",
                           "viper","sparseMatrixStats")

for (i in bioconductor.packages) {
  if (!requireNamespace(i, quietly = TRUE)) {
    BiocManager::install(i, ask = F, update = F)
  }
}

# install packages from Github
if (!requireNamespace("irGSEA", quietly = TRUE)) { 
    devtools::install_github("chuiqin/irGSEA", force =T)
}

To use gene sets from clusterProfiler (like KEGG or GO) with irGSEA, you must convert the clusterProfiler output into a named list of gene symbols.

Workflow for KEGG:

1. Download KEGG pathways using clusterProfiler::gson_KEGG(species = "hsa").
2. Convert the resulting object to a data frame.
3. Convert Entrez IDs to Gene Symbols using clusterProfiler::bitr() and org.Hs.eg.db.
4. Clean the data (remove NAs).
5. Group the data by pathway name and create a named list where each element is a character vector of unique gene symbols.

Workflow for GO:

1. Download GO terms using clusterProfiler::gson_GO(OrgDb = "org.Hs.eg.db", ont = "BP").
2. Convert to a data frame and map Entrez IDs to Symbols.
3. Group by term name and create a named list of unique gene symbols.

# Example: Converting KEGG to irGSEA compatible list
library(clusterProfiler)
library(tidyverse)

# 1. Get KEGG data
kk <- clusterProfiler::gson_KEGG(species = "hsa")

# 2. Convert to data frame
kegg.list <- dplyr::left_join(kk@gsid2name, kk@gsid2gene, by = "gsid")

# 3. Convert ID to Symbol
gene_name <- clusterProfiler::bitr(kegg.list$gene, fromType = "ENTREZID", toType = "SYMBOL", OrgDb = "org.Hs.eg.db")
kegg.list <- dplyr::full_join(kegg.list, gene_name, by = c("gene"="ENTREZID"))
kegg.list <- kegg.list[complete.cases(kegg.list[, c("gene", "SYMBOL")]), ]

# 4. Format as named list for irGSEA
kegg.list$name <- factor(kegg.list$name)
kegg.list <- kegg.list %>% 
  dplyr::group_split(name, .keep = F) %>%
  purrr::map( ~.x %>% dplyr::pull(SYMBOL) %>% unique(.)) %>%
  purrr::set_names(levels(kegg.list$name))

Use the following functions to generate global visualizations of gene set enrichment across clusters. These functions typically take an enrichment result object (e.g., result.dge) and a method (such as "RRA").

• irGSEA.heatmap(): Generates a heatmap showing co-upregulated or co-downregulated gene sets per cluster.
• irGSEA.bubble(): Generates a bubble plot for co-upregulated or co-downregulated gene sets. Note: If you encounter an error regarding caller_env, uninstall ggtree and reinstall it from GitHub using remotes::install_github("YuLab-SMU/ggtree").
• irGSEA.upset(): Generates an UpSet plot showing intersections of significant gene sets among clusters.
• irGSEA.barplot(): Generates a stacked bar plot showing intersections of significant gene sets across multiple methods (e.g., "AUCell", "UCell", "singscore", "ssgsea", "JASMINE", "viper", "RRA").

# Heatmap
irGSEA.heatmap(object = result.dge, method = "RRA", top = 50)

# Bubble plot
irGSEA.bubble(object = result.dge, method = "RRA", top = 50)

# Upset plot
irGSEA.upset(object = result.dge, method = "RRA")

# Stacked bar plot
irGSEA.barplot(object = result.dge, method = c("AUCell", "UCell", "singscore", "ssgsea", "JASMINE", "viper", "RRA"))

Use these functions to inspect the expression and distribution of specific gene sets within a single enrichment method on a specific object (e.g., a Seurat object).

• irGSEA.density.scatterplot(): Shows expression and distribution on a reduction plot (e.g., reduction = "umap").
• irGSEA.halfvlnplot(): Shows expression and distribution among clusters using a half violin plot.
• irGSEA.vlnplot(): Compares expression and distribution of a gene set across multiple methods (e.g., method = c("AUCell", "UCell", ...)).
• irGSEA.ridgeplot(): Shows expression and distribution among clusters using a ridge plot.
• irGSEA.densityheatmap(): Shows expression and distribution among clusters using a density heatmap.

# Density scatterplot on UMAP
irGSEA.density.scatterplot(object = pbmc3k.final, method = "UCell", show.geneset = "HALLMARK-INFLAMMATORY-RESPONSE", reduction = "umap")

# Half violin plot
irGSEA.halfvlnplot(object = pbmc3k.final, method = "UCell", show.geneset = "HALLMARK-INFLAMMATORY-RESPONSE")

# Violin plot comparing methods
irGSEA.vlnplot(object = pbmc3k.final, method = c("AUCell", "UCell", "singscore", "ssgsea", "JASMINE", "viper"), show.geneset = "HALLMARK-INFLAMMATORY-RESPONSE")

# Ridge plot
irGSEA.ridgeplot(object = pbmc3k.final, method = "UCell", show.geneset = "HALLMARK-INFLAMMATORY-RESPONSE")

# Density heatmap
irGSEA.densityheatmap(object = pbmc3k.final, method = "UCell", show.geneset = "HALLMARK-INFLAMMATORY-RESPONSE")

The irGSEA package requires gene sets to be provided as a named list of gene symbols. Since the latest MSigDB releases are often distributed as .db.zip (SQLite) or .gmt files, you must convert them into the appropriate list format.

There are two primary methods to prepare these sets:

Method 1: Using the MSigDB SQLite database (db.zip)

This method involves downloading the database, connecting via RSQLite, joining several tables (like gene_set, gene_set_details, and gene_symbol) to map IDs to symbols, and then splitting the resulting tibble into a named list.

Method 2: Using the MSigDB GMT file

This is a simpler approach where you download the .gmt file and use clusterProfiler::read.gmt() to load it. You then filter the terms (e.g., using HALLMARK_, GOBP_, or KEGG_ prefixes) and convert them into a named list.

Common Collection Filters

When preparing lists for irGSEA, you typically filter for specific collections:

• Hallmark: Filter for collection == "H" or terms containing HALLMARK_.
• GO Biological Process: Filter for collection == "C5:GO:BP" or terms containing GOBP_.
• KEGG: Filter for collection == "C2:CP:KEGG_MEDICUS" or terms containing KEGG_.

# Example: Converting a GMT file to the list format required by irGSEA
library(clusterProfiler)
library(tidyverse)

# 1. Download and read GMT
gmt_url <- "https://data.broadinstitute.org/gsea-msigdb/msigdb/release/2023.2.Hs/msigdb.v2023.2.Hs.symbols.gmt"
local_gmt <- "./msigdb.v2023.2.Hs.symbols.gmt"
download.file(gmt_url, local_gmt)
msigdb <- clusterProfiler::read.gmt(local_gmt)

# 2. Convert Hallmark sets to a named list
msigdb.h <- msigdb %>% 
  dplyr::filter(str_detect(term, pattern = regex("HALLMARK_", ignore_case=TRUE))) 
msigdb.h$term <- factor(msigdb.h$term)
msigdb.h <- msigdb.h %>% 
  dplyr::group_split(term, .keep = F) %>%
  purrr::map( ~.x %>% dplyr::pull(gene) %>% unique(.)) %>%
  purrr::set_names(levels(msigdb.h$term))

If you intend to use specific methods like VISION, gficf, or Sargent, you need to install additional packages.

Requirements for specific methods:

• VISION: Requires VISION from GitHub.
• mdt: Requires ranger from GitHub.
• gficf: Requires RcppML (version > 0.3.7) from GitHub and then library(RcppML) before running gficf.
• GSVApy / ssGSEApy: Requires SeuratDisk from GitHub.
• sargent: Requires sargent from Gitee.
• pagoda2: Requires scde from GitHub.

Python Environment Setup

irGSEA uses reticulate to manage a Conda environment. If you encounter errors finding the conda binary, run reticulate::install_miniconda().

To ensure Python dependencies are available, create the irGSEA environment and install the following via pip:

• anndata
• scanpy
• argparse
• gseapy
• decoupler

# Create conda env
if (! "irGSEA" %in% reticulate::conda_list()$name) {
      reticulate::conda_create("irGSEA")
}

# Install python packages
require.package <- c("anndata", "scanpy", "argparse", "gseapy", "decoupler")
for (i in require.package) {
  reticulate::conda_install(envname = "irGSEA", packages = i, pip = T)
}

# VISION
if (!requireNamespace("VISION", quietly = TRUE)) { 
    devtools::install_github("YosefLab/VISION", force =T)
}

# mdt need ranger
if (!requireNamespace("ranger", quietly = TRUE)) { 
  devtools::install_github("imbs-hl/ranger", force =T)
}

# gficf need RcppML (version > 0.3.7) package
if (utils::packageVersion("RcppML") > "0.3.7") {
  message("The version of RcppML should greater than 0.3.7 and install RcppML package from Github")
  devtools::install_github("zdebruine/RcppML", force =T)
}

# please first `library(RcppML)` if you want to perform gficf
if (!requireNamespace("gficf", quietly = TRUE)) { 
    devtools::install_github("gambalab/gficf", force =T)
}

# GSVApy and ssGSEApy need SeuratDisk package
if (!requireNamespace("SeuratDisk", quietly = TRUE)) {
    devtools::install_github("mojaveazure/seurat-disk", force =T)
}

# sargent
if (!requireNamespace("sargent", quietly = TRUE)) {
    devtools::install_git("https://gitee.com/fan_chuiqin/PMCB-Sargent.git", force =T)
}

# pagoda2 need scde package
if (!requireNamespace("scde", quietly = TRUE)) {
    devtools::install_github("hms-dbmi/scde", force =T)
}

#### create conda env
if (! "irGSEA" %in% reticulate::conda_list()$name) {
      reticulate::conda_create("irGSEA")
}

# if python package exist
python.package <- reticulate::py_list_packages(envname = "irGSEA")$package
require.package <- c("anndata", "scanpy", "argparse", "gseapy", "decoupler")
for (i in seq_along(require.package)) {
  if (i %in% python.package) {
    reticulate::conda_install(envname = "irGSEA", packages = i, pip = T)
  }
}

The irGSEA.score() function calculates enrichment scores for multiple gene sets using various methods and returns a Seurat object containing the resulting score matrices as new assays.

Important Considerations:

• Performance: Methods like AUCell or ssgsea can be slow with large numbers of genes or cells. It is recommended to use high-quality subsets.
• Parallelization: If you encounter an error Valid ‘mctype’: ‘snow’ or ‘doMC’ when ncores > 1, ensure AUCell is version >= 1.14 or set ncores = 1.
• Warnings: You may see warnings regarding constant expression values, small gene sets, or experimental dgCMatrix usage; these can generally be ignored.
• Dependencies: To access the full MSigDB dataset, ensure the msigdbdf package is installed.

Parameters:

• object: A Seurat object.
• assay: The assay to use (e.g., "RNA").
• slot: The data slot (e.g., "data").
• method: A vector of methods to use. Supported: "AUCell", "UCell", "singscore", "ssgsea", "JASMINE", "viper".
• species: The organism (e.g., "Homo sapiens").
• geneid: The type of gene identifier (e.g., "symbol").
• msigdb: Boolean, whether to use MSigDB (default TRUE).

pbmc3k.final <- irGSEA.score(object = pbmc3k.final, assay = "RNA", 
                             slot = "data", seeds = 123, ncores = 4,
                             min.cells = 3, min.feature = 0,
                             custom = F, geneset = NULL, msigdb = T, 
                             species = "Homo sapiens", category = "H",  
                             subcategory = NULL, geneid = "symbol",
                             method = c("AUCell", "UCell", "singscore", 
                                        "ssgsea", "JASMINE", "viper"),
                             aucell.MaxRank = NULL, ucell.MaxRank = NULL, 
                             kcdf = 'Gaussian')

The irGSEA.hub() function calculates hub genes for a specified gene set based on the correlation between the gene set's score and the expression (or rank) of the genes included in that set.

Parameters:

• object: The input object.
• assay: The assay to use (e.g., "RNA").
• slot: The slot to use (e.g., "data").
• method: A vector of enrichment methods.
• show.geneset: A vector of gene set names.
• type: Either "expression" or "rank" (default is "rank").
• ncores: Number of cores to use.
• top: Number of top genes to return.
• correlation.color: Colors for the correlation results.

Returns a list containing $hub_result (a data frame with method, geneset, gene, correlation, and p.value) and $hub_plot (a list of plots for each gene set).

hub.result <- irGSEA.hub(object = pbmc3k.final, assay = "RNA", slot = "data",
                         method = c("AUCell","UCell","singscore", "ssgsea",
                                     "JASMINE", "viper"),
                         show.geneset = c("HALLMARK-INFLAMMATORY-RESPONSE",
                                          "HALLMARK-APOPTOSIS"),
                         ncores = 4, type = "rank", maxRank = 2000, top = 5,
                         correlation.color = c("#0073c2","white","#efc000"),
                         method.color = NULL)

# Access results
head(hub.result$hub_result)
hub.result$hub_plot$`HALLMARK-APOPTOSIS`

The irGSEA.integrate() function identifies common differentially expressed gene sets across multiple enrichment analysis methods.

How it works:

1. It performs a Wilcoxon test on all enrichment score matrices and gene sets.
2. It selects gene sets with an adjusted p-value < 0.05.
3. It uses RRA (Robust Rank Aggregation) to integrate these results.
4. It identifies gene sets that are statistically significant and common across all selected enrichment methods.

Returns: A list containing the integrated results.

Parameters:

• object: A Seurat object containing enrichment score assays.
• group.by: The metadata column used for grouping (e.g., cell type annotations).
• method: A vector of the methods used during the scoring step (e.g., c("AUCell", "UCell", ...)).

result.dge <- irGSEA.integrate(object = pbmc3k.final, 
                               group.by = "seurat_annotations",
                               metadata = NULL, col.name = NULL,
                               method = c("AUCell","UCell","singscore",
                                          "ssgsea", "JASMINE", "viper"))

class(result.dge) # Returns [1] "list"