nx_phred.nx source
↩ module page · 138 lines · 5796 B
1// nx_phred.nx -- Phred quality score encode / decode helpers.
2//
3// license_tier: INDEPENDENT_REDERIVE
4// genealogy_id: international-research-sources/sam-spec-v1.4.6 + ewing-green-1998-phred
5//
6// G2.0a of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. Bridges to
7// FASTQ / SAM real-world data: every sequencer-emitted base carries
8// a Phred quality score Q = -10 * log10(P_error), encoded in FASTQ
9// as ASCII byte (Q + offset). This primitive ships ONLY the
10// encode/decode layer; probability conversions (Q -> P_error in
11// fixed-point) require log/exp primitives queued for G2.0b.
12//
13// Two encodings:
14// Phred+33 (Sanger / Illumina 1.8+): the modern universal default
15// - ASCII = Q + 33
16// - Q range [0..93] maps to ASCII [33..126] = '!' to '~'
17// - Used by ALL contemporary sequencers + SAM/BAM
18// Phred+64 (legacy Illumina 1.3 / 1.5 / 1.7): historical only
19// - ASCII = Q + 64
20// - Q range [0..62] maps to ASCII [64..126] = '@' to '~'
21// - Substrate supports for archived data + cross-compatibility
22//
23// Universal Q-value benchmarks:
24// Q10 = 90% accuracy (10% error)
25// Q20 = 99% accuracy (1% error) <- variant-calling minimum threshold
26// Q30 = 99.9% accuracy (0.1% error) <- typical Illumina target
27// Q40 = 99.99% accuracy <- typical PacBio HiFi
28// Q60 = 99.9999% accuracy <- consensus / theoretical max
29//
30// API:
31// nx_phred33_encode(q) -> i64 Q [0..93] -> ASCII byte; -1 out of range
32// nx_phred33_decode(ascii) -> i64 ASCII [33..126] -> Q; -1 out of range
33// nx_phred64_encode(q) -> i64 Q [0..62] -> ASCII byte; -1 out of range
34// nx_phred64_decode(ascii) -> i64 ASCII [64..126] -> Q; -1 out of range
35// nx_phred33_encode_string(qs, n, out_ascii) -> i64 bulk encode
36// nx_phred33_decode_string(ascii, n, out_qs) -> i64 bulk decode
37//
38// What G2.0a does NOT do (deferred):
39// - Q -> P_error fixed-point conversion -- needs nx_fixed_point + log10 (G2.0b)
40// - Per-position quality aggregation (consensus quality scoring) -- G2.0c
41// - Quality recalibration (BQSR-class) -- G2.5
42// - Solexa-encoded ancient data (pre-2009) -- not in scope
43//
44// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml)
45// intended_use: "FASTQ + SAM Phred quality encoding/decoding;
46// the bridge between sequencer-emitted ASCII
47// and Q-value arithmetic"
48// sil_target: SIL2
49// asil_target: QM
50// dal_target: DAL C
51// iec_62304_class: B
52// evidence: [no_floating_point, deterministic,
53// bit_equal_reproducible,
54// sam_spec_v1_4_6_phred33_canonical,
55// composes_nx_const_phred_constants,
56// round_trip_q_to_ascii_to_q_KAT,
57// license_tier_INDEPENDENT_REDERIVE]
58// hazard_register: [bug-tape-phred33-vs-phred64-misdetect,
59// bug-tape-phred-out-of-range-silent-wraparound,
60// bug-tape-bulk-encode-decode-length-mismatch,
61// bug-tape-q-value-cap-93-not-enforced]
62// residual_risk: "G2.0a does not auto-detect encoding;
63// callers must know which offset their data
64// uses. Detection heuristics (presence of
65// ASCII < 64) live in G2.0a.1 + are noisy
66// for low-Q legacy data. Variant-quality
67// ARITHMETIC requires G2.0b log primitives."
68// verdict: NOT_YET_EVALUATED
69
70import "nx_syscalls.nx"
71import "nx_const.nx"
72
73// Phred+33 single-value encode.
74// Q [0..93] -> ASCII [33..126]. Returns -1 if Q out of range.
75func nx_phred33_encode(q: i64) -> i64 {
76 if q < NX_PHRED_MIN_Q { return -1 }
77 if q > NX_PHRED_MAX_Q { return -1 }
78 return q + NX_PHRED33_OFFSET
79}
80
81// Phred+33 single-value decode.
82// ASCII [33..126] -> Q [0..93]. Returns -1 if ASCII out of range.
83func nx_phred33_decode(ascii: i64) -> i64 {
84 let a: i64 = ascii & 0xff
85 if a < NX_PHRED33_OFFSET { return -1 }
86 let q: i64 = a - NX_PHRED33_OFFSET
87 if q > NX_PHRED_MAX_Q { return -1 }
88 return q
89}
90
91// Phred+64 single-value encode.
92// Q [0..62] -> ASCII [64..126]. Returns -1 if Q out of range.
93// Cap at 62 ensures encoded ASCII does not exceed 126 ('~').
94func nx_phred64_encode(q: i64) -> i64 {
95 if q < NX_PHRED_MIN_Q { return -1 }
96 if q > 62 { return -1 }
97 return q + NX_PHRED64_OFFSET
98}
99
100// Phred+64 single-value decode.
101// ASCII [64..126] -> Q [0..62]. Returns -1 if ASCII out of range.
102func nx_phred64_decode(ascii: i64) -> i64 {
103 let a: i64 = ascii & 0xff
104 if a < NX_PHRED64_OFFSET { return -1 }
105 let q: i64 = a - NX_PHRED64_OFFSET
106 if q > 62 { return -1 }
107 return q
108}
109
110// Bulk Phred+33 encode. qs[i] is the Q value for position i;
111// out_ascii[i] gets the encoded byte. Returns 0 on success or
112// -1 if any q is out of range (and stops at the bad index).
113func nx_phred33_encode_string(qs: *i64, n: i64, out_ascii: *u8) -> i64 {
114 if n < 0 { return -1 }
115 var i: i64 = 0
116 while i < n {
117 let a: i64 = nx_phred33_encode(qs[i])
118 if a < 0 { return -1 }
119 out_ascii[i] = a & 0xff
120 i = i + 1
121 }
122 return 0
123}
124
125// Bulk Phred+33 decode. ascii[i] is the FASTQ byte; out_qs[i] gets
126// the decoded Q. Returns 0 on success or -1 if any byte is out of
127// range (and stops at the bad index).
128func nx_phred33_decode_string(ascii: *u8, n: i64, out_qs: *i64) -> i64 {
129 if n < 0 { return -1 }
130 var i: i64 = 0
131 while i < n {
132 let q: i64 = nx_phred33_decode(ascii[i] & 0xff)
133 if q < 0 { return -1 }
134 out_qs[i] = q
135 i = i + 1
136 }
137 return 0
138}