nx_variant.nx source
↩ module page · 326 lines · 12924 B
1// nx_variant.nx -- first SNV caller from pileup base counts.
2//
3// license_tier: INDEPENDENT_REDERIVE
4// genealogy_id: international-research-sources/li-2009-samtools-mpileup + depristo-2011-gatk
5//
6// G2.2 of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. The first variant-
7// calling primitive: given a reference base and a pileup-base-count
8// histogram, classify the position as no-variant / heterozygous /
9// homozygous and identify the alt allele.
10//
11// Decision logic (G2.2 reference impl):
12// total_acgt = counts[A] + counts[C] + counts[G] + counts[T]
13// (excludes N and DEL -- the canonical "informative
14// coverage" used by most callers)
15// if total_acgt < min_coverage: NONE (insufficient evidence)
16// alt_code = argmax over c != ref of counts[c]
17// (leftmost-max wins on tie for determinism)
18// if counts[alt_code] == 0: NONE
19// alt_freq_pct = counts[alt_code] * 100 / total_acgt
20// if alt_freq_pct < min_alt_freq_pct: NONE (below threshold)
21// if alt_freq_pct >= min_hom_freq_pct: HOM
22// else: HET
23//
24// Why this simple rule first:
25// Real SNV callers (HaplotypeCaller, Mutect2, FreeBayes) use
26// genotype likelihoods + Bayesian priors + strand-bias filters +
27// base-quality recalibration + indel-realignment + ML scoring.
28// All require log/exp + nx_phred Q->P conversion + statistical
29// distributions -- queued for G2.2b through G2.5.
30// G2.2 ships the IRREDUCIBLE COUNT-BASED LOGIC that every fancy
31// caller's threshold step ultimately reduces to. Composable;
32// sufficient for high-coverage clean data; honest about its limits.
33//
34// What G2.2 does NOT do (deferred):
35// - Quality-weighted counts (Phred-aware) -- G2.2b
36// - Genotype likelihoods + Bayesian posteriors -- G2.2c
37// - Strand-bias filters (skew between fwd + rev reads) -- G2.2d
38// - Indel calling (insertions / deletions vs reference) -- G2.3
39// - Multi-allelic site reporting (2+ alt alleles) -- G2.4
40// - Phasing (linking variants across reads) -- G2.5
41// - Somatic tumor-vs-normal -- G2.6
42// - Population priors (gnomAD-class) -- G2.7
43//
44// API:
45// snv_call(ref_code, counts,
46// min_coverage, min_alt_freq_pct, min_hom_freq_pct,
47// out_alt_code, out_alt_count) -> i64
48//
49// ref_code : NX_DNA_A/C/G/T (0..3)
50// counts : 6-slot histogram from pileup_aggregate_counts
51// min_coverage : minimum total ACGT count to attempt a call
52// min_alt_freq_pct : minimum alt-allele frequency (0..100) to call
53// min_hom_freq_pct : threshold above which call is HOM not HET
54// out_alt_code : 2-bit code of the most-supported non-ref base
55// out_alt_count : raw count for that alt allele
56// Returns: NX_VARIANT_NONE / NX_VARIANT_HET / NX_VARIANT_HOM (from nx_const)
57// or -1 on bad input (ref_code out of range, etc.)
58//
59// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml)
60// intended_use: "First SNV variant caller; the count-based
61// decision step that every probabilistic caller
62// ultimately reduces to. Composes pileup
63// aggregation into actionable variant calls."
64// sil_target: SIL2
65// asil_target: QM
66// dal_target: DAL C
67// iec_62304_class: B
68// evidence: [no_floating_point, deterministic,
69// bit_equal_reproducible,
70// composes_nx_pileup_aggregate_KAT,
71// composes_nx_const_variant_categories,
72// leftmost_alt_tie_break_stable,
73// documented_simple_count_rule_with_honest_limits,
74// license_tier_INDEPENDENT_REDERIVE]
75// hazard_register: [bug-tape-snv-counts-include-ref-double-counted,
76// bug-tape-snv-alt-freq-divide-by-zero,
77// bug-tape-snv-tie-break-rightmost-drift,
78// bug-tape-snv-N-bases-skewing-freq-calculation,
79// bug-tape-snv-hom-het-threshold-inverted]
80// residual_risk: "Pure count rule with no quality weighting;
81// sequencer-error positions at high coverage
82// will produce false-HET calls at the error
83// rate. Quality-weighted version (G2.2b)
84// and probabilistic genotyping (G2.2c) are
85// the upgrade path."
86// verdict: NOT_YET_EVALUATED
87
88import "nx_syscalls.nx"
89import "nx_const.nx"
90import "nx_pileup.nx"
91
92// Indel variant caller (G2.3c) -- mirrors snv_call's count-based
93// decision rule applied to an insertion-length histogram.
94//
95// Decision:
96// - Find argmax over length L > 0 of ins_hist[L].
97// (Length 0 = "no insertion", excluded from alt search.
98// Leftmost-tie-break: smaller L wins on equal count.)
99// - alt_count = ins_hist[best_len]
100// - alt_freq_pct = alt_count * 100 / total_coverage
101// - Apply min_coverage / min_alt_freq_pct / min_hom_freq_pct
102// thresholds same as snv_call.
103//
104// Args:
105// ins_hist : output of indel_count_insertions_by_length
106// max_len : highest bucket in ins_hist
107// total_coverage : total reads spanning this position
108// (includes ins_hist[0] = no-insertion reads)
109// min_coverage / min_alt_freq_pct / min_hom_freq_pct : thresholds
110// out_event_len : dominant insertion length on a call (>= 1)
111// out_event_count : count of reads supporting that length
112//
113// Returns NX_VARIANT_NONE / NX_VARIANT_HET / NX_VARIANT_HOM, or -1
114// on bad input.
115//
116// What G2.3c does NOT do (deferred):
117// - Deletion-variant calling (parallel primitive; same pattern
118// using DEL counts from pileup_aggregate_counts) -- G2.3d
119// - Insertion ALLELE distinguishing (different inserted bases of
120// same length) -- G2.3e
121// - Quality-weighted indel calling -- G2.3f
122func indel_call(ins_hist: *i64,
123 max_len: i64,
124 total_coverage: i64,
125 min_coverage: i64,
126 min_alt_freq_pct: i64,
127 min_hom_freq_pct: i64,
128 out_event_len: *i64,
129 out_event_count: *i64) -> i64 {
130 if max_len < 0 { return -1 }
131 if total_coverage < 0 { return -1 }
132
133 if total_coverage < min_coverage {
134 out_event_len[0] = 0
135 out_event_count[0] = 0
136 return NX_VARIANT_NONE
137 }
138
139 var best_len: i64 = 0
140 var best_count: i64 = 0
141 var L: i64 = 1
142 while L <= max_len {
143 let n: i64 = ins_hist[L]
144 if n > best_count {
145 best_count = n
146 best_len = L
147 }
148 L = L + 1
149 }
150
151 out_event_len[0] = best_len
152 out_event_count[0] = best_count
153
154 if best_count == 0 { return NX_VARIANT_NONE }
155
156 let freq_pct: i64 = (best_count * 100) / total_coverage
157 if freq_pct < min_alt_freq_pct { return NX_VARIANT_NONE }
158 if freq_pct >= min_hom_freq_pct { return NX_VARIANT_HOM }
159 return NX_VARIANT_HET
160}
161
162// Multi-allelic SNV caller (G2.4c) -- returns ALL non-ref alleles
163// with frequency >= min_alt_freq_pct, in descending count order.
164// Enables VCF emission of multi-allelic sites like "A -> G,T" via
165// vcf_write_var_line with comma-separated ALT bytes.
166//
167// Decision:
168// - Scan all 4 ACGT codes; skip ref_code
169// - For each non-ref base with count > 0:
170// alt_freq_pct = count * 100 / total_acgt
171// if alt_freq_pct >= min_alt_freq_pct: emit
172// - Sort emitted alts by count descending (ties: leftmost-tie-break
173// preserved -- smaller code wins on equal count)
174//
175// Args:
176// ref_code, counts, min_alt_freq_pct : same as snv_call
177// out_alts, out_alt_counts : arrays of up to max_alts entries
178// max_alts : capacity (typically 3 since there are 3 non-ref bases)
179// Returns:
180// number of alts emitted (0..3), or -1 on bad input
181//
182// Below min_coverage: returns 0 (no variants reported). Caller
183// applies min_coverage check separately if needed.
184func snv_call_multi(ref_code: i64,
185 counts: *i64,
186 min_alt_freq_pct: i64,
187 out_alts: *i64,
188 out_alt_counts: *i64,
189 max_alts: i64) -> i64 {
190 if ref_code < 0 { return -1 }
191 if ref_code > 3 { return -1 }
192 if max_alts <= 0 { return -1 }
193
194 let total_acgt: i64 = counts[NX_DNA_A] + counts[NX_DNA_C] +
195 counts[NX_DNA_G] + counts[NX_DNA_T]
196 if total_acgt <= 0 { return 0 }
197
198 // Collect candidates that pass min_alt_freq_pct, in ascending
199 // code order (so leftmost-tie-break holds when we sort by count).
200 let cand_codes: *i64 = sys_mmap(64) as *i64
201 let cand_counts: *i64 = sys_mmap(64) as *i64
202 var n_cand: i64 = 0
203 var c: i64 = 0
204 while c < 4 {
205 if c != ref_code {
206 let n: i64 = counts[c]
207 if n > 0 {
208 let freq_pct: i64 = (n * 100) / total_acgt
209 if freq_pct >= min_alt_freq_pct {
210 cand_codes[n_cand] = c
211 cand_counts[n_cand] = n
212 n_cand = n_cand + 1
213 }
214 }
215 }
216 c = c + 1
217 }
218
219 // Insertion-sort by count descending (stable on ties so leftmost
220 // code wins on equal count).
221 var i: i64 = 1
222 while i < n_cand {
223 let key_code: i64 = cand_codes[i]
224 let key_count: i64 = cand_counts[i]
225 var j: i64 = i - 1
226 var keep_going: i64 = 1
227 while keep_going == 1 {
228 if j < 0 {
229 keep_going = 0
230 } else {
231 if cand_counts[j] < key_count {
232 cand_codes[j + 1] = cand_codes[j]
233 cand_counts[j + 1] = cand_counts[j]
234 j = j - 1
235 } else {
236 keep_going = 0
237 }
238 }
239 }
240 cand_codes[j + 1] = key_code
241 cand_counts[j + 1] = key_count
242 i = i + 1
243 }
244
245 // Emit up to max_alts.
246 var emit: i64 = n_cand
247 if emit > max_alts { emit = max_alts }
248 var k: i64 = 0
249 while k < emit {
250 out_alts[k] = cand_codes[k]
251 out_alt_counts[k] = cand_counts[k]
252 k = k + 1
253 }
254 return emit
255}
256
257// Deletion variant caller (G2.3d) -- identical logic to indel_call
258// but operating on the deletion-length histogram (from a separate
259// aggregator over pileup_deletion_after_ref outputs).
260//
261// Args + return mirror indel_call exactly. The insertion and
262// deletion histogram arrays share the same shape (length-keyed
263// count buckets), so callers can reuse indel_count_insertions_by_length
264// with deletion lengths passed in. This primitive exists as a
265// named entry point for symmetry; logic delegates to indel_call's
266// pattern.
267func del_call(del_hist: *i64,
268 max_len: i64,
269 total_coverage: i64,
270 min_coverage: i64,
271 min_alt_freq_pct: i64,
272 min_hom_freq_pct: i64,
273 out_event_len: *i64,
274 out_event_count: *i64) -> i64 {
275 // Logic is structurally identical to indel_call; no special-case
276 // distinction at the count-rule layer. Distinct primitive name
277 // preserves call-site readability + lets G2.3f quality-weighted
278 // and G2.3e allele-distinguishing variants diverge cleanly later.
279 return indel_call(del_hist, max_len, total_coverage,
280 min_coverage, min_alt_freq_pct, min_hom_freq_pct,
281 out_event_len, out_event_count)
282}
283
284func snv_call(ref_code: i64,
285 counts: *i64,
286 min_coverage: i64,
287 min_alt_freq_pct: i64,
288 min_hom_freq_pct: i64,
289 out_alt_code: *i64,
290 out_alt_count: *i64) -> i64 {
291 if ref_code < 0 { return -1 }
292 if ref_code > 3 { return -1 }
293
294 let total_acgt: i64 = counts[NX_DNA_A] + counts[NX_DNA_C] +
295 counts[NX_DNA_G] + counts[NX_DNA_T]
296 if total_acgt < min_coverage {
297 out_alt_code[0] = -1
298 out_alt_count[0] = 0
299 return NX_VARIANT_NONE
300 }
301
302 // Find leftmost-max alt code (most-supported non-ref base).
303 var best_code: i64 = -1
304 var best_count: i64 = 0
305 var c: i64 = 0
306 while c < 4 {
307 if c != ref_code {
308 let n: i64 = counts[c]
309 if n > best_count {
310 best_count = n
311 best_code = c
312 }
313 }
314 c = c + 1
315 }
316
317 out_alt_code[0] = best_code
318 out_alt_count[0] = best_count
319
320 if best_count == 0 { return NX_VARIANT_NONE }
321
322 let freq_pct: i64 = (best_count * 100) / total_acgt
323 if freq_pct < min_alt_freq_pct { return NX_VARIANT_NONE }
324 if freq_pct >= min_hom_freq_pct { return NX_VARIANT_HOM }
325 return NX_VARIANT_HET
326}