diff --git a/.Rbuildignore b/.Rbuildignore index 0e5321c..2e3efa5 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -8,3 +8,4 @@ ^_pkgdown\.yml$ ^docs$ ^pkgdown$ +^\.github$ diff --git a/.github/.gitignore b/.github/.gitignore new file mode 100644 index 0000000..2d19fc7 --- /dev/null +++ b/.github/.gitignore @@ -0,0 +1 @@ +*.html diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml index cea5906..bfc9f4d 100644 --- a/.github/workflows/pkgdown.yaml +++ b/.github/workflows/pkgdown.yaml @@ -2,7 +2,7 @@ # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help on: push: - branches: main + branches: [main, master] pull_request: release: types: [published] diff --git a/DESCRIPTION b/DESCRIPTION index 4a6edcd..551fec2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: leapR Title: Layered enrichment analysis of pathways R -Version: 0.99.7 +Version: 0.99.9 Authors@R: c( person("Sara", "Gosline", email = "sara.gosline@pnnl.gov", role = c('aut',"cre"), comment = c(ORCID = "0000-0002-6534-4774")), person("Jason", "McDermott", email = "jason.mcdermott@pnnl.gov", role = "aut"), @@ -17,7 +17,8 @@ biocViews: Proteomics, Pathways, GeneExpression, - Transcriptomics + Transcriptomics, + Software Imports: stats, gplots, @@ -30,7 +31,8 @@ Imports: stringr, tidyr, SummarizedExperiment, - BiocStyle + BiocStyle, + BiocFileCache Suggests: knitr, rmarkdown, @@ -38,3 +40,4 @@ Suggests: VignetteBuilder: knitr License: MIT + file LICENSE Config/testthat/edition: 3 +URL: https://pnnl.github.io/leapR/ diff --git a/R/calcTTest.R b/R/calcTTest.R index 095c35c..ca0df02 100644 --- a/R/calcTTest.R +++ b/R/calcTTest.R @@ -18,11 +18,15 @@ #' @examples #' #' library(leapR) +#' library(BiocFileCache) +#' +#' path <- tools::R_user_dir("leapR", which = "cache") +#' bfc <- BiocFileCache(path, ask = FALSE) +#' #' url <- "https://api.figshare.com/v2/file/download/56536214" -#' tdata <- download.file(url,method='libcurl',destfile='transData.rda') -#' load('transData.rda') -#' p <- file.remove("transData.rda") -#' +#' tc <- bfcadd(bfc, "tdat", fpath = url) +#' load(tc) +#' #' # read in the pathways #' data("ncipid") #' diff --git a/R/cluster_enrichment.R b/R/cluster_enrichment.R index d1f3120..7e76adf 100644 --- a/R/cluster_enrichment.R +++ b/R/cluster_enrichment.R @@ -19,12 +19,14 @@ #' @export #' @examples #' library(leapR) -#' -#' # read in the example transcriptomic data +#' library(BiocFileCache) +#' +#' path <- tools::R_user_dir("leapR", which = "cache") +#' bfc <- BiocFileCache(path, ask = FALSE) +#' #' url <- "https://api.figshare.com/v2/file/download/56536214" -#' tdata <- download.file(url,method='libcurl',destfile='transData.rda') -#' load('transData.rda') -#' p <- file.remove("transData.rda") +#' tc <- bfcadd(bfc, "tdat", fpath = url) +#' load(tc) #' #' # read in the pathways #' data("ncipid") diff --git a/R/combine_omics.R b/R/combine_omics.R index 101d595..78dbbcc 100644 --- a/R/combine_omics.R +++ b/R/combine_omics.R @@ -18,23 +18,21 @@ #' #' @examples #' library(leapR) -#' url <- 'https://api.figshare.com/v2/file/download/56536217' -#' -#' pdata <- download.file(url,method='libcurl',destfile='protData.rda') -#' load('protData.rda') -#' p <- file.remove("protData.rda") -#' +#' library(BiocFileCache) +#' path <- tools::R_user_dir("leapR", which = "cache") +#' bfc <- BiocFileCache(path, ask = FALSE) +#' +#' url <- "https://api.figshare.com/v2/file/download/56536217" +#' pc <- bfcadd(bfc, "pdat", fpath = url) +#' load(pc) +#' #' url <- "https://api.figshare.com/v2/file/download/56536214" -#' tdata <- download.file(url,method='libcurl',destfile='transData.rda') -#' load('transData.rda') -#' p <- file.remove("transData.rda") -#' -#' url <- 'https://api.figshare.com/v2/file/download/56536211' -#' phdata<-download.file(url,method='libcurl',destfile = 'phosData.rda') -#' #phosphodata<-read.csv("phdata",check.names=FALSE,row.names=1) -#' load('phosData.rda') -#' p <- file.remove('phosData.rda')# read in the example protein data -#' +#' tc <- bfcadd(bfc, "tdat", fpath = url) +#' load(tc) +#' +#' url <- "https://api.figshare.com/v2/file/download/56536211" +#' phc <- bfcadd(bfc, "phdat", fpath = url) +#' load(phc) #' #' # merge the three datasets by rows and add prefix tags for #' # different omics types diff --git a/R/enrichment_in_groups.R b/R/enrichment_in_groups.R index b82cc45..4d92940 100644 --- a/R/enrichment_in_groups.R +++ b/R/enrichment_in_groups.R @@ -18,6 +18,7 @@ #' to log your data before calling. NOTE: if you do not call `suppressWarnings` then #' the KS test will warn you about ties. #' @param minsize minimum size of set +#' @param log_transformed Set to TRUE if data is already log-transformed #' @param mapping_column column name of mapping identifiers #' @param abundance_column columns mapping abundance, either in the `assay` #' matrix or `rowData` @@ -114,7 +115,7 @@ enrichment_in_groups <- function(geneset, names(backvals) <- backlist#[-group_ind] in_back <- length(backvals) - outgroup_mean = mean(backvals[-group_ind], na.rm = T) + outgroup_mean = mean(backvals[-group_ind], na.rm = TRUE) in_path <- length(in_group) #how many left after na.rm if ((in_path > minsize) & (any(!is.na(in_path))) & diff --git a/R/leapR-package.R b/R/leapR-package.R index ad09996..f82de0d 100644 --- a/R/leapR-package.R +++ b/R/leapR-package.R @@ -143,14 +143,14 @@ #' \cr #' @examples #' library(leapR) -#' -#' # read in the example abundance data -#' # read in the example transcriptomic data -#' tdata <- download.file("https://api.figshare.com/v2/file/download/56536214", -#' method='libcurl',destfile='transData.rda') -#' load('transData.rda') -#' p <- file.remove("transData.rda") -#' +#' library(BiocFileCache) +#' +#' path <- tools::R_user_dir("leapR", which = "cache") +#' bfc <- BiocFileCache(path, ask = FALSE) +#' +#' url <- "https://api.figshare.com/v2/file/download/56536214" +#' tc <- bfcadd(bfc, "tdat", fpath = url) +#' load(tc) #' # read in the pathways #' data("ncipid") #' diff --git a/R/leapR.R b/R/leapR.R index 05318c4..4b53700 100644 --- a/R/leapR.R +++ b/R/leapR.R @@ -140,13 +140,14 @@ #' \cr #' @examples #' library(leapR) -#' -#' # read in the example abundance data -#' # read in the example transcriptomic data -#' tdata <- download.file("https://api.figshare.com/v2/file/download/56536214", -#' method='libcurl',destfile='transData.rda') -#' load('transData.rda') -#' p <- file.remove("transData.rda") +#' library(BiocFileCache) +#' +#' path <- tools::R_user_dir("leapR", which = "cache") +#' bfc <- BiocFileCache(path, ask = FALSE) +#' +#' url <- "https://api.figshare.com/v2/file/download/56536214" +#' tc <- bfcadd(bfc, "tdat", fpath = url) +#' load(tc) #' #' # read in the pathways #' data("ncipid") diff --git a/_pkgdown.yml b/_pkgdown.yml index 69b0979..2ecd09d 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -1,4 +1,4 @@ -url: https://pnnl-github.io/leapR +url: https://pnnl.github.io/leapR/ template: bootstrap: 5 diff --git a/docs/404.html b/docs/404.html deleted file mode 100644 index a4be44e..0000000 --- a/docs/404.html +++ /dev/null @@ -1,82 +0,0 @@ - - -
- - - - -Copyright 2025 Battelle Memorial Institute - - -Redistribution and use in source and binary forms, with or without -modification, are permitted provided that the following conditions -are met: - - -1. Redistributions of source code must retain the above copyright -notice, this list of conditions and the following disclaimer. - - -2. Redistributions in binary form must reproduce the above copyright -notice, this list of conditions and the following disclaimer in the -documentation and/or other materials provided with the distribution. - - -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, -INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, -BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; -LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT -LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN -ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -POSSIBILITY OF SUCH DAMAGE. -- -
examples.RmdA sample data set is included that is from the CPTAC study of 169
-ovarian tumors. We include the dataset as a object, containing three
-assays (transcriptomics, global proteomics, and phosphoproteomics) to
-enable interoperability with other tools, and store example file as
-rda on Figshare
-as example.
This data can be loaded as follows:
-
-url <- "https://api.figshare.com/v2/file/download/56536217"
-pdata <- download.file(url, method = "libcurl", destfile = "protData.rda")
-# as.matrix()
-load("protData.rda")
-
-p <- file.remove("protData.rda")
-
-url <- "https://api.figshare.com/v2/file/download/56536214"
-tdata <- download.file(url, method = "libcurl", destfile = "transData.rda")
-load("transData.rda")
-p <- file.remove("transData.rda")
-
-url <- "https://api.figshare.com/v2/file/download/56536211"
-phdata <- download.file(url, method = "libcurl", destfile = "phosData.rda")
-load("phosData.rda")
-p <- file.remove("phosData.rda")We also have local data we can load
- -We will compare the ability of transcriptomics, proteomics, and
-phosphoproteomics to inform about differences between short and long
-surviving patient groups. In addition to other methods, we also employ a
-calcTTest function that takes two sets of samples from the
-SummarizedExperiment object and computes the t-test between
-them. The results are then stored in the rowData of the
-same object, so that they can be used for enrichment later on.
This spans the multiple enrichment methods in leapR and also includes -multi-omics.
-The resulting heatmap is presented as Figure 2 in the paper.
-
-# load the single omic and multi-omic pathway databases
-data("krbpaths")
-data("mo_krbpaths")
-
-# comparison enrichment in transcriptional data
-transdata.comp.enrichment.svl <- leapR::leapR(
- geneset = krbpaths,
- enrichment_method = "enrichment_comparison",
- eset = tset,
- assay_name = "transcriptomics",
- primary_columns = shortlist,
- secondary_columns = longlist
-)
-
-# comparison enrichment in proteomics data
-# this is the same code used above, just repeated here for clarity
-protdata.comp.enrichment.svl <- leapR::leapR(
- geneset = krbpaths,
- enrichment_method = "enrichment_comparison",
- eset = pset,
- assay_name = "proteomics",
- primary_columns = shortlist,
- secondary_columns = longlist
-)
-
-# comparison enrichment in phosphoproteomics data
-phosphodata.comp.enrichment.svl <- leapR::leapR(
- geneset = krbpaths,
- enrichment_method = "enrichment_comparison",
- eset = phset,
- assay_name = "phosphoproteomics",
- primary_columns = shortlist,
- secondary_columns = longlist, id_column = "hgnc_id"
-)
-
-
-# set enrichment in transcriptomics data
-# perform the comparison t-test
-tset <- leapR::calcTTest(tset, assay_name = "transcriptomics",
- shortlist, longlist)
-
-
-## now we need to run enrichment in sets with target list, not eset
-transdata.set.enrichment.svl <- leapR::leapR(
- geneset = krbpaths,
- eset = tset,
- assay_name = "transcriptomics",
- enrichment_method = "enrichment_in_sets",
- primary_columns = "pvalue",
- greaterthan = FALSE, threshold = 0.05
-)
-
-
-pset <- leapR::calcTTest(pset, assay_name = "proteomics",
- shortlist, longlist)
-
-protdata.set.enrichment.svl <- leapR::leapR(
- geneset = krbpaths,
- eset = pset,
- assay_name = "proteomics",
- enrichment_method = "enrichment_in_sets",
- primary_columns = "pvalue",
- greaterthan = FALSE, threshold = 0.05
-)
-
-
-phset <- leapR::calcTTest(phset, assay_name = "phosphoproteomics",
- shortlist, longlist)
-
-phosphodata.set.enrichment.svl <- leapR::leapR(
- geneset = krbpaths,
- enrichment_method = "enrichment_in_sets",
- id_column = "hgnc_id",
- assay_name = "phosphoproteomics",
- eset = phset, primary_columns = "pvalue",
- greaterthan = FALSE, threshold = 0.05
-)
-
-# order enrichment in transcriptomics data
-transdata.order.enrichment.svl <- leapR::leapR(
- geneset = krbpaths,
- enrichment_method = "enrichment_in_order",
- eset = tset,
- assay_name = "transcriptomics",
- primary_columns = "difference"
-)
-
-# order enrichment in proteomics data
-protdata.order.enrichment.svl <- leapR::leapR(
- geneset = krbpaths,
- enrichment_method = "enrichment_in_order",
- eset = pset,
- assay_name = "proteomics",
- primary_columns = "difference"
-)
-
-# order enrichment in phosphoproteomics data
-
-
-phosphodata.order.enrichment.svl <- leapR::leapR(
- geneset = krbpaths,
- enrichment_method = "enrichment_in_order",
- id_column = "hgnc_id",
- method = 'ztest',
- eset = phset,
- assay_name = "phosphoproteomics",
- primary_columns = "difference"
-)
-
-# correlation difference in transcriptomics data
-transdata.corr.enrichment.svl <- leapR::leapR(
- geneset = krbpaths,
- enrichment_method = "correlation_comparison",
- eset = tset,
- assay_name = "transcriptomics",
- primary_columns = shortlist,
- secondary_columns = longlist
-)
-# correlation difference in proteomics data
-protdata.corr.enrichment.svl <- leapR::leapR(
- geneset = krbpaths,
- enrichment_method = "correlation_comparison",
- eset = pset,
- assay_name = "proteomics",
- primary_columns = shortlist,
- secondary_columns = longlist
-)
-# correlation difference in phosphoproteomics data
-phosphodata.corr.enrichment.svl <- leapR::leapR(
- geneset = krbpaths,
- enrichment_method = "correlation_comparison",
- eset = phset,
- assay_name = "phosphoproteomics",
- primary_columns = shortlist,
- secondary_columns = longlist, id_column = "hgnc_id"
-)
-
-# combine the omics data into one with prefix tags
-comboset <- leapR::combine_omics(list(pset, phset, tset),
- c(NA, "hgnc_id", NA))
-
-# comparison enrichment for combodata
-# when we use expression set, we do not need to use the mo_krbpaths
-#since the id mapping column is used
-combodata.enrichment.svl <- leapR::leapR(
- geneset = krbpaths, # mo_krbpaths,
- enrichment_method = "enrichment_comparison",
- eset = comboset,
- assay_name = "combined",
- primary_columns = shortlist,
- secondary_columns = longlist, id_column = "id"
-)
-
-
-# set enrichment in combo data
-# perform the comparison t test
-comboset <- leapR::calcTTest(comboset,
- assay_name = "combined",
- shortlist, longlist)
-
-combodata.set.enrichment.svl <- leapR::leapR(
- geneset = krbpaths,
- enrichment_method = "enrichment_in_sets",
- eset = comboset, primary_columns = "pvalue",
- assay_name = "combined",
- id_column = "id",
- greaterthan = FALSE, threshold = 0.05
-)
-
-# order enrichment in combo data
-combodata.order.enrichment.svl <- leapR::leapR(
- geneset = krbpaths,
- enrichment_method = "enrichment_in_order",
- assay_name = "combined",
- eset = comboset, primary_columns = "difference",
- id_column = "id"
-)
-
-# correlation difference in combo data
-combodata.corr.enrichment.svl <- leapR::leapR(
- geneset = krbpaths,
- enrichment_method = "correlation_comparison",
- eset = comboset,
- assay_name = "combined",
- primary_columns = shortlist,
- id_column = "id",
- secondary_columns = longlist
-)
-
-
-# now take all these results and combine them into one figure
-all_results <- list(
- transdata.comp.enrichment.svl,
- protdata.comp.enrichment.svl,
- phosphodata.comp.enrichment.svl,
- combodata.enrichment.svl,
- transdata.set.enrichment.svl,
- protdata.set.enrichment.svl,
- phosphodata.set.enrichment.svl,
- combodata.set.enrichment.svl,
- transdata.order.enrichment.svl,
- protdata.order.enrichment.svl,
- phosphodata.order.enrichment.svl,
- combodata.order.enrichment.svl,
- transdata.corr.enrichment.svl,
- protdata.corr.enrichment.svl,
- phosphodata.corr.enrichment.svl,
- combodata.corr.enrichment.svl
-)
-
-pathways_of_interest <- c(
- "KEGG_APOPTOSIS",
- "KEGG_CELL_CYCLE",
- "KEGG_ERBB_SIGNALING_PATHWAY",
- "KEGG_FOCAL_ADHESION",
- "KEGG_INSULIN_SIGNALING_PATHWAY",
- "KEGG_MAPK_SIGNALING_PATHWAY",
- "KEGG_MISMATCH_REPAIR",
- "KEGG_MTOR_SIGNALING_PATHWAY",
- "KEGG_OXIDATIVE_PHOSPHORYLATION",
- "KEGG_P53_SIGNALING_PATHWAY",
- "KEGG_PATHWAYS_IN_CANCER",
- "KEGG_PROTEASOME",
- "KEGG_RIBOSOME",
- "KEGG_VEGF_SIGNALING_PATHWAY",
- "KEGG_WNT_SIGNALING_PATHWAY"
-)
-
-
-results.frame <- data.frame(
- pathway = pathways_of_interest,
- td.comp = all_results[[1]][pathways_of_interest, "BH_pvalue"] < 0.05,
- pd.comp = all_results[[2]][pathways_of_interest, "BH_pvalue"] < 0.05,
- fd.comp = all_results[[3]][pathways_of_interest, "BH_pvalue"] < 0.05,
- cd.comp = all_results[[4]][pathways_of_interest, "BH_pvalue"] < 0.05,
- td.set = all_results[[5]][pathways_of_interest, "BH_pvalue"] < 0.05,
- pd.set = all_results[[6]][pathways_of_interest, "BH_pvalue"] < 0.05,
- fd.set = all_results[[7]][pathways_of_interest, "BH_pvalue"] < 0.05,
- cd.set = all_results[[8]][pathways_of_interest, "BH_pvalue"] < 0.05,
- td.order = all_results[[9]][pathways_of_interest, "BH_pvalue"] < 0.05,
- pd.order = all_results[[10]][pathways_of_interest, "BH_pvalue"] < 0.05,
- fd.order = all_results[[11]][pathways_of_interest, "BH_pvalue"] < 0.05,
- cd.order = all_results[[12]][pathways_of_interest, "BH_pvalue"] < 0.05,
- td.corr = all_results[[13]][pathways_of_interest, "BH_pvalue"] < 0.05,
- pd.corr = all_results[[14]][pathways_of_interest, "BH_pvalue"] < 0.05,
- fd.corr = all_results[[15]][pathways_of_interest, "BH_pvalue"] < 0.05,
- cd.corr = all_results[[16]][pathways_of_interest, "BH_pvalue"] < 0.05
-)
-
-results.frame.or <- data.frame(
- pathway = pathways_of_interest,
- td.comp = all_results[[1]][pathways_of_interest, "oddsratio"],
- pd.comp = all_results[[2]][pathways_of_interest, "oddsratio"],
- fd.comp = all_results[[3]][pathways_of_interest, "oddsratio"],
- cd.comp = all_results[[4]][pathways_of_interest, "oddsratio"],
- td.set = log(all_results[[5]][pathways_of_interest, "oddsratio"], 2),
- pd.set = log(all_results[[6]][pathways_of_interest, "oddsratio"], 2),
- fd.set = log(all_results[[7]][pathways_of_interest, "oddsratio"], 2),
- cd.set = log(all_results[[8]][pathways_of_interest, "oddsratio"], 2),
- td.order = all_results[[9]][pathways_of_interest, "oddsratio"],
- pd.order = all_results[[10]][pathways_of_interest, "oddsratio"],
- fd.order = all_results[[11]][pathways_of_interest, "oddsratio"],
- cd.order = all_results[[12]][pathways_of_interest, "oddsratio"],
- td.corr = all_results[[13]][pathways_of_interest, "oddsratio"],
- pd.corr = all_results[[14]][pathways_of_interest, "oddsratio"],
- fd.corr = all_results[[15]][pathways_of_interest, "oddsratio"],
- cd.corr = all_results[[16]][pathways_of_interest, "oddsratio"]
-)
-
-rownames(results.frame) <- results.frame[, 1]
-rownames(results.frame.or) <- results.frame.or[, 1]
-results.frame.sig <- results.frame[, 2:17] * results.frame.or[, 2:17]
-
-heatmap.2(as.matrix(results.frame.sig[, c(1:4, 9:16)]), Colv = NA,
- trace = "none", breaks = c(-1, -.1, -0.0001, 0, 0.1, 1),
- col = c("blue", "lightblue", "grey", "pink", "red"),
- dendrogram = "none")
Figure 3. An application of the KSEA-like approach in leapR as -applied to our example data. In this example we are looking for known -substrate sets of kinases (from Phosphosite Plus) that are enriched in -the short vs long comparison of phosphopeptides.
-
-# this comparison of abundance in substrates between case and control
-# is lopsided in the sense that phosphorylation levels were previously
-# reported to be overall higher in the short survivors. Thus the
-# results are not terribly interesting (all kinases are in the same
-# direction)
-phosphodata.ksea.comp.svl <- leapR::leapR(
- geneset = kinasesubstrates,
- enrichment_method = "enrichment_comparison",
- eset = phset,
- assay_name = "phosphoproteomics",
- primary_columns = shortlist, secondary_columns = longlist
-)
-
-
-# thus for the example we'll look at correlation between known substrates
-# in the case v control conditions
-phosphodata.ksea.corr.svl <- leapR::leapR(
- geneset = kinasesubstrates,
- enrichment_method = "correlation_comparison",
- eset = phset,
- assay_name = "phosphoproteomics",
- primary_columns = shortlist,
- secondary_columns = longlist
-)
-
-# for the example we are using an UNCORRECTED PVALUE
-# which will allow us to plot more values, but
-# for real applications it's necessary to use the
-# CORRECTED PVALUE
-
-# here are all the kinases *significant (*uncorrected) from the analysis
-or <- order(phosphodata.ksea.corr.svl[, "pvalue"])
-ksea_result <- phosphodata.ksea.corr.svl[or, ][1:9, ]
-ksea_cols <- rep("grey", 9)
-ksea_cols[which(ksea_result[, "oddsratio"] > 0)] <- "black"
-
-# plot left panel: correlation significance of top most significant kinases
-barplot(ksea_result[, "oddsratio"],
- horiz = TRUE, xlim = c(-1, 0.5),
- names.arg = rownames(ksea_result), las = 1, col = ksea_cols
-)
-
-# plot right panel: abundance comparison results of the same kinases
-barplot(phosphodata.ksea.comp.svl[rownames(ksea_result), "oddsratio"],
- horiz = TRUE, names.arg = rownames(ksea_result), las = 1, col = "black"
-)
leapR.RmdThis is intended to be a short introduction to the leapR
-package. First we need to load the required libraries:
-# install from bioconductor
-if (!require(BiocManager)) {
- install.packages('BiocManager')
- BiocManager::install('leapR')
-}Dataset - an expression dataset, contained in the -Bioconductor object, that at the bare minimum has a matrix of components -(rows) measured in the same system under multiple different conditions -(columns)
-Component - the things being measured, genes, proteins, -methylation site, phosphosite, etc. For functional (currently) the -component must be associated with a gene name. That is, there’s not -currently a way to calculate pathway enrichment using lipids.
-Pathway - a set of components that works together to -accomplish something or are related to each other in some other way. -This includes classic signaling and metabolic pathways, but also -molecular function and localization categories and other groups of -related components, like genome location, conservation, etc.
-Condition - a sample where the treatment, environmental -conditions, patient, time point or some combination of those is -varied.
-The overall idea for functional enrichment is to determine which -pathways are statistically over-represented in one group versus another, -display statistically differential abundance from one group to another, -or are statistically differentially distributed in a ranked list based -on the abundance of one sample. Each of these purposes has a different -underlying statistical test (or family of tests) and the results of each -can be interpreted in somewhat different ways. The purpose of this -vignette is to give the user a very brief introduction on how to use the -package, not to discuss the underlying statistical choices that need to -be made when analyzing such data.
-There are a number of caveats (probably non-exhaustive) with doing -this kind of analysis.
-One important point is to use data that’s been normalized in a -particular way to do these analyses. Data here has been normalized as a -Z score by row (gene/protein/etc.). So, for each row, calculate -the mean and standard deviation across all the conditions (columns) and -then express as a Z score.
-Here’s why. All high-throughput technologies (microarray, RNAseq, -MS-assisted proteomics, metabolomics, lipidomics, etc.) suffer from the -same limitation. The detectability of each molecule being detected -(protein, RNA, etc.) is different and, in general, it’s impossible to -accurately determine how detectable each one is. The multi-omic -functional enrichment process lumps together measurements from different -components (proteins, genes, etc.) to summarize a pathway. If the -component measurements aren’t directly comparable (they aren’t) then -this can and will introduce significant systematic errors and won’t -produce the results you’re looking for. Careful consideration must be -given that the results of the analysis reflect the question being asked -and that the normalization method hasn’t obscured the desired -results.
-The background of comparison for functional enrichment is always -important, but it mainly impacts the Fisher’s exact tests in the -examples below. The background answers the question: “My functional -group of interest is statistically enriched relative to what?” For -Fisher’s exact tests this is critical. Generally, it is best to compare -enrichment against the components observed in the data (the experiment’s -universe) rather than the universe of all possible components. For -example, a proteomics dataset from plasma may contain a limited set of -proteins compared with all possible human proteins; using the observed -proteins as the background usually yields more meaningful results. Using -all possible proteins will result in substantially different -findings.
-When testing the statistical significance of differences in a lot of -pathways it’s necessary to correct for multiple hypotheses. This -essentially accounts for the possibility you might see SOMETHING -significant by chance if you just test enough things- so it moves p -values in a less significant direction. The more things you test, the -greater this move will be. So pathway databases with lots of pathways -are affected more by this correction, making it harder to get a -significant result (which is a good thing actually).
-Two ‘databases’ (organized text files) are included for pathways. The -example is taken from the NCI’s Pathway Interaction Database (PID) and -covers signaling pathways in human - but is no longer being actively -maintained. They can be loaded as follows:
-
-data(ncipid)The identifiers (gene names, e.g.) for the data input MUST match the -identifiers used in the pathway database. The two included human -databases use the HGNC-approved gene names. Which means your data has to -use the same identifiers.
-A sample data set is included that is from the CPTAC study of 169
-ovarian tumors. We include the dataset as a object, containing three
-assays (transcriptomics, global proteomics, and phosphoproteomics) to
-enable interoperability with other tools, and store example file as
-rda on Figshare
-as example.
This data can be loaded as follows:
-
-# Inspect the `pset` SummarizedExperiment.
-str(pset)
-#> Formal class 'SummarizedExperiment' [package "SummarizedExperiment"] with 5 slots
-#> ..@ colData :Formal class 'DFrame' [package "S4Vectors"] with 6 slots
-#> .. .. ..@ rownames : chr [1:174] "TCGA-09-1664" "TCGA-09-2056" "TCGA-13-1404" "TCGA-13-1409" ...
-#> .. .. ..@ nrows : int 174
-#> .. .. ..@ elementType : chr "ANY"
-#> .. .. ..@ elementMetadata: NULL
-#> .. .. ..@ metadata : list()
-#> .. .. ..@ listData : Named list()
-#> ..@ assays :Formal class 'SimpleAssays' [package "SummarizedExperiment"] with 1 slot
-#> .. .. ..@ data:Formal class 'SimpleList' [package "S4Vectors"] with 4 slots
-#> .. .. .. .. ..@ listData :List of 1
-#> .. .. .. .. .. ..$ proteomics: num [1:18632, 1:174] -4.0649 -0.1398 -0.0366 0.768 0.2437 ...
-#> .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
-#> .. .. .. .. .. .. .. ..$ : chr [1:18632] "C9orf152" "ELMO2" "RPS11" "CREB3L1" ...
-#> .. .. .. .. .. .. .. ..$ : chr [1:174] "TCGA-09-1664" "TCGA-09-2056" "TCGA-13-1404" "TCGA-13-1409" ...
-#> .. .. .. .. ..@ elementType : chr "ANY"
-#> .. .. .. .. ..@ elementMetadata: NULL
-#> .. .. .. .. ..@ metadata : list()
-#> ..@ NAMES : chr [1:18632] "C9orf152" "ELMO2" "RPS11" "CREB3L1" ...
-#> ..@ elementMetadata:Formal class 'DFrame' [package "S4Vectors"] with 6 slots
-#> .. .. ..@ rownames : NULL
-#> .. .. ..@ nrows : int 18632
-#> .. .. ..@ elementType : chr "ANY"
-#> .. .. ..@ elementMetadata: NULL
-#> .. .. ..@ metadata : list()
-#> .. .. ..@ listData : Named list()
-#> ..@ metadata : list()
-dim(SummarizedExperiment::assay(pset, "proteomics"))
-#> [1] 18632 174
-head(rownames(pset))
-#> [1] "C9orf152" "ELMO2" "RPS11" "CREB3L1" "PNMA1" "MMP2"
-head(colnames(pset))
-#> [1] "TCGA-09-1664" "TCGA-09-2056" "TCGA-13-1404" "TCGA-13-1409" "TCGA-13-1410"
-#> [6] "TCGA-13-1482"We also include some groups of patients to compare stored as R data -objects:
-
-data(shortlist)
-data(longlist)
-
-## columns that we want to use for results
-
-cols_to_display <- c("ingroup_n", "outgroup_n", "background_n",
- "pvalue", "BH_pvalue")The data are now loaded and ready to go through some of the -examples.
-We include five examples of how to use this tool, depending on the -analysis at hand. ## Comparison of one condition/group versus another -condition/group.
-There are a number of ways to do this. I generally use a simple -approach which assesses the statistical difference in distributions -between the abundance values from all the members of a pathway in all -the group members from one group with those from the other group using a -t test.
-This is a ‘bag of values’ approach and it does not pay attention to -the relationships between values in different groups (i.e. that each -group has measurements for the same component). There are likely issues -that rise because of this and caveats associated with it. However, it -works fairly well.
-In this example we are assessing the enrichment of pathways in a
-group of short surviving patients versus in a group of long
-surviving
-patients. We can also do a single patient-to-patient comparison or
-compare a single patient to a group of patients.
Better corrected p-values are more enriched. However, you can get
-good p-values when the algorithm only considers a limited number of
-components from a pathway. That is, the pathway may have 30 members and
-the p-value is coming from values from just 3 members. You can look at
-the ingroup_n column from the result matrix to see this
-(and screen out if desired).
It is VERY important to also consider the effect size. That is, the
-difference between the mean of one group and the mean of the other
-group. If there are large numbers of components in the pathway being
-compared it is relatively easy to get a significant p value with small
-effect size. Though this may be a real difference it is often not as
-interesting as a smaller group with worse p value and greater effect
-size. You can look at the effect size by comparing the
-ingroup_mean and outgroup_mean columns.
-# in this example we lump a bunch of patients together (the 'short survivors')
-# and compare them to another group (the 'long survivors')
-
-### using enrichment_wrapper function
-protdata.enrichment.svl <- leapR::leapR(
- geneset = ncipid,
- enrichment_method = "enrichment_comparison",
- eset = pset,
- assay_name = "proteomics",
- primary_columns = shortlist,
- secondary_columns = longlist
-)
-
-or <- order(unlist(protdata.enrichment.svl[, "pvalue"]))
-rmarkdown::paged_table(protdata.enrichment.svl[or, cols_to_display])
-# another application is to compare just one patient against another
-# (this would be the equivalent of comparing one time point to another)
-
-### using enrichment_wrapper function
-protdata.enrichment.svl.ovo <- leapR::leapR(
- geneset = ncipid,
- enrichment_method = "enrichment_comparison",
- eset = pset,
- assay_name = "proteomics",
- primary_columns = shortlist[1],
- secondary_columns = longlist[1]
-)
-or <- order(unlist(protdata.enrichment.svl.ovo[, "pvalue"]))
-rmarkdown::paged_table(protdata.enrichment.svl.ovo[or, cols_to_display])When we only compare one sample to another, we get no enriched -pathways.
-For this test I use Fisher’s exact which is a simple comparison of -the overlap of two sets (think of it like a statistical Venn diagram -with two groups). It’s also referred to as a hypergeometric test.
-Caveat 1. Fisher’s exact does not consider abundance values -but only lists of components. Generally this requires some separation of -a group of interest using differential expression, module membership -(from a network for example), or some other method.
-Caveat 2. The choice of background for comparison can make a -big difference on outcome. For example, in a proteomics experiment where -you’re looking at enrichment in a group of highly differentially -expressed proteins, you could choose to use all possible proteins as a -background, or you could use just those proteins that were observed by -proteomics (generally a much more limited set). The second option is -generally the best since the first options will result in (partly to -mostly) functions that are enriched in proteins that are seen in -proteomics. That is, the most abundant proteins, which is generally not -the desired outcome.
-In the example below I construct a genelist of interest using a -simple abundance threshold on the data then use a background of all the -genes in the example dataset (which is a limited number). I then do a -simple hierarchical clustering on the data, extract modules, and step -through each module to calculate enrichment for them, outputting the -results into a separate text file.
-As with the t test comparison above it is important to look at the -number of pathway members included in the comparison (look at the -in_path column). There is no ‘effect size’ problem with Fisher’s exact -since it’s just a set comparison, but it’s important to note that -significant p values can arise from a pathway being -underrepresented in the genelist, which often times is not the -desired result. The foldx column gives a ratio of in versus not in the -genelist, values > 1 being enriched and <1 being depleted.
-
-# for this example we will construct a list of genes from the expression data
-# to emulate what you might be inputting
-genelist <- rownames(pset)[which(SummarizedExperiment::assay(pset,
- "proteomics")[, 1] > 0.5)]
-
-protdata.enrichment.sets.test <- leapR::leapR(
- geneset = ncipid,
- enrichment_method = "enrichment_in_sets",
- eset = pset,
- assay_name = "proteomics",
- targets = genelist
-)
-or <- order(protdata.enrichment.sets.test[, "pvalue"])
-rmarkdown::paged_table(protdata.enrichment.sets.test[or, cols_to_display])
-
-
-
-# in this example we construct some modules from the hierarchical clustering
-# of the data
-protdata_naf <- SummarizedExperiment::assay(pset, "proteomics")
-
-# hierarchical clustering is not too happy with lots of missing values
-# so we'll do a zero fill on this to get the modules
-protdata_naf[which(is.na(protdata_naf))] <- 0
-
-# construct the hierarchical clustering using the 'wardD' method, which
-# seems to give more even sized modules
-protdata_hc <- hclust(dist(protdata_naf), method = "ward.D2")
-
-# arbitrarily we'll chop the clusters into 5 modules
-modules <- cutree(protdata_hc, k = 5)
-
-## sara: created list
-clusters <- lapply(unique(modules), function(x) names(which(modules == x)))
-
-# modules is a named list of values where each value is a module
-# number and the name is the gene name
-
-# To do enrichment for one module (module 1 in this case) do this
-protdata.enrichment.sets.module_1 <- leapR::leapR(
- geneset = ncipid,
- enrichment_method = "enrichment_in_sets",
- eset = pset,
- assay_name = "proteomics",
- targets = names(modules[which(modules == 1)])
-)
-
-# To do enrichment on all modules and return the list of enrichment results
-protdata.enrichment.sets.modules <- do.call(rbind,
- leapR::cluster_enrichment(
- eset = pset,
- assay_name= 'proteomics',
- geneset = ncipid,
- clusters = clusters,
- sigfilter = 0.25))
-## nothing is enriched
-rmarkdown::paged_table(protdata.enrichment.sets.modules[, cols_to_display])
-# Plot the top enriched gene sets from Fisher's exact test
-# Stars indicate significance (None seen here)
-plot_leapr_bar(
- protdata.enrichment.sets.test,
- title = "Fisher's Exact Test: Top Enriched Pathways",
- top_n = 10,
- star_thresholds = c(0.05, 0.01, 0.001),
- wrap = 40
-)
Similar to the popular GSEA, KS tests whether a group of components -(the pathway) is distributed in a statistically significant manner in a -ranked list of components. That is, if all the members of the pathway -are clustered together at the top of the list (highly abundant, e.g.) or -at the bottom of the list (low abundance, e.g.) this will return good p -values. I should note that GSEA uses a more sophisticated approach than -this and their application has a lot of bells and whistles.
-In the example below I’m simply calculating enrichment for one of the -patients in the list (arbitrarily selected). The ranking value is -relative protein abundance in this case, but can be any continuous -measure or derived value. For example, you could calculate the topology -of all proteins in a network and use the topology measure (degree) as -the measure.
-Similar to the other examples be cautious of pathways with good p -values that consider a small number of pathway numbers (in_path column). -The MeanPath column gives a measure that shows how far above or below -the median the mean rank of the pathway is (normalized to -1,1). The -Zscore column is a Zscore calculated on the basis of the mean percentage -rank of the pathway relative to the mean of the entire list divided by -the standard deviation of the pathway rank. The foldx column expresses -the mean percentage rank of the pathway relative to the entire list - -closer to 0 is higher in the list and closer to 1 is closer to the -bottom of the list. Each of these should give consistent results, but -will be somewhat different.
-
-# This is how you calculate enrichment in a ranked list
-# (for example from topology)
-### using enrichment_wrapper function
-protdata.enrichment.order <- leapR::leapR(
- geneset = ncipid, "enrichment_in_order",
- eset = pset,
- method = 'ks',
- assay_name = "proteomics",
- primary_columns = shortlist[1]
-)
-
-
-or <- order(protdata.enrichment.order[, "pvalue"])
-rmarkdown::paged_table(protdata.enrichment.order[or, cols_to_display])
-# Plot the ranked enrichment results
-plot_leapr_bar(
- protdata.enrichment.order,
- title = "Kolmogorov-Smirnov Test: Ranked Enrichment",
- top_n = 12,
- fill_sig = "#2166AC",
- fill_ns = "#B2DFDB",
- wrap = 38
-)
Given that the KS test has not always been the best for -gene set enrichment we also implement the one-sample z test.
-
-# This is how you calculate enrichment in a ranked list
-# (for example from topology)
-### using enrichment_wrapper function
-protdata.enrichment.order <- leapR::leapR(
- geneset = ncipid, "enrichment_in_order",
- eset = pset,
- method = 'ztest',
- assay_name = "proteomics",
- primary_columns = shortlist[1]
-)
-
-
-or <- order(protdata.enrichment.order[, "pvalue"])
-rmarkdown::paged_table(protdata.enrichment.order[or, cols_to_display])
-
-plot_leapr_bar(
- protdata.enrichment.order,
- title = "Z Test: Ranked Enrichment",
- top_n = 12,
- fill_sig = "#2166AC",
- fill_ns = "#B2DFDB",
- wrap = 38
-)
The idea here is to use the correlation of pathway members to each -other versus to non-pathway members as a way to assess functional -enrichment. This idea seems sound- pathways that are varying in a -correlated way across a bunch of conditions (say time points or -patients) may be more active and more important than others. However, -more testing and validation is needed to show that this is the case.
-The ingroup_mean gives the mean correlation of the
-pathway members to each other and outgroup_mean gives the correlation of
-the pathway members to non-pathway members. Background_mean gives the
-mean correlation of all non-pathway members. The pvalue and
-BH_pvalue are for the pathway members to each other versus
-those pathway members to non-pathway components. The
-pvalue_background and BH_pvalue_background are
-for the pathway member correlation relative to non-pathway member
-correlation (which is similar but slightly different than the other
-p-values).
-### using enrichment_wrapper function
-protdata.enrichment.correlation <- leapR::leapR(
- geneset = ncipid,
- enrichment_method = "correlation_enrichment",
- assay_name = "proteomics",
- eset = pset
-)
-
-or <- order(protdata.enrichment.correlation[, "pvalue"])
-rmarkdown::paged_table(head(protdata.enrichment.correlation[or,
- cols_to_display]))
-
-protdata.enrichment.correlation.short <- leapR::leapR(
- geneset = ncipid,
- enrichment_method = "correlation_enrichment",
- assay_name = "proteomics",
- eset = pset[, shortlist]
-)
-or <- order(protdata.enrichment.correlation.short[, "pvalue"])
-rmarkdown::paged_table(head(protdata.enrichment.correlation.short[or,
- cols_to_display]))
-
-protdata.enrichment.correlation.long <- leapR::leapR(
- geneset = ncipid,
- enrichment_method = "correlation_enrichment",
- assay_name = "proteomics",
- eset = pset[, longlist]
-)
-or <- order(protdata.enrichment.correlation.long[, "pvalue"])
-rmarkdown::paged_table(head(protdata.enrichment.correlation.long[or,
- cols_to_display]))
-# Compare correlation patterns across conditions
-plot_leapr_bar(
- protdata.enrichment.correlation.short,
- title = "Correlation Enrichment: Short Survivors",
- top_n = 10,
- fill_sig = "#1B9E77",
- fill_ns = "#D8F0E8",
- wrap = 36
-)
In this example we will use phosphoproteomics data to assess the -enrichment in known kinase substrates (a proxy for kinase activity)
-
-data("kinasesubstrates")
-
-# for an individual tumor calculate the Kinase-Substrate
-# Enrichment (similar to KSEA)
-# This uses the site-specific phosphorylation data to determine
-# which kinases
-# might be active by assessing the enrichment of the
-# phosphorylation of their known substrates
-
-phosphodata.ksea.order <- leapR::leapR(
- geneset = kinasesubstrates,
- enrichment_method = "enrichment_in_order",
- assay_name = "phosphoproteomics",
- eset = phset,
- method = 'ztest',
- primary_columns = "TCGA-13-1484")
-
-or <- order(phosphodata.ksea.order[, "pvalue"])
-rmarkdown::paged_table(phosphodata.ksea.order[or, cols_to_display])
-
-
-# now do the same thing but use a threshold
-phosphodata.sets.order <- leapR::leapR(
- geneset = kinasesubstrates,
- enrichment_method = "enrichment_in_sets",
- eset = phset,
- assay_name = "phosphoproteomics",
- threshold = 0.5,
- primary_columns = "TCGA-13-1484"
-)
-or <- order(phosphodata.sets.order[, "pvalue"])
-rmarkdown::paged_table(phosphodata.sets.order[or, cols_to_display])
-plot <- plot_leapr_bar(
- phosphodata.sets.order,
- title = "Kinase Substrate Enrichment (Phosphoproteomics)",
- top_n = 15,
- star_thresholds = c(0.05, 0.01, 1e-3),
- wrap = 36,
- fill_sig = "#0F766E", # dark teal for significant
- fill_ns = "#99F6E4", # light teal for non-significant
- outline = "grey30",
- axis_text_y_size = 7,
- axis_text_x_size = 8
-)
-
-plot
-
-# You can also modify the plot further using standard ggplot2 arguments
-plot + ggplot2::labs(
- y = expression(-log[10]("adjusted p-value")),
- caption = "* p<0.05, ** p<0.01, *** p<0.001"
-)
Lastly we print out the session info!
-
-sessionInfo()
-#> R version 4.5.2 (2025-10-31)
-#> Platform: aarch64-apple-darwin20
-#> Running under: macOS Tahoe 26.2
-#>
-#> Matrix products: default
-#> BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.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] stringr_1.6.0 tibble_3.3.1 dplyr_1.1.4 ggplot2_4.0.1
-#> [5] rmarkdown_2.30 gplots_3.3.0 leapR_0.99.7 BiocStyle_2.38.0
-#>
-#> loaded via a namespace (and not attached):
-#> [1] SummarizedExperiment_1.40.0 gtable_0.3.6
-#> [3] xfun_0.56 bslib_0.10.0
-#> [5] htmlwidgets_1.6.4 caTools_1.18.3
-#> [7] Biobase_2.70.0 lattice_0.22-7
-#> [9] tzdb_0.5.0 vctrs_0.7.1
-#> [11] tools_4.5.2 bitops_1.0-9
-#> [13] generics_0.1.4 stats4_4.5.2
-#> [15] pkgconfig_2.0.3 Matrix_1.7-4
-#> [17] KernSmooth_2.23-26 RColorBrewer_1.1-3
-#> [19] S7_0.2.1 desc_1.4.3
-#> [21] S4Vectors_0.48.0 lifecycle_1.0.5
-#> [23] compiler_4.5.2 farver_2.1.2
-#> [25] textshaping_1.0.4 Seqinfo_1.0.0
-#> [27] htmltools_0.5.9 sass_0.4.10
-#> [29] yaml_2.3.12 pkgdown_2.2.0
-#> [31] pillar_1.11.1 jquerylib_0.1.4
-#> [33] DelayedArray_0.36.0 cachem_1.1.0
-#> [35] abind_1.4-8 gtools_3.9.5
-#> [37] tidyselect_1.2.1 digest_0.6.39
-#> [39] stringi_1.8.7 bookdown_0.46
-#> [41] labeling_0.4.3 fastmap_1.2.0
-#> [43] grid_4.5.2 cli_3.6.5
-#> [45] SparseArray_1.10.8 magrittr_2.0.4
-#> [47] S4Arrays_1.10.1 readr_2.1.6
-#> [49] withr_3.0.2 scales_1.4.0
-#> [51] XVector_0.50.0 matrixStats_1.5.0
-#> [53] otel_0.2.0 ragg_1.5.0
-#> [55] hms_1.1.4 evaluate_1.0.5
-#> [57] knitr_1.51 GenomicRanges_1.62.1
-#> [59] IRanges_2.44.0 rlang_1.1.7
-#> [61] glue_1.8.0 BiocManager_1.30.27
-#> [63] BiocGenerics_0.56.0 rstudioapi_0.18.0
-#> [65] jsonlite_2.0.0 R6_2.6.1
-#> [67] MatrixGenerics_1.22.0 systemfonts_1.3.1
-#> [69] fs_1.6.6Layered Enrichment Analysis of Pathways in R (leapR) a tool that carries out statistical enrichment analysis on single- or multi-omics data.
-leapR is available through Bioconductor repository here
-if (!require("BiocManager", quietly = TRUE))
- install.packages("BiocManager")
-BiocManager::install(version = "3.21")
-BiocManager::install('BiocStyle')
-BiocManager::install('leapR')
-Once you have successfully installed the package you can load the vignette to read examples using the vignette('leapR') command.
The primary function of the leapR package is the leapR function itself. This function serves a wrapper to run different styles of enrichment functions on the data. The package contains other functions to support pathway information and multi-omics datasets.
Here is a list of enrichment arguments that can be called with the leapR command.
| Argument | -Description | -
|---|---|
enrichment_in_sets |
-Calculates enrichment in pathway membership in a list (e.g. highly differential proteins) relative to background using Fisher’s exact test. | -
enrichment_in_order |
-Calculates enrichment of pathways based on a ranked list using the Kolmogorov-Smirnov test | -
enrichment_comparison |
-Compares the distribution of abundances between two sets of conditions for each pathway using a t test | -
enrichment_in_pathways |
-Compares the distribution of abundances in a pathway with the background distribution of abundances using a t test | -
correlation_enrichment |
-Calculates the enrichment of a pathway based on correlation between pathway members across conditions versus correlation between members not in the pathway | -
enrichment_in_relationships |
-Calculates the enrichment of a pathway in specified interactions relative to non-pathway members | -
calcTTest.Rdcalculates a t-test for two distributions of data on a per-gene basis -append results to ExpressionSet with two extra columns: `pvalue` and -`difference` for each feature
-
- library(leapR)
- url <- "https://api.figshare.com/v2/file/download/56536214"
- tdata <- download.file(url,method='libcurl',destfile='transData.rda')
- load('transData.rda')
- p <- file.remove("transData.rda")
-
- # read in the pathways
- data("ncipid")
-
- # read in the patient groups
- data("shortlist")
- data("longlist")
- calcTTest(tset, 'transcriptomics', shortlist, longlist)
-#> class: SummarizedExperiment
-#> dim: 1999 174
-#> metadata(0):
-#> assays(1): transcriptomics
-#> rownames(1999): NOC2L ISG15 ... ARL6 MINA
-#> rowData names(2): pvalue difference
-#> colnames(174): TCGA-13-1484 TCGA-13-1495 ... TCGA-61-1995 TCGA-61-2008
-#> colData names(0):
-cluster_enrichment.RdCluster enrichment Run enrichment (Fisher's exact) on clusters (lists of -identifier groups)
-is an SummarizedExperiment containing data that is clustered
is the name of the assay
is a GeneSet object for pathway annotation
is a list of clusters (gene lists) to calculate enrichment -on, generally the result of the `cutree` function
minimum significance threshold default is .05
This function will calculate enrichment (Fisher's exact test for - membership overlap) on
-a series of lists of genes, such as from a set of clusters. The -results are returned as
-a list of results matrices in the order of the input clusters.
- library(leapR)
-
- # read in the example transcriptomic data
- url <- "https://api.figshare.com/v2/file/download/56536214"
- tdata <- download.file(url,method='libcurl',destfile='transData.rda')
- load('transData.rda')
- p <- file.remove("transData.rda")
-
- # read in the pathways
- data("ncipid")
-
- # for the example we will limit the number of transcripts considered
- #- arbitrarily in this case
- transdata <- SummarizedExperiment::assay(tset,'transcriptomics')
- transdata[which(is.na(transdata),arr.ind=TRUE)]<-0.0
- # perform heirarchical clustering on the data
- transdata.hc <- hclust(dist(transdata), method="ward.D2")
-
- transdata.hc.clusters <- cutree(transdata.hc, k=5)
- clust.list <- lapply(seq_len(5), function(x) {
- return(names(which(transdata.hc.clusters==x)))})
- #calculates enrichment for each of the clusters individually a
- #and returns a list of enrichment results
- transdata.hc.enrichment <- leapR::cluster_enrichment(eset=tset,
- assay_name='transcriptomics',
- geneset=ncipid,
- clusters=clust.list)
-
-
-
-combine_omics.Rdcombine_omics -Combine two or more omics matrices into one multi-omics matrix with -'tagged' ids.
-This combines matrices of different omics types together and -adds prefix tags to the ids.
- library(leapR)
- url <- 'https://api.figshare.com/v2/file/download/56536217'
-
- pdata <- download.file(url,method='libcurl',destfile='protData.rda')
- load('protData.rda')
- p <- file.remove("protData.rda")
-
- url <- "https://api.figshare.com/v2/file/download/56536214"
- tdata <- download.file(url,method='libcurl',destfile='transData.rda')
- load('transData.rda')
- p <- file.remove("transData.rda")
-
- url <- 'https://api.figshare.com/v2/file/download/56536211'
- phdata<-download.file(url,method='libcurl',destfile = 'phosData.rda')
- #phosphodata<-read.csv("phdata",check.names=FALSE,row.names=1)
- load('phosData.rda')
- p <- file.remove('phosData.rda')# read in the example protein data
-
-
- # merge the three datasets by rows and add prefix tags for
- # different omics types
- multi_omics <- combine_omics(list(pset, tset, phset),
- list(NA,NA,'hgnc_id'))
-
-
-correlation_comparison_enrichment.Rd# internal function to calculate enrichment in differences in correlation -# between two groups -# access through the leapr wrapper
-correlation_comparison_enrichment(
- geneset,
- eset,
- assay_name,
- set1,
- set2,
- mapping_column = NA
-)enrichment_in_abundance.RdEnrichment in abundance calculates enrichment in pathways by the difference -in abundance of the pathway members.
-enrichment_in_abundance(
- geneset,
- eset,
- assay_name,
- mapping_column = NULL,
- abundance_column = NULL,
- fdr = 0,
- matchset = NULL,
- sample_comparison = NULL,
- min_p_threshold = NULL,
- sample_n = NULL,
- silence_try_errors = TRUE
-)Gene set to calculate enrichment
Molecular abundance data in `SummarizedExperiment` format
Name of assay to compare
Column to use to map identifiers
Columns to use to quantify abundance
number of times to sample for FDR value
Name of a set to use for enrichment
list of samples to use as comparison. if missing -background (eset) is used
Only include p-values lower than this
size of sample to use
set to true to silence try errors
enrichment_in_groups.RdCalculate the enrichment in pathways using Fisher's exact or -Kolmogorov-Smirnov test, using either the abundance column to identify -feature or the targets list. access through leapr wrapper
-enrichment_in_groups(
- geneset,
- targets = c(),
- background = NULL,
- assay_name = NULL,
- method = "fishers",
- minsize = 5,
- mapping_column = NULL,
- abundance_column = NULL,
- randomize = FALSE,
- silence_try_errors = TRUE
-)geneset to use for enrichment
targets to use for enrichment
`SummarizedExperiment` describing background to use
is the name of the assay to use from the background
method to use for statistical test, options are -'fishers', 'ks', 'ztest', or 'chisq'. Remember that KS test assumes normality, so it would be good -to log your data before calling. NOTE: if you do not call `suppressWarnings` then -the KS test will warn you about ties.
minimum size of set
column name of mapping identifiers
columns mapping abundance, either in the `assay` -matrix or `rowData`
true/false whether to randomize
true/false to silence errors
enrichment_in_relationships.Rdenrichment_in_relationships function description is a general way to -determine if a pathway -is enriched in relationships (interactions, correlation) between its members -# access through leapr wrapper
-enrichment_in_relationships(
- geneset,
- relationships,
- idmap = NA,
- silence_try_errors = TRUE
-)calcTTest()
-
- cluster_enrichment()
-
- combine_omics()
-
- correlation_comparison_enrichment()
-
- correlation_enrichment()
-
- enrichment_in_abundance()
-
- enrichment_in_groups()
-
- enrichment_in_relationships()
-
- get_pathway_information()
-
- kinasesubstrates
-
- krbpaths
-
- leapR()
-
- longlist
-
- ncipid
-
- plot_leapr_bar()
-
- read_gene_sets()
-
- shortlist
-
- krbpaths.RdKEGG, Reactome, BioCarta Pathways
-leapR-package.RdleapR is a package that identifies pathways that are enriched across diverse 'omics experiments. It leverages any tabular expression data (proteomics, transcriptomics) using the `SummarizedExperiment` object. It works with any pathway in the .gct file format.
-Maintainer: Sara Gosline sara.gosline@pnnl.gov (ORCID)
-Authors:
Jason McDermott jason.mcdermott@pnnl.gov
Other contributors:
Vincent Danna vincent.danna@pnnl.gov [contributor]
National Institutes of Health [funder]
leapR.RdleapR is a wrapper function that consolidates multiple enrichment methods.
-is a list of four vectors, gene names, gene descriptions, gene -sizes and a matrix of genes. It represents .gmt format pathway files.
is a character string specifying the method of -enrichment to be performed, one of: "enrichment_comparison", -"enrichment_in_order", "enrichment_in_sets", "enrichment_in_pathway", -"correlation_enrichment".
is a `SummarizedExperiment` object containing expression data, -with features as rows and n sample/conditions as columns.
is the assay to be analyzed within the `eset`. -Recommended to describe the data type (e.g. transcriptomics, proteomics) -so that it can be integrated in `combine_omics`
further arguments
Further arguments and enrichment method optional argument information:
| id_column | Is a character string, present in the rowData slot,
-that is used to specify a column for identifiers to map to enrichment
-libraries.
-If missing, the rownames of the SummarizedExperiment assay will be used. |
| primary_columns | Is a character vector composed of column names from
-eset (either in the `assay` or in the `rowData`),
-that specifies a set of primary columns to calculate enrichment on.
-The meaning of this varies according to the enrichment method used - see
-the descriptions for each method below.
-This is an optional argument used with 'enrichment_in_order',
-'enrichment_in_sets', and 'enrichment_comparison' methods. |
| secondary_columns | |
| Is a character vector of column names for comparison, -pulled from the `assay` of the SummarizedExperiment. This is an -optional argument used with 'enrichment_comparison' methods. | |
| threshold | Is a numeric value, an optional argument used with -'enrichment_in sets' method which filters out abundance values or p-values -(depending on what `primary_columns` is used) -either above or below it. |
| greaterthan | |
Is a logical value that defaults to TRUE, it's used with
-'enrichment_in_sets' method.
-When set to TRUE, genes with `primary_columns` value above the
-threshold argument are kept.
-When set to FALSE genes with `primary_columns` value below the
-threshold argument are kept.
-This is an optional argument used with 'enrichment_in_sets' method. | |
| minsize | Is a numeric value, an optional argument used with -'enrichment_in_sets' and 'enrichment_in_order". |
| fdr | |
| A numerical value which specifies how many times to randomly -sample genes to calculate an empirical false discovery rate, is an optional -argument used with 'enrichment_comparison' method. | |
| min_p_threshold | Is a numeric value, a lower p-value threshold and is an -optional argument used with 'enrichment_comparison' method. |
| sample_n | |
| Is a way to subsample the number of components considered for -each calculation randomly. This is an optional argument used with -'enrichment_comparison' method. |
Enrichment Methods:
-
-enrichment_comparison
-
-Compares the distribution of abundances between two sets of
-conditions for each pathway using a t test. For each pathway in
-geneset uses a t test to compare the distribution of abundance
-values/numbers in eset primary_columns with those in
-eset secondary_columns. Lower p-values for pathways indicate
-that the expression of the pathway is significantly different between the
-set of conditions in primary_columns and the set of conditions in
-secondary_columns.
-Optionally, users can specify fdr which will calculate an empirical
-p-value by randomizing abundances fdr number of times. If the
-min_p_threshold is specified the method will only return pathways
-with an adjusted p-value lower than the specified threshold. If
-sample_n is specified the method will subsample the
-pathway members to the specified number of components.
-
-enrichment_in_order
-
-Calculates enrichment of pathways based on a ranked list using the
-Kolmogorov-Smirnov test. For each pathway in geneset uses a
-Kolmogorov-Smirnov test for rank order to test if the distribution of ranked
-abundance values in the eset primary_columns is significant
-relative to a random distribution. Note that currently
-primary_columns only accepts a single column for this method.
-
-enrichment_in_sets
-
-Calculates enrichment in pathway membership in a list (e.g. highly
-differential proteins) relative to background using Fisher's exact test. For
-each pathway in geneset uses a Fisher's exact test over- or under-
-representation of a list of components specified. If targets are
-specified this must be a vector of identifiers to serve as the target list
-for comparison. If eset and primary_columns are specified then
-threshold specifies a threshold value for determining the target list
-of components to test. Specifying greaterthan to be False
-will result in components with values lower than the specified
-threshold. If eset is a data frame or matrix, the background
-used for calculation will be taken as the rownames of eset
-
-enrichment_in_pathway
-
-Compares the distribution of abundances in a pathway with the background
-distribution of abundances using a t test. For each pathway in
-geneset calculates the significance of the difference between the
-abundances from pathway members versus abundance of non-pathway members in
-the set of conditions specified by primary_columns. Optionally, users
-can specify fdr which will calculate an empirical p-value by
-randomizing abundances fdr number of times. If the
-min_p_threshold is specified the method will only return pathways
-with an adjusted p-value lower than the specified threshold. If
-sample_n is specified the method will subsample the
-pathway members to the specified number of components.
-
-correlation_enrichment
-
-Calculates the enrichment of a pathway based on correlation between pathway
-members across conditions versus correlation between members not in the
-pathway. For each pathway in geneset calculates the pairwise
-correlation between all pathway members and non-pathway members
-across the specified primary_columns conditions in eset. Note
-that for large matrices this can take a long time. A p-value is calculated
-based on comparing the correlation within the members of a pathway with the
-correlation values between members of the pathway and non-members of the
-pathway.
-
library(leapR)
-
- # read in the example abundance data
- # read in the example transcriptomic data
- tdata <- download.file("https://api.figshare.com/v2/file/download/56536214",
- method='libcurl',destfile='transData.rda')
- load('transData.rda')
- p <- file.remove("transData.rda")
-
- # read in the pathways
- data("ncipid")
-
- # read in the patient groups
- data("shortlist")
- data("longlist")
-
- # use enrichment_comparison to calculate enrichment in one set of
- # conditions (shortlist) and another (longlist)
- short_v_long = leapR(geneset=ncipid, assay_name='transcriptomics',
- enrichment_method='enrichment_comparison',
- eset=tset, primary_columns=shortlist,
- secondary_columns=longlist)
-
- # use enrichment_in_sets to calculate the most enriched pathways
- # from the highest abundance proteins
- # from one condition
- onept_sets = leapR(geneset=ncipid, assay_name='transcriptomics',
- enrichment_method='enrichment_in_sets',
- eset=tset, primary_columns="TCGA-13-1484", threshold=1.5)
-
- # use enrichment_in_order to calculate the most enriched pathways from the
- # same condition
- # Note: that this uses the entire set of abundance values and their order -
- # whereas the previous example uses a hard threshold to get a short list of
- # most abundant proteins and calculates enrichment based on set overlap.
- # The results are likely to be similar - but with some notable differences.
- onept_order = leapR(geneset=ncipid, assay_name='transcriptomics',
- enrichment_method='enrichment_in_order',
- eset=tset, primary_columns="TCGA-13-1484")
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-#> Warning: p-value will be approximate in the presence of ties
-
- # use enrichment_in_pathway to calculate the most enriched pathways in a
- # set of conditions based on abundance in the pathway members versus
- # abundance in non-pathway members
- short_pathways = leapR(geneset=ncipid, assay_name='transcriptomics',
- enrichment_method='enrichment_in_pathway',
- eset=tset, primary_columns=shortlist)
-
- # use correlation_enrichment to calculate the most enriched pathways in
- # correlation across the shortlist conditions
- short_correlation_pathways = leapR(geneset=ncipid,
- assay_name='transcriptomics',
- enrichment_method='correlation_enrichment',
- eset=tset, primary_columns=shortlist)
-
-
-ncipid.RdA list of pathways and the genes that comprise these pathways
-plot_leapr_bar.RdThis plotting helper expects leapR generated results to plot. -It will use BH_pvalue if present, otherwise pvalue.
-plot_leapr_bar(
- res_df,
- title = NULL,
- top_n = 15,
- star_thresholds = c(0.05, 0.01, 0.001),
- wrap = 42,
- max_stars = 5L,
- fill_sig = "#2C7BB6",
- fill_ns = "#BFD7FF",
- outline = NA,
- axis_text_y_size = 8,
- axis_text_x_size = 9
-)A leapR df containing BH_pvalue (or pvalue) and a pathway/term label column.
Plot title.
Number of top pathways/genes to display.
list of numeric significance thresholds for star annotations.
Wrap width for pathway labels (helps formatting).
Maximum number of stars to draw per bar (default 5).
Fill color for significant bars
Fill color for non-significant bars.
Bar border color.
Font size for y-axis (category) labels.
Font size for x-axis (numeric) labels.
read_gene_sets.Rdread_gene_sets is a function to import external pathway -database files in .gmt format
-read_gene_sets(
- gsfile,
- gene.labels = NA,
- gs.size.threshold.min = 5,
- gs.size.threshold.max = 15000
-)gfile <- system.file('extdata','h.all.v2024.1.Hs.symbols.gmt',
- package='leapR')
-glist <- read_gene_sets(gfile)
-
-