nx_pileup.nx source
↩ module page · 433 lines · 17495 B
1// nx_pileup.nx -- CIGAR walking primitives for variant-calling pileup.
2//
3// license_tier: INDEPENDENT_REDERIVE
4// genealogy_id: international-research-sources/sam-spec-v1.4.6 + li-2009-samtools-pileup
5//
6// G2.1a of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. The math foundation
7// for variant calling: given an aligned read's CIGAR string and its
8// alignment ref-start, compute "what query base contributes to a
9// specific reference position?" Aggregating this across many reads
10// at a single ref position IS the pileup, and the pileup IS the
11// input to germline + somatic variant callers.
12//
13// CIGAR walking model (SAM-spec consumption flags from nx_cigar):
14// For each op in alignment order:
15// - if consumes_ref: advance ref-position counter by op_len
16// - if consumes_query: advance query-position counter by op_len
17// Until ref-position counter reaches the target. The resulting
18// query-position counter is the read offset of the base that
19// aligns to the target ref position (or in_gap=1 if the target
20// falls inside a D / N op).
21//
22// Why this primitive matters:
23// - Variant calling: aggregate per-position base counts -> SNVs
24// - Indel detection: I/D ops at specific ref positions
25// - Coverage depth: count reads spanning a position
26// - Quality recalibration (BQSR): per-position quality stats
27// - Pileup visualisation (IGV / samtools tview)
28//
29// API:
30// cigar_query_pos_at_ref(
31// cigar_codes, cigar_lens, n_groups,
32// r_start, target_r,
33// out_q_pos, out_in_gap
34// ) -> i64
35//
36// r_start : reference start position of this alignment (the
37// SAM POS field minus 1 if interpreting SAM 1-indexed;
38// substrate uses 0-indexed throughout)
39// target_r : the reference position we want the read-base for
40// out_q_pos : query offset (0-indexed into read sequence)
41// written if op covering target_r consumes query;
42// set to -1 if target falls in a D/N (deletion)
43// out_in_gap : 1 if target_r falls inside a D / N op; 0 otherwise
44// Returns: 0 success
45// -1 target_r is outside the alignment region
46// (or n_groups < 0 / target before start)
47//
48// What G2.1a does NOT do (deferred):
49// - Multi-read pileup aggregation (sum of bases at one position) -- G2.1b
50// - Quality-weighted pileup -- G2.1c (needs nx_phred composition)
51// - Indel realignment around homopolymers -- G2.5
52// - Variant calling proper (HaplotypeCaller-class) -- G2.2
53// - Tumor-vs-normal somatic calling (Mutect2-class) -- G2.3
54//
55// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml)
56// intended_use: "CIGAR walking for variant-calling pileup;
57// the math primitive that turns aligned-read
58// + position queries into base contributions"
59// sil_target: SIL2
60// asil_target: QM
61// dal_target: DAL C
62// iec_62304_class: B
63// evidence: [no_floating_point, deterministic,
64// bit_equal_reproducible,
65// composes_nx_cigar_consume_helpers,
66// sam_spec_v1_4_6_canonical_consumption,
67// deletion_in_gap_KAT,
68// insertion_shifts_query_KAT,
69// license_tier_INDEPENDENT_REDERIVE]
70// hazard_register: [bug-tape-pileup-1-indexed-vs-0-indexed-off-by-one,
71// bug-tape-pileup-soft-clip-double-counted,
72// bug-tape-pileup-target-equals-end-position-edge,
73// bug-tape-pileup-insertion-base-claimed-at-wrong-pos]
74// residual_risk: "All coordinates 0-indexed throughout the
75// substrate; SAM-spec is 1-indexed for POS
76// so SAM readers must subtract 1 before
77// passing r_start. Documented in API but
78// a recurring footgun across genomics tools."
79// verdict: NOT_YET_EVALUATED
80
81import "nx_syscalls.nx"
82import "nx_const.nx"
83import "nx_cigar.nx"
84import "nx_sequence.nx"
85import "nx_sequence_ambig.nx"
86
87// Pileup contribution flags from pileup_contribute_one_read.
88const NX_PILEUP_FLAG_BASE: i64 = 0 // normal ACGT base contributed
89const NX_PILEUP_FLAG_N: i64 = 1 // ambiguous base (N) contributed
90const NX_PILEUP_FLAG_DEL: i64 = 2 // position in deletion gap, no base
91const NX_PILEUP_FLAG_OUT: i64 = 3 // target_r outside alignment region
92
93// Walk the CIGAR to find which query position corresponds to a target
94// reference position. See header for return codes + flag semantics.
95func cigar_query_pos_at_ref(cigar_codes: *u8, cigar_lens: *i64, n_groups: i64,
96 r_start: i64, target_r: i64,
97 out_q_pos: *i64, out_in_gap: *i64) -> i64 {
98 if n_groups < 0 { return -1 }
99 if target_r < r_start { return -1 }
100
101 var cur_r: i64 = r_start
102 var cur_q: i64 = 0
103
104 var g: i64 = 0
105 while g < n_groups {
106 let op: i64 = cigar_codes[g] & 0xff
107 let len: i64 = cigar_lens[g]
108 let cq: i64 = nx_cigar_consumes_query(op)
109 let cr: i64 = nx_cigar_consumes_ref(op)
110
111 if cr == 1 {
112 // This op spans ref [cur_r, cur_r + len). Does it cover target?
113 if target_r < cur_r + len {
114 let r_offset_in_op: i64 = target_r - cur_r
115 if cq == 1 {
116 // M / = / X -- target is matched / mismatched; emit q-pos.
117 out_q_pos[0] = cur_q + r_offset_in_op
118 out_in_gap[0] = 0
119 } else {
120 // D / N -- target is in a deletion / skip; no read base.
121 out_q_pos[0] = -1
122 out_in_gap[0] = 1
123 }
124 return 0
125 }
126 cur_r = cur_r + len
127 }
128
129 if cq == 1 {
130 // M / I / S / = / X consume query. Already handled match
131 // case above; for I / S we advance query without touching ref.
132 cur_q = cur_q + len
133 }
134
135 g = g + 1
136 }
137
138 // Target is past the end of the alignment.
139 return -1
140}
141
142// Single-read pileup contribution. Composes cigar_query_pos_at_ref
143// + nx_sequence base accessors + N sidecar check into one primitive
144// suitable for accumulation across many reads at a target ref pos.
145//
146// out_base_code: when flag is BASE or N, holds the 2-bit DNA code (0..3)
147// (for N positions the code is whatever bases[] stored,
148// typically 0 = A as a placeholder per nx_sequence_ambig
149// convention; caller should rely on out_flags not the
150// code value when flag==N)
151// out_flags : NX_PILEUP_FLAG_BASE / _N / _DEL / _OUT
152// Returns: 0 contributed (BASE or N or DEL)
153// -1 outside alignment region (caller skips this read)
154func pileup_contribute_one_read(
155 cigar_codes: *u8, cigar_lens: *i64, cigar_n: i64,
156 r_start: i64, target_r: i64,
157 read_bases: *u8, read_nbits: *u8,
158 out_base_code: *i64, out_flags: *i64) -> i64 {
159
160 let q_pos_buf: *i64 = sys_mmap(16) as *i64
161 let in_gap_buf: *i64 = sys_mmap(16) as *i64
162
163 let rc: i64 = cigar_query_pos_at_ref(cigar_codes, cigar_lens, cigar_n,
164 r_start, target_r,
165 q_pos_buf, in_gap_buf)
166 if rc < 0 {
167 out_base_code[0] = -1
168 out_flags[0] = NX_PILEUP_FLAG_OUT
169 return -1
170 }
171
172 if in_gap_buf[0] == 1 {
173 out_base_code[0] = -1
174 out_flags[0] = NX_PILEUP_FLAG_DEL
175 return 0
176 }
177
178 let q_pos: i64 = q_pos_buf[0]
179 let code: i64 = dna_get_base(read_bases, q_pos)
180 out_base_code[0] = code
181
182 if dna_is_n(read_nbits, q_pos) == 1 {
183 out_flags[0] = NX_PILEUP_FLAG_N
184 } else {
185 out_flags[0] = NX_PILEUP_FLAG_BASE
186 }
187 return 0
188}
189
190// Indices into the 6-slot counts array produced by
191// pileup_aggregate_counts. First 4 align with NX_DNA_A/C/G/T so
192// callers can index by base-code directly.
193const NX_PILEUP_COUNT_N: i64 = 4
194const NX_PILEUP_COUNT_DEL: i64 = 5
195const NX_PILEUP_COUNT_LEN: i64 = 6
196
197// Detect an INSERTION event immediately after target_r in a read's CIGAR.
198//
199// VCF convention: an insertion is reported at the reference position
200// IMMEDIATELY BEFORE the inserted bases (i.e., the last reference-
201// consuming base before the I op). This primitive walks the CIGAR
202// and returns the insertion event at target_r if one exists.
203//
204// Deletions are already detected by pileup_contribute_one_read via
205// NX_PILEUP_FLAG_DEL; this primitive completes the indel-detection
206// pair.
207//
208// Args:
209// cigar_codes, cigar_lens, n_groups : compressed CIGAR from nx_cigar
210// r_start, target_r : alignment ref-start + target ref position (0-indexed)
211// out_ins_length : bp count of the inserted run (>=1 on hit)
212// out_q_pos : query offset where the inserted bases start
213// Returns:
214// 0 = insertion event found at target_r (out_* written)
215// 1 = no insertion at target_r (target_r is in alignment but not pre-I)
216// -1 = target_r outside alignment region
217func pileup_insertion_after_ref(cigar_codes: *u8, cigar_lens: *i64, n_groups: i64,
218 r_start: i64, target_r: i64,
219 out_ins_length: *i64,
220 out_q_pos: *i64) -> i64 {
221 if n_groups < 0 { return -1 }
222 if target_r < r_start { return -1 }
223
224 var cur_r: i64 = r_start
225 var cur_q: i64 = 0
226
227 var g: i64 = 0
228 while g < n_groups {
229 let op: i64 = cigar_codes[g] & 0xff
230 let len: i64 = cigar_lens[g]
231 let cq: i64 = nx_cigar_consumes_query(op)
232 let cr: i64 = nx_cigar_consumes_ref(op)
233
234 if cr == 1 {
235 if target_r < cur_r + len {
236 // target_r lives inside this op's ref span.
237 if cq == 1 {
238 // M / = / X consume both -- candidate for "last
239 // base before an I op".
240 if target_r == cur_r + len - 1 {
241 // Check next op for I.
242 let next_g: i64 = g + 1
243 if next_g < n_groups {
244 let next_op: i64 = cigar_codes[next_g] & 0xff
245 if next_op == NX_CIGAR_I {
246 out_ins_length[0] = cigar_lens[next_g]
247 out_q_pos[0] = cur_q + (target_r - cur_r) + 1
248 return 0
249 }
250 }
251 }
252 }
253 // Not an insertion at this position.
254 out_ins_length[0] = 0
255 out_q_pos[0] = -1
256 return 1
257 }
258 cur_r = cur_r + len
259 }
260 if cq == 1 {
261 cur_q = cur_q + len
262 }
263 g = g + 1
264 }
265 return -1
266}
267
268// Encode an insertion event's bases + length into a single i64 hash
269// suitable for allele-aggregation (G2.3e). Same inserted-base sequence
270// + same length -> same hash; any difference -> different hash.
271//
272// Layout:
273// bits 56..63 = length (8 bits, max 255 -- realistic indels <= 50bp)
274// bits 0..55 = 2-bit packed bases, base 0 in highest position
275// within the field (consistent with nx_sequence k-mer
276// packing convention). Capacity: up to 28 bases.
277//
278// Returns the encoded i64 hash, or -1 for invalid input (ins_len <= 0
279// or > 28, or q_pos < 0). Insertions longer than 28bp use the longer-
280// allele primitive (G2.3e.2, deferred).
281//
282// To decode: high 8 bits = length, low bits = packed bases.
283// Caller uses dna_get_base-style extraction at bit positions
284// 2*(L-1-i) .. 2*(L-1-i)+1 for base i.
285func ins_event_to_hash(read_bases: *u8, q_pos: i64, ins_len: i64) -> i64 {
286 if q_pos < 0 { return -1 }
287 if ins_len <= 0 { return -1 }
288 if ins_len > 28 { return -1 }
289
290 var packed: i64 = 0
291 var i: i64 = 0
292 while i < ins_len {
293 let code: i64 = dna_get_base(read_bases, q_pos + i)
294 let bit_pos: i64 = (ins_len - 1 - i) * 2
295 packed = packed | ((code & 3) << bit_pos)
296 i = i + 1
297 }
298 return (ins_len << 56) | packed
299}
300
301// Decode an insertion hash back into length + packed bases.
302// out_packed receives the 2-bit packed bases (as in_event_to_hash).
303// Returns the length, or -1 if hash has length out of range (1..28).
304func ins_event_decode_hash(hash: i64, out_packed: *i64) -> i64 {
305 let len: i64 = (hash >> 56) & 0xff
306 if len < 1 { return -1 }
307 if len > 28 { return -1 }
308 // Mask off the length field; keep only the bases bits.
309 let mask: i64 = (1 << 56) - 1
310 out_packed[0] = hash & mask
311 return len
312}
313
314// Detect a DELETION event immediately after target_r in a read's CIGAR.
315// Symmetric to pileup_insertion_after_ref: deletion is reported at
316// the reference position IMMEDIATELY BEFORE the deleted bases (the
317// last reference-consuming M/=/X base before the D op), matching the
318// VCF convention.
319//
320// For example, CIGAR "3M3D2M" starting at ref 100:
321// 3M consumes ref 100,101,102 (target_r=102 is the last M)
322// 3D consumes ref 103,104,105 (no query advance)
323// 2M consumes ref 106,107
324//
325// Returns 0 if a deletion event starts after target_r (with
326// out_del_length written), 1 if no deletion event at target_r,
327// -1 if target_r is outside the alignment region.
328func pileup_deletion_after_ref(cigar_codes: *u8, cigar_lens: *i64, n_groups: i64,
329 r_start: i64, target_r: i64,
330 out_del_length: *i64) -> i64 {
331 if n_groups < 0 { return -1 }
332 if target_r < r_start { return -1 }
333
334 var cur_r: i64 = r_start
335
336 var g: i64 = 0
337 while g < n_groups {
338 let op: i64 = cigar_codes[g] & 0xff
339 let len: i64 = cigar_lens[g]
340 let cq: i64 = nx_cigar_consumes_query(op)
341 let cr: i64 = nx_cigar_consumes_ref(op)
342
343 if cr == 1 {
344 if target_r < cur_r + len {
345 if cq == 1 {
346 if target_r == cur_r + len - 1 {
347 // Last base of M/=/X op. Check next op for D.
348 let next_g: i64 = g + 1
349 if next_g < n_groups {
350 let next_op: i64 = cigar_codes[next_g] & 0xff
351 if next_op == NX_CIGAR_D {
352 out_del_length[0] = cigar_lens[next_g]
353 return 0
354 }
355 }
356 }
357 }
358 out_del_length[0] = 0
359 return 1
360 }
361 cur_r = cur_r + len
362 }
363 g = g + 1
364 }
365 return -1
366}
367
368// Aggregate per-read insertion lengths into a length histogram.
369//
370// Per-read input: ins_lens[i] is the insertion length detected at
371// the position for read i. 0 means no insertion at this position
372// (either rc==1 from pileup_insertion_after_ref OR read outside).
373// Length L > max_len gets bucketed into out_counts[max_len].
374//
375// Args:
376// ins_lens, n_reads : per-read insertion-length array
377// out_counts : pre-zeroed array of size max_len+1
378// max_len : highest length bucket (lengths > max_len
379// collapse into out_counts[max_len])
380// Returns:
381// total count of reads with ins_len > 0 (i.e., reads with any
382// insertion at this position). -1 on bad input.
383func indel_count_insertions_by_length(ins_lens: *i64, n_reads: i64,
384 out_counts: *i64, max_len: i64) -> i64 {
385 if n_reads < 0 { return -1 }
386 if max_len < 0 { return -1 }
387
388 var with_ins: i64 = 0
389 var i: i64 = 0
390 while i < n_reads {
391 let len: i64 = ins_lens[i]
392 if len < 0 { return -1 } // defensive; ins_len should never be negative
393 var bucket: i64 = len
394 if bucket > max_len { bucket = max_len }
395 out_counts[bucket] = out_counts[bucket] + 1
396 if len > 0 { with_ins = with_ins + 1 }
397 i = i + 1
398 }
399 return with_ins
400}
401
402// Aggregate per-read pileup contributions into a 6-slot count
403// histogram [A, C, G, T, N, DEL]. out_counts must have capacity
404// for 6 i64s and be pre-zeroed (sys_mmap pages are zeroed).
405// Reads with flag == OUT are skipped (they did not align here).
406// Returns the total number of contributing reads (BASE + N + DEL).
407func pileup_aggregate_counts(base_codes: *i64, flags: *i64, n_reads: i64,
408 out_counts: *i64) -> i64 {
409 if n_reads < 0 { return -1 }
410 var total: i64 = 0
411 var i: i64 = 0
412 while i < n_reads {
413 let f: i64 = flags[i]
414 if f == NX_PILEUP_FLAG_BASE {
415 let c: i64 = base_codes[i] & 3
416 out_counts[c] = out_counts[c] + 1
417 total = total + 1
418 } else {
419 if f == NX_PILEUP_FLAG_N {
420 out_counts[NX_PILEUP_COUNT_N] = out_counts[NX_PILEUP_COUNT_N] + 1
421 total = total + 1
422 } else {
423 if f == NX_PILEUP_FLAG_DEL {
424 out_counts[NX_PILEUP_COUNT_DEL] = out_counts[NX_PILEUP_COUNT_DEL] + 1
425 total = total + 1
426 }
427 // OUT flag is skipped silently.
428 }
429 }
430 i = i + 1
431 }
432 return total
433}