1.3.3

.Rbuildignore

@@ -3,3 +3,4 @@
 LICENSE
 inst/hg38*
 README.md
+tests/*

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: SCORPION
 Type: Package
 Title: Single Cell Oriented Reconstruction of PANDA Individually Optimized Networks
-Version: 1.3.2
+Version: 1.3.3
 Authors@R: c(person("Daniel","Osorio", email = "daniecos@uio.no", role = c("aut","cre"), comment = c(ORCID = "0000-0003-4424-8422")),
   person("Marieke L.","Kuijjer", email = "marieke.kuijjer@ncmbm.uio.no", role = c("aut"), comment = c(ORCID = "0000-0001-6280-3130")))
 Description: Constructs cell-type–specific gene regulatory networks from single-cell RNA-sequencing data. The method implements the SCORPION algorithm, which first aggregates individual cells into super-cells and then applies PANDA (Passing Attributes between Networks for Data Assimilation) to infer transcription factor–target regulatory relationships. It also provides statistical methods for differential edge analysis.
@@ -10,7 +10,7 @@ Encoding: UTF-8
 LazyData: true
 Depends: R (>= 3.5.0)
 Imports: cli, methods, irlba, igraph, RANN, Matrix, pbapply, dplyr, furrr, future
-Suggests: RhpcBLASctl, testthat
+Suggests: RhpcBLASctl, testthat, mori
 URL: https://github.com/kuijjerlab/SCORPION
 BugReports: https://github.com/kuijjerlab/SCORPION/issues
 RoxygenNote: 7.3.3

R/fastCorrelation.R

@@ -21,6 +21,10 @@ fastCorrelation <- function(X, Y, method = 'pearson'){
   nx <- sqrt(colSums(RX * RX))
   ny <- sqrt(colSums(RY * RY))
 
-  # Correlation
-  crossprod(RX, RY) / (nx %o% ny)
+  # Correlation — column-by-column division avoids full-size outer product
+  result <- crossprod(RX, RY)
+  for (j in seq_len(ncol(result))) {
+    result[, j] <- result[, j] / (nx * ny[j])
+  }
+  result
 }

R/makeSuperCells.R

@@ -137,6 +137,7 @@ makeSuperCells <- function(X,
       PCA.presampled$u %*% diag(PCA.presampled$d)
     PCA.presampled$rotation <- PCA.presampled$v
   }
+  rm(X.for.pca)
 
   sc.nw <-
     buildKNN(
@@ -147,6 +148,7 @@ makeSuperCells <- function(X,
       dist_method = "euclidean",
       directed = directed
     )
+  rm(PCA.presampled)
 
   # simplify
 
@@ -179,6 +181,7 @@ makeSuperCells <- function(X,
     igraph::contract(sc.nw$graph.knn, membership.presampled)
   SC.NW <-
     igraph::simplify(SC.NW, remove.loops = T, edge.attr.comb = "sum")
+  rm(g.s, sc.nw, SC.NW)
 
   if (do.approx) {
     # PCA.averaged.SC      <-

R/normalizeNetwork.R

@@ -6,30 +6,42 @@ normalizeNetwork <- function(X) {
   mu0 <- mean(X)
   std0 <- sd(X)
 
-  # Row normalization (R's column-major recycling handles row-wise broadcast)
+  # Row stats
   mu1 <- rowMeans(X)
   std1 <- rowSds(X) * sqrt((nc - 1) / nc)
-  Z1 <- (X - mu1) / std1
 
-  # Column normalization (needs explicit broadcast)
+  # Column stats (vectors only — no full-matrix broadcast)
   mu2 <- colMeans(X)
   std2 <- colSds(X) * sqrt((nr - 1) / nr)
-  mu2 <- rep(mu2, each = nr)
-  dim(mu2) <- c(nr, nc)
-  std2 <- rep(std2, each = nr)
-  dim(std2) <- c(nr, nc)
-  Z2 <- (X - mu2) / std2
-
-  # combine and return
-  normMat <- (Z1 + Z2) / sqrt(2)
-
-  Z0 <- (X - mu0) / std0
-  f1 <- is.na(Z1)
-  f2 <- is.na(Z2)
-
-  normMat[f1] <- (Z2[f1] + Z0[f1]) / sqrt(2)
-  normMat[f2] <- (Z1[f2] + Z0[f2]) / sqrt(2)
-  normMat[f1 & f2] <- 2 * Z0[f1 & f2] / sqrt(2)
+
+  has_zero_var <- any(std1 == 0) || any(std2 == 0)
+
+  if (!has_zero_var) {
+    # Fast path: accumulate Z1 + Z2 in a single matrix
+    normMat <- (X - mu1) / std1
+    for (j in seq_len(nc)) {
+      normMat[, j] <- normMat[, j] + (X[, j] - mu2[j]) / std2[j]
+    }
+    normMat <- normMat / sqrt(2)
+  } else {
+    # Slow path: handle NaN from zero-variance rows/columns
+    Z1 <- (X - mu1) / std1
+
+    Z2 <- X
+    for (j in seq_len(nc)) {
+      Z2[, j] <- (X[, j] - mu2[j]) / std2[j]
+    }
+
+    normMat <- (Z1 + Z2) / sqrt(2)
+
+    Z0 <- (X - mu0) / std0
+    f1 <- is.na(Z1)
+    f2 <- is.na(Z2)
+
+    normMat[f1] <- (Z2[f1] + Z0[f1]) / sqrt(2)
+    normMat[f2] <- (Z1[f2] + Z0[f2]) / sqrt(2)
+    normMat[f1 & f2] <- 2 * Z0[f1 & f2] / sqrt(2)
+  }
 
   normMat
 }

R/removeBatch.R

@@ -1,18 +1,34 @@
 remove_batch <- function(X, batch) {
-    batch <- as.factor(batch)
+    batch <- droplevels(as.factor(batch))
+
+    # With fewer than two batch levels there is nothing to correct
+    if (nlevels(batch) < 2) {
+        return(X)
+    }
 
     # Design matrix
     H <- model.matrix(~batch)
 
     # Solve coefficients: beta = (H'H)^(-1) H' X'
-    HtH_inv <- solve(crossprod(H))
-    HtX <- crossprod(H, t(X))
-    beta <- HtH_inv %*% HtX
+    # Fall back to a Moore-Penrose pseudo-inverse when H'H is rank-deficient
+    # (e.g. collinear or singleton batches) so the run is not aborted with a
+    # Lapack 'system is exactly singular' error.
+    HtH <- crossprod(H)
+    HtH_inv <- tryCatch(
+        solve(HtH),
+        error = function(e) {
+            s <- svd(HtH)
+            tol <- max(dim(HtH)) * .Machine$double.eps * max(s$d)
+            keep <- s$d > tol
+            s$v[, keep, drop = FALSE] %*%
+                ((1 / s$d[keep]) * t(s$u[, keep, drop = FALSE]))
+        }
+    )
+    beta <- HtH_inv %*% t(X %*% H)
 
-    # Correction: H * beta, then subtract
-    X <- t(t(X) - H %*% beta)
+    # Correction: subtract H %*% beta from X
+    X <- X - tcrossprod(t(beta), H)
 
-    # Clear memory
     gc()
     return(X)
 }

R/runPANDA.R

@@ -151,6 +151,7 @@ runPANDA <- function(motif = NULL, expr = NULL, ppi = NULL, alpha = 0.1, hamming
       cli::cli_alert_success("Verified sufficient samples")
     }
   }
+  rm(expr)
   if (any(is.na(geneCoreg))) {
     # check for NA and replace them by zero
     diag(geneCoreg) <- 1
@@ -197,7 +198,6 @@ runPANDA <- function(motif = NULL, expr = NULL, ppi = NULL, alpha = 0.1, hamming
     tanimoto_fn <- tanimoto
   }
 
-  minusAlpha <- 1 - alpha
   step <- 0
   hamming_cur <- 1
   if (progress) {
@@ -209,27 +209,33 @@ runPANDA <- function(motif = NULL, expr = NULL, ppi = NULL, alpha = 0.1, hamming
       cli::cli_alert_warning(paste0("Reached maximum iterations, iter =", iter))
       break
     }
-    # Precompute squared norms for regulatoryNetwork (reused across tanimoto calls)
-    reg_row_sq <- rowSums(regulatoryNetwork * regulatoryNetwork)
-    reg_col_sq <- colSums(regulatoryNetwork * regulatoryNetwork)
+    # Precompute squared norms once (shared across tanimoto calls)
+    reg_sq <- regulatoryNetwork * regulatoryNetwork
+    reg_row_sq <- rowSums(reg_sq)
+    reg_col_sq <- colSums(reg_sq)
+    rm(reg_sq)
 
-    Responsibility <- tanimoto_fn(tfCoopNetwork, regulatoryNetwork,
+    # Combine Responsibility + Availability directly into RA
+    RA <- tanimoto_fn(tfCoopNetwork, regulatoryNetwork,
                                y_norm_sq = reg_col_sq)
-    Availability <- tanimoto_fn(regulatoryNetwork, geneCoreg,
+    RA <- RA + tanimoto_fn(regulatoryNetwork, geneCoreg,
                              x_norm_sq = reg_row_sq)
-    RA <- 0.5 * (Responsibility + Availability)
+    RA <- 0.5 * RA
 
     hamming_cur <- sum(abs(regulatoryNetwork - RA)) / (num.TFs * num.genes)
-    regulatoryNetwork <- minusAlpha * regulatoryNetwork + alpha * RA
+    regulatoryNetwork <- regulatoryNetwork + alpha * (RA - regulatoryNetwork)
+    rm(RA)
 
     # tcrossprod/crossprod use optimized BLAS and avoid explicit t()
     ppi <- tanimoto_fn(regulatoryNetwork, type = "tcrossprod")
     ppi <- update.diagonal(ppi, num.TFs, alpha, step)
-    tfCoopNetwork <- minusAlpha * tfCoopNetwork + alpha * ppi
+    tfCoopNetwork <- tfCoopNetwork + alpha * (ppi - tfCoopNetwork)
+    rm(ppi)
 
     CoReg2 <- tanimoto_fn(regulatoryNetwork, type = "crossprod")
     CoReg2 <- update.diagonal(CoReg2, num.genes, alpha, step)
-    geneCoreg <- minusAlpha * geneCoreg + alpha * CoReg2
+    geneCoreg <- geneCoreg + alpha * (CoReg2 - geneCoreg)
+    rm(CoReg2)
 
     if (progress) {
       # message("Iteration ", step,": hamming distance = ", round(hamming_cur,5))

R/runSCORPION.R

@@ -239,6 +239,13 @@ runSCORPION <- function(gexMatrix,
     cli::cli_h1("SCORPION")
   }
 
+  # Control BLAS threading to respect nCores
+  if (requireNamespace("RhpcBLASctl", quietly = TRUE)) {
+    old_blas <- RhpcBLASctl::blas_get_num_procs()
+    RhpcBLASctl::blas_set_num_threads(nCores)
+    on.exit(RhpcBLASctl::blas_set_num_threads(old_blas), add = TRUE)
+  }
+
   # Normalizing data
   if (normalizeData) {
     if (showProgress) {
@@ -259,6 +266,12 @@ runSCORPION <- function(gexMatrix,
     gexMatrix <- remove_batch(X = gexMatrix, batch = batch)
     gexMatrix <- gexMatrix + mean_expr
   }
+  rm(batch)
+
+  # Pre-convert to data.frame once so workers receive the converted objects
+  # instead of converting on every call
+  tfMotifs <- as.data.frame(tfMotifs)
+  ppiNet <- as.data.frame(ppiNet)
 
   # Setting min number of cells to construct network
   min_cells <- max(minCells, 30)
@@ -272,6 +285,7 @@ runSCORPION <- function(gexMatrix,
   } else {
     cli::cli_abort('groupBy must match cellsMetadata column name')
   }
+  rm(cellsMetadata)
 
   total_net <- length(unique(metadata$network_id))
   metadata <- metadata %>% filter(.data$n_cells >= min_cells)
@@ -282,17 +296,25 @@ runSCORPION <- function(gexMatrix,
     cli::cli_alert_success(paste0(filtered_net, " networks meet the minimum cell requirement (", min_cells, ")"))
   }
 
-  compute_network <- function(idx) {
-    selected_network <- network_ids[idx]
+  if (filtered_net == 0) {
+    cli::cli_abort("No groups have enough cells (>= {min_cells}) to build a network")
+  }
 
-    selected_cells <- metadata %>%
-      filter(.data$network_id %in% selected_network)
-    selected_cells <- gexMatrix[, selected_cells$cell_id]
+  network_ids <- unique(metadata$network_id)
 
+  # Pre-split gexMatrix into per-group chunks so each worker only receives
+  # the subset it needs, instead of the full matrix.
+  gex_chunks <- lapply(network_ids, function(nid) {
+    cells <- metadata$cell_id[metadata$network_id == nid]
+    gexMatrix[, cells, drop = FALSE]
+  })
+  rm(gexMatrix, metadata)
+
+  compute_network <- function(gex_chunk) {
     network <- scorpion(
-      gexMatrix = selected_cells,
-      tfMotifs = as.data.frame(tfMotifs),
-      ppiNet = as.data.frame(ppiNet),
+      gexMatrix = gex_chunk,
+      tfMotifs = tfMotifs,
+      ppiNet = ppiNet,
       computingEngine = computingEngine,
       nCores = 1,
       gammaValue = gammaValue,
@@ -308,44 +330,78 @@ runSCORPION <- function(gexMatrix,
       scaleByPresent = scaleByPresent,
       filterExpr = filterExpr
     )[[outNet]]
-    
+
     return(network)
   }
 
-  network_ids <- unique(metadata$network_id)
+  # The TF-motif prior and PPI network are broadcast to every parallel worker.
+  # If the optional 'mori' package is available, place them in OS-backed shared
+  # memory so workers map the same physical pages instead of each receiving a
+  # full serialized copy. Falls back to standard serialization when absent, or
+  # when disabled via options(scorpion.use_mori = FALSE).
+  use_mori <- nCores > 1 &&
+    isTRUE(getOption("scorpion.use_mori", TRUE)) &&
+    requireNamespace("mori", quietly = TRUE)
+  if (use_mori) {
+    if (!is.null(tfMotifs)) tfMotifs <- mori::share(tfMotifs)
+    if (!is.null(ppiNet)) ppiNet <- mori::share(ppiNet)
+  }
+
+  furrr_opts <- furrr::furrr_options(
+    seed = TRUE,
+    packages = if (use_mori) "mori" else NULL
+  )
 
   if (nCores > 1) {
+    # Allow arbitrarily large globals to be exported to workers.
+    old_maxsize <- getOption("future.globals.maxSize")
+    options(future.globals.maxSize = Inf)
+    on.exit(options(future.globals.maxSize = old_maxsize), add = TRUE)
     old_plan <- future::plan(future::multisession, workers = nCores)
     on.exit(future::plan(old_plan), add = TRUE)
   } else {
     old_plan <- future::plan(future::sequential)
     on.exit(future::plan(old_plan), add = TRUE)
   }
 
+  n_total <- length(gex_chunks)
+
   if (showProgress) {
-    cli::cli_alert_info("Computing networks")
+    cli::cli_alert_info(paste0("Computing ", n_total, " networks"))
     if (nCores > 1) {
       cli::cli_alert_info(paste0("Using ", nCores, " cores for parallel processing"))
+      network_matrices <- furrr::future_map(gex_chunks, compute_network, .options = furrr_opts, .progress = FALSE)
+    } else {
+      network_matrices <- vector("list", n_total)
+      for (i in seq_len(n_total)) {
+        cli::cli_alert_info(paste0("Network ", i, "/", n_total, ": ", network_ids[i]))
+        network_matrices[[i]] <- compute_network(gex_chunks[[i]])
+      }
     }
-    network_matrices <- furrr::future_map(seq_along(network_ids), compute_network, .options = furrr::furrr_options(seed = TRUE), .progress = TRUE)
     cli::cli_alert_success("Networks successfully constructed")
   } else {
-    network_matrices <- furrr::future_map(seq_along(network_ids), compute_network, .options = furrr::furrr_options(seed = TRUE), .progress = TRUE)
+    network_matrices <- furrr::future_map(gex_chunks, compute_network, .options = furrr_opts, .progress = FALSE)
   }
+  rm(gex_chunks)
 
-  # Build TF-target pairs from first network using same method as before
+  # Build TF-target pairs from first network
   first_net <- network_matrices[[1]]
   tf_target_df <- as.data.frame(as.table(first_net))[, 1:2]
   colnames(tf_target_df) <- c("tf", "target")
-
-  # Extract weights as vectors (column-major order matches as.table order)
-  weight_matrix <- vapply(
-    network_matrices,
-    as.vector,
-    numeric(length(first_net))
-  )
+  n_edges <- length(first_net)
+  rm(first_net)
+
+  # Stream extraction: pull each network's weights into pre-allocated matrix,
+  # then NULL out the list element immediately to free memory.
+  n_nets <- length(network_matrices)
+  weight_matrix <- matrix(NA_real_, nrow = n_edges, ncol = n_nets)
   colnames(weight_matrix) <- network_ids
-
+  for (k in seq_len(n_nets)) {
+    weight_matrix[, k] <- as.vector(network_matrices[[k]])
+    network_matrices[k] <- list(NULL)
+  }
+  rm(network_matrices)
+
   # Combine into final data frame
   networks <- data.frame(
     tf = as.character(tf_target_df$tf),

R/scorpion.R

@@ -148,6 +148,13 @@ scorpion <- function(tfMotifs = NULL,
     cli::cli_h1("SCORPION")
   }
 
+  # Control BLAS threading to respect nCores
+  if (requireNamespace("RhpcBLASctl", quietly = TRUE)) {
+    old_blas <- RhpcBLASctl::blas_get_num_procs()
+    RhpcBLASctl::blas_set_num_threads(nCores)
+    on.exit(RhpcBLASctl::blas_set_num_threads(old_blas), add = TRUE)
+  }
+
   if (isTRUE(filterExpr)) {
     gexMatrix <- gexMatrix[rowSums(gexMatrix) > 0, ]
   }