Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 24 additions & 10 deletions src/intronIC/scoring/background.py
Original file line number Diff line number Diff line change
Expand Up @@ -464,20 +464,27 @@ def _build_final_pwm_sets(
if human_ps is None:
continue

# Iterate subtypes in sorted order so that when multiple
# subtypes map to the same pwm_subtype (e.g., non-canonical 5'
# dnts all default to 'gtag'), the last-writer-wins assignment
# to ``matrices[human_key]`` is deterministic. Without this,
# streaming and in-memory pipelines can produce different
# corrected PWMs for the same data because the underlying
# ``_data`` dict's insertion order differs (sequential
# accumulate vs parallel merge_worker_counts).
for subtype_dnt in sorted(acc.get_all_subtypes()):
# Multiple intron dnts map to the same pwm_subtype (all non-canonical
# 5' dnts default to 'gtag'). The HIGHEST-n dnt must define the slot.
# The previous code iterated sorted(dnts) and did last-writer-wins,
# which let a rare low-n non-canonical dnt (e.g. TT n=3) clobber the
# real GT background (n~96k) with a ~human matrix (w=n/(n+n0)~=0),
# silently disabling the species correction on any annotation that
# carries non-canonical dnts. Sort by (-n, dnt) [order-independent =>
# streaming/in-memory deterministic, the original goal] + first-writer-wins.
_u2_seen = set()
for subtype_dnt in sorted(
acc.get_all_subtypes(),
key=lambda d: (-acc._data.get(d, {}).get('n', 0), d),
):
pwm_subtype = self.FIVE_DNT_TO_SUBTYPE.get(
subtype_dnt,
self.THREE_DNT_TO_SUBTYPE.get(subtype_dnt, 'gtag')
)
human_key = ('u2', pwm_subtype)
if human_key in _u2_seen:
continue # a higher-n dnt already defined this background
_u2_seen.add(human_key)
human_pwm = human_ps.matrices.get(human_key)
if human_pwm is None:
# Try fallback
Expand Down Expand Up @@ -558,9 +565,16 @@ def _build_final_pwm_sets(

human_bp = self.human_u2.get('bp')
if human_bp:
for subtype_dnt in sorted(self._bp_acc.get_all_subtypes()):
_bp_u2_seen = set()
for subtype_dnt in sorted(
self._bp_acc.get_all_subtypes(),
key=lambda d: (-self._bp_acc._data.get(d, {}).get('n', 0), d),
):
pwm_subtype = self.FIVE_DNT_TO_SUBTYPE.get(subtype_dnt, 'gtag')
human_key = ('u2', pwm_subtype)
if human_key in _bp_u2_seen:
continue # highest-n dnt wins (see 5'/3' loop; clobber-bug fix)
_bp_u2_seen.add(human_key)
human_pwm = human_bp.matrices.get(human_key)
if human_pwm is None:
human_pwm = human_bp.select_best('u2', pwm_subtype)
Expand Down
32 changes: 32 additions & 0 deletions tests/unit/test_scoring/test_background.py
Original file line number Diff line number Diff line change
Expand Up @@ -209,3 +209,35 @@ def test_trim_excludes_high_scorers(self, human_u2_pwm_sets):
a_freq = u2_five.matrix[0, 17]
g_freq = u2_five.matrix[2, 17]
assert a_freq > g_freq

def test_low_n_noncanonical_dnt_does_not_clobber(
self, human_u2_pwm_sets, default_config
):
"""Regression (U2 subtype-clobber bug): a rare non-canonical 5' dnt that
defaults to 'gtag' (e.g. TT) must NOT overwrite the high-n GT-derived gtag
background. The previous code sorted dnts and did last-writer-wins, so a
handful of TT introns (n=3, blend weight ~0 => ~human prior) clobbered the
GT background (n=2000), silently disabling the species correction on any
annotation carrying non-canonical dnts. Fix: highest-n dnt defines the slot."""
bg = SpeciesBackground(
human_u2_pwm_sets, default_config,
five_len=12, three_len=10, bp_len=12,
)
# 2000 canonical GT-AG introns, 5' window all-A => strongly species-shifted
# away from the uniform (0.25) human prior.
for i in range(2000):
bg.accumulate(f'gt_{i}', 'GT', 'AG', 'A' * 12, 'C' * 10, 'A' * 50)
# 3 non-canonical TT introns (5' dnt 'TT' -> defaults to 'gtag'), 5' all-G.
# 'TT' > 'GT' alphabetically, so last-writer-wins would let these win.
for i in range(3):
bg.accumulate(f'tt_{i}', 'TT', 'AG', 'G' * 12, 'C' * 10, 'A' * 50)

five = bg.build_corrected_pwm_sets()['five'].matrices[('u2', 'gtag')].matrix
# Scored window is cols 17-28, A is row 0. Under the fix the gtag background
# reflects the 2000 GT introns (A ~0.75); under the clobber bug it reverts
# to ~human uniform (~0.25).
a_scored = float(five[0, 17:29].mean())
assert a_scored > 0.5, (
f"gtag 5' A-freq={a_scored:.3f} in scored window: expected species-shifted "
f"(~0.75 from 2000 GT introns); ~0.25 means 3 TT introns clobbered it to human"
)
Loading