From 29371be8e15a4975892879932fe6a567166ba0e3 Mon Sep 17 00:00:00 2001 From: zhanghao-njmu <542370159@qq.com> Date: Sun, 19 Nov 2023 23:21:39 +0800 Subject: [PATCH] The calculation strategy for p.adjust in RunDEtest has been modified so that p.adjust is conducted within each group of comparisons instead of globally. --- R/SCP-analysis.R | 24 +++++++--------- R/SCP-app.R | 8 ++++-- R/SCP-cellqc.R | 17 ++++++----- R/SCP-plot.R | 65 ++++++++++++++++++++++++------------------ R/Seurat-function.R | 18 ++++++------ man/CellCorHeatmap.Rd | 8 +++--- man/DynamicHeatmap.Rd | 8 +++--- man/FeatureStatPlot.Rd | 8 ++++++ man/GroupHeatmap.Rd | 8 +++--- man/RunCellQC.Rd | 10 +++---- man/RunUMAP2.Rd | 2 ++ 11 files changed, 99 insertions(+), 77 deletions(-) diff --git a/R/SCP-analysis.R b/R/SCP-analysis.R index a8b3df7e..b0c1bc87 100644 --- a/R/SCP-analysis.R +++ b/R/SCP-analysis.R @@ -1925,6 +1925,9 @@ RunDEtest <- function(srt, group_by = NULL, group1 = NULL, group2 = NULL, cells1 markers[, "gene"] <- rownames(markers) markers[, "group1"] <- as.character(group) markers[, "group2"] <- "others" + if ("p_val" %in% colnames(markers)) { + markers[, "p_val_adj"] <- p.adjust(markers[, "p_val"], method = p.adjust.method) + } return(markers) } else { return(NULL) @@ -1935,9 +1938,6 @@ RunDEtest <- function(srt, group_by = NULL, group1 = NULL, group2 = NULL, cells1 if (!is.null(AllMarkers) && nrow(AllMarkers) > 0) { rownames(AllMarkers) <- NULL AllMarkers[, "group1"] <- factor(AllMarkers[, "group1"], levels = levels(cell_group)) - if ("p_val" %in% colnames(AllMarkers)) { - AllMarkers[, "p_val_adj"] <- p.adjust(AllMarkers[, "p_val"], method = p.adjust.method) - } AllMarkers[, "test_group_number"] <- as.integer(table(AllMarkers[["gene"]])[AllMarkers[, "gene"]]) AllMarkersMatrix <- as.data.frame.matrix(table(AllMarkers[, c("gene", "group1")])) AllMarkers[, "test_group"] <- apply(AllMarkersMatrix, 1, function(x) { @@ -1965,6 +1965,9 @@ RunDEtest <- function(srt, group_by = NULL, group1 = NULL, group2 = NULL, cells1 markers[, "gene"] <- rownames(markers) markers[, "group1"] <- as.character(pair[i, 1]) markers[, "group2"] <- as.character(pair[i, 2]) + if ("p_val" %in% colnames(markers)) { + markers[, "p_val_adj"] <- p.adjust(markers[, "p_val"], method = p.adjust.method) + } return(markers) } else { return(NULL) @@ -1975,9 +1978,6 @@ RunDEtest <- function(srt, group_by = NULL, group1 = NULL, group2 = NULL, cells1 if (!is.null(PairedMarkers) && nrow(PairedMarkers) > 0) { rownames(PairedMarkers) <- NULL PairedMarkers[, "group1"] <- factor(PairedMarkers[, "group1"], levels = levels(cell_group)) - if ("p_val" %in% colnames(PairedMarkers)) { - PairedMarkers[, "p_val_adj"] <- p.adjust(PairedMarkers[, "p_val"], method = p.adjust.method) - } PairedMarkers[, "test_group_number"] <- as.integer(table(PairedMarkers[["gene"]])[PairedMarkers[, "gene"]]) PairedMarkersMatrix <- as.data.frame.matrix(table(PairedMarkers[, c("gene", "group1")])) PairedMarkers[, "test_group"] <- apply(PairedMarkersMatrix, 1, function(x) { @@ -2011,19 +2011,19 @@ RunDEtest <- function(srt, group_by = NULL, group1 = NULL, group2 = NULL, cells1 markers[, "gene"] <- rownames(markers) markers[, "group1"] <- as.character(group) markers[, "group2"] <- "others" + if ("p_val" %in% colnames(markers)) { + markers[, "p_val_adj"] <- p.adjust(markers[, "p_val"], method = p.adjust.method) + } return(markers) } else { return(NULL) } } }, BPPARAM = BPPARAM) - ConservedMarkers <- do.call(rbind.data.frame, lapply(ConservedMarkers, function(x) x[, c("avg_log2FC", "pct.1", "pct.2", "max_pval", "p_val", "gene", "group1", "group2")])) + ConservedMarkers <- do.call(rbind.data.frame, lapply(ConservedMarkers, function(x) x[, c("avg_log2FC", "pct.1", "pct.2", "max_pval", "p_val", "p_val_adj", "gene", "group1", "group2")])) if (!is.null(ConservedMarkers) && nrow(ConservedMarkers) > 0) { rownames(ConservedMarkers) <- NULL ConservedMarkers[, "group1"] <- factor(ConservedMarkers[, "group1"], levels = levels(cell_group)) - if ("p_val" %in% colnames(ConservedMarkers)) { - ConservedMarkers[, "p_val_adj"] <- p.adjust(ConservedMarkers[, "p_val"], method = p.adjust.method) - } ConservedMarkers[, "test_group_number"] <- as.integer(table(ConservedMarkers[["gene"]])[ConservedMarkers[, "gene"]]) ConservedMarkersMatrix <- as.data.frame.matrix(table(ConservedMarkers[, c("gene", "group1")])) ConservedMarkers[, "test_group"] <- apply(ConservedMarkersMatrix, 1, function(x) { @@ -2085,9 +2085,6 @@ RunDEtest <- function(srt, group_by = NULL, group1 = NULL, group2 = NULL, cells1 if (!is.null(DisturbedMarkers) && nrow(DisturbedMarkers) > 0) { rownames(DisturbedMarkers) <- NULL DisturbedMarkers[, "group1"] <- factor(DisturbedMarkers[, "group1"], levels = levels(cell_group)) - if ("p_val" %in% colnames(DisturbedMarkers)) { - DisturbedMarkers[, "p_val_adj"] <- p.adjust(DisturbedMarkers[, "p_val"], method = p.adjust.method) - } DisturbedMarkers[, "test_group_number"] <- as.integer(table(unique(DisturbedMarkers[, c("gene", "group1")])[["gene"]])[DisturbedMarkers[, "gene"]]) DisturbedMarkersMatrix <- as.data.frame.matrix(table(DisturbedMarkers[, c("gene", "group1")])) DisturbedMarkers[, "test_group"] <- apply(DisturbedMarkersMatrix, 1, function(x) { @@ -5398,6 +5395,7 @@ srt_to_adata <- function(srt, features = NULL, var[["highly_variable"]] <- features %in% VariableFeatures(srt, assay = assay_X) } + # X <- t(as_matrix(slot(srt@assays[[assay_X]], slot_X)[features, , drop = FALSE])) X <- t(GetAssayData(srt, assay = assay_X, slot = slot_X)[features, , drop = FALSE]) adata <- sc$AnnData( X = np_array(X, dtype = np$float32), diff --git a/R/SCP-app.R b/R/SCP-app.R index d2003679..af18df30 100644 --- a/R/SCP-app.R +++ b/R/SCP-app.R @@ -145,7 +145,7 @@ CreateMetaFile <- function(srt, MetaFile, name = NULL, write_tools = FALSE, writ meta_asfeatures <- c(meta_asfeatures, var) } else { if (length(unique(meta)) > ignore_nlevel) { - warning("The number of categories in ", var, " is greater than ", ignore_nlevel, ", it will be ignored.", immediate. = TRUE) + warning("The number of categories in ", var, " is greater than ", ignore_nlevel, ". It will be ignored.", immediate. = TRUE) } else { meta_asgroups <- c(meta_asgroups, var) } @@ -156,7 +156,11 @@ CreateMetaFile <- function(srt, MetaFile, name = NULL, write_tools = FALSE, writ if (paste0(name, "/metadata/", var) %in% paste0(h5ls(MetaFile)$group, "/", h5ls(MetaFile)$name)) { message("Group ", paste0(name, "/metadata/", var), " already exists in the ", MetaFile) } else { - h5write(obj = meta, file = MetaFile, name = paste0(name, "/metadata/", var), write.attributes = write.attributes, level = compression_level) + if (all(is.na(meta))) { + warning("All of values in ", var, " is NA. It will be ignored.", immediate. = TRUE) + } else { + h5write(obj = meta, file = MetaFile, name = paste0(name, "/metadata/", var), write.attributes = write.attributes, level = compression_level) + } } } diff --git a/R/SCP-cellqc.R b/R/SCP-cellqc.R index c716d3b8..7ceb2b6b 100644 --- a/R/SCP-cellqc.R +++ b/R/SCP-cellqc.R @@ -259,10 +259,10 @@ isOutlier <- function(x, nmads = 2.5, constant = 1.4826, type = c("both", "lower #' #' @inheritParams RunDoubletCalling #' @param split.by Name of the sample variable to split the Seurat object. Default is NULL. +#' @param return_filtered Logical indicating whether to return a cell-filtered Seurat object. Default is FALSE. #' @param qc_metrics A character vector specifying the quality control metrics to be applied. Default is #' `c("doublets", "outlier", "umi", "gene", "mito", "ribo", "ribo_mito_ratio", "species")`. -#' @param return_filtered Logical indicating whether to return a cell-filtered Seurat object. Default is FALSE. -#' @param outlier_cutoff A character vector specifying the outlier cutoff values. Default is +#' @param outlier_threshold A character vector specifying the outlier threshold. Default is #' `c("log10_nCount:lower:2.5", "log10_nCount:higher:5", "log10_nFeature:lower:2.5", "log10_nFeature:higher:5", "featurecount_dist:lower:2.5")`. See \link[scuttle]{isOutlier}. #' @param db_coefficient The coefficient used to calculate the doublet rate. Default is 0.01. Doublet rate is calculated as`ncol(srt) / 1000 * db_coefficient` #' @param outlier_n Minimum number of outlier metrics that meet the conditions for determining outlier cells. Default is 1. @@ -311,11 +311,10 @@ isOutlier <- function(x, nmads = 2.5, constant = 1.4826, type = c("both", "lower #' @importFrom Matrix colSums t #' @export #' -RunCellQC <- function(srt, assay = "RNA", split.by = NULL, +RunCellQC <- function(srt, assay = "RNA", split.by = NULL, return_filtered = FALSE, qc_metrics = c("doublets", "outlier", "umi", "gene", "mito", "ribo", "ribo_mito_ratio", "species"), - return_filtered = FALSE, db_method = "scDblFinder", db_rate = NULL, db_coefficient = 0.01, - outlier_cutoff = c( + outlier_threshold = c( "log10_nCount:lower:2.5", "log10_nCount:higher:5", "log10_nFeature:lower:2.5", @@ -344,7 +343,7 @@ RunCellQC <- function(srt, assay = "RNA", split.by = NULL, } status <- check_DataType(srt, slot = "counts", assay = assay) if (status != "raw_counts") { - warning("Data type is not raw counts!") + warning("Data type is not raw counts!", immediate. = TRUE) } if (!paste0("nCount_", assay) %in% colnames(srt@meta.data)) { srt@meta.data[[paste0("nCount_", assay)]] <- colSums(srt[[assay]]@counts) @@ -429,15 +428,15 @@ RunCellQC <- function(srt, assay = "RNA", split.by = NULL, # theme(panel.background = element_rect(fill = "grey")) # nrow(lower_df) - var <- sapply(strsplit(outlier_cutoff, ":"), function(x) x[[1]]) + var <- sapply(strsplit(outlier_threshold, ":"), function(x) x[[1]]) var_valid <- var %in% colnames(srt@meta.data) | sapply(var, FUN = function(x) exists(x, where = environment())) if (any(!var_valid)) { stop("Variable ", paste0(names(var_valid)[!var_valid], collapse = ","), " is not found in the srt object.") } - outlier <- lapply(strsplit(outlier_cutoff, ":"), function(m) { + outlier <- lapply(strsplit(outlier_threshold, ":"), function(m) { colnames(srt)[isOutlier(get(m[1]), nmads = as.numeric(m[3]), type = m[2])] }) - names(outlier) <- outlier_cutoff + names(outlier) <- outlier_threshold # print(unlist(lapply(outlier, length))) outlier_tb <- table(unlist(outlier)) outlier_qc <- c(outlier_qc, names(outlier_tb)[outlier_tb >= outlier_n]) diff --git a/R/SCP-plot.R b/R/SCP-plot.R index 602d2a61..4d003fc2 100644 --- a/R/SCP-plot.R +++ b/R/SCP-plot.R @@ -1590,11 +1590,11 @@ CellDimPlot <- function(srt, group.by, reduction = NULL, dims = c(1, 2), split.b p <- p + scattermore::geom_scattermore( data = dat[is.na(dat[, "group.by"]), , drop = FALSE], mapping = aes(x = .data[["x"]], y = .data[["y"]]), color = bg_color, - pointsize = floor(pt.size), alpha = pt.alpha, pixels = raster.dpi + pointsize = ceiling(pt.size), alpha = pt.alpha, pixels = raster.dpi ) + scattermore::geom_scattermore( data = dat[!is.na(dat[, "group.by"]), , drop = FALSE], mapping = aes(x = .data[["x"]], y = .data[["y"]], color = .data[["group.by"]]), - pointsize = floor(pt.size), alpha = pt.alpha, pixels = raster.dpi + pointsize = ceiling(pt.size), alpha = pt.alpha, pixels = raster.dpi ) } else if (isTRUE(hex)) { check_R("hexbin") @@ -2325,12 +2325,12 @@ FeatureDimPlot <- function(srt, features, reduction = NULL, dims = c(1, 2), spli p <- p + scattermore::geom_scattermore( data = dat[dat[, "color_blend"] == bg_color, , drop = FALSE], mapping = aes(x = .data[["x"]], y = .data[["y"]], color = .data[["color_blend"]]), - pointsize = floor(pt.size), alpha = pt.alpha, pixels = raster.dpi + pointsize = ceiling(pt.size), alpha = pt.alpha, pixels = raster.dpi ) + scattermore::geom_scattermore( data = dat[dat[, "color_blend"] != bg_color, , drop = FALSE], mapping = aes(x = .data[["x"]], y = .data[["y"]], color = .data[["color_blend"]]), - pointsize = floor(pt.size), alpha = pt.alpha, pixels = raster.dpi + pointsize = ceiling(pt.size), alpha = pt.alpha, pixels = raster.dpi ) + scale_color_identity() + new_scale_color() @@ -2599,12 +2599,12 @@ FeatureDimPlot <- function(srt, features, reduction = NULL, dims = c(1, 2), spli p <- p + scattermore::geom_scattermore( data = dat[is.na(dat[, "value"]), , drop = FALSE], mapping = aes(x = .data[["x"]], y = .data[["y"]], color = .data[["value"]]), - pointsize = floor(pt.size), alpha = pt.alpha, pixels = raster.dpi + pointsize = ceiling(pt.size), alpha = pt.alpha, pixels = raster.dpi ) + scattermore::geom_scattermore( data = dat[!is.na(dat[, "value"]), , drop = FALSE], mapping = aes(x = .data[["x"]], y = .data[["y"]], color = .data[["value"]]), - pointsize = floor(pt.size), alpha = pt.alpha, pixels = raster.dpi + pointsize = ceiling(pt.size), alpha = pt.alpha, pixels = raster.dpi ) } else if (isTRUE(hex)) { check_R("hexbin") @@ -3399,6 +3399,8 @@ FeatureDimPlot3D <- function(srt, features, reduction = NULL, dims = c(1, 2, 3), #' FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "SubCellType", add_box = TRUE) #' FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "SubCellType", add_point = TRUE) #' FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "SubCellType", add_trend = TRUE) +#' FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "SubCellType", add_stat = "mean") +#' FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "SubCellType", add_line = 0.2, line_type = 2) #' FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "SubCellType", split.by = "Phase") #' FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "SubCellType", split.by = "Phase", add_box = TRUE, add_trend = TRUE) #' FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "SubCellType", split.by = "Phase", comparisons = TRUE) @@ -3471,7 +3473,8 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b add_box = FALSE, box_color = "black", box_width = 0.1, box_ptsize = 2, add_point = FALSE, pt.color = "grey30", pt.size = NULL, pt.alpha = 1, jitter.width = 0.4, jitter.height = 0.1, add_trend = FALSE, trend_color = "black", trend_linewidth = 1, trend_ptsize = 2, - add_stat = c("none", "mean", "median"), stat_color = "black", stat_size = 1, + add_stat = c("none", "mean", "median"), stat_color = "black", stat_size = 1, stat_stroke = 1, stat_shape = 25, + add_line = NULL, line_color = "red", line_size = 1, line_type = 1, cells.highlight = NULL, cols.highlight = "red", sizes.highlight = 1, alpha.highlight = 1, calculate_coexp = FALSE, same.y.lims = FALSE, y.min = NULL, y.max = NULL, y.trans = "identity", y.nbreaks = 5, @@ -3524,7 +3527,8 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b add_box = add_box, box_color = box_color, box_width = box_width, box_ptsize = box_ptsize, add_point = add_point, pt.color = pt.color, pt.size = pt.size, pt.alpha = pt.alpha, jitter.width = jitter.width, jitter.height = jitter.height, add_trend = add_trend, trend_color = trend_color, trend_linewidth = trend_linewidth, trend_ptsize = trend_ptsize, - add_stat = add_stat, stat_color = stat_color, stat_size = stat_size, + add_stat = add_stat, stat_color = stat_color, stat_size = stat_size, stat_stroke = stat_stroke, stat_shape = stat_shape, + add_line = add_line, line_color = line_color, line_size = line_size, line_type = line_type, cells.highlight = cells.highlight, cols.highlight = cols.highlight, sizes.highlight = sizes.highlight, alpha.highlight = alpha.highlight, calculate_coexp = calculate_coexp, same.y.lims = same.y.lims, y.min = y.min, y.max = y.max, y.trans = y.trans, y.nbreaks = y.nbreaks, @@ -3551,7 +3555,8 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b add_box = add_box, box_color = box_color, box_width = box_width, box_ptsize = box_ptsize, add_point = add_point, pt.color = pt.color, pt.size = pt.size, pt.alpha = pt.alpha, jitter.width = jitter.width, jitter.height = jitter.height, add_trend = add_trend, trend_color = trend_color, trend_linewidth = trend_linewidth, trend_ptsize = trend_ptsize, - add_stat = add_stat, stat_color = stat_color, stat_size = stat_size, + add_stat = add_stat, stat_color = stat_color, stat_size = stat_size, stat_stroke = stat_stroke, stat_shape = stat_shape, + add_line = add_line, line_color = line_color, line_size = line_size, line_type = line_type, cells.highlight = cells.highlight, cols.highlight = cols.highlight, sizes.highlight = sizes.highlight, alpha.highlight = alpha.highlight, calculate_coexp = calculate_coexp, same.y.lims = same.y.lims, y.min = y.min, y.max = y.max, y.trans = y.trans, y.nbreaks = y.nbreaks, @@ -3667,7 +3672,8 @@ ExpressionStatPlot <- function(exp.data, meta.data, stat.by, group.by = NULL, sp add_box = FALSE, box_color = "black", box_width = 0.1, box_ptsize = 2, add_point = FALSE, pt.color = "grey30", pt.size = NULL, pt.alpha = 1, jitter.width = 0.4, jitter.height = 0.1, add_trend = FALSE, trend_color = "black", trend_linewidth = 1, trend_ptsize = 2, - add_stat = c("none", "mean", "median"), stat_color = "black", stat_size = 1, + add_stat = c("none", "mean", "median"), stat_color = "black", stat_size = 1, stat_stroke = 1, stat_shape = 25, + add_line = NULL, line_color = "red", line_size = 1, line_type = 1, cells.highlight = NULL, cols.highlight = "red", sizes.highlight = 1, alpha.highlight = 1, calculate_coexp = FALSE, same.y.lims = FALSE, y.min = NULL, y.max = NULL, y.trans = "identity", y.nbreaks = 5, @@ -3686,6 +3692,9 @@ ExpressionStatPlot <- function(exp.data, meta.data, stat.by, group.by = NULL, sp fill.by <- match.arg(fill.by) sig_label <- match.arg(sig_label) add_stat <- match.arg(add_stat) + if (!is.null(add_line)) { + stopifnot(is.numeric(add_line)) + } if (missing(exp.data)) { exp.data <- matrix(0, nrow = 1, ncol = nrow(meta.data), dimnames = list("", rownames(meta.data))) @@ -4245,10 +4254,16 @@ ExpressionStatPlot <- function(exp.data, meta.data, stat.by, group.by = NULL, sp } if (add_stat != "none") { p <- p + stat_summary( - fun = add_stat, geom = "point", mapping = aes(group = .data[["split.by"]], shape = 95), - position = position_dodge(width = 0.9), color = stat_color, fill = "white", size = stat_size, stroke = 10, + fun = add_stat, geom = "point", mapping = aes(group = .data[["split.by"]], shape = stat_shape), + position = position_dodge(width = 0.9), color = stat_color, fill = stat_color, size = stat_size, stroke = stat_stroke, ) + scale_shape_identity() } + if (!is.null(add_line)) { + p <- p + geom_hline( + yintercept = add_line, + color = line_color, linetype = line_type, linewidth = line_size + ) + } if (nrow(dat) == 0) { p <- p + facet_null() @@ -5356,7 +5371,7 @@ FeatureCorPlot <- function(srt, features, group.by = NULL, split.by = NULL, cell if (isTRUE(raster)) { p <- p + scattermore::geom_scattermore( mapping = aes(x = .data[[f1]], y = .data[[f2]], color = .data[[group.by]]), - pointsize = floor(pt.size), alpha = pt.alpha, pixels = raster.dpi + pointsize = ceiling(pt.size), alpha = pt.alpha, pixels = raster.dpi ) } else { p <- p + geom_point(aes(x = .data[[f1]], y = .data[[f2]], color = .data[[group.by]]), @@ -7786,7 +7801,7 @@ mestimate <- function(data) { #' features = de_top$gene, feature_split = de_top$group1, group.by = "CellType", #' heatmap_palette = "YlOrRd", #' cell_annotation = c("Phase", "G2M_score", "Neurod2"), cell_annotation_palette = c("Dark2", "Paired", "Paired"), -#' cell_annotation_params = list(height = grid::unit(20, "mm")), +#' cell_annotation_params = list(height = grid::unit(10, "mm")), #' feature_annotation = c("TF", "CSPA"), #' feature_annotation_palcolor = list(c("gold", "steelblue"), c("forestgreen")), #' add_dot = TRUE, add_bg = TRUE, nlabel = 0, show_row_names = TRUE @@ -7797,7 +7812,7 @@ mestimate <- function(data) { #' features = de_top$gene, feature_split = de_top$group1, group.by = "CellType", #' heatmap_palette = "YlOrRd", #' cell_annotation = c("Phase", "G2M_score", "Neurod2"), cell_annotation_palette = c("Dark2", "Paired", "Paired"), -#' cell_annotation_params = list(width = grid::unit(20, "mm")), +#' cell_annotation_params = list(width = grid::unit(10, "mm")), #' feature_annotation = c("TF", "CSPA"), #' feature_annotation_palcolor = list(c("gold", "steelblue"), c("forestgreen")), #' add_dot = TRUE, add_bg = TRUE, @@ -7862,7 +7877,7 @@ GroupHeatmap <- function(srt, features = NULL, group.by = NULL, split.by = NULL, add_violin = FALSE, fill.by = "feature", fill_palette = "Dark2", fill_palcolor = NULL, heatmap_palette = "RdBu", heatmap_palcolor = NULL, group_palette = "Paired", group_palcolor = NULL, cell_split_palette = "simspec", cell_split_palcolor = NULL, feature_split_palette = "simspec", feature_split_palcolor = NULL, - cell_annotation = NULL, cell_annotation_palette = "Paired", cell_annotation_palcolor = NULL, cell_annotation_params = if (flip) list(width = grid::unit(20, "mm")) else list(height = grid::unit(20, "mm")), + cell_annotation = NULL, cell_annotation_palette = "Paired", cell_annotation_palcolor = NULL, cell_annotation_params = if (flip) list(width = grid::unit(10, "mm")) else list(height = grid::unit(10, "mm")), feature_annotation = NULL, feature_annotation_palette = "Dark2", feature_annotation_palcolor = NULL, feature_annotation_params = if (flip) list(height = grid::unit(5, "mm")) else list(width = grid::unit(5, "mm")), use_raster = NULL, raster_device = "png", raster_by_magick = FALSE, height = NULL, width = NULL, units = "inch", seed = 11, ht_params = list()) { @@ -8878,9 +8893,7 @@ GroupHeatmap <- function(srt, features = NULL, group.by = NULL, split.by = NULL, if ((!is.null(row_split) && length(index) > 0) || any(c(anno_terms, anno_keys, anno_features)) || !is.null(width) || !is.null(height)) { fix <- TRUE if (is.null(width) || is.null(height)) { - message("The size of the heatmap is fixed because certain elements are not scalable.\n - The width and height of the heatmap are determined by the size of the current viewport.\n - If you want to have more control over the size, you can manually set the parameters 'width' and 'height'.") + message("\nThe size of the heatmap is fixed because certain elements are not scalable.\nThe width and height of the heatmap are determined by the size of the current viewport.\nIf you want to have more control over the size, you can manually set the parameters 'width' and 'height'.\n") } } else { fix <- FALSE @@ -9822,9 +9835,7 @@ FeatureHeatmap <- function(srt, features = NULL, cells = NULL, group.by = NULL, if ((!is.null(row_split) && length(index) > 0) || any(c(anno_terms, anno_keys, anno_features)) || !is.null(width) || !is.null(height)) { fix <- TRUE if (is.null(width) || is.null(height)) { - message("The size of the heatmap is fixed because certain elements are not scalable.\n - The width and height of the heatmap are determined by the size of the current viewport.\n - If you want to have more control over the size, you can manually set the parameters 'width' and 'height'.") + message("\nThe size of the heatmap is fixed because certain elements are not scalable.\nThe width and height of the heatmap are determined by the size of the current viewport.\nIf you want to have more control over the size, you can manually set the parameters 'width' and 'height'.\n") } } else { fix <- FALSE @@ -10064,8 +10075,8 @@ CellCorHeatmap <- function(srt_query, srt_ref = NULL, bulk_ref = NULL, heatmap_palette = "RdBu", heatmap_palcolor = NULL, query_group_palette = "Paired", query_group_palcolor = NULL, ref_group_palette = "simspec", ref_group_palcolor = NULL, - query_cell_annotation = NULL, query_cell_annotation_palette = "Paired", query_cell_annotation_palcolor = NULL, query_cell_annotation_params = if (flip) list(height = grid::unit(20, "mm")) else list(width = grid::unit(20, "mm")), - ref_cell_annotation = NULL, ref_cell_annotation_palette = "Paired", ref_cell_annotation_palcolor = NULL, ref_cell_annotation_params = if (flip) list(width = grid::unit(20, "mm")) else list(height = grid::unit(20, "mm")), + query_cell_annotation = NULL, query_cell_annotation_palette = "Paired", query_cell_annotation_palcolor = NULL, query_cell_annotation_params = if (flip) list(height = grid::unit(10, "mm")) else list(width = grid::unit(10, "mm")), + ref_cell_annotation = NULL, ref_cell_annotation_palette = "Paired", ref_cell_annotation_palcolor = NULL, ref_cell_annotation_params = if (flip) list(width = grid::unit(10, "mm")) else list(height = grid::unit(10, "mm")), use_raster = NULL, raster_device = "png", raster_by_magick = FALSE, height = NULL, width = NULL, units = "inch", seed = 11, ht_params = list()) { set.seed(seed) @@ -10892,7 +10903,7 @@ CellCorHeatmap <- function(srt_query, srt_ref = NULL, bulk_ref = NULL, #' cell_annotation_palette = c("Paired", "simspec", "Purples"), #' separate_annotation = list("SubCellType", c("Nnat", "Irx1")), #' separate_annotation_palette = c("Paired", "Set1"), -#' separate_annotation_params = list(height = grid::unit(20, "mm")), +#' separate_annotation_params = list(height = grid::unit(10, "mm")), #' feature_annotation = c("TF", "CSPA"), #' feature_annotation_palcolor = list(c("gold", "steelblue"), c("forestgreen")), #' pseudotime_label = 25, @@ -10912,7 +10923,7 @@ CellCorHeatmap <- function(srt_query, srt_ref = NULL, bulk_ref = NULL, #' cell_annotation_palette = c("Paired", "simspec", "Purples"), #' separate_annotation = list("SubCellType", c("Nnat", "Irx1")), #' separate_annotation_palette = c("Paired", "Set1"), -#' separate_annotation_params = list(width = grid::unit(20, "mm")), +#' separate_annotation_params = list(width = grid::unit(10, "mm")), #' feature_annotation = c("TF", "CSPA"), #' feature_annotation_palcolor = list(c("gold", "steelblue"), c("forestgreen")), #' pseudotime_label = 25, @@ -10969,7 +10980,7 @@ DynamicHeatmap <- function(srt, lineages, features = NULL, use_fitted = FALSE, b feature_split_palette = "simspec", feature_split_palcolor = NULL, cell_annotation = NULL, cell_annotation_palette = "Paired", cell_annotation_palcolor = NULL, cell_annotation_params = if (flip) list(width = grid::unit(5, "mm")) else list(height = grid::unit(5, "mm")), feature_annotation = NULL, feature_annotation_palette = "Dark2", feature_annotation_palcolor = NULL, feature_annotation_params = if (flip) list(height = grid::unit(5, "mm")) else list(width = grid::unit(5, "mm")), - separate_annotation = NULL, separate_annotation_palette = "Paired", separate_annotation_palcolor = NULL, separate_annotation_params = if (flip) list(width = grid::unit(20, "mm")) else list(height = grid::unit(20, "mm")), + separate_annotation = NULL, separate_annotation_palette = "Paired", separate_annotation_palcolor = NULL, separate_annotation_params = if (flip) list(width = grid::unit(10, "mm")) else list(height = grid::unit(10, "mm")), reverse_ht = NULL, use_raster = NULL, raster_device = "png", raster_by_magick = FALSE, height = NULL, width = NULL, units = "inch", seed = 11, ht_params = list()) { set.seed(seed) diff --git a/R/Seurat-function.R b/R/Seurat-function.R index 460b7788..b984d1f9 100644 --- a/R/Seurat-function.R +++ b/R/Seurat-function.R @@ -587,7 +587,7 @@ RunUMAP2 <- function(object, ...) { RunUMAP2.Seurat <- function(object, reduction = "pca", dims = NULL, features = NULL, neighbor = NULL, graph = NULL, assay = NULL, slot = "data", - umap.method = "uwot", reduction.model = NULL, + umap.method = "uwot", reduction.model = NULL, n_threads = NULL, return.model = FALSE, n.neighbors = 30L, n.components = 2L, metric = "cosine", n.epochs = 200L, spread = 1, min.dist = 0.3, set.op.mix.ratio = 1, local.connectivity = 1L, negative.sample.rate = 5L, @@ -667,7 +667,7 @@ RunUMAP2.Seurat <- function(object, #' @importFrom Matrix sparseMatrix #' @export RunUMAP2.default <- function(object, assay = NULL, - umap.method = "uwot", reduction.model = NULL, + umap.method = "uwot", reduction.model = NULL, n_threads = NULL, return.model = FALSE, n.neighbors = 30L, n.components = 2L, metric = "cosine", n.epochs = 200L, spread = 1, min.dist = 0.3, set.op.mix.ratio = 1, local.connectivity = 1L, negative.sample.rate = 5L, @@ -879,7 +879,7 @@ RunUMAP2.default <- function(object, assay = NULL, if (umap.method == "uwot") { if (inherits(x = object, what = "dist")) { embeddings <- uwot::umap( - X = object, n_neighbors = n.neighbors, n_threads = 1, n_components = n.components, + X = object, n_neighbors = n.neighbors, n_threads = n_threads, n_components = n.components, metric = metric, n_epochs = n.epochs, learning_rate = learning.rate, min_dist = min.dist, spread = spread, set_op_mix_ratio = set.op.mix.ratio, local_connectivity = local.connectivity, repulsion_strength = repulsion.strength, @@ -904,7 +904,7 @@ RunUMAP2.default <- function(object, assay = NULL, # object <- srt@neighbors$Seurat.nn # object <- list(idx = Indices(object), dist = Distances(object)) out <- uwot::umap( - X = NULL, nn_method = object, n_threads = 1, n_components = n.components, + X = NULL, nn_method = object, n_threads = n_threads, n_components = n.components, metric = metric, n_epochs = n.epochs, learning_rate = learning.rate, min_dist = min.dist, spread = spread, set_op_mix_ratio = set.op.mix.ratio, local_connectivity = local.connectivity, repulsion_strength = repulsion.strength, @@ -962,7 +962,7 @@ RunUMAP2.default <- function(object, assay = NULL, # idx <- t(as_matrix(apply(object, 2, function(x) order(x, decreasing = TRUE)[1:n.neighbors]))) # connectivity <- t(as_matrix(apply(object, 2, function(x) x[order(x, decreasing = TRUE)[1:n.neighbors]]))) out <- uwot::umap( - X = NULL, nn_method = nn, n_threads = 1, n_components = n.components, + X = NULL, nn_method = nn, n_threads = n_threads, n_components = n.components, metric = metric, n_epochs = n.epochs, learning_rate = learning.rate, min_dist = min.dist, spread = spread, set_op_mix_ratio = set.op.mix.ratio, local_connectivity = local.connectivity, repulsion_strength = repulsion.strength, @@ -1004,7 +1004,7 @@ RunUMAP2.default <- function(object, assay = NULL, if (inherits(x = object, what = "matrix") || inherits(x = object, what = "Matrix")) { # object <- srt@reductions$Seuratpca@cell.embeddings out <- uwot::umap( - X = object, n_neighbors = n.neighbors, n_threads = 1, n_components = n.components, + X = object, n_neighbors = n.neighbors, n_threads = n_threads, n_components = n.components, metric = metric, n_epochs = n.epochs, learning_rate = learning.rate, min_dist = min.dist, spread = spread, set_op_mix_ratio = set.op.mix.ratio, local_connectivity = local.connectivity, repulsion_strength = repulsion.strength, @@ -1067,7 +1067,7 @@ RunUMAP2.default <- function(object, assay = NULL, } embeddings <- uwot::umap_transform( X = NULL, nn_method = object, model = model, n_epochs = n.epochs, - n_threads = 1, verbose = verbose + n_threads = n_threads, verbose = verbose ) rownames(x = embeddings) <- row.names(object[["idx"]]) colnames(x = embeddings) <- paste0(reduction.key, 1:n.components) @@ -1088,7 +1088,7 @@ RunUMAP2.default <- function(object, assay = NULL, } embeddings <- uwot::umap_transform( X = NULL, nn_method = object, model = model, n_epochs = n.epochs, - n_threads = 1, verbose = verbose + n_threads = n_threads, verbose = verbose ) rownames(x = embeddings) <- row.names(object[["idx"]]) colnames(x = embeddings) <- paste0(reduction.key, 1:n.components) @@ -1103,7 +1103,7 @@ RunUMAP2.default <- function(object, assay = NULL, if (inherits(x = object, what = "matrix") || inherits(x = object, what = "Matrix")) { embeddings <- uwot::umap_transform( X = object, model = model, n_epochs = n.epochs, - n_threads = 1, verbose = verbose + n_threads = n_threads, verbose = verbose ) rownames(x = embeddings) <- row.names(object) colnames(x = embeddings) <- paste0(reduction.key, 1:n.components) diff --git a/man/CellCorHeatmap.Rd b/man/CellCorHeatmap.Rd index 9c12ba94..85ca5dbd 100644 --- a/man/CellCorHeatmap.Rd +++ b/man/CellCorHeatmap.Rd @@ -59,13 +59,13 @@ CellCorHeatmap( query_cell_annotation = NULL, query_cell_annotation_palette = "Paired", query_cell_annotation_palcolor = NULL, - query_cell_annotation_params = if (flip) list(height = grid::unit(20, "mm")) else - list(width = grid::unit(20, "mm")), + query_cell_annotation_params = if (flip) list(height = grid::unit(10, "mm")) else + list(width = grid::unit(10, "mm")), ref_cell_annotation = NULL, ref_cell_annotation_palette = "Paired", ref_cell_annotation_palcolor = NULL, - ref_cell_annotation_params = if (flip) list(width = grid::unit(20, "mm")) else - list(height = grid::unit(20, "mm")), + ref_cell_annotation_params = if (flip) list(width = grid::unit(10, "mm")) else + list(height = grid::unit(10, "mm")), use_raster = NULL, raster_device = "png", raster_by_magick = FALSE, diff --git a/man/DynamicHeatmap.Rd b/man/DynamicHeatmap.Rd index c3a248c9..57410020 100644 --- a/man/DynamicHeatmap.Rd +++ b/man/DynamicHeatmap.Rd @@ -110,8 +110,8 @@ DynamicHeatmap( separate_annotation = NULL, separate_annotation_palette = "Paired", separate_annotation_palcolor = NULL, - separate_annotation_params = if (flip) list(width = grid::unit(20, "mm")) else - list(height = grid::unit(20, "mm")), + separate_annotation_params = if (flip) list(width = grid::unit(10, "mm")) else + list(height = grid::unit(10, "mm")), reverse_ht = NULL, use_raster = NULL, raster_device = "png", @@ -399,7 +399,7 @@ ht4 <- DynamicHeatmap( cell_annotation_palette = c("Paired", "simspec", "Purples"), separate_annotation = list("SubCellType", c("Nnat", "Irx1")), separate_annotation_palette = c("Paired", "Set1"), - separate_annotation_params = list(height = grid::unit(20, "mm")), + separate_annotation_params = list(height = grid::unit(10, "mm")), feature_annotation = c("TF", "CSPA"), feature_annotation_palcolor = list(c("gold", "steelblue"), c("forestgreen")), pseudotime_label = 25, @@ -419,7 +419,7 @@ ht5 <- DynamicHeatmap( cell_annotation_palette = c("Paired", "simspec", "Purples"), separate_annotation = list("SubCellType", c("Nnat", "Irx1")), separate_annotation_palette = c("Paired", "Set1"), - separate_annotation_params = list(width = grid::unit(20, "mm")), + separate_annotation_params = list(width = grid::unit(10, "mm")), feature_annotation = c("TF", "CSPA"), feature_annotation_palcolor = list(c("gold", "steelblue"), c("forestgreen")), pseudotime_label = 25, diff --git a/man/FeatureStatPlot.Rd b/man/FeatureStatPlot.Rd index 7697c3f3..58bc795e 100644 --- a/man/FeatureStatPlot.Rd +++ b/man/FeatureStatPlot.Rd @@ -41,6 +41,12 @@ FeatureStatPlot( add_stat = c("none", "mean", "median"), stat_color = "black", stat_size = 1, + stat_stroke = 1, + stat_shape = 25, + add_line = NULL, + line_color = "red", + line_size = 1, + line_type = 1, cells.highlight = NULL, cols.highlight = "red", sizes.highlight = 1, @@ -235,6 +241,8 @@ FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "SubCe FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "SubCellType", add_box = TRUE) FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "SubCellType", add_point = TRUE) FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "SubCellType", add_trend = TRUE) +FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "SubCellType", add_stat = "mean") +FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "SubCellType", add_line = 0.2, line_type = 2) FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "SubCellType", split.by = "Phase") FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "SubCellType", split.by = "Phase", add_box = TRUE, add_trend = TRUE) FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "SubCellType", split.by = "Phase", comparisons = TRUE) diff --git a/man/GroupHeatmap.Rd b/man/GroupHeatmap.Rd index 25e57f65..93586a4d 100644 --- a/man/GroupHeatmap.Rd +++ b/man/GroupHeatmap.Rd @@ -105,8 +105,8 @@ GroupHeatmap( cell_annotation = NULL, cell_annotation_palette = "Paired", cell_annotation_palcolor = NULL, - cell_annotation_params = if (flip) list(width = grid::unit(20, "mm")) else list(height - = grid::unit(20, "mm")), + cell_annotation_params = if (flip) list(width = grid::unit(10, "mm")) else list(height + = grid::unit(10, "mm")), feature_annotation = NULL, feature_annotation_palette = "Dark2", feature_annotation_palcolor = NULL, @@ -405,7 +405,7 @@ ht4 <- GroupHeatmap(pancreas_sub, features = de_top$gene, feature_split = de_top$group1, group.by = "CellType", heatmap_palette = "YlOrRd", cell_annotation = c("Phase", "G2M_score", "Neurod2"), cell_annotation_palette = c("Dark2", "Paired", "Paired"), - cell_annotation_params = list(height = grid::unit(20, "mm")), + cell_annotation_params = list(height = grid::unit(10, "mm")), feature_annotation = c("TF", "CSPA"), feature_annotation_palcolor = list(c("gold", "steelblue"), c("forestgreen")), add_dot = TRUE, add_bg = TRUE, nlabel = 0, show_row_names = TRUE @@ -416,7 +416,7 @@ ht5 <- GroupHeatmap(pancreas_sub, features = de_top$gene, feature_split = de_top$group1, group.by = "CellType", heatmap_palette = "YlOrRd", cell_annotation = c("Phase", "G2M_score", "Neurod2"), cell_annotation_palette = c("Dark2", "Paired", "Paired"), - cell_annotation_params = list(width = grid::unit(20, "mm")), + cell_annotation_params = list(width = grid::unit(10, "mm")), feature_annotation = c("TF", "CSPA"), feature_annotation_palcolor = list(c("gold", "steelblue"), c("forestgreen")), add_dot = TRUE, add_bg = TRUE, diff --git a/man/RunCellQC.Rd b/man/RunCellQC.Rd index e239d7ed..16b64a88 100644 --- a/man/RunCellQC.Rd +++ b/man/RunCellQC.Rd @@ -8,13 +8,13 @@ RunCellQC( srt, assay = "RNA", split.by = NULL, + return_filtered = FALSE, qc_metrics = c("doublets", "outlier", "umi", "gene", "mito", "ribo", "ribo_mito_ratio", "species"), - return_filtered = FALSE, db_method = "scDblFinder", db_rate = NULL, db_coefficient = 0.01, - outlier_cutoff = c("log10_nCount:lower:2.5", "log10_nCount:higher:5", + outlier_threshold = c("log10_nCount:lower:2.5", "log10_nCount:higher:5", "log10_nFeature:lower:2.5", "log10_nFeature:higher:5", "featurecount_dist:lower:2.5"), outlier_n = 1, UMI_threshold = 3000, @@ -40,18 +40,18 @@ RunCellQC( \item{split.by}{Name of the sample variable to split the Seurat object. Default is NULL.} +\item{return_filtered}{Logical indicating whether to return a cell-filtered Seurat object. Default is FALSE.} + \item{qc_metrics}{A character vector specifying the quality control metrics to be applied. Default is `c("doublets", "outlier", "umi", "gene", "mito", "ribo", "ribo_mito_ratio", "species")`.} -\item{return_filtered}{Logical indicating whether to return a cell-filtered Seurat object. Default is FALSE.} - \item{db_method}{Doublet-calling methods used. Can be one of \code{scDblFinder}, \code{Scrublet}, \code{DoubletDetection}, \code{scds_cxds}, \code{scds_bcds}, \code{scds_hybrid}} \item{db_rate}{The expected doublet rate. By default this is assumed to be 1\% per thousand cells captured (so 4\% among 4000 thousand cells), which is appropriate for 10x datasets.} \item{db_coefficient}{The coefficient used to calculate the doublet rate. Default is 0.01. Doublet rate is calculated as`ncol(srt) / 1000 * db_coefficient`} -\item{outlier_cutoff}{A character vector specifying the outlier cutoff values. Default is +\item{outlier_threshold}{A character vector specifying the outlier threshold. Default is `c("log10_nCount:lower:2.5", "log10_nCount:higher:5", "log10_nFeature:lower:2.5", "log10_nFeature:higher:5", "featurecount_dist:lower:2.5")`. See \link[scuttle]{isOutlier}.} \item{outlier_n}{Minimum number of outlier metrics that meet the conditions for determining outlier cells. Default is 1.} diff --git a/man/RunUMAP2.Rd b/man/RunUMAP2.Rd index 3d759b29..fb32bb47 100644 --- a/man/RunUMAP2.Rd +++ b/man/RunUMAP2.Rd @@ -19,6 +19,7 @@ RunUMAP2(object, ...) slot = "data", umap.method = "uwot", reduction.model = NULL, + n_threads = NULL, return.model = FALSE, n.neighbors = 30L, n.components = 2L, @@ -45,6 +46,7 @@ RunUMAP2(object, ...) assay = NULL, umap.method = "uwot", reduction.model = NULL, + n_threads = NULL, return.model = FALSE, n.neighbors = 30L, n.components = 2L,