feat: pwm.grad and pwm.grad.ism per-bp gradient vtracks - #111
Open
aviezerl wants to merge 15 commits into
Open
Conversation
Collaborator
Author
Code reviewFound 2 issues:
tests/testthat/test-vtrack-pwm-grad.R
|
Add API skeleton (Task 1 of the PWM gradient vtrack plan):
- R: .vtrack_params_pwm_grad validator (rejects score.thresh, defaults
aggregate to 'lse', validates 'lse'|'max'); register pwm.grad and
pwm.grad.ism in dispatch list and sourceless funcs; doc table rows.
- C++ PWMScorer: 4 new ScoringMode entries (GRAD_LSE, GRAD_MAX,
GRAD_LSE_ISM, GRAD_MAX_ISM); early-out stub in score_interval throws
rdb::verror until subsequent tasks implement them.
- C++ Track_var: PWM_GRAD / PWM_GRAD_ISM enum entries + matching
FUNC_NAMES strings; is_seq_variable, is_sequence_based_function and
is_pwm_function extended so SequenceVarProcessor classifies the new
vtracks under the PWM path and routes through pwm_scorer.
- C++ add_vtrack_var: dispatch parses pwm.grad / pwm.grad.ism (with
optional 'aggregate' rparam) and constructs PWMScorer with the right
ScoringMode.
End state: gvtrack.create('g', NULL, 'pwm.grad', pssm = ..., aggregate = 'lse')
parses, gvtrack.ls() shows the vtrack, validators reject score.thresh and
bad aggregate values, and gextract on a real DB throws
"PWMScorer: gradient modes not yet implemented" until Tasks 3-9 land.
Implements the linearized MAX-mode gradient (Task 3 of the pwm.grad plan) for the simplest case: bidirect=FALSE, strand=1, no spatial weighting. - Precompute m_worst_col0_fwd / m_worst_col_last_fwd at construction. - Route GRAD_MAX through a dedicated score_grad_max() that scans the iterator-clamped anchor range [i_min, i_max] (matching pwm.max's argmax convention) and returns M[0, b_p] - worst_col0 at the head anchor, 0 elsewhere. - Reject bidirect=TRUE / strand=-1 with verror BEFORE the try/catch in score_interval so the error reaches R instead of being squashed to NaN. - Skip the sliding-window dispatch for all four gradient modes; defer spatial weighting to Task 9 (NaN for now under m_use_spat). Tests verify both head=argmax (positive gradient matching the oracle) and head!=argmax (engine returns 0, oracle agrees), plus the two error paths.
Unifies the GRAD_MAX and GRAD_LSE answer paths under score_grad_linearized, using score_forward_original / score_reverse_original so the rc'd target (when strand_mode == -1) is canonicalized to fwd-genome semantics. Bidirect combines per-strand head scores via softmax. Tie-break in argmax favors the head anchor to match the oracle's R which.max convention.
Extends score_grad to handle ism=true: at the head anchor, only the head's per-anchor score changes under a base flip, so 3 alternative aggregates can be computed in O(1) using the pre-scanned best_no_head and lse_no_head. Covers both LSE and MAX aggregations and all strand modes (bidirect, fwd, rc). Renames score_grad_linearized to score_grad and unifies the four grad modes through one function. Tests: ISM-vs-linearized divergence (argmax shift on a periodic test seq), ISM bidirect on asymmetric PSSMs, ISM strand=-1, both LSE and MAX.
Per-anchor spat_log is added to each anchor's score in the scan; ISM uses the head's spat_log to keep the flipped per-anchor score on the same scale as best_no_head and lse_no_head. Tests: spat_factor=1 invariance and non-trivial spat with both lin and ism.
New vignettes/PWM-Functions.Rmd covering the existing PWM aggregating tracks (pwm, pwm.max, pwm.max.pos, pwm.count), edit-distance variants, and the new per-position attribution tracks pwm.grad / pwm.grad.ism. Math sections evaluate as pure-R demos (no genome DB needed); genome examples use eval=FALSE. Adds an articles: section to _pkgdown.yml listing the new vignette alongside the existing ones.
Cache per-anchor (fwd, rc) raw scores in a deque. On consecutive iterator steps with matching chrom/strand/window, slide the deque (pop trailing, score+push leading) instead of doing a full O(W*L) rescan. The aggregate is still recomputed from the deque per pivot (O(W)), so ISM_MAX which needs best_no_head doesn't require a second monotonic deque. Bench (HOMER.CTCF, mm10 chr1:3M, 100 kb, W=500, iter=1): pwm 0.11s (1.00x) pwm.grad 1.47s (13.25x; was 66.58x) Full O(1) running aggregates (RunningLogSumExp + dual RunningMaxDeque) would buy another ~5-10x but require a sliding-window 'max excluding the head' deque that's tricky for ISM_MAX.
- Add a dedicated 'Strands and bidirect' section that lays out the per-anchor score formulas, the bidirect=TRUE strand-union, and the position-reporting convention (pwm.max.pos sign, grad pivot at interval.start in fwd-genome coords). - Drop the 'Why LSE rather than just MAX?' subsection. - Trim gradient-directed framing in the intro and pwm/pwm.max sections. - Shorten the gradient section's motivation; one-line strand handling reference points to the new section.
The PWM is a linear log-likelihood model, so the 'gradient' is just a PSSM column lookup; integrated-gradients / DeepLIFT contribution rules don't apply (they collapse to gradient × input for linear models). ISM is standard in motif analysis (e.g. motifbreakR for SNPs). Renaming the user-facing description to 'per-bp PSSM-column contribution' (with explicit 'softmax-weighted' for LSE and 'argmax-conditioned' for MAX) is more accurate than 'gradient / saliency / DeepLIFT-style attribution'. The tracks themselves stay named pwm.grad / pwm.grad.ism.
PWM_GRAD / PWM_GRAD_ISM are anchored at seq_interval.start, so the multi-part filter aggregation that summed scores across disconnected unmasked parts mixed gradients at different genomic anchors. The single- part fallback had the same problem when the mask covered the pivot. Score only the unmasked part that starts at seq_interval.start; return NaN when the pivot is masked. Add filter tests for both pwm.grad and pwm.grad.ism (mask pivot, multi-part mask, non-intersecting mask).
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Two new sequence-based virtual tracks expose per-bp gradient / saliency of the PWM aggregate (LSE or MAX) at the iterator interval start:
pwm.grad— linearized gradient (DeepLIFT-style). For LSE:g(p) = w_p * (M[0, b_p] - min_b M[0, b])wherew_p = exp(score_p - f_LSE). For MAX: same diff if argmax is the head anchor, else 0.pwm.grad.ism— in-silico mutagenesis:g(p) = f(actual) - min_b' f(seq with seq[p] := b'). Computed in O(1) per pivot using pre-scannedbest_no_head/lse_no_head(only the head anchor's per-anchor score changes under a flip).Both support
aggregate = "lse" | "max", all strand modes (bidirect / fwd / rc), spatial weighting, prior, and extend. Pivot is fixed at the iterator interval start for v1.Spec:
dev/notes/features/2026-05-07_pwm-gradient-vtrack-design.mdPlan:
dev/notes/features/2026-05-07_pwm-gradient-vtrack-plan.mdTest plan
tests/testthat/helper-pwm-grad-oracle.R) for both linearized and ISM, validated on trivial inputs.tests/testthat/test-vtrack-pwm-grad.R).spat_factor,spat_bin) + invariance under all-1 spat.pwm,gscreen,gsummary.pwm,pwm.max,pwm.count, sliding-window paths unchanged.vignettes/PWM-Functions.Rmdbuilds (verified viarmarkdown::render; fulldevtools::build_vignettesto be exercised by CI).dev/benchmarks/2026-05-07_pwm-grad-bench.Rto be run on a real genome DB.Notes
spat_factorindex follows the rc'd-target order (matches existing engine convention; same forpwm.count).