diff --git a/NEWS.md b/NEWS.md index 036dfbcbf..9a3f6b744 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 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 * `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. diff --git a/R/vtrack.R b/R/vtrack.R index 36d628555..65e6fca24 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 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 #' } #' #' \strong{Edit distance summarizers} 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/man/gvtrack.create.Rd b/man/gvtrack.create.Rd index 9b8404f0a..14aba2d73 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 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 } \strong{Edit distance summarizers} diff --git a/src/PWMScorer.cpp b/src/PWMScorer.cpp index 16efb0a19..a0dbc24be 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() @@ -147,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 @@ -337,6 +367,411 @@ float PWMScorer::score_without_spatial(const std::string& target, int64_t motif_ return compute_position_result(pos_idx, target.length(), motif_length, best_dir); } +// Per-bp gradient at the iterator interval start. +// +// 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. +// +// 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. +// +// Strand handling unifies bidirect / fwd-only / rc-only via score_*_original, +// which canonicalize the rc'd target back to fwd-genome semantics. +// +// 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, + 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 + // re-checked here. Caller is score_interval, which handles iterator + // clamping and parameter validation upstream. + + 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(); + } + + 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; + } + + // 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]]; + } + 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 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. 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(); + float best_no_head = -std::numeric_limits::infinity(); + size_t argmax = i_min; + 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 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; + if (bidirect) { + if (!std::isfinite(fwd) && !std::isfinite(rc)) { + s_raw = -std::numeric_limits::infinity(); + } else if (!std::isfinite(fwd)) { + s_raw = rc; + } else if (!std::isfinite(rc)) { + s_raw = fwd; + } else { + s_raw = fwd; + log_sum_log(s_raw, rc); + } + } else { + s_raw = use_fwd ? fwd : rc; + } + + // 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) { + head_fwd = fwd; + head_rc = rc; + head_score = s; + head_spat_log = spat_log; + } + + if (std::isfinite(s)) { + any_finite = true; + 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; + } + } + } + } + + // 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; + } + + if (!any_finite) { + return std::numeric_limits::quiet_NaN(); + } + + // 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 (!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 (!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; + } + + // ---- 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(); + } + + // 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 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_raw = -std::numeric_limits::infinity(); + } else if (!std::isfinite(fwd_alt)) { + s_alt_raw = rc_alt; + } else if (!std::isfinite(rc_alt)) { + s_alt_raw = fwd_alt; + } else { + s_alt_raw = fwd_alt; + log_sum_log(s_alt_raw, rc_alt); + } + } else { + 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; + 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 float PWMScorer::score_with_spatial(const std::string& target, int64_t motif_length) { @@ -741,6 +1176,30 @@ float PWMScorer::score_with_sliding_window(const std::string& target, float PWMScorer::score_interval(const GInterval& interval, const GenomeChromKey& chromkey) { + // 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_MAX: + case GRAD_LSE: + 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) { + rdb::verror("pwm.grad: strand must be -1, 0, or 1"); + } + break; + 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); @@ -803,6 +1262,18 @@ float PWMScorer::score_interval(const GInterval& interval, const GenomeChromKey& i_min = 0; } + // 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_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, + interval.chromid, interval.start, interval.end, + lse_aggregate, ism); + } + // Try spatial sliding window optimization if (m_use_spat) { size_t stride = 0; @@ -842,8 +1313,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 021247281..7479e15d0 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, // 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, @@ -62,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 @@ -167,11 +203,28 @@ class PWMScorer : public GenomeSeqScorer // Utilities inline float get_spatial_log_factor(size_t pos_index) const; + // 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, 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, + int chromid, int64_t interval_start, int64_t interval_end, + bool lse_aggregate, bool ism); + // 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; @@ -180,6 +233,7 @@ class PWMScorer : public GenomeSeqScorer // Cache SlideCache m_slide; SpatSlideCache m_spat_slide; + GradSlideCache m_grad_slide; }; #endif // PWM_SCORER_H_ 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/src/TrackExpressionVars.cpp b/src/TrackExpressionVars.cpp index 40e85e3fc..4fe5c9ce1 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,48 @@ 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)) { + 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") { + 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) { diff --git a/tests/testthat/helper-pwm-grad-oracle.R b/tests/testthat/helper-pwm-grad-oracle.R new file mode 100644 index 000000000..791118601 --- /dev/null +++ b/tests/testthat/helper-pwm-grad-oracle.R @@ -0,0 +1,215 @@ +# 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") +} + +# 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, + 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 (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) + + # 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, + 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) + } + + 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..328fa2885 --- /dev/null +++ b/tests/testthat/test-vtrack-pwm-grad.R @@ -0,0 +1,1105 @@ +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) + } + } +}) + +# ---- 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 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))) + + # 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 = iter_W - 1) + res <- gextract("g_bidir", iterator = 1, intervals = pivot) + + 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 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.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 = iter_W - 1) + res <- gextract("g_minus", iterator = 1, intervals = pivot) + + 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 + +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) +}) + +# 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) +}) + +# 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) +}) + +# 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) +}) + +# ============================================================ +# 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) +}) diff --git a/vignettes/PWM-Functions.Rmd b/vignettes/PWM-Functions.Rmd new file mode 100644 index 000000000..68aa0e5ea --- /dev/null +++ b/vignettes/PWM-Functions.Rmd @@ -0,0 +1,445 @@ +--- +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 aggregating +scorers (`pwm`, `pwm.max`, `pwm.max.pos`, `pwm.count`), the +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 +they always evaluate. 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 forward-strand per-anchor score is + +$$ +\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. + +```{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 +``` + +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) +``` + +## 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 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. + +```{r} +max(all_anchor_scores(M, seq)) +``` + +### `pwm` (LSE): integrated likelihood + +`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 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) { + 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 similar-strength hits `pwm` can be substantially +larger. + +### `pwm.max.pos`: argmax position + +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 per-anchor +score is the strand union, 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) +``` + +## Common parameters + +These mirror what `?gvtrack.create` documents; one-liner each: + +- **`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; + 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 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 +**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. + +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`, +and only column 0 of the PSSM is involved. + +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}}). +$$ + +`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 anchors. The diff term equals the column-0 information +content under an "actual minus worst" baseline. + +For the MAX aggregate: + +$$ +g_{\text{MAX}}(p) \;=\; +\begin{cases} +M[0,\, b_p] - \min_b M[0,\, b] & \text{if } a^* = p,\\ +0 & \text{otherwise.} +\end{cases} +$$ + +i.e. the column-0 diff iff the head anchor is the argmax of the window. + +### 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). +$$ + +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 the two definitions differ + +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 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. + +`pwm.grad` is cheaper; `pwm.grad.ism` is more conservative around +argmax-shift cases. + +### Strands + +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( + 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) +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) +b_p <- "A" +diff_col0 <- M[1, b_p] - min(M[1, ]) +g_lse <- w_p * diff_col0 +c(w_p = w_p, diff_col0 = diff_col0, g_lse = g_lse) +``` + +ISM gradient: flip `seq[1]` to each of C/G/T and recompute the LSE. + +```{r} +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)) +ism_lse(M, seq) - min(f_alt) +``` + +The two values agree to within numerical noise; this is the typical +case. + +### Genome-track recipe + +```{r, eval = FALSE} +W <- 500 # window length + +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") + +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) +``` + +`out$grad - out$grad_ism` flags positions where the argmax flips under +single substitution. Multiplying `out$grad` by another per-bp track +gives a contribution-weighted version of that track. + +## Performance notes + +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 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.