nx_sequence.nx source

1// nx_sequence.nx -- foundational DNA / RNA sequence algebra primitive. 2// 3// license_tier: INDEPENDENT_REDERIVE 4// genealogy_id: international-research-sources/biology/dna_2bit_encoding_field_standard 5// 6// G0.1 of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. Encoding + 7// reverse-complement + k-mer enumeration + canonical k-mer. Built 8// bits-up against existing nx_syscalls; no other genomics dependency. 9// 10// 2-bit packed DNA encoding: 11// A = 00, C = 01, G = 10, T = 11 12// Complement = code XOR 0b11 (A<->T, C<->G). 13// 4 bases per byte, packed big-endian within byte: 14// byte bits 7-6 = base position 0 within byte 15// byte bits 5-4 = base position 1 16// byte bits 3-2 = base position 2 17// byte bits 1-0 = base position 3 18// Bytes ordered first-to-last. A trailing partial byte zero-pads 19// the low bits (so "ACG" = 0x18 = 00 01 10 00 in one byte). 20// 21// k-mer encoding (k <= 32): 22// Single i64 holding the 2*k low bits. Base position 0 of the k-mer 23// sits in the highest 2 bits of the 2*k-bit window. This means 24// lexicographic ordering of bases (A < C < G < T) matches arithmetic 25// ordering of the packed i64, so canonical k-mer = min(fwd, rc) is 26// a simple arithmetic min. 27// 28// Why 2-bit: 29// - Minimum storage for an alphabet of 4 (human reference ~750 MB 30// uncompressed FASTA -> ~190 MB packed) 31// - Reverse-complement of a packed code is a single XOR 32// - SIMD-friendly: 16 bases per 32-bit word, 32 bases per i64 33// - k-mer hashing fits naturally 34// 35// N / ambiguity handling: 36// Out of scope for G0.1. See nx_sequence_ambig.nx (G0.5) for the 37// IUPAC sidecar bitset. In G0.1, callers must filter Ns upstream 38// or the encoder returns -1 (refused) on any non-ACGT input. 39// 40// API: 41// dna_encode_base(c) -> i64 ASCII -> 2-bit code or -1 42// dna_decode_base(code) -> i64 2-bit -> ASCII 43// dna_complement_code(code) -> i64 2-bit -> complemented 2-bit 44// dna_pack(ascii, len, out_bases) -> i64 pack ASCII into 2-bit; 0 ok / -1 refused 45// dna_unpack(bases, len, out) -> i64 unpack 2-bit into ASCII; 0 46// dna_get_base(bases, pos) -> i64 read 2-bit code at base pos 47// dna_set_base(bases, pos, code) -> i64 write 2-bit code at base pos 48// dna_revcomp(in, len, out) -> i64 reverse-complement packed seq; 0 49// dna_kmer_at(bases, pos, k) -> i64 extract packed k-mer at pos 50// dna_kmer_revcomp(kmer, k) -> i64 reverse-complement a packed k-mer 51// dna_kmer_canonical(kmer, k) -> i64 min(fwd, revcomp) -- canonical form 52// 53// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml) 54// intended_use: "DNA / RNA sequence algebra -- foundation for 55// alignment, variant calling, k-mer indexing, 56// CRISPR guide design, polygenic scoring" 57// sil_target: SIL2 (data-type primitive; correctness 58// errors here cascade into wrong 59// variant calls + wrong guide designs) 60// asil_target: QM 61// dal_target: DAL C 62// iec_62304_class: B 63// evidence: [no_floating_point, no_table_lookup, 64// bit_equal_reproducible, 65// 2bit_encoding_field_standard, 66// revcomp_round_trip_verified, 67// kmer_canonical_palindrome_verified, 68// constant_time_per_base, 69// license_tier_INDEPENDENT_REDERIVE] 70// hazard_register: [bug-tape-iupac-N-silently-dropped, 71// bug-tape-revcomp-partial-byte-misalign, 72// bug-tape-kmer-k-greater-than-32-overflow, 73// bug-tape-canonical-kmer-not-symmetric] 74// residual_risk: "G0.1 refuses non-ACGT input; callers using 75// long-read data with Ns must pre-filter or 76// wait for nx_sequence_ambig (G0.5). Partial- 77// byte alignment handled per-base, not packed 78// SIMD; G0.4 perf gate may require a packed 79// fast path." 80// verdict: NOT_YET_EVALUATED 81 82import "nx_syscalls.nx" 83import "nx_const.nx" 84 85// All DNA + ASCII + complement constants are sourced from nx_const 86// (NishiLang's `const` requires integer literals, so we use the 87// NX_ constants directly in function bodies rather than locally 88// re-aliasing them). Substrate-wide source of truth lives in 89// nx_const.nx per docs/NISHI_CONSTANT_REGISTRY.md. 90 91// ---- single-base conversions ---- 92 93// ASCII base -> 2-bit code. Returns -1 for any non-ACGT input 94// (refused per G0.1 scope; Ns handled in G0.5). 95func dna_encode_base(c: i64) -> i64 { 96 let cc: i64 = c & 0xff 97 if cc == NX_ASCII_DNA_A { return NX_DNA_A } 98 if cc == NX_ASCII_DNA_a { return NX_DNA_A } 99 if cc == NX_ASCII_DNA_C { return NX_DNA_C } 100 if cc == NX_ASCII_DNA_c { return NX_DNA_C } 101 if cc == NX_ASCII_DNA_G { return NX_DNA_G } 102 if cc == NX_ASCII_DNA_g { return NX_DNA_G } 103 if cc == NX_ASCII_DNA_T { return NX_DNA_T } 104 if cc == NX_ASCII_DNA_t { return NX_DNA_T } 105 return -1 106} 107 108// 2-bit code -> uppercase ASCII base. Caller responsibility to pass 109// a valid code 0..3; out-of-range returns 0x4E ('N') as a fail-safe. 110func dna_decode_base(code: i64) -> i64 { 111 let cc: i64 = code & 3 112 if cc == NX_DNA_A { return NX_ASCII_DNA_A } 113 if cc == NX_DNA_C { return NX_ASCII_DNA_C } 114 if cc == NX_DNA_G { return NX_ASCII_DNA_G } 115 if cc == NX_DNA_T { return NX_ASCII_DNA_T } 116 return 0x4E 117} 118 119// Complement of a 2-bit code: A<->T, C<->G. 120func dna_complement_code(code: i64) -> i64 { 121 return (code & 3) ^ NX_DNA_COMPLEMENT_MASK 122} 123 124// ---- per-base packed-array accessors ---- 125 126// Read 2-bit code at base position pos from a packed array. 127// byte index = pos / 4; sub-index = pos % 4 (0..3, 0 = high bits). 128// Shift amount = 6 - 2 * subindex. 129func dna_get_base(bases: *u8, pos: i64) -> i64 { 130 let byte_idx: i64 = pos >> 2 131 let sub: i64 = pos & 3 132 let shift: i64 = 6 - (sub << 1) 133 let b: i64 = bases[byte_idx] & 0xff 134 return (b >> shift) & 3 135} 136 137// Write 2-bit code at base position pos into a packed array. 138// Reads current byte, clears the target 2-bit slot, ORs in the new code. 139func dna_set_base(bases: *u8, pos: i64, code: i64) -> i64 { 140 let byte_idx: i64 = pos >> 2 141 let sub: i64 = pos & 3 142 let shift: i64 = 6 - (sub << 1) 143 let mask: i64 = 3 << shift 144 let cur: i64 = bases[byte_idx] & 0xff 145 let cleared: i64 = cur & (mask ^ 0xff) 146 let written: i64 = cleared | ((code & 3) << shift) 147 bases[byte_idx] = written & 0xff 148 return 0 149} 150 151// ---- whole-sequence pack / unpack ---- 152 153// Pack ASCII bases into 2-bit form. out_bases must have capacity 154// >= (len + 3) / 4 bytes, pre-zeroed (callers using sys_mmap get 155// zeroed pages). Returns 0 on success, -1 if any base is non-ACGT. 156func dna_pack(ascii: *u8, len: i64, out_bases: *u8) -> i64 { 157 var i: i64 = 0 158 while i < len { 159 let c: i64 = ascii[i] & 0xff 160 let code: i64 = dna_encode_base(c) 161 if code < 0 { return -1 } 162 dna_set_base(out_bases, i, code) 163 i = i + 1 164 } 165 return 0 166} 167 168// Unpack 2-bit bases back to ASCII. out must have capacity >= len. 169// No terminator written; caller wraps in length-prefix or sentinel. 170func dna_unpack(bases: *u8, len: i64, out: *u8) -> i64 { 171 var i: i64 = 0 172 while i < len { 173 let code: i64 = dna_get_base(bases, i) 174 out[i] = dna_decode_base(code) & 0xff 175 i = i + 1 176 } 177 return 0 178} 179 180// ---- reverse-complement ---- 181 182// Reverse-complement of a packed sequence. out_bases must have 183// capacity >= (len + 3) / 4 bytes, pre-zeroed. out_bases may NOT 184// alias in_bases. Walks base-by-base for correctness across the 185// partial-byte boundary; a packed SIMD fast path lives in G0.4 perf. 186func dna_revcomp(in_bases: *u8, len: i64, out_bases: *u8) -> i64 { 187 var i: i64 = 0 188 while i < len { 189 let fwd: i64 = dna_get_base(in_bases, i) 190 let rc: i64 = dna_complement_code(fwd) 191 dna_set_base(out_bases, len - 1 - i, rc) 192 i = i + 1 193 } 194 return 0 195} 196 197// ---- k-mer extraction ---- 198 199// Extract a k-mer at base position pos as a packed i64. k must be 200// 1..32. Base 0 of the k-mer sits in the highest 2 bits of the 201// 2*k-bit window: kmer = (b0 << (2k-2)) | (b1 << (2k-4)) | ... | bN. 202// Caller responsibility: pos + k <= len. 203func dna_kmer_at(bases: *u8, pos: i64, k: i64) -> i64 { 204 var kmer: i64 = 0 205 var i: i64 = 0 206 var shift: i64 = (k - 1) << 1 207 while i < k { 208 let b: i64 = dna_get_base(bases, pos + i) 209 kmer = kmer | ((b & 3) << shift) 210 shift = shift - 2 211 i = i + 1 212 } 213 return kmer 214} 215 216// Reverse-complement a packed k-mer in-register. Reverses the 217// base order AND complements each 2-bit code. Used by the canonical 218// k-mer function and by minimizer / hash flows. 219func dna_kmer_revcomp(kmer: i64, k: i64) -> i64 { 220 var rc: i64 = 0 221 var i: i64 = 0 222 var fwd_shift: i64 = (k - 1) << 1 223 var rc_shift: i64 = 0 224 while i < k { 225 let fwd_code: i64 = (kmer >> fwd_shift) & 3 226 let comp: i64 = fwd_code ^ NX_DNA_COMPLEMENT_MASK 227 rc = rc | (comp << rc_shift) 228 fwd_shift = fwd_shift - 2 229 rc_shift = rc_shift + 2 230 i = i + 1 231 } 232 return rc 233} 234 235// Canonical k-mer: lexicographic min of the k-mer and its reverse 236// complement. Because our packing puts base 0 in the highest 2 bits, 237// lexicographic order matches arithmetic order, so this is min(a, b). 238// Palindromic k-mers (kmer == revcomp(kmer)) are their own canonical. 239func dna_kmer_canonical(kmer: i64, k: i64) -> i64 { 240 let rc: i64 = dna_kmer_revcomp(kmer, k) 241 if kmer < rc { return kmer } 242 return rc 243}