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