Add code to run scran::scoreMarkers

code/method-pars.R

@@ -11,12 +11,37 @@ suppressMessages({
 # Read the overall config file.
 config <- yaml::read_yaml(here::here("config.yaml"))
 
-scran_pars <- expand_grid(
+findMarkers_pars <- expand_grid(
   method = "scran",
+  func = "findMarkers",
   test.type = c("t", "wilcox", "binom"),
   pval.type = c("any", "all", "some")
 ) %>%
-  mutate(file_name = paste0(method, "_", test.type, "_", pval.type))
+  mutate(file_name = paste0(method, "_", func, "_", test.type, "_", pval.type))
+
+scoreMarkers_pars <- expand_grid(
+  method = "scran",
+  func = "scoreMarkers",
+  metric = c(
+    "mean.logFC.cohen",
+    "min.logFC.cohen",
+    "median.logFC.cohen",
+    "max.logFC.cohen",
+    "rank.logFC.cohen",
+    "mean.AUC",
+    "min.AUC",
+    "max.AUC",
+    "rank.AUC",
+    "mean.logFC.detected",
+    "min.logFC.detected",
+    "median.logFC.detected",
+    "max.logFC.detected",
+    "rank.logFC.detected"
+  )
+) %>%
+  mutate(metric_readable = gsub("\\.", "_", metric)) %>%
+  mutate(file_name = paste0(method, "_", func, "_", metric_readable)) %>%
+  select(-metric_readable)
 
 seurat_pars <- expand_grid(
   method = "seurat",
@@ -104,7 +129,8 @@ cosg_pars <- expand_grid(
   mutate(file_name = paste0(method))
 
 method_pars_data <- list(
-  scran_pars,
+  findMarkers_pars,
+  scoreMarkers_pars,
   seurat_pars,
   lm_pars,
   edger_pars,

code/run-scran.R

@@ -10,6 +10,8 @@
 #'
 run_scran <- function(sce, pars) {
 
+  if (pars$func == "findMarkers") {
+
   times <- system.time({
     raw_result <- scran::findMarkers(
       sce,
@@ -48,26 +50,59 @@ run_scran <- function(sce, pars) {
         p_value_adj = FDR,
         gene = gene
       )
-
     }
 
-    # We do know the log fold change when not using the Wilcox test, but it
-    # doesn't add any extra information.
-    x <- dplyr::mutate(
-      x,
-      log_fc = 0,
-      scaled_statistic = 0
-    )
+    x <- dplyr::mutate(x, scaled_statistic = 0)
 
     x
   })
 
+  } else if (pars$func == "scoreMarkers") {
+
+    metric <- pars$metric
+    rank_by <- pars$rank_by
+
+    times <- system.time({
+      raw_result <- scran::scoreMarkers(sce)
+    })
+
+    result <- lapply(raw_result, function(x) {
+
+      x <- as.data.frame(x)
+      x <- x[, metric, drop = FALSE]
+      x <- tibble::rownames_to_column(x, var = "gene")
+      x <- tibble::as_tibble(x)
+
+      x <- dplyr::select(
+        x,
+        raw_statistic = !!metric,
+        gene = gene,
+      )
+
+      x <- dplyr::mutate(
+        x,
+        scaled_statistic = 0,
+        p_value = 0,
+        p_value_adj = 0
+      )
+
+      x <- dplyr::arrange(x, dplyr::desc(raw_statistic))
+    })
+  }
+
   cluster_names <- names(result)
   for (i in seq_along(result)) {
     result[[i]] <- dplyr::mutate(result[[i]], cluster = cluster_names[[i]])
   }
   result <- dplyr::bind_rows(!!!result)
 
+  # Add in one-vs-rest log fold changes.
+  result$log_fc <- calculate_log_fc(
+    sce,
+    genes = result$gene,
+    clusters = result$cluster
+  )
+
   validate_mgs_result(result)
 
   # Get the elapsed time.
@@ -75,5 +110,3 @@ run_scran <- function(sce, pars) {
 
   list(time = time, result = result, raw_result = raw_result, pars = pars)
 }
-
-