From a42122cc5a9903a1170f8a93ed034fed49de4c66 Mon Sep 17 00:00:00 2001 From: aviezerl Date: Thu, 7 May 2026 22:29:25 +0300 Subject: [PATCH 01/15] feat(pwm.grad): plumbing for pwm.grad / pwm.grad.ism vtracks 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. --- R/vtrack.R | 35 +++++++++++++++++++++++++++++++++ man/gvtrack.create.Rd | 2 ++ src/PWMScorer.cpp | 14 +++++++++++++ src/PWMScorer.h | 6 +++++- src/TrackExpressionVars.cpp | 39 +++++++++++++++++++++++++++++-------- src/TrackExpressionVars.h | 8 +++++++- 6 files changed, 94 insertions(+), 10 deletions(-) diff --git a/R/vtrack.R b/R/vtrack.R index 36d628555..c55a3d6f1 100644 --- a/R/vtrack.R +++ b/R/vtrack.R @@ -246,6 +246,36 @@ kmer_params } +#' Validate and process PWM gradient function parameters +#' @noRd +.vtrack_params_pwm_grad <- function(func, params, dots) { + if (!is.null(params)) { + if (!is.list(params) || !("pssm" %in% names(params))) { + stop(sprintf("%s requires a list with at least 'pssm' matrix parameter", func)) + } + merged_dots <- params + } else { + merged_dots <- dots + } + + if (!is.null(merged_dots$score.thresh)) { + stop(sprintf("%s does not accept score.thresh", func)) + } + aggregate <- if (!is.null(merged_dots$aggregate)) merged_dots$aggregate else "lse" + if (!aggregate %in% c("lse", "max")) { + stop(sprintf( + "%s: aggregate must be 'lse' or 'max', got %s", + func, aggregate + )) + } + out <- .vtrack_params_pwm(func, params, dots) + out$aggregate <- aggregate + # gradient functions never take score.thresh; force-clear the default 0 + # introduced by .vtrack_params_pwm so it cannot be misinterpreted downstream. + out$score.thresh <- NULL + out +} + #' Validate and process PWM edit distance function parameters #' @noRd .vtrack_params_pwm_edit_distance <- function(func, params, dots) { @@ -560,6 +590,8 @@ pwm.max = .vtrack_params_pwm, pwm.max.pos = .vtrack_params_pwm, pwm.count = .vtrack_params_pwm, + pwm.grad = .vtrack_params_pwm_grad, + pwm.grad.ism = .vtrack_params_pwm_grad, kmer.count = .vtrack_params_kmer, kmer.frac = .vtrack_params_kmer, masked.count = .vtrack_params_masked, @@ -576,6 +608,7 @@ # Functions that don't require a source track .VTRACK_SOURCELESS_FUNCS <- c( "pwm", "pwm.max", "pwm.max.pos", "pwm.count", + "pwm.grad", "pwm.grad.ism", "kmer.count", "kmer.frac", "masked.count", "masked.frac", "pwm.edit_distance", "pwm.edit_distance.pos", "pwm.max.edit_distance", @@ -672,6 +705,8 @@ #' NULL (sequence) \tab pwm.max \tab pssm, bidirect, prior, extend, spat_* \tab Maximum log-likelihood score among all anchors (per-position union across strands). \cr #' NULL (sequence) \tab pwm.max.pos \tab pssm, bidirect, prior, extend, spat_* \tab 1-based position of the best-scoring anchor (signed by strand when \code{bidirect = TRUE}); coordinates are always relative to the iterator interval after any \code{gvtrack.iterator()} shifts/extensions. \cr #' NULL (sequence) \tab pwm.count \tab pssm, score.thresh, bidirect, prior, extend, strand, spat_* \tab Count of anchors whose score exceeds \code{score.thresh} (per-position union). \cr +#' NULL (sequence) \tab pwm.grad \tab pssm, aggregate, bidirect, prior, extend, spat_* \tab Linearized per-bp gradient (DeepLIFT-style) of the PWM aggregate at the interval start. \cr +#' NULL (sequence) \tab pwm.grad.ism \tab pssm, aggregate, bidirect, prior, extend, spat_* \tab In-silico mutagenesis per-bp gradient at the interval start. \cr #' } #' #' \strong{Edit distance summarizers} diff --git a/man/gvtrack.create.Rd b/man/gvtrack.create.Rd index 9b8404f0a..02a35d6d8 100644 --- a/man/gvtrack.create.Rd +++ b/man/gvtrack.create.Rd @@ -138,6 +138,8 @@ interval after all modifier adjustments. NULL (sequence) \tab pwm.max \tab pssm, bidirect, prior, extend, spat_* \tab Maximum log-likelihood score among all anchors (per-position union across strands). \cr NULL (sequence) \tab pwm.max.pos \tab pssm, bidirect, prior, extend, spat_* \tab 1-based position of the best-scoring anchor (signed by strand when \code{bidirect = TRUE}); coordinates are always relative to the iterator interval after any \code{gvtrack.iterator()} shifts/extensions. \cr NULL (sequence) \tab pwm.count \tab pssm, score.thresh, bidirect, prior, extend, strand, spat_* \tab Count of anchors whose score exceeds \code{score.thresh} (per-position union). \cr + NULL (sequence) \tab pwm.grad \tab pssm, aggregate, bidirect, prior, extend, spat_* \tab Linearized per-bp gradient (DeepLIFT-style) of the PWM aggregate at the interval start. \cr + NULL (sequence) \tab pwm.grad.ism \tab pssm, aggregate, bidirect, prior, extend, spat_* \tab In-silico mutagenesis per-bp gradient at the interval start. \cr } \strong{Edit distance summarizers} diff --git a/src/PWMScorer.cpp b/src/PWMScorer.cpp index 16efb0a19..90e5064ea 100644 --- a/src/PWMScorer.cpp +++ b/src/PWMScorer.cpp @@ -741,6 +741,20 @@ float PWMScorer::score_with_sliding_window(const std::string& target, float PWMScorer::score_interval(const GInterval& interval, const GenomeChromKey& chromkey) { + // Gradient modes are not yet implemented. Stub here so that any code path + // (score_without_spatial / score_with_spatial / sliding window / spatial sliding + // window) that would otherwise silently fall through to a NaN default produces + // a clear error instead. + switch (m_mode) { + case GRAD_LSE: + case GRAD_MAX: + case GRAD_LSE_ISM: + case GRAD_MAX_ISM: + rdb::verror("PWMScorer: gradient modes not yet implemented"); + default: + break; + } + // Calculate expanded interval to include full motif coverage int64_t motif_length = m_pssm.size(); GInterval expanded_interval = calculate_expanded_interval(interval, chromkey, motif_length); diff --git a/src/PWMScorer.h b/src/PWMScorer.h index 021247281..b7ae5b0b0 100644 --- a/src/PWMScorer.h +++ b/src/PWMScorer.h @@ -18,7 +18,11 @@ class PWMScorer : public GenomeSeqScorer TOTAL_LIKELIHOOD, // Integrated log-likelihood across all positions MAX_LIKELIHOOD, // Maximum log-likelihood score MAX_LIKELIHOOD_POS, // Position of maximum log-likelihood - MOTIF_COUNT // Count of positions exceeding threshold + MOTIF_COUNT, // Count of positions exceeding threshold + GRAD_LSE, // Linearized per-bp gradient (DeepLIFT-style) under LSE aggregate + GRAD_MAX, // Linearized per-bp gradient under MAX aggregate + GRAD_LSE_ISM, // In-silico mutagenesis per-bp gradient under LSE aggregate + GRAD_MAX_ISM // In-silico mutagenesis per-bp gradient under MAX aggregate }; PWMScorer(const DnaPSSM &pssm, const std::string &genome_root, bool extend = true, diff --git a/src/TrackExpressionVars.cpp b/src/TrackExpressionVars.cpp index 40e85e3fc..447a333d9 100644 --- a/src/TrackExpressionVars.cpp +++ b/src/TrackExpressionVars.cpp @@ -58,7 +58,9 @@ PWMEditDistanceScorer::Direction parse_direction_param(SEXP rparams, const std:: const char *TrackExpressionVars::Track_var::FUNC_NAMES[TrackExpressionVars::Track_var::NUM_FUNCS] = { "avg", "min", "max", "nearest", "stddev", "sum", "lse", "quantile", "global.percentile", "global.percentile.min", "global.percentile.max", - "weighted.sum", "area", "pwm", "pwm.max", "pwm.max.pos", "pwm.count", "kmer.count", "kmer.frac", + "weighted.sum", "area", "pwm", "pwm.max", "pwm.max.pos", "pwm.count", + "pwm.grad", "pwm.grad.ism", + "kmer.count", "kmer.frac", "masked.count", "masked.frac", "max.pos.abs", "max.pos.relative", "min.pos.abs", "min.pos.relative", "exists", "size", "sample", "sample.pos.abs", "sample.pos.relative", @@ -455,7 +457,8 @@ void TrackExpressionVars::add_vtrack_var(const string &vtrack, SEXP rvtrack) string func = CHAR(STRING_ELT(rfunc, 0)); transform(func.begin(), func.end(), func.begin(), ::tolower); - if (func == "pwm" || func == "pwm.max" || func == "pwm.max.pos" || func == "pwm.count") { + if (func == "pwm" || func == "pwm.max" || func == "pwm.max.pos" || func == "pwm.count" || + func == "pwm.grad" || func == "pwm.grad.ism") { // Create the Track_var without a Track_n_imdf m_track_vars.push_back(Track_var()); Track_var &var = m_track_vars.back(); @@ -463,24 +466,44 @@ void TrackExpressionVars::add_vtrack_var(const string &vtrack, SEXP rvtrack) var.val_func = (func == "pwm" ? Track_var::PWM : func == "pwm.max" ? Track_var::PWM_MAX : func == "pwm.max.pos" ? Track_var::PWM_MAX_POS : - Track_var::PWM_COUNT); + func == "pwm.count" ? Track_var::PWM_COUNT : + func == "pwm.grad" ? Track_var::PWM_GRAD : + Track_var::PWM_GRAD_ISM); var.track_n_imdf = nullptr; // No track needed for PWM var.seq_imdf1d = nullptr; - + SEXP rparams = get_rvector_col(rvtrack, "params", vtrack.c_str(), false); // Parse PWM parameters using helper struct TrackExprParams::PWMParams pwm_params = TrackExprParams::PWMParams::parse(rparams, vtrack); + // Gradient functions read aggregate ("lse" or "max") from params + std::string aggregate = "lse"; + if (func == "pwm.grad" || func == "pwm.grad.ism") { + SEXP r_agg = get_rvector_col(rparams, "aggregate", vtrack.c_str(), false); + if (!Rf_isNull(r_agg) && Rf_isString(r_agg)) { + aggregate = CHAR(STRING_ELT(r_agg, 0)); + } + if (aggregate != "lse" && aggregate != "max") { + verror("Virtual track %s: aggregate must be 'lse' or 'max', got '%s'", + vtrack.c_str(), aggregate.c_str()); + } + } + + PWMScorer::ScoringMode mode = + func == "pwm" ? PWMScorer::TOTAL_LIKELIHOOD : + func == "pwm.max" ? PWMScorer::MAX_LIKELIHOOD : + func == "pwm.max.pos" ? PWMScorer::MAX_LIKELIHOOD_POS : + func == "pwm.count" ? PWMScorer::MOTIF_COUNT : + func == "pwm.grad" ? (aggregate == "lse" ? PWMScorer::GRAD_LSE : PWMScorer::GRAD_MAX) : + (aggregate == "lse" ? PWMScorer::GRAD_LSE_ISM : PWMScorer::GRAD_MAX_ISM); + // Construct scorer with shared sequence fetcher for caching var.pwm_scorer = std::make_unique( pwm_params.core.pssm, &m_shared_seqfetch, pwm_params.extend_flag, - func == "pwm" ? PWMScorer::TOTAL_LIKELIHOOD : - func == "pwm.max" ? PWMScorer::MAX_LIKELIHOOD : - func == "pwm.max.pos" ? PWMScorer::MAX_LIKELIHOOD_POS : - PWMScorer::MOTIF_COUNT, + mode, static_cast(pwm_params.core.strand_mode), pwm_params.core.spat_factor, pwm_params.core.spat_bin_size, diff --git a/src/TrackExpressionVars.h b/src/TrackExpressionVars.h index b59a529db..5448acce3 100644 --- a/src/TrackExpressionVars.h +++ b/src/TrackExpressionVars.h @@ -133,6 +133,8 @@ class TrackExpressionVars { PWM_MAX, PWM_MAX_POS, PWM_COUNT, + PWM_GRAD, + PWM_GRAD_ISM, KMER_COUNT, KMER_FRAC, MASKED_COUNT, @@ -429,6 +431,8 @@ inline bool TrackExpressionVars::is_seq_variable(unsigned ivar) const { m_track_vars[ivar].val_func == Track_var::PWM_MAX || m_track_vars[ivar].val_func == Track_var::PWM_MAX_POS || m_track_vars[ivar].val_func == Track_var::PWM_COUNT || + m_track_vars[ivar].val_func == Track_var::PWM_GRAD || + m_track_vars[ivar].val_func == Track_var::PWM_GRAD_ISM || m_track_vars[ivar].val_func == Track_var::KMER_COUNT || m_track_vars[ivar].val_func == Track_var::KMER_FRAC || m_track_vars[ivar].val_func == Track_var::MASKED_COUNT || @@ -445,6 +449,7 @@ inline bool TrackExpressionVars::is_seq_variable(unsigned ivar) const { inline bool TrackExpressionVars::is_sequence_based_function(Track_var::Val_func func) { return func == Track_var::PWM || func == Track_var::PWM_MAX || func == Track_var::PWM_MAX_POS || func == Track_var::PWM_COUNT || + func == Track_var::PWM_GRAD || func == Track_var::PWM_GRAD_ISM || func == Track_var::KMER_COUNT || func == Track_var::KMER_FRAC || func == Track_var::MASKED_COUNT || func == Track_var::MASKED_FRAC || func == Track_var::PWM_EDIT_DISTANCE || @@ -457,7 +462,8 @@ inline bool TrackExpressionVars::is_sequence_based_function(Track_var::Val_func inline bool TrackExpressionVars::is_pwm_function(Track_var::Val_func func) { return func == Track_var::PWM || func == Track_var::PWM_MAX || - func == Track_var::PWM_MAX_POS || func == Track_var::PWM_COUNT; + func == Track_var::PWM_MAX_POS || func == Track_var::PWM_COUNT || + func == Track_var::PWM_GRAD || func == Track_var::PWM_GRAD_ISM; } inline bool TrackExpressionVars::is_pwm_edit_distance_function(Track_var::Val_func func) { From 85e9529519c7ac9c7ced3a7eeb3c276934f7dc03 Mon Sep 17 00:00:00 2001 From: aviezerl Date: Thu, 7 May 2026 22:37:35 +0300 Subject: [PATCH 02/15] fix(pwm.grad): require aggregate to be length-1 string --- src/TrackExpressionVars.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/TrackExpressionVars.cpp b/src/TrackExpressionVars.cpp index 447a333d9..4fe5c9ce1 100644 --- a/src/TrackExpressionVars.cpp +++ b/src/TrackExpressionVars.cpp @@ -481,7 +481,11 @@ void TrackExpressionVars::add_vtrack_var(const string &vtrack, SEXP rvtrack) std::string aggregate = "lse"; if (func == "pwm.grad" || func == "pwm.grad.ism") { SEXP r_agg = get_rvector_col(rparams, "aggregate", vtrack.c_str(), false); - if (!Rf_isNull(r_agg) && Rf_isString(r_agg)) { + if (!Rf_isNull(r_agg)) { + if (!Rf_isString(r_agg) || Rf_length(r_agg) != 1) { + verror("Virtual track %s: aggregate must be a single string", + vtrack.c_str()); + } aggregate = CHAR(STRING_ELT(r_agg, 0)); } if (aggregate != "lse" && aggregate != "max") { From 175a1cce0ff981f596d05211fabfda4417928b60 Mon Sep 17 00:00:00 2001 From: aviezerl Date: Thu, 7 May 2026 22:41:51 +0300 Subject: [PATCH 03/15] test(pwm.grad): R oracle for linearized + ISM gradients --- tests/testthat/helper-pwm-grad-oracle.R | 181 ++++++++++++++++++++++++ tests/testthat/test-vtrack-pwm-grad.R | 143 +++++++++++++++++++ 2 files changed, 324 insertions(+) create mode 100644 tests/testthat/helper-pwm-grad-oracle.R create mode 100644 tests/testthat/test-vtrack-pwm-grad.R diff --git a/tests/testthat/helper-pwm-grad-oracle.R b/tests/testthat/helper-pwm-grad-oracle.R new file mode 100644 index 000000000..0de6904e3 --- /dev/null +++ b/tests/testthat/helper-pwm-grad-oracle.R @@ -0,0 +1,181 @@ +# Oracle helpers for PWM gradient vtrack tests. +# +# These are pure R reference implementations for the linearized gradient (B) +# and ISM gradient (A) defined in the PWM gradient vtrack design. +# They are intentionally simple (correctness over speed) and reuse the same +# prior + per-row normalization as `manual_pwm_scores_single_strand` in +# helper-pwm.R, so engine outputs can be validated against these in later +# tasks. + +# Map base character -> 1-based index (A=1, C=2, G=3, T=4). NA for non-ACGT. +.base_to_idx <- function(base) { + switch(base, + "A" = 1L, + "C" = 2L, + "G" = 3L, + "T" = 4L, + NA_integer_ + ) +} + +# Complement of an index (1<->4, 2<->3). +.complement_idx <- function(idx) c(4L, 3L, 2L, 1L)[idx] + +# Reverse-complement of a sequence string. +.rc_seq <- function(seq) { + comp <- chartr("ACGT", "TGCA", seq) + paste(rev(strsplit(comp, "")[[1]]), collapse = "") +} + +# Apply prior + per-row normalization, then take log. Matches the path used in +# manual_pwm_scores_single_strand (same prior treatment). +.normalize_pssm <- function(pssm, prior = 0.01) { + if (prior > 0) { + normalized <- matrix(0, nrow = nrow(pssm), ncol = ncol(pssm)) + for (i in seq_len(nrow(pssm))) { + row <- pssm[i, ] + prior + normalized[i, ] <- row / sum(row) + } + pssm <- normalized + } + pssm +} + +# Per-anchor scores in a window. +# $fwd[a] : log-prob of fwd PSSM applied to seq[a..a+L-1] +# $rc[a] : log-prob of fwd PSSM applied to RC of seq[a..a+L-1] +# (== rc PSSM applied to seq[a..a+L-1]) +# $comb[a]: log(exp(fwd[a]) + exp(rc[a])) +.oracle_anchor_scores <- function(seq, pssm, prior = 0.01) { + L <- nrow(pssm) + n <- nchar(seq) + stopifnot(n >= L) + + # Fwd uses the same code path as manual_pwm_scores_single_strand. + fwd <- manual_pwm_scores_single_strand(seq, pssm, prior = prior) + + # Rc: apply fwd PSSM to RC of each window. Easiest: build rc PSSM + # = pssm[L:1, c(4,3,2,1)] and run on the original seq. This is + # equivalent (since prior + normalization is row-wise) to running the + # fwd PSSM on the RC of the seq. + rc_pssm <- pssm[L:1, c(4L, 3L, 2L, 1L), drop = FALSE] + if (!is.null(colnames(pssm))) colnames(rc_pssm) <- colnames(pssm) + rc <- manual_pwm_scores_single_strand(seq, rc_pssm, prior = prior) + + comb <- vapply(seq_along(fwd), function(i) { + log_sum_exp(c(fwd[i], rc[i])) + }, numeric(1)) + + list(fwd = fwd, rc = rc, comb = comb) +} + +# Aggregate a vector of per-anchor scores via "lse" or "max". +.aggregate_scores <- function(s, aggregate) { + aggregate <- match.arg(aggregate, c("lse", "max")) + if (aggregate == "max") { + return(max(s)) + } + log_sum_exp(s) +} + +# Pick the per-anchor score vector to use given (bidirect, strand). +# bidirect = TRUE -> comb +# bidirect = FALSE, strand = +1 -> fwd +# bidirect = FALSE, strand = -1 -> rc +.pick_scores <- function(scores, bidirect, strand) { + if (bidirect) { + return(scores$comb) + } + if (strand == 1L || strand == 1) { + return(scores$fwd) + } + if (strand == -1L || strand == -1) { + return(scores$rc) + } + stop("strand must be +1 or -1 when bidirect = FALSE") +} + +# Linearized gradient at p = interval start, head anchor, column 0 (= base b_p). +# +# Returns a single non-negative scalar. For an all-N seq or otherwise invalid +# head base, returns 0. +manual_pwm_grad <- function(pssm, seq, aggregate = c("lse", "max"), + bidirect = TRUE, strand = 1L, prior = 0.01) { + aggregate <- match.arg(aggregate) + L <- nrow(pssm) + + # Normalized log PSSM (used for diff lookups). + M <- log(.normalize_pssm(pssm, prior = prior)) + + # Per-anchor scores. + scores <- .oracle_anchor_scores(seq, pssm, prior = prior) + s <- .pick_scores(scores, bidirect, strand) + f <- .aggregate_scores(s, aggregate) + + # Head base. + b_p_char <- substr(seq, 1L, 1L) + b_p <- .base_to_idx(b_p_char) + if (is.na(b_p)) { + return(0) + } + + # Diffs at column 0 of head anchor for each strand. + diff_fwd <- M[1L, b_p] - min(M[1L, ]) + diff_rc <- M[L, .complement_idx(b_p)] - min(M[L, ]) + + if (!bidirect) { + diff <- if (strand == 1L || strand == 1) diff_fwd else diff_rc + if (aggregate == "max") { + # 1-based index of head anchor is 1. + return(if (which.max(s) == 1L) diff else 0) + } + # LSE: weight by softmax of head anchor. + w_p <- exp(s[1L] - f) + return(w_p * diff) + } + + # Bidirect: per-strand softmax at the head anchor, combined diff. + fwd_p <- scores$fwd[1L] + rc_p <- scores$rc[1L] + comb_p <- scores$comb[1L] + sm_fwd <- exp(fwd_p - comb_p) + sm_rc <- 1 - sm_fwd + combined_diff <- sm_fwd * diff_fwd + sm_rc * diff_rc + + if (aggregate == "max") { + return(if (which.max(s) == 1L) combined_diff else 0) + } + w_p_comb <- exp(comb_p - f) + w_p_comb * combined_diff +} + +# In-silico mutagenesis at p = interval start. +# g = f_actual - min over b' != b_p of f(seq with seq[1] := b'). +# If b_p is not ACGT, returns 0. +manual_pwm_grad_ism <- function(pssm, seq, aggregate = c("lse", "max"), + bidirect = TRUE, strand = 1L, prior = 0.01) { + aggregate <- match.arg(aggregate) + + aggregate_for_seq <- function(s_seq) { + sc <- .oracle_anchor_scores(s_seq, pssm, prior = prior) + s <- .pick_scores(sc, bidirect, strand) + .aggregate_scores(s, aggregate) + } + + f_actual <- aggregate_for_seq(seq) + + b_p_char <- substr(seq, 1L, 1L) + b_p <- .base_to_idx(b_p_char) + if (is.na(b_p)) { + return(0) + } + + bases <- c("A", "C", "G", "T") + others <- bases[-b_p] + f_alt <- vapply(others, function(b) { + modified <- paste0(b, substr(seq, 2L, nchar(seq))) + aggregate_for_seq(modified) + }, numeric(1)) + + f_actual - min(f_alt) +} diff --git a/tests/testthat/test-vtrack-pwm-grad.R b/tests/testthat/test-vtrack-pwm-grad.R new file mode 100644 index 000000000..93e785dee --- /dev/null +++ b/tests/testthat/test-vtrack-pwm-grad.R @@ -0,0 +1,143 @@ +test_that("oracle: flat PSSM gives g = 0 (lse, max, ism)", { + pssm <- matrix(0.25, 3, 4, dimnames = list(NULL, c("A", "C", "G", "T"))) + expect_equal(manual_pwm_grad(pssm, "ACGTAA", "lse", bidirect = TRUE), 0) + expect_equal(manual_pwm_grad(pssm, "ACGTAA", "max", bidirect = TRUE), 0) + expect_equal(manual_pwm_grad_ism(pssm, "ACGTAA", "lse", bidirect = TRUE), 0) +}) + +test_that("oracle: L=1 PSSM, single strand, g_MAX = M[1,b_p] - min(M[1,])", { + pssm <- matrix(c(0.7, 0.1, 0.1, 0.1), 1, 4, + byrow = TRUE, + dimnames = list(NULL, c("A", "C", "G", "T")) + ) + M <- log((pssm + 0.01) / sum(pssm + 0.01)) + expected <- unname(M[1, 1] - min(M[1, ])) + expect_equal( + manual_pwm_grad(pssm, "ACGT", "max", bidirect = FALSE, strand = 1L), + expected, + tolerance = 1e-10 + ) +}) + +test_that("oracle: linearized matches ISM when argmax is stable under flips (single strand)", { + pssm <- matrix( + c( + 0.97, 0.01, 0.01, 0.01, + 0.01, 0.97, 0.01, 0.01, + 0.01, 0.01, 0.97, 0.01 + ), + 3, 4, + byrow = TRUE, + dimnames = list(NULL, c("A", "C", "G", "T")) + ) + seq <- "ACGAAAAA" # head anchor "ACG" perfect, others weak. + g_lin <- manual_pwm_grad(pssm, seq, "max", bidirect = FALSE, strand = 1L) + g_ism <- manual_pwm_grad_ism(pssm, seq, "max", bidirect = FALSE, strand = 1L) + expect_equal(g_lin, g_ism, tolerance = 1e-8) +}) + +test_that("oracle: linearized lse with concentrated softmax matches max", { + # Longer PSSM + uniform A-tail in seq pushes the non-head anchor scores + # so far below the head anchor that the softmax concentrates at >= 1-1e-4 + # on the head, making g_LSE numerically indistinguishable from g_MAX. + pssm <- matrix( + c( + 0.97, 0.01, 0.01, 0.01, + 0.01, 0.97, 0.01, 0.01, + 0.01, 0.01, 0.97, 0.01, + 0.01, 0.01, 0.97, 0.01, + 0.01, 0.01, 0.97, 0.01 + ), + 5, 4, + byrow = TRUE, + dimnames = list(NULL, c("A", "C", "G", "T")) + ) + seq <- "ACGGGAAAAAAAA" # head anchor "ACGGG" perfect, others weak. + g_max <- manual_pwm_grad(pssm, seq, "max", bidirect = FALSE, strand = 1L) + g_lse <- manual_pwm_grad(pssm, seq, "lse", bidirect = FALSE, strand = 1L) + expect_equal(g_max, g_lse, tolerance = 1e-3) +}) + +test_that("oracle: bidirect on palindromic PSSM matches single-strand on palindrome", { + pssm <- matrix( + c( + 0.7, 0.1, 0.1, 0.1, + 0.1, 0.4, 0.4, 0.1, + 0.1, 0.1, 0.1, 0.7 + ), + 3, 4, + byrow = TRUE, + dimnames = list(NULL, c("A", "C", "G", "T")) + ) + rc <- pssm[3:1, c(4, 3, 2, 1)] + dimnames(rc) <- dimnames(pssm) + expect_equal(rc, pssm) + g <- manual_pwm_grad(pssm, "ATCAT", "lse", bidirect = TRUE) + expect_true(is.finite(g) && g >= 0) +}) + +test_that("oracle: rc by RC-of-seq equals rc by RC-of-PSSM", { + pssm <- matrix( + c( + 0.7, 0.1, 0.1, 0.1, + 0.1, 0.6, 0.2, 0.1, + 0.1, 0.1, 0.6, 0.2 + ), + 3, 4, + byrow = TRUE, + dimnames = list(NULL, c("A", "C", "G", "T")) + ) + seq <- "ACGCAT" + g_strand_minus <- manual_pwm_grad(pssm, seq, "max", bidirect = FALSE, strand = -1L) + expect_true(is.finite(g_strand_minus) && g_strand_minus >= 0) +}) + +test_that("oracle: rc-of-seq vs rc-of-PSSM give same per-anchor scores", { + # Self-consistency: scoring fwd PSSM on RC(seq) at anchor a should equal + # scoring rc PSSM on seq at the mirrored anchor (n - L - a + 1). + pssm <- matrix( + c( + 0.6, 0.2, 0.1, 0.1, + 0.1, 0.5, 0.3, 0.1, + 0.2, 0.1, 0.5, 0.2 + ), + 3, 4, + byrow = TRUE, + dimnames = list(NULL, c("A", "C", "G", "T")) + ) + seq <- "ACGTACGTAC" + L <- nrow(pssm) + n <- nchar(seq) + + sc <- .oracle_anchor_scores(seq, pssm, prior = 0.01) + + # Independent path: rc score at anchor a == fwd score on RC(seq) at + # position (n - L - (a - 1) + 1) - 1 + 1 = n - L - a + 2 (1-based). + rc_seq <- .rc_seq(seq) + fwd_on_rc <- manual_pwm_scores_single_strand(rc_seq, pssm, prior = 0.01) + + n_anchors <- n - L + 1L + for (a in seq_len(n_anchors)) { + mirror <- n_anchors - a + 1L + expect_equal(sc$rc[a], fwd_on_rc[mirror], tolerance = 1e-12) + } +}) + +test_that("oracle: ISM is always >= 0", { + pssm <- matrix( + c( + 0.7, 0.1, 0.1, 0.1, + 0.1, 0.6, 0.2, 0.1, + 0.1, 0.1, 0.6, 0.2 + ), + 3, 4, + byrow = TRUE, + dimnames = list(NULL, c("A", "C", "G", "T")) + ) + for (seq in c("ACGCAT", "TTTTTT", "ACGACG", "GGGGGG")) { + for (agg in c("lse", "max")) { + g <- manual_pwm_grad_ism(pssm, seq, agg, bidirect = TRUE) + expect_gte(g, -1e-10) + } + } +}) From ca04082409a2941cc362004e73ce6c3c9b87ac54 Mon Sep 17 00:00:00 2001 From: aviezerl Date: Thu, 7 May 2026 22:55:20 +0300 Subject: [PATCH 04/15] feat(pwm.grad): GRAD_MAX single-strand fwd implementation 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. --- src/PWMScorer.cpp | 154 ++++++++++++++++++++++++-- src/PWMScorer.h | 10 ++ tests/testthat/test-vtrack-pwm-grad.R | 134 ++++++++++++++++++++++ 3 files changed, 290 insertions(+), 8 deletions(-) diff --git a/src/PWMScorer.cpp b/src/PWMScorer.cpp index 90e5064ea..b38efef58 100644 --- a/src/PWMScorer.cpp +++ b/src/PWMScorer.cpp @@ -108,6 +108,31 @@ static inline float pos_value_with_dir(const DnaPSSM& pssm, return best_val; } +// Precompute the worst-base log-prob at the head and tail PSSM columns. These +// are used as the "neutral" baseline for the linearized gradient +// (g = M[col, b_p] - min_b M[col, b]). Computed once at construction so the +// per-call gradient path is O(1) on top of the existing scan. +static inline void compute_worst_col_logprobs(const DnaPSSM& pssm, + float& worst_col0, + float& worst_col_last) +{ + if (pssm.size() == 0) { + worst_col0 = 0.0f; + worst_col_last = 0.0f; + return; + } + const DnaProbVec& col0 = pssm[0]; + worst_col0 = std::min({col0.get_log_prob_from_code(0), + col0.get_log_prob_from_code(1), + col0.get_log_prob_from_code(2), + col0.get_log_prob_from_code(3)}); + const DnaProbVec& col_last = pssm[pssm.size() - 1]; + worst_col_last = std::min({col_last.get_log_prob_from_code(0), + col_last.get_log_prob_from_code(1), + col_last.get_log_prob_from_code(2), + col_last.get_log_prob_from_code(3)}); +} + PWMScorer::PWMScorer(const DnaPSSM& pssm, const std::string& genome_root, bool extend, ScoringMode mode, char strand, const std::vector& spat_factor, int spat_bin_size, float score_thresh) @@ -123,6 +148,8 @@ PWMScorer::PWMScorer(const DnaPSSM& pssm, const std::string& genome_root, bool e m_spat_log_factors[i] = std::log(std::max(1e-30f, spat_factor[i])); } } + + compute_worst_col_logprobs(m_pssm, m_worst_col0_fwd, m_worst_col_last_fwd); } PWMScorer::PWMScorer(const DnaPSSM& pssm, GenomeSeqFetch* shared_seqfetch, bool extend, @@ -140,6 +167,8 @@ PWMScorer::PWMScorer(const DnaPSSM& pssm, GenomeSeqFetch* shared_seqfetch, bool m_spat_log_factors[i] = std::log(std::max(1e-30f, spat_factor[i])); } } + + compute_worst_col_logprobs(m_pssm, m_worst_col0_fwd, m_worst_col_last_fwd); } void PWMScorer::invalidate_cache() @@ -337,6 +366,82 @@ float PWMScorer::score_without_spatial(const std::string& target, int64_t motif_ return compute_position_result(pos_idx, target.length(), motif_length, best_dir); } +// Linearized per-bp gradient (DeepLIFT-style) under MAX aggregation, fwd-only. +// +// Pivot is fixed at the iterator interval start (target[0] == genomic +// interval.start). The argmax is taken over anchors with starts in +// [i_min, i_max] inside the fetched/expanded target, which mirrors the clamping +// applied to pwm.max so the two functions agree on which anchor "wins". When +// the head anchor (i = 0) is the argmax, the gradient is +// M[0, b_p] - min_b M[0, b] +// where M[0, .] are the (prior-adjusted) log-probs at column 0 of the PSSM and +// b_p is the head base. Otherwise the gradient is 0. +// +// Restrictions for Task 3: single-strand fwd only. bidirect=TRUE (m_pssm +// bidirect) and strand=-1 are deferred to a later task (full bidirect / rc +// formula needs both fwd and rc head scores at the head anchor). +float PWMScorer::score_grad_max(const std::string& target, + size_t i_min, size_t i_max, + size_t motif_length) +{ + // Precondition: bidirect/strand restrictions have already been enforced + // by the dispatch in score_interval (those verrors live outside the + // try/catch block so they reach R; calling verror here would be silently + // converted to NaN by the surrounding handler). + + // Iterator-interval was empty after clamping, or interval shorter than motif. + if (i_min > i_max) { + return std::numeric_limits::quiet_NaN(); + } + if (target.size() < (size_t)motif_length) { + return std::numeric_limits::quiet_NaN(); + } + if (i_max + (size_t)motif_length > target.size()) { + return std::numeric_limits::quiet_NaN(); + } + + // The pivot ("p" in the spec) is the iterator interval start, i.e. target[0]. + // It is the *head anchor* only when i_min == 0; otherwise the head anchor + // sits outside the iterator-interval scan and the gradient at p is 0 by + // definition (column 0 of every scanned anchor uses some other base). + if (i_min > 0) { + return 0.0f; + } + + // Scan fwd-only over [i_min, i_max] and find the argmax. + float best = -std::numeric_limits::infinity(); + size_t argmax = i_min; + bool any_finite = false; + for (size_t i = i_min; i <= i_max; ++i) { + float v = 0.0f; + std::string::const_iterator it = target.begin() + i; + m_pssm.calc_like(it, v); + if (std::isfinite(v) && v > best) { + best = v; + argmax = i; + any_finite = true; + } + } + + if (!any_finite) { + // All anchors hit non-ACGT bases (-Inf): cannot define a gradient. + return std::numeric_limits::quiet_NaN(); + } + + // Head base (target[0]). N or other non-ACGT -> NA. + int b_p = DnaLookupTables::BASE_ENCODE[(unsigned char)target[0]]; + if (b_p < 0) { + return std::numeric_limits::quiet_NaN(); + } + + if (argmax != 0) { + return 0.0f; + } + + // Argmax is at the head anchor: gradient = M[0, b_p] - worst_col0. + return m_pssm[0].get_log_prob_from_code(b_p) - m_worst_col0_fwd; +} + // Score with spatial weighting float PWMScorer::score_with_spatial(const std::string& target, int64_t motif_length) { @@ -741,16 +846,28 @@ float PWMScorer::score_with_sliding_window(const std::string& target, float PWMScorer::score_interval(const GInterval& interval, const GenomeChromKey& chromkey) { - // Gradient modes are not yet implemented. Stub here so that any code path - // (score_without_spatial / score_with_spatial / sliding window / spatial sliding - // window) that would otherwise silently fall through to a NaN default produces - // a clear error instead. + // Gradient modes: GRAD_MAX (single-strand fwd) is implemented below via a + // dedicated answer path. The remaining gradient modes are not yet wired up + // and must produce a clear error instead of falling through into the + // generic scan paths. + // + // NOTE: these checks are intentionally placed BEFORE the try/catch block + // below so that verror surfaces as an R-level error rather than being + // swallowed and converted to NaN by the TGLException handler. switch (m_mode) { case GRAD_LSE: - case GRAD_MAX: case GRAD_LSE_ISM: case GRAD_MAX_ISM: - rdb::verror("PWMScorer: gradient modes not yet implemented"); + rdb::verror("PWMScorer: gradient mode not yet implemented"); + case GRAD_MAX: + // Task 3 supports only single-strand forward; reject bidirect + // PSSMs and strand=-1 here so the caller sees an error rather than + // a silent NaN. + if (m_pssm.is_bidirect() || m_strand != 1) { + rdb::verror("pwm.grad currently supports only bidirect=FALSE, strand=1; " + "bidirect=TRUE and strand=-1 will be added in a future task"); + } + break; default: break; } @@ -817,6 +934,21 @@ float PWMScorer::score_interval(const GInterval& interval, const GenomeChromKey& i_min = 0; } + // Spatial weighting for gradient modes is deferred to a later task. + // Surface a clear error rather than silently using non-spatial code. + if (m_use_spat && + (m_mode == GRAD_LSE || m_mode == GRAD_MAX || + m_mode == GRAD_LSE_ISM || m_mode == GRAD_MAX_ISM)) { + rdb::verror("PWMScorer: spatial weighting not yet supported for gradient modes"); + } + + // Gradient modes (single-strand fwd) take a dedicated answer path. + // They share the iterator clamping (i_min..i_max) above so that the + // argmax convention matches pwm.max exactly. + if (!m_use_spat && m_mode == GRAD_MAX) { + return score_grad_max(target, i_min, i_max, motif_len); + } + // Try spatial sliding window optimization if (m_use_spat) { size_t stride = 0; @@ -856,8 +988,14 @@ float PWMScorer::score_interval(const GInterval& interval, const GenomeChromKey& return score_with_spatial(target, motif_length); } - // Non-spatial sliding window optimization (original code) - if (!m_use_spat && m_mode != MAX_LIKELIHOOD_POS) { + // Non-spatial sliding window optimization (original code). + // Gradient modes are not supported by the sliding window machinery + // (they need a per-call full scan), so force them through the + // dedicated gradient answer paths above (or, if not yet + // implemented, the verror dispatch at the top of score_interval). + if (!m_use_spat && m_mode != MAX_LIKELIHOOD_POS && + m_mode != GRAD_LSE && m_mode != GRAD_MAX && + m_mode != GRAD_LSE_ISM && m_mode != GRAD_MAX_ISM) { return score_with_sliding_window(target, interval, expanded_interval, i_min, i_max, motif_len); } } diff --git a/src/PWMScorer.h b/src/PWMScorer.h index b7ae5b0b0..71852067a 100644 --- a/src/PWMScorer.h +++ b/src/PWMScorer.h @@ -171,11 +171,21 @@ class PWMScorer : public GenomeSeqScorer // Utilities inline float get_spatial_log_factor(size_t pos_index) const; + // Gradient-mode answer paths (single-strand fwd; bidirect / rc deferred) + float score_grad_max(const std::string& target, size_t i_min, size_t i_max, + size_t motif_length); + // Core members DnaPSSM m_pssm; ScoringMode m_mode; float m_score_thresh = 0.0f; + // Worst per-base log-prob at the head and tail PSSM columns. Used as the + // baseline for the linearized gradient ("actual minus worst" per spec). + // Computed once at construction. + float m_worst_col0_fwd = 0.0f; + float m_worst_col_last_fwd = 0.0f; + // Spatial weighting bool m_use_spat = false; std::vector m_spat_log_factors; diff --git a/tests/testthat/test-vtrack-pwm-grad.R b/tests/testthat/test-vtrack-pwm-grad.R index 93e785dee..8be9b6318 100644 --- a/tests/testthat/test-vtrack-pwm-grad.R +++ b/tests/testthat/test-vtrack-pwm-grad.R @@ -141,3 +141,137 @@ test_that("oracle: ISM is always >= 0", { } } }) + +# ---- Engine tests for GRAD_MAX (Task 3) ---- + +create_isolated_test_db() + +# Build a strong-consensus PSSM (0.97/0.01/0.01/0.01) whose column j has its +# peak at the j-th base of `bases`. `bases` must be a length-L character vector +# of A/C/G/T (no Ns). +.consensus_pssm <- function(bases) { + L <- length(bases) + stopifnot(all(bases %in% c("A", "C", "G", "T"))) + pssm <- matrix(0.01, + nrow = L, ncol = 4, + dimnames = list(NULL, c("A", "C", "G", "T")) + ) + for (j in seq_len(L)) { + idx <- match(bases[j], c("A", "C", "G", "T")) + pssm[j, idx] <- 0.97 + } + pssm +} + +test_that("pwm.grad max: head=argmax => engine matches oracle col-0 diff", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 20L # number of iterator-interval positions per pivot + # Use a 1bp iterator interval; the iterator extension makes the vtrack + # see a window of length iter_W (sshift=0, eshift=iter_W-1). + pivot <- gintervals(1, 200, 201) + # The engine fetches sequence over [pivot.start + sshift, + # pivot.end + eshift + (L-1)) + # because extend=TRUE adds (L-1) at the end. That's [200, 200 + iter_W + L - 1). + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + # Pick PSSM consensus = first L bases of seq_ext, so the head anchor wins. + head_bases <- strsplit(substr(seq_ext, 1, L), "")[[1]] + pssm <- .consensus_pssm(head_bases) + + gvtrack.create("g", NULL, "pwm.grad", + pssm = pssm, + aggregate = "max", bidirect = FALSE, strand = 1, + extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g", sshift = 0, eshift = iter_W - 1) + + res <- gextract("g", iterator = 1, intervals = pivot) + + expected <- manual_pwm_grad(pssm, seq_ext, "max", + bidirect = FALSE, strand = 1L, prior = 0.01 + ) + expect_equal(res$g[1], expected, tolerance = 1e-5, ignore_attr = TRUE) + # Sanity: the head-wins case must produce a strictly positive gradient + # (the head base is the consensus and column 0 is non-uniform). + expect_gt(res$g[1], 0) +}) + +test_that("pwm.grad max: head not argmax => engine returns 0 (matches oracle)", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 30L + offset <- 10L # interior anchor index whose match wins + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + # Construct a PSSM whose consensus matches the interior anchor at `offset`, + # NOT the head anchor at 0. Argmax should land at the interior, so head + # gradient is 0 (provided the head bases are not coincidentally also a + # consensus run, which we verify below). + interior_bases <- strsplit(substr(seq_ext, offset + 1, offset + L), "")[[1]] + pssm <- .consensus_pssm(interior_bases) + + # Guard the test setup: interior pattern must differ from head bases. + head_bases <- strsplit(substr(seq_ext, 1, L), "")[[1]] + skip_if( + identical(interior_bases, head_bases), + "test fixture: interior PSSM consensus collides with head bases; pick a different interval" + ) + + gvtrack.create("g", NULL, "pwm.grad", + pssm = pssm, + aggregate = "max", bidirect = FALSE, strand = 1, + extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g", sshift = 0, eshift = iter_W - 1) + + res <- gextract("g", iterator = 1, intervals = pivot) + + expected <- manual_pwm_grad(pssm, seq_ext, "max", + bidirect = FALSE, strand = 1L, prior = 0.01 + ) + expect_equal(res$g[1], expected, tolerance = 1e-5, ignore_attr = TRUE) + # And it should be exactly 0 in this construction (argmax interior). + expect_equal(expected, 0, tolerance = 1e-10) +}) + +test_that("pwm.grad max: bidirect=TRUE errors out (Task 3 scope)", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + pssm <- create_test_pssm() + gvtrack.create("g_bidir", NULL, "pwm.grad", + pssm = pssm, + aggregate = "max", bidirect = TRUE, + extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_bidir", sshift = 0, eshift = 19) + + expect_error( + gextract("g_bidir", iterator = 1, intervals = gintervals(1, 200, 201)), + "bidirect=FALSE" + ) +}) + +test_that("pwm.grad max: strand=-1 errors out (Task 3 scope)", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + pssm <- create_test_pssm() + gvtrack.create("g_minus", NULL, "pwm.grad", + pssm = pssm, + aggregate = "max", bidirect = FALSE, strand = -1, + extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_minus", sshift = 0, eshift = 19) + + expect_error( + gextract("g_minus", iterator = 1, intervals = gintervals(1, 200, 201)), + "strand=-1" + ) +}) From 2d2402eb62692caacf952421dc7c4cf16f5ef0cd Mon Sep 17 00:00:00 2001 From: aviezerl Date: Thu, 7 May 2026 23:00:30 +0300 Subject: [PATCH 05/15] feat(pwm.grad): GRAD_LSE single-strand fwd implementation --- src/PWMScorer.cpp | 85 ++++++++++++++++- src/PWMScorer.h | 2 + tests/testthat/test-vtrack-pwm-grad.R | 132 ++++++++++++++++++++++++++ 3 files changed, 215 insertions(+), 4 deletions(-) diff --git a/src/PWMScorer.cpp b/src/PWMScorer.cpp index b38efef58..025b34c9c 100644 --- a/src/PWMScorer.cpp +++ b/src/PWMScorer.cpp @@ -442,6 +442,80 @@ float PWMScorer::score_grad_max(const std::string& target, return m_pssm[0].get_log_prob_from_code(b_p) - m_worst_col0_fwd; } +// Linearized per-bp gradient (DeepLIFT-style) under LSE aggregation, fwd-only. +// +// Same pivot and clamping conventions as score_grad_max: pivot is target[0], +// argmax/LSE is taken over anchors with starts in [i_min, i_max], so the +// resulting f_LSE matches what `pwm` (TOTAL_LIKELIHOOD) would report on the +// same window with the same strand restriction. +// +// At the head anchor (i = 0), with the prior-adjusted log-probs M: +// w_p = exp(score_p - f_LSE) (softmax weight of head anchor) +// g = w_p * (M[0, b_p] - min_b M[0, b]) +// +// The "integration over W" comes through f_LSE's denominator, which spans all +// scanned anchors. When the head anchor is not in the scan ([i_min > 0]) the +// gradient is 0 by construction. +float PWMScorer::score_grad_lse(const std::string& target, + size_t i_min, size_t i_max, + size_t motif_length) +{ + // Precondition: bidirect/strand restrictions enforced upstream by the + // dispatch in score_interval (those verrors live outside the try/catch). + + if (i_min > i_max) { + return std::numeric_limits::quiet_NaN(); + } + if (target.size() < (size_t)motif_length) { + return std::numeric_limits::quiet_NaN(); + } + if (i_max + (size_t)motif_length > target.size()) { + return std::numeric_limits::quiet_NaN(); + } + + // Pivot is target[0]; head anchor only exists when i_min == 0. + if (i_min > 0) { + return 0.0f; + } + + // Head base. Non-ACGT -> NA (gradient is undefined when the pivot's base + // is unknown; matches score_grad_max). + int b_p = DnaLookupTables::BASE_ENCODE[(unsigned char)target[0]]; + if (b_p < 0) { + return std::numeric_limits::quiet_NaN(); + } + + // Scan and accumulate LSE; remember head score at i = 0. + float lse = -std::numeric_limits::infinity(); + float head = -std::numeric_limits::infinity(); + bool any_finite = false; + for (size_t i = i_min; i <= i_max; ++i) { + float v = 0.0f; + std::string::const_iterator it = target.begin() + i; + m_pssm.calc_like(it, v); + if (i == 0) { + head = v; + } + if (std::isfinite(v)) { + log_sum_log(lse, v); + any_finite = true; + } + } + + if (!any_finite) { + return std::numeric_limits::quiet_NaN(); + } + + // If the head anchor is itself -Inf, w_p = 0 and the gradient is 0. + if (!std::isfinite(head)) { + return 0.0f; + } + + const float w_p = std::exp(head - lse); + const float diff = m_pssm[0].get_log_prob_from_code(b_p) - m_worst_col0_fwd; + return w_p * diff; +} + // Score with spatial weighting float PWMScorer::score_with_spatial(const std::string& target, int64_t motif_length) { @@ -855,14 +929,14 @@ float PWMScorer::score_interval(const GInterval& interval, const GenomeChromKey& // below so that verror surfaces as an R-level error rather than being // swallowed and converted to NaN by the TGLException handler. switch (m_mode) { - case GRAD_LSE: case GRAD_LSE_ISM: case GRAD_MAX_ISM: rdb::verror("PWMScorer: gradient mode not yet implemented"); case GRAD_MAX: - // Task 3 supports only single-strand forward; reject bidirect - // PSSMs and strand=-1 here so the caller sees an error rather than - // a silent NaN. + case GRAD_LSE: + // GRAD_MAX (Task 3) and GRAD_LSE (Task 4) currently support only + // single-strand forward; reject bidirect PSSMs and strand=-1 here so + // the caller sees an error rather than a silent NaN. if (m_pssm.is_bidirect() || m_strand != 1) { rdb::verror("pwm.grad currently supports only bidirect=FALSE, strand=1; " "bidirect=TRUE and strand=-1 will be added in a future task"); @@ -948,6 +1022,9 @@ float PWMScorer::score_interval(const GInterval& interval, const GenomeChromKey& if (!m_use_spat && m_mode == GRAD_MAX) { return score_grad_max(target, i_min, i_max, motif_len); } + if (!m_use_spat && m_mode == GRAD_LSE) { + return score_grad_lse(target, i_min, i_max, motif_len); + } // Try spatial sliding window optimization if (m_use_spat) { diff --git a/src/PWMScorer.h b/src/PWMScorer.h index 71852067a..a381b625b 100644 --- a/src/PWMScorer.h +++ b/src/PWMScorer.h @@ -174,6 +174,8 @@ class PWMScorer : public GenomeSeqScorer // Gradient-mode answer paths (single-strand fwd; bidirect / rc deferred) float score_grad_max(const std::string& target, size_t i_min, size_t i_max, size_t motif_length); + float score_grad_lse(const std::string& target, size_t i_min, size_t i_max, + size_t motif_length); // Core members DnaPSSM m_pssm; diff --git a/tests/testthat/test-vtrack-pwm-grad.R b/tests/testthat/test-vtrack-pwm-grad.R index 8be9b6318..7b94b62e9 100644 --- a/tests/testthat/test-vtrack-pwm-grad.R +++ b/tests/testthat/test-vtrack-pwm-grad.R @@ -275,3 +275,135 @@ test_that("pwm.grad max: strand=-1 errors out (Task 3 scope)", { "strand=-1" ) }) + +# Task 4: GRAD_LSE single-strand fwd + +test_that("pwm.grad lse: concentrated softmax (head=consensus) matches oracle", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 20L + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + head_bases <- strsplit(substr(seq_ext, 1, L), "")[[1]] + pssm <- .consensus_pssm(head_bases) + + gvtrack.create("g_lse", NULL, "pwm.grad", + pssm = pssm, + aggregate = "lse", bidirect = FALSE, strand = 1, + extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_lse", sshift = 0, eshift = iter_W - 1) + res <- gextract("g_lse", iterator = 1, intervals = pivot) + + expected <- manual_pwm_grad(pssm, seq_ext, "lse", + bidirect = FALSE, strand = 1L, prior = 0.01 + ) + expect_equal(res$g_lse[1], expected, tolerance = 1e-5, ignore_attr = TRUE) + expect_gt(res$g_lse[1], 0) +}) + +test_that("pwm.grad lse: head not argmax => w_p small, g_lse < g_max", { + # Same construction as the head-not-argmax test for max, but compare lse + # vs max instead of just checking equality. + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 30L + offset <- 10L + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + interior_bases <- strsplit(substr(seq_ext, offset + 1, offset + L), "")[[1]] + pssm <- .consensus_pssm(interior_bases) + + head_bases <- strsplit(substr(seq_ext, 1, L), "")[[1]] + skip_if( + identical(interior_bases, head_bases), + "test fixture: interior PSSM consensus collides with head bases" + ) + + gvtrack.create("g_lse", NULL, "pwm.grad", + pssm = pssm, + aggregate = "lse", bidirect = FALSE, strand = 1, + extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_lse", sshift = 0, eshift = iter_W - 1) + res <- gextract("g_lse", iterator = 1, intervals = pivot) + + expected <- manual_pwm_grad(pssm, seq_ext, "lse", + bidirect = FALSE, strand = 1L, prior = 0.01 + ) + expect_equal(res$g_lse[1], expected, tolerance = 1e-5, ignore_attr = TRUE) + # When the interior anchor is the LSE-dominant one, w_p (head) is small, + # so g_lse stays close to 0 even though M[0, b_p] - worst is large. + expect_lt(res$g_lse[1], 0.5) +}) + +test_that("pwm.grad lse: w_p computed correctly via softmax denominator", { + # Hand-verifiable case: build a 2-anchor situation where head and one + # interior anchor have known scores, then compare the engine's w_p * diff + # against a direct computation. + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 5L # exactly 5 anchors + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + head_bases <- strsplit(substr(seq_ext, 1, L), "")[[1]] + pssm <- .consensus_pssm(head_bases) + + gvtrack.create("g_lse", NULL, "pwm.grad", + pssm = pssm, + aggregate = "lse", bidirect = FALSE, strand = 1, + extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_lse", sshift = 0, eshift = iter_W - 1) + res <- gextract("g_lse", iterator = 1, intervals = pivot) + + expected <- manual_pwm_grad(pssm, seq_ext, "lse", + bidirect = FALSE, strand = 1L, prior = 0.01 + ) + expect_equal(res$g_lse[1], expected, tolerance = 1e-5, ignore_attr = TRUE) + + # Cross-check: w_p must be in (0, 1] and g_lse <= g_max. + g_max_oracle <- manual_pwm_grad(pssm, seq_ext, "max", + bidirect = FALSE, strand = 1L, prior = 0.01 + ) + expect_lte(res$g_lse[1], g_max_oracle + 1e-8) + expect_gt(res$g_lse[1], 0) +}) + +test_that("pwm.grad lse: bidirect=TRUE / strand=-1 error out (Task 4 scope)", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + pssm <- create_test_pssm() + + gvtrack.create("g_bidir_lse", NULL, "pwm.grad", + pssm = pssm, + aggregate = "lse", bidirect = TRUE, + extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_bidir_lse", sshift = 0, eshift = 19) + expect_error( + gextract("g_bidir_lse", iterator = 1, intervals = gintervals(1, 200, 201)), + "bidirect=FALSE" + ) + + gvtrack.create("g_minus_lse", NULL, "pwm.grad", + pssm = pssm, + aggregate = "lse", bidirect = FALSE, strand = -1, + extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_minus_lse", sshift = 0, eshift = 19) + expect_error( + gextract("g_minus_lse", iterator = 1, intervals = gintervals(1, 200, 201)), + "strand=-1" + ) +}) From 0ac91556310a0af6fbf4c96ac411d6e29dc9c7c2 Mon Sep 17 00:00:00 2001 From: aviezerl Date: Thu, 7 May 2026 23:18:16 +0300 Subject: [PATCH 06/15] feat(pwm.grad): bidirect and rc-only support for linearized gradient 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. --- src/PWMScorer.cpp | 285 +++++++++++++++----------- src/PWMScorer.h | 12 +- tests/testthat/test-vtrack-pwm-grad.R | 189 +++++++++++++---- 3 files changed, 319 insertions(+), 167 deletions(-) diff --git a/src/PWMScorer.cpp b/src/PWMScorer.cpp index 025b34c9c..abe72d3be 100644 --- a/src/PWMScorer.cpp +++ b/src/PWMScorer.cpp @@ -366,30 +366,43 @@ float PWMScorer::score_without_spatial(const std::string& target, int64_t motif_ return compute_position_result(pos_idx, target.length(), motif_length, best_dir); } -// Linearized per-bp gradient (DeepLIFT-style) under MAX aggregation, fwd-only. +// Linearized per-bp gradient (DeepLIFT-style) at the iterator interval start. // -// Pivot is fixed at the iterator interval start (target[0] == genomic -// interval.start). The argmax is taken over anchors with starts in -// [i_min, i_max] inside the fetched/expanded target, which mirrors the clamping -// applied to pwm.max so the two functions agree on which anchor "wins". When -// the head anchor (i = 0) is the argmax, the gradient is -// M[0, b_p] - min_b M[0, b] -// where M[0, .] are the (prior-adjusted) log-probs at column 0 of the PSSM and -// b_p is the head base. Otherwise the gradient is 0. +// Pivot is fixed at target[0] (== genomic interval.start). Anchors are scanned +// over [i_min, i_max]; this matches the clamping applied to pwm/pwm.max so the +// gradient's denominator/argmax agrees with the un-gradient aggregate. // -// Restrictions for Task 3: single-strand fwd only. bidirect=TRUE (m_pssm -// bidirect) and strand=-1 are deferred to a later task (full bidirect / rc -// formula needs both fwd and rc head scores at the head anchor). -float PWMScorer::score_grad_max(const std::string& target, - size_t i_min, size_t i_max, - size_t motif_length) +// Strand handling: +// - Bidirect (m_pssm.is_bidirect() == true): per-anchor combined score is +// comb_a = log(exp(fwd_a) + exp(rc_a)). At the head anchor we split the +// diff by strand softmax: +// sm_fwd = exp(fwd_p - comb_p) +// sm_rc = 1 - sm_fwd +// diff_fwd = M[0, b_p] - min_b M[0, b] +// diff_rc = M[L-1, b_p_comp] - min_b M[L-1, b] (b_p_comp = 3 - b_p) +// combined_diff = sm_fwd * diff_fwd + sm_rc * diff_rc +// - Single-strand fwd (!is_bidirect, m_strand == 1): scan fwd only; +// diff = diff_fwd; no strand softmax. +// - Single-strand rc (!is_bidirect, m_strand == -1): scan rc only; +// diff = diff_rc; no strand softmax. +// +// Aggregation: +// - MAX (lse_aggregate == false): returns combined_diff (or single-strand +// diff) iff the argmax anchor index is 0 (head); else 0. +// - LSE (lse_aggregate == true): returns w_p * (combined_diff or +// single-strand diff), with w_p = exp(head_score - f_LSE). +// +// The "integration over W" enters through the LSE denominator, which spans all +// scanned anchors. +float PWMScorer::score_grad_linearized(const std::string& target, + size_t i_min, size_t i_max, + size_t motif_length, + bool lse_aggregate) { - // Precondition: bidirect/strand restrictions have already been enforced - // by the dispatch in score_interval (those verrors live outside the - // try/catch block so they reach R; calling verror here would be silently - // converted to NaN by the surrounding handler). + // Precondition: only the score_thresh / spatial restrictions need not be + // re-checked here. Caller is score_interval, which handles iterator + // clamping and parameter validation upstream. - // Iterator-interval was empty after clamping, or interval shorter than motif. if (i_min > i_max) { return std::numeric_limits::quiet_NaN(); } @@ -400,120 +413,155 @@ float PWMScorer::score_grad_max(const std::string& target, return std::numeric_limits::quiet_NaN(); } - // The pivot ("p" in the spec) is the iterator interval start, i.e. target[0]. - // It is the *head anchor* only when i_min == 0; otherwise the head anchor - // sits outside the iterator-interval scan and the gradient at p is 0 by - // definition (column 0 of every scanned anchor uses some other base). - if (i_min > 0) { + const bool bidirect = m_pssm.is_bidirect(); + const bool use_fwd = bidirect || m_strand == 1; + const bool use_rc = bidirect || m_strand == -1; + const size_t L = (size_t)motif_length; + + // The "head anchor" in fwd-genome (= the anchor at genomic interval.start) + // sits at a different target index depending on strand_mode: + // - strand_mode != -1 (target is fwd): head_idx = 0 + // - strand_mode == -1 (target is rc'd): head_idx = tlen - L + // If the head anchor is outside the iterator-clamped scan, the gradient at + // the pivot is 0 by construction (column 0 of every scanned anchor uses + // some other genomic base). + const size_t head_idx = (m_strand == -1) + ? (target.size() - (size_t)motif_length) + : 0; + if (head_idx < i_min || head_idx > i_max) { return 0.0f; } - // Scan fwd-only over [i_min, i_max] and find the argmax. - float best = -std::numeric_limits::infinity(); - size_t argmax = i_min; - bool any_finite = false; - for (size_t i = i_min; i <= i_max; ++i) { - float v = 0.0f; - std::string::const_iterator it = target.begin() + i; - m_pssm.calc_like(it, v); - if (std::isfinite(v) && v > best) { - best = v; - argmax = i; - any_finite = true; - } - } - - if (!any_finite) { - // All anchors hit non-ACGT bases (-Inf): cannot define a gradient. - return std::numeric_limits::quiet_NaN(); + // Pivot's fwd-genome base, in fwd PSSM coords (b_p_fwd) plus its complement + // for rc-strand contributions. + int b_p_fwd; + if (m_strand == -1) { + const unsigned char tail = (unsigned char)target[target.size() - 1]; + b_p_fwd = DnaLookupTables::COMPLEMENT_ENCODE[tail]; + } else { + b_p_fwd = DnaLookupTables::BASE_ENCODE[(unsigned char)target[0]]; } - - // Head base (target[0]). N or other non-ACGT -> NA. - int b_p = DnaLookupTables::BASE_ENCODE[(unsigned char)target[0]]; - if (b_p < 0) { + if (b_p_fwd < 0) { return std::numeric_limits::quiet_NaN(); } + const int b_p_fwd_comp = 3 - b_p_fwd; + + // Per-strand diffs at the head anchor (in fwd-genome semantics). + const float diff_fwd = m_pssm[0].get_log_prob_from_code(b_p_fwd) - m_worst_col0_fwd; + const float diff_rc = m_pssm[L - 1].get_log_prob_from_code(b_p_fwd_comp) - m_worst_col_last_fwd; + + // Scan and accumulate. We track: + // - LSE over per-anchor scores (used iff lse_aggregate) + // - argmax index of per-anchor scores (used iff !lse_aggregate) + // - head fwd / rc / comb scores (for the strand softmax at the head anchor) + // - head per-anchor score (for tie-break: head wins ties) + float lse = -std::numeric_limits::infinity(); + float best = -std::numeric_limits::infinity(); + size_t argmax = i_min; + float head_fwd = -std::numeric_limits::infinity(); + float head_rc = -std::numeric_limits::infinity(); + float head_score = -std::numeric_limits::infinity(); + bool any_finite = false; - if (argmax != 0) { - return 0.0f; - } + for (size_t i = i_min; i <= i_max; ++i) { + // score_forward_original / score_reverse_original return scores in + // fwd-genome semantics regardless of whether `target` is rc'd. + float fwd = -std::numeric_limits::infinity(); + float rc = -std::numeric_limits::infinity(); + if (use_fwd) { + fwd = score_forward_original(m_pssm, target, i, m_strand); + } + if (use_rc) { + rc = score_reverse_original(m_pssm, target, i, m_strand); + } - // Argmax is at the head anchor: gradient = M[0, b_p] - worst_col0. - return m_pssm[0].get_log_prob_from_code(b_p) - m_worst_col0_fwd; -} + // Per-anchor score: comb (bidirect) or single-strand. + float s; + if (bidirect) { + if (!std::isfinite(fwd) && !std::isfinite(rc)) { + s = -std::numeric_limits::infinity(); + } else if (!std::isfinite(fwd)) { + s = rc; + } else if (!std::isfinite(rc)) { + s = fwd; + } else { + s = fwd; + log_sum_log(s, rc); + } + } else { + s = use_fwd ? fwd : rc; + } -// Linearized per-bp gradient (DeepLIFT-style) under LSE aggregation, fwd-only. -// -// Same pivot and clamping conventions as score_grad_max: pivot is target[0], -// argmax/LSE is taken over anchors with starts in [i_min, i_max], so the -// resulting f_LSE matches what `pwm` (TOTAL_LIKELIHOOD) would report on the -// same window with the same strand restriction. -// -// At the head anchor (i = 0), with the prior-adjusted log-probs M: -// w_p = exp(score_p - f_LSE) (softmax weight of head anchor) -// g = w_p * (M[0, b_p] - min_b M[0, b]) -// -// The "integration over W" comes through f_LSE's denominator, which spans all -// scanned anchors. When the head anchor is not in the scan ([i_min > 0]) the -// gradient is 0 by construction. -float PWMScorer::score_grad_lse(const std::string& target, - size_t i_min, size_t i_max, - size_t motif_length) -{ - // Precondition: bidirect/strand restrictions enforced upstream by the - // dispatch in score_interval (those verrors live outside the try/catch). + if (i == head_idx) { + head_fwd = fwd; + head_rc = rc; + head_score = s; + } - if (i_min > i_max) { - return std::numeric_limits::quiet_NaN(); - } - if (target.size() < (size_t)motif_length) { - return std::numeric_limits::quiet_NaN(); - } - if (i_max + (size_t)motif_length > target.size()) { - return std::numeric_limits::quiet_NaN(); + if (std::isfinite(s)) { + any_finite = true; + if (lse_aggregate) { + log_sum_log(lse, s); + } + if (s > best) { + best = s; + argmax = i; + } + } } - // Pivot is target[0]; head anchor only exists when i_min == 0. - if (i_min > 0) { - return 0.0f; + // Tie-break: if the head anchor matches `best`, prefer it. This matches + // the oracle's R `which.max` semantics, which returns the first index of + // the max in fwd-anchor order. For strand=-1 the rc'd-target scan visits + // fwd-anchors in reverse order, so without this tie-break ties land on + // the wrong anchor. + if (std::isfinite(head_score) && head_score == best) { + argmax = head_idx; } - // Head base. Non-ACGT -> NA (gradient is undefined when the pivot's base - // is unknown; matches score_grad_max). - int b_p = DnaLookupTables::BASE_ENCODE[(unsigned char)target[0]]; - if (b_p < 0) { + if (!any_finite) { return std::numeric_limits::quiet_NaN(); } - // Scan and accumulate LSE; remember head score at i = 0. - float lse = -std::numeric_limits::infinity(); - float head = -std::numeric_limits::infinity(); - bool any_finite = false; - for (size_t i = i_min; i <= i_max; ++i) { - float v = 0.0f; - std::string::const_iterator it = target.begin() + i; - m_pssm.calc_like(it, v); - if (i == 0) { - head = v; + // Combined diff at the head anchor. + float combined_diff; + if (bidirect) { + // If the head's comb score is -Inf (both strands hit -Inf at i=0), + // gradient is 0 (head doesn't participate in the aggregate). + const bool fwd_finite = std::isfinite(head_fwd); + const bool rc_finite = std::isfinite(head_rc); + if (!fwd_finite && !rc_finite) { + return 0.0f; } - if (std::isfinite(v)) { - log_sum_log(lse, v); - any_finite = true; + if (!fwd_finite) { + combined_diff = diff_rc; + } else if (!rc_finite) { + combined_diff = diff_fwd; + } else { + float comb_p = head_fwd; + log_sum_log(comb_p, head_rc); + const float sm_fwd = std::exp(head_fwd - comb_p); + const float sm_rc = 1.0f - sm_fwd; + combined_diff = sm_fwd * diff_fwd + sm_rc * diff_rc; } + } else { + const float head_single = use_fwd ? head_fwd : head_rc; + if (!std::isfinite(head_single)) { + return 0.0f; + } + combined_diff = use_fwd ? diff_fwd : diff_rc; } - if (!any_finite) { - return std::numeric_limits::quiet_NaN(); + if (!lse_aggregate) { + return (argmax == head_idx) ? combined_diff : 0.0f; } - // If the head anchor is itself -Inf, w_p = 0 and the gradient is 0. - if (!std::isfinite(head)) { + // LSE: weight by softmax of the head anchor's per-anchor score. + if (!std::isfinite(head_score)) { return 0.0f; } - - const float w_p = std::exp(head - lse); - const float diff = m_pssm[0].get_log_prob_from_code(b_p) - m_worst_col0_fwd; - return w_p * diff; + const float w_p = std::exp(head_score - lse); + return w_p * combined_diff; } // Score with spatial weighting @@ -934,12 +982,11 @@ float PWMScorer::score_interval(const GInterval& interval, const GenomeChromKey& rdb::verror("PWMScorer: gradient mode not yet implemented"); case GRAD_MAX: case GRAD_LSE: - // GRAD_MAX (Task 3) and GRAD_LSE (Task 4) currently support only - // single-strand forward; reject bidirect PSSMs and strand=-1 here so - // the caller sees an error rather than a silent NaN. - if (m_pssm.is_bidirect() || m_strand != 1) { - rdb::verror("pwm.grad currently supports only bidirect=FALSE, strand=1; " - "bidirect=TRUE and strand=-1 will be added in a future task"); + // Linearized gradient supports bidirect, single-strand fwd, and + // single-strand rc. m_strand must be in {-1, 0, 1} (the parser + // already enforces this; defensive check). + if (m_strand != -1 && m_strand != 0 && m_strand != 1) { + rdb::verror("pwm.grad: strand must be -1, 0, or 1"); } break; default: @@ -1019,11 +1066,9 @@ float PWMScorer::score_interval(const GInterval& interval, const GenomeChromKey& // Gradient modes (single-strand fwd) take a dedicated answer path. // They share the iterator clamping (i_min..i_max) above so that the // argmax convention matches pwm.max exactly. - if (!m_use_spat && m_mode == GRAD_MAX) { - return score_grad_max(target, i_min, i_max, motif_len); - } - if (!m_use_spat && m_mode == GRAD_LSE) { - return score_grad_lse(target, i_min, i_max, motif_len); + if (!m_use_spat && (m_mode == GRAD_MAX || m_mode == GRAD_LSE)) { + const bool lse = (m_mode == GRAD_LSE); + return score_grad_linearized(target, i_min, i_max, motif_len, lse); } // Try spatial sliding window optimization diff --git a/src/PWMScorer.h b/src/PWMScorer.h index a381b625b..2b46f48a5 100644 --- a/src/PWMScorer.h +++ b/src/PWMScorer.h @@ -171,11 +171,13 @@ class PWMScorer : public GenomeSeqScorer // Utilities inline float get_spatial_log_factor(size_t pos_index) const; - // Gradient-mode answer paths (single-strand fwd; bidirect / rc deferred) - float score_grad_max(const std::string& target, size_t i_min, size_t i_max, - size_t motif_length); - float score_grad_lse(const std::string& target, size_t i_min, size_t i_max, - size_t motif_length); + // Linearized per-bp gradient (DeepLIFT-style) at the iterator interval start. + // Handles fwd-only, rc-only, and bidirect strand modes with either MAX + // (`lse_aggregate = false`) or LSE (`lse_aggregate = true`) aggregation. + float score_grad_linearized(const std::string& target, + size_t i_min, size_t i_max, + size_t motif_length, + bool lse_aggregate); // Core members DnaPSSM m_pssm; diff --git a/tests/testthat/test-vtrack-pwm-grad.R b/tests/testthat/test-vtrack-pwm-grad.R index 7b94b62e9..1cd6578c7 100644 --- a/tests/testthat/test-vtrack-pwm-grad.R +++ b/tests/testthat/test-vtrack-pwm-grad.R @@ -240,40 +240,174 @@ test_that("pwm.grad max: head not argmax => engine returns 0 (matches oracle)", expect_equal(expected, 0, tolerance = 1e-10) }) -test_that("pwm.grad max: bidirect=TRUE errors out (Task 3 scope)", { +test_that("pwm.grad max: bidirect=TRUE matches oracle on asymmetric PSSM", { remove_all_vtracks() withr::defer(remove_all_vtracks()) - pssm <- create_test_pssm() + L <- 3L + iter_W <- 20L + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + # Asymmetric PSSM (rc != self) so bidirect's strand-softmax is non-trivial. + pssm <- matrix( + c( + 0.7, 0.1, 0.1, 0.1, + 0.1, 0.6, 0.2, 0.1, + 0.1, 0.1, 0.5, 0.3 + ), + L, 4, + byrow = TRUE, + dimnames = list(NULL, c("A", "C", "G", "T")) + ) + gvtrack.create("g_bidir", NULL, "pwm.grad", pssm = pssm, aggregate = "max", bidirect = TRUE, extend = TRUE, prior = 0.01 ) - gvtrack.iterator("g_bidir", sshift = 0, eshift = 19) + gvtrack.iterator("g_bidir", sshift = 0, eshift = iter_W - 1) + res <- gextract("g_bidir", iterator = 1, intervals = pivot) - expect_error( - gextract("g_bidir", iterator = 1, intervals = gintervals(1, 200, 201)), - "bidirect=FALSE" + expected <- manual_pwm_grad(pssm, seq_ext, "max", + bidirect = TRUE, prior = 0.01 ) + expect_equal(res$g_bidir[1], expected, tolerance = 1e-5, ignore_attr = TRUE) }) -test_that("pwm.grad max: strand=-1 errors out (Task 3 scope)", { +test_that("pwm.grad lse: bidirect=TRUE matches oracle on asymmetric PSSM", { remove_all_vtracks() withr::defer(remove_all_vtracks()) - pssm <- create_test_pssm() - gvtrack.create("g_minus", NULL, "pwm.grad", + L <- 3L + iter_W <- 20L + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + pssm <- matrix( + c( + 0.7, 0.1, 0.1, 0.1, + 0.1, 0.6, 0.2, 0.1, + 0.1, 0.1, 0.5, 0.3 + ), + L, 4, + byrow = TRUE, + dimnames = list(NULL, c("A", "C", "G", "T")) + ) + + gvtrack.create("g_bidir_lse", NULL, "pwm.grad", + pssm = pssm, + aggregate = "lse", bidirect = TRUE, + extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_bidir_lse", sshift = 0, eshift = iter_W - 1) + res <- gextract("g_bidir_lse", iterator = 1, intervals = pivot) + + expected <- manual_pwm_grad(pssm, seq_ext, "lse", + bidirect = TRUE, prior = 0.01 + ) + expect_equal(res$g_bidir_lse[1], expected, tolerance = 1e-5, ignore_attr = TRUE) +}) + +test_that("pwm.grad: bidirect on palindromic PSSM is consistent", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + # Palindromic: rc(pssm) == pssm. The bidirect comb score is fwd + log(2) + # (offset, irrelevant for argmax). Strand softmax at any anchor is 0.5/0.5. + pssm <- matrix( + c( + 0.7, 0.1, 0.1, 0.1, + 0.1, 0.4, 0.4, 0.1, + 0.1, 0.1, 0.1, 0.7 + ), + 3, 4, + byrow = TRUE, + dimnames = list(NULL, c("A", "C", "G", "T")) + ) + rc_test <- pssm[3:1, c(4, 3, 2, 1)] + dimnames(rc_test) <- dimnames(pssm) + expect_equal(rc_test, pssm) # confirm palindromic + + iter_W <- 20L + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + nrow(pssm) - 1))) + + gvtrack.create("g_bidir", NULL, "pwm.grad", pssm = pssm, + aggregate = "lse", bidirect = TRUE, + extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_bidir", sshift = 0, eshift = iter_W - 1) + res <- gextract("g_bidir", iterator = 1, intervals = pivot) + + expected <- manual_pwm_grad(pssm, seq_ext, "lse", + bidirect = TRUE, prior = 0.01 + ) + expect_equal(res$g_bidir[1], expected, tolerance = 1e-5, ignore_attr = TRUE) +}) + +test_that("pwm.grad max: bidirect=FALSE strand=-1 matches oracle", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 20L + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + # Build PSSM whose rc-consensus matches the first L bases of seq_ext, so + # rc-PSSM scoring at the head anchor is strong. + head_bases <- strsplit(substr(seq_ext, 1, L), "")[[1]] + # rc-PSSM at column 0 sees complement(head[0]); we want a consensus PSSM + # whose rc orientation has strong matches at our head bases. Simplest: use + # the consensus PSSM for the actual head bases; the rc scoring will see + # complement-of-head at L-1, L-2, ... so it will be a weak match. To make + # rc scoring strong, use the rc of the consensus PSSM instead. + fwd_pssm <- .consensus_pssm(head_bases) + rc_target_pssm <- fwd_pssm[L:1, c(4, 3, 2, 1)] + dimnames(rc_target_pssm) <- dimnames(fwd_pssm) + + gvtrack.create("g_minus", NULL, "pwm.grad", + pssm = rc_target_pssm, aggregate = "max", bidirect = FALSE, strand = -1, extend = TRUE, prior = 0.01 ) - gvtrack.iterator("g_minus", sshift = 0, eshift = 19) + gvtrack.iterator("g_minus", sshift = 0, eshift = iter_W - 1) + res <- gextract("g_minus", iterator = 1, intervals = pivot) - expect_error( - gextract("g_minus", iterator = 1, intervals = gintervals(1, 200, 201)), - "strand=-1" + expected <- manual_pwm_grad(rc_target_pssm, seq_ext, "max", + bidirect = FALSE, strand = -1L, prior = 0.01 ) + expect_equal(res$g_minus[1], expected, tolerance = 1e-5, ignore_attr = TRUE) +}) + +test_that("pwm.grad lse: bidirect=FALSE strand=-1 matches oracle", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 20L + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + head_bases <- strsplit(substr(seq_ext, 1, L), "")[[1]] + fwd_pssm <- .consensus_pssm(head_bases) + rc_target_pssm <- fwd_pssm[L:1, c(4, 3, 2, 1)] + dimnames(rc_target_pssm) <- dimnames(fwd_pssm) + + gvtrack.create("g_minus_lse", NULL, "pwm.grad", + pssm = rc_target_pssm, + aggregate = "lse", bidirect = FALSE, strand = -1, + extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_minus_lse", sshift = 0, eshift = iter_W - 1) + res <- gextract("g_minus_lse", iterator = 1, intervals = pivot) + + expected <- manual_pwm_grad(rc_target_pssm, seq_ext, "lse", + bidirect = FALSE, strand = -1L, prior = 0.01 + ) + expect_equal(res$g_minus_lse[1], expected, tolerance = 1e-5, ignore_attr = TRUE) }) # Task 4: GRAD_LSE single-strand fwd @@ -378,32 +512,3 @@ test_that("pwm.grad lse: w_p computed correctly via softmax denominator", { expect_lte(res$g_lse[1], g_max_oracle + 1e-8) expect_gt(res$g_lse[1], 0) }) - -test_that("pwm.grad lse: bidirect=TRUE / strand=-1 error out (Task 4 scope)", { - remove_all_vtracks() - withr::defer(remove_all_vtracks()) - - pssm <- create_test_pssm() - - gvtrack.create("g_bidir_lse", NULL, "pwm.grad", - pssm = pssm, - aggregate = "lse", bidirect = TRUE, - extend = TRUE, prior = 0.01 - ) - gvtrack.iterator("g_bidir_lse", sshift = 0, eshift = 19) - expect_error( - gextract("g_bidir_lse", iterator = 1, intervals = gintervals(1, 200, 201)), - "bidirect=FALSE" - ) - - gvtrack.create("g_minus_lse", NULL, "pwm.grad", - pssm = pssm, - aggregate = "lse", bidirect = FALSE, strand = -1, - extend = TRUE, prior = 0.01 - ) - gvtrack.iterator("g_minus_lse", sshift = 0, eshift = 19) - expect_error( - gextract("g_minus_lse", iterator = 1, intervals = gintervals(1, 200, 201)), - "strand=-1" - ) -}) From f28031dcf1907365dab8ed61347cd15208725f3f Mon Sep 17 00:00:00 2001 From: aviezerl Date: Thu, 7 May 2026 23:29:58 +0300 Subject: [PATCH 07/15] feat(pwm.grad): in-silico mutagenesis (Tasks 6/7/8 unified) 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. --- src/PWMScorer.cpp | 211 +++++++++++++++++------ src/PWMScorer.h | 14 +- tests/testthat/test-vtrack-pwm-grad.R | 231 ++++++++++++++++++++++++++ 3 files changed, 400 insertions(+), 56 deletions(-) diff --git a/src/PWMScorer.cpp b/src/PWMScorer.cpp index abe72d3be..1fa29e6a6 100644 --- a/src/PWMScorer.cpp +++ b/src/PWMScorer.cpp @@ -366,38 +366,44 @@ float PWMScorer::score_without_spatial(const std::string& target, int64_t motif_ return compute_position_result(pos_idx, target.length(), motif_length, best_dir); } -// Linearized per-bp gradient (DeepLIFT-style) at the iterator interval start. +// Per-bp gradient at the iterator interval start. // -// Pivot is fixed at target[0] (== genomic interval.start). Anchors are scanned -// over [i_min, i_max]; this matches the clamping applied to pwm/pwm.max so the -// gradient's denominator/argmax agrees with the un-gradient aggregate. +// Two definitions are available: +// - Linearized (ism=false): "actual minus worst" weighted by softmax (LSE) +// or argmax indicator (MAX). Closed-form, O(1) on top of the scan. +// - In-silico mutagenesis (ism=true): f(actual) - min_b' f(seq with seq[p] +// flipped to b'). With the pivot at the head anchor, only the head's +// per-anchor score changes under a flip, so 3 alternative aggregates can +// be computed in O(1) each. // -// Strand handling: -// - Bidirect (m_pssm.is_bidirect() == true): per-anchor combined score is -// comb_a = log(exp(fwd_a) + exp(rc_a)). At the head anchor we split the -// diff by strand softmax: -// sm_fwd = exp(fwd_p - comb_p) -// sm_rc = 1 - sm_fwd -// diff_fwd = M[0, b_p] - min_b M[0, b] -// diff_rc = M[L-1, b_p_comp] - min_b M[L-1, b] (b_p_comp = 3 - b_p) -// combined_diff = sm_fwd * diff_fwd + sm_rc * diff_rc -// - Single-strand fwd (!is_bidirect, m_strand == 1): scan fwd only; -// diff = diff_fwd; no strand softmax. -// - Single-strand rc (!is_bidirect, m_strand == -1): scan rc only; -// diff = diff_rc; no strand softmax. +// Pivot is the iterator interval start (target[0] in fwd-target, target[tlen-1] +// in rc'd target). Anchors are scanned over [i_min, i_max]; this matches the +// clamping applied to pwm/pwm.max so the denominator/argmax agree with the +// un-gradient aggregate. // -// Aggregation: -// - MAX (lse_aggregate == false): returns combined_diff (or single-strand -// diff) iff the argmax anchor index is 0 (head); else 0. -// - LSE (lse_aggregate == true): returns w_p * (combined_diff or -// single-strand diff), with w_p = exp(head_score - f_LSE). +// Strand handling unifies bidirect / fwd-only / rc-only via score_*_original, +// which canonicalize the rc'd target back to fwd-genome semantics. // -// The "integration over W" enters through the LSE denominator, which spans all -// scanned anchors. -float PWMScorer::score_grad_linearized(const std::string& target, - size_t i_min, size_t i_max, - size_t motif_length, - bool lse_aggregate) +// Linearized result: +// - MAX: return combined_diff if argmax == head_idx, else 0 +// - LSE: return w_p * combined_diff, w_p = exp(head_score - f_LSE) +// where combined_diff = sm_fwd * diff_fwd + sm_rc * diff_rc (bidirect) or +// the single-strand diff. diff_fwd = M[0,b_p] - min_b M[0,b], +// diff_rc = M[L-1,b_p_comp] - min_b M[L-1,b]. +// +// ISM result: flip the pivot's base to each of 3 alternatives. Only the head +// anchor's per-anchor score changes (column 0 fwd, last col rc, or both for +// bidirect). For each alternative b': +// - fwd_alt = head_fwd - M[0, b_p_fwd] + M[0, b'] +// - rc_alt = head_rc - M[L-1, b_p_fwd_comp] + M[L-1, comp(b')] +// - s_alt = log(exp(fwd_alt) + exp(rc_alt)) (bidirect) or single-strand +// - MAX: f_alt = max(s_alt, max_excluding_head) +// - LSE: f_alt = log(exp(f_LSE) - exp(head_score) + exp(s_alt)) [stable form] +// The ISM gradient is f_actual - min_b' f_alt. +float PWMScorer::score_grad(const std::string& target, + size_t i_min, size_t i_max, + size_t motif_length, + bool lse_aggregate, bool ism) { // Precondition: only the score_thresh / spatial restrictions need not be // re-checked here. Caller is score_interval, which handles iterator @@ -450,13 +456,13 @@ float PWMScorer::score_grad_linearized(const std::string& target, const float diff_fwd = m_pssm[0].get_log_prob_from_code(b_p_fwd) - m_worst_col0_fwd; const float diff_rc = m_pssm[L - 1].get_log_prob_from_code(b_p_fwd_comp) - m_worst_col_last_fwd; - // Scan and accumulate. We track: - // - LSE over per-anchor scores (used iff lse_aggregate) - // - argmax index of per-anchor scores (used iff !lse_aggregate) - // - head fwd / rc / comb scores (for the strand softmax at the head anchor) - // - head per-anchor score (for tie-break: head wins ties) + // Scan and accumulate. We track all-anchor stats plus excluding-head + // stats; the latter let ISM compute max/LSE after a single-base flip in + // O(1) without a second pass. float lse = -std::numeric_limits::infinity(); + float lse_no_head = -std::numeric_limits::infinity(); float best = -std::numeric_limits::infinity(); + float best_no_head = -std::numeric_limits::infinity(); size_t argmax = i_min; float head_fwd = -std::numeric_limits::infinity(); float head_rc = -std::numeric_limits::infinity(); @@ -500,13 +506,17 @@ float PWMScorer::score_grad_linearized(const std::string& target, if (std::isfinite(s)) { any_finite = true; - if (lse_aggregate) { - log_sum_log(lse, s); - } + log_sum_log(lse, s); if (s > best) { best = s; argmax = i; } + if (i != head_idx) { + log_sum_log(lse_no_head, s); + if (s > best_no_head) { + best_no_head = s; + } + } } } @@ -552,16 +562,115 @@ float PWMScorer::score_grad_linearized(const std::string& target, combined_diff = use_fwd ? diff_fwd : diff_rc; } - if (!lse_aggregate) { - return (argmax == head_idx) ? combined_diff : 0.0f; + if (!ism) { + // Linearized gradient. + if (!lse_aggregate) { + return (argmax == head_idx) ? combined_diff : 0.0f; + } + if (!std::isfinite(head_score)) { + return 0.0f; + } + const float w_p = std::exp(head_score - lse); + return w_p * combined_diff; } - // LSE: weight by softmax of the head anchor's per-anchor score. - if (!std::isfinite(head_score)) { - return 0.0f; + // ---- ISM --------------------------------------------------------------- + // + // Flip the pivot's base to each of 3 alternatives, recompute the head + // anchor's per-anchor score in O(1), recompute the aggregate in O(1) using + // the excluding-head stats, take the min, return f_actual - min_alt. + // + // If the head's per-anchor score is -Inf (head doesn't participate), the + // gradient is f_actual - min_b' f_alt where every f_alt equals the + // excluding-head aggregate (head still won't participate after flip if + // the alt also becomes -Inf, which won't happen since alt is always ACGT). + // Fall through to the general formula below. + const float f_actual = lse_aggregate ? lse : best; + if (!std::isfinite(f_actual)) { + return std::numeric_limits::quiet_NaN(); } - const float w_p = std::exp(head_score - lse); - return w_p * combined_diff; + + // Per-column log-prob deltas: subtract the actual base's contribution, + // add the alternative's. fwd uses M[0, .], rc uses M[L-1, complement(.)]. + const float fwd_minus_actual = use_fwd + ? m_pssm[0].get_log_prob_from_code(b_p_fwd) : 0.0f; + const float rc_minus_actual = use_rc + ? m_pssm[L - 1].get_log_prob_from_code(b_p_fwd_comp) : 0.0f; + + float min_f_alt = std::numeric_limits::infinity(); + for (int b_alt = 0; b_alt < 4; ++b_alt) { + if (b_alt == b_p_fwd) continue; + + // Flipped per-strand head scores (in fwd-genome semantics). + float fwd_alt = -std::numeric_limits::infinity(); + float rc_alt = -std::numeric_limits::infinity(); + if (use_fwd) { + const float add = m_pssm[0].get_log_prob_from_code(b_alt); + if (std::isfinite(head_fwd) && std::isfinite(add)) { + fwd_alt = head_fwd - fwd_minus_actual + add; + } else if (std::isfinite(add)) { + // head_fwd was -Inf because actual base hit a -Inf prob col. + // The flipped score is just the column-by-column sum with the + // new base at column 0. We can't recover that without a full + // re-eval; treat as -Inf (the alt also can't fix unrelated + // -Inf columns, which only happen for non-ACGT bases at non- + // pivot positions). + } + } + if (use_rc) { + const int b_alt_comp = 3 - b_alt; + const float add = m_pssm[L - 1].get_log_prob_from_code(b_alt_comp); + if (std::isfinite(head_rc) && std::isfinite(add)) { + rc_alt = head_rc - rc_minus_actual + add; + } + } + + // Combined per-anchor score after flip. + float s_alt; + if (bidirect) { + if (!std::isfinite(fwd_alt) && !std::isfinite(rc_alt)) { + s_alt = -std::numeric_limits::infinity(); + } else if (!std::isfinite(fwd_alt)) { + s_alt = rc_alt; + } else if (!std::isfinite(rc_alt)) { + s_alt = fwd_alt; + } else { + s_alt = fwd_alt; + log_sum_log(s_alt, rc_alt); + } + } else { + s_alt = use_fwd ? fwd_alt : rc_alt; + } + + // Aggregate after flip: only the head anchor's score changes. + float f_alt; + if (lse_aggregate) { + // f_alt = log(exp(lse_no_head) + exp(s_alt)) + if (!std::isfinite(s_alt)) { + f_alt = lse_no_head; + } else if (!std::isfinite(lse_no_head)) { + f_alt = s_alt; + } else { + f_alt = s_alt; + log_sum_log(f_alt, lse_no_head); + } + } else { + // f_alt = max(s_alt, best_no_head) + f_alt = std::max(s_alt, best_no_head); + } + + if (f_alt < min_f_alt) { + min_f_alt = f_alt; + } + } + + if (!std::isfinite(min_f_alt)) { + // All alternatives produced -Inf; gradient undefined. + return std::numeric_limits::quiet_NaN(); + } + const float g = f_actual - min_f_alt; + // Clamp tiny numerical underflow to zero. + return (g < 0.0f && g > -1e-5f) ? 0.0f : g; } // Score with spatial weighting @@ -977,12 +1086,11 @@ float PWMScorer::score_interval(const GInterval& interval, const GenomeChromKey& // below so that verror surfaces as an R-level error rather than being // swallowed and converted to NaN by the TGLException handler. switch (m_mode) { - case GRAD_LSE_ISM: - case GRAD_MAX_ISM: - rdb::verror("PWMScorer: gradient mode not yet implemented"); case GRAD_MAX: case GRAD_LSE: - // Linearized gradient supports bidirect, single-strand fwd, and + case GRAD_MAX_ISM: + case GRAD_LSE_ISM: + // All gradient modes support bidirect, single-strand fwd, and // single-strand rc. m_strand must be in {-1, 0, 1} (the parser // already enforces this; defensive check). if (m_strand != -1 && m_strand != 0 && m_strand != 1) { @@ -1066,9 +1174,12 @@ float PWMScorer::score_interval(const GInterval& interval, const GenomeChromKey& // Gradient modes (single-strand fwd) take a dedicated answer path. // They share the iterator clamping (i_min..i_max) above so that the // argmax convention matches pwm.max exactly. - if (!m_use_spat && (m_mode == GRAD_MAX || m_mode == GRAD_LSE)) { - const bool lse = (m_mode == GRAD_LSE); - return score_grad_linearized(target, i_min, i_max, motif_len, lse); + if (!m_use_spat && + (m_mode == GRAD_MAX || m_mode == GRAD_LSE || + m_mode == GRAD_MAX_ISM || m_mode == GRAD_LSE_ISM)) { + const bool lse_aggregate = (m_mode == GRAD_LSE || m_mode == GRAD_LSE_ISM); + const bool ism = (m_mode == GRAD_MAX_ISM || m_mode == GRAD_LSE_ISM); + return score_grad(target, i_min, i_max, motif_len, lse_aggregate, ism); } // Try spatial sliding window optimization diff --git a/src/PWMScorer.h b/src/PWMScorer.h index 2b46f48a5..79b063876 100644 --- a/src/PWMScorer.h +++ b/src/PWMScorer.h @@ -171,13 +171,15 @@ class PWMScorer : public GenomeSeqScorer // Utilities inline float get_spatial_log_factor(size_t pos_index) const; - // Linearized per-bp gradient (DeepLIFT-style) at the iterator interval start. + // Per-bp gradient at the iterator interval start. // Handles fwd-only, rc-only, and bidirect strand modes with either MAX - // (`lse_aggregate = false`) or LSE (`lse_aggregate = true`) aggregation. - float score_grad_linearized(const std::string& target, - size_t i_min, size_t i_max, - size_t motif_length, - bool lse_aggregate); + // (lse_aggregate=false) or LSE (lse_aggregate=true) aggregation, and + // either linearized (ism=false) or in-silico-mutagenesis (ism=true) + // gradient definitions. + float score_grad(const std::string& target, + size_t i_min, size_t i_max, + size_t motif_length, + bool lse_aggregate, bool ism); // Core members DnaPSSM m_pssm; diff --git a/tests/testthat/test-vtrack-pwm-grad.R b/tests/testthat/test-vtrack-pwm-grad.R index 1cd6578c7..954238cff 100644 --- a/tests/testthat/test-vtrack-pwm-grad.R +++ b/tests/testthat/test-vtrack-pwm-grad.R @@ -512,3 +512,234 @@ test_that("pwm.grad lse: w_p computed correctly via softmax denominator", { expect_lte(res$g_lse[1], g_max_oracle + 1e-8) expect_gt(res$g_lse[1], 0) }) + +# Task 6: GRAD_MAX_ISM tests + +test_that("pwm.grad.ism max: argmax-shift case differs from linearized", { + # The test-DB seq at chr1:200-222 is periodic ("CCCTAACCC..."), so even + # with a strong consensus PSSM matching the head, flipping the head base + # doesn't drop the max because anchor 7 also matches the consensus. + # Linearized would report ~3.9 (the col-0 diff); ISM correctly returns 0 + # because the argmax shifts under the flip. This is the canonical + # ISM-vs-linearized divergence case. + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 20L + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + head_bases <- strsplit(substr(seq_ext, 1, L), "")[[1]] + pssm <- .consensus_pssm(head_bases) + + gvtrack.create("g_ism", NULL, "pwm.grad.ism", + pssm = pssm, aggregate = "max", + bidirect = FALSE, strand = 1, extend = TRUE, prior = 0.01 + ) + gvtrack.create("g_lin", NULL, "pwm.grad", + pssm = pssm, aggregate = "max", + bidirect = FALSE, strand = 1, extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_ism", sshift = 0, eshift = iter_W - 1) + gvtrack.iterator("g_lin", sshift = 0, eshift = iter_W - 1) + res <- gextract(c("g_ism", "g_lin"), iterator = 1, intervals = pivot) + + expected_ism <- manual_pwm_grad_ism(pssm, seq_ext, "max", + bidirect = FALSE, strand = 1L, prior = 0.01 + ) + expect_equal(res$g_ism[1], expected_ism, tolerance = 1e-5, ignore_attr = TRUE) + # ISM = 0 (argmax shift); linearized > 0; they differ. + expect_equal(expected_ism, 0, tolerance = 1e-10) + expect_gt(res$g_lin[1], 1) +}) + +test_that("pwm.grad.ism max: matches oracle on weak-consensus PSSM", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 15L + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + # Mid-strength PSSM so flipping can plausibly shift argmax. + pssm <- matrix( + c( + 0.5, 0.3, 0.1, 0.1, + 0.1, 0.5, 0.3, 0.1, + 0.1, 0.1, 0.5, 0.3 + ), + L, 4, + byrow = TRUE, + dimnames = list(NULL, c("A", "C", "G", "T")) + ) + + gvtrack.create("g_ism", NULL, "pwm.grad.ism", + pssm = pssm, aggregate = "max", + bidirect = FALSE, strand = 1, extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_ism", sshift = 0, eshift = iter_W - 1) + res <- gextract("g_ism", iterator = 1, intervals = pivot) + + expected <- manual_pwm_grad_ism(pssm, seq_ext, "max", + bidirect = FALSE, strand = 1L, prior = 0.01 + ) + expect_equal(res$g_ism[1], expected, tolerance = 1e-5, ignore_attr = TRUE) + expect_gte(res$g_ism[1], -1e-7) +}) + +# Task 7: GRAD_LSE_ISM tests + +test_that("pwm.grad.ism lse: matches oracle (concentrated softmax)", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 20L + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + head_bases <- strsplit(substr(seq_ext, 1, L), "")[[1]] + pssm <- .consensus_pssm(head_bases) + + gvtrack.create("g_ism_lse", NULL, "pwm.grad.ism", + pssm = pssm, aggregate = "lse", + bidirect = FALSE, strand = 1, extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_ism_lse", sshift = 0, eshift = iter_W - 1) + res <- gextract("g_ism_lse", iterator = 1, intervals = pivot) + + expected <- manual_pwm_grad_ism(pssm, seq_ext, "lse", + bidirect = FALSE, strand = 1L, prior = 0.01 + ) + expect_equal(res$g_ism_lse[1], expected, tolerance = 1e-5, ignore_attr = TRUE) +}) + +test_that("pwm.grad.ism lse: matches oracle on mid-strength PSSM", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 15L + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + pssm <- matrix( + c( + 0.5, 0.3, 0.1, 0.1, + 0.1, 0.5, 0.3, 0.1, + 0.1, 0.1, 0.5, 0.3 + ), + L, 4, + byrow = TRUE, + dimnames = list(NULL, c("A", "C", "G", "T")) + ) + + gvtrack.create("g_ism_lse", NULL, "pwm.grad.ism", + pssm = pssm, aggregate = "lse", + bidirect = FALSE, strand = 1, extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_ism_lse", sshift = 0, eshift = iter_W - 1) + res <- gextract("g_ism_lse", iterator = 1, intervals = pivot) + + expected <- manual_pwm_grad_ism(pssm, seq_ext, "lse", + bidirect = FALSE, strand = 1L, prior = 0.01 + ) + expect_equal(res$g_ism_lse[1], expected, tolerance = 1e-5, ignore_attr = TRUE) +}) + +# Task 8: bidirect for ISM + +test_that("pwm.grad.ism max: bidirect=TRUE matches oracle on asymmetric PSSM", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 20L + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + pssm <- matrix( + c( + 0.6, 0.2, 0.1, 0.1, + 0.1, 0.5, 0.3, 0.1, + 0.1, 0.1, 0.5, 0.3 + ), + L, 4, + byrow = TRUE, + dimnames = list(NULL, c("A", "C", "G", "T")) + ) + + gvtrack.create("g_ism_bid", NULL, "pwm.grad.ism", + pssm = pssm, aggregate = "max", + bidirect = TRUE, extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_ism_bid", sshift = 0, eshift = iter_W - 1) + res <- gextract("g_ism_bid", iterator = 1, intervals = pivot) + + expected <- manual_pwm_grad_ism(pssm, seq_ext, "max", + bidirect = TRUE, prior = 0.01 + ) + expect_equal(res$g_ism_bid[1], expected, tolerance = 1e-5, ignore_attr = TRUE) +}) + +test_that("pwm.grad.ism lse: bidirect=TRUE matches oracle on asymmetric PSSM", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 20L + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + pssm <- matrix( + c( + 0.6, 0.2, 0.1, 0.1, + 0.1, 0.5, 0.3, 0.1, + 0.1, 0.1, 0.5, 0.3 + ), + L, 4, + byrow = TRUE, + dimnames = list(NULL, c("A", "C", "G", "T")) + ) + + gvtrack.create("g_ism_bid_lse", NULL, "pwm.grad.ism", + pssm = pssm, aggregate = "lse", + bidirect = TRUE, extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_ism_bid_lse", sshift = 0, eshift = iter_W - 1) + res <- gextract("g_ism_bid_lse", iterator = 1, intervals = pivot) + + expected <- manual_pwm_grad_ism(pssm, seq_ext, "lse", + bidirect = TRUE, prior = 0.01 + ) + expect_equal(res$g_ism_bid_lse[1], expected, tolerance = 1e-5, ignore_attr = TRUE) +}) + +test_that("pwm.grad.ism max: strand=-1 matches oracle", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 20L + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + head_bases <- strsplit(substr(seq_ext, 1, L), "")[[1]] + fwd_pssm <- .consensus_pssm(head_bases) + rc_target_pssm <- fwd_pssm[L:1, c(4, 3, 2, 1)] + dimnames(rc_target_pssm) <- dimnames(fwd_pssm) + + gvtrack.create("g_ism_minus", NULL, "pwm.grad.ism", + pssm = rc_target_pssm, aggregate = "max", + bidirect = FALSE, strand = -1, extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_ism_minus", sshift = 0, eshift = iter_W - 1) + res <- gextract("g_ism_minus", iterator = 1, intervals = pivot) + + expected <- manual_pwm_grad_ism(rc_target_pssm, seq_ext, "max", + bidirect = FALSE, strand = -1L, prior = 0.01 + ) + expect_equal(res$g_ism_minus[1], expected, tolerance = 1e-5, ignore_attr = TRUE) +}) From b047c04f4416dd0b05499d9ebcdca3783bc85aee Mon Sep 17 00:00:00 2001 From: aviezerl Date: Thu, 7 May 2026 23:36:19 +0300 Subject: [PATCH 08/15] feat(pwm.grad): spatial weighting support 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. --- src/PWMScorer.cpp | 71 ++++++++-------- tests/testthat/helper-pwm-grad-oracle.R | 40 ++++++++- tests/testthat/test-vtrack-pwm-grad.R | 105 ++++++++++++++++++++++++ 3 files changed, 178 insertions(+), 38 deletions(-) diff --git a/src/PWMScorer.cpp b/src/PWMScorer.cpp index 1fa29e6a6..fbd763eb0 100644 --- a/src/PWMScorer.cpp +++ b/src/PWMScorer.cpp @@ -458,15 +458,18 @@ float PWMScorer::score_grad(const std::string& target, // Scan and accumulate. We track all-anchor stats plus excluding-head // stats; the latter let ISM compute max/LSE after a single-base flip in - // O(1) without a second pass. + // O(1) without a second pass. When spatial weighting is on, each anchor's + // score is augmented by a position-dependent log factor; head_fwd/head_rc + // are kept *raw* (no spat) so the ISM flip math can recompute them. float lse = -std::numeric_limits::infinity(); float lse_no_head = -std::numeric_limits::infinity(); float best = -std::numeric_limits::infinity(); float best_no_head = -std::numeric_limits::infinity(); size_t argmax = i_min; - float head_fwd = -std::numeric_limits::infinity(); - float head_rc = -std::numeric_limits::infinity(); - float head_score = -std::numeric_limits::infinity(); + float head_fwd = -std::numeric_limits::infinity(); // raw, no spat + float head_rc = -std::numeric_limits::infinity(); // raw, no spat + float head_score = -std::numeric_limits::infinity(); // with spat + float head_spat_log = 0.0f; bool any_finite = false; for (size_t i = i_min; i <= i_max; ++i) { @@ -481,27 +484,32 @@ float PWMScorer::score_grad(const std::string& target, rc = score_reverse_original(m_pssm, target, i, m_strand); } - // Per-anchor score: comb (bidirect) or single-strand. - float s; + // Per-anchor raw score: comb (bidirect) or single-strand. + float s_raw; if (bidirect) { if (!std::isfinite(fwd) && !std::isfinite(rc)) { - s = -std::numeric_limits::infinity(); + s_raw = -std::numeric_limits::infinity(); } else if (!std::isfinite(fwd)) { - s = rc; + s_raw = rc; } else if (!std::isfinite(rc)) { - s = fwd; + s_raw = fwd; } else { - s = fwd; - log_sum_log(s, rc); + s_raw = fwd; + log_sum_log(s_raw, rc); } } else { - s = use_fwd ? fwd : rc; + s_raw = use_fwd ? fwd : rc; } + // Apply spatial log factor (= 0 when m_use_spat is false). + const float spat_log = get_spatial_log_factor(i - i_min); + const float s = std::isfinite(s_raw) ? (s_raw + spat_log) : s_raw; + if (i == head_idx) { - head_fwd = fwd; - head_rc = rc; - head_score = s; + head_fwd = fwd; + head_rc = rc; + head_score = s; + head_spat_log = spat_log; } if (std::isfinite(s)) { @@ -625,22 +633,24 @@ float PWMScorer::score_grad(const std::string& target, } } - // Combined per-anchor score after flip. - float s_alt; + // Combined per-anchor raw score after flip, then apply head's spat_log. + float s_alt_raw; if (bidirect) { if (!std::isfinite(fwd_alt) && !std::isfinite(rc_alt)) { - s_alt = -std::numeric_limits::infinity(); + s_alt_raw = -std::numeric_limits::infinity(); } else if (!std::isfinite(fwd_alt)) { - s_alt = rc_alt; + s_alt_raw = rc_alt; } else if (!std::isfinite(rc_alt)) { - s_alt = fwd_alt; + s_alt_raw = fwd_alt; } else { - s_alt = fwd_alt; - log_sum_log(s_alt, rc_alt); + s_alt_raw = fwd_alt; + log_sum_log(s_alt_raw, rc_alt); } } else { - s_alt = use_fwd ? fwd_alt : rc_alt; + s_alt_raw = use_fwd ? fwd_alt : rc_alt; } + const float s_alt = std::isfinite(s_alt_raw) + ? (s_alt_raw + head_spat_log) : s_alt_raw; // Aggregate after flip: only the head anchor's score changes. float f_alt; @@ -1163,20 +1173,11 @@ float PWMScorer::score_interval(const GInterval& interval, const GenomeChromKey& i_min = 0; } - // Spatial weighting for gradient modes is deferred to a later task. - // Surface a clear error rather than silently using non-spatial code. - if (m_use_spat && - (m_mode == GRAD_LSE || m_mode == GRAD_MAX || - m_mode == GRAD_LSE_ISM || m_mode == GRAD_MAX_ISM)) { - rdb::verror("PWMScorer: spatial weighting not yet supported for gradient modes"); - } - - // Gradient modes (single-strand fwd) take a dedicated answer path. + // Gradient modes route through score_grad regardless of spatial. // They share the iterator clamping (i_min..i_max) above so that the // argmax convention matches pwm.max exactly. - if (!m_use_spat && - (m_mode == GRAD_MAX || m_mode == GRAD_LSE || - m_mode == GRAD_MAX_ISM || m_mode == GRAD_LSE_ISM)) { + if (m_mode == GRAD_MAX || m_mode == GRAD_LSE || + m_mode == GRAD_MAX_ISM || m_mode == GRAD_LSE_ISM) { const bool lse_aggregate = (m_mode == GRAD_LSE || m_mode == GRAD_LSE_ISM); const bool ism = (m_mode == GRAD_MAX_ISM || m_mode == GRAD_LSE_ISM); return score_grad(target, i_min, i_max, motif_len, lse_aggregate, ism); diff --git a/tests/testthat/helper-pwm-grad-oracle.R b/tests/testthat/helper-pwm-grad-oracle.R index 0de6904e3..791118601 100644 --- a/tests/testthat/helper-pwm-grad-oracle.R +++ b/tests/testthat/helper-pwm-grad-oracle.R @@ -95,20 +95,48 @@ stop("strand must be +1 or -1 when bidirect = FALSE") } +# Apply a per-anchor spatial log-factor to a vector of per-anchor scores. +# `spat_factor` is a vector of multiplicative factors (in linear scale); we +# convert to log and add. `spat_bin_size` controls how many anchors share each +# factor (matches the engine's get_spatial_log_factor convention). When +# `spat_factor` is NULL or empty, scores are returned unchanged. +.apply_spat <- function(scores, spat_factor, spat_bin_size = 1L) { + if (is.null(spat_factor) || length(spat_factor) == 0) { + return(scores) + } + n <- length(scores) + bin_sz <- max(1L, as.integer(spat_bin_size)) + bins <- ((seq_len(n) - 1L) %/% bin_sz) + 1L + bins <- pmin(bins, length(spat_factor)) + log_factors <- log(pmax(1e-30, spat_factor))[bins] + scores + log_factors +} + # Linearized gradient at p = interval start, head anchor, column 0 (= base b_p). # # Returns a single non-negative scalar. For an all-N seq or otherwise invalid # head base, returns 0. manual_pwm_grad <- function(pssm, seq, aggregate = c("lse", "max"), - bidirect = TRUE, strand = 1L, prior = 0.01) { + bidirect = TRUE, strand = 1L, prior = 0.01, + spat_factor = NULL, spat_bin_size = 1L) { aggregate <- match.arg(aggregate) L <- nrow(pssm) # Normalized log PSSM (used for diff lookups). M <- log(.normalize_pssm(pssm, prior = prior)) - # Per-anchor scores. + # Per-anchor scores (with spatial weighting applied). scores <- .oracle_anchor_scores(seq, pssm, prior = prior) + scores$fwd <- .apply_spat(scores$fwd, spat_factor, spat_bin_size) + scores$rc <- .apply_spat(scores$rc, spat_factor, spat_bin_size) + # Recompute comb after spatial weighting (spatial is the same for fwd and + # rc at a given anchor, so log(exp(fwd_spat) + exp(rc_spat)) = log(exp(fwd) + # + exp(rc)) + spat, but recomputing keeps the code symmetric and handles + # the all-Inf cases cleanly). + scores$comb <- vapply(seq_along(scores$fwd), function(i) { + log_sum_exp(c(scores$fwd[i], scores$rc[i])) + }, numeric(1)) + s <- .pick_scores(scores, bidirect, strand) f <- .aggregate_scores(s, aggregate) @@ -153,11 +181,17 @@ manual_pwm_grad <- function(pssm, seq, aggregate = c("lse", "max"), # g = f_actual - min over b' != b_p of f(seq with seq[1] := b'). # If b_p is not ACGT, returns 0. manual_pwm_grad_ism <- function(pssm, seq, aggregate = c("lse", "max"), - bidirect = TRUE, strand = 1L, prior = 0.01) { + bidirect = TRUE, strand = 1L, prior = 0.01, + spat_factor = NULL, spat_bin_size = 1L) { aggregate <- match.arg(aggregate) aggregate_for_seq <- function(s_seq) { sc <- .oracle_anchor_scores(s_seq, pssm, prior = prior) + sc$fwd <- .apply_spat(sc$fwd, spat_factor, spat_bin_size) + sc$rc <- .apply_spat(sc$rc, spat_factor, spat_bin_size) + sc$comb <- vapply(seq_along(sc$fwd), function(i) { + log_sum_exp(c(sc$fwd[i], sc$rc[i])) + }, numeric(1)) s <- .pick_scores(sc, bidirect, strand) .aggregate_scores(s, aggregate) } diff --git a/tests/testthat/test-vtrack-pwm-grad.R b/tests/testthat/test-vtrack-pwm-grad.R index 954238cff..e39458cba 100644 --- a/tests/testthat/test-vtrack-pwm-grad.R +++ b/tests/testthat/test-vtrack-pwm-grad.R @@ -743,3 +743,108 @@ test_that("pwm.grad.ism max: strand=-1 matches oracle", { ) expect_equal(res$g_ism_minus[1], expected, tolerance = 1e-5, ignore_attr = TRUE) }) + +# Task 9: spatial weighting + +test_that("pwm.grad: spat_factor = all-1 matches non-spatial", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 20L + pivot <- gintervals(1, 200, 201) + + head_bases <- strsplit( + toupper(gseq.extract(gintervals(1, 200, 200 + L))), "" + )[[1]] + pssm <- .consensus_pssm(head_bases) + + spat_one <- rep(1, iter_W) + + gvtrack.create("g_no_spat", NULL, "pwm.grad", + pssm = pssm, aggregate = "lse", + bidirect = FALSE, strand = 1, extend = TRUE, prior = 0.01 + ) + gvtrack.create("g_spat_one", NULL, "pwm.grad", + pssm = pssm, spat_factor = spat_one, spat_bin = 1L, + aggregate = "lse", bidirect = FALSE, strand = 1, + extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_no_spat", sshift = 0, eshift = iter_W - 1) + gvtrack.iterator("g_spat_one", sshift = 0, eshift = iter_W - 1) + + res <- gextract(c("g_no_spat", "g_spat_one"), + iterator = 1, intervals = pivot + ) + expect_equal(res$g_no_spat[1], res$g_spat_one[1], + tolerance = 1e-5, ignore_attr = TRUE + ) +}) + +test_that("pwm.grad lse: non-trivial spat_factor matches oracle", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 10L + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + head_bases <- strsplit(substr(seq_ext, 1, L), "")[[1]] + pssm <- .consensus_pssm(head_bases) + + # Down-weight late anchors so the LSE concentrates more on early ones, + # which boosts the head's softmax weight. + spat_factor <- c(rep(1, 5), rep(0.01, 5)) + + gvtrack.create("g_spat", NULL, "pwm.grad", + pssm = pssm, spat_factor = spat_factor, spat_bin = 1L, + aggregate = "lse", bidirect = FALSE, strand = 1, + extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_spat", sshift = 0, eshift = iter_W - 1) + res <- gextract("g_spat", iterator = 1, intervals = pivot) + + expected <- manual_pwm_grad(pssm, seq_ext, "lse", + bidirect = FALSE, strand = 1L, prior = 0.01, + spat_factor = spat_factor, spat_bin_size = 1L + ) + expect_equal(res$g_spat[1], expected, tolerance = 1e-5, ignore_attr = TRUE) +}) + +test_that("pwm.grad.ism lse: spat_factor matches oracle", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 10L + pivot <- gintervals(1, 200, 201) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + + pssm <- matrix( + c( + 0.5, 0.3, 0.1, 0.1, + 0.1, 0.5, 0.3, 0.1, + 0.1, 0.1, 0.5, 0.3 + ), + L, 4, + byrow = TRUE, + dimnames = list(NULL, c("A", "C", "G", "T")) + ) + + spat_factor <- c(rep(1, 5), rep(0.05, 5)) + + gvtrack.create("g_ism_spat", NULL, "pwm.grad.ism", + pssm = pssm, spat_factor = spat_factor, spat_bin = 1L, + aggregate = "lse", bidirect = FALSE, strand = 1, + extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_ism_spat", sshift = 0, eshift = iter_W - 1) + res <- gextract("g_ism_spat", iterator = 1, intervals = pivot) + + expected <- manual_pwm_grad_ism(pssm, seq_ext, "lse", + bidirect = FALSE, strand = 1L, prior = 0.01, + spat_factor = spat_factor, spat_bin_size = 1L + ) + expect_equal(res$g_ism_spat[1], expected, tolerance = 1e-5, ignore_attr = TRUE) +}) From e5ac674b3e077f85adf05ed7ce4087b30c9ed38c Mon Sep 17 00:00:00 2001 From: aviezerl Date: Thu, 7 May 2026 23:37:38 +0300 Subject: [PATCH 09/15] test(pwm.grad): NA, composition, multi-chrom, gscreen / gsummary edge cases --- tests/testthat/test-vtrack-pwm-grad.R | 117 ++++++++++++++++++++++++++ 1 file changed, 117 insertions(+) diff --git a/tests/testthat/test-vtrack-pwm-grad.R b/tests/testthat/test-vtrack-pwm-grad.R index e39458cba..1679e9e1f 100644 --- a/tests/testthat/test-vtrack-pwm-grad.R +++ b/tests/testthat/test-vtrack-pwm-grad.R @@ -848,3 +848,120 @@ test_that("pwm.grad.ism lse: spat_factor matches oracle", { ) expect_equal(res$g_ism_spat[1], expected, tolerance = 1e-5, ignore_attr = TRUE) }) + +# Task 10: NA / composition / edge-case tests + +test_that("pwm.grad: interval shorter than L returns NA", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 5L + pivot <- gintervals(1, 200, 201) + pssm <- .consensus_pssm(c("A", "C", "G", "T", "A")) + + # eshift only adds 2 bp -> total interval length = 3, < L = 5 -> NA. + gvtrack.create("g_short", NULL, "pwm.grad", + pssm = pssm, aggregate = "lse", + bidirect = FALSE, strand = 1, extend = FALSE, prior = 0.01 + ) + gvtrack.iterator("g_short", sshift = 0, eshift = 2) + res <- gextract("g_short", iterator = 1, intervals = pivot) + expect_true(is.na(res$g_short[1])) +}) + +test_that("pwm.grad: composes with pwm in gextract", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 20L + pivot_block <- gintervals(1, 200, 210) + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + head_bases <- strsplit(substr(seq_ext, 1, L), "")[[1]] + pssm <- .consensus_pssm(head_bases) + + gvtrack.create("score", NULL, "pwm", + pssm = pssm, bidirect = FALSE, strand = 1, + extend = TRUE, prior = 0.01 + ) + gvtrack.create("g", NULL, "pwm.grad", + pssm = pssm, aggregate = "lse", + bidirect = FALSE, strand = 1, extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("score", sshift = 0, eshift = iter_W - 1) + gvtrack.iterator("g", sshift = 0, eshift = iter_W - 1) + + res <- gextract(c("score", "g"), iterator = 1, intervals = pivot_block) + expect_equal(nrow(res), 10L) # 10 1-bp positions + expect_true(all(res$g >= -1e-7, na.rm = TRUE)) + # Score and gradient should both vary across positions. + expect_gt(length(unique(res$score)), 1L) +}) + +test_that("pwm.grad: gscreen filter on gradient", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 20L + seq_ext <- toupper(gseq.extract(gintervals(1, 200, 200 + iter_W + L - 1))) + head_bases <- strsplit(substr(seq_ext, 1, L), "")[[1]] + pssm <- .consensus_pssm(head_bases) + + gvtrack.create("g", NULL, "pwm.grad", + pssm = pssm, aggregate = "lse", + bidirect = FALSE, strand = 1, extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g", sshift = 0, eshift = iter_W - 1) + + out <- gscreen("g > 0.1", iterator = 1, intervals = gintervals(1, 200, 220)) + expect_s3_class(out, "data.frame") + # All returned intervals must satisfy the predicate. + if (nrow(out) > 0) { + vals <- gextract("g", iterator = 1, intervals = out) + expect_true(all(vals$g > 0.1, na.rm = TRUE)) + } +}) + +test_that("pwm.grad: multi-chromosome intervals", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 10L + pssm <- .consensus_pssm(c("A", "C", "G")) + + gvtrack.create("g", NULL, "pwm.grad", + pssm = pssm, aggregate = "lse", + bidirect = FALSE, strand = 1, extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g", sshift = 0, eshift = iter_W - 1) + + intervals <- rbind( + gintervals(1, 200, 210), + gintervals(2, 200, 210) + ) + res <- gextract("g", iterator = 1, intervals = intervals) + expect_equal(nrow(res), 20L) + expect_true(all(c(1, 2) %in% res$chrom1 | + c("chr1", "chr2") %in% res$chrom)) +}) + +test_that("pwm.grad: gsummary aggregation runs without error", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + + L <- 3L + iter_W <- 10L + pssm <- .consensus_pssm(c("A", "C", "G")) + + gvtrack.create("g", NULL, "pwm.grad", + pssm = pssm, aggregate = "max", + bidirect = FALSE, strand = 1, extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g", sshift = 0, eshift = iter_W - 1) + + s <- gsummary("g", iterator = 1, intervals = gintervals(1, 200, 250)) + expect_true(is.numeric(s)) + expect_true(length(s) >= 5L) +}) From 14cf460da35a9b41e3514d63476ed695cc8365eb Mon Sep 17 00:00:00 2001 From: aviezerl Date: Thu, 7 May 2026 23:39:51 +0300 Subject: [PATCH 10/15] docs(pwm.grad): NEWS entry for pwm.grad / pwm.grad.ism --- NEWS.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/NEWS.md b/NEWS.md index 036dfbcbf..706c87f12 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # misha (development version) +* Added `pwm.grad` and `pwm.grad.ism` virtual tracks: per-bp gradient / saliency of the PWM aggregate (LSE or MAX) at the iterator interval start. `pwm.grad` is the linearized DeepLIFT-style attribution; `pwm.grad.ism` is in-silico mutagenesis. + # misha 5.6.28 * `gtrack.copy()` gained a `db` argument to copy tracks across databases, and an `overwrite` argument to replace existing destinations. Format conversion (per-chromosome <-> indexed) and chromosome-order remap are handled automatically. Multi-track input is also supported. From 797990eab2d8436dd06a12049df0c21ca523bbde Mon Sep 17 00:00:00 2001 From: aviezerl Date: Thu, 7 May 2026 23:44:48 +0300 Subject: [PATCH 11/15] docs(vignette): PWM virtual tracks (existing family + pwm.grad) 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. --- _pkgdown.yml | 9 + vignettes/PWM-Functions.Rmd | 465 ++++++++++++++++++++++++++++++++++++ 2 files changed, 474 insertions(+) create mode 100644 vignettes/PWM-Functions.Rmd diff --git a/_pkgdown.yml b/_pkgdown.yml index 010ab1188..3d99a8a86 100755 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -1,3 +1,12 @@ +articles: + - title: Articles + contents: + - Misha-Basics + - PWM-Functions + - Manual + - Genomes + - Database-Formats + reference: - title: Track functions desc: Functions to manipulate and create tracks diff --git a/vignettes/PWM-Functions.Rmd b/vignettes/PWM-Functions.Rmd new file mode 100644 index 000000000..781e073d7 --- /dev/null +++ b/vignettes/PWM-Functions.Rmd @@ -0,0 +1,465 @@ +--- +title: "PWM Virtual Tracks" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{PWM Virtual Tracks} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.width = 7, fig.height = 4 +) +``` + +This vignette covers misha's PWM virtual track family: the existing +aggregating scorers (`pwm`, `pwm.max`, `pwm.max.pos`, `pwm.count`), the +edit-distance variants, and the per-position attribution tracks +(`pwm.grad`, `pwm.grad.ism`). + +The aggregating sections (`pwm`, `pwm.max`, ...) are well documented on +the `gvtrack.create` help page, so the focus here is on the mental +model and on examples a reader can adapt to their own database. The +math for the gradient tracks is worked out in more detail because they +are new and there is no existing reference. + +The math chunks below run as plain R (no genome database required) so +they always evaluate. The genome-track examples are shown with +`eval = FALSE` because they need a database; the patterns transfer +unchanged to any misha DB. + +## What is a PSSM? + +A position-specific scoring matrix (PSSM, also "PWM") is a per-column +distribution over A/C/G/T. In misha the input is an `L x 4` matrix: + +- **Rows**: motif columns, in 5' -> 3' order (one row per nucleotide + position in the motif). +- **Columns**: bases, named `A`, `C`, `G`, `T` (any extra columns are + ignored). + +The raw matrix typically holds frequencies or counts. Internally misha +adds a Dirichlet pseudocount (`prior`, default `0.01`), normalizes each +row to sum to 1, and takes the natural log. The result is a matrix of +log-probabilities `M[c, b]`: the log-likelihood that base `b` appears +at column `c` of the motif. + +```{r} +# A 3-column PSSM with consensus ACG. +pssm <- matrix( + c( + 0.97, 0.01, 0.01, 0.01, + 0.01, 0.97, 0.01, 0.01, + 0.01, 0.01, 0.97, 0.01 + ), + nrow = 3, ncol = 4, byrow = TRUE, + dimnames = list(NULL, c("A", "C", "G", "T")) +) +pssm +``` + +```{r} +# Apply prior + per-row normalization, then log. +prior <- 0.01 +M <- log((pssm + prior) / rowSums(pssm + prior)) +round(M, 3) +``` + +`M[1, "A"]` is close to 0 (high probability), while every off-consensus +entry is a large negative number. The `prior` keeps log-zero rows +finite: with `prior = 0` and a strict frequency matrix, missing bases +become `-Inf` and dominate any sum. `0.01` is a conservative default; a +raw count matrix often wants something larger. + +## Per-anchor score + +Given a sequence `s`, an *anchor* is a position `a` at which the motif +is laid down. The per-anchor (single-strand, forward) score is + +$$ +\text{score}_a \;=\; \sum_{c=0}^{L-1} M[c,\, s_{a+c}]. +$$ + +It is just a column-by-column lookup followed by a sum. Units are log +probability; values are <= 0 by construction (each term is the log of a +probability in [0, 1]). + +```{r} +# Score the anchor at position 1 of the sequence "ACGTTT". +seq <- "ACGTTT" +bases <- strsplit(seq, "")[[1]] +score_anchor <- function(M, bases, anchor) { + L <- nrow(M) + sum(vapply(seq_len(L), function(c) M[c, bases[anchor + c - 1]], numeric(1))) +} +score_anchor(M, bases, 1) # "ACG" - perfect consensus +score_anchor(M, bases, 2) # "CGT" - off-consensus +``` + +The first anchor matches the consensus and scores close to 0; the +second is well below. + +For a window of length `W >= L`, we get `W - L + 1` per-anchor scores. + +```{r} +all_anchor_scores <- function(M, seq) { + bases <- strsplit(seq, "")[[1]] + W <- length(bases) + L <- nrow(M) + vapply(seq_len(W - L + 1), function(a) score_anchor(M, bases, a), numeric(1)) +} +all_anchor_scores(M, seq) +``` + +## Aggregating over an interval + +The PWM virtual tracks each reduce these per-anchor scores to a single +scalar per iterator interval. + +### `pwm.max`: best anchor + +`pwm.max` returns the maximum per-anchor score in the window. Easiest +to interpret ("how good is the best motif hit here?"), but discards +information about secondary hits. + +```{r} +max(all_anchor_scores(M, seq)) +``` + +### `pwm` (LSE): smooth aggregate + +`pwm` returns the log-sum-exp of all per-anchor scores: + +$$ +f_{\text{LSE}} \;=\; \log \sum_a \exp(\text{score}_a). +$$ + +This is dominated by the best anchor when one anchor clearly wins, but +also adds support for many mediocre anchors (an "OR" over noisy +matches). LSE is also smooth: small changes in the input give small +changes in the output, which is friendlier for downstream filtering and +gradient-based attribution. + +```{r} +log_sum_exp <- function(x) { + m <- max(x) + m + log(sum(exp(x - m))) +} +s <- all_anchor_scores(M, seq) +log_sum_exp(s) +``` + +For a single dominant hit `pwm` and `pwm.max` are nearly equal; for a +window with many weak similar-strength hits `pwm` can be substantially +larger. + +### `pwm.max.pos`: argmax position + +`pwm.max.pos` returns the 1-based offset of the argmax anchor inside +the iterator interval, signed by strand when `bidirect = TRUE`. Useful +when you want to know *where* the best hit sits within a region. + +### `pwm.count`: thresholded hit count + +`pwm.count` requires `score.thresh`. It counts anchors whose score is +at or above `score.thresh`. With `bidirect = TRUE` the union over +strands is taken at each genomic start, so palindromic hits count once. + +```{r, eval = FALSE} +gvtrack.create("motif_count", NULL, "pwm.count", + pssm = pssm, score.thresh = -2, bidirect = TRUE +) +gvtrack.iterator("motif_count", sshift = 0, eshift = 500) +gextract("motif_count", iterator = 500, intervals = my_intervals) +``` + +### Why LSE rather than just MAX? + +For any aggregation over a long window, LSE has three practical +advantages over MAX: + +1. **Integrates many anchors.** A region with five mediocre hits has a + higher LSE than one with a single mediocre hit, but the same MAX. +2. **Smooth.** Differentiable everywhere; small input changes -> small + output changes. This matters for gradient attribution + (next section). +3. **Numerically stable when implemented with the standard + `m + log(sum(exp(x - m)))` trick.** + +When you only care about "is there a strong hit anywhere in this +window?", `pwm.max` is fine. For typical "motif energy" scoring across +500 bp - 1 kb windows, `pwm` (LSE) is the better default. + +## Common parameters + +These mirror what `?gvtrack.create` documents; one-liner each: + +- **`bidirect`** (default `TRUE`): scan both strands and combine per + genomic start. With `bidirect = FALSE`, only the strand given by + `strand` is scanned. +- **`strand`** (`+1` / `-1`): only used when `bidirect = FALSE`. +- **`prior`** (default `0.01`): Dirichlet pseudocount added to PSSM + rows before normalization. Raise it for noisy / count-based PSSMs; + set to 0 to disable. +- **`extend`** (default `TRUE`): pad the fetched sequence so anchors + near the right edge of the iterator interval still get a full + window. Without extension, an interval of length `W` only yields + `W - L + 1` anchors. +- **`spat_factor`** + **`spat_bin`**: optional position-dependent log + weights. `spat_factor` is a positive numeric vector; `spat_bin` is + the bin width in bp. Each anchor's score is multiplied by the log + factor for its bin. Use this to weight central bases of a window + more than flanks. +- **`score.thresh`**: required by `pwm.count` and the edit-distance + variants. Anchors with score >= threshold are "hits". Rejected by + `pwm.grad` and `pwm.grad.ism`. + +## Edit-distance variants + +These ask "how many single-base edits would it take to make / break a +motif match?". They share the same PSSM machinery but answer a +discrete combinatorial question instead of returning a continuous +score. + +- **`pwm.edit_distance`**: minimum number of edits (substitutions, and + optionally indels via `max_indels`) needed across all windows in the + iterator interval to push the score across `score.thresh`. + `direction = "above"` (default) targets making a match; + `direction = "below"` targets disrupting one. +- **`pwm.max.edit_distance`**: same idea, but always evaluated at the + argmax window of `pwm.max`. +- **`pwm.n_mutations`**: number of independent single-base + substitutions that each cross `score.thresh` on their own. Returns 0 + if the threshold is already met, NA if no single substitution + suffices. + +There is also an LSE-aware family (`pwm.edit_distance.lse`, +`pwm.edit_distance.lse.pos`) that targets the LSE aggregate rather +than per-window scores. + +See `?gvtrack.create` for the full parameter list and worked examples. + +## Per-position attribution: `pwm.grad` and `pwm.grad.ism` + +The aggregating tracks above all collapse a window into one scalar. +`pwm.grad` and `pwm.grad.ism` answer a different question: **at this +specific base pair, how much does the local sequence support the PWM +aggregate?** The output is one non-negative number per iterator +interval, computed at the anchor that starts at `interval.start`. + +Combined with a `iterator = 1` step and `gvtrack.iterator(eshift = +W - 1)`, this gives a per-bp saliency / attribution track you can plot +or combine with conservation, ATAC, methylation, etc. + +### Motivation + +Two use cases: + +1. **Saliency.** "Of the base pairs in this window, which ones are + actually carrying the motif signal?" - a poor man's + integrated-gradient / DeepLIFT. +2. **Single-bp coupling.** Combine PWM evidence with another per-bp + signal (conservation, ChIP, etc.) at single-base resolution. With + `pwm` you only get one number per window; with `pwm.grad` you get + one per bp. + +### Math intuition (linearized, `pwm.grad`) + +Think of each base as a one-hot vector `x_p` over A/C/G/T. The score +`score_a` is linear in `x_p` whenever `a` contains `p`, with gradient +`M[p - a, b]`. With pivot fixed at the iterator start `p = +interval.start`, only one anchor (`a = p`) contains `p`, and only +column 0 of the PSSM is involved. + +For the LSE aggregate, the gradient is the softmax-weighted column-0 +score, contracted against `(x_p - x_p^{\text{worst}})` where the +reference is the worst possible base at column 0: + +$$ +g_{\text{LSE}}(p) \;=\; w_p \cdot \big(M[0,\, b_p] - \min_b M[0,\, b]\big), +\qquad w_p = \exp(\text{score}_p - f_{\text{LSE}}). +$$ + +Three useful sanity points: + +- **Non-negative.** `M[0, b_p] - min_b M[0, b]` >= 0 by construction; + `w_p` in [0, 1]. +- **Bounded by per-column information.** The diff term equals the + column's information content under the "actual minus worst" baseline. +- **Modulated by softmax weight.** If the head anchor dominates the + LSE (`w_p ≈ 1`), the gradient saturates to the diff. If the LSE + spreads over many anchors (`w_p` small), the gradient is small even + if the head anchor is itself a good match - the *attribution* of + this particular base to the aggregate is diluted. + +For the MAX aggregate, the gradient picks up the diff iff the head +anchor is the argmax: + +$$ +g_{\text{MAX}}(p) \;=\; +\begin{cases} +M[0,\, b_p] - \min_b M[0,\, b] & \text{if } a^* = p,\\ +0 & \text{otherwise.} +\end{cases} +$$ + +### Math intuition (in-silico mutagenesis, `pwm.grad.ism`) + +The ISM gradient is "actual minus best replacement": + +$$ +g_{\text{ISM}}(p) \;=\; f(\text{seq}) \;-\; \min_{b' \neq b_p} f\big(\text{seq with } s_p := b'\big). +$$ + +Always >= 0. With pivot fixed at `interval.start`, only `score_p` +changes when we flip `s_p`, so the update is O(1) per flip and ISM +costs only ~3-4x more than the linearized form. + +### When `pwm.grad` and `pwm.grad.ism` differ + +The two definitions agree exactly when the argmax is stable under +flips of `s_p`, and agree to high precision when LSE concentrates on +the head anchor. They diverge in two cases: + +1. **Argmax shifts.** For MAX, if flipping `s_p` makes a *different* + anchor become the argmax, ISM accounts for that ("the second-best + hit kicks in") and the linearized form doesn't. +2. **Curvature in LSE.** Linearization is a first-order approximation + around the actual sequence; ISM evaluates the aggregate at three + discrete neighbours and takes the worst. For sharply non-linear + regions of the LSE landscape they differ slightly, but in practice + the disagreement is small. + +In short: **`pwm.grad` is cheaper and almost always good enough; +`pwm.grad.ism` is more conservative around argmax-shift edge cases.** + +### Worked numerical example + +Consider an `L = 3` PSSM and a 6-bp window starting with consensus +"ACG": + +```{r} +pssm <- matrix( + c( + 0.97, 0.01, 0.01, 0.01, + 0.01, 0.97, 0.01, 0.01, + 0.01, 0.01, 0.97, 0.01 + ), + nrow = 3, ncol = 4, byrow = TRUE, + dimnames = list(NULL, c("A", "C", "G", "T")) +) +M <- log((pssm + 0.01) / rowSums(pssm + 0.01)) + +seq <- "ACGTTT" +s <- all_anchor_scores(M, seq) # 4 per-anchor scores +s +``` + +Linearized LSE gradient at `p = 1`: + +```{r} +f_lse <- log_sum_exp(s) +score_p <- s[1] +w_p <- exp(score_p - f_lse) # softmax weight at the head +b_p <- "A" # head base +diff_col0 <- M[1, b_p] - min(M[1, ]) +w_p +diff_col0 +g_lse <- w_p * diff_col0 +g_lse +``` + +`w_p ≈ 1` because the head anchor "ACG" is the dominant LSE +contributor; `diff_col0` is the column-0 information content. The +result is essentially the column's entropy gap. + +ISM gradient: flip `seq[1]` to each of C/G/T and recompute the LSE. + +```{r} +ism_lse <- function(M, seq) { + s <- all_anchor_scores(M, seq) + log_sum_exp(s) +} +flips <- setdiff(c("A", "C", "G", "T"), substr(seq, 1, 1)) +f_alt <- vapply(flips, function(b) { + ism_lse(M, paste0(b, substr(seq, 2, nchar(seq)))) +}, numeric(1)) +f_alt +ism_lse(M, seq) - min(f_alt) +``` + +The two values agree to within numerical noise; this is the typical +case. + +### Genome-track recipe + +Here is the canonical pattern: extract `pwm` and `pwm.grad` over a 1 +kb iterator window, then walk the per-bp gradient with `iterator = 1`. + +```{r, eval = FALSE} +W <- 500 # window length (bp) used for both aggregate and gradient + +pssm <- matrix(...) # your 4-col PSSM, dimnames A/C/G/T + +gvtrack.create("score", NULL, "pwm", pssm = pssm) +gvtrack.create("grad", NULL, "pwm.grad", pssm = pssm, aggregate = "lse") +gvtrack.create("grad_ism", NULL, "pwm.grad.ism", pssm = pssm, aggregate = "lse") + +# Each per-bp value uses a window starting at this bp and extending W bp right. +gvtrack.iterator("score", sshift = 0, eshift = W - 1) +gvtrack.iterator("grad", sshift = 0, eshift = W - 1) +gvtrack.iterator("grad_ism", sshift = 0, eshift = W - 1) + +regions <- gintervals(1, 100000, 110000) +out <- gextract(c("score", "grad", "grad_ism"), + intervals = regions, iterator = 1 +) +head(out) +``` + +Common downstream patterns: + +- **Per-bp coupling**: `out$grad * out$other_track` to get an + attribution-weighted version of `other_track`. +- **Argmax-shift hunting**: `out$grad_ism - out$grad` highlights bp + where ISM and linearization diverge, i.e. where the argmax flips + under a single substitution. +- **Smoothing**: a `gtrack.smooth` over `pwm.grad` gives a soft + saliency band over the canonical motif location. + +### `bidirect = TRUE` for asymmetric PSSMs + +For palindromic PSSMs the rc strand contributes column-0 of the rc +PSSM (= column `L-1` of the fwd PSSM, with complement). With +`bidirect = TRUE` the gradient uses a per-anchor strand softmax to +combine the two contributions; this is handled internally and you +don't need to touch anything beyond `bidirect = TRUE`. For asymmetric +PSSMs, the bidirect gradient at a bp can pick up signal from either +strand depending on which one wins locally. + +## Performance notes + +In rough numbers for an `L = 12` PSSM, `W = 500` bp window, scanned at +1 bp resolution: + +- `pwm.grad` is roughly comparable to `pwm` (both do the same scan; + the gradient post-step is O(1) per anchor). +- `pwm.grad.ism` is ~3-4x `pwm.grad` (three flips x O(1) per + position, with a small constant for log-space arithmetic). + +For genome-wide attribution, you almost never want `iterator = 1` +across the whole genome. Common patterns: + +- Use a sparse intervals set as iterator (e.g. peaks, TSS windows) and + `iterator = 1` only inside those; gives single-bp attribution where + it matters and skips empty space. +- Use a coarser iterator step (e.g. `iterator = 10` or `iterator = + 100`) for whole-genome scans; gradient is still meaningful, just + averaged over the step. + +When in doubt, profile on one chromosome before committing to a +genome-wide scan: the expensive part is reading sequence from disk, +not the gradient computation itself. From 873ef76b7382eaff20453c9f5686144b42d7c4fd Mon Sep 17 00:00:00 2001 From: aviezerl Date: Fri, 8 May 2026 07:57:26 +0300 Subject: [PATCH 12/15] perf(pwm.grad): per-anchor sliding cache (~5x speedup) 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. --- src/PWMScorer.cpp | 119 ++++++++++++++++++++++++++++++++++++++++------ src/PWMScorer.h | 34 +++++++++++++ 2 files changed, 139 insertions(+), 14 deletions(-) diff --git a/src/PWMScorer.cpp b/src/PWMScorer.cpp index fbd763eb0..a0dbc24be 100644 --- a/src/PWMScorer.cpp +++ b/src/PWMScorer.cpp @@ -176,6 +176,7 @@ void PWMScorer::invalidate_cache() m_slide.valid = false; m_slide.stride = 0; m_spat_slide.valid = false; + m_grad_slide.clear(); } // Get spatial log factor for a given position index @@ -403,6 +404,7 @@ float PWMScorer::score_without_spatial(const std::string& target, int64_t motif_ float PWMScorer::score_grad(const std::string& target, size_t i_min, size_t i_max, size_t motif_length, + int chromid, int64_t interval_start, int64_t interval_end, bool lse_aggregate, bool ism) { // Precondition: only the score_thresh / spatial restrictions need not be @@ -461,6 +463,98 @@ float PWMScorer::score_grad(const std::string& target, // O(1) without a second pass. When spatial weighting is on, each anchor's // score is augmented by a position-dependent log factor; head_fwd/head_rc // are kept *raw* (no spat) so the ISM flip math can recompute them. + // ---- Slide-or-seed phase -------------------------------------------------- + // + // Populate m_grad_slide.fwd / m_grad_slide.rc so they hold per-anchor RAW + // scores in target-index order [i_min, i_max]. On consecutive iterator + // steps with the same chrom / strand / window geometry, slide the deque + // (pop trailing anchors, score & push leading anchors) instead of doing a + // full O(W*L) rescan; otherwise seed from scratch. + const size_t W = (i_min <= i_max) ? (i_max - i_min + 1) : 0; + bool can_slide = false; + size_t stride = 0; + if (m_grad_slide.valid && + m_grad_slide.chromid == chromid && + m_grad_slide.strand_mode == m_strand && + m_grad_slide.i_min == i_min && + m_grad_slide.i_max == i_max && + m_grad_slide.fwd.size() == W && + m_grad_slide.rc.size() == W) { + const int64_t step = interval_start - m_grad_slide.last_interval_start; + const int64_t step_end = interval_end - m_grad_slide.last_interval_end; + if (step > 0 && step == step_end && (size_t)step <= W) { + stride = (size_t)step; + can_slide = true; + } + } + + if (can_slide) { + if (m_strand != -1) { + // Target is fwd-orientation: a 1bp genomic forward step drops the + // anchor at the LEFT of the target (low target index) and adds a + // new anchor at the RIGHT (high target index). The fetched window + // for the new pivot has its left edge advanced by `stride`, so + // target indices for the SAME genomic anchor have shifted by + // -stride; in target-local terms we pop_front and push_back. + for (size_t k = 0; k < stride; ++k) { + m_grad_slide.fwd.pop_front(); + m_grad_slide.rc.pop_front(); + const size_t new_anchor_i = i_max - stride + 1 + k; + float fwd_v = -std::numeric_limits::infinity(); + float rc_v = -std::numeric_limits::infinity(); + if (use_fwd) fwd_v = score_forward_original(m_pssm, target, new_anchor_i, m_strand); + if (use_rc) rc_v = score_reverse_original(m_pssm, target, new_anchor_i, m_strand); + m_grad_slide.fwd.push_back(fwd_v); + m_grad_slide.rc.push_back(rc_v); + } + } else { + // Target is rc'd: the rc'd buffer's HIGH target index corresponds + // to the LOW genomic coordinate. A genomic forward step still adds + // new anchors at the LOW genomic coordinates (because the iterator + // step extends the END of the genomic window), which in rc'd-target + // space means LOW target indices. Pop trailing tail anchors + // (high target indices) and push new anchors at the front to + // preserve target-index ordering. + for (size_t k = 0; k < stride; ++k) { + m_grad_slide.fwd.pop_back(); + m_grad_slide.rc.pop_back(); + } + for (size_t k = stride; k > 0; --k) { + const size_t new_anchor_i = i_min + (k - 1); + float fwd_v = -std::numeric_limits::infinity(); + float rc_v = -std::numeric_limits::infinity(); + if (use_fwd) fwd_v = score_forward_original(m_pssm, target, new_anchor_i, m_strand); + if (use_rc) rc_v = score_reverse_original(m_pssm, target, new_anchor_i, m_strand); + m_grad_slide.fwd.push_front(fwd_v); + m_grad_slide.rc.push_front(rc_v); + } + } + m_grad_slide.last_interval_start = interval_start; + m_grad_slide.last_interval_end = interval_end; + } else { + // Seed: full scan once. + m_grad_slide.fwd.clear(); + m_grad_slide.rc.clear(); + for (size_t k = 0; k < W; ++k) { + const size_t i = i_min + k; + float fwd_v = -std::numeric_limits::infinity(); + float rc_v = -std::numeric_limits::infinity(); + if (use_fwd) fwd_v = score_forward_original(m_pssm, target, i, m_strand); + if (use_rc) rc_v = score_reverse_original(m_pssm, target, i, m_strand); + m_grad_slide.fwd.push_back(fwd_v); + m_grad_slide.rc.push_back(rc_v); + } + m_grad_slide.valid = true; + m_grad_slide.chromid = chromid; + m_grad_slide.strand_mode = m_strand; + m_grad_slide.last_interval_start = interval_start; + m_grad_slide.last_interval_end = interval_end; + m_grad_slide.i_min = i_min; + m_grad_slide.i_max = i_max; + } + + // ---- Aggregation phase ---------------------------------------------------- + float lse = -std::numeric_limits::infinity(); float lse_no_head = -std::numeric_limits::infinity(); float best = -std::numeric_limits::infinity(); @@ -472,17 +566,10 @@ float PWMScorer::score_grad(const std::string& target, float head_spat_log = 0.0f; bool any_finite = false; - for (size_t i = i_min; i <= i_max; ++i) { - // score_forward_original / score_reverse_original return scores in - // fwd-genome semantics regardless of whether `target` is rc'd. - float fwd = -std::numeric_limits::infinity(); - float rc = -std::numeric_limits::infinity(); - if (use_fwd) { - fwd = score_forward_original(m_pssm, target, i, m_strand); - } - if (use_rc) { - rc = score_reverse_original(m_pssm, target, i, m_strand); - } + for (size_t k = 0; k < W; ++k) { + const size_t i = i_min + k; + const float fwd = m_grad_slide.fwd[k]; + const float rc = m_grad_slide.rc[k]; // Per-anchor raw score: comb (bidirect) or single-strand. float s_raw; @@ -501,8 +588,10 @@ float PWMScorer::score_grad(const std::string& target, s_raw = use_fwd ? fwd : rc; } - // Apply spatial log factor (= 0 when m_use_spat is false). - const float spat_log = get_spatial_log_factor(i - i_min); + // Apply spatial log factor (= 0 when m_use_spat is false). Spatial is + // iterator-local (indexed by k = i - i_min), so it's recomputed per + // pivot even though the underlying fwd/rc are cached. + const float spat_log = get_spatial_log_factor(k); const float s = std::isfinite(s_raw) ? (s_raw + spat_log) : s_raw; if (i == head_idx) { @@ -1180,7 +1269,9 @@ float PWMScorer::score_interval(const GInterval& interval, const GenomeChromKey& m_mode == GRAD_MAX_ISM || m_mode == GRAD_LSE_ISM) { const bool lse_aggregate = (m_mode == GRAD_LSE || m_mode == GRAD_LSE_ISM); const bool ism = (m_mode == GRAD_MAX_ISM || m_mode == GRAD_LSE_ISM); - return score_grad(target, i_min, i_max, motif_len, lse_aggregate, ism); + return score_grad(target, i_min, i_max, motif_len, + interval.chromid, interval.start, interval.end, + lse_aggregate, ism); } // Try spatial sliding window optimization diff --git a/src/PWMScorer.h b/src/PWMScorer.h index 79b063876..ce743d242 100644 --- a/src/PWMScorer.h +++ b/src/PWMScorer.h @@ -66,6 +66,38 @@ class PWMScorer : public GenomeSeqScorer int hit_count = 0; }; + // Per-anchor (fwd, rc) raw score cache for grad modes. Sliding by 1bp pops + // one anchor at the trailing edge and adds one at the leading edge; the + // per-pivot aggregate is recomputed from the cache (O(W)) but the per-anchor + // scoring is amortized to O(L) across slides, replacing an O(W*L) full scan. + // + // Cache holds RAW (no spatial weighting) per-strand scores so the per-pivot + // spatial assignment can be re-applied: different pivots assign different + // spatial bins to the same anchor (spatial is iterator-local, not anchor- + // genomic). + struct GradSlideCache { + bool valid = false; + int chromid = -1; + char strand_mode = 0; + int64_t last_interval_start = -1; + int64_t last_interval_end = -1; + size_t i_min = 0; + size_t i_max = 0; // i_max - i_min + 1 == window size + std::deque fwd; // raw fwd score per anchor in target-index order + std::deque rc; // raw rc score per anchor in target-index order + + void clear() { + valid = false; + chromid = -1; + last_interval_start = -1; + last_interval_end = -1; + i_min = 0; + i_max = 0; + fwd.clear(); + rc.clear(); + } + }; + // Spatial sliding window cache with bin-aware delta updates struct SpatSlideCache { // Geometry / semantics @@ -179,6 +211,7 @@ class PWMScorer : public GenomeSeqScorer float score_grad(const std::string& target, size_t i_min, size_t i_max, size_t motif_length, + int chromid, int64_t interval_start, int64_t interval_end, bool lse_aggregate, bool ism); // Core members @@ -200,6 +233,7 @@ class PWMScorer : public GenomeSeqScorer // Cache SlideCache m_slide; SpatSlideCache m_spat_slide; + GradSlideCache m_grad_slide; }; #endif // PWM_SCORER_H_ From 78ead3b789bd4325fe3e0106ba583d1104e12b3b Mon Sep 17 00:00:00 2001 From: aviezerl Date: Fri, 8 May 2026 11:47:26 +0300 Subject: [PATCH 13/15] docs(vignette): rework strands section, drop LSE-vs-MAX motivation - 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. --- vignettes/PWM-Functions.Rmd | 319 ++++++++++++++++-------------------- 1 file changed, 142 insertions(+), 177 deletions(-) diff --git a/vignettes/PWM-Functions.Rmd b/vignettes/PWM-Functions.Rmd index 781e073d7..c0530b6d3 100644 --- a/vignettes/PWM-Functions.Rmd +++ b/vignettes/PWM-Functions.Rmd @@ -15,19 +15,13 @@ knitr::opts_chunk$set( ) ``` -This vignette covers misha's PWM virtual track family: the existing -aggregating scorers (`pwm`, `pwm.max`, `pwm.max.pos`, `pwm.count`), the +This vignette covers misha's PWM virtual track family: the aggregating +scorers (`pwm`, `pwm.max`, `pwm.max.pos`, `pwm.count`), the edit-distance variants, and the per-position attribution tracks (`pwm.grad`, `pwm.grad.ism`). -The aggregating sections (`pwm`, `pwm.max`, ...) are well documented on -the `gvtrack.create` help page, so the focus here is on the mental -model and on examples a reader can adapt to their own database. The -math for the gradient tracks is worked out in more detail because they -are new and there is no existing reference. - The math chunks below run as plain R (no genome database required) so -they always evaluate. The genome-track examples are shown with +they always evaluate. Genome-track examples are shown with `eval = FALSE` because they need a database; the patterns transfer unchanged to any misha DB. @@ -77,15 +71,14 @@ raw count matrix often wants something larger. ## Per-anchor score Given a sequence `s`, an *anchor* is a position `a` at which the motif -is laid down. The per-anchor (single-strand, forward) score is +is laid down. The forward-strand per-anchor score is $$ -\text{score}_a \;=\; \sum_{c=0}^{L-1} M[c,\, s_{a+c}]. +\text{score}_a^{\text{fwd}} \;=\; \sum_{c=0}^{L-1} M[c,\, s_{a+c}]. $$ It is just a column-by-column lookup followed by a sum. Units are log -probability; values are <= 0 by construction (each term is the log of a -probability in [0, 1]). +probability; values are <= 0 by construction. ```{r} # Score the anchor at position 1 of the sequence "ACGTTT". @@ -99,9 +92,6 @@ score_anchor(M, bases, 1) # "ACG" - perfect consensus score_anchor(M, bases, 2) # "CGT" - off-consensus ``` -The first anchor matches the consensus and scores close to 0; the -second is well below. - For a window of length `W >= L`, we get `W - L + 1` per-anchor scores. ```{r} @@ -114,22 +104,71 @@ all_anchor_scores <- function(M, seq) { all_anchor_scores(M, seq) ``` +## Strands and `bidirect` + +Anchors are always indexed in **forward-genome coordinates**: anchor +`a` covers genomic positions `[a, a + L)`. What changes with the +strand setting is the *score* attached to each anchor. + +For each genomic anchor `a` there are two scores: + +- **Forward score** `score_a^{fwd}`: PSSM `M` applied directly to the + bases at `[a, a + L)`. Column `c` of the PSSM looks up base `s_{a+c}`. +- **Reverse score** `score_a^{rc}`: PSSM applied to the + reverse-complement of `[a, a + L)`. Equivalently, column `c` of the + PSSM looks up the *complement* of `s_{a + L - 1 - c}`. + +`bidirect` and `strand` decide which of these the per-anchor aggregate +uses: + +| Setting | Per-anchor score | +|---|---| +| `bidirect = TRUE` (default) | `log(exp(score_a^{fwd}) + exp(score_a^{rc}))` (per-genomic-position union over strands) | +| `bidirect = FALSE, strand = 1` | `score_a^{fwd}` | +| `bidirect = FALSE, strand = -1` | `score_a^{rc}` | + +The `bidirect = TRUE` case is the most common. The combined score is +the LSE union of the two strands per anchor; equivalent to "what's the +best motif support at this genomic position, considering both +orientations". This is computed at *every* anchor, then aggregated +(MAX / LSE / COUNT) over the iterator window. + +The `strand` argument is **only** used when `bidirect = FALSE`. + +### Position reporting + +Several tracks return positions: + +- **`pwm.max.pos`**: 1-based offset of the argmax anchor inside the + iterator interval (after `gvtrack.iterator()` shifts). With + `bidirect = TRUE`, the value is **signed**: positive when the + forward strand wins at the argmax, negative when the reverse strand + wins. The absolute value is always the offset, never a genomic + coordinate. +- **`pwm.grad`, `pwm.grad.ism`**: report the gradient at a single + position - the iterator interval **start**, in forward-genome + coordinates, regardless of strand. + +So if you want a per-bp gradient track at every position `p` in a +region, you call `gextract` with `iterator = 1` and an iterator +extension `gvtrack.iterator(eshift = W - 1)` so each per-bp evaluation +sees a window `[p, p + W)`. The reported value is at `p`, in +fwd-genome coords, exactly the position you iterated to. + ## Aggregating over an interval -The PWM virtual tracks each reduce these per-anchor scores to a single -scalar per iterator interval. +The aggregating PWM virtual tracks reduce the per-anchor scores in an +iterator window to a single scalar. ### `pwm.max`: best anchor -`pwm.max` returns the maximum per-anchor score in the window. Easiest -to interpret ("how good is the best motif hit here?"), but discards -information about secondary hits. +`pwm.max` returns the maximum per-anchor score in the window. ```{r} max(all_anchor_scores(M, seq)) ``` -### `pwm` (LSE): smooth aggregate +### `pwm` (LSE): integrated likelihood `pwm` returns the log-sum-exp of all per-anchor scores: @@ -138,10 +177,9 @@ f_{\text{LSE}} \;=\; \log \sum_a \exp(\text{score}_a). $$ This is dominated by the best anchor when one anchor clearly wins, but -also adds support for many mediocre anchors (an "OR" over noisy -matches). LSE is also smooth: small changes in the input give small -changes in the output, which is friendlier for downstream filtering and -gradient-based attribution. +also picks up support from many mediocre anchors (an "OR" over noisy +matches). Useful when you care about cumulative motif evidence in a +region rather than just the single best hit. ```{r} log_sum_exp <- function(x) { @@ -153,20 +191,19 @@ log_sum_exp(s) ``` For a single dominant hit `pwm` and `pwm.max` are nearly equal; for a -window with many weak similar-strength hits `pwm` can be substantially +window with many similar-strength hits `pwm` can be substantially larger. ### `pwm.max.pos`: argmax position -`pwm.max.pos` returns the 1-based offset of the argmax anchor inside -the iterator interval, signed by strand when `bidirect = TRUE`. Useful -when you want to know *where* the best hit sits within a region. +Returns the offset of the argmax anchor inside the iterator interval. +See "Position reporting" above for the sign convention. ### `pwm.count`: thresholded hit count `pwm.count` requires `score.thresh`. It counts anchors whose score is -at or above `score.thresh`. With `bidirect = TRUE` the union over -strands is taken at each genomic start, so palindromic hits count once. +at or above `score.thresh`. With `bidirect = TRUE` the per-anchor +score is the strand union, so palindromic hits count once. ```{r, eval = FALSE} gvtrack.create("motif_count", NULL, "pwm.count", @@ -176,30 +213,11 @@ gvtrack.iterator("motif_count", sshift = 0, eshift = 500) gextract("motif_count", iterator = 500, intervals = my_intervals) ``` -### Why LSE rather than just MAX? - -For any aggregation over a long window, LSE has three practical -advantages over MAX: - -1. **Integrates many anchors.** A region with five mediocre hits has a - higher LSE than one with a single mediocre hit, but the same MAX. -2. **Smooth.** Differentiable everywhere; small input changes -> small - output changes. This matters for gradient attribution - (next section). -3. **Numerically stable when implemented with the standard - `m + log(sum(exp(x - m)))` trick.** - -When you only care about "is there a strong hit anywhere in this -window?", `pwm.max` is fine. For typical "motif energy" scoring across -500 bp - 1 kb windows, `pwm` (LSE) is the better default. - ## Common parameters These mirror what `?gvtrack.create` documents; one-liner each: -- **`bidirect`** (default `TRUE`): scan both strands and combine per - genomic start. With `bidirect = FALSE`, only the strand given by - `strand` is scanned. +- **`bidirect`** (default `TRUE`): see "Strands and `bidirect`" above. - **`strand`** (`+1` / `-1`): only used when `bidirect = FALSE`. - **`prior`** (default `0.01`): Dirichlet pseudocount added to PSSM rows before normalization. Raise it for noisy / count-based PSSMs; @@ -244,59 +262,31 @@ See `?gvtrack.create` for the full parameter list and worked examples. ## Per-position attribution: `pwm.grad` and `pwm.grad.ism` -The aggregating tracks above all collapse a window into one scalar. -`pwm.grad` and `pwm.grad.ism` answer a different question: **at this -specific base pair, how much does the local sequence support the PWM -aggregate?** The output is one non-negative number per iterator -interval, computed at the anchor that starts at `interval.start`. - -Combined with a `iterator = 1` step and `gvtrack.iterator(eshift = -W - 1)`, this gives a per-bp saliency / attribution track you can plot -or combine with conservation, ATAC, methylation, etc. +The aggregating tracks above collapse a window into one scalar. +`pwm.grad` and `pwm.grad.ism` return one scalar per iterator interval +**evaluated at a specific bp** - the anchor that starts at +`interval.start`. With `iterator = 1` and `gvtrack.iterator(eshift = +W - 1)`, this gives a per-bp track over a region. -### Motivation +### Linearized gradient (`pwm.grad`) -Two use cases: +Treat each base as a one-hot vector `x_p` over A/C/G/T. With pivot +fixed at `p = interval.start`, only one anchor (`a = p`) contains `p`, +and only column 0 of the PSSM is involved. -1. **Saliency.** "Of the base pairs in this window, which ones are - actually carrying the motif signal?" - a poor man's - integrated-gradient / DeepLIFT. -2. **Single-bp coupling.** Combine PWM evidence with another per-bp - signal (conservation, ChIP, etc.) at single-base resolution. With - `pwm` you only get one number per window; with `pwm.grad` you get - one per bp. - -### Math intuition (linearized, `pwm.grad`) - -Think of each base as a one-hot vector `x_p` over A/C/G/T. The score -`score_a` is linear in `x_p` whenever `a` contains `p`, with gradient -`M[p - a, b]`. With pivot fixed at the iterator start `p = -interval.start`, only one anchor (`a = p`) contains `p`, and only -column 0 of the PSSM is involved. - -For the LSE aggregate, the gradient is the softmax-weighted column-0 -score, contracted against `(x_p - x_p^{\text{worst}})` where the -reference is the worst possible base at column 0: +For the LSE aggregate: $$ g_{\text{LSE}}(p) \;=\; w_p \cdot \big(M[0,\, b_p] - \min_b M[0,\, b]\big), \qquad w_p = \exp(\text{score}_p - f_{\text{LSE}}). $$ -Three useful sanity points: +`w_p` is the softmax weight of the head anchor in the window: it +captures how much the head anchor contributes to the LSE relative to +the other `W - L` anchors. The diff term equals the column-0 +information content under an "actual minus worst" baseline. -- **Non-negative.** `M[0, b_p] - min_b M[0, b]` >= 0 by construction; - `w_p` in [0, 1]. -- **Bounded by per-column information.** The diff term equals the - column's information content under the "actual minus worst" baseline. -- **Modulated by softmax weight.** If the head anchor dominates the - LSE (`w_p ≈ 1`), the gradient saturates to the diff. If the LSE - spreads over many anchors (`w_p` small), the gradient is small even - if the head anchor is itself a good match - the *attribution* of - this particular base to the aggregate is diluted. - -For the MAX aggregate, the gradient picks up the diff iff the head -anchor is the argmax: +For the MAX aggregate: $$ g_{\text{MAX}}(p) \;=\; @@ -306,40 +296,51 @@ M[0,\, b_p] - \min_b M[0,\, b] & \text{if } a^* = p,\\ \end{cases} $$ -### Math intuition (in-silico mutagenesis, `pwm.grad.ism`) +i.e. the column-0 diff iff the head anchor is the argmax of the window. -The ISM gradient is "actual minus best replacement": +### In-silico mutagenesis (`pwm.grad.ism`) $$ g_{\text{ISM}}(p) \;=\; f(\text{seq}) \;-\; \min_{b' \neq b_p} f\big(\text{seq with } s_p := b'\big). $$ -Always >= 0. With pivot fixed at `interval.start`, only `score_p` -changes when we flip `s_p`, so the update is O(1) per flip and ISM -costs only ~3-4x more than the linearized form. +Flip the head base to each of the three other letters, recompute the +aggregate, take the worst, and subtract from the unperturbed +aggregate. With pivot at `interval.start`, only `score_p` changes +under the flip, so the update is O(1) per alternative. -### When `pwm.grad` and `pwm.grad.ism` differ +### When the two definitions differ -The two definitions agree exactly when the argmax is stable under -flips of `s_p`, and agree to high precision when LSE concentrates on -the head anchor. They diverge in two cases: +The two agree exactly when the argmax is stable under flips of `s_p`, +and agree to high precision when LSE concentrates on the head anchor. +They diverge when: -1. **Argmax shifts.** For MAX, if flipping `s_p` makes a *different* - anchor become the argmax, ISM accounts for that ("the second-best - hit kicks in") and the linearized form doesn't. -2. **Curvature in LSE.** Linearization is a first-order approximation - around the actual sequence; ISM evaluates the aggregate at three - discrete neighbours and takes the worst. For sharply non-linear - regions of the LSE landscape they differ slightly, but in practice - the disagreement is small. +1. **Argmax shifts under flip.** For MAX, flipping `s_p` may make a + *different* anchor become the argmax; ISM picks that up and the + linearized form doesn't. +2. **LSE curvature.** Linearization is first-order around the actual + sequence; ISM evaluates the aggregate at three discrete neighbours. + The disagreement is small in practice. -In short: **`pwm.grad` is cheaper and almost always good enough; -`pwm.grad.ism` is more conservative around argmax-shift edge cases.** +`pwm.grad` is cheaper; `pwm.grad.ism` is more conservative around +argmax-shift cases. -### Worked numerical example +### Strands -Consider an `L = 3` PSSM and a 6-bp window starting with consensus -"ACG": +The pivot is at `interval.start` in fwd-genome coordinates regardless +of strand. What changes is which PSSM column matters: + +- `bidirect = FALSE, strand = 1`: column 0 of `M`, base `s_p`. +- `bidirect = FALSE, strand = -1`: column `L-1` of `M`, complement of + `s_p` (because the rc-strand PSSM applied at `p` reads its + *first* column from `s_p` complemented, which corresponds to fwd + PSSM column `L-1`). +- `bidirect = TRUE`: per-strand softmax at the head anchor combines + the two contributions. You don't need to think about it explicitly; + pass `bidirect = TRUE` and the engine computes the appropriate + combined gradient. + +### Worked numerical example ```{r} pssm <- matrix( @@ -354,7 +355,7 @@ pssm <- matrix( M <- log((pssm + 0.01) / rowSums(pssm + 0.01)) seq <- "ACGTTT" -s <- all_anchor_scores(M, seq) # 4 per-anchor scores +s <- all_anchor_scores(M, seq) s ``` @@ -363,31 +364,21 @@ Linearized LSE gradient at `p = 1`: ```{r} f_lse <- log_sum_exp(s) score_p <- s[1] -w_p <- exp(score_p - f_lse) # softmax weight at the head -b_p <- "A" # head base +w_p <- exp(score_p - f_lse) +b_p <- "A" diff_col0 <- M[1, b_p] - min(M[1, ]) -w_p -diff_col0 g_lse <- w_p * diff_col0 -g_lse +c(w_p = w_p, diff_col0 = diff_col0, g_lse = g_lse) ``` -`w_p ≈ 1` because the head anchor "ACG" is the dominant LSE -contributor; `diff_col0` is the column-0 information content. The -result is essentially the column's entropy gap. - ISM gradient: flip `seq[1]` to each of C/G/T and recompute the LSE. ```{r} -ism_lse <- function(M, seq) { - s <- all_anchor_scores(M, seq) - log_sum_exp(s) -} +ism_lse <- function(M, seq) log_sum_exp(all_anchor_scores(M, seq)) flips <- setdiff(c("A", "C", "G", "T"), substr(seq, 1, 1)) f_alt <- vapply(flips, function(b) { ism_lse(M, paste0(b, substr(seq, 2, nchar(seq)))) }, numeric(1)) -f_alt ism_lse(M, seq) - min(f_alt) ``` @@ -396,11 +387,8 @@ case. ### Genome-track recipe -Here is the canonical pattern: extract `pwm` and `pwm.grad` over a 1 -kb iterator window, then walk the per-bp gradient with `iterator = 1`. - ```{r, eval = FALSE} -W <- 500 # window length (bp) used for both aggregate and gradient +W <- 500 # window length pssm <- matrix(...) # your 4-col PSSM, dimnames A/C/G/T @@ -408,7 +396,6 @@ gvtrack.create("score", NULL, "pwm", pssm = pssm) gvtrack.create("grad", NULL, "pwm.grad", pssm = pssm, aggregate = "lse") gvtrack.create("grad_ism", NULL, "pwm.grad.ism", pssm = pssm, aggregate = "lse") -# Each per-bp value uses a window starting at this bp and extending W bp right. gvtrack.iterator("score", sshift = 0, eshift = W - 1) gvtrack.iterator("grad", sshift = 0, eshift = W - 1) gvtrack.iterator("grad_ism", sshift = 0, eshift = W - 1) @@ -420,46 +407,24 @@ out <- gextract(c("score", "grad", "grad_ism"), head(out) ``` -Common downstream patterns: - -- **Per-bp coupling**: `out$grad * out$other_track` to get an - attribution-weighted version of `other_track`. -- **Argmax-shift hunting**: `out$grad_ism - out$grad` highlights bp - where ISM and linearization diverge, i.e. where the argmax flips - under a single substitution. -- **Smoothing**: a `gtrack.smooth` over `pwm.grad` gives a soft - saliency band over the canonical motif location. - -### `bidirect = TRUE` for asymmetric PSSMs - -For palindromic PSSMs the rc strand contributes column-0 of the rc -PSSM (= column `L-1` of the fwd PSSM, with complement). With -`bidirect = TRUE` the gradient uses a per-anchor strand softmax to -combine the two contributions; this is handled internally and you -don't need to touch anything beyond `bidirect = TRUE`. For asymmetric -PSSMs, the bidirect gradient at a bp can pick up signal from either -strand depending on which one wins locally. +`out$grad - out$grad_ism` flags positions where the argmax flips under +single substitution. Multiplying `out$grad` by another per-bp track +gives an attribution-weighted version of that track. ## Performance notes -In rough numbers for an `L = 12` PSSM, `W = 500` bp window, scanned at -1 bp resolution: - -- `pwm.grad` is roughly comparable to `pwm` (both do the same scan; - the gradient post-step is O(1) per anchor). -- `pwm.grad.ism` is ~3-4x `pwm.grad` (three flips x O(1) per - position, with a small constant for log-space arithmetic). - -For genome-wide attribution, you almost never want `iterator = 1` -across the whole genome. Common patterns: - -- Use a sparse intervals set as iterator (e.g. peaks, TSS windows) and - `iterator = 1` only inside those; gives single-bp attribution where - it matters and skips empty space. -- Use a coarser iterator step (e.g. `iterator = 10` or `iterator = - 100`) for whole-genome scans; gradient is still meaningful, just - averaged over the step. - -When in doubt, profile on one chromosome before committing to a -genome-wide scan: the expensive part is reading sequence from disk, -not the gradient computation itself. +For an `L = 12` PSSM, `W = 500` bp window, scanned at 1 bp resolution: + +- `pwm` and `pwm.max` use a sliding-window cache and are very fast + (O(L) per pivot amortized). +- `pwm.grad` and `pwm.grad.ism` use a per-anchor sliding cache (the + per-anchor scoring is amortized to O(L) across pivots) but + recompute the per-window aggregate per pivot (O(W)). Expect + roughly 10-15x slower than `pwm` at iterator = 1. +- `pwm.grad.ism` is essentially the same cost as `pwm.grad`: the + three flips contribute O(1) each on top of the per-pivot scan. + +For genome-wide attribution, a sparse iterator (peaks, TSS windows) +plus `iterator = 1` *inside* those is usually the right shape. Coarser +iterator steps (`iterator = 10` or `100`) work for whole-genome scans +when single-bp resolution isn't required. From 0ccd5b09c80d1629fa03f5622481870d20ccc078 Mon Sep 17 00:00:00 2001 From: aviezerl Date: Fri, 8 May 2026 13:12:25 +0300 Subject: [PATCH 14/15] docs(pwm.grad): tone down DeepLIFT / saliency framing MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- NEWS.md | 2 +- R/vtrack.R | 2 +- man/gvtrack.create.Rd | 2 +- src/PWMScorer.h | 8 ++++---- vignettes/PWM-Functions.Rmd | 29 ++++++++++++++++++++++------- 5 files changed, 29 insertions(+), 14 deletions(-) diff --git a/NEWS.md b/NEWS.md index 706c87f12..9a3f6b744 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,6 @@ # misha (development version) -* Added `pwm.grad` and `pwm.grad.ism` virtual tracks: per-bp gradient / saliency of the PWM aggregate (LSE or MAX) at the iterator interval start. `pwm.grad` is the linearized DeepLIFT-style attribution; `pwm.grad.ism` is in-silico mutagenesis. +* Added `pwm.grad` and `pwm.grad.ism` virtual tracks: per-bp PSSM-column contribution at the iterator interval start, derived from the LSE or MAX aggregate. `pwm.grad` is the softmax-weighted column contribution (linearized form); `pwm.grad.ism` is in-silico mutagenesis (`motifbreakR`-style flip-and-rescore, packaged as a sliding-window track). # misha 5.6.28 diff --git a/R/vtrack.R b/R/vtrack.R index c55a3d6f1..65e6fca24 100644 --- a/R/vtrack.R +++ b/R/vtrack.R @@ -705,7 +705,7 @@ #' NULL (sequence) \tab pwm.max \tab pssm, bidirect, prior, extend, spat_* \tab Maximum log-likelihood score among all anchors (per-position union across strands). \cr #' NULL (sequence) \tab pwm.max.pos \tab pssm, bidirect, prior, extend, spat_* \tab 1-based position of the best-scoring anchor (signed by strand when \code{bidirect = TRUE}); coordinates are always relative to the iterator interval after any \code{gvtrack.iterator()} shifts/extensions. \cr #' NULL (sequence) \tab pwm.count \tab pssm, score.thresh, bidirect, prior, extend, strand, spat_* \tab Count of anchors whose score exceeds \code{score.thresh} (per-position union). \cr -#' NULL (sequence) \tab pwm.grad \tab pssm, aggregate, bidirect, prior, extend, spat_* \tab Linearized per-bp gradient (DeepLIFT-style) of the PWM aggregate at the interval start. \cr +#' NULL (sequence) \tab pwm.grad \tab pssm, aggregate, bidirect, prior, extend, spat_* \tab Softmax-weighted per-bp PSSM-column contribution at the interval start (linearized form). \cr #' NULL (sequence) \tab pwm.grad.ism \tab pssm, aggregate, bidirect, prior, extend, spat_* \tab In-silico mutagenesis per-bp gradient at the interval start. \cr #' } #' diff --git a/man/gvtrack.create.Rd b/man/gvtrack.create.Rd index 02a35d6d8..14aba2d73 100644 --- a/man/gvtrack.create.Rd +++ b/man/gvtrack.create.Rd @@ -138,7 +138,7 @@ interval after all modifier adjustments. NULL (sequence) \tab pwm.max \tab pssm, bidirect, prior, extend, spat_* \tab Maximum log-likelihood score among all anchors (per-position union across strands). \cr NULL (sequence) \tab pwm.max.pos \tab pssm, bidirect, prior, extend, spat_* \tab 1-based position of the best-scoring anchor (signed by strand when \code{bidirect = TRUE}); coordinates are always relative to the iterator interval after any \code{gvtrack.iterator()} shifts/extensions. \cr NULL (sequence) \tab pwm.count \tab pssm, score.thresh, bidirect, prior, extend, strand, spat_* \tab Count of anchors whose score exceeds \code{score.thresh} (per-position union). \cr - NULL (sequence) \tab pwm.grad \tab pssm, aggregate, bidirect, prior, extend, spat_* \tab Linearized per-bp gradient (DeepLIFT-style) of the PWM aggregate at the interval start. \cr + NULL (sequence) \tab pwm.grad \tab pssm, aggregate, bidirect, prior, extend, spat_* \tab Softmax-weighted per-bp PSSM-column contribution at the interval start (linearized form). \cr NULL (sequence) \tab pwm.grad.ism \tab pssm, aggregate, bidirect, prior, extend, spat_* \tab In-silico mutagenesis per-bp gradient at the interval start. \cr } diff --git a/src/PWMScorer.h b/src/PWMScorer.h index ce743d242..7479e15d0 100644 --- a/src/PWMScorer.h +++ b/src/PWMScorer.h @@ -19,10 +19,10 @@ class PWMScorer : public GenomeSeqScorer MAX_LIKELIHOOD, // Maximum log-likelihood score MAX_LIKELIHOOD_POS, // Position of maximum log-likelihood MOTIF_COUNT, // Count of positions exceeding threshold - GRAD_LSE, // Linearized per-bp gradient (DeepLIFT-style) under LSE aggregate - GRAD_MAX, // Linearized per-bp gradient under MAX aggregate - GRAD_LSE_ISM, // In-silico mutagenesis per-bp gradient under LSE aggregate - GRAD_MAX_ISM // In-silico mutagenesis per-bp gradient under MAX aggregate + GRAD_LSE, // Softmax-weighted per-bp PSSM-column contribution under LSE aggregate + GRAD_MAX, // Argmax-conditioned per-bp PSSM-column contribution under MAX aggregate + GRAD_LSE_ISM, // In-silico mutagenesis per-bp contribution under LSE aggregate + GRAD_MAX_ISM // In-silico mutagenesis per-bp contribution under MAX aggregate }; PWMScorer(const DnaPSSM &pssm, const std::string &genome_root, bool extend = true, diff --git a/vignettes/PWM-Functions.Rmd b/vignettes/PWM-Functions.Rmd index c0530b6d3..68aa0e5ea 100644 --- a/vignettes/PWM-Functions.Rmd +++ b/vignettes/PWM-Functions.Rmd @@ -17,7 +17,7 @@ knitr::opts_chunk$set( This vignette covers misha's PWM virtual track family: the aggregating scorers (`pwm`, `pwm.max`, `pwm.max.pos`, `pwm.count`), the -edit-distance variants, and the per-position attribution tracks +edit-distance variants, and the per-position contribution tracks (`pwm.grad`, `pwm.grad.ism`). The math chunks below run as plain R (no genome database required) so @@ -260,7 +260,7 @@ than per-window scores. See `?gvtrack.create` for the full parameter list and worked examples. -## Per-position attribution: `pwm.grad` and `pwm.grad.ism` +## Per-position contribution: `pwm.grad` and `pwm.grad.ism` The aggregating tracks above collapse a window into one scalar. `pwm.grad` and `pwm.grad.ism` return one scalar per iterator interval @@ -268,7 +268,22 @@ The aggregating tracks above collapse a window into one scalar. `interval.start`. With `iterator = 1` and `gvtrack.iterator(eshift = W - 1)`, this gives a per-bp track over a region. -### Linearized gradient (`pwm.grad`) +A PWM is a linear log-likelihood model, so "the gradient" is just a +PSSM column. These tracks expose two flavors of per-base contribution: + +- **`pwm.grad`**: a softmax-weighted column-0 contribution. The math + is the gradient of the LSE / MAX aggregate w.r.t. the one-hot input + at the head bp, contracted with `actual_base - worst_base`. This + reduces to "how much does this base's PSSM column contribute, + weighted by how much the head anchor matters in the window". +- **`pwm.grad.ism`**: in-silico mutagenesis. Same idea as + `motifbreakR`-style SNP-impact scoring, exposed as a sliding-window + track: flip the base, recompute the aggregate, take the worst. + +These are not new methods - they are standard motif-scanning math +packaged as misha vtracks composable with other 1D tracks. + +### Softmax-weighted contribution (`pwm.grad`) Treat each base as a one-hot vector `x_p` over A/C/G/T. With pivot fixed at `p = interval.start`, only one anchor (`a = p`) contains `p`, @@ -283,8 +298,8 @@ $$ `w_p` is the softmax weight of the head anchor in the window: it captures how much the head anchor contributes to the LSE relative to -the other `W - L` anchors. The diff term equals the column-0 -information content under an "actual minus worst" baseline. +the other anchors. The diff term equals the column-0 information +content under an "actual minus worst" baseline. For the MAX aggregate: @@ -409,7 +424,7 @@ head(out) `out$grad - out$grad_ism` flags positions where the argmax flips under single substitution. Multiplying `out$grad` by another per-bp track -gives an attribution-weighted version of that track. +gives a contribution-weighted version of that track. ## Performance notes @@ -424,7 +439,7 @@ For an `L = 12` PSSM, `W = 500` bp window, scanned at 1 bp resolution: - `pwm.grad.ism` is essentially the same cost as `pwm.grad`: the three flips contribute O(1) each on top of the per-pivot scan. -For genome-wide attribution, a sparse iterator (peaks, TSS windows) +For genome-wide per-bp scans, a sparse iterator (peaks, TSS windows) plus `iterator = 1` *inside* those is usually the right shape. Coarser iterator steps (`iterator = 10` or `100`) work for whole-genome scans when single-bp resolution isn't required. From 9085d2a995fd30669c87133797e929f4810d0acc Mon Sep 17 00:00:00 2001 From: aviezerl Date: Mon, 11 May 2026 13:20:27 +0300 Subject: [PATCH 15/15] fix(pwm.grad): respect gvtrack.filter at the pivot 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). --- src/SequenceVarProcessor.cpp | 19 ++++ tests/testthat/test-vtrack-pwm-grad.R | 138 ++++++++++++++++++++++++++ 2 files changed, 157 insertions(+) diff --git a/src/SequenceVarProcessor.cpp b/src/SequenceVarProcessor.cpp index d3aeed14e..18f2e9faa 100644 --- a/src/SequenceVarProcessor.cpp +++ b/src/SequenceVarProcessor.cpp @@ -402,9 +402,28 @@ void SequenceVarProcessor::process_individual_sequence_vars( vector unmasked_parts; ivar->filter->subtract(seq_interval, unmasked_parts); + const bool is_grad = + ivar->val_func == TrackExpressionVars::Track_var::PWM_GRAD || + ivar->val_func == TrackExpressionVars::Track_var::PWM_GRAD_ISM; + if (unmasked_parts.empty()) { // Completely masked ivar->var[idx] = numeric_limits::quiet_NaN(); + } else if (is_grad) { + // Gradient modes are anchored at seq_interval.start; we can only + // score the unmasked part that begins at that pivot. Any other + // part is at a different genomic anchor and must not contribute. + bool scored = false; + for (const auto& part : unmasked_parts) { + if (part.start == seq_interval.start) { + ivar->var[idx] = ivar->pwm_scorer->score_interval(part, m_iu.get_chromkey()); + scored = true; + break; + } + } + if (!scored) { + ivar->var[idx] = numeric_limits::quiet_NaN(); + } } else if (unmasked_parts.size() == 1) { // Single unmasked part ivar->var[idx] = ivar->pwm_scorer->score_interval(unmasked_parts[0], m_iu.get_chromkey()); diff --git a/tests/testthat/test-vtrack-pwm-grad.R b/tests/testthat/test-vtrack-pwm-grad.R index 1679e9e1f..328fa2885 100644 --- a/tests/testthat/test-vtrack-pwm-grad.R +++ b/tests/testthat/test-vtrack-pwm-grad.R @@ -965,3 +965,141 @@ test_that("pwm.grad: gsummary aggregation runs without error", { expect_true(is.numeric(s)) expect_true(length(s) >= 5L) }) + +# ============================================================ +# Filters +# ============================================================ + +# Helpers shared by the filter tests below. They all use a 1-bp pivot at +# chr1:200 and the same head-consensus PSSM, so the gradient is anchored at +# genomic position 200 and is positive when the pivot is unmasked. +.grad_filter_setup <- function(aggregate = "lse", L = 3L, iter_W = 20L, + pivot_start = 200L) { + pivot <- gintervals(1, pivot_start, pivot_start + 1L) + seq_ext <- toupper( + gseq.extract(gintervals(1, pivot_start, pivot_start + iter_W + L - 1L)) + ) + head_bases <- strsplit(substr(seq_ext, 1, L), "")[[1]] + pssm <- .consensus_pssm(head_bases) + list( + pivot = pivot, seq_ext = seq_ext, pssm = pssm, L = L, + iter_W = iter_W, pivot_start = pivot_start + ) +} + +test_that("pwm.grad: gvtrack.filter masking the pivot returns NA", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + s <- .grad_filter_setup(aggregate = "lse") + + gvtrack.create("g", NULL, "pwm.grad", + pssm = s$pssm, aggregate = "lse", + bidirect = FALSE, strand = 1, extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g", sshift = 0, eshift = s$iter_W - 1L) + # Mask covers the pivot at chr1:200. + gvtrack.filter("g", filter = gintervals(1, s$pivot_start, s$pivot_start + 5L)) + + res <- gextract("g", iterator = 1, intervals = s$pivot) + expect_true(is.na(res$g[1])) +}) + +test_that("pwm.grad: gvtrack.filter masking the pivot returns NA (ism)", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + s <- .grad_filter_setup(aggregate = "lse") + + gvtrack.create("g", NULL, "pwm.grad.ism", + pssm = s$pssm, aggregate = "lse", + bidirect = FALSE, strand = 1, extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g", sshift = 0, eshift = s$iter_W - 1L) + gvtrack.filter("g", filter = gintervals(1, s$pivot_start, s$pivot_start + 5L)) + + res <- gextract("g", iterator = 1, intervals = s$pivot) + expect_true(is.na(res$g[1])) +}) + +test_that("pwm.grad: multi-part filter scores only the pivot-anchored part", { + # Mask splits seq_interval into two unmasked parts. The pivot is in the + # first part; the second part lives at a different genomic anchor. The + # gradient must be computed at the pivot-anchored part only - summing + # across parts (the pre-fix behavior) would mix gradients at different + # anchor positions. + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + s <- .grad_filter_setup(aggregate = "lse") + + gvtrack.create("g", NULL, "pwm.grad", + pssm = s$pssm, aggregate = "lse", + bidirect = FALSE, strand = 1, extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g", sshift = 0, eshift = s$iter_W - 1L) + + # seq_interval = [200, 200 + iter_W + L - 1) = [200, 222) for default args. + # Mask [210, 215) -> unmasked parts [200, 210) and [215, 222). + mask_lo <- s$pivot_start + 10L + mask_hi <- s$pivot_start + 15L + gvtrack.filter("g", filter = gintervals(1, mask_lo, mask_hi)) + + res <- gextract("g", iterator = 1, intervals = s$pivot) + + # Expected: score on the first unmasked part [pivot, mask_lo). The engine + # expands that part by L-1 at the end for motif scanning, so the oracle + # sees (mask_lo - pivot) + (L - 1) bases. + part_len <- mask_lo - s$pivot_start + (s$L - 1L) + part_seq <- substr(s$seq_ext, 1, part_len) + expected <- manual_pwm_grad(s$pssm, part_seq, "lse", + bidirect = FALSE, strand = 1L, prior = 0.01 + ) + expect_equal(res$g[1], expected, tolerance = 1e-5, ignore_attr = TRUE) +}) + +test_that("pwm.grad.ism: multi-part filter scores only the pivot-anchored part", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + s <- .grad_filter_setup(aggregate = "lse") + + gvtrack.create("g", NULL, "pwm.grad.ism", + pssm = s$pssm, aggregate = "lse", + bidirect = FALSE, strand = 1, extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g", sshift = 0, eshift = s$iter_W - 1L) + + mask_lo <- s$pivot_start + 10L + mask_hi <- s$pivot_start + 15L + gvtrack.filter("g", filter = gintervals(1, mask_lo, mask_hi)) + + res <- gextract("g", iterator = 1, intervals = s$pivot) + part_len <- mask_lo - s$pivot_start + (s$L - 1L) + part_seq <- substr(s$seq_ext, 1, part_len) + expected <- manual_pwm_grad_ism(s$pssm, part_seq, "lse", + bidirect = FALSE, strand = 1L, prior = 0.01 + ) + expect_equal(res$g[1], expected, tolerance = 1e-5, ignore_attr = TRUE) +}) + +test_that("pwm.grad: gvtrack.filter not intersecting seq_interval matches unfiltered", { + remove_all_vtracks() + withr::defer(remove_all_vtracks()) + s <- .grad_filter_setup(aggregate = "lse") + + gvtrack.create("g_filt", NULL, "pwm.grad", + pssm = s$pssm, aggregate = "lse", + bidirect = FALSE, strand = 1, extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_filt", sshift = 0, eshift = s$iter_W - 1L) + gvtrack.create("g_nofilt", NULL, "pwm.grad", + pssm = s$pssm, aggregate = "lse", + bidirect = FALSE, strand = 1, extend = TRUE, prior = 0.01 + ) + gvtrack.iterator("g_nofilt", sshift = 0, eshift = s$iter_W - 1L) + + # Mask far from seq_interval; should not affect the result. + gvtrack.filter("g_filt", + filter = gintervals(1, s$pivot_start + 1000L, s$pivot_start + 1100L) + ) + + res <- gextract(c("g_filt", "g_nofilt"), iterator = 1, intervals = s$pivot) + expect_equal(res$g_filt[1], res$g_nofilt[1], tolerance = 1e-9) +})