diff --git a/.Rbuildignore b/.Rbuildignore
index f977d10..703fee7 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -3,3 +3,4 @@
 LICENSE
 inst/hg38*
 README.md
+tests/*
\ No newline at end of file
diff --git a/DESCRIPTION b/DESCRIPTION
index 00b3e24..95aa9de 100644
--- a/DESCRIPTION
+++ b/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
diff --git a/R/fastCorrelation.R b/R/fastCorrelation.R
index a73d633..75d0db0 100644
--- a/R/fastCorrelation.R
+++ b/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
 }
diff --git a/R/makeSuperCells.R b/R/makeSuperCells.R
index 42397d6..9ab9801 100644
--- a/R/makeSuperCells.R
+++ b/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      <-
diff --git a/R/normalizeNetwork.R b/R/normalizeNetwork.R
index 19374e7..4372d6f 100644
--- a/R/normalizeNetwork.R
+++ b/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
 }
diff --git a/R/removeBatch.R b/R/removeBatch.R
index 587a79a..a1f84fe 100644
--- a/R/removeBatch.R
+++ b/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)
 }
diff --git a/R/runPANDA.R b/R/runPANDA.R
index 6440193..5cc1cba 100644
--- a/R/runPANDA.R
+++ b/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))
diff --git a/R/runSCORPION.R b/R/runSCORPION.R
index 9229585..a495769 100644
--- a/R/runSCORPION.R
+++ b/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,13 +330,33 @@ 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 {
@@ -322,30 +364,44 @@ runSCORPION <- function(gexMatrix,
     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),
diff --git a/R/scorpion.R b/R/scorpion.R
index c98b9f6..2929a46 100644
--- a/R/scorpion.R
+++ b/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, ]
   }
diff --git a/R/tanimotoSimilarity.R b/R/tanimotoSimilarity.R
index 67a61c5..1dfab6b 100644
--- a/R/tanimotoSimilarity.R
+++ b/R/tanimotoSimilarity.R
@@ -19,8 +19,13 @@ tanimoto <- function(X, Y = NULL,
     if (is.null(x_norm_sq)) x_norm_sq <- rowSums(X * X)
   }
 
-  den <- outer(x_norm_sq, y_norm_sq, "+") - abs(Amat)
-  Amat / sqrt(den)
+  # Column-by-column normalization to avoid full-size outer() temporary
+  nc <- ncol(Amat)
+  for (j in seq_len(nc)) {
+    col_j <- Amat[, j]
+    Amat[, j] <- col_j / sqrt(x_norm_sq + y_norm_sq[j] - abs(col_j))
+  }
+  Amat
 }
 
 #' @importFrom utils getFromNamespace
@@ -49,6 +54,10 @@ tanimoto_gpu <- function(X, Y = NULL,
     if (is.null(x_norm_sq)) x_norm_sq <- rowSums(X * X)
   }
 
-  den <- outer(x_norm_sq, y_norm_sq, "+") - abs(Amat)
-  Amat / sqrt(den)
+  nc <- ncol(Amat)
+  for (j in seq_len(nc)) {
+    col_j <- Amat[, j]
+    Amat[, j] <- col_j / sqrt(x_norm_sq + y_norm_sq[j] - abs(col_j))
+  }
+  Amat
 }
diff --git a/R/testEdges.R b/R/testEdges.R
index 622903a..22f9338 100644
--- a/R/testEdges.R
+++ b/R/testEdges.R
@@ -55,7 +55,7 @@
 #'     \item{tStatistic: Test statistic}
 #'     \item{pValue: Raw p-value}
 #'     \item{pAdj: Adjusted p-value}
-#'     \item{For two-sample tests: meanGroup1, meanGroup2, diffMean (Group1 - Group2), cohensD, log2FoldChange}
+#'     \item{For two-sample tests: meanGroup1, meanGroup2, cohensD, log2FoldChange (Group1 - Group2)}
 #'   }
 #' @details
 #' For single-sample tests, the function tests whether the mean edge weight across
@@ -278,6 +278,9 @@ testEdges <- function(networksDF,
     
     # Set up parallel plan; use sequential reset as safety net on exit
     # to guarantee worker processes are killed even if an error occurs
+    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(future::sequential)
@@ -295,6 +298,7 @@ testEdges <- function(networksDF,
     # then restore the caller's original plan
     future::plan(future::sequential)
     future::plan(old_plan)
+    options(future.globals.maxSize = old_maxsize)
     on.exit()  # cancel the on.exit guard since cleanup is done
   }
   
@@ -452,10 +456,10 @@ testEdgesTwoSample <- function(networksDF, group1, group2, alternative, minLog2F
   meanEdge <- (meanEdge1 + meanEdge2) / 2
   diffMean <- meanEdge1 - meanEdge2
   
-  # Calculate log2 fold change from quantiles (needed for filtering)
-  # Convert z-scores to quantiles using pnorm with log.p=TRUE for precision
-  # log2(q1/q2) = (log(q1) - log(q2)) / log(2)
-  log2FC <- (pnorm(meanEdge1, log.p = TRUE) - pnorm(meanEdge2, log.p = TRUE)) / log(2)
+  # limma-style log2 fold change: limma's logFC is the model coefficient, i.e. the
+  # difference of group means on log2-scale input (defined for all reals, incl.
+  # negative means). Treating edge weights as log2-scale, this equals diffMean.
+  log2FC <- meanEdge1 - meanEdge2
   
   # Filter by minimum log2 fold change
   keep_idx <- abs(log2FC) >= minLog2FC
@@ -525,7 +529,6 @@ testEdgesTwoSample <- function(networksDF, group1, group2, alternative, minLog2F
     target = tf_target$target,
     meanGroup1 = meanEdge1,
     meanGroup2 = meanEdge2,
-    diffMean = diffMean,
     cohensD = cohensD,
     log2FoldChange = log2FC,
     meanEdge = meanEdge,
@@ -556,10 +559,10 @@ testEdgesPaired <- function(networksDF, group1, group2, alternative, minLog2FC,
   diff_matrix <- edge_matrix1 - edge_matrix2
   diffMean <- rowMeans(diff_matrix, na.rm = TRUE)
   
-  # Calculate log2 fold change from quantiles (needed for filtering)
-  # Convert z-scores to quantiles using pnorm with log.p=TRUE for precision
-  # log2(q1/q2) = (log(q1) - log(q2)) / log(2)
-  log2FC <- (pnorm(meanEdge1, log.p = TRUE) - pnorm(meanEdge2, log.p = TRUE)) / log(2)
+  # limma-style log2 fold change: limma's logFC is the model coefficient, i.e. the
+  # difference of group means on log2-scale input (defined for all reals, incl.
+  # negative means). Treating edge weights as log2-scale, this equals diffMean.
+  log2FC <- meanEdge1 - meanEdge2
   
   # Filter by minimum log2 fold change
   keep_idx <- abs(log2FC) >= minLog2FC
@@ -626,7 +629,6 @@ testEdgesPaired <- function(networksDF, group1, group2, alternative, minLog2FC,
     target = tf_target$target,
     meanGroup1 = meanEdge1,
     meanGroup2 = meanEdge2,
-    diffMean = diffMean,
     cohensD = cohensD,
     log2FoldChange = log2FC,
     meanEdge = meanEdge,
diff --git a/README.md b/README.md
index d4b2c27..bfcb7d0 100644
--- a/README.md
+++ b/README.md
@@ -298,9 +298,8 @@ A data frame containing:
 | `tf`, `target` | TF-target pair identifiers |
 | `meanEdge` | Mean edge weight (single-sample) |
 | `meanGroup1`, `meanGroup2` | Group means (two-sample) |
-| `diffMean` | Difference in means, Group1 − Group2 (two-sample) |
 | `cohensD` | Cohen's d effect size (two-sample and paired tests) |
-| `log2FoldChange` | Log2 fold change (two-sample) |
+| `log2FoldChange` | Log2 fold change, Group1 − Group2 (two-sample) |
 | `tStatistic` | t-statistic |
 | `pValue` | Raw p-value |
 | `pAdj` | Adjusted p-value |
diff --git a/man/testEdges.Rd b/man/testEdges.Rd
index 50e0719..b9ee5e3 100644
--- a/man/testEdges.Rd
+++ b/man/testEdges.Rd
@@ -84,7 +84,7 @@ A data.frame containing:
     \item{tStatistic: Test statistic}
     \item{pValue: Raw p-value}
     \item{pAdj: Adjusted p-value}
-    \item{For two-sample tests: meanGroup1, meanGroup2, diffMean (Group1 - Group2), cohensD, log2FoldChange}
+    \item{For two-sample tests: meanGroup1, meanGroup2, cohensD, log2FoldChange (Group1 - Group2)}
   }
 }
 \description{
diff --git a/tests/testthat/test_edge_statistics.R b/tests/testthat/test_edge_statistics.R
index 881eaba..bcac885 100644
--- a/tests/testthat/test_edge_statistics.R
+++ b/tests/testthat/test_edge_statistics.R
@@ -282,7 +282,7 @@ test_that("testEdges two-sample mean difference is correct", {
     result_row <- results[results$tf == mock$df$tf[i] & 
                            results$target == mock$df$target[i], ]
     
-    expect_equal(result_row$diffMean, expected_diff, 
+    expect_equal(result_row$log2FoldChange, expected_diff, 
                  tolerance = 1e-10,
                  label = paste("mean difference for edge", i))
     expect_equal(result_row$meanGroup1, mean(edge_vals_g1), 
@@ -449,18 +449,18 @@ test_that("testEdges paired mean difference is correct", {
     empiricalNull = FALSE
   )
   
-  # Check that diffMean equals mean of differences (not difference of means)
+  # Check that log2FoldChange equals mean of differences (not difference of means)
   for (i in 1:10) {
     edge_vals_g1 <- as.numeric(mock$df[i, mock$group1])
     edge_vals_g2 <- as.numeric(mock$df[i, mock$group2])
     
-    # For paired test, diffMean should be mean(g1 - g2)
+    # For paired test, log2FoldChange should be mean(g1 - g2)
     expected_diff <- mean(edge_vals_g1 - edge_vals_g2)
     
     result_row <- results[results$tf == mock$df$tf[i] & 
                            results$target == mock$df$target[i], ]
     
-    expect_equal(result_row$diffMean, expected_diff, 
+    expect_equal(result_row$log2FoldChange, expected_diff, 
                  tolerance = 1e-10,
                  label = paste("paired mean difference for edge", i))
   }
@@ -867,7 +867,7 @@ test_that("testEdges parallel two-sample matches serial and t.test()", {
   
   expect_equal(results_serial$tStatistic, results_parallel$tStatistic, tolerance = 1e-10)
   expect_equal(results_serial$pValue, results_parallel$pValue, tolerance = 1e-10)
-  expect_equal(results_serial$diffMean, results_parallel$diffMean, tolerance = 1e-10)
+  expect_equal(results_serial$log2FoldChange, results_parallel$log2FoldChange, tolerance = 1e-10)
   expect_equal(results_serial$pAdj, results_parallel$pAdj, tolerance = 1e-10)
   
   # Parallel results should still match t.test()
@@ -915,7 +915,7 @@ test_that("testEdges parallel paired matches serial and t.test()", {
   
   expect_equal(results_serial$tStatistic, results_parallel$tStatistic, tolerance = 1e-10)
   expect_equal(results_serial$pValue, results_parallel$pValue, tolerance = 1e-10)
-  expect_equal(results_serial$diffMean, results_parallel$diffMean, tolerance = 1e-10)
+  expect_equal(results_serial$log2FoldChange, results_parallel$log2FoldChange, tolerance = 1e-10)
   expect_equal(results_serial$pAdj, results_parallel$pAdj, tolerance = 1e-10)
   
   # Parallel results should still match t.test(paired = TRUE)
diff --git a/tests/testthat/test_mori_sharing.R b/tests/testthat/test_mori_sharing.R
new file mode 100644
index 0000000..def8eb3
--- /dev/null
+++ b/tests/testthat/test_mori_sharing.R
@@ -0,0 +1,113 @@
+# Tests for optional 'mori' shared-memory acceleration of runSCORPION().
+# All tests skip unless mori (and the parallel stack) is installed.
+
+test_that("runSCORPION() with mori matches results without mori", {
+  skip_on_cran()
+  skip_if_not_installed("mori")
+  skip_if_not_installed("furrr")
+  skip_if_not_installed("future")
+
+  data(scorpionTest)
+
+  # Groups must each have enough cells to build a network (>= 30 cells).
+  groups <- table(scorpionTest$metadata$region)
+  skip_if(sum(groups >= 30) < 2, "Need >= 2 groups with >= 30 cells for parallel test")
+
+  old_opt <- getOption("scorpion.use_mori")
+  on.exit(options(scorpion.use_mori = old_opt), add = TRUE)
+
+  run <- function() {
+    set.seed(1)
+    runSCORPION(
+      gexMatrix = scorpionTest$gex,
+      tfMotifs = scorpionTest$tf,
+      ppiNet = scorpionTest$ppi,
+      cellsMetadata = scorpionTest$metadata,
+      groupBy = "region",
+      alphaValue = 0.8,
+      nCores = 2L,
+      showProgress = FALSE
+    )
+  }
+
+  options(scorpion.use_mori = FALSE)
+  res_plain <- run()
+
+  options(scorpion.use_mori = TRUE)
+  res_mori <- run()
+
+  expect_equal(dim(res_plain), dim(res_mori))
+  expect_equal(colnames(res_plain), colnames(res_mori))
+  expect_equal(res_plain$tf, res_mori$tf)
+  expect_equal(res_plain$target, res_mori$target)
+
+  weight_cols <- setdiff(colnames(res_plain), c("tf", "target"))
+  for (col in weight_cols) {
+    expect_equal(res_plain[[col]], res_mori[[col]],
+                 tolerance = 1e-10,
+                 label = paste("weights for network", col))
+  }
+})
+
+test_that("mori::share() emits a compact serialized reference", {
+  skip_on_cran()
+  skip_if_not_installed("mori")
+
+  # The whole point of mori: a shared object serializes to a tiny reference
+  # instead of a full copy, so it is cheap to broadcast to every worker.
+  x <- as.data.frame(matrix(rnorm(2e5), ncol = 4))
+  shared <- mori::share(x)
+
+  size_plain <- length(serialize(x, NULL))
+  size_shared <- length(serialize(shared, NULL))
+
+  expect_lt(size_shared, size_plain)
+  expect_lt(size_shared, 10000)  # a few bytes / KB, not the full payload
+})
+
+test_that("runSCORPION() parallel performance with vs without mori", {
+  skip_on_cran()
+  skip_if_not_installed("mori")
+  skip_if_not_installed("furrr")
+  skip_if_not_installed("future")
+
+  data(scorpionTest)
+
+  groups <- table(scorpionTest$metadata$region)
+  skip_if(sum(groups >= 30) < 2, "Need >= 2 groups with >= 30 cells for parallel test")
+
+  old_opt <- getOption("scorpion.use_mori")
+  on.exit(options(scorpion.use_mori = old_opt), add = TRUE)
+
+  run <- function() {
+    runSCORPION(
+      gexMatrix = scorpionTest$gex,
+      tfMotifs = scorpionTest$tf,
+      ppiNet = scorpionTest$ppi,
+      cellsMetadata = scorpionTest$metadata,
+      groupBy = "region",
+      alphaValue = 0.8,
+      nCores = 2L,
+      showProgress = FALSE
+    )
+  }
+
+  # Warm up worker processes so plan startup cost is not attributed to a run.
+  invisible(run())
+
+  options(scorpion.use_mori = FALSE)
+  t_plain <- system.time(run())[["elapsed"]]
+
+  options(scorpion.use_mori = TRUE)
+  t_mori <- system.time(run())[["elapsed"]]
+
+  message(sprintf(
+    "runSCORPION parallel timing: without mori = %.2fs, with mori = %.2fs (%.2fx)",
+    t_plain, t_mori, t_plain / t_mori
+  ))
+
+  # Both paths must complete successfully; timing is reported, not asserted,
+  # since wall-clock speedup depends on prior size, worker count and hardware.
+  expect_true(is.finite(t_plain) && t_plain >= 0)
+  expect_true(is.finite(t_mori) && t_mori >= 0)
+})