Analysing single-cell RNA-sequencing Data
Introduction
The ReactomeGSA package is a client to the web-based Reactome Analysis System. Essentially, it performs a gene set analysis using the latest version of the Reactome pathway database as a backend.
This vignette shows how the ReactomeGSA package can be used to perform a pathway analysis of cell clusters in single-cell RNA-sequencing data.
Citation
To cite this package, use
Griss J. ReactomeGSA, https://github.com/reactome/ReactomeGSA (2019)Installation
The ReactomeGSA package can be directly installed from
Bioconductor:
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
if (!require(ReactomeGSA))
BiocManager::install("ReactomeGSA")
# install the ReactomeGSA.data package for the example data
if (!require(ReactomeGSA.data))
BiocManager::install("ReactomeGSA.data")For more information, see https://bioconductor.org/install/.
Example data
As an example we load single-cell RNA-sequencing data of B cells extracted from the dataset published by Jerby-Arnon et al. (Cell, 2018).
Note: This is not a complete Seurat object. To decrease the size, the object only contains gene expression values and cluster annotations.
library(ReactomeGSA.data)
#> Loading required package: limma
#> Loading required package: edgeR
#> Loading required package: ReactomeGSA
#> Loading required package: Seurat
#> Loading required package: SeuratObject
#> Loading required package: sp
#> 'SeuratObject' was built under R 4.4.0 but the current version is
#> 4.4.1; it is recomended that you reinstall 'SeuratObject' as the ABI
#> for R may have changed
#> 'SeuratObject' was built with package 'Matrix' 1.7.0 but the current
#> version is 1.7.1; it is recomended that you reinstall 'SeuratObject' as
#> the ABI for 'Matrix' may have changed
#>
#> Attaching package: 'SeuratObject'
#> The following objects are masked from 'package:base':
#>
#> intersect, t
data(jerby_b_cells)
jerby_b_cells
#> An object of class Seurat
#> 23686 features across 920 samples within 1 assay
#> Active assay: RNA (23686 features, 0 variable features)
#> 2 layers present: counts, dataPathway analysis of cell clusters
The pathway analysis is at the very end of a scRNA-seq workflow. This means, that any Q/C was already performed, the data was normalized and cells were already clustered.
The ReactomeGSA package can now be used to get pathway-level expression values for every cell cluster. This is achieved by calculating the mean gene expression for every cluster and then submitting this data to a gene set variation analysis.
All of this is wrapped in the single analyse_sc_clusters
function.
library(ReactomeGSA)
gsva_result <- analyse_sc_clusters(jerby_b_cells, verbose = TRUE)
#> Calculating average cluster expression...
#> Converting expression data to string... (This may take a moment)
#> Conversion complete
#> Submitting request to Reactome API...
#> Compressing request data...
#> Reactome Analysis submitted succesfully
#> Converting dataset Seurat...
#> Mapping identifiers...
#> Performing gene set analysis using ssGSEA
#> Analysing dataset 'Seurat' using ssGSEA
#> Retrieving result...The resulting object is a standard
ReactomeAnalysisResult object.
gsva_result
#> ReactomeAnalysisResult object
#> Reactome Release: 90
#> Results:
#> - Seurat:
#> 1741 pathways
#> 11759 fold changes for genes
#> No Reactome visualizations available
#> ReactomeAnalysisResultpathways returns the pathway-level expression values per
cell cluster:
pathway_expression <- pathways(gsva_result)
# simplify the column names by removing the default dataset identifier
colnames(pathway_expression) <- gsub("\\.Seurat", "", colnames(pathway_expression))
pathway_expression[1:3,]
#> Name Cluster_1 Cluster_10 Cluster_11
#> R-HSA-1059683 Interleukin-6 signaling -0.03556065 -0.04521415 0.14931444
#> R-HSA-109703 PKB-mediated events 0.32848879 -0.19771237 0.05087568
#> R-HSA-109704 PI3K Cascade -0.12080150 -0.13450596 0.16088061
#> Cluster_12 Cluster_13 Cluster_2 Cluster_3 Cluster_4
#> R-HSA-1059683 0.06075330 0.0004459001 0.14896349 -0.09445591 -0.11884117
#> R-HSA-109703 0.09556102 0.0274677841 -0.05142555 0.15462891 -0.16856699
#> R-HSA-109704 -0.04114424 0.0680222582 0.05924518 0.05295837 -0.02208554
#> Cluster_5 Cluster_6 Cluster_7 Cluster_8 Cluster_9
#> R-HSA-1059683 -0.13882420 -0.07577576 -0.06560393 0.174814091 -0.04388071
#> R-HSA-109703 0.10554980 -0.06251424 -0.30506281 -0.007670405 -0.09929333
#> R-HSA-109704 0.02008053 -0.12603775 -0.01371145 0.101810119 -0.04142070A simple approach to find the most relevant pathways is to assess the maximum difference in expression for every pathway:
# find the maximum differently expressed pathway
max_difference <- do.call(rbind, apply(pathway_expression, 1, function(row) {
values <- as.numeric(row[2:length(row)])
return(data.frame(name = row[1], min = min(values), max = max(values)))
}))
max_difference$diff <- max_difference$max - max_difference$min
# sort based on the difference
max_difference <- max_difference[order(max_difference$diff, decreasing = T), ]
head(max_difference)
#> name
#> R-HSA-140180 COX reactions
#> R-HSA-1296067 Potassium transport channels
#> R-HSA-392023 Adrenaline signalling through Alpha-2 adrenergic receptor
#> R-HSA-8964540 Alanine metabolism
#> R-HSA-3248023 Regulation by TREX1
#> R-HSA-350864 Regulation of thyroid hormone activity
#> min max diff
#> R-HSA-140180 -0.9647899 0.9841810 1.948971
#> R-HSA-1296067 -1.0000000 0.8858649 1.885865
#> R-HSA-392023 -0.9008335 0.9738051 1.874639
#> R-HSA-8964540 -0.8731077 0.9938765 1.866984
#> R-HSA-3248023 -0.9185236 0.9421670 1.860691
#> R-HSA-350864 -0.9134206 0.9394455 1.852866Plotting the results
The ReactomeGSA package contains two functions to visualize these pathway results. The first simply plots the expression for a selected pathway:
For a better overview, the expression of multiple pathways can be
shown as a heatmap using gplots heatmap.2
function:
# Additional parameters are directly passed to gplots heatmap.2 function
plot_gsva_heatmap(gsva_result, max_pathways = 15, margins = c(6,20))
The plot_gsva_heatmap function can also be used to only
display specific pahtways:
# limit to selected B cell related pathways
relevant_pathways <- c("R-HSA-983170", "R-HSA-388841", "R-HSA-2132295", "R-HSA-983705", "R-HSA-5690714")
plot_gsva_heatmap(gsva_result,
pathway_ids = relevant_pathways, # limit to these pathways
margins = c(6,30), # adapt the figure margins in heatmap.2
dendrogram = "col", # only plot column dendrogram
scale = "row", # scale for each pathway
key = FALSE, # don't display the color key
lwid=c(0.1,4)) # remove the white space on the left
#> Warning in plot_gsva_heatmap(gsva_result, pathway_ids = relevant_pathways, :
#> Warning: No results for the following pathways: R-HSA-983705
This analysis shows us that cluster 8 has a marked up-regulation of B Cell receptor signalling, which is linked to a co-stimulation of the CD28 family. Additionally, there is a gradient among the cluster with respect to genes releated to antigen presentation.
Therefore, we are able to further classify the observed B cell subtypes based on their pathway activity.
Pathway-level PCA
The pathway-level expression analysis can also be used to run a
Principal Component Analysis on the samples. This is simplified through
the function plot_gsva_pca:
In this analysis, cluster 11 is a clear outlier from the other B cell subtypes and therefore might be prioritised for further evaluation.
Session Info
sessionInfo()
#> R version 4.4.1 (2024-06-14)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] ReactomeGSA.data_1.19.0 Seurat_5.1.0 SeuratObject_5.0.2
#> [4] sp_2.1-4 ReactomeGSA_1.21.0 edgeR_4.4.0
#> [7] limma_3.63.0 rmarkdown_2.28
#>
#> loaded via a namespace (and not attached):
#> [1] RColorBrewer_1.1-3 sys_3.4.3 jsonlite_1.8.9
#> [4] magrittr_2.0.3 spatstat.utils_3.1-0 farver_2.1.2
#> [7] vctrs_0.6.5 ROCR_1.0-11 spatstat.explore_3.3-3
#> [10] progress_1.2.3 htmltools_0.5.8.1 curl_5.2.3
#> [13] sass_0.4.9 sctransform_0.4.1 parallelly_1.38.0
#> [16] KernSmooth_2.23-24 bslib_0.8.0 htmlwidgets_1.6.4
#> [19] ica_1.0-3 plyr_1.8.9 plotly_4.10.4
#> [22] zoo_1.8-12 cachem_1.1.0 buildtools_1.0.0
#> [25] igraph_2.1.1 mime_0.12 lifecycle_1.0.4
#> [28] pkgconfig_2.0.3 Matrix_1.7-1 R6_2.5.1
#> [31] fastmap_1.2.0 fitdistrplus_1.2-1 future_1.34.0
#> [34] shiny_1.9.1 digest_0.6.37 colorspace_2.1-1
#> [37] patchwork_1.3.0 tensor_1.5 RSpectra_0.16-2
#> [40] irlba_2.3.5.1 labeling_0.4.3 progressr_0.15.0
#> [43] fansi_1.0.6 spatstat.sparse_3.1-0 httr_1.4.7
#> [46] polyclip_1.10-7 abind_1.4-8 compiler_4.4.1
#> [49] withr_3.0.2 fastDummies_1.7.4 highr_0.11
#> [52] gplots_3.2.0 MASS_7.3-61 gtools_3.9.5
#> [55] caTools_1.18.3 tools_4.4.1 lmtest_0.9-40
#> [58] httpuv_1.6.15 future.apply_1.11.3 goftest_1.2-3
#> [61] glue_1.8.0 nlme_3.1-166 promises_1.3.0
#> [64] grid_4.4.1 Rtsne_0.17 cluster_2.1.6
#> [67] reshape2_1.4.4 generics_0.1.3 gtable_0.3.6
#> [70] spatstat.data_3.1-2 tidyr_1.3.1 hms_1.1.3
#> [73] data.table_1.16.2 utf8_1.2.4 BiocGenerics_0.53.0
#> [76] spatstat.geom_3.3-3 RcppAnnoy_0.0.22 ggrepel_0.9.6
#> [79] RANN_2.6.2 pillar_1.9.0 stringr_1.5.1
#> [82] spam_2.11-0 RcppHNSW_0.6.0 later_1.3.2
#> [85] splines_4.4.1 dplyr_1.1.4 lattice_0.22-6
#> [88] survival_3.7-0 deldir_2.0-4 tidyselect_1.2.1
#> [91] locfit_1.5-9.10 maketools_1.3.1 miniUI_0.1.1.1
#> [94] pbapply_1.7-2 knitr_1.48 gridExtra_2.3
#> [97] scattermore_1.2 xfun_0.48 Biobase_2.67.0
#> [100] statmod_1.5.0 matrixStats_1.4.1 stringi_1.8.4
#> [103] lazyeval_0.2.2 yaml_2.3.10 evaluate_1.0.1
#> [106] codetools_0.2-20 tibble_3.2.1 cli_3.6.3
#> [109] uwot_0.2.2 xtable_1.8-4 reticulate_1.39.0
#> [112] munsell_0.5.1 jquerylib_0.1.4 Rcpp_1.0.13
#> [115] globals_0.16.3 spatstat.random_3.3-2 png_0.1-8
#> [118] spatstat.univar_3.0-1 parallel_4.4.1 ggplot2_3.5.1
#> [121] prettyunits_1.2.0 dotCall64_1.2 bitops_1.0-9
#> [124] listenv_0.9.1 viridisLite_0.4.2 scales_1.3.0
#> [127] ggridges_0.5.6 crayon_1.5.3 leiden_0.4.3.1
#> [130] purrr_1.0.2 rlang_1.1.4 cowplot_1.1.3