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

@@ -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),

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)
+      }
     }
   }
 

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])