diff --git a/.gitignore b/.gitignore index cdf82c5..53d3b2c 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,6 @@ .Rbuildignore SPAROscore.Rproj inst/doc +/vignettes/*.html +/vignettes/*.R + diff --git a/DESCRIPTION b/DESCRIPTION index 8682a84..323961c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,24 +1,26 @@ -Package: SPAROscore Type: Package -Title: A package to compute gene signature scores from transcriptomics data +Package: SPAROscore +Title: A package to compute gene signature scores from transcriptomics + data Version: 0.99.0 -Authors@R: c( - person( - "Venkatesh", "Kamaraj", - email = "venkatesh.kamaraj@adelaide.edu.au", - role = c("aut", "cre"), - comment = c(ORCID = "0009-0008-7309-9810") - ) - ) -Description: SPAROscore is a gene signature scoring method designed to be - robust across diverse gene expression datasets. SPAROscore adapts to - varying levels of sparsity, allowing signature scores to be efficiently - computed across bulk, single-cell, and spatial transcriptomic datasets. -biocViews: GeneExpression, SingleCell, Spatial, Transcriptomics, - GeneSetEnrichment, Pathways, Software +Authors@R: + person("Venkatesh", "Kamaraj", , "venkatesh.kamaraj@adelaide.edu.au", role = c("aut", "cre"), + comment = c(ORCID = "0009-0008-7309-9810")) +Description: SPAROscore is a gene signature scoring method designed to be + robust across diverse gene expression datasets. SPAROscore adapts to + varying levels of sparsity, allowing signature scores to be + efficiently computed across bulk, single-cell, and spatial + transcriptomic datasets. The resulting scores quantify the relative + expression of a gene signature compared with the background expression + of each sample/cell/domain, making the scores straightforward to + interpret biologically. Internally, SPAROscore ranks genes within + each column and computes the Spearman footrule distance between the + observed ranks of the signature genes and a background gene expression + rank estimated from the geometric mean expression of the sample, cell, + or spatial domain. License: GPL-3 -Encoding: UTF-8 -RoxygenNote: 7.3.3 +URL: https://github.com/MangiolaLaboratory/SPAROscore +BugReports: https://github.com/MangiolaLaboratory/SPAROscore/issues Imports: DelayedArray, DelayedMatrixStats, @@ -37,18 +39,20 @@ Suggests: GSEABase, knitr, msigdbr, + muscData, rmarkdown, scater, Seurat, - SeuratData, SingleCellExperiment, SpatialExperiment, spatialLIBD, SummarizedExperiment, testthat (>= 3.0.0), tidyr -VignetteBuilder: knitr +VignetteBuilder: + knitr +biocViews: GeneExpression, SingleCell, Spatial, Transcriptomics, + GeneSetEnrichment, Pathways, Software Config/testthat/edition: 3 -Remotes: satijalab/seurat-data -URL: https://github.com/MangiolaLaboratory/SPAROscore -BugReports: https://github.com/MangiolaLaboratory/SPAROscore/issues +Encoding: UTF-8 +RoxygenNote: 7.3.3 diff --git a/R/get_scores_generic.R b/R/get_scores_generic.R index f4abe61..a9ff6b2 100644 --- a/R/get_scores_generic.R +++ b/R/get_scores_generic.R @@ -59,7 +59,7 @@ #' Defaults to NULL. #' #' When down_signatures is not NULL, -#' Final Score = Score(signatures) - Score(down_signatures) +#' Final Score <- Score(signatures) - Score(down_signatures) #' #' #' @param prefix Character string to be appended before the headers of the diff --git a/R/helper_functions.R b/R/helper_functions.R index 328c45a..839e86f 100644 --- a/R/helper_functions.R +++ b/R/helper_functions.R @@ -229,7 +229,7 @@ append_to_matrix_like_object <- function(matrix_like_object, numeric_vector){ #' #' #' @examples -#' \dontrun{ +#' \donttest{ #' counts <- matrix( #' c(0, 5, 10, #' 2, 0, 20, diff --git a/R/sparoscore_generic.R b/R/sparoscore_generic.R index f5fc845..92e23c6 100644 --- a/R/sparoscore_generic.R +++ b/R/sparoscore_generic.R @@ -111,7 +111,7 @@ #' Defaults to NULL. #' #' When down_signatures is not NULL, -#' Final Score = Score(signatures) - Score(down_signatures) +#' Final Score <- Score(signatures) - Score(down_signatures) #' #' #' @param prefix Character string to be appended before the headers of the @@ -218,7 +218,7 @@ #' # Seurat object #' # ------------------------------------------------------------------ #' -#' \dontrun{ +#' \donttest{ #' seurat_object <- sparoscore( #' data = seurat_object, #' signatures = c("CCR7", "IL7R", "LTB") @@ -257,7 +257,7 @@ #' # SingleCellExperiment / SummarizedExperiment #' # ------------------------------------------------------------------ #' -#' \dontrun{ +#' \donttest{ #' sce <- sparoscore( #' data = sce, #' signatures = c("CCR7", "IL7R", "LTB") diff --git a/vignettes/.gitignore b/vignettes/.gitignore deleted file mode 100644 index 097b241..0000000 --- a/vignettes/.gitignore +++ /dev/null @@ -1,2 +0,0 @@ -*.html -*.R diff --git a/vignettes/sparoscore_data_structures.Rmd b/vignettes/sparoscore_data_structures.Rmd index 2ddb5df..4e656cb 100644 --- a/vignettes/sparoscore_data_structures.Rmd +++ b/vignettes/sparoscore_data_structures.Rmd @@ -10,7 +10,7 @@ vignette: > %\VignetteEncoding{UTF-8} --- -```{r, include = FALSE} +```{r global-options, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", @@ -97,7 +97,7 @@ as well as the hallmark inflammatory response genes taken from the Molecular Signatures Database (MSigDB). ### Load the gene expression dataset -```{r} +```{r load-data} # load the airway dataset data(airway) @@ -109,7 +109,7 @@ metadata <- as.data.frame(colData(airway)) ``` ### Load the gene signatures -```{r} +```{r load-signature} # set the glucocorticoid response markers from the original study as a signature response_genes <- c("ENSG00000112936", "ENSG00000170606", "ENSG00000120129", "ENSG00000152795", "ENSG00000211445", "ENSG00000163884", @@ -131,7 +131,7 @@ inflammation_genes <- msigdbr(species = "Homo sapiens", collection = "H") %>% # Gene expression as matrix-like data structures ## Dense matrix -```{r} +```{r matrix} # load gene expression data as dense matrix counts_matrix <- as.matrix(counts_data) @@ -147,7 +147,7 @@ Note that one signature gene is absent from the input expression data, and SPAROscore reports this as a warning during scoring. ## SparseMatrix -```{r} +```{r smatrix} # load gene expression data as sparseMatrix counts_sparse <- as(counts_matrix, "sparseMatrix") @@ -164,7 +164,7 @@ them to dense matrices, reducing memory usage and improving computational efficiency. ## DelayedArrays -```{r} +```{r dmatrix} # load gene expression data as delayedArray counts_delayed <- DelayedArray(counts_matrix) @@ -180,7 +180,7 @@ Note that the gene ranking is performed natively on the DelayedArray object, enabling improved computational efficiency. ## Data frame -```{r} +```{r dataframe} # load gene expression data as delayedArray counts_df <- as.data.frame(counts_matrix) @@ -197,7 +197,7 @@ because matrix operations are more efficient for this step. ### Verify if the output scores are the same with differing input formats -```{r} +```{r check-results} identical(scores_matrix, scores_sparse) identical(scores_matrix, scores_delayed) identical(scores_matrix, scores_df) @@ -206,7 +206,7 @@ identical(scores_matrix, scores_df) # Gene signature representations ## Single signature as a character vector -```{r} +```{r character} # compute scores from a character vector signature scores_for_vector <- sparoscore(data = counts_matrix, signatures = response_genes) @@ -217,20 +217,21 @@ scores_for_vector ## Single signature as a GeneSet object -```{r} +```{r geneset} # get a GeneSet object response_geneset <- GeneSet(response_genes) # compute scores from a GeneSet object scores_for_geneset <- sparoscore(data = counts_sparse, - signatures = response_geneset, prefix = "hello") + signatures = response_geneset, + prefix = "hello") # view scores scores_for_geneset ``` ## Multiple signatures as a named list -```{r} +```{r list} # create a named list of vectors to score for both signatures simultaneously gene_signatures <- list("glucocorticoid_response" = response_genes, "inflammation_response" = inflammation_genes) @@ -244,7 +245,7 @@ scores_for_list ``` ## Multiple signatures as a GeneSetCollection object -```{r} +```{r genesetcollection} # create a GeneSetCollection object to score for both signatures simultaneously gene_set_collection <- GeneSetCollection( mapply(function(g, n) GeneSet(unique(g), setName = n), @@ -268,7 +269,7 @@ named `ranks`. In this example, let us work with a SingleCellExperiment object -```{r} +```{r sce} # convert airways to a SingleCellExperiment object airway_sce <- as(airway, "SingleCellExperiment") @@ -309,7 +310,7 @@ For Seurat objects, after scoring, SPAROscore returns a Seurat object with the scores and the rank caps appended to the metadata. The ranks are stored in a layer named `ranks` in the same assay as the counts. -```{r} +```{r seurat} # convert airway to a seurat object airway_seurat <- CreateSeuratObject( counts = assay(airway), @@ -364,6 +365,6 @@ For examples of in-depth transcriptomics analyses using SPAROscore, refer here: For examples of additional options and advanced usage, see: [Advanced options](sparoscore_advanced_options.html) -```{r} +```{r session-info} sessionInfo() ``` diff --git a/vignettes/sparoscore_getting_started.Rmd b/vignettes/sparoscore_getting_started.Rmd index ffb3c3b..496b502 100644 --- a/vignettes/sparoscore_getting_started.Rmd +++ b/vignettes/sparoscore_getting_started.Rmd @@ -10,7 +10,7 @@ vignette: > %\VignetteEncoding{UTF-8} --- -```{r, include = FALSE} +```{r global-options, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", @@ -61,7 +61,7 @@ reported in the [original study](https://doi.org/10.1371/journal.pone.0099625). The expression matrix contains genes in rows and samples in columns. The accompanying metadata indicates whether each sample received dexamethasone treatment (trt) or served as an untreated control (untrt). -```{r} +```{r load-dataset} # load the airway dataset data(airway) @@ -87,7 +87,7 @@ head(signature_genes %in% rownames(counts_data)) ``` ## Compute signature scores using SPAROscore -```{r} +```{r compute-scores} scores <- sparoscore(data = counts_data, signatures = signature_genes) # view the scores @@ -99,7 +99,7 @@ scores Since dexamethasone activates glucocorticoid signalling, treated samples are expected to exhibit higher glucocorticoid response scores than untreated controls. -```{r} +```{r visualise-results} # add the signature scores to the metadata metadata$SPAROscore <- scores[, 1] @@ -159,6 +159,6 @@ For examples of additional options and more advanced usage, see: [Advanced options](sparoscore_advanced_options.html) -```{r} +```{r session-info} sessionInfo() ``` diff --git a/vignettes/sparoscore_transcriptomics_technologies.Rmd b/vignettes/sparoscore_transcriptomics_technologies.Rmd index d20017c..88ed0a0 100644 --- a/vignettes/sparoscore_transcriptomics_technologies.Rmd +++ b/vignettes/sparoscore_transcriptomics_technologies.Rmd @@ -10,7 +10,7 @@ vignette: > %\VignetteEncoding{UTF-8} --- -```{r, include = FALSE} +```{r global-options, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", @@ -29,7 +29,7 @@ library(ggplot2) library(dplyr) library(edgeR) library(Seurat) -library(SeuratData) +library(muscData) library(msigdbr) library(scater) library(spatialLIBD) @@ -103,7 +103,7 @@ The expression matrix contains genes in rows and samples in columns. The accompanying metadata indicates whether each sample received dexamethasone treatment (trt) or served as an untreated control (untrt). -```{r} +```{r bulk-1} # load the airway dataset data(airway) @@ -132,7 +132,7 @@ The number of signature missing from the dataset will also be reported by SPAROscore as a warning. ## Compute SPAROscore -```{r} +```{r bulk-2} bulk_scores <- sparoscore( data = counts_data, signatures = gc_response_genes) bulk_scores @@ -143,7 +143,7 @@ signature genes have higher expression relative to the background expression of the sample, consistent with stronger pathway activation. ## Combine scores with the sample metadata -```{r} +```{r bulk-3} metadata$SPAROscore <- bulk_scores[,1] metadata[, c("Run", "dex", "SPAROscore")] @@ -157,7 +157,7 @@ Since dexamethasone activates glucocorticoid signalling, treated samples are expected to exhibit higher glucocorticoid response scores than untreated controls. -```{r} +```{r bulk-4} metadata %>% ggplot(aes(x = dex, y = SPAROscore, fill = dex)) + geom_boxplot(width = 0.6, outlier.shape = NA) + @@ -177,7 +177,7 @@ successfully captures the biological response induced by the treatment. The observed difference can be formally evaluated using the non-parametric Wilcoxon ranksum test -```{r} +```{r bulk-5} wilcox.test(SPAROscore ~ dex, data = metadata) ``` @@ -191,7 +191,7 @@ sample-to-sample variation in transcriptomic datasets. By colouring samples according to their SPAROscore rather than their treatment group, we can determine whether pathway activity explains the dominant source of variation. -```{r, warning=FALSE} +```{r bulk-6, warning=FALSE} # Create the DGEList object bulk_dge <- DGEList(counts = counts_data) @@ -243,7 +243,7 @@ To illustrate this, we compare the SPAROscore with the expression of most differentially expressed glucocorticoid response gene (C7) from the [original study](https://doi.org/10.1371/journal.pone.0099625). -```{r} +```{r bulk-7} # assign the top gene from the original study top_gene <- "ENSG00000112936" @@ -295,7 +295,8 @@ cell types and experimental conditions to be explored. In this example, the IFN-β stimulation dataset from [Kang et al.](https://doi.org/10.1038/nbt.4042), which is available through the -SeuratData package will be used. The dataset consists of ~14,000 +muscData package will be used. We will work with a seurat object of the same to +showcase the versatality of SPAROscore. The dataset consists of ~29,000 peripheral blood mononuclear cells (PBMCs) collected from lupus patients before and after stimulation with interferon-β (IFN-β). @@ -305,17 +306,22 @@ Although the cells were stimulated with IFN-β, the Hallmark collection contains a single canonical type I interferon response signature, which captures the conserved transcriptional response induced by both IFN-α and IFN-β signalling. -## Load the dataset from SeuratData +## Load the dataset from muscData -```{r, warning=FALSE} -# install the dataset for first-time use -InstallData("ifnb") +```{r sc-1, warning=FALSE} +# load the ifnb data from the muscData library +ifnb_sce_data <- Kang18_8vs8() -# load the ifnb dataset -ifnb <- LoadData("ifnb") +# working with a seurat object +ifnb <- CreateSeuratObject( + counts = counts(ifnb_sce_data), + meta.data = as.data.frame(colData(ifnb_sce_data)), + assay = "RNA" +) # check the layer in the RNA assay Layers(ifnb[["RNA"]]) +#> [1] "counts" "data" # view the metadata of the ifnb seurat object head(ifnb[[]]) @@ -328,7 +334,7 @@ are also given under the `seurat_annotations` metadata column. ## Load the gene signature from MSigDB -```{r} +```{r sc-2} ifn_signature <- msigdbr( species = "Homo sapiens", collection = "H" ) %>% filter(gs_name == "HALLMARK_INTERFERON_ALPHA_RESPONSE") %>% pull(gene_symbol) %>% @@ -342,7 +348,7 @@ head(ifn_signature %in% rownames(ifnb)) ## Compute SPAROscore The interferon-beta response for every cell can be quantified using SPAROscore -```{r} +```{r sc-3} ifnb <- sparoscore(data = ifnb, assay = "RNA", layer = "counts", @@ -370,7 +376,7 @@ Since interferon signalling is expected to increase following stimulation, the distribution of SPAROscore between the two experimental conditions should reflect that. -```{r} +```{r sc-4} VlnPlot(ifnb, features = "IFN_SPAROscore", group.by = "stim", @@ -388,7 +394,7 @@ scores than unstimulated controls. The difference between the two experimental groups can also be assessed statistically using the non-parametric Wilcoxon ranksum test. -```{r} +```{r sc-5} wilcox.test(IFN_SPAROscore ~ stim, ifnb[[]]) ``` @@ -401,10 +407,10 @@ to quantify that biological phenomenon. Because pathway scores are calculated for every individual cell, they can be compared across annotated immune cell populations. -```{r, warning=FALSE} +```{r sc-6, warning=FALSE} VlnPlot(ifnb, features = "IFN_SPAROscore", - group.by = "seurat_annotations", + group.by = "cell", pt.size = 0, flip = TRUE ) @@ -420,7 +426,7 @@ transcriptional response following interferon stimulation. One of the most useful ways to explore pathway activity is to visualise the scores on the UMAP embeddings. -```{r, warning=FALSE, message=FALSE} +```{r sc-7, warning=FALSE, message=FALSE} # get umap embedding from the gene expression data # normalise the gene expression @@ -461,7 +467,7 @@ computationally intensive ranking step. For example, an inflammatory response signature could be scored as follows. -```{r} +```{r sc-8} # load the inflammatory response hallmark genes from msigdb inflammation_signature <- msigdbr(species = "Homo sapiens",collection = "H") %>% filter(gs_name == "HALLMARK_INFLAMMATORY_RESPONSE") %>% @@ -500,9 +506,9 @@ quantify interferon pathway activity across each cell population. ### Aggregate cells into pseudo-bulk profiles -```{r} +```{r pb-1} pseudobulk_ifn <- AggregateExpression(ifnb, - group.by = c("seurat_annotations", "stim"), + group.by = c("cell", "stim"), assays = "RNA", return.seurat = TRUE, slot = "counts") @@ -514,7 +520,7 @@ cell type and treatment (for example, "B_CTRL" or "B_STIM"), while rows represent genes. ### Compute scores at pseudo-bulk level -```{r} +```{r pb-2} pseudobulk_ifn <- sparoscore(data = pseudobulk_ifn, signatures = ifn_signature, assay = "RNA", @@ -522,9 +528,9 @@ pseudobulk_ifn <- sparoscore(data = pseudobulk_ifn, ``` ### Visualise pseudo-bulk level scores -```{r} +```{r pb-3} pseudobulk_ifn[[]] %>% - ggplot(aes(y = seurat_annotations, x = SPAROscore, fill = stim)) + + ggplot(aes(y = cell, x = SPAROscore, fill = stim)) + geom_col(position = position_dodge(width = 0.8)) + theme_classic() + scale_fill_manual(values = c("#2166AC", "#B2182B"))+ @@ -590,7 +596,7 @@ lower in the cortex. For computational efficiency, we analyse one representative tissue section. -```{r} +```{r sp-1} # load the DLPFC Visium dataset dlpfc_spe <- fetch_data(type = "spe") @@ -602,7 +608,7 @@ The object is a SpatialExperiment, containing spot-level gene expression counts, spatial coordinates, tissue images, and manually annotated cortical layers. ## Visualise the tissue morphology -```{r} +```{r sp-2} # for simplicity of the vignette, work with one sample dlpfc_spe <- dlpfc_spe[ , dlpfc_spe$sample_id == "151675"] @@ -617,7 +623,7 @@ vis_clus(spe = dlpfc_spe, For demonstration purposes, we use a collection of well-established oligodendrocyte marker genes. -```{r} +```{r sp-3} oligo_signature <- c( "ENSG00000167748", "ENSG00000166923", "ENSG00000174607", "ENSG00000135678", "ENSG00000089123", "ENSG00000183018", "ENSG00000129538", "ENSG00000136937", @@ -636,7 +642,7 @@ head(oligo_signature %in% rownames(dlpfc_spe)) Since SpatialExperiment objects are supported directly by SPAROscore, pathway scores can be computed without extracting the expression matrix. -```{r} +```{r sp-4} dlpfc_spe <- sparoscore(data = dlpfc_spe, assay = "counts", signatures = oligo_signature, @@ -658,7 +664,7 @@ visualise biological processes within their anatomical context. Using the spot coordinates stored in the SpatialExperiment object, the SPAROscore can be projected directly onto the tissue section. -```{r} +```{r sp-5} vis_gene(spe = dlpfc_spe, geneid = "Oligo_SPAROscore", spatial = TRUE, @@ -679,7 +685,7 @@ of oligodendrocyte gene expression. The manual annotations for each cortical layer enables quantitative comparison of oligodendrocyte pathway activity. -```{r} +```{r sp-6} # differentiate white matter from cortex based on annotations dlpfc_metadata <- as.data.frame(colData(dlpfc_spe)) %>% filter(layer_guess != "NA") %>% # remove missing annotations @@ -710,7 +716,7 @@ matter while the cortical grey matter predominantly contains neuronal cells. To further illustrate the difference, oligodendrocyte signature scores between cortical spots and white matter can be compared statistically. -```{r} +```{r sp-7} wilcox.test(Oligo_SPAROscore ~ Region, data = dlpfc_metadata) ``` The significant reduction in oligodendrocyte signature scores within grey matter @@ -763,6 +769,6 @@ usage vignette can be found here: [Advanced options](sparoscore_advanced_options.html) -```{r} +```{r session-info} sessionInfo() ```