Skip to content

Handle IUPAC ambiguity codes in the reference (#291) - #309

Merged
joshfactorial merged 2 commits into
developfrom
fix/291-iupac-handling
Jun 8, 2026
Merged

Handle IUPAC ambiguity codes in the reference (#291)#309
joshfactorial merged 2 commits into
developfrom
fix/291-iupac-handling

Conversation

@joshfactorial

@joshfactorial joshfactorial commented Jun 8, 2026

Copy link
Copy Markdown
Collaborator

Fixes #291. Release v4.5.3 (version bump + ChangeLog entry included).

Problem

GRCh38 and other assemblies carry IUPAC ambiguity codes (R Y S W K M B D H V) alongside A/C/G/T/N. NEAT had no handling:

  • They survived reference load — Read.convert_masking() only neutralized non-ACGT bases when the segment also contained an N, and skipped bases before the first N.
  • A surviving code at a sequencing-error site hit NUC_IND[snv_reference] (error_models.py:256); NUC_IND only has A/C/G/T → KeyError, aborting the run.

(Mutation generation was already safe — skips non-ALLOWED_NUCL positions — and the BAM SEQ encoder already supports IUPAC nibbles. The crash was specific to the sequencing-error path.)

Fix

Resolve ambiguity codes to a concrete base at reference load, so every downstream consumer (reads, BAM, golden VCF, error/trinucleotide lookups) only ever sees A/C/G/T/N.

  • common: IUPAC_CODES map + resolve_iupac_bases() — replaces each code with one of the bases it represents using the run's seeded RNG (reproducible). Vectorized; the no-codes common case is a single scan, so genome-scale references aren't penalized. N is left for its existing low-quality masking.
  • split_inputs: resolves right after the existing .upper(); logs a warning with the count.
  • error_models: defense-in-depth — NUC_IND lookup skips a non-ACGT base instead of raising, so nothing can crash even if a code slips through.

Release

(Branch is rebased on current develop at 4.5.2, so the bump applies cleanly.)

Tests

  • Resolver unit tests: every code resolves into its allowed set, counts, reproducibility, N left alone, length preserved, lowercase passthrough.
  • End-to-end runner test on a reference with all ten codes + an N run, asserting FASTQ output is clean.

Full suite: 790 passed, 3 skipped. Verified end-to-end (single- and multi-threaded) on a synthetic IUPAC reference: no crash; FASTQ/BAM/VCF outputs contain only A/C/G/T/N.

Known limitation

The input-VCF ref-match check (vcf_func.py:137) reads the original reference file, not the resolved chunks — it won't crash (string comparison), but at an ambiguous position an input variant is matched against the original code while reads show the resolved base. Rare edge case; left as follow-up.

🤖 Generated with Claude Code

joshfactorial and others added 2 commits June 8, 2026 16:50
GRCh38 and other assemblies carry IUPAC ambiguity codes (R, Y, S, W, K,
M, B, D, H, V). NEAT previously had no handling for them:
- they survived reference load (convert_masking only neutralized non-ACGT
  bases when an N was also present in the segment, and skipped bases
  before the first N), and
- a surviving code at a sequencing-error site hit NUC_IND[base] and
  raised KeyError, crashing the run.

Fix: resolve ambiguity codes to a concrete base at reference load.
- common: add IUPAC_CODES map + resolve_iupac_bases(), which replaces
  each code with one of the bases it represents using the run's seeded
  RNG (reproducible). Vectorized; the no-codes common case is a single
  scan, so genome-scale references aren't penalized. 'N' is left for its
  existing low-quality masking.
- split_inputs: resolve right after the existing upper(), so every
  downstream consumer (reads, BAM, golden VCF, error/trinucleotide
  lookups) only ever sees A/C/G/T/N. Logs a warning with the count.
- error_models: defense-in-depth — NUC_IND lookup now skips a non-ACGT
  base instead of raising, so nothing can crash even if a code slips
  through.

Tests: resolver unit tests (every code, counts, reproducibility, N left
alone, length preserved) and an end-to-end runner test on a reference
containing all ten codes plus an N run, asserting FASTQ output is clean.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Bump version to 4.5.3 and add the ChangeLog entry for the IUPAC
ambiguity-code handling fix.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
@joshfactorial
joshfactorial force-pushed the fix/291-iupac-handling branch from fa28705 to ee55d91 Compare June 8, 2026 21:51
@joshfactorial
joshfactorial merged commit fcbf898 into develop Jun 8, 2026
1 check passed
@joshfactorial
joshfactorial deleted the fix/291-iupac-handling branch June 8, 2026 21:56
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant