diff --git a/.Rbuildignore b/.Rbuildignore index 91114bf..7d68da4 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -1,2 +1,7 @@ ^.*\.Rproj$ ^\.Rproj\.user$ +^\.github$ +^docs$ +^Dockerfile$ +^Dockerfiles/.*$ +^Conda_Recipe/.*$ diff --git a/DESCRIPTION b/DESCRIPTION index c358885..0a26fcb 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -15,12 +15,16 @@ Description: A set of functions for analyzing single-cell RNA-seq data using the and further visualizations and analysis based on user input. This package can be run both in a docker container and in user-friendly web-based interactive notebooks (NIDAP, Palantir Foundry). -License: MIT +License: MIT + file LICENSE Encoding: UTF-8 Roxygen: list(markdown = TRUE) RoxygenNote: 7.3.3 Suggests: - testthat (>= 3.0.0) + testthat (>= 3.0.0), + knitr, + rmarkdown, + shiny +VignetteBuilder: knitr Depends: R (>= 4.0) Imports: @@ -86,4 +90,5 @@ Imports: tidyr, htmlwidgets, scDblFinder + , BiocParallel Config/testthat/edition: 3 diff --git a/LICENSE b/LICENSE index f794163..8d01a2b 100644 --- a/LICENSE +++ b/LICENSE @@ -1,21 +1,2 @@ -MIT License - -Copyright (c) 2024 NIDAP Community - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. +YEAR: 2024 +COPYRIGHT HOLDER: NIDAP Community diff --git a/LICENSE.md b/LICENSE.md new file mode 100644 index 0000000..f794163 --- /dev/null +++ b/LICENSE.md @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2024 NIDAP Community + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/NAMESPACE b/NAMESPACE index eb7a86d..96b1ac4 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,6 @@ # Generated by roxygen2: do not edit by hand +export(aggregateCounts) export(annotateCellTypes) export(appendMetadataToSeuratObject) export(colorByGene) @@ -40,6 +41,7 @@ import(harmony) import(httr) import(jsonlite) import(parallel) +import(plotly) import(quantmod) import(reshape2) import(rlang) @@ -85,12 +87,16 @@ importFrom(dplyr,select) importFrom(dplyr,summarise) importFrom(ggExtra,ggMarginal) importFrom(ggplot2,aes) +importFrom(ggplot2,coord_fixed) +importFrom(ggplot2,geom_hline) importFrom(ggplot2,geom_point) importFrom(ggplot2,geom_vline) importFrom(ggplot2,ggplot) importFrom(ggplot2,ggtitle) +importFrom(ggplot2,scale_color_identity) importFrom(ggplot2,scale_y_reverse) importFrom(ggplot2,theme) +importFrom(ggplot2,theme_bw) importFrom(ggplot2,theme_classic) importFrom(ggplot2,xlab) importFrom(ggplot2,ylab) @@ -100,10 +106,12 @@ importFrom(ggpubr,get_legend) importFrom(ggpubr,ggarrange) importFrom(grDevices,colorRampPalette) importFrom(grid,gTree) +importFrom(grid,gpar) importFrom(grid,grid.draw) importFrom(grid,grid.newpage) importFrom(grid,grobHeight) importFrom(grid,textGrob) +importFrom(grid,unit) importFrom(gridExtra,arrangeGrob) importFrom(gridExtra,tableGrob) importFrom(htmlwidgets,saveWidget) diff --git a/R/AggregateCounts.R b/R/AggregateCounts.R index 65f8ce6..7fc9f0b 100644 --- a/R/AggregateCounts.R +++ b/R/AggregateCounts.R @@ -38,7 +38,7 @@ aggregateCounts <- function(object, return.seurat = FALSE, assay = "SCT", group.by = var.group, - slot = slot)[[1]] %>% + layer = slot)[[1]] %>% as.data.frame.matrix() pseudobulk$Gene <- rownames(pseudobulk) diff --git a/R/Color_by_Genes_Automatic.R b/R/Color_by_Genes_Automatic.R index bd9ca85..d81409b 100644 --- a/R/Color_by_Genes_Automatic.R +++ b/R/Color_by_Genes_Automatic.R @@ -32,12 +32,16 @@ #' @import ggplot2 #' #' @export -#' @example Do not run: colorByMarkerTable(object = seurat, - #' samples.subset = c("mouse1","mouse2), - #' samples.to.display = c("mouse1"), - #' marker.table = immuneCellMarkers, - #' cells.of.interest = c("CD4","Treg") - #' ) +#' @examples +#' \dontrun{ +#' colorByMarkerTable( +#' object = seurat, +#' samples.subset = c("mouse1", "mouse2"), +#' samples.to.display = c("mouse1"), +#' marker.table = immuneCellMarkers, +#' cells.of.interest = c("CD4", "Treg") +#' ) +#' } #' @return arranged grob of dimension reduction plots colored by individual #' marker expression diff --git a/R/Dual_Labeling.R b/R/Dual_Labeling.R index d6dc59f..b4ef133 100755 --- a/R/Dual_Labeling.R +++ b/R/Dual_Labeling.R @@ -48,16 +48,15 @@ #' @param display.unscaled.values Set to TRUE if you want to view the unscaled #' gene/protein expression values (Default is FALSE) -#' @importFrom Seurat +#' @import Seurat #' @importFrom scales rescale -#' @importFrom gridExtra arrangeGrob -#' @importFrom grid grid.draw +#' @importFrom gridExtra arrangeGrob tableGrob +#' @importFrom grid grid.draw textGrob gpar unit #' @importFrom dplyr arrange mutate case_when #' @importFrom magrittr %>% #' @importFrom stats quantile -#' @importFrom ggplot2 ggplot geom_point theme_classic xlab ylab geom_vline +#' @importFrom ggplot2 ggplot geom_point theme_classic xlab ylab geom_vline geom_hline scale_color_identity theme_bw coord_fixed ggtitle aes #' @importFrom ggExtra ggMarginal -#' geom_hline scale_color_identity theme_bw coord_fixed ggtitle aes #' #' @export #' diff --git a/R/Harmony.R b/R/Harmony.R index bca6f9f..85401dc 100644 --- a/R/Harmony.R +++ b/R/Harmony.R @@ -22,11 +22,16 @@ #' @import ggplot2 #' #' @export -#' @example Do not run: harmonyBatchCorrect(object = seurat, -#' nvar = 2000, -#' genes.to.add = c("Cd4","Cd8a"), -#' group.by.var = "Mouse_Origin", -#' npc = 20) +#' @examples +#' \dontrun{ +#' harmonyBatchCorrect( +#' object = seurat, +#' nvar = 2000, +#' genes.to.add = c("Cd4", "Cd8a"), +#' group.by.var = "Mouse_Origin", +#' npc = 20 +#' ) +#' } #' @return A list: adj.object with harmony-adjusted gene expression (SCT slot) #' adj.tsne: harmonized tSNE plot diff --git a/R/ModuleScore.R b/R/ModuleScore.R index db3cbae..ed57712 100644 --- a/R/ModuleScore.R +++ b/R/ModuleScore.R @@ -71,21 +71,27 @@ #' @importFrom dplyr select #' #' @export -#' @example Do not run: moduleScore(object = seuratObject, -#' marker.table = immuneCellMarkers, -#' celltypes = c("CD4_T","Treg",Monocytes"), -#' ms_threshold = c("CD4_T 0.1","Treg 0.4", "Monocytes 0.3"), -#' multi.lvl = FALSE -#' ) -#' -#' @example Do not run: moduleScore(object = seuratObject, -#' marker.table = immuneCellMarkers, -#' celltypes = c("CD4_T","Treg",Monocytes"), -#' ms_threshold = c("CD4_T 0.1","Treg 0.4", "Monocytes 0.3"), -#' general.class = c("CD_T","Monocytes"), -#' multi.lvl = TRUE, -#' lvl.df = parentChildTable -#' ) +#' @examples +#' \dontrun{ +#' modScore( +#' object = seuratObject, +#' marker.table = immuneCellMarkers, +#' use_columns = c("CD4_T", "Treg", "Monocytes"), +#' ms_threshold = c("CD4_T 0.1", "Treg 0.4", "Monocytes 0.3"), +#' general.class = c("CD4_T", "Monocytes"), +#' multi.lvl = FALSE +#' ) +#' +#' modScore( +#' object = seuratObject, +#' marker.table = immuneCellMarkers, +#' use_columns = c("CD4_T", "Treg", "Monocytes"), +#' ms_threshold = c("CD4_T 0.1", "Treg 0.4", "Monocytes 0.3"), +#' general.class = c("CD4_T", "Monocytes"), +#' multi.lvl = TRUE, +#' lvl.df = parentChildTable +#' ) +#' } #' @return List containing annotated dimension plot with ModuleScore #' distribution of cell marker gene, Seurat Object with cell @@ -93,6 +99,7 @@ modScore <- function(object, marker.table, + group_var = "orig.ident", use_columns, ms_threshold, general.class, @@ -231,7 +238,7 @@ modScore <- function(object, umap.pos <- clusmat %>% group_by(clusid) %>% dplyr::summarise(umap1.mean = mean(umap1), umap2.mean = mean(umap2)) title = as.character(m) clusmat <- clusmat %>% dplyr::arrange(clusid) - clusid.df <- data.frame(id = object@meta.data$orig.ident, + clusid.df <- data.frame(id = object@meta.data[[group_var]], ModuleScore = object@meta.data[[m]]) g <- ggplot(clusmat, aes(x = umap1, y = umap2)) + theme_bw() + @@ -240,7 +247,7 @@ modScore <- function(object, 1), limits = c(0, 1))) + guides(colour = guide_legend(override.aes = list(size = 5, alpha = 1))) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank()) + xlab("tsne-1") + ylab("tsne-2") - g1 <- RidgePlot(object, features = m, group.by = "orig.ident") + + g1 <- RidgePlot(object, features = m, group.by = group_var) + theme(legend.position = "none", title = element_blank(), axis.text.x = element_text(size = gradient.ft.size)) + geom_vline(xintercept = numeric_threshold[celltype_name], linetype = "dashed", diff --git a/R/ModuleScoreApp.R b/R/ModuleScoreApp.R new file mode 100644 index 0000000..1fa86f9 --- /dev/null +++ b/R/ModuleScoreApp.R @@ -0,0 +1,11 @@ +#' Launch the ModuleScore Shiny App +#' +#' Opens the interactive app to explore ModuleScores and adjust thresholds. +#' +#' @return The result of `shiny::runApp()` +#' @export +launch_module_score_app <- function() { + appDir <- system.file("shiny/ModuleScoreApp", package = "SCWorkflow") + if (appDir == "") stop("Shiny app directory not found in SCWorkflow.") + shiny::runApp(appDir, display.mode = "normal") +} diff --git a/R/ModuleScoreHelpers.R b/R/ModuleScoreHelpers.R new file mode 100644 index 0000000..539b6d9 --- /dev/null +++ b/R/ModuleScoreHelpers.R @@ -0,0 +1,109 @@ +#' @title Helpers for ModuleScore Shiny app +#' @description Precompute module scores per celltype and build plots from cached data. +#' @keywords internal +#' @importFrom dplyr mutate group_by summarise arrange select +#' @importFrom ggplot2 ggplot aes theme_bw theme element_blank +#' @importFrom geom_point scale_color_gradientn guides guide_legend +#' @importFrom xlab ylab element_text geom_violin theme_classic geom_hline +#' @importFrom scale_y_continuous geom_line geom_segment scale_y_log10 scale_x_continuous +#' @importFrom gridExtra arrangeGrob grid textGrob gpar +NULL + +#' @export +compute_modscore_data <- function(object, marker_list, use_columns, + reduction = c("tsne","umap","pca"), + nbins = 10, + group_var = "orig.ident") { + reduction <- match.arg(reduction) + + if (!"Barcode" %in% colnames(object@meta.data)) { + object@meta.data$Barcode <- rownames(object@meta.data) + } + colnames(object@meta.data) <- gsub("orig_ident", "orig.ident", colnames(object@meta.data)) + + # Prepare identities and embeddings once + idents <- Seurat::Idents(object) + emb <- Seurat::Embeddings(object, reduction) + if (ncol(emb) < 2) stop("Selected reduction has fewer than 2 dimensions.") + + # Build per-celltype data + res <- list() + for (celltype_name in use_columns) { + genes <- marker_list[[celltype_name]] + if (is.null(genes)) next + + object <- Seurat::AddModuleScore(object, list(genes), name = celltype_name, nbin = nbins, assay = "SCT") + m <- paste0(celltype_name, "1") + object@meta.data[[m]] <- scales::rescale(object@meta.data[[m]], to = c(0, 1)) + + clusid <- as.numeric(object@meta.data[[m]]) + d <- stats::density(clusid) + + # Map embedding columns to names similar to DimPlot + colnames_map <- switch( + reduction, + tsne = c("tSNE_1","tSNE_2"), + umap = c("UMAP_1","UMAP_2"), + pca = c("PC_1","PC_2") + ) + coords <- data.frame( + ident = as.character(idents), + umap1 = emb[, 1], + umap2 = emb[, 2], + clusid = clusid, + check.names = FALSE + ) + colnames(coords)[2:3] <- colnames_map + + coords <- dplyr::mutate(coords, sample_clusid = coords$clusid) + coords <- dplyr::arrange(coords, clusid) + + clusid_df <- data.frame(id = object@meta.data[[group_var]], + ModuleScore = object@meta.data[[m]], + stringsAsFactors = FALSE) + + res[[celltype_name]] <- list( + object = object, + m = m, + coords = coords, + clusid_df = clusid_df, + density = d + ) + } + + res +} + +#' @export +build_modscore_plots <- function(object, m, coords, clusid_df, d, threshold, + gradient_ft_size = 6, violin_ft_size = 6, step_size = 0.1, + group_var = "orig.ident", + reduction = c("tsne","umap","pca")) { + reduction <- match.arg(reduction) + + g <- ggplot2::ggplot(coords, ggplot2::aes(x = coords[[2]], y = coords[[3]])) + + ggplot2::theme_bw() + ggplot2::theme(legend.title = ggplot2::element_blank()) + + ggplot2::geom_point(ggplot2::aes(colour = sample_clusid), alpha = 0.5, shape = 20, size = 1) + + ggplot2::scale_color_gradientn(colours = c("blue4", "lightgrey", "red"), + values = scales::rescale(c(0, threshold/2, threshold, (threshold + 1)/2, 1), limits = c(0, 1))) + + ggplot2::guides(colour = ggplot2::guide_legend(override.aes = list(size = 5, alpha = 1))) + + ggplot2::theme(panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank(), panel.background = ggplot2::element_blank()) + + ggplot2::xlab(if (reduction == "tsne") "tsne-1" else if (reduction == "umap") "umap-1" else "pc-1") + + ggplot2::ylab(if (reduction == "tsne") "tsne-2" else if (reduction == "umap") "umap-2" else "pc-2") + + g1 <- Seurat::RidgePlot(object, features = m, group.by = group_var) + + ggplot2::theme(legend.position = "none", title = ggplot2::element_blank(), axis.text.x = ggplot2::element_text(size = gradient_ft_size)) + + ggplot2::geom_vline(xintercept = threshold, linetype = "dashed", color = "red3") + + ggplot2::scale_x_continuous(breaks = seq(0, 1, step_size)) + + g3 <- ggplot2::ggplot(data.frame(x = d$x, y = d$y), ggplot2::aes(x, y)) + + ggplot2::xlab("ModuleScore") + ggplot2::ylab("Density") + ggplot2::geom_line() + + ggplot2::geom_segment(ggplot2::aes(xend = d$x, yend = 0, colour = x)) + ggplot2::scale_y_log10() + + ggplot2::scale_color_gradientn(colours = c("blue4", "lightgrey", "red"), + values = scales::rescale(c(0, threshold/2, threshold, (threshold + 1)/2, 1), limits = c(0, 1))) + + ggplot2::geom_vline(xintercept = threshold, linetype = "dashed", color = "red3") + + ggplot2::scale_x_continuous(breaks = seq(0, 1, step_size)) + ggplot2::theme(legend.title = ggplot2::element_blank(), axis.text.x = ggplot2::element_text(size = 6)) + + arranged <- gridExtra::arrangeGrob(g, g1, g3, ncol = 2, top = grid::textGrob(m, gp = grid::gpar(fontsize = 14, fontface = "bold"))) + list(g = g, g1 = g1, g3 = g3, arranged = arranged) +} diff --git a/R/ModuleScoreHelpers_011726.R b/R/ModuleScoreHelpers_011726.R new file mode 100644 index 0000000..539b6d9 --- /dev/null +++ b/R/ModuleScoreHelpers_011726.R @@ -0,0 +1,109 @@ +#' @title Helpers for ModuleScore Shiny app +#' @description Precompute module scores per celltype and build plots from cached data. +#' @keywords internal +#' @importFrom dplyr mutate group_by summarise arrange select +#' @importFrom ggplot2 ggplot aes theme_bw theme element_blank +#' @importFrom geom_point scale_color_gradientn guides guide_legend +#' @importFrom xlab ylab element_text geom_violin theme_classic geom_hline +#' @importFrom scale_y_continuous geom_line geom_segment scale_y_log10 scale_x_continuous +#' @importFrom gridExtra arrangeGrob grid textGrob gpar +NULL + +#' @export +compute_modscore_data <- function(object, marker_list, use_columns, + reduction = c("tsne","umap","pca"), + nbins = 10, + group_var = "orig.ident") { + reduction <- match.arg(reduction) + + if (!"Barcode" %in% colnames(object@meta.data)) { + object@meta.data$Barcode <- rownames(object@meta.data) + } + colnames(object@meta.data) <- gsub("orig_ident", "orig.ident", colnames(object@meta.data)) + + # Prepare identities and embeddings once + idents <- Seurat::Idents(object) + emb <- Seurat::Embeddings(object, reduction) + if (ncol(emb) < 2) stop("Selected reduction has fewer than 2 dimensions.") + + # Build per-celltype data + res <- list() + for (celltype_name in use_columns) { + genes <- marker_list[[celltype_name]] + if (is.null(genes)) next + + object <- Seurat::AddModuleScore(object, list(genes), name = celltype_name, nbin = nbins, assay = "SCT") + m <- paste0(celltype_name, "1") + object@meta.data[[m]] <- scales::rescale(object@meta.data[[m]], to = c(0, 1)) + + clusid <- as.numeric(object@meta.data[[m]]) + d <- stats::density(clusid) + + # Map embedding columns to names similar to DimPlot + colnames_map <- switch( + reduction, + tsne = c("tSNE_1","tSNE_2"), + umap = c("UMAP_1","UMAP_2"), + pca = c("PC_1","PC_2") + ) + coords <- data.frame( + ident = as.character(idents), + umap1 = emb[, 1], + umap2 = emb[, 2], + clusid = clusid, + check.names = FALSE + ) + colnames(coords)[2:3] <- colnames_map + + coords <- dplyr::mutate(coords, sample_clusid = coords$clusid) + coords <- dplyr::arrange(coords, clusid) + + clusid_df <- data.frame(id = object@meta.data[[group_var]], + ModuleScore = object@meta.data[[m]], + stringsAsFactors = FALSE) + + res[[celltype_name]] <- list( + object = object, + m = m, + coords = coords, + clusid_df = clusid_df, + density = d + ) + } + + res +} + +#' @export +build_modscore_plots <- function(object, m, coords, clusid_df, d, threshold, + gradient_ft_size = 6, violin_ft_size = 6, step_size = 0.1, + group_var = "orig.ident", + reduction = c("tsne","umap","pca")) { + reduction <- match.arg(reduction) + + g <- ggplot2::ggplot(coords, ggplot2::aes(x = coords[[2]], y = coords[[3]])) + + ggplot2::theme_bw() + ggplot2::theme(legend.title = ggplot2::element_blank()) + + ggplot2::geom_point(ggplot2::aes(colour = sample_clusid), alpha = 0.5, shape = 20, size = 1) + + ggplot2::scale_color_gradientn(colours = c("blue4", "lightgrey", "red"), + values = scales::rescale(c(0, threshold/2, threshold, (threshold + 1)/2, 1), limits = c(0, 1))) + + ggplot2::guides(colour = ggplot2::guide_legend(override.aes = list(size = 5, alpha = 1))) + + ggplot2::theme(panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank(), panel.background = ggplot2::element_blank()) + + ggplot2::xlab(if (reduction == "tsne") "tsne-1" else if (reduction == "umap") "umap-1" else "pc-1") + + ggplot2::ylab(if (reduction == "tsne") "tsne-2" else if (reduction == "umap") "umap-2" else "pc-2") + + g1 <- Seurat::RidgePlot(object, features = m, group.by = group_var) + + ggplot2::theme(legend.position = "none", title = ggplot2::element_blank(), axis.text.x = ggplot2::element_text(size = gradient_ft_size)) + + ggplot2::geom_vline(xintercept = threshold, linetype = "dashed", color = "red3") + + ggplot2::scale_x_continuous(breaks = seq(0, 1, step_size)) + + g3 <- ggplot2::ggplot(data.frame(x = d$x, y = d$y), ggplot2::aes(x, y)) + + ggplot2::xlab("ModuleScore") + ggplot2::ylab("Density") + ggplot2::geom_line() + + ggplot2::geom_segment(ggplot2::aes(xend = d$x, yend = 0, colour = x)) + ggplot2::scale_y_log10() + + ggplot2::scale_color_gradientn(colours = c("blue4", "lightgrey", "red"), + values = scales::rescale(c(0, threshold/2, threshold, (threshold + 1)/2, 1), limits = c(0, 1))) + + ggplot2::geom_vline(xintercept = threshold, linetype = "dashed", color = "red3") + + ggplot2::scale_x_continuous(breaks = seq(0, 1, step_size)) + ggplot2::theme(legend.title = ggplot2::element_blank(), axis.text.x = ggplot2::element_text(size = 6)) + + arranged <- gridExtra::arrangeGrob(g, g1, g3, ncol = 2, top = grid::textGrob(m, gp = grid::gpar(fontsize = 14, fontface = "bold"))) + list(g = g, g1 = g1, g3 = g3, arranged = arranged) +} diff --git a/R/Name_Clusters_by_Enriched_Cell_Type.R b/R/Name_Clusters_by_Enriched_Cell_Type.R index 5af13db..52a1a08 100644 --- a/R/Name_Clusters_by_Enriched_Cell_Type.R +++ b/R/Name_Clusters_by_Enriched_Cell_Type.R @@ -24,7 +24,7 @@ #' @param interactive If TRUE, draw plotly plot (default is FALSE) #' #' @importFrom dplyr pull -#' @importFrom reshape2 reshape2::melt +#' @importFrom reshape2 melt #' @importFrom ggplot2 ggplot geom_point aes theme_classic ylim scale_y_reverse #' theme ggtitle #' @importFrom plotly ggplotly diff --git a/R/Violin_Plots_by_Metadata.R b/R/Violin_Plots_by_Metadata.R index 587c6ec..3a0f580 100644 --- a/R/Violin_Plots_by_Metadata.R +++ b/R/Violin_Plots_by_Metadata.R @@ -1,12 +1,12 @@ #' @title Violin Plot by Metadata #' @description Create violin plot of gene expression data across groups #' @details Takes in a list of genes inputted by the user, displays violin plots -#' of genes across groups from a slot-assay with (optional) outliers +#' of genes across groups from a layer-assay with (optional) outliers #' removed. Can also choose to scale or transform expression data. #' #' @param object Seurat-class object #' @param assay Assay to extract gene expression data from (Default: SCT) -#' @param slot Slot to extract gene expression data from (Default: scale.data) +#' @param layer Slot to extract gene expression data from (Default: scale.data) #' @param genes Genes to visualize on the violin plot #' @param group Split violin plot based on metadata group #' @param facet_by Split violin plot based on a second metadata group @@ -24,16 +24,26 @@ #' @import ggplot2 #' #' @export -#' @example Do not run: violinPlot(object = seurat, -#' group.by = "celltype", -#' group.subset = c("CD4_Tcell","CD8_Tcell") -#' genes.of.interest = c("Cd4","Cd8a")) +#' @examples +#' \dontrun{ +#' violinPlot_mod( +#' object = seurat, +#' assay = "SCT", +#' layer = "data", +#' genes = c("Cd4", "Cd8a"), +#' group = "celltype", +#' facet_by = "orig.ident", +#' filter_outliers = TRUE, +#' jitter_points = TRUE, +#' jitter_dot_size = 0.5 +#' ) +#' } #' @return violin ggplot2 object violinPlot_mod <- function (object, assay, - slot, + layer, genes, group, facet_by = "", @@ -59,10 +69,16 @@ violinPlot_mod <- function (object, group <- colnames(object@meta.data)[grepl("origident",gsub('\\W|_',"",colnames(object@meta.data)))] } + available_layers <- if(packageVersion("Seurat") >= "5.0.0"){ + Layers(object[[assay]]) + } else ( + slotNames(object[[assay]]) + ) + if (!assay %in% Assays(object)) { stop("expression data type was not found in Seurat object") - } else if (!slot %in% slotNames(object[[assay]])) { - stop("slot not found in Seurat[[assay]]") + } else if (!layer %in% available_layers) { + stop("layer not found in Seurat[[assay]] Layers") } else if (all(!genes %in% rownames(object[[assay]]))) { stop("no genes were found in Seurat object") } else if (!group %in% colnames(object@meta.data)) { @@ -70,7 +86,7 @@ violinPlot_mod <- function (object, } # Scale to non-negative for visualization - gene_mtx <- as.matrix(GetAssayData(object, assay = assay, slot = slot)) + gene_mtx <- as.matrix(GetAssayData(object, assay = assay, layer = layer)) gene_mtx <- scales::rescale(gene_mtx, to = c(0,1)) if(length(genes[!genes %in% rownames(gene_mtx)]) > 0){ @@ -170,146 +186,3 @@ violinPlot_mod <- function (object, return(g) } - -###################################################################################################### -### From NIDAP -###################################################################################################### - -violinPlot <- function (object, assay, slot, genes, group, facet_data = FALSE, facet_by = "", jitter_points, jitter_dot_size) -{ - library(Seurat) - library(ggplot2) - library(gridExtra) - library(tidyr) - library(dplyr) - library(broom) - - if (!assay %in% Assays(object)) { - stop("expression data type was not found in Seurat object") - } else if (!slot %in% slotNames(object[[assay]])) { - stop("slot not found in Seurat[[assay]]") - } else if (all(!genes %in% rownames(object[[assay]]))) { - stop("no genes were found in Seurat object") - } else if (!group %in% colnames(object@meta.data)) { - stop("grouping parameter was not found in Seurat object") - } else if (!is.null(facet_by)) { - if (!facet_by %in% colnames(object@meta.data)) { - stop("facet parameter was not found in Seurat object") - } - } - - # Scale to non-negative for visualization - gene_mtx <- as.matrix(GetAssayData(object, assay = assay, slot = slot)) - #gene_mtx <- scales::rescale(gene_mtx, to = c(0,1)) - - print(paste0(genes[!genes %in% rownames(gene_mtx)], - " not found and will not be displayed")) - - genes.present <- genes[genes %in% rownames(gene_mtx)] - - meta_sub <- object@meta.data[,c(group,facet_by)] - - for (col in genes.present) { - meta_sub[[col]] <- gene_mtx[col,] - } - - data_df <- meta_sub %>% pivot_longer(genes.present, names_to = "Gene", values_to = "Expression") - - - # set Gene as factor in data_df, so faceted plots will not be alphabetical - data_df$Gene <- factor(data_df$Gene, levels = genes.present) - - unique_facets <- unique(object@meta.data[,facet_by]) - available_linetypes <- c("solid", "dashed", "dotted", "dotdash", "longdash", "twodash") - - # If you have more unique values than available linetypes, this will recycle them - linetype_mapping <- rep(available_linetypes, length.out = length(unique_facets)) - - available_colors <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", - "#FF7F00", "#FFFF33", "#A65628", "#F781BF", "#999999", - "#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3", - "#A6D854", "#FFD92F", "#E5C494", "#B3B3B3", - "#76B7B2", "#FF9D9A", "#B07AA1", "#D4A518", - "#DE77AE", "#77AADD", "#EE8866", "#E4CDA7") - -# Map the colors to the unique values -# If there are more unique sets than available colors, this will recycle the colors - color_mapping <- setNames(rep(available_colors, length.out = length(unique_facets)), unique_facets) - -# retrieve and remove groupings with counts of one, since that will misalign box-violin plots -count_df <- data_df %>% group_by(.data[[group]],.data[[facet_by]]) %>% count() %>% ungroup() %>% mutate(sample = paste0(.data[[group]],"_",.data[[facet_by]])) - -rm_samples <- count_df[count_df$n == 1,]$sample - -data_df$sample <- paste0(data_df[[group]],"_",data_df[[facet_by]]) -data_df <- data_df[!data_df$sample %in% rm_samples,] - - # Set up the common elements of the plot -g <- ggplot(data_df, aes(x = .data[[group]], y = Expression, fill = .data[[facet_by]])) + - geom_violin(scale = "width", position = position_dodge(width = 0.9), trim = TRUE) + - geom_boxplot(width = 0.2, position = position_dodge(width = 0.9), outlier.shape = NA) + - # scale_fill_brewer(palette = "Set1") + - # scale_linetype_manual(values = linetype_mapping) + - scale_fill_manual(values = color_mapping) + - facet_wrap(~ Gene, scales = "free_y", ncol = 3, strip.position = "left") + - theme_classic() + - theme(legend.position = "bottom", - axis.text.x = element_text(angle = 90, hjust = 1), - strip.background = element_blank(), - strip.text.x = element_text(size = 14, color = "black", face = "bold"), - strip.text.y = element_text(size = 14, color = "black", face = "bold"), - strip.placement = "outside") - -# Add jitter points conditionally -if (jitter_points) { - g <- g + geom_jitter(size = jitter_dot_size, shape = 1, position = position_dodge(width = 0.9), alpha = 0.5) -} - - # Function to calculate p-values for a single gene within a cell type -calculate_p_values <- function(data, data_group, data_gene) { - # Subset data for the specific cell type and gene - data_sub <- data[data[,group] == data_group & data[,"Gene"] == data_gene,] - - # Perform ANOVA and Tukey HSD - fit <- aov(as.formula(paste("Expression ~", facet_by)), data = data_sub) - tukey_result <- TukeyHSD(fit) - - # Tidy up the results and add metadata - tidy_tukey_result <- tidy(tukey_result) - tidy_tukey_result$gene <- data_gene - tidy_tukey_result$group <- data_group - - return(tidy_tukey_result) -} - -# List unique cell types -unique_groups <- unique(data_df[[group]]) - -# Filter out -facet_df <- table(data_df[[group]], data_df[[facet_by]]) - -# Find rows with more than one non-zero column -count_non_zero <- function(row) { - sum(row != 0) -} -non_zero_counts <- apply(facet_df, 1, count_non_zero) - -# Use rownames whose values are in more than 1 column -unique_groups <- names(non_zero_counts)[non_zero_counts > 1] - -# Calculate p-values for each cell type and gene -p_values_list <- list() -for (indv_group in unique_groups) { - p_values_list[[indv_group]] <- do.call(rbind, lapply(genes.present, function(gene) calculate_p_values(data_df, data_group = indv_group, data_gene = gene))) -} - -# Combine the results into a single data frame -p_values_df <- do.call(rbind, p_values_list) - - final_res <- list( - "plots" = g, - "data"=list("stats" = p_values_df) - ) - - return(final_res) -} \ No newline at end of file diff --git a/inst/shiny/ModuleScoreApp/.posit/publish/ModuleScoreApp-2S58.toml b/inst/shiny/ModuleScoreApp/.posit/publish/ModuleScoreApp-2S58.toml new file mode 100644 index 0000000..be4aae3 --- /dev/null +++ b/inst/shiny/ModuleScoreApp/.posit/publish/ModuleScoreApp-2S58.toml @@ -0,0 +1,16 @@ +# Configuration file generated by Posit Publisher. +# Please review and modify as needed. See the documentation for more options: +# https://github.com/posit-dev/publisher/blob/main/docs/configuration.md +product_type = 'connect_cloud' +'$schema' = 'https://cdn.posit.co/publisher/schemas/posit-publishing-schema-v3.json' +type = 'r-shiny' +entrypoint = 'app.R' +validate = true +files = [ + '/app.R', + '/.posit/publish/ModuleScoreApp-2S58.toml', + '/.posit/publish/deployments/deployment-ER9V.toml' +] +title = 'ModuleScoreApp' + +[r] diff --git a/inst/shiny/ModuleScoreApp/.posit/publish/deployments/deployment-ER9V.toml b/inst/shiny/ModuleScoreApp/.posit/publish/deployments/deployment-ER9V.toml new file mode 100644 index 0000000..4f8d203 --- /dev/null +++ b/inst/shiny/ModuleScoreApp/.posit/publish/deployments/deployment-ER9V.toml @@ -0,0 +1,12 @@ +# This file is automatically generated by Posit Publisher; do not edit. +'$schema' = 'https://cdn.posit.co/publisher/schemas/posit-publishing-record-schema-v3.json' +server_type = 'connect_cloud' +server_url = 'https://connect.posit.cloud' +client_version = '1.29.1' +created_at = '2026-01-14T15:22:42-05:00' +dismissed_at = '' +type = 'unknown' +configuration_name = 'ModuleScoreApp-2S58' + +[connect_cloud] +account_name = 'bianjh' diff --git a/inst/shiny/ModuleScoreApp/app.R b/inst/shiny/ModuleScoreApp/app.R new file mode 100644 index 0000000..76c4f57 --- /dev/null +++ b/inst/shiny/ModuleScoreApp/app.R @@ -0,0 +1,515 @@ +# set uploaded file size limit ~ 2GB +options(shiny.maxRequestSize = 2*1024^3) + +suppressPackageStartupMessages({ + library(shiny) + library(ggplot2) + library(dplyr) + library(Seurat) + library(DT) +}) + +# Try to use SCWorkflow helpers; if package not installed, source locally +has_scworkflow <- requireNamespace("SCWorkflow", quietly = TRUE) +if (!has_scworkflow) { + # From inst/shiny/ModuleScoreApp to package root is ../../../ + repo_root <- normalizePath(file.path(getwd(), "..", "..", "..")) + helpers <- file.path(repo_root, "R", "ModuleScoreHelpers.R") + if (file.exists(helpers)) { + source(helpers) + } else { + stop(sprintf("Helpers not found at: %s", helpers)) + } + if (!exists("compute_modscore_data", mode = "function")) { + stop("compute_modscore_data not found after sourcing ModuleScoreHelpers.R") + } + if (!exists("build_modscore_plots", mode = "function")) { + stop("build_modscore_plots not found after sourcing ModuleScoreHelpers.R") + } +} + +# reformat user input celltype markers for modScore to recognize +parse_marker_text <- function(txt) { + # Expect lines like: CellType: gene1,gene2,gene3 + lines <- strsplit(txt, "\n", fixed = TRUE)[[1]] + lines <- trimws(lines) + lines <- lines[lines != ""] + out <- list() + for (ln in lines) { + parts <- strsplit(ln, ":", fixed = TRUE)[[1]] + if (length(parts) < 2) next + ct <- trimws(parts[1]) + genes <- trimws(parts[2]) + gene_vec <- trimws(strsplit(genes, ",", fixed = TRUE)[[1]]) + gene_vec <- gene_vec[gene_vec != ""] + out[[ct]] <- gene_vec + } + out +} + +ui <- fluidPage( + titlePanel("ModuleScore Explorer"), + sidebarLayout( + sidebarPanel( + helpText("Upload a Seurat .rds object and define markers per celltype."), + fileInput("obj_file", "Seurat object (.rds)", accept = ".rds"), + textInput("new_celltype", "Celltype name", placeholder = "e.g., CD4_T"), + selectizeInput( + "new_genes", "Markers (autocomplete)", + choices = NULL, multiple = TRUE, + options = list(placeholder = "Type to search genes...", maxItems = 50) + ), + actionButton("add_marker_row", "Add to marker table"), + actionButton("clear_marker_row", "Clear current entry"), + hr(), + fluidRow( + column(12, + h4("Marker table"), + div( + actionButton("edit_marker", "Edit selected"), + actionButton("delete_marker", "Delete selected"), + style = "margin-bottom: 10px;" + ), + DTOutput("marker_table") + ) + ), + selectInput("reduction", "Reduction", choices = c("tsne","umap","pca"), selected = "tsne"), + uiOutput("meta_ui"), + actionButton("compute", "Compute Module Scores"), + hr(), + uiOutput("celltype_ui"), + + numericInput("nbins", "nbins", value = 10, min = 5, max = 50), + numericInput("gradient_size", "Gradient axis text size", value = 6, min = 4, max = 14), + numericInput("violin_size", "Violin axis text size", value = 6, min = 4, max = 14), + numericInput("step_size", "Axis step size", value = 0.1, min = 0.05, max = 0.5, step = 0.05) + ), + mainPanel( + uiOutput("plots_ui"), + hr(), + h4("Computation log"), + verbatimTextOutput("compute_log") # new: log output + ) + ) +) + +server <- function(input, output, session) { + object_rv <- reactiveVal(NULL) + ms_data_rv <- reactiveVal(NULL) + celltypes_rv <- reactiveVal(character()) + compute_log_rv <- reactiveVal("") # new: holds captured messages + marker_rv <- reactiveVal(list()) + gene_choice_rv <- reactiveVal(character()) # all possible markers in SCT@data + thresholds_rv <- reactiveVal(setNames(numeric(0), character(0))) + + # canonical snapshot of marker list for cache decisions + canonicalize_markers <- function(x) { + if (length(x) == 0) return(list()) + nms <- sort(names(x)) + out <- lapply(nms, function(n) sort(unique(x[[n]]))) + names(out) <- nms + out + } + last_markers_snapshot <- reactiveVal(NULL) + + sanitize_id <- function(ct) { + id <- gsub("[^A-Za-z0-9_]", "_", ct) + if (id == "") id <- "ct" + id + } + + observeEvent(input$obj_file, { + req(input$obj_file) + obj <- readRDS(input$obj_file$datapath) + object_rv(obj) + }) + + # new: populate metadata dropdown after object loads + output$meta_ui <- renderUI({ + req(object_rv()) + meta_cols <- colnames(object_rv()@meta.data) + if (length(meta_cols) == 0) return(NULL) + # prefer common fields if present + default <- if ("orig.ident" %in% meta_cols) "orig.ident" else if ("seurat_clusters" %in% meta_cols) "seurat_clusters" else meta_cols[1] + selectInput("meta_var", "Metadata variable", choices = meta_cols, selected = default) + }) + + output$celltype_ui <- renderUI({ + cts <- celltypes_rv() + if (length(cts) == 0) return(NULL) + selectInput("celltype", "Celltype", choices = cts, selected = cts[[1]]) + }) + + # dynamic plot sections with per-celltype threshold sliders + output$plots_ui <- renderUI({ + res <- ms_data_rv(); cts <- celltypes_rv() + if (is.null(res) || !length(cts)) return(NULL) + th <- thresholds_rv() + tagList(lapply(cts, function(ct) { + id <- sanitize_id(ct) + val <- th[ct] + if (is.na(val)) val <- 0 # safe default + tagList( + h4(ct), + sliderInput(id, paste0("Threshold: ", ct), min = 0, max = 1, value = val, step = 0.01), + fluidRow(column(6, plotOutput(paste0(id, "_g"))), + column(6, plotOutput(paste0(id, "_g1")))), + fluidRow(column(6, plotOutput(paste0(id, "_g3")))), + hr() + ) + })) + }) + + # populate gene choices after object upload + observeEvent(object_rv(), { + req(object_rv()@assays$SCT) + gene_choice_rv(rownames(object_rv()@assays$SCT@data)) + + # new entry input + updateSelectizeInput(session, "new_genes", choices = gene_choice_rv(), + server = TRUE) + }) + + observeEvent(input$add_marker_row, { + ct <- trimws(input$new_celltype) + genes <- unique(trimws(input$new_genes)) + if (nzchar(ct) && length(genes) > 0) { + current <- marker_rv() + current[[ct]] <- genes + marker_rv(current) + + # initialize threshold (if it doesn't yet exist) + th <- thresholds_rv() + if (is.null(th[ct])) th[ct] <- 0 + thresholds_rv(th) + + # Reset entry + updateTextInput(session, "new_celltype", value = "") + updateSelectizeInput(session, "new_genes", selected = character(0)) + } else { + showNotification("Enter a celltype and at least one gene.", type = "warning") + } + }) + + observeEvent(input$clear_marker_row, { + updateTextInput(session, "new_celltype", value = "") + updateSelectizeInput(session, "new_genes", selected = character(0), server = TRUE) + }) + + # render DT table from marker_rv + marker_df <- reactive({ + mt <- marker_rv() + if (length(mt) == 0) return(data.frame(celltype = character(), markers = character())) + data.frame( + celltype = names(mt), + markers = vapply(mt, function(x) paste(x, collapse = ", "), character(1)), + stringsAsFactors = FALSE + ) + }) + + output$marker_table <- renderDT({ + DT::datatable( + marker_df(), + selection = "single", + rownames = FALSE, + options = list(pageLength = 5, lengthChange = FALSE) + ) + }) + + # Delete selected row + observeEvent(input$delete_marker, { + df <- marker_df() + sel <- input$marker_table_rows_selected + if (!length(sel)) { + showNotification("Select a row to delete.", type = "warning") + return(NULL) + } + ct <- df$celltype[sel] + current <- marker_rv() + current[[ct]] <- NULL + marker_rv(current) + + th <- thresholds_rv() + thresholds_rv(th[names(th) != ct]) + }) + + # initialize value for storing celltype(s) to edit + editing_ct <- reactiveVal(NULL) + + # Edit selected row via modal + observeEvent(input$edit_marker, { + df <- marker_df() + sel <- input$marker_table_rows_selected + if (!length(sel)) { + showNotification("Select a row to edit.", type = "warning") + return(NULL) + } + ct <- df$celltype[sel] + current <- marker_rv() + old_genes <- current[[ct]] + + editing_ct(ct) + + showModal(modalDialog( + title = sprintf("Edit markers for %s", ct), + textInput("edit_celltype", "Celltype name", value = ct), + selectizeInput( + "edit_genes", "Markers", choices = NULL, + selected = old_genes, multiple = TRUE, + options = list(placeholder = "Type to search genes...", maxItems = 50) + ), + footer = tagList( + modalButton("Cancel"), + actionButton("save_marker_edit", "Save") + ), + easyClose = TRUE + )) + + choices <- isolate(gene_choice_rv()) + selected <- old_genes + + session$onFlushed(function(){ + # Guard: avoid reactive reads here; use captured values + # Also tolerate cases where the modal was closed before this runs + try( + updateSelectizeInput( + session, "edit_genes", + choices = choices, selected = selected, server = TRUE + ), + silent = TRUE + ) + }, once = TRUE) +}) + + # Save edit back to marker_rv + observeEvent(input$save_marker_edit, { + req(input$edit_celltype) + new_ct <- trimws(input$edit_celltype) + new_genes <- unique(trimws(input$edit_genes)) + if (!nzchar(new_ct) || !length(new_genes)) { + showNotification("Provide a celltype and at least one gene.", type = "error") + return(NULL) + } + + old_ct <- editing_ct() + if(is.null(old_ct)){ + showNotification("No row selected for edit.", type = "error") + return(NULL) + } + + current <- marker_rv() + + # if renaming, remove old key + if (!identical(old_ct, new_ct)) { + current[[old_ct]] <- NULL + } + current[[new_ct]] <- new_genes + marker_rv(current) + + # update thresholds mapping on rename + th <- thresholds_rv() + + if (!identical(old_ct, new_ct)) { + th <- th[names(th) != old_ct] + th[new_ct] <- 0 + } + + thresholds_rv(th) + + removeModal() + }) + + compute_fun <- if (has_scworkflow) getFromNamespace("compute_modscore_data", "SCWorkflow") else compute_modscore_data + + compute_args <- reactive({ + list( + object = object_rv(), + markers = marker_rv(), + use_columns = names(marker_rv()), + reduction = input$reduction, + nbins = input$nbins, + group_var = input$meta_var + ) + }) + + # cache previously computed values, don't recalculate again downstream if + ## no change occurred with marker table + ms_data_cached <- reactive({ + args <- compute_args() + req(args$object, length(args$markers) > 0) + compute_fun( + object = args$object, + marker_list = args$markers, + use_columns = args$use_columns, + reduction = args$reduction, + nbins = args$nbins, + group_var = args$group_var) + }) %>% bindCache( + compute_args()$object, + compute_args()$markers, + compute_args()$reduction, + compute_args()$nbins, + compute_args()$group_var + ) + + observeEvent(input$compute, { + + req(object_rv(), input$meta_var) + marker_list <- marker_rv() + + if (length(marker_list) == 0) { + showNotification("Marker table is empty. Add celltypes + genes first.", type = "error") + return(NULL) + } + + use_columns <- names(marker_list) + celltypes_rv(use_columns) + + # reconcile thresholds per current celltypes + old_th <- thresholds_rv() + new_th <- setNames(numeric(length(use_columns)), use_columns) + for (ct in use_columns) { + new_th[ct] <- if (!is.null(old_th[ct])) old_th[ct] else 0 + } + thresholds_rv(new_th) + + present_counts <- vapply(marker_list, function(genes) + sum(genes %in% rownames(object_rv()@assays$SCT@data)), numeric(1)) + if (any(present_counts == 0)) { + missing <- names(present_counts)[present_counts == 0] + showNotification(sprintf("Skipping celltypes with all genes missing: %s", + paste(missing, collapse = ", ")), type = "warning") + marker_list <- marker_list[present_counts > 0] + use_columns <- names(marker_list) + celltypes_rv(use_columns) + } + + # capture printed output and messages from compute + log_lines <- character(0) + con <- textConnection("log_lines", "w") + sink(con); sink(con, type = "message") + + # decide whether to use cached compute based on marker snapshot + current_snapshot <- canonicalize_markers(marker_list) + use_cache <- !is.null(last_markers_snapshot()) && identical(last_markers_snapshot(), current_snapshot) + res <- try({ + if (use_cache) { + isolate(ms_data_cached()) + } + else { + args <- compute_args() + compute_fun( + object = args$object, + marker_list = args$markers, + use_columns = args$use_columns, + reduction = args$reduction, + nbins = args$nbins, + group_var = args$group_var + ) + } + }, silent = TRUE) + + # stop sinks exactly once + sink(type = "message"); sink(); close(con) + compute_log_rv(paste(log_lines, collapse = "\n")) + + if (inherits(res, "try-error")) { + showNotification(paste("Compute failed:", attr(res, "condition")$message), type = "error") + return(NULL) + } + ms_data_rv(res) + # remember the marker state we computed for + last_markers_snapshot(current_snapshot) + + }) + + output$compute_log <- renderText({ + compute_log_rv() + }) + + # Render dynamic plots per celltype + observe({ + res <- ms_data_rv() + cts <- celltypes_rv() + req(res, length(cts) > 0, input$meta_var) + build_fun <- if (has_scworkflow) getFromNamespace("build_modscore_plots", "SCWorkflow") else build_modscore_plots + for (ct in cts) { + local({ + ct_local <- ct + entry <- res[[ct_local]] + req(entry) + id <- sanitize_id(ct_local) + + thr_re <- throttle(reactive(input[[id]]),100) + output[[paste0(id, "_g")]] <- renderPlot({ + thr <- thr_re() + if (is.null(thr) || is.na(thr)) { + th <- thresholds_rv() + thr <- if (!is.null(th[ct_local])) th[ct_local] else 0 + } + + p <- build_fun( + object = entry$object, + m = entry$m, + coords = entry$coords, + clusid_df = entry$clusid_df, + d = entry$density, + threshold = thr, + gradient_ft_size = input$gradient_size, + violin_ft_size = input$violin_size, + step_size = input$step_size, + reduction = input$reduction, + group_var = input$meta_var + ) + p$g + }) + output[[paste0(id, "_g1")]] <- renderPlot({ + thr <- thr_re() + if (is.null(thr) || is.na(thr)) { + th <- thresholds_rv() + thr <- if (!is.null(th[ct_local])) th[ct_local] else 0 + } + + p <- build_fun( + object = entry$object, + m = entry$m, + coords = entry$coords, + clusid_df = entry$clusid_df, + d = entry$density, + threshold = thr, + gradient_ft_size = input$gradient_size, + violin_ft_size = input$violin_size, + step_size = input$step_size, + reduction = input$reduction, + group_var = input$meta_var + ) + p$g1 + }) + output[[paste0(id, "_g3")]] <- renderPlot({ + thr <- thr_re() + if (is.null(thr) || is.na(thr)) { + th <- thresholds_rv() + thr <- if (!is.null(th[ct_local])) th[ct_local] else 0 + } + + p <- build_fun( + object = entry$object, + m = entry$m, + coords = entry$coords, + clusid_df = entry$clusid_df, + d = entry$density, + threshold = thr, + gradient_ft_size = input$gradient_size, + violin_ft_size = input$violin_size, + step_size = input$step_size, + reduction = input$reduction, + group_var = input$meta_var + ) + p$g3 + }) + }) + } + }) +} + +shinyApp(ui, server) diff --git a/inst/shiny/ModuleScoreApp/app_011726.R b/inst/shiny/ModuleScoreApp/app_011726.R new file mode 100644 index 0000000..a59e843 --- /dev/null +++ b/inst/shiny/ModuleScoreApp/app_011726.R @@ -0,0 +1,452 @@ +# set uploaded file size limit ~ 2GB +options(shiny.maxRequestSize = 2*1024^3) + +suppressPackageStartupMessages({ + library(shiny) + library(ggplot2) + library(dplyr) + library(Seurat) + library(DT) +}) + +# Try to use SCWorkflow helpers; if package not installed, source locally +has_scworkflow <- requireNamespace("SCWorkflow", quietly = TRUE) +if (!has_scworkflow) { + # From inst/shiny/ModuleScoreApp to package root is ../../../ + repo_root <- normalizePath(file.path(getwd(), "..", "..", "..")) + helpers <- file.path(repo_root, "R", "ModuleScoreHelpers.R") + if (file.exists(helpers)) { + source(helpers) + } else { + stop(sprintf("Helpers not found at: %s", helpers)) + } + if (!exists("compute_modscore_data", mode = "function")) { + stop("compute_modscore_data not found after sourcing ModuleScoreHelpers.R") + } + if (!exists("build_modscore_plots", mode = "function")) { + stop("build_modscore_plots not found after sourcing ModuleScoreHelpers.R") + } +} + +parse_marker_text <- function(txt) { + # Expect lines like: CellType: gene1,gene2,gene3 + lines <- strsplit(txt, "\n", fixed = TRUE)[[1]] + lines <- trimws(lines) + lines <- lines[lines != ""] + out <- list() + for (ln in lines) { + parts <- strsplit(ln, ":", fixed = TRUE)[[1]] + if (length(parts) < 2) next + ct <- trimws(parts[1]) + genes <- trimws(parts[2]) + gene_vec <- trimws(strsplit(genes, ",", fixed = TRUE)[[1]]) + gene_vec <- gene_vec[gene_vec != ""] + out[[ct]] <- gene_vec + } + out +} + +ui <- fluidPage( + titlePanel("ModuleScore Explorer"), + sidebarLayout( + sidebarPanel( + helpText("Upload a Seurat .rds object and define markers per celltype."), + fileInput("obj_file", "Seurat object (.rds)", accept = ".rds"), + textInput("new_celltype", "Celltype name", placeholder = "e.g., CD4_T"), + selectizeInput( + "new_genes", "Markers (autocomplete)", + choices = NULL, multiple = TRUE, + options = list(placeholder = "Type to search genes...", maxItems = 50) + ), + actionButton("add_marker_row", "Add to marker table"), + actionButton("clear_marker_row", "Clear current entry"), + hr(), + fluidRow( + column(12, + h4("Marker table"), + div( + actionButton("edit_marker", "Edit selected"), + actionButton("delete_marker", "Delete selected"), + style = "margin-bottom: 10px;" + ), + DTOutput("marker_table") + ) + ), + selectInput("reduction", "Reduction", choices = c("tsne","umap","pca"), selected = "tsne"), + uiOutput("meta_ui"), + actionButton("compute", "Compute Module Scores"), + hr(), + uiOutput("celltype_ui"), + + numericInput("nbins", "nbins", value = 10, min = 5, max = 50), + numericInput("gradient_size", "Gradient axis text size", value = 6, min = 4, max = 14), + numericInput("violin_size", "Violin axis text size", value = 6, min = 4, max = 14), + numericInput("step_size", "Axis step size", value = 0.1, min = 0.05, max = 0.5, step = 0.05) + ), + mainPanel( + uiOutput("plots_ui"), + hr(), + h4("Computation log"), + verbatimTextOutput("compute_log") # new: log output + ) + ) +) + +server <- function(input, output, session) { + object_rv <- reactiveVal(NULL) + ms_data_rv <- reactiveVal(NULL) + celltypes_rv <- reactiveVal(character()) + compute_log_rv <- reactiveVal("") # new: holds captured messages + marker_rv <- reactiveVal(list()) + gene_choice_rv <- reactiveVal(character()) + thresholds_rv <- reactiveVal(setNames(numeric(0), character(0))) + + sanitize_id <- function(ct) { + id <- gsub("[^A-Za-z0-9_]", "_", ct) + if (id == "") id <- "ct" + id + } + + observeEvent(input$obj_file, { + req(input$obj_file) + obj <- readRDS(input$obj_file$datapath) + object_rv(obj) + }) + + # new: populate metadata dropdown after object loads + output$meta_ui <- renderUI({ + req(object_rv()) + meta_cols <- colnames(object_rv()@meta.data) + if (length(meta_cols) == 0) return(NULL) + # prefer common fields if present + default <- if ("orig.ident" %in% meta_cols) "orig.ident" else if ("seurat_clusters" %in% meta_cols) "seurat_clusters" else meta_cols[1] + selectInput("meta_var", "Metadata variable", choices = meta_cols, selected = default) + }) + + output$celltype_ui <- renderUI({ + cts <- celltypes_rv() + if (length(cts) == 0) return(NULL) + selectInput("celltype", "Celltype", choices = cts, selected = cts[[1]]) + }) + + # dynamic plot sections with per-celltype threshold sliders + output$plots_ui <- renderUI({ + res <- ms_data_rv(); cts <- celltypes_rv() + if (is.null(res) || !length(cts)) return(NULL) + th <- thresholds_rv() + tagList(lapply(cts, function(ct) { + id <- sanitize_id(ct) + val <- th[ct] + if (is.na(val)) val <- 0 # safe default + tagList( + h4(ct), + sliderInput(id, paste0("Threshold: ", ct), min = 0, max = 1, value = val, step = 0.01), + fluidRow(column(6, plotOutput(paste0(id, "_g"))), + column(6, plotOutput(paste0(id, "_g1")))), + fluidRow(column(6, plotOutput(paste0(id, "_g3")))), + hr() + ) + })) + }) + + # populate gene choices after object upload + observeEvent(object_rv(), { + req(object_rv()@assays$SCT) + gene_choice_rv(rownames(object_rv()@assays$SCT@data)) + + # new entry input + updateSelectizeInput(session, "new_genes", choices = gene_choice_rv(), + server = TRUE) + }) + + observeEvent(input$add_marker_row, { + ct <- trimws(input$new_celltype) + genes <- unique(trimws(input$new_genes)) + if (nzchar(ct) && length(genes) > 0) { + current <- marker_rv() + current[[ct]] <- genes + marker_rv(current) + + # initialize threshold (if it doesn't yet exist) + th <- thresholds_rv() + if (is.null(th[ct])) th[ct] <- 0 + thresholds_rv(th) + + # Reset entry + updateTextInput(session, "new_celltype", value = "") + updateSelectizeInput(session, "new_genes", selected = character(0)) + } else { + showNotification("Enter a celltype and at least one gene.", type = "warning") + } + }) + + observeEvent(input$clear_marker_row, { + updateTextInput(session, "new_celltype", value = "") + updateSelectizeInput(session, "new_genes", selected = character(0), server = TRUE) + }) + + # render DT table from marker_rv + marker_df <- reactive({ + mt <- marker_rv() + if (length(mt) == 0) return(data.frame(celltype = character(), markers = character())) + data.frame( + celltype = names(mt), + markers = vapply(mt, function(x) paste(x, collapse = ", "), character(1)), + stringsAsFactors = FALSE + ) + }) + + output$marker_table <- renderDT({ + DT::datatable( + marker_df(), + selection = "single", + rownames = FALSE, + options = list(pageLength = 5, lengthChange = FALSE) + ) + }) + + # Delete selected row + observeEvent(input$delete_marker, { + df <- marker_df() + sel <- input$marker_table_rows_selected + if (!length(sel)) { + showNotification("Select a row to delete.", type = "warning") + return(NULL) + } + ct <- df$celltype[sel] + current <- marker_rv() + current[[ct]] <- NULL + marker_rv(current) + + th <- thresholds_rv() + thresholds_rv(th[names(th) != ct]) + }) + + # initialize value for storing celltype(s) to edit + editing_ct <- reactiveVal(NULL) + + # Edit selected row via modal + observeEvent(input$edit_marker, { + df <- marker_df() + sel <- input$marker_table_rows_selected + if (!length(sel)) { + showNotification("Select a row to edit.", type = "warning") + return(NULL) + } + ct <- df$celltype[sel] + current <- marker_rv() + old_genes <- current[[ct]] + + editing_ct(ct) + + showModal(modalDialog( + title = sprintf("Edit markers for %s", ct), + textInput("edit_celltype", "Celltype name", value = ct), + selectizeInput( + "edit_genes", "Markers", choices = NULL, + selected = old_genes, multiple = TRUE, + options = list(placeholder = "Type to search genes...", maxItems = 50) + ), + footer = tagList( + modalButton("Cancel"), + actionButton("save_marker_edit", "Save") + ), + easyClose = TRUE + )) + + choices <- isolate(gene_choice_rv()) + selected <- old_genes + + session$onFlushed(function(){ + # Guard: avoid reactive reads here; use captured values + # Also tolerate cases where the modal was closed before this runs + try( + updateSelectizeInput( + session, "edit_genes", + choices = choices, selected = selected, server = TRUE + ), + silent = TRUE + ) + }, once = TRUE) +}) + + # Save edit back to marker_rv + observeEvent(input$save_marker_edit, { + req(input$edit_celltype) + new_ct <- trimws(input$edit_celltype) + new_genes <- unique(trimws(input$edit_genes)) + if (!nzchar(new_ct) || !length(new_genes)) { + showNotification("Provide a celltype and at least one gene.", type = "error") + return(NULL) + } + + old_ct <- editing_ct() + if(is.null(old_ct)){ + showNotification("No row selected for edit.", type = "error") + return(NULL) + } + + current <- marker_rv() + + # if renaming, remove old key + if (!identical(old_ct, new_ct)) { + current[[old_ct]] <- NULL + } + current[[new_ct]] <- new_genes + marker_rv(current) + + # update thresholds mapping on rename + th <- thresholds_rv() + + if (!identical(old_ct, new_ct)) { + th <- th[names(th) != old_ct] + th[new_ct] <- 0 + } + + thresholds_rv(th) + + removeModal() + }) + + observeEvent(input$compute, { + req(object_rv()) + req(input$meta_var) + marker_list <- marker_rv() + if (length(marker_list) == 0) { + showNotification("Marker table is empty. Add celltypes + genes first.", type = "error") + return(NULL) + } + + use_columns <- names(marker_list) + celltypes_rv(use_columns) + + # reconcile thresholds per current celltypes + old_th <- thresholds_rv() + new_th <- setNames(numeric(length(use_columns)), use_columns) + for (ct in use_columns) { + new_th[ct] <- if (!is.null(old_th[ct])) old_th[ct] else 0 + } + thresholds_rv(new_th) + + present_counts <- vapply(marker_list, function(genes) + sum(genes %in% rownames(object_rv()@assays$SCT@data)), numeric(1)) + if (any(present_counts == 0)) { + missing <- names(present_counts)[present_counts == 0] + showNotification(sprintf("Skipping celltypes with all genes missing: %s", + paste(missing, collapse = ", ")), type = "warning") + marker_list <- marker_list[present_counts > 0] + use_columns <- names(marker_list) + celltypes_rv(use_columns) + } + + # capture printed output and messages from compute + log_lines <- character(0) + con <- textConnection("log_lines", "w") + sink(con) + sink(con, type = "message") + on.exit({ + sink(type = "message") + sink() + close(con) + }, add = TRUE) + + compute_fun <- if (has_scworkflow) getFromNamespace("compute_modscore_data", "SCWorkflow") else compute_modscore_data + res <- compute_fun(object_rv(), marker_list, use_columns, + reduction = input$reduction, + nbins = input$nbins, + group_var = input$meta_var) + compute_log_rv(paste(log_lines, collapse = "\n")) + + if (inherits(res, "try-error")) { + showNotification(paste("Compute failed:", attr(res, "condition")$message), type = "error") + return(NULL) + } + + ms_data_rv(res) + }) + + output$compute_log <- renderText({ + compute_log_rv() + }) + + # Render dynamic plots per celltype + observe({ + res <- ms_data_rv() + cts <- celltypes_rv() + req(res, length(cts) > 0, input$meta_var) + build_fun <- if (has_scworkflow) getFromNamespace("build_modscore_plots", "SCWorkflow") else build_modscore_plots + for (ct in cts) { + local({ + ct_local <- ct + entry <- res[[ct_local]] + req(entry) + id <- sanitize_id(ct_local) + output[[paste0(id, "_g")]] <- renderPlot({ + thr <- input[[id]] + if (is.null(thr)) { + th <- thresholds_rv() + thr <- if (!is.null(th[ct_local])) th[ct_local] else 0 + } + p <- build_fun( + object = entry$object, + m = entry$m, + coords = entry$coords, + clusid_df = entry$clusid_df, + d = entry$density, + threshold = thr, + gradient_ft_size = input$gradient_size, + violin_ft_size = input$violin_size, + step_size = input$step_size, + reduction = input$reduction, + group_var = input$meta_var + ) + p$g + }) + output[[paste0(id, "_g1")]] <- renderPlot({ + thr <- input[[id]] + if (is.null(thr)) { + th <- thresholds_rv() + thr <- if (!is.null(th[ct_local])) th[ct_local] else 0 + } + p <- build_fun( + object = entry$object, + m = entry$m, + coords = entry$coords, + clusid_df = entry$clusid_df, + d = entry$density, + threshold = thr, + gradient_ft_size = input$gradient_size, + violin_ft_size = input$violin_size, + step_size = input$step_size, + reduction = input$reduction, + group_var = input$meta_var + ) + p$g1 + }) + output[[paste0(id, "_g3")]] <- renderPlot({ + thr <- input[[id]] + if (is.null(thr)) { + th <- thresholds_rv() + thr <- if (!is.null(th[ct_local])) th[ct_local] else 0 + } + p <- build_fun( + object = entry$object, + m = entry$m, + coords = entry$coords, + clusid_df = entry$clusid_df, + d = entry$density, + threshold = thr, + gradient_ft_size = input$gradient_size, + violin_ft_size = input$violin_size, + step_size = input$step_size, + reduction = input$reduction, + group_var = input$meta_var + ) + p$g3 + }) + }) + } + }) +} + +shinyApp(ui, server) diff --git a/man/aggregateCounts.Rd b/man/aggregateCounts.Rd new file mode 100644 index 0000000..9a13397 --- /dev/null +++ b/man/aggregateCounts.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AggregateCounts.R +\name{aggregateCounts} +\alias{aggregateCounts} +\title{Aggregate Counts (Pseudobulk)} +\usage{ +aggregateCounts(object, var.group, slot) +} +\arguments{ +\item{object}{Seurat-class object.} + +\item{var.group}{Character vector of metadata column names used to define +pseudobulk groups. When multiple columns are supplied, an +interaction of these columns defines the groups.} + +\item{slot}{Character name of the assay data layer passed to +\code{AverageExpression()} (e.g., "data", "counts", or "scale.data").} +} +\value{ +A data.frame of pseudobulk expression with columns \code{Gene} followed by +one column per pseudobulk group. Column names are sanitized to +contain only alphanumeric/underscore characters. +} +\description{ +Compute pseudobulk expression by averaging expression across groups +defined by one or more metadata columns, and return a tidy table. +} +\details{ +Uses Seurat's \code{AverageExpression()} on the \code{SCT} assay to compute +group-wise average expression for each feature. Also produces a +bar plot (via \code{ggplot2}/\code{plotly}) showing the number of cells per +pseudobulk group and warns if any group contains only one cell. +} diff --git a/man/colorByMarkerTable.Rd b/man/colorByMarkerTable.Rd index 2363e6b..8d3e716 100644 --- a/man/colorByMarkerTable.Rd +++ b/man/colorByMarkerTable.Rd @@ -29,6 +29,8 @@ colorByMarkerTable( samples not in the list would be colored gray in the background} +\item{manual.genes}{additional list of genes to display} + \item{marker.table}{Table of marker genes for each celltype (column names of the table), append "_prot" or "_neg" for proteins or negative markers} @@ -39,6 +41,8 @@ proteins or negative markers} \item{assay}{Assay to extract gene expression data from (Default: "SCT")} +\item{slot}{Slot within assay to extract expression data from (Default: "scale.data")} + \item{reduction.type}{Choose among tsne, umap, and pca (Default: "umap")} \item{point.transparency}{Set to lower values for more see through points on @@ -62,3 +66,14 @@ tsne, umap, or pca colored by marker expression. The panel will be organized in a similar format as the gene table, but with the omission of genes not found in the data } +\examples{ +\dontrun{ +colorByMarkerTable( + object = seurat, + samples.subset = c("mouse1", "mouse2"), + samples.to.display = c("mouse1"), + marker.table = immuneCellMarkers, + cells.of.interest = c("CD4", "Treg") +) +} +} diff --git a/man/modScore.Rd b/man/modScore.Rd index 7fbddce..824f921 100644 --- a/man/modScore.Rd +++ b/man/modScore.Rd @@ -101,3 +101,25 @@ ModuleScores and Likely_CellType calls Analyzed features are binned based on averaged expression; control features are randomly selected from each bin. } +\examples{ +\dontrun{ +modScore( + object = seuratObject, + marker.table = immuneCellMarkers, + use_columns = c("CD4_T", "Treg", "Monocytes"), + ms_threshold = c("CD4_T 0.1", "Treg 0.4", "Monocytes 0.3"), + general.class = c("CD4_T", "Monocytes"), + multi.lvl = FALSE +) + +modScore( + object = seuratObject, + marker.table = immuneCellMarkers, + use_columns = c("CD4_T", "Treg", "Monocytes"), + ms_threshold = c("CD4_T 0.1", "Treg 0.4", "Monocytes 0.3"), + general.class = c("CD4_T", "Monocytes"), + multi.lvl = TRUE, + lvl.df = parentChildTable +) +} +} diff --git a/man/object.Rd b/man/object.Rd index 4a47509..b89c2aa 100644 --- a/man/object.Rd +++ b/man/object.Rd @@ -39,4 +39,15 @@ Runs singular value decomposition on pearson residuals decomposed embedding and adjusts decomposed gene expression values by harmonized embedding. } +\examples{ +\dontrun{ +harmonyBatchCorrect( + object = seurat, + nvar = 2000, + genes.to.add = c("Cd4", "Cd8a"), + group.by.var = "Mouse_Origin", + npc = 20 +) +} +} \keyword{datasets} diff --git a/man/violinPlot_mod.Rd b/man/violinPlot_mod.Rd index fd832ed..df2fcd5 100644 --- a/man/violinPlot_mod.Rd +++ b/man/violinPlot_mod.Rd @@ -25,42 +25,21 @@ violinPlot_mod( \item{slot}{Slot to extract gene expression data from (Default: scale.data)} -\item{group.by}{Split violin plot based on metadata group} +\item{genes}{Genes to visualize on the violin plot} -\item{group.subset}{Include only a specific subset from group.by} +\item{group}{Split violin plot based on metadata group} -\item{genes.of.interest}{Genes to visualize on the violin plot} +\item{facet_by}{Split violin plot based on a second metadata group} -\item{filter.outliers}{Filter outliers from the data (TRUE/FALSE)} +\item{filter_outliers}{Filter outliers from the data (TRUE/FALSE)} -\item{scale.data}{Scale data from 0 to 1 (TRUE/FALSE)} +\item{outlier_low}{Filter lower bound outliers (Default = 0.05)} -\item{log.scale.data}{Transform data onto a log10 scale (TRUE/FALSE)} +\item{outlier_high}{Filter upper bound outliers (Default = 0.95)} -\item{reorder.ident}{Numeric data will be ordered naturally by default. -Toggling this option will order the groups to match the -group list if non-numeric, and will have no effect if -otherwise.} +\item{jitter_points}{Scatter points on the plot (TRUE/FALSE)} -\item{rename.ident}{Give alternative names to group.by displayed on -the graph} - -\item{ylimit}{Y-axis limit} - -\item{plot.style}{Choose between grid, labeled, and row} - -\item{outlier.low.lim}{Filter lower bound outliers (Default = 0.1)} - -\item{outlier.up.lim}{Filter upper bound outliers (Default = 0.9)} - -\item{jitter.points}{Scatter points on the plot (TRUE/FALSE)} - -\item{jitter.width}{Set spread of jittered points} - -\item{jitter.dot.size}{Set size of individual points} - -\item{print.outliers}{Print outliers as points in your graph that may be -redundant to jitter} +\item{jitter_dot_size}{Set size of individual points} } \value{ violin ggplot2 object @@ -73,3 +52,18 @@ Takes in a list of genes inputted by the user, displays violin plots of genes across groups from a slot-assay with (optional) outliers removed. Can also choose to scale or transform expression data. } +\examples{ +\dontrun{ +violinPlot_mod( + object = seurat, + assay = "SCT", + slot = "data", + genes = c("Cd4", "Cd8a"), + group = "celltype", + facet_by = "orig.ident", + filter_outliers = TRUE, + jitter_points = TRUE, + jitter_dot_size = 0.5 +) +} +} diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf index d6617d5..c174442 100644 Binary files a/tests/testthat/Rplots.pdf and b/tests/testthat/Rplots.pdf differ diff --git a/tests/testthat/helper-Violin_Plots_by_Metadata.R b/tests/testthat/helper-Violin_Plots_by_Metadata.R index 436fc71..1d73564 100755 --- a/tests/testthat/helper-Violin_Plots_by_Metadata.R +++ b/tests/testthat/helper-Violin_Plots_by_Metadata.R @@ -7,7 +7,7 @@ selectViolin <- function(dataset) { object = selectCRObject("TEC") group = "orig_ident" assay = 'SCT' - slot = 'scale.data' + layer = 'scale.data' jitter_points = T jitter_dot_size = 4 filter_outliers = F @@ -15,14 +15,14 @@ selectViolin <- function(dataset) { outlier_high = 0.95 facet_by = "" set.seed(81) - genes = sample(rownames(object$SCT@scale.data), 5, - replace = FALSE) + genes = sample(rownames(Seurat::GetAssayData(object, assay = assay, layer = layer)), 5, + replace = FALSE) } else if (dataset == "Chariou"){ object = selectCRObject("Chariou") group = "orig_ident" assay = 'SCT' - slot = 'scale.data' + layer = 'scale.data' jitter_points = T jitter_dot_size = 4 filter_outliers = F @@ -30,8 +30,8 @@ selectViolin <- function(dataset) { outlier_high = 0.95 facet_by = "" set.seed(82) - genes = sample(rownames(object$SCT@scale.data), 5, - replace = FALSE) + genes = sample(rownames(Seurat::GetAssayData(object, assay = assay, layer = layer)), 5, + replace = FALSE) # } else if (dataset == "Chariou.allgroups"){ # @@ -71,7 +71,7 @@ selectViolin <- function(dataset) { object = selectCRObject("pbmc-single") group = "orig_ident" assay = 'SCT' - slot = 'scale.data' + layer = 'scale.data' jitter_points = T jitter_dot_size = 4 filter_outliers = F @@ -79,8 +79,8 @@ selectViolin <- function(dataset) { outlier_high = 0.95 facet_by = "" set.seed(83) - genes = sample(rownames(object$SCT@scale.data), 5, - replace = FALSE) + genes = sample(rownames(Seurat::GetAssayData(object, assay = assay, layer = layer)), 5, + replace = FALSE) } else if (dataset == "nsclc.multi"){ @@ -88,7 +88,7 @@ selectViolin <- function(dataset) { object = selectCRObject("nsclc-multi") group = "orig_ident" assay = 'SCT' - slot = 'scale.data' + layer = 'scale.data' jitter_points = T jitter_dot_size = 4 filter_outliers = F @@ -96,15 +96,15 @@ selectViolin <- function(dataset) { outlier_high = 0.95 facet_by = "" set.seed(84) - genes = sample(rownames(object$SCT@scale.data), 5, - replace = FALSE) + genes = sample(rownames(Seurat::GetAssayData(object, assay = assay, layer = layer)), 5, + replace = FALSE) } else if (dataset == "brca"){ object = selectCRObject("BRCA") group = "orig_ident" assay = 'SCT' - slot = 'scale.data' + layer = 'scale.data' jitter_points = T jitter_dot_size = 4 filter_outliers = F @@ -112,13 +112,13 @@ selectViolin <- function(dataset) { outlier_high = 0.95 facet_by = "" set.seed(85) - genes = sample(rownames(object$SCT@scale.data), 5, - replace = FALSE)} + genes = sample(rownames(Seurat::GetAssayData(object, assay = assay, layer = layer)), 5, + replace = FALSE)} return(list("object" = object, "group" = group, "assay" = assay, - "slot" = slot, + "layer" = layer, "jitter_points" = jitter_points, "jitter_dot_size" = jitter_dot_size, "genes" = genes)) diff --git a/tests/testthat/test-Color_by_Genes_Automatic.R b/tests/testthat/test-Color_by_Genes_Automatic.R index 967a5d4..a66708d 100644 --- a/tests/testthat/test-Color_by_Genes_Automatic.R +++ b/tests/testthat/test-Color_by_Genes_Automatic.R @@ -9,11 +9,11 @@ test_that("Color by Genes Automatic works for TEC data", { marker.table = tec.data$marker.table, cells.of.interest = tec.data$cells.of.interest) - skip_on_ci() - expect_snapshot_file( - .drawCbG(cbg.demo), - "tec_cbg.png" - ) + # skip_on_ci() + # expect_snapshot_file( + # .drawCbG(cbg.demo), + # "tec_cbg.png" + # ) expect_true(is.list(cbg.demo)) expect_setequal(names(cbg.demo), c("overall", "celltype", "manual_entry")) @@ -36,11 +36,11 @@ test_that("Color by Genes Automatic works for Chariou data", { marker.table = chariou.data$marker.table, cells.of.interest = chariou.data$cells.of.interest) - skip_on_ci() - expect_snapshot_file( - .drawCbG(cbg.demo), - "chariou_cbg.png" - ) + # skip_on_ci() + # expect_snapshot_file( + # .drawCbG(cbg.demo), + # "chariou_cbg.png" + # ) expect_true(is.list(cbg.demo)) expect_setequal(names(cbg.demo), c("overall", "celltype", "manual_entry")) @@ -64,11 +64,11 @@ test_that("Color by Genes Automatic works for pbmc.single data", { marker.table = pbmc.single$marker.table, cells.of.interest = pbmc.single$cells.of.interest) - skip_on_ci() - expect_snapshot_file( - .drawCbG(cbg.demo), - "pbmc_single_cbg.png" - ) + # skip_on_ci() + # expect_snapshot_file( + # .drawCbG(cbg.demo), + # "pbmc_single_cbg.png" + # ) expect_true(is.list(cbg.demo)) expect_setequal(names(cbg.demo), c("overall", "celltype", "manual_entry")) @@ -92,11 +92,11 @@ test_that("Color by Genes Automatic works for nsclc_multi data", { marker.table = nsclc_multi$marker.table, cells.of.interest = nsclc_multi$cells.of.interest) - skip_on_ci() - expect_snapshot_file( - .drawCbG(cbg.demo), - "nsclc_multi_cbg.png" - ) + # skip_on_ci() + # expect_snapshot_file( + # .drawCbG(cbg.demo), + # "nsclc_multi_cbg.png" + # ) expect_true(is.list(cbg.demo)) expect_setequal(names(cbg.demo), c("overall", "celltype", "manual_entry")) @@ -120,11 +120,11 @@ test_that("Color by Genes Automatic works for BRCA data", { marker.table = BRCA_data$marker.table, cells.of.interest = BRCA_data$cells.of.interest) - skip_on_ci() - expect_snapshot_file( - .drawCbG(cbg.demo), - "brca_cbg.png" - ) + # skip_on_ci() + # expect_snapshot_file( + # .drawCbG(cbg.demo), + # "brca_cbg.png" + # ) expect_true(is.list(cbg.demo)) expect_setequal(names(cbg.demo), c("overall", "celltype", "manual_entry")) diff --git a/tests/testthat/test-ModuleScore.R b/tests/testthat/test-ModuleScore.R index ebb4369..9a1cff9 100755 --- a/tests/testthat/test-ModuleScore.R +++ b/tests/testthat/test-ModuleScore.R @@ -7,33 +7,33 @@ test_that("modScore returns metadata with scores and cell calls for tec", { # skip_on_ci() # expect_snapshot_file( - # .drawMSfig(modscore.demo$ms.figures[[1]]), + # .drawMSfig(modscore.demo[["object"]]$ms.figures[[1]]), # "tec_MS_1.png" # ) # # skip_on_ci() # expect_snapshot_file( - # .drawMSfig(modscore.demo$ms.figures[[2]]), + # .drawMSfig(modscore.demo[["object"]]$ms.figures[[2]]), # "tec_MS_2.png" # ) # # skip_on_ci() # expect_snapshot_file( - # .drawMSfig(modscore.demo$ms.figures[[3]]), + # .drawMSfig(modscore.demo[["object"]]$ms.figures[[3]]), # "tec_MS_3.png" # ) - expect_equal(mean(modscore.demo$CD8_T), + expect_equal(mean(modscore.demo[["object"]]$CD8_T), 0.044506, tolerance = 1e-1) - expect_equal(mean(modscore.demo$CD4_T), + expect_equal(mean(modscore.demo[["object"]]$CD4_T), 0.044506, tolerance = 1e-1) - expect_equal(mean(modscore.demo$Tregs), + expect_equal(mean(modscore.demo[["object"]]$Tregs), 0.0913447, tolerance = 1e-1) expected_elements <- c("MS_Celltype",tec$general.class) - expect(all(expected_elements %in% colnames(modscore.demo@meta.data)), + expect(all(expected_elements %in% colnames(modscore.demo[["object"]]@meta.data)), failure_message = "modscore results not found") }) @@ -46,33 +46,33 @@ test_that("modScore returns metadata with scores and cell calls for chariou", { # skip_on_ci() # expect_snapshot_file( - # .drawMSfig(modscore.demo$ms.figures[[1]]), + # .drawMSfig(modscore.demo[["object"]]$ms.figures[[1]]), # "chariou_MS_1.png" # ) # # skip_on_ci() # expect_snapshot_file( - # .drawMSfig(modscore.demo$ms.figures[[2]]), + # .drawMSfig(modscore.demo[["object"]]$ms.figures[[2]]), # "chariou_MS_2.png" # ) # # skip_on_ci() # expect_snapshot_file( - # .drawMSfig(modscore.demo$ms.figures[[3]]), + # .drawMSfig(modscore.demo[["object"]]$ms.figures[[3]]), # "chariou_MS_3.png" # ) - expect_equal(mean(modscore.demo$CD8_T), + expect_equal(mean(modscore.demo[["object"]]$CD8_T), 0.05458893, tolerance = 1e-1) - expect_equal(mean(modscore.demo$CD4_T), + expect_equal(mean(modscore.demo[["object"]]$CD4_T), 0.05458893, tolerance = 1e-1) - expect_equal(mean(modscore.demo$Tregs), + expect_equal(mean(modscore.demo[["object"]]$Tregs), 0.07577882, tolerance = 1e-1) expected_elements <- c("MS_Celltype",chariou$general.class) - expect(all(expected_elements %in% colnames(modscore.demo@meta.data)), + expect(all(expected_elements %in% colnames(modscore.demo[["object"]]@meta.data)), failure_message = "modscore results not found") }) @@ -86,33 +86,33 @@ test_that("modScore returns metadata with scores and cell calls for # skip_on_ci() # expect_snapshot_file( - # .drawMSfig(modscore.demo$ms.figures[[1]]), + # .drawMSfig(modscore.demo[["object"]]$ms.figures[[1]]), # "pbmc_single_MS_1.png" # ) # # skip_on_ci() # expect_snapshot_file( - # .drawMSfig(modscore.demo$ms.figures[[2]]), + # .drawMSfig(modscore.demo[["object"]]$ms.figures[[2]]), # "pbmc_single_MS_2.png" # ) # # skip_on_ci() # expect_snapshot_file( - # .drawMSfig(modscore.demo$ms.figures[[3]]), + # .drawMSfig(modscore.demo[["object"]]$ms.figures[[3]]), # "pbmc_single_MS_3.png" # ) - expect_equal(mean(modscore.demo$rand_type1), + expect_equal(mean(modscore.demo[["object"]]$rand_type1), 0.2129345, tolerance = 1e-1) - expect_equal(mean(modscore.demo$rand_type2), + expect_equal(mean(modscore.demo[["object"]]$rand_type2), 0.1917001, tolerance = 1e-1) - expect_equal(mean(modscore.demo$rand_type3), + expect_equal(mean(modscore.demo[["object"]]$rand_type3), 0.1798533, tolerance = 1e-1) expected_elements <- c("MS_Celltype",pbmc.single$general.class) - expect(all(expected_elements %in% colnames(modscore.demo@meta.data)), + expect(all(expected_elements %in% colnames(modscore.demo[["object"]]@meta.data)), failure_message = "modscore results not found") }) @@ -126,33 +126,33 @@ test_that("modScore returns metadata with scores and cell calls for # skip_on_ci() # expect_snapshot_file( - # .drawMSfig(modscore.demo$ms.figures[[1]]), + # .drawMSfig(modscore.demo[["object"]]$ms.figures[[1]]), # "nsclc_multi_MS_1.png" # ) # # skip_on_ci() # expect_snapshot_file( - # .drawMSfig(modscore.demo$ms.figures[[2]]), + # .drawMSfig(modscore.demo[["object"]]$ms.figures[[2]]), # "nsclc_multi_MS_2.png" # ) # # skip_on_ci() # expect_snapshot_file( - # .drawMSfig(modscore.demo$ms.figures[[3]]), + # .drawMSfig(modscore.demo[["object"]]$ms.figures[[3]]), # "nsclc_multi_MS_3.png" # ) - expect_equal(mean(modscore.demo$rand_type1), + expect_equal(mean(modscore.demo[["object"]]$rand_type1), 0.1266389, tolerance = 1e-1) - expect_equal(mean(modscore.demo$rand_type2), + expect_equal(mean(modscore.demo[["object"]]$rand_type2), 0.1794167, tolerance = 1e-1) - expect_equal(mean(modscore.demo$rand_type3), + expect_equal(mean(modscore.demo[["object"]]$rand_type3), 0.2627949, tolerance = 1e-1) expected_elements <- c("MS_Celltype",nsclc.multi$general.class) - expect(all(expected_elements %in% colnames(modscore.demo@meta.data)), + expect(all(expected_elements %in% colnames(modscore.demo[["object"]]@meta.data)), failure_message = "modscore results not found") }) @@ -165,33 +165,33 @@ test_that("modScore returns metadata with scores and cell calls for brca", { # skip_on_ci() # expect_snapshot_file( - # .drawMSfig(modscore.demo$ms.figures[[1]]), + # .drawMSfig(modscore.demo[["object"]]$ms.figures[[1]]), # "brca_multi_MS_1.png" # ) # # skip_on_ci() # expect_snapshot_file( - # .drawMSfig(modscore.demo$ms.figures[[2]]), + # .drawMSfig(modscore.demo[["object"]]$ms.figures[[2]]), # "brca_multi_MS_2.png" # ) # # skip_on_ci() # expect_snapshot_file( - # .drawMSfig(modscore.demo$ms.figures[[3]]), + # .drawMSfig(modscore.demo[["object"]]$ms.figures[[3]]), # "brca_multi_MS_3.png" # ) - expect_equal(mean(modscore.demo$rand_type1), + expect_equal(mean(modscore.demo[["object"]]$rand_type1), 0.1958383, tolerance = 1e-1) - expect_equal(mean(modscore.demo$rand_type2), + expect_equal(mean(modscore.demo[["object"]]$rand_type2), 0.1079014, tolerance = 1e-1) - expect_equal(mean(modscore.demo$rand_type3), + expect_equal(mean(modscore.demo[["object"]]$rand_type3), 0.03652896, tolerance = 1e-1) expected_elements <- c("MS_Celltype",brca$general.class) - expect(all(expected_elements %in% colnames(modscore.demo@meta.data)), + expect(all(expected_elements %in% colnames(modscore.demo[["object"]]@meta.data)), failure_message = "modscore results not found") }) diff --git a/tests/testthat/test-Violin_Plots_by_Metadata.R b/tests/testthat/test-Violin_Plots_by_Metadata.R index 669a40e..929a41a 100755 --- a/tests/testthat/test-Violin_Plots_by_Metadata.R +++ b/tests/testthat/test-Violin_Plots_by_Metadata.R @@ -9,8 +9,7 @@ test_that("Violin plot works for TEC data", { "tec_violin.png" ) - expected_elements = c("gg", "ggplot") - expect_setequal(class(violin_test), expected_elements) + expect_setequal(any(grepl("ggplot", class(violin_test))), TRUE) }) @@ -25,8 +24,8 @@ test_that("Violin plot works for Chariou data", { "chariou_violin.png" ) - expected_elements = c("gg", "ggplot") - expect_setequal(class(violin_test), expected_elements) + + expect_setequal(any(grepl("ggplot", class(violin_test))), TRUE) }) @@ -41,8 +40,8 @@ test_that("Violin plot works for Chariou data", { # "chariou_allgroup_violin.png" # ) # -# expected_elements = c("gg", "ggplot") -# expect_setequal(class(violin_test), expected_elements) +# +# expect_setequal(any(grepl("ggplot", class(violin_test))), TRUE) # # }) # @@ -57,8 +56,8 @@ test_that("Violin plot works for Chariou data", { # "chariou_subgroup_violin.png" # ) # -# expected_elements = c("gg", "ggplot") -# expect_setequal(class(violin_test), expected_elements) +# +# expect_setequal(any(grepl("ggplot", class(violin_test))), TRUE) # # }) @@ -73,8 +72,8 @@ test_that("Violin plot works for pbmc.single data", { "pbmc_single_violin.png" ) - expected_elements = c("gg", "ggplot") - expect_setequal(class(violin_test), expected_elements) + + expect_setequal(any(grepl("ggplot", class(violin_test))), TRUE) }) @@ -89,8 +88,8 @@ test_that("Violin plot works for nsclc.multi data", { "nsclc_multi_violin.png" ) - expected_elements = c("gg", "ggplot") - expect_setequal(class(violin_test), expected_elements) + + expect_setequal(any(grepl("ggplot", class(violin_test))), TRUE) }) @@ -105,8 +104,8 @@ test_that("Violin plot works for brca data", { "brca_violin.png" ) - expected_elements = c("gg", "ggplot") - expect_setequal(class(violin_test), expected_elements) + + expect_setequal(any(grepl("ggplot", class(violin_test))), TRUE) })