Hairpin behavior: 25C fold, multiloop splitting, degenerate-dS/Tm guards - #14
Open
KowalskiBio wants to merge 12 commits into
Open
Hairpin behavior: 25C fold, multiloop splitting, degenerate-dS/Tm guards#14KowalskiBio wants to merge 12 commits into
KowalskiBio wants to merge 12 commits into
Conversation
Lib/strider._native gains dimer_mfe_candidates, a statement-by-statement Rust port of strider.thermo.dimer_thermo._dimer_mfe_candidates (DNA): the inter-strand DP fill, dangles, terminal penalties, and the full Mathews-Turner interior-energy set (exact 1x1/1x2/2x2 tables included). Sequence-independent table data is NOT hand-copied: native/codegen_tables.py packs strider.thermo.parameters_dna keys (2-bit-per-base) and emits sorted static arrays (native/src/tables_dna.rs, 16 tables / 12694 entries) that the Rust side binary-searches. strider.thermo.dimer_thermo prefers the native DP when material='dna' and no param override is active (param_context or engine custom ParameterSet); all other cases keep the pure-Python implementation, which also stays as the test oracle (_dimer_mfe_candidates_py). Parity (tests/test_native_dp_parity.py): * candidate-list fuzz: 600 random pairs, 0/600 energy or pair mismatches * end-to-end: dimer_thermo/dimer_thermo_subopt (incl. ValueError behavior) 300/300 identical with the DP toggled between native and Python Timing on the Oligool IDT-analysis workload (2 hairpins + 3 dimer pairs): 41.5 ms -> 32.1 ms wall; the dimer DP segment (32 ms cumulative) drops to ~2 ms. Full suite unchanged: 649 passed, 4 pre-existing ViennaRNA-env failures, no regressions.
…ram-override guard test Addresses review feedback on the DP port: Alphabet divergence (EmilioVenegas#1): * tables.rs: CODE_TABLE defaults unknown bytes to u32::MAX (not 3=T), so any degenerate base (N/R/Y/...) produces a packed code that never matches a generated entry — binary-search misses, lookup returns None, the caller uses the Python default, exactly like dict.get(key, default) * pack(): early-return on sentinel so the check costs one branch per key * dimer_thermo.py: alphabet guard in the native wrapper rejects non-ACGT sequences and falls back to Python — the primary gate * tests: degenerate-base sequences now in the fuzz corpus (15%); explicit N/R cases in test_fallback_guard_degenerate_alphabet; test_native_raw_degenerate_matches_python proves the sentinel fix works even for direct (wrapper-bypassing) raw calls Codegen staleness (EmilioVenegas#2): * codegen_tables.py docstring now correctly references the test that validates the committed tables (it previously claimed a check that did not exist) * test_tables_dna_regenerated_matches_committed: regenerates from parameters_dna and compares byte-for-byte; catches silent drift when the Python tables are edited but the Rust ones were not regenerated Param-override guard (EmilioVenegas#3): * test_fallback_guard_param_override: _param_override with an empty dG dict forces the Python path (mirroring lookup_table fallback) * test_fallback_guard_custom_engine_params: engine with _uses_custom_params() = True forces the Python path 8/8 new + existing postgres tests pass; full suite unchanged (654 passed, 4 pre-existing ViennaRNA-environment failures).
The native-parity job only ran tests/test_native_parity.py; the DP fuzz, the guard tests and the codegen staleness check live in tests/test_native_dp_parity.py, which never executed in CI. Run both.
…s + dH interior tables - convert_stack: the .par matrix's column pair is written in reversed orientation (from the j-side); keys must read outer[0] + inner[1] + inner[0] + outer[1]. Verified against RNA.eval_structure_verbose: GCGC/CGCG = -2.2, CCGG/GGCC = -1.8, AATT/TTAA = -1.0, ATAT = -0.9, TATA = -0.6. - convert_bulge_size: drop the leading 0.0 'sentinel' — index 0 IS the 1-nt bulge entry (2.9). The bug made single-base bulges free and shifted every larger bulge by one size — the root cause of bogus single-bulge hairpin structures the DP over-stabilized. - convert_interior_size: same off-by-one fix (index i = VR internal[i+1]); size-table sentinels are now 30.0 (INF), matching the native set. - generate_paramset: emit the previously computed-but-dropped interior_1_1 / interior_2_2 dH tables, so structure_enthalpy no longer silently mixes native dG values into the enthalpy walk. - verification: add bulge-containing hairpins to the test set — the old straight-stem-only set passed with the bug present. 7/7 structures now match ViennaRNA eval_structure exactly at dangles=0.
…2 exterior stacks
- dimer_thermo / dimer_thermo_subopt / dimer_tm accept paramset (instance
or name) and use it to predict AND score the inter-strand helix,
matching the hairpin API.
- param_context + structure_* walkers resolve parameter-set NAMES via a
memoized loader instead of crashing on str ('str' has no .dG).
- hairpin_thermo folds with the same paramset it scores with
(param_context wrapping fold_mfe); fraction_folded returns 0.0 for
non-folding sequences instead of raising ValueError.
- NEW: dangles option (default 0, bit-identical energies) on fold_mfe /
subopt_structures / fold_complex / ThermoEngine(dangles=...) — exterior
dangling-end stacks at W-matrix stem placements: best negative single
dangle per exterior stem. Verified against ViennaRNA eval on a 25-case
single-stem flank grid: single-flank cases match exactly, two-flank
cases track the better of the two flanks; 38/42 end-to-end sequences
within 0.4 kcal/mol of VR(dangles=2) (remaining are G.T-wobble stems).
- hairpin dG/dH structure walk applies the same exterior-tail dangle, so
hairpin_thermo stays energy-consistent with the dangle-aware MFE.
…2004 paramset scripts/generate_mathews2004_params.py uses ViennaRNA's distributed dna_mathews2004.par as a development-time extraction input to produce strider/thermo/parameters/mathews2004-dna.json; the transcribed values are primary literals of Mathews et al. 1999 (JMB 288:911-940) — physical measurements, not copyrightable expression — and the JSON is self-contained static data. Neither the .par file nor ViennaRNA is redistributed or required at runtime.
…er flag, tests, notices * fraction_folded: zero case derives from MFE pairs; caller mistakes raise again (bad salt model, degenerate entropy, multiloop MFE) instead of being flattened to 0.0 by the blanket ValueError catch. * dangles scope: engine trim to MFE/suboptimal - pfunc/ensemble/defect/ differentiable/equilibrium explicitly ignore the flag; dropped from the pfunc cache key. Naming note at the API surface for the dangles=2 vs VR dangles=2 divergence (best single negative dangle per exterior stem). * dimer path: dimer_thermo/_subopt/_tm accept dangles=0|2 threaded into _sum_dimer_elements (each terminus gathers adjacent negative dangle stacks). Also a latent-domain fix: the enthalpy walk no longer ingests DG-valued dangle tables (dangle tables are DG-only). * ext_e 5' branch: drop the redundant (ko - 1) not in nicks check that discarded valid strand-start 5' dangles in fold_mfe/fold_complex. * notices: value-level confirmation of the dna_mathews2004.par identity (ATAT -0.9 vs -0.8, GT/CG +1.2 vs +1.3 markers vs the 1999 par file). * _resolve_name: documented the shared-instance read-only contract (no in-place writes exist in the tree today). * tests/test_mathews2004.py (12 tests): pinned mathews2004-dna folded energies and table markers, 42-case dangles=0 bit-identity grid, name-path regression, fraction_folded zero + error propagation, dimer dangles flag, pfunc/cache-key scoping, strand-start flank regression pin. * changelog: behavior-changes and documentation entries for the above.
ViennaRNA marks 8 hairpin-mismatch dG entries as DEF but resolves them to -0.50 kcal/mol at runtime. strider was treating DEF as 0.0, causing dG mismatches of 0.50 kcal/mol on any hairpin with those mismatch combinations (verified against eval_structure on 10k oligos). The dH values for these entries ARE defined in dna_mathews2004.par but were also discarded. Now always stored when defined (96 entries, up from 73). Hardcode the 8 known DEF dG entries to -0.50 via a lookup set, since the .par file does not encode the default and ViennaRNA resolves it internally. All other DEF entries remain 0.0 (the correct default for non-WC pair types). Regenerated mathews2004-dna.json passes all verification checks against ViennaRNA eval_structure at dangles=0.
ensemble_dg and multistrand_pairs now accept pair_probs=True (default, preserves existing behavior). When False, skip _external_Q_nodangle and _pair_probs_outside and return a zero matrix; dG is byte-identical because the inside partition (_fill_dp_nicks + _apply_coaxial_external) fully determines Q[0][n-1] before the outside pass runs, and the outside pass writes only to local arrays. Thread the flag through ThermoEngine.pfunc -> _pfunc_dispatch -> _pfunc_native -> _pfunc_native_inner. Use a distinct 'pfunc_dg' cache-key op for the no-probs variant so cached results never cross-contaminate, and extend the dangles-exclusion check to cover it.
…splitting, dS threshold, and 25C fold Three fixes for hairpin_thermo: 1. dS threshold: raise the degenerate-entropy guard from exact-zero to |dS| < 0.5 cal/mol/K, and add a Tm > 200C ceiling. Both return None (via ValueError) instead of absurd values like 528C. The two-state model is not applicable when dS is near-zero or Tm exceeds 200C. 2. Multiloop splitting: when parse_hairpin_pairs returns None (multiloop), try _split_stem_groups to decompose into individual stem-loops and recurse on the most stable one. Mirrors the splitStemGroups algorithm in Oligool's HairpinSVG.tsx. Recovers Tm for GC-rich primers whose MFE is a multiloop but contain a valid hairpin stem. 3. Fold at 25C via ThermoEngine: replace fold_mfe(seq, 37C, ...) with ThermoEngine(celsius=25.0) for the internal MFE. At 37C, marginally stable hairpins (Tm 25-37C) unfold and are missed, yielding NA Tm. ThermoEngine applies dG(T) = dH - T*dS, which at 25C finds weak structures. Salt is NOT passed (folds at 1M Na+ reference) so the MFE structure is salt-independent, as the Owczarzy model requires. Results on 10k oligo parity table: - Extreme Tms (>200C): eliminated (was ~15 sequences) - Mathews NA rate: 51.3% -> 26.5% - SantaLucia NA rate: 49.8% -> 35.9% - Mathews mean abs diff vs IDT: 7.36 -> 6.50C - All 660 strider tests pass (excluding pre-existing dimer failure)
…split bug Fixes issues a reviewer would catch on the previous commit (70ec8e9): 1. Dead import: remove unused fold_mfe import from hairpin_thermo (ThermoEngine handles folding now). 2. fraction_folded consistency: update fraction_folded to use ThermoEngine at 25C (FOLD_CELSIUS) instead of fold_mfe at 37C, so both functions find the same MFE structure. 3. Redundant exception: except (ValueError, Exception) -> except Exception. 4. Named constants: replace magic numbers 0.0005 and 200.0 with MIN_DS_CAL (0.5 cal/mol/K) and MAX_TM_CELSIUS (200C), with physical motivation documented. Add FOLD_CELSIUS (25C) constant. 5. _split_stem_groups bug: end position used g[-1][1] (innermost pair closing) instead of g[0][1] (outermost pair closing), truncating sub-structures. Fixed. 6. Tests: add 7 new tests covering dS threshold, Tm ceiling, 25C fold behavior, multiloop splitting correctness, and edge cases for _split_stem_groups (single hairpin, unbalanced, no pairs). All 667 tests pass (was 660 + 7 new).
1. except Exception: continue -> except ValueError: continue in the multiloop stem-loop candidate loop. A bare Exception catch masks real errors (e.g. ImportError, TypeError) behind 'no valid stem-loop found in multiloop'. ValueError is what hairpin_thermo raises for degenerate dS, extreme Tm, or non-hairpin structures. 2. test_degenerate_ds_raises now tests behavior: patches the energy walkers to force near-zero dS and asserts ValueError is raised. Previously it just asserted MIN_DS_CAL == 0.5 and MAX_TM_CELSIUS == 200.0, restating the constants rather than testing the guard. 3. Expose FOLD_CELSIUS as a named fold_celsius parameter on hairpin_thermo / hairpin_tm / fraction_folded (default 25.0C). The 25C fold matches IDT OligoAnalyzer's reporting temperature. Pass fold_celsius=37.0 to restore the pre-25C behavior.
Owner
|
Read the parity report, solid work on the methodology. You were right on the outliers: breaking down the three worst cases clearly shows distinct causes across implementations rather than a shared comparison flaw. The common-subset result also settles #13 cleanly (Mathews +0.29 vs SantaLucia +1.19 ΔG error on 4,309 shared oligos), showing the gain isn't an NA filtering artifact. A few points on the PR:
Next steps on the gate: because this benchmark compares predictions against predictions, I will validate the old vs new defaults against experimental Tm data (Vallone set) locally and post the results here. Defaulting |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Stacked on #13 (
mathews2004-dangles). Until #13 merges, this PR's diff also shows #13's commits; the effective diff here is the 3 commits at the tip.What changed
Three behavioral changes to
hairpin_thermo, split out from #13 per review:1. Fold at 25 C instead of 37 C
hairpin_thermonow folds the internal MFE at 25 C (matching IDT OligoAnalyzer's reporting temperature) instead of the 37 C dG-table reference. At 37 C, marginally stable hairpins (Tm 25-37 C) unfold and are missed, yielding NA Tm. The 25 C choice is because IDT reports hairpin analysis at 25 C, which is a stronger justification than the NA-rate argument in the original commit message.FOLD_CELSIUSis now exposed as a namedfold_celsiusparameter onhairpin_thermo,hairpin_tm, andfraction_folded(default 25.0). Passfold_celsius=37.0to restore the pre-25 C behavior.2. Multiloop splitting
When the MFE is a multiloop,
hairpin_thermonow decomposes it into individual stem-loops and scores the most stable one, instead of raisingValueError. Recovers Tm for GC-rich primers whose MFE is a multiloop but contain a valid hairpin stem.3. Degenerate dS / extreme Tm guards
Near-zero |dS| (< 0.5 cal/mol/K) raises
ValueErrorinstead of returning an absurd Tm. Tm above 200 C raisesValueError(eliminates ~15 sequences with extreme Tms on the 10k parity table).Review notes addressed
except Exception: continuein the multiloop candidate loop changed toexcept ValueError: continue. The bareExceptioncatch masked real errors (ImportError, TypeError) behind "no valid stem-loop found in multiloop".test_degenerate_ds_raisesnow tests behavior: patches the energy walkers to force near-zero dS and assertsValueErroris raised. Previously it assertedMIN_DS_CAL == 0.5andMAX_TM_CELSIUS == 200.0, restating constants.Evidence gap
The 10k oligo parity table measures agreement with other prediction tools (IDT, primer3, ViennaRNA). Tuning defaults toward IDT can move away from experimental agreement. The Vallone hairpin set (experimentally measured Tm) is not in the repo and should be run under old (37 C fold, no guards) and new (25 C fold, guards) behavior side by side before this PR lands. If the new defaults hold up there, they land with evidence behind them.
Test plan
_split_stem_groupsedge cases.