Analysing single-cell RNA-sequencing Data

Johannes Griss

2023-04-25

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")
#> Loading required package: ReactomeGSA

# install the ReactomeGSA.data package for the example data
if (!require(ReactomeGSA.data))
  BiocManager::install("ReactomeGSA.data")
#> Loading required package: ReactomeGSA.data
#> Loading required package: limma
#> Loading required package: edgeR
#> Loading required package: Seurat
#> Attaching SeuratObject

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)
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)

Pathway 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: 84
#>   Results:
#>   - Seurat:
#>     1752 pathways
#>     11028 fold changes for genes
#>   No Reactome visualizations available
#> ReactomeAnalysisResult

pathways 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.1068604 0.09743804  0.1420590
#> R-HSA-109606  Intrinsic Pathway for Apoptosis 0.1132631 0.10959578  0.1128281
#> R-HSA-109703              PKB-mediated events 0.1273525 0.05283679  0.1066424
#>               Cluster_12 Cluster_13 Cluster_2  Cluster_3  Cluster_4  Cluster_5
#> R-HSA-1059683 0.11074994 0.10043978 0.1146915 0.11282191 0.11003725 0.10353440
#> R-HSA-109606  0.11912879 0.12836408 0.1081178 0.11094937 0.11310679 0.10564916
#> R-HSA-109703  0.09574043 0.07376632 0.0835158 0.08429958 0.05591271 0.04644355
#>                Cluster_6  Cluster_7  Cluster_8  Cluster_9
#> R-HSA-1059683 0.09556152 0.12114396 0.13516241 0.10128279
#> R-HSA-109606  0.10996000 0.11850204 0.12124599 0.11416078
#> R-HSA-109703  0.12405405 0.07729952 0.07839453 0.01445622

A 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        min
#> R-HSA-350864           Regulation of thyroid hormone activity -0.4876651
#> R-HSA-8964540                              Alanine metabolism -0.5061335
#> R-HSA-190374  FGFR1c and Klotho ligand binding and activation -0.3432262
#> R-HSA-140180                                    COX reactions -0.3451557
#> R-HSA-9024909           BDNF activates NTRK2 (TRKB) signaling -0.3744760
#> R-HSA-5263617       Metabolism of ingested MeSeO2H into MeSeH -0.1936180
#>                     max      diff
#> R-HSA-350864  0.3757525 0.8634177
#> R-HSA-8964540 0.2563056 0.7624391
#> R-HSA-190374  0.4160889 0.7593151
#> R-HSA-140180  0.3726285 0.7177843
#> R-HSA-9024909 0.3236443 0.6981203
#> R-HSA-5263617 0.4938665 0.6874845

Plotting the results

The ReactomeGSA package contains two functions to visualize these pathway results. The first simply plots the expression for a selected pathway:

plot_gsva_pathway(gsva_result, pathway_id = rownames(max_difference)[1])

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

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:

plot_gsva_pca(gsva_result)

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.3.0 RC (2023-04-13 r84269)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.2 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.17-bioc/R/lib/libRblas.so 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_GB              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: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ReactomeGSA.data_1.13.0 SeuratObject_4.1.3      Seurat_4.3.0           
#> [4] edgeR_3.42.0            limma_3.56.0            ReactomeGSA_1.14.0     
#> 
#> loaded via a namespace (and not attached):
#>   [1] RColorBrewer_1.1-3     jsonlite_1.8.4         magrittr_2.0.3        
#>   [4] spatstat.utils_3.0-2   farver_2.1.1           rmarkdown_2.21        
#>   [7] vctrs_0.6.2            ROCR_1.0-11            spatstat.explore_3.1-0
#>  [10] htmltools_0.5.5        progress_1.2.2         curl_5.0.0            
#>  [13] sass_0.4.5             sctransform_0.3.5      parallelly_1.35.0     
#>  [16] KernSmooth_2.23-20     bslib_0.4.2            htmlwidgets_1.6.2     
#>  [19] ica_1.0-3              plyr_1.8.8             plotly_4.10.1         
#>  [22] zoo_1.8-12             cachem_1.0.7           igraph_1.4.2          
#>  [25] mime_0.12              lifecycle_1.0.3        pkgconfig_2.0.3       
#>  [28] Matrix_1.5-4           R6_2.5.1               fastmap_1.1.1         
#>  [31] fitdistrplus_1.1-11    future_1.32.0          shiny_1.7.4           
#>  [34] digest_0.6.31          colorspace_2.1-0       patchwork_1.1.2       
#>  [37] tensor_1.5             irlba_2.3.5.1          labeling_0.4.2        
#>  [40] progressr_0.13.0       fansi_1.0.4            spatstat.sparse_3.0-1 
#>  [43] httr_1.4.5             polyclip_1.10-4        abind_1.4-5           
#>  [46] compiler_4.3.0         withr_2.5.0            highr_0.10            
#>  [49] gplots_3.1.3           MASS_7.3-59            gtools_3.9.4          
#>  [52] caTools_1.18.2         tools_4.3.0            lmtest_0.9-40         
#>  [55] httpuv_1.6.9           future.apply_1.10.0    goftest_1.2-3         
#>  [58] glue_1.6.2             nlme_3.1-162           promises_1.2.0.1      
#>  [61] grid_4.3.0             Rtsne_0.16             cluster_2.1.4         
#>  [64] reshape2_1.4.4         generics_0.1.3         gtable_0.3.3          
#>  [67] spatstat.data_3.0-1    tidyr_1.3.0            data.table_1.14.8     
#>  [70] hms_1.1.3              sp_1.6-0               utf8_1.2.3            
#>  [73] spatstat.geom_3.1-0    RcppAnnoy_0.0.20       ggrepel_0.9.3         
#>  [76] RANN_2.6.1             pillar_1.9.0           stringr_1.5.0         
#>  [79] later_1.3.0            splines_4.3.0          dplyr_1.1.2           
#>  [82] lattice_0.21-8         survival_3.5-5         deldir_1.0-6          
#>  [85] tidyselect_1.2.0       locfit_1.5-9.7         miniUI_0.1.1.1        
#>  [88] pbapply_1.7-0          knitr_1.42             gridExtra_2.3         
#>  [91] scattermore_0.8        xfun_0.39              matrixStats_0.63.0    
#>  [94] stringi_1.7.12         lazyeval_0.2.2         yaml_2.3.7            
#>  [97] evaluate_0.20          codetools_0.2-19       tibble_3.2.1          
#> [100] BiocManager_1.30.20    cli_3.6.1              uwot_0.1.14           
#> [103] xtable_1.8-4           reticulate_1.28        munsell_0.5.0         
#> [106] jquerylib_0.1.4        Rcpp_1.0.10            globals_0.16.2        
#> [109] spatstat.random_3.1-4  png_0.1-8              parallel_4.3.0        
#> [112] ellipsis_0.3.2         ggplot2_3.4.2          prettyunits_1.1.1     
#> [115] bitops_1.0-7           listenv_0.9.0          viridisLite_0.4.1     
#> [118] scales_1.2.1           ggridges_0.5.4         leiden_0.4.3          
#> [121] purrr_1.0.1            crayon_1.5.2           rlang_1.1.0           
#> [124] cowplot_1.1.1