fix: score each sequence independently in compute_pwm

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).
 

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

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)
+    }
+})