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.

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))
  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: Seurat
#> Attaching SeuratObject
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...
#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.
#> Compressing request data...
#> Reactome Analysis submitted succesfully
#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.
#> Converting dataset Seurat...
#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.
#> Mapping identifiers...
#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.
#> Performing gene set analysis using ssGSEA
#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.
#> Analysing dataset 'Seurat' using ssGSEA
#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.
#> Retrieving result...
#> Warning in value[[3L]](cond): Failed to process certificate. This is generally
#> caused by a problem with Windows and R. Reverting to non-encrypted connection to
#> the ReactomeGSA system.

The resulting object is a standard ReactomeAnalysisResult object.

gsva_result
#> ReactomeAnalysisResult object
#>   Reactome Release: 75
#>   Results:
#>   - Seurat:
#>     1726 pathways
#>     10877 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.1032472 0.10008086  0.1497649
#> R-HSA-109606  Intrinsic Pathway for Apoptosis 0.1095918 0.10716671  0.1123212
#> R-HSA-109703              PKB-mediated events 0.1271271 0.05247245  0.1060551
#>               Cluster_12 Cluster_13  Cluster_2  Cluster_3  Cluster_4  Cluster_5
#> R-HSA-1059683 0.10709318  0.1079596 0.11536591 0.11046754 0.11324535 0.10860458
#> R-HSA-109606  0.11484536  0.1264604 0.10426341 0.10737779 0.11095973 0.10312804
#> R-HSA-109703  0.09532549  0.0729352 0.08280237 0.08390533 0.05531344 0.04627873
#>                Cluster_6  Cluster_7  Cluster_8  Cluster_9
#> R-HSA-1059683 0.09265093 0.11630862 0.13375275 0.10534315
#> R-HSA-109606  0.10819825 0.11376882 0.11859846 0.11439947
#> R-HSA-109703  0.12356494 0.07705906 0.07798356 0.01401364

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.4872121
#> R-HSA-8964540                              Alanine metabolism -0.5062198
#> R-HSA-190374  FGFR1c and Klotho ligand binding and activation -0.3440953
#> R-HSA-140180                                    COX reactions -0.3451202
#> R-HSA-9024909           BDNF activates NTRK2 (TRKB) signaling -0.3747605
#> R-HSA-9025046           NTF3 activates NTRK2 (TRKB) signaling -0.3957010
#>                     max      diff
#> R-HSA-350864  0.3746567 0.8618688
#> R-HSA-8964540 0.2550200 0.7612398
#> R-HSA-190374  0.4152810 0.7593763
#> R-HSA-140180  0.3719546 0.7170748
#> R-HSA-9024909 0.3228439 0.6976044
#> R-HSA-9025046 0.2986478 0.6943488

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.0.5 (2021-03-31)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows Server 2012 R2 x64 (build 9600)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=C                          
#> [2] LC_CTYPE=English_United States.1252   
#> [3] LC_MONETARY=English_United States.1252
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.1252    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ReactomeGSA.data_1.4.0 SeuratObject_4.0.0     Seurat_4.0.1          
#> [4] edgeR_3.32.1           limma_3.46.0           ReactomeGSA_1.4.2     
#> 
#> loaded via a namespace (and not attached):
#>   [1] Rtsne_0.15            colorspace_2.0-0      deldir_0.2-10        
#>   [4] ellipsis_0.3.1        ggridges_0.5.3        spatstat.data_2.1-0  
#>   [7] farver_2.1.0          leiden_0.3.7          listenv_0.8.0        
#>  [10] ggrepel_0.9.1         fansi_0.4.2           codetools_0.2-18     
#>  [13] splines_4.0.5         knitr_1.32            polyclip_1.10-0      
#>  [16] jsonlite_1.7.2        ica_1.0-2             cluster_2.1.1        
#>  [19] png_0.1-7             uwot_0.1.10           shiny_1.6.0          
#>  [22] sctransform_0.3.2     spatstat.sparse_2.0-0 BiocManager_1.30.12  
#>  [25] compiler_4.0.5        httr_1.4.2            assertthat_0.2.1     
#>  [28] Matrix_1.3-2          fastmap_1.1.0         lazyeval_0.2.2       
#>  [31] later_1.1.0.1         prettyunits_1.1.1     htmltools_0.5.1.1    
#>  [34] tools_4.0.5           igraph_1.2.6          gtable_0.3.0         
#>  [37] glue_1.4.2            RANN_2.6.1            reshape2_1.4.4       
#>  [40] dplyr_1.0.5           Rcpp_1.0.6            scattermore_0.7      
#>  [43] jquerylib_0.1.3       vctrs_0.3.7           nlme_3.1-152         
#>  [46] lmtest_0.9-38         xfun_0.22             stringr_1.4.0        
#>  [49] globals_0.14.0        mime_0.10             miniUI_0.1.1.1       
#>  [52] lifecycle_1.0.0       irlba_2.3.3           gtools_3.8.2         
#>  [55] goftest_1.2-2         future_1.21.0         MASS_7.3-53.1        
#>  [58] zoo_1.8-9             scales_1.1.1          spatstat.core_2.0-0  
#>  [61] hms_1.0.0             promises_1.2.0.1      spatstat.utils_2.1-0 
#>  [64] parallel_4.0.5        RColorBrewer_1.1-2    curl_4.3             
#>  [67] yaml_2.2.1            reticulate_1.18       pbapply_1.4-3        
#>  [70] gridExtra_2.3         ggplot2_3.3.3         sass_0.3.1           
#>  [73] rpart_4.1-15          stringi_1.5.3         highr_0.9            
#>  [76] caTools_1.18.2        rlang_0.4.10          pkgconfig_2.0.3      
#>  [79] matrixStats_0.58.0    bitops_1.0-6          evaluate_0.14        
#>  [82] lattice_0.20-41       ROCR_1.0-11           purrr_0.3.4          
#>  [85] tensor_1.5            labeling_0.4.2        patchwork_1.1.1      
#>  [88] htmlwidgets_1.5.3     cowplot_1.1.1         tidyselect_1.1.0     
#>  [91] parallelly_1.24.0     RcppAnnoy_0.0.18      plyr_1.8.6           
#>  [94] magrittr_2.0.1        R6_2.5.0              gplots_3.1.1         
#>  [97] generics_0.1.0        DBI_1.1.1             mgcv_1.8-34          
#> [100] pillar_1.6.0          fitdistrplus_1.1-3    survival_3.2-10      
#> [103] abind_1.4-5           tibble_3.1.0          future.apply_1.7.0   
#> [106] crayon_1.4.1          KernSmooth_2.23-18    utf8_1.2.1           
#> [109] spatstat.geom_2.1-0   plotly_4.9.3          rmarkdown_2.7        
#> [112] progress_1.2.2        locfit_1.5-9.4        grid_4.0.5           
#> [115] data.table_1.14.0     digest_0.6.27         xtable_1.8-4         
#> [118] tidyr_1.1.3           httpuv_1.5.5          munsell_0.5.0        
#> [121] viridisLite_0.4.0     bslib_0.2.4