From 0b7a728a9e47175db10dec2d2786740ee64f60fa Mon Sep 17 00:00:00 2001
From: Robrecht Cannoodt <rcannood@gmail.com>
Date: Tue, 29 Sep 2026 14:15:52 +0200
Subject: [PATCH] Select ATAC peaks by hvg_score and drop empty cells in
 process_dataset

---
 CHANGELOG.md                                 |  6 ++++++
 src/data_processors/process_dataset/script.R | 19 ++++++++++++++++---
 2 files changed, 22 insertions(+), 3 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1cf2a082..9d121145 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,12 @@
 
 * Added `src/utils/exit_codes.py`, so components can mark themselves non-applicable for a dataset (PR #32).
 
+## MAJOR CHANGES
+
+* `process_dataset`: When the target modality is ATAC, keep the 10k peaks with the highest `hvg_score` instead of 10k random peaks, so the target is the part of the chromatin accessibility that actually differs between cells.
+
+* `process_dataset`: Remove cells without counts in either modality, after the peak selection. Empty cells have no profile to predict, and made `cellmapper_scvi` fail on `pbmc_cite` (43 cells without ADT counts).
+
 ## MINOR CHANGES
 
 * `file_train_mod1`, `file_train_mod2`, `file_test_mod1`, `file_test_mod2`: Declare `uns["modality"]`, which `process_dataset` already writes and which six methods and the `run_benchmark` workflow already read (PR #29).
diff --git a/src/data_processors/process_dataset/script.R b/src/data_processors/process_dataset/script.R
index 1ab7a2c1..336714d7 100644
--- a/src/data_processors/process_dataset/script.R
+++ b/src/data_processors/process_dataset/script.R
@@ -89,9 +89,8 @@ if (ad1_mod == "ATAC" && "gene_activity" %in% names(ad1$obsm)) {
 if (ad2_mod == "ATAC") {
   # subset to make the task computationally feasible
   if (ncol(ad2) > 10000) {
-    poss_ix <- which(Matrix::colSums(ad2$layers[["normalized"]]) > 0)
-    # sample by position -- sample(x, n) treats a length-1 x as seq_len(x)
-    sel_ix <- sort(poss_ix[sample.int(length(poss_ix), min(10000, length(poss_ix)))])
+    # keep the most variable peaks
+    sel_ix <- sort(order(ad2$var[["hvg_score"]], decreasing = TRUE)[seq_len(10000)])
     ad2 <- ad2[, sel_ix]$copy()
     ad2_var <- ad2_var[sel_ix, , drop = FALSE]
   }
@@ -103,6 +102,20 @@ if (ad2_mod == "ATAC") {
   }
 }
 
+# drop cells without counts in either modality -- after the peak selection above,
+# which can leave cells empty in the ATAC modality
+keep <- which(
+  Matrix::rowSums(ad1$layers[["counts"]]) > 0 &
+    Matrix::rowSums(ad2$layers[["counts"]]) > 0
+)
+if (length(keep) < nrow(ad1)) {
+  cat("Removing ", nrow(ad1) - length(keep), " cells without counts in either modality\n", sep = "")
+  ad1 <- ad1[keep, ]$copy()
+  ad2 <- ad2[keep, ]$copy()
+  ad1_obsm <- lapply(ad1_obsm, function(x) x[keep, , drop = FALSE])
+  ad2_obsm <- lapply(ad2_obsm, function(x) x[keep, , drop = FALSE])
+}
+
 cat("Creating train/test split\n")
 if ("is_train" %in% colnames(ad1$obs)) {
   is_train <- which(ad1$obs[["is_train"]] == "train")