From 9e2e3d8422a41c025463e958a94c8b983d33da7d Mon Sep 17 00:00:00 2001 From: aviezerl Date: Tue, 21 Jul 2026 19:31:25 +0300 Subject: [PATCH] fix: symmetric N/* handling across strands in PWM energy scoring compute_pwm() scored N (and the * wildcard) inconsistently between strands: integrate_energy() and integrate_energy_max() used the column's average log-probability (get_avg_log_prob) on the forward strand but a flat log(0.25) (c_log_quarter) on the reverse strand. With bidirect=TRUE this broke strand symmetry - a window and its reverse-complement scored differently (~1.6 nats) whenever an N fell on an informative position. Use get_avg_log_prob() on both strands, matching the forward branch and the other likelihood routines (max_like_match, integrate_like_seg). Only affects sequences containing N/*. Claude-Session: https://claude.ai/code/session_01PK3qefBGDoBwd9w5226FEq --- NEWS.md | 7 +++++++ src/DnaPSSM.cpp | 8 ++++---- tests/testthat/test-compute_pwm.R | 22 ++++++++++++++++++++++ 3 files changed, 33 insertions(+), 4 deletions(-) diff --git a/NEWS.md b/NEWS.md index 6c8813c..e03a4b3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,12 @@ # prego 0.0.10 +* Fix: `compute_pwm()` scored `N` (and the `*` wildcard) inconsistently between + strands - the forward strand used the column's average log-probability while + the reverse strand used a flat `log(0.25)`. With `bidirect = TRUE` this made a + sequence and its reverse-complement score differently whenever an `N` fell on + an informative position. Both strands now use the column average + (`get_avg_log_prob()`), matching the other likelihood routines, so scoring is + strand-symmetric again. Only affects sequences containing `N`/`*`. * Fix: `calc_seq_pwm()` / `extract_pwm()` now error clearly when given sequences of unequal length instead of silently recycling the shorter ones (via `rbind`) and returning wrong scores. These functions build a single rectangular one-hot diff --git a/src/DnaPSSM.cpp b/src/DnaPSSM.cpp index 1c09011..9afdbf9 100755 --- a/src/DnaPSSM.cpp +++ b/src/DnaPSSM.cpp @@ -872,10 +872,10 @@ void DnaPSSM::integrate_energy_max(const string &target, float &energy, vectorget_log_prob('C'); break; case '*': - logp_rev += c_log_quarter; + logp_rev += p->get_avg_log_prob(); break; case 'N': - logp_rev += c_log_quarter; + logp_rev += p->get_avg_log_prob(); break; default: break; @@ -947,10 +947,10 @@ void DnaPSSM::integrate_energy(const string &target, float &energy, vectorget_log_prob('C'); break; case '*': - logp += c_log_quarter; + logp += p->get_avg_log_prob(); break; case 'N': - logp += c_log_quarter; + logp += p->get_avg_log_prob(); break; default: break; diff --git a/tests/testthat/test-compute_pwm.R b/tests/testthat/test-compute_pwm.R index adb7fe8..b1164f6 100644 --- a/tests/testthat/test-compute_pwm.R +++ b/tests/testthat/test-compute_pwm.R @@ -59,6 +59,28 @@ test_that("compute_pwm works with 'max' func", { expect_true(all(windows_r == windows_log_sum_exp_r)) }) +test_that("compute_pwm handles N/* symmetrically on both strands", { + # For a bidirectional motif, a window and its reverse-complement must score + # identically. This held for plain ACGT but broke when the window contained + # an N at an informative position: the forward strand used the column's + # average log-prob while the reverse strand used a flat log(0.25). + win <- "AAATAAAAAAAAAAA" # single motif-length window (nchar == nrow(pssm) == 15) + win_n <- "AAANAAAAAAAAAAA" # N at position 4 (an informative, T-dominated column) + # tolerance is well above float32 traversal noise (~1e-6) but far below the + # ~1.6 nat gap the asymmetry produced. + for (f in c("max", "logSumExp")) { + expect_lt(abs( + compute_pwm(win_n, pssm, func = f, bidirect = TRUE) - + compute_pwm(rc(win_n), pssm, func = f, bidirect = TRUE) + ), 1e-4) + # sanity: the plain (N-free) window is strand-symmetric too + expect_lt(abs( + compute_pwm(win, pssm, func = f, bidirect = TRUE) - + compute_pwm(rc(win), pssm, func = f, bidirect = TRUE) + ), 1e-4) + } +}) + 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