Overview

This package performs gene set enrichment analyses including functional/pathway enrichment and fold change GSEA. Checkout BIGpicture for ways to plot these results.

Setup

Software and packages

This pipeline should be completed in R version 4.5.0 (2025-04-11) or newer. Please also download the following packages.

# CRAN packages
install.packages(c("tidyverse"))

# Bioconductor packages
## These are dependencies for kimma that are beneficial to download separately
install.packages("BiocManager")
BiocManager::install(c("limma", "edgeR"))

# GitHub packages
install.packages("devtools")
devtools::install_github("BIGslu/kimma")
devtools::install_github("BIGslu/SEARchways")

And load them into your current R session.

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(kimma)
library(SEARchways)

To ensure reproducibility with this document, set the random seed value.

set.seed(651)

Data background

Example data were obtained from media controls and human rhinovirus (HRV) infected human plasmacytoid dendritic cells. You can learn more about these data in

Dill-McFarland KA, Schwartz JT, Zhao H, Shao B, Fulkerson PC, Altman MC, Gill MA. 2022. Eosinophil-mediated suppression and Anti-IL-5 enhancement of plasmacytoid dendritic cell interferon responses in asthma. J Allergy Clin Immunol. Epub ahead of print. doi: 10.1016/j.jaci.2022.03.025. — GitHub

Specifically, this vignette uses bulk human RNA-seq processed from fastq to counts using SEAsnake and counts to voom using edgeR/limma. This results in voom-normalized, log2 counts per million (CPM) gene expression and associated sample and gene metadata. In total, 6 donors and 1000 random genes were selected for this vignette. Actual donor identifiers have been altered for privacy.

Load data

Gene fold change estimates.

example.gene.list    <- SEARchways::example.gene.list
str(example.gene.list)
## List of 2
##  $ HRV1: Named num [1:100] 0.1636 0.26787 0.53233 -0.09021 0.00173 ...
##   ..- attr(*, "names")= chr [1:100] "ENSG00000198898" "ENSG00000121410" "ENSG00000092096" "ENSG00000076928" ...
##  $ HRV2: Named num [1:100] -4.663 -0.367 -0.876 0.898 -0.601 ...
##   ..- attr(*, "names")= chr [1:100] "ENSG00000174123" "ENSG00000165675" "ENSG00000103024" "ENSG00000151067" ...

Extract gene names only.

example.gene.list2 <- list(
  "HRV1"=names(example.gene.list$HRV1),
  "HRV2"=names(SEARchways::example.gene.list$HRV2))
str(example.gene.list2)
## List of 2
##  $ HRV1: chr [1:100] "ENSG00000198898" "ENSG00000121410" "ENSG00000092096" "ENSG00000076928" ...
##  $ HRV2: chr [1:100] "ENSG00000174123" "ENSG00000165675" "ENSG00000103024" "ENSG00000151067" ...

Create data frame formats as well.

example.gene.df <- data.frame(
  group = "HRV1",
  gene = names(example.gene.list$HRV1),
  logFC = example.gene.list$HRV1
)
head(example.gene.df)
##                 group            gene       logFC
## ENSG00000198898  HRV1 ENSG00000198898  0.16360114
## ENSG00000121410  HRV1 ENSG00000121410  0.26786954
## ENSG00000092096  HRV1 ENSG00000092096  0.53232684
## ENSG00000076928  HRV1 ENSG00000076928 -0.09020896
## ENSG00000171307  HRV1 ENSG00000171307  0.00172912
## ENSG00000116863  HRV1 ENSG00000116863  0.90057251
example.gene.df2 <- example.gene.df[,1:2]
head(example.gene.df2)
##                 group            gene
## ENSG00000198898  HRV1 ENSG00000198898
## ENSG00000121410  HRV1 ENSG00000121410
## ENSG00000092096  HRV1 ENSG00000092096
## ENSG00000076928  HRV1 ENSG00000076928
## ENSG00000171307  HRV1 ENSG00000171307
## ENSG00000116863  HRV1 ENSG00000116863

Shared parameters

All SEARchways functions except iterEnrich can take either a named list (gene_list) or data frame (gene_df) as input (see examples above). Gene names can be HGNC/MGI symbols, ENSEMBL, or ENTREZ ID. Use ID to specify based on your data.

Additionally, all SEARchways functions can utilize the Broad MSigDB gene set data base. Select the species, collection, and subcollection as needed. See function help pages for the lists of options.

Alternatively, you may also provide a custom data base such as WCGNA modules with db.

With the exception of data input, these parameters will not be demoed below since they work similarly across all functions.

Functional / pathway / hypergeometric enrichment

Fisher’s exact test of selected genes.

flexEnrich: Flexible enrichment

flexEnrich allows more customization than BIGprofiler but can run slower. By default, enrichment is performed against a protein-coding background and gene sets must be between 10 and 1E10 genes.

flexEnrich(
  gene_list = example.gene.list2,
  ID = "ENSEMBL", species = "human", 
  collection = "H")
## # A tibble: 51 × 11
##    group n_query_genes n_background_genes gs_collection pathway  n_pathway_genes
##    <chr>         <int>              <int> <chr>         <chr>              <int>
##  1 HRV1            100              21940 H             HALLMAR…             200
##  2 HRV1            100              21940 H             HALLMAR…             139
##  3 HRV1            100              21940 H             HALLMAR…             201
##  4 HRV1            100              21940 H             HALLMAR…             200
##  5 HRV1            100              21940 H             HALLMAR…             200
##  6 HRV1            100              21940 H             HALLMAR…             200
##  7 HRV1            100              21940 H             HALLMAR…             158
##  8 HRV1            100              21940 H             HALLMAR…             200
##  9 HRV1            100              21940 H             HALLMAR…             200
## 10 HRV1            100              21940 H             HALLMAR…             200
## # ℹ 41 more rows
## # ℹ 5 more variables: n_query_genes_in_pathway <int>, `k/K` <dbl>, pval <dbl>,
## #   FDR <dbl>, genes <list>

Use a data frame input.

flexEnrich(
  gene_df = example.gene.df2,
  ID = "ENSEMBL", species = "human", 
  collection = "H")

Change the background to a custom one (custom_bg) or include all genes by turning off protein_coding.

flexEnrich(
  gene_list = example.gene.list2,
  ID = "ENSEMBL", species = "human", 
  collection = "H",
  protein_coding = FALSE)

Restrict the gene sets included to those with at least 3 query genes and total size between 100 and 1000 genes.

flexEnrich(
  gene_list = example.gene.list2,
  ID = "ENSEMBL", species = "human", 
  collection = "H",
  minOverlap = 3,
  minGeneSetSize = 100,
  maxGeneSetSize = 1000)

Speed up computation by not saving the query genes found in each gene set.

flexEnrich(
  gene_list = example.gene.list2,
  ID = "ENSEMBL", species = "human", 
  collection = "H",
  print_genes = FALSE)

BIGprofiler: Less flexible but faster enrichment

BIGprofiler performs the same enrichment as flexEnrich but with fewer options. It does not allow you to alter the background, always utilizing protein-coding only. Thus, when this background is desired, BIGprofiler may be utilized since it runs faster than flexEnrich on large queries.

BIGprofiler(
  gene_list = example.gene.list2,
  ID = "ENSEMBL", species = "human", 
  collection = "H")
## [1] "HRV1"
## 
## [1] "HRV2"
## # A tibble: 51 × 13
##    group n_query_genes gs_collection gs_subcollection n_background_genes pathway
##    <chr>         <int> <chr>         <chr>                         <dbl> <chr>  
##  1 HRV1            100 H             ""                             4386 HALLMA…
##  2 HRV1            100 H             ""                             4386 HALLMA…
##  3 HRV1            100 H             ""                             4386 HALLMA…
##  4 HRV1            100 H             ""                             4386 HALLMA…
##  5 HRV1            100 H             ""                             4386 HALLMA…
##  6 HRV1            100 H             ""                             4386 HALLMA…
##  7 HRV1            100 H             ""                             4386 HALLMA…
##  8 HRV1            100 H             ""                             4386 HALLMA…
##  9 HRV1            100 H             ""                             4386 HALLMA…
## 10 HRV1            100 H             ""                             4386 HALLMA…
## # ℹ 41 more rows
## # ℹ 7 more variables: n_pathway_genes <dbl>, n_query_genes_in_pathway <dbl>,
## #   `k/K` <dbl>, pval <dbl>, FDR <dbl>, qvalue <dbl>, genes <list>

Restrict the gene sets included to those with total size between 100 and 1000 genes (similar to flexEnrich).

BIGprofiler(
  gene_list = example.gene.list2,
  ID = "ENSEMBL", species = "human", 
  collection = "H",
  minGSSize = 100, maxGSSize = 1000)

iterEnrich: Iterative enrichment

iterEnrich allows you to run enrichment on inputs other than genes. For example, if you have SNPs, each SNP may annotate to > 1 gene. Regular enrichment would result in over-representation of pathways associated these multiple annotation SNPs. Thus, iterEnrich subsamples gene lists to 1 annotation per feature (e.g. SNP, probe, etc) and runs enrichment across many subsamples to determine a more representative result.

Here, we create a fake gene annotation where each SNP is associated with 2 genes.

example.gene.anno <- data.frame(
  SNP = c(paste0("SNP",1:50), paste0("SNP",1:50)),
  gene = sample(example.gene.list2[["HRV1"]], 
                replace = TRUE, size = 100)
)
head(example.gene.anno)
##    SNP            gene
## 1 SNP1 ENSG00000076928
## 2 SNP2 ENSG00000213995
## 3 SNP3 ENSG00000073803
## 4 SNP4 ENSG00000184232
## 5 SNP5 ENSG00000107018
## 6 SNP6 ENSG00000105289

Run iterative enrichment. We reduce the number of iters and increase the processors to run quickly for this vignette.

iter <- iterEnrich(
  anno_df = example.gene.anno,
  ID = "ENSEMBL", species = "human", 
  collection = "H",
  niter = 10, ncores = 6)

# Summary results across iters
iter$summary
##      group                                  pathway   pval_min   pval_max
## 1  anno_df                     HALLMARK_COAGULATION 0.02908510 1.00000000
## 2  anno_df                      HALLMARK_COMPLEMENT 0.32085006 1.00000000
## 3  anno_df                     HALLMARK_E2F_TARGETS 0.30054682 0.31953545
## 4  anno_df         HALLMARK_ESTROGEN_RESPONSE_EARLY 0.30054682 0.31953545
## 5  anno_df          HALLMARK_ESTROGEN_RESPONSE_LATE 0.30054682 0.31953545
## 6  anno_df                  HALLMARK_G2M_CHECKPOINT 0.04915364 0.05611768
## 7  anno_df                      HALLMARK_GLYCOLYSIS 0.30054682 1.00000000
## 8  anno_df                 HALLMARK_HEME_METABOLISM 0.31953545 1.00000000
## 9  anno_df           HALLMARK_INFLAMMATORY_RESPONSE 0.30180159 0.32085006
## 10 anno_df               HALLMARK_KRAS_SIGNALING_UP 0.31953545 1.00000000
## 11 anno_df       HALLMARK_OXIDATIVE_PHOSPHORYLATION 0.04915364 0.31953545
## 12 anno_df                     HALLMARK_P53_PATHWAY 0.05611768 0.30054682
## 13 anno_df         HALLMARK_PI3K_AKT_MTOR_SIGNALING 0.17076862 0.18263882
## 14 anno_df HALLMARK_REACTIVE_OXYGEN_SPECIES_PATHWAY 0.08357430 0.08971189
## 15 anno_df                 HALLMARK_SPERMATOGENESIS 0.21409813 1.00000000
## 16 anno_df                  HALLMARK_UV_RESPONSE_DN 0.22666112 0.24181594
## 17 anno_df                  HALLMARK_UV_RESPONSE_UP 0.24581556 0.26202152
##    pval_median       FDR k_median   K  k/K_median        genes
## 1   0.02908510 0.3180002        2 139 0.014388489 ENSG0000....
## 2   0.32085006 0.3636301        1 201 0.004975124 ENSG0000....
## 3   0.31953545 0.3636301        1 200 0.005000000 ENSG0000....
## 4   0.31953545 0.3636301        1 200 0.005000000 ENSG0000....
## 5   0.31953545 0.3636301        1 200 0.005000000 ENSG0000....
## 6   0.05611768 0.3180002        2 200 0.010000000 ENSG0000....
## 7   1.00000000 1.0000000        0 200 0.000000000 ENSG0000....
## 8   0.31953545 0.3636301        1 200 0.005000000 ENSG0000....
## 9   0.32085006 0.3636301        1 201 0.004975124 ENSG0000....
## 10  0.31953545 0.3636301        1 200 0.005000000 ENSG0000....
## 11  0.31953545 0.3636301        1 200 0.005000000 ENSG0000....
## 12  0.05611768 0.3180002        2 200 0.010000000 ENSG0000....
## 13  0.18263882 0.3636301        1 105 0.009523810 ENSG0000....
## 14  0.08971189 0.3636301        1  49 0.020408163 ENSG0000....
## 15  1.00000000 1.0000000        0 135 0.000000000 ENSG0000....
## 16  0.24181594 0.3636301        1 144 0.006944444 ENSG0000....
## 17  0.26202152 0.3636301        1 158 0.006329114 ENSG0000....
# All individual iter results
head(iter$`k/K_iterations`)
##                             pathway k/K_1       k/K_2       k/K_3       k/K_4
## 10             HALLMARK_COAGULATION 0.000 0.014388489 0.014388489 0.014388489
## 11              HALLMARK_COMPLEMENT 0.000 0.004975124 0.004975124 0.004975124
## 13             HALLMARK_E2F_TARGETS 0.005 0.005000000 0.005000000 0.005000000
## 15 HALLMARK_ESTROGEN_RESPONSE_EARLY 0.005 0.005000000 0.005000000 0.005000000
## 16  HALLMARK_ESTROGEN_RESPONSE_LATE 0.005 0.005000000 0.005000000 0.005000000
## 18          HALLMARK_G2M_CHECKPOINT 0.010 0.010000000 0.010000000 0.010000000
##          k/K_5       k/K_6       k/K_7 k/K_8       k/K_9 k/K_10
## 10 0.014388489 0.014388489 0.014388489 0.000 0.014388489  0.000
## 11 0.004975124 0.004975124 0.004975124 0.000 0.004975124  0.000
## 13 0.005000000 0.005000000 0.005000000 0.005 0.005000000  0.005
## 15 0.005000000 0.005000000 0.005000000 0.005 0.005000000  0.005
## 16 0.005000000 0.005000000 0.005000000 0.005 0.005000000  0.005
## 18 0.010000000 0.010000000 0.010000000 0.010 0.010000000  0.010
head(iter$p_iterations)
##                             pathway     pval_1     pval_2     pval_3     pval_4
## 10             HALLMARK_COAGULATION 1.00000000 0.02908510 0.02908510 0.02908510
## 11              HALLMARK_COMPLEMENT 1.00000000 0.32085006 0.32085006 0.32085006
## 13             HALLMARK_E2F_TARGETS 0.30054682 0.31953545 0.31953545 0.31953545
## 15 HALLMARK_ESTROGEN_RESPONSE_EARLY 0.30054682 0.31953545 0.31953545 0.31953545
## 16  HALLMARK_ESTROGEN_RESPONSE_LATE 0.30054682 0.31953545 0.31953545 0.31953545
## 18          HALLMARK_G2M_CHECKPOINT 0.04915364 0.05611768 0.05611768 0.05611768
##        pval_5     pval_6     pval_7     pval_8     pval_9    pval_10
## 10 0.02908510 0.02908510 0.02908510 1.00000000 0.02908510 1.00000000
## 11 0.32085006 0.32085006 0.32085006 1.00000000 0.32085006 1.00000000
## 13 0.31953545 0.31953545 0.31953545 0.30054682 0.31953545 0.30054682
## 15 0.31953545 0.31953545 0.31953545 0.30054682 0.31953545 0.30054682
## 16 0.31953545 0.31953545 0.31953545 0.30054682 0.31953545 0.30054682
## 18 0.05611768 0.05611768 0.05611768 0.04915364 0.05611768 0.04915364

Change the background to a custom one (custom_bg) or include all genes by turning off protein_coding.

iterEnrich(
  anno_df = example.gene.anno,
  ID = "ENSEMBL", species = "human", 
  collection = "H",
  niter = 10, ncores = 6,
  protein_coding = FALSE)

Restrict the gene sets included to those with at least 3 query genes and total size between 100 and 1000 genes.

iterEnrich(
  anno_df = example.gene.anno,
  ID = "ENSEMBL", species = "human", 
  collection = "H",
  niter = 10, ncores = 6,
  minOverlap = 3,
  minGeneSetSize = 100,
  maxGeneSetSize = 1000)

Speed up computation by not saving the query genes found in each gene set.

iterEnrich(
  anno_df = example.gene.anno,
  ID = "ENSEMBL", species = "human", 
  collection = "H",
  niter = 10, ncores = 6,
  print_genes = FALSE)

Other parameters

  • Change the column names to match your annotation data with anno_featCol and anno_annotationCol
  • Change the FDR method with p_adjust. See p.adjust.methods for options

GSEA

Gene set enrichment analysis (GSEA) as described by Broad. This method utilizes the fgsea framework.

BIGsea(
  gene_list = example.gene.list,
  ID = "ENSEMBL", species = "human", 
  collection = "C2", subcollection = "CP")
## Running HRV1
## Running HRV2
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (1% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
##   group gs_collection gs_subcollection method
## 1  HRV1            C2               CP  multi
## 2  HRV2            C2               CP  multi
## 3  HRV2            C2               CP  multi
##                                            pathway      pval       FDR
## 1            REACTOME_TRANSPORT_OF_SMALL_MOLECULES 0.8736264 0.8736264
## 2 REACTOME_POST_TRANSLATIONAL_PROTEIN_MODIFICATION 0.5366430 0.5366430
## 3         REACTOME_RNA_POLYMERASE_II_TRANSCRIPTION 0.1465721 0.2931442
##      log2err        ES       NES size  leadingEdge
## 1 0.04782996 0.2925544 0.6594805   10 ENSG0000....
## 2 0.08431443 0.4127165 0.9346453   10 ENSG0000....
## 3 0.17821987 0.6103486 1.3822065   10 ENSG0000....

Use a data frame input.

BIGsea(
  gene_list = example.gene.df,
  ID = "ENSEMBL", species = "human", 
  collection = "C2", subcollection = "CP")

Restrict the gene sets included to those with total size between 100 and 1000 genes.

BIGsea(
  gene_list = example.gene.list,
  ID = "ENSEMBL", species = "human", 
  collection = "C2", subcollection = "CP",
  minGeneSetSize = 100, maxGeneSetSize = 1000)

Change the GSEA method. By default, rand = "multi", which runs an adaptive multilevel splitting Monte Carlo approach.

For very small sample groups, it is recommended to run “simple” randomization, which randomizes genes in gene sets only.

BIGsea(
  gene_list = example.gene.list,
  ID = "ENSEMBL", species = "human", 
  collection = "C2", subcollection = "CP",
  rand = "simple")

Alternatively, choose the “label” method to randomize samples in groups and calculate new fold changes for each randomization. This method runs kimma under the hood to determine fold changes. Thus, you have all the model flexibility of kimma, not just a simple A vs B fold change.

The “label” method requires dat and rand_var in order to calculate fold changes from the original expression data. You also must provide kimma::kmFit model parameters model, run_, use_weights, etc.

Finally, we strongly recommend only running on select gene sets with pw and reducing computational time with multiple processors. Here, we run only 3 permutations for speed as well.

#Force original fold change list name to match example.voom$targets variable
example.gene.list.virus <- list(
  "virus" = example.gene.list$HRV1
)

gsea_label <- BIGsea(
  gene_list = example.gene.list.virus,
  ID = "ENSEMBL", species = "human", 
  collection = "C2", subcollection = "CP",
  rand = "label",
  #Data for fold change calculations
  dat = kimma::example.voom, 
  rand_var = "virus",
  #Parameters for kimma models
  model="~virus",
  run_lm=TRUE, use_weights=TRUE,
  #Recommend only running on select gene sets of interest
  pw = "REACTOME_TRANSPORT_OF_SMALL_MOLECULES",
  nperm = 3, processors = 3)
## Running virus
## Loading required namespace: limma
## Label randomization runs multiple kimma models. This may take a long time.
## [1] "kimma permutation 1"
## [1] "kimma permutation 2"
## [1] "kimma permutation 3"
##   |                                                                              |                                                                      |   0%  |                                                                              |=======================                                               |  33%  |                                                                              |===============================================                       |  67%  |                                                                              |======================================================================| 100%
#GSEA results
gsea_label$gsea
##   group gs_collection gs_subcollection method
## 1 virus            C2               CP  label
##                                 pathway pval FDR        ES NES nMoreExtreme
## 1 REACTOME_TRANSPORT_OF_SMALL_MOLECULES    1   1 0.2925544 NaN            0
##   size  leadingEdge
## 1   10 ENSG0000....
#Fold change estimates for each permutation
head(gsea_label$estimates)
##                          V1         V2          V3
## ENSG00000000460 -0.48411710 -0.3572232 -0.22709577
## ENSG00000001460  0.69172192 -0.4946608 -0.06399949
## ENSG00000002587 -0.02961869  0.1458342  0.37694009
## ENSG00000005189  1.72116622  0.6523883  0.07678363
## ENSG00000005801  0.02720169  0.2159433 -0.19241759
## ENSG00000005884  0.89349003  0.1290022 -1.49447436

Once you’ve run the “label” method once, you can reuse the fold change permutations for additional gene sets with rand_est.

gsea_label2 <- BIGsea(
  gene_list = example.gene.list.virus,
  ID = "ENSEMBL", species = "human", 
  collection = "C2", subcollection = "CP",
  rand = "label",
  #Data for fold change calculations
  dat = kimma::example.voom,
  rand_var = "virus",
  #Pre-made fold change results
  rand_est = gsea_label$estimates,
  #Recommend only running on select gene sets of interest
  pw = "REACTOME_POST_TRANSLATIONAL_PROTEIN_MODIFICATION",
  nperm = 3, processors = 3)
## Running virus

R session

sessionInfo()
## R version 4.5.0 (2025-04-11)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sequoia 15.4.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Los_Angeles
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] SEARchways_1.5.0 kimma_1.6.3      lubridate_1.9.4  forcats_1.0.0   
##  [5] stringr_1.5.1    dplyr_1.1.4      purrr_1.0.4      readr_2.1.5     
##  [9] tidyr_1.3.1      tibble_3.2.1     ggplot2_3.5.2    tidyverse_2.0.0 
## 
## loaded via a namespace (and not attached):
##   [1] DBI_1.2.3               gson_0.1.0              rlang_1.1.6            
##   [4] magrittr_2.0.3          DOSE_4.2.0              msigdbr_10.0.2         
##   [7] compiler_4.5.0          RSQLite_2.3.9           png_0.1-8              
##  [10] vctrs_0.6.5             reshape2_1.4.4          pkgconfig_2.0.3        
##  [13] crayon_1.5.3            fastmap_1.2.0           XVector_0.48.0         
##  [16] utf8_1.2.4              rmarkdown_2.29          tzdb_0.5.0             
##  [19] enrichplot_1.28.1       UCSC.utils_1.4.0        bit_4.6.0              
##  [22] xfun_0.52               cachem_1.1.0            aplot_0.2.5            
##  [25] GenomeInfoDb_1.44.0     jsonlite_2.0.0          blob_1.2.4             
##  [28] BiocParallel_1.42.0     parallel_4.5.0          R6_2.6.1               
##  [31] bslib_0.9.0             stringi_1.8.7           RColorBrewer_1.1-3     
##  [34] limma_3.64.0            msigdbdf_24.1.1         jquerylib_0.1.4        
##  [37] GOSemSim_2.34.0         Rcpp_1.0.14             assertthat_0.2.1       
##  [40] iterators_1.0.14        knitr_1.50              ggtangle_0.0.6         
##  [43] R.utils_2.13.0          IRanges_2.42.0          igraph_2.1.4           
##  [46] Matrix_1.7-3            splines_4.5.0           timechange_0.3.0       
##  [49] tidyselect_1.2.1        qvalue_2.40.0           rstudioapi_0.17.1      
##  [52] yaml_2.3.10             doParallel_1.0.17       codetools_0.2-20       
##  [55] lattice_0.22-7          plyr_1.8.9              treeio_1.32.0          
##  [58] Biobase_2.68.0          withr_3.0.2             KEGGREST_1.48.0        
##  [61] evaluate_1.0.3          gridGraphics_0.5-1      Biostrings_2.76.0      
##  [64] ggtree_3.16.0           pillar_1.10.2           foreach_1.5.2          
##  [67] stats4_4.5.0            clusterProfiler_4.16.0  ggfun_0.1.8            
##  [70] generics_0.1.3          S4Vectors_0.46.0        hms_1.1.3              
##  [73] tidytree_0.4.6          scales_1.4.0            glue_1.8.0             
##  [76] lazyeval_0.2.2          tools_4.5.0             data.table_1.17.0      
##  [79] fgsea_1.34.0            babelgene_22.9          fs_1.6.6               
##  [82] fastmatch_1.1-6         cowplot_1.1.3           grid_4.5.0             
##  [85] ape_5.8-1               AnnotationDbi_1.70.0    colorspace_2.1-1       
##  [88] nlme_3.1-168            patchwork_1.3.0         GenomeInfoDbData_1.2.14
##  [91] cli_3.6.5               gtable_0.3.6            R.methodsS3_1.8.2      
##  [94] yulab.utils_0.2.0       sass_0.4.10             digest_0.6.37          
##  [97] BiocGenerics_0.54.0     ggrepel_0.9.6           ggplotify_0.1.2        
## [100] farver_2.1.2            memoise_2.0.1           htmltools_0.5.8.1      
## [103] R.oo_1.27.0             lifecycle_1.0.4         httr_1.4.7             
## [106] statmod_1.5.0           GO.db_3.21.0            bit64_4.6.0-1