Using the ReactomeGSA package
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.
The main advantages of using the Reactome Analysis System are:
- Simultaneous analysis and visualization of different types of ’omics data
- Support for direct comparison across different species
- Directly linked with Reactome’s powerful pathway browser producing publication-ready figures of your gene set analysis
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")For more information, see https://bioconductor.org/install/.
Getting available methods
The Reactome Analysis System will be continuously updated. Before starting your analysis it is therefore a good approach to check which methods are available.
This can simply be done by using:
library(ReactomeGSA)
available_methods <- get_reactome_methods(print_methods = FALSE, return_result = TRUE)
# only show the names of the available methods
available_methods$name
#> [1] "PADOG" "Camera" "ssGSEA"To get more information about a specific method, set
print_details to TRUE and specify the
method:
# Use this command to print the description of the specific method to the console
# get_reactome_methods(print_methods = TRUE, print_details = TRUE, method = "PADOG", return_result = FALSE)
# show the parameter names for the method
padog_params <- available_methods$parameters[available_methods$name == "PADOG"][[1]]
paste0(padog_params$name, " (", padog_params$type, ", ", padog_params$default, ")")
#> [1] "use_interactors (bool, False)"
#> [2] "include_disease_pathways (bool, True)"
#> [3] "max_missing_values (float, 0.5)"
#> [4] "create_reactome_visualization (bool, True)"
#> [5] "create_reports (bool, False)"
#> [6] "email (string, )"
#> [7] "reactome_server (string, production)"
#> [8] "sample_groups (string, )"
#> [9] "discrete_norm_function (string, TMM)"
#> [10] "continuous_norm_function (string, none)"Creating an analysis request
To start a gene set analysis, you first have to create an analysis request. This is a simple S4 class that takes care of submitting multiple datasets simultaneously to the analysis system.
When creating the request object, you already have to specify the analysis method you want to use:
Setting parameters
To get a list of supported parameters for each method, use the
get_reactome_methods function (see above).
Parameters are simply set using the set_parameters
function:
# set the maximum number of allowed missing values to 50%
my_request <- set_parameters(request = my_request, max_missing_values = 0.5)
my_request
#> ReactomeAnalysisRequestObject
#> Method = Camera
#> Parameters:
#> - max_missing_values: 0.5
#> Datasets: none
#> ReactomeAnalysisRequestMultiple parameters can by set simulataneously by simply adding more name-value pairs to the function call.
Adding datasets
One analysis request can contain multiple datasets. This can be used to, for example, visualize the results of an RNA-seq and Proteomics experiment (of the same / similar samples) side by side:
This is a limma EList object with the sample data
already added
class(griss_melanoma_proteomics)
#> [1] "EList"
#> attr(,"package")
#> [1] "limma"
head(griss_melanoma_proteomics$samples)
#> patient.id condition cell.type
#> M-D MOCK PBMCB P3 MOCK PBMCB
#> M-D MCM PBMCB P3 MCM PBMCB
#> M-K MOCK PBMCB P4 MOCK PBMCB
#> M-K MCM PBMCB P4 MCM PBMCB
#> P-A MOCK PBMCB P1 MOCK PBMCB
#> P-A MCM PBMCB P1 MCM PBMCBThe dataset can now simply be added to the request using the
add_dataset function:
my_request <- add_dataset(request = my_request,
expression_values = griss_melanoma_proteomics,
name = "Proteomics",
type = "proteomics_int",
comparison_factor = "condition",
comparison_group_1 = "MOCK",
comparison_group_2 = "MCM",
additional_factors = c("cell.type", "patient.id"))
my_request
#> ReactomeAnalysisRequestObject
#> Method = Camera
#> Parameters:
#> - max_missing_values: 0.5
#> Datasets:
#> - Proteomics (proteomics_int)
#> No parameters set.
#> ReactomeAnalysisRequestSeveral datasets (of the same experiment) can be added to one
request. This RNA-seq data is stored as an edgeR
DGEList object:
data("griss_melanoma_rnaseq")
# only keep genes with >= 100 reads in total
total_reads <- rowSums(griss_melanoma_rnaseq$counts)
griss_melanoma_rnaseq <- griss_melanoma_rnaseq[total_reads >= 100, ]
# this is a edgeR DGEList object
class(griss_melanoma_rnaseq)
#> [1] "DGEList"
#> attr(,"package")
#> [1] "edgeR"
head(griss_melanoma_rnaseq$samples)
#> group lib.size norm.factors patient cell_type treatment
#> 195-13 MOCK 29907534 1.0629977 P1 TIBC MOCK
#> 195-14 MCM 26397322 0.9927768 P1 TIBC MCM
#> 195-19 MOCK 18194834 1.0077827 P2 PBMCB MOCK
#> 195-20 MCM 24282215 1.0041410 P2 PBMCB MCM
#> 197-11 MOCK 22628117 0.9522869 P1 PBMCB MOCK
#> 197-12 MCM 23319849 1.0115732 P1 PBMCB MCMAgain, the dataset can simply be added using
add_dataset. Here, we added an additional parameter to the
add_dataset call. Such additional parameters are treated as
additional dataset-level parameters.
# add the dataset
my_request <- add_dataset(request = my_request,
expression_values = griss_melanoma_rnaseq,
name = "RNA-seq",
type = "rnaseq_counts",
comparison_factor = "treatment",
comparison_group_1 = "MOCK",
comparison_group_2 = "MCM",
additional_factors = c("cell_type", "patient"),
# This adds the dataset-level parameter 'discrete_norm_function' to the request
discrete_norm_function = "TMM")
#> Converting expression data to string... (This may take a moment)
#> Conversion complete
my_request
#> ReactomeAnalysisRequestObject
#> Method = Camera
#> Parameters:
#> - max_missing_values: 0.5
#> Datasets:
#> - Proteomics (proteomics_int)
#> No parameters set.
#> - RNA-seq (rnaseq_counts)
#> discrete_norm_function: TMM
#> ReactomeAnalysisRequestSample annotations
Datasets can be passed as limma EList, edgeR
DGEList, any implementation of the Bioconductor
ExpressionSet, or simply a data.frame.
For the first three, sample annotations are simply read from the
respective slot. When supplying the expression values as a
data.frame, the sample_data parameter has to
be set using a data.frame where each row represents one
sample and each column one proptery. If the the sample_data
option is set while providing the expression data as an
EList, DGEList, or ExpressionSet,
the data in sample_data will be used instead of the sample
annotations in the expression data object.
Name
Each dataset has to have a name. This can be anything but has to be unique within one analysis request.
Type
The ReactomeAnalysisSystem supports different types of ’omics data.
To get a list of supported types, use the
get_reactome_data_types function:
get_reactome_data_types()
#> rnaseq_counts:
#> RNA-seq (raw counts)
#> Raw RNA-seq based read counts per gene (recommended).
#> rnaseq_norm:
#> RNA-seq (normalized)
#> log2 transformed, normalized RNA-seq based read counts per gene (f.e. RPKM, TPM)
#> proteomics_int:
#> Proteomics (intensity)
#> Intensity-based quantitative proteomics data (for example, iTRAQ/TMT or intensity-based label-free quantitation). Values must be log2 transformed.
#> proteomics_sc:
#> Proteomics (spectral counts)
#> Raw spectral-counts of label-free proteomics experiments
#> microarray_norm:
#> Microarray (normalized)
#> Normalized and log2 transformed microarray-based gene expression values.Defining the experimental design
Defining the experimental design for a
ReactomeAnalysisRequest is very simple. Basically, it only
takes three parameters:
comparison_factor: Name of the property within the sample data to usecomparison_group_1: The first group to comparecomparison_group_2: The second group to compare
The value set in comparison_factor must match a column
name in the sample data (either the slot in an Elist,
DGEList, or ExpressionSet object or in the
sample_data parameter).
Additionally, it is possible to define blocking factors. These are
supported by all methods that rely on linear models in the backend. Some
methods though might simply ignore this parameter. For more information
on whether a method supports blocking factors, please use
get_reactome_methods.
Blocking factors can simply be set additional_factors to
a vector of names. These should again reference properties (or columns)
in the sample data.
Submitting the request
Once the ReactomeAnalysisRequest is created, the
complete analysis can be run using
perform_reactome_analysis:
result <- perform_reactome_analysis(request = my_request, compress = F)
#> Submitting request to Reactome API...
#> Reactome Analysis submitted succesfully
#> Updating to new REACTOME version...
#> Converting dataset Proteomics...
#> Mapping identifiers...
#> Performing gene set analysis using Camera
#> Analysing dataset 'Proteomics' using Camera
#> Analysing dataset 'RNA-seq' using Camera
#> Creating REACTOME visualization
#> Retrieving result...Investigating the result
The result object is a ReactomeAnalysisResult S4 class
with several helper functions to access the data.
To retrieve the names of all available results (generally one per
dataset), use the names function:
For every dataset, different result types may be available. These can
be shown using the result_types function:
The Camera analysis method returns two types of results,
pathway-level data and gene- / protein-level fold changes.
A specific result can be retrieved using the get_result
method:
# retrieve the fold-change data for the proteomics dataset
proteomics_fc <- get_result(result, type = "fold_changes", name = "Proteomics")
head(proteomics_fc)
#> Identifier logFC AveExpr t P.Value adj.P.Val B
#> 1 Q14526 0.4937650 -3.346909 14.504637 1.489621e-10 8.846861e-07 13.962655
#> 2 Q6VY07 0.2981411 -3.330347 13.506868 4.246464e-10 1.260987e-06 13.071842
#> 3 P07093 1.7950301 -3.648968 12.297371 1.657293e-09 3.280887e-06 11.882006
#> 4 P10124 1.0758634 -3.436961 10.333180 1.946522e-08 2.890099e-05 9.650706
#> 5 P55210 0.5018522 -3.347932 9.513709 6.059130e-08 7.197035e-05 8.593544
#> 6 O43683 -0.4754083 -3.345551 -9.364551 7.507346e-08 7.431021e-05 8.392279Additionally, it is possible to directly merge the pathway level data
for all result sets using the pathways function:
combined_pathways <- pathways(result)
head(combined_pathways)
#> Name
#> R-HSA-1428517 Aerobic respiration and respiratory electron transport
#> R-HSA-611105 Respiratory electron transport
#> R-HSA-6799198 Complex I biogenesis
#> R-HSA-72649 Translation initiation complex formation
#> R-HSA-72662 Activation of the mRNA upon binding of the cap-binding complex and eIFs, and subsequent binding to 43S
#> R-HSA-72702 Ribosomal scanning and start codon recognition
#> Direction.Proteomics FDR.Proteomics PValue.Proteomics
#> R-HSA-1428517 Up 2.872443e-10 1.785314e-13
#> R-HSA-611105 Up 2.872443e-10 2.536374e-13
#> R-HSA-6799198 Up 1.092147e-08 1.446552e-11
#> R-HSA-72649 Down 7.888749e-08 1.393157e-10
#> R-HSA-72662 Down 1.493054e-07 3.295924e-10
#> R-HSA-72702 Down 1.642546e-07 4.961295e-10
#> NGenes.Proteomics av_foldchange.Proteomics sig.Proteomics
#> R-HSA-1428517 199 0.10614760 TRUE
#> R-HSA-611105 131 0.11620244 TRUE
#> R-HSA-6799198 64 0.12741142 TRUE
#> R-HSA-72649 57 -0.09502561 TRUE
#> R-HSA-72662 58 -0.09089517 TRUE
#> R-HSA-72702 57 -0.09160251 TRUE
#> Direction.RNA-seq FDR.RNA-seq PValue.RNA-seq NGenes.RNA-seq
#> R-HSA-1428517 Down 8.043354e-05 3.565708e-06 239
#> R-HSA-611105 Down 5.129250e-04 3.279983e-05 149
#> R-HSA-6799198 Down 4.485027e-03 3.747786e-04 68
#> R-HSA-72649 Down 8.107325e-02 1.463072e-02 58
#> R-HSA-72662 Down 1.085202e-01 2.154226e-02 59
#> R-HSA-72702 Down 7.988876e-02 1.429160e-02 58
#> av_foldchange.RNA-seq sig.RNA-seq
#> R-HSA-1428517 -0.14949659 TRUE
#> R-HSA-611105 -0.15378463 TRUE
#> R-HSA-6799198 -0.15727714 TRUE
#> R-HSA-72649 -0.11783700 FALSE
#> R-HSA-72662 -0.09156838 FALSE
#> R-HSA-72702 -0.11805469 FALSEVisualising results
The ReactomeGSA package includes several basic plotting functions to
visualise the pathway results. For comparative gene set analysis like
the one presented here, two functions are available:
plot_correlations and plot_volcano.
plot_correlations can be used to quickly assess how
similar two datasets are on the pathway level:
plot_correlations(result)
#> Comparing 1 vs 2
#> [[1]]
#> Warning: Removed 290 rows containing missing values or values outside the scale range
#> (`geom_point()`).
Individual datasets can further be visualised using volcano plots of the pathway data:
Finally, it is possible to view the result as a heatmap:
plot_heatmap(result) +
# reduce text size to create a better HTML rendering
ggplot2::theme(text = ggplot2::element_text(size = 6))
By default, only 30 pathways are shown in the heatmap. It is also possible to easily manually select pathways of interest to plot:
# get the data ready to plot with ggplot2
plot_data <- plot_heatmap(result, return_data = TRUE)
# select the pathways of interest - here all pathways
# with "Interleukin" in their name
interleukin_pathways <- grepl("Interleukin", plot_data$Name)
interesting_data <- plot_data[interleukin_pathways, ]
# create the heatmap
ggplot2::ggplot(interesting_data, ggplot2::aes(x = dataset, y = Name, fill = direction)) +
ggplot2::geom_tile() +
ggplot2::scale_fill_brewer(palette = "RdYlBu") +
ggplot2::labs(x = "Dataset", fill = "Direction") +
ggplot2::theme(text = ggplot2::element_text(size = 6))
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] Biobase_2.67.0 BiocGenerics_0.53.0 ReactomeGSA.data_1.19.0
#> [4] Seurat_5.1.0 SeuratObject_5.0.2 sp_2.1-4
#> [7] ReactomeGSA_1.21.0 edgeR_4.4.0 limma_3.63.0
#> [10] 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 spatstat.geom_3.3-3
#> [76] RcppAnnoy_0.0.22 ggrepel_0.9.6 RANN_2.6.2
#> [79] pillar_1.9.0 stringr_1.5.1 spam_2.11-0
#> [82] RcppHNSW_0.6.0 later_1.3.2 splines_4.4.1
#> [85] dplyr_1.1.4 lattice_0.22-6 survival_3.7-0
#> [88] deldir_2.0-4 tidyselect_1.2.1 locfit_1.5-9.10
#> [91] maketools_1.3.1 miniUI_0.1.1.1 pbapply_1.7-2
#> [94] knitr_1.48 gridExtra_2.3 scattermore_1.2
#> [97] xfun_0.48 statmod_1.5.0 matrixStats_1.4.1
#> [100] stringi_1.8.4 lazyeval_0.2.2 yaml_2.3.10
#> [103] evaluate_1.0.1 codetools_0.2-20 tibble_3.2.1
#> [106] cli_3.6.3 uwot_0.2.2 xtable_1.8-4
#> [109] reticulate_1.39.0 munsell_0.5.1 jquerylib_0.1.4
#> [112] Rcpp_1.0.13 globals_0.16.3 spatstat.random_3.3-2
#> [115] png_0.1-8 spatstat.univar_3.0-1 parallel_4.4.1
#> [118] ggplot2_3.5.1 prettyunits_1.2.0 dotCall64_1.2
#> [121] bitops_1.0-9 listenv_0.9.1 viridisLite_0.4.2
#> [124] scales_1.3.0 ggridges_0.5.6 crayon_1.5.3
#> [127] leiden_0.4.3.1 purrr_1.0.2 rlang_1.1.4
#> [130] cowplot_1.1.3