From cc6e7ba127baac548d81830843ecfbc13bab75da Mon Sep 17 00:00:00 2001 From: aviezerl Date: Tue, 14 Jul 2026 00:20:27 +0300 Subject: [PATCH] fix: score each sequence independently in compute_pwm compute_pwm() capped the motif-scan window to nchar(sequences[1]) for the whole batch, so a sequence longer than the first was only scanned over its first few positions and could miss its real motif hit - making a sequence's score depend on its batch companions. Use max(nchar(sequences)) for both the scan range and the flat spatial bin; the C++ scan already clamps per-sequence to each sequence's own end, so every sequence is now scored on its full length. Equal-length batches (regression / screen_kmers) are unchanged since max == first. Claude-Session: https://claude.ai/code/session_01PK3qefBGDoBwd9w5226FEq --- NEWS.md | 6 ++++++ R/pssm-utils.R | 8 ++++++-- tests/testthat/test-compute_pwm.R | 15 +++++++++++++++ 3 files changed, 27 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index c940bd9..824dcab 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,11 @@ # prego 0.0.10 +* Fix: `compute_pwm()` gave a sequence a different score depending on the other + sequences in the batch. The motif-scan window was capped to the length of the + *first* sequence for the whole batch, so a longer sequence sitting behind a + shorter one was only scanned over its first `nchar(sequences[1])` positions and + missed motif hits further along. Each sequence is now scored on its own full + length, independently of its batch companions. * Added `return_all` parameter to `regress_pwm` (multi-kmer path). When TRUE, returns every candidate-kmer regression (sorted by validation score) instead of just the best one - useful for getting N independent motifs without the residual-rounds approach used by `motif_num > 1`. When `sample_for_kmers = TRUE`, each candidate is refit on the full data. * Improved docs for `regress_pwm` (clarified the three operating modes, fixed `n_motifs`/`comb_modle` typos in the return-value section). diff --git a/R/pssm-utils.R b/R/pssm-utils.R index 654fd56..1374956 100644 --- a/R/pssm-utils.R +++ b/R/pssm-utils.R @@ -26,7 +26,9 @@ compute_pwm <- function(sequences, pssm, spat = NULL, spat_min = 1, spat_max = NULL, bidirect = TRUE, prior = 0.01, func = "logSumExp") { if (is.null(spat)) { spat <- data.frame(bin = 0, spat_factor = 1) - binsize <- nchar(sequences[[1]]) + # One flat bin spanning the longest sequence, so every sequence is + # scored on its own full length (bin index stays 0 regardless of length). + binsize <- max(nchar(sequences)) } else { validate_spat(spat) binsize <- unique(diff(spat$bin)) @@ -73,7 +75,9 @@ compute_pwm <- function(sequences, pssm, spat = NULL, spat_min = 1, spat_max = N pssm_mat = pssm_mat, is_bidirect = bidirect, spat_min = 0, - spat_max = nchar(sequences[1]), + # Longest sequence: the C++ scan clamps per sequence to its own end, + # so each is scored independently of others in the batch. + spat_max = max(nchar(sequences)), spat_factor = spat$spat_factor, bin_size = binsize, use_max = use_max diff --git a/tests/testthat/test-compute_pwm.R b/tests/testthat/test-compute_pwm.R index 56e3741..adb7fe8 100644 --- a/tests/testthat/test-compute_pwm.R +++ b/tests/testthat/test-compute_pwm.R @@ -58,3 +58,18 @@ test_that("compute_pwm works with 'max' func", { windows_log_sum_exp_r <- purrr::map_dbl(windows, compute_pwm, pssm, func = "logSumExp") expect_true(all(windows_r == windows_log_sum_exp_r)) }) + +test_that("compute_pwm score of a sequence is independent of batch companions", { + # A long sequence must score the same whether computed alone or in a batch + # with a shorter sequence. Regression test for the scan window being capped + # to nchar(sequences[1]) for every sequence. + long_seq <- substr(s, 1, 100) + short_seq <- substr(s, 1, 20) + for (f in c("max", "logSumExp")) { + alone <- compute_pwm(long_seq, pssm, func = f) + with_short_first <- compute_pwm(c(short_seq, long_seq), pssm, func = f)[2] + with_short_last <- compute_pwm(c(long_seq, short_seq), pssm, func = f)[1] + expect_equal(with_short_first, alone, info = f) + expect_equal(with_short_last, alone, info = f) + } +})