Skip to content

Select ATAC peaks by hvg_score and drop empty cells in process_dataset #78

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
rcannood merged 1 commit into from
Sep 30, 2026
Merged

Select ATAC peaks by hvg_score and drop empty cells in process_dataset #78

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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).
Expand Down
19 changes: 16 additions & 3 deletions src/data_processors/process_dataset/script.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]
}
Expand All @@ -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")
Expand Down
Loading