nx_sequence_ambig.nx source
↩ module page · 151 lines · 6484 B
1// nx_sequence_ambig.nx -- N / IUPAC ambiguity sidecar for nx_sequence.
2//
3// license_tier: INDEPENDENT_REDERIVE
4// genealogy_id: international-research-sources/iupac-1985-nucleic-acid-nomenclature
5//
6// G0.5a of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. Closes the
7// N-handling gap in nx_sequence: G0.1 dna_pack refuses any non-ACGT
8// input (returns -1), which is correct for clean fixtures but
9// blocks real-world FASTQ reads where N is the universal "unknown
10// base" marker (Illumina low-quality cycles, ONT undetermined,
11// PacBio CCS hard-mask, reference-genome assembly gaps).
12//
13// Design -- parallel sidecar bitset:
14// The 2-bit packed bases array stays the same as G0.1 -- existing
15// primitives (revcomp / k-mer / FM-index / alignment) read it
16// unchanged. N positions are encoded as code-A (0b00) in the bases
17// array AND have the corresponding bit set in n_bits[]. Caller
18// checks dna_is_n at any position to know whether the A code is
19// real or a stand-in for N.
20//
21// Why a sidecar bitset, not 4-bit IUPAC codes throughout:
22// - Preserves bit-for-bit compatibility with G0.1 primitives
23// - Memory: 1 bit per base vs 4 bits = same 4x compression factor
24// applied to ambiguity tracking
25// - Most positions are NOT ambiguous; sidecar is mostly zeros
26// and compresses well in any downstream blob storage
27// - Full IUPAC (R Y K M S W B D H V) is rare outside reference
28// ambiguity codes and lives in G0.5b nx_sequence_iupac (a
29// separate primitive that uses 4-bit codes per position)
30//
31// Bitset layout (matches 2-bit DNA packing convention):
32// bit (7 - pos%8) of n_bits[pos/8] holds the N-flag for base pos.
33// MSB-first within byte to mirror the bases packing.
34// Position 0 is in the high bit (bit 7) of n_bits[0].
35//
36// API:
37// dna_pack_ambig(ascii, len, out_bases, out_nbits) -> i64
38// dna_unpack_ambig(bases, nbits, len, out_ascii) -> i64
39// dna_is_n(nbits, pos) -> i64
40// dna_set_n(nbits, pos, is_n) -> i64
41//
42// out_bases capacity: (len + 3) / 4 bytes (same as G0.1)
43// out_nbits capacity: (len + 7) / 8 bytes (1 bit per position)
44// Both must be pre-zeroed (sys_mmap returns zeroed pages).
45//
46// What G0.5a does NOT do (deferred):
47// - Full IUPAC ambiguity codes R/Y/K/M/S/W/B/D/H/V -- G0.5b
48// - U (uracil, RNA) -- folded into T via existing dna_encode_base path
49// - Soft-masked lowercase (repeat-masker convention) -- G0.5c
50// - Quality-score sidecar (Phred per position) -- nx_qphred (G2.0a)
51//
52// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml)
53// intended_use: "Real-world FASTQ + reference-genome N
54// handling; the gap-fill that lets every
55// upstream primitive (G0.1 / G0.3 / G1.x)
56// ingest sequencer output without rejection"
57// sil_target: SIL2
58// asil_target: QM
59// dal_target: DAL C
60// iec_62304_class: B
61// evidence: [no_floating_point, deterministic,
62// bit_equal_reproducible,
63// msb_first_bitset_matches_2bit_packing,
64// pack_unpack_round_trip_with_N_KAT,
65// non_IUPAC_refused_explicit_minus_one,
66// license_tier_INDEPENDENT_REDERIVE]
67// hazard_register: [bug-tape-nbits-bit-order-mismatch-with-bases,
68// bug-tape-N-treated-as-A-without-sidecar-check,
69// bug-tape-iupac-silently-accepted-as-N,
70// bug-tape-sidecar-too-small-overrun]
71// residual_risk: "G0.5a treats N as 'unknown' only. IUPAC
72// partial-match codes (Y = C or T, etc.)
73// are REFUSED here and require G0.5b's
74// 4-bit-per-position layout for correctness
75// in variant calling and primer design."
76// verdict: NOT_YET_EVALUATED
77
78import "nx_syscalls.nx"
79import "nx_const.nx"
80import "nx_sequence.nx"
81
82// Read the N-flag for a given base position.
83// Returns 1 if N (ambiguous), 0 otherwise.
84func dna_is_n(nbits: *u8, pos: i64) -> i64 {
85 let byte_idx: i64 = pos >> 3
86 let sub: i64 = pos & 7
87 let shift: i64 = 7 - sub
88 let b: i64 = nbits[byte_idx] & 0xff
89 return (b >> shift) & 1
90}
91
92// Write the N-flag for a given base position.
93// is_n = 1 sets the bit; is_n = 0 clears it.
94func dna_set_n(nbits: *u8, pos: i64, is_n: i64) -> i64 {
95 let byte_idx: i64 = pos >> 3
96 let sub: i64 = pos & 7
97 let shift: i64 = 7 - sub
98 let mask: i64 = 1 << shift
99 let cur: i64 = nbits[byte_idx] & 0xff
100 if (is_n & 1) == 1 {
101 nbits[byte_idx] = (cur | mask) & 0xff
102 } else {
103 let cleared: i64 = cur & (mask ^ 0xff)
104 nbits[byte_idx] = cleared & 0xff
105 }
106 return 0
107}
108
109// Pack ASCII bases (allowing A/C/G/T + N) into 2-bit form + N sidecar.
110// Returns 0 on success, -1 if any base is not in {A,C,G,T,N,a,c,g,t,n}.
111// N positions: bases[pos] gets code A (00); nbits[pos] gets bit set.
112func dna_pack_ambig(ascii: *u8, len: i64,
113 out_bases: *u8, out_nbits: *u8) -> i64 {
114 if len < 0 { return -1 }
115 var i: i64 = 0
116 while i < len {
117 let c: i64 = ascii[i] & 0xff
118 var is_n: i64 = 0
119 if c == NX_ASCII_DNA_N { is_n = 1 }
120 if c == 0x6E { is_n = 1 } // 'n' lowercase
121
122 if is_n == 1 {
123 dna_set_base(out_bases, i, NX_DNA_A)
124 dna_set_n(out_nbits, i, 1)
125 } else {
126 let code: i64 = dna_encode_base(c)
127 if code < 0 { return -1 } // refuses other IUPAC + bad bytes
128 dna_set_base(out_bases, i, code)
129 // dna_set_n with 0 is a no-op on a zeroed buffer; skipped
130 // intentionally for the hot path (assumes caller-zeroed).
131 }
132 i = i + 1
133 }
134 return 0
135}
136
137// Unpack 2-bit bases + N sidecar back to ASCII, restoring 'N' at
138// flagged positions. out_ascii capacity >= len; no terminator written.
139func dna_unpack_ambig(bases: *u8, nbits: *u8, len: i64, out_ascii: *u8) -> i64 {
140 var i: i64 = 0
141 while i < len {
142 if dna_is_n(nbits, i) == 1 {
143 out_ascii[i] = NX_ASCII_DNA_N & 0xff
144 } else {
145 let code: i64 = dna_get_base(bases, i)
146 out_ascii[i] = dna_decode_base(code) & 0xff
147 }
148 i = i + 1
149 }
150 return 0
151}