nx_align_backtrace.nx source
↩ module page · 474 lines · 17340 B
1// nx_align_backtrace.nx -- Smith-Waterman backtrace with CIGAR output.
2//
3// license_tier: INDEPENDENT_REDERIVE
4// genealogy_id: international-research-sources/smith-waterman-1981 + sam-spec-v1.4.6
5//
6// G1.0c of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. Walks the SW DP
7// table backward from the max-score cell to recover the operation
8// string of the optimal local alignment. Output uses SAM/BAM
9// CIGAR codes (NX_CIGAR_M / _I / _D from nx_const).
10//
11// Why a separate primitive rather than extending smith_waterman_linear:
12// The existing G1.0 SW returns score + end-coordinates and discards
13// the DP table internally. Backtrace needs the full table. Rather
14// than complicate the G1.0 API with optional-output flags, this
15// primitive owns its own table allocation + DP fill + backtrace.
16// Some DP-code duplication is the price for clean API separation.
17// (G1.5 perf path will share a sw_dp helper used by both.)
18//
19// Tie-break for predecessor selection:
20// When multiple predecessors produce the current cell's value
21// (rare but possible), the order checked is: diagonal -> up -> left.
22// This favours matches/mismatches over gaps, matching the
23// minimap2 + BWA-MEM convention. Deterministic across hosts.
24//
25// Stop condition:
26// Backtrace stops when H[i][j] == 0 -- the local-alignment 0-floor
27// marks the start of the optimal region. Out-coordinates record
28// where the alignment STARTED (i_start, j_start), so the alignment
29// spans a[i_start .. i_max) vs b[j_start .. j_max) using ops as
30// the operation stream from left to right.
31//
32// What G1.0c does NOT do (deferred):
33// - Affine-gap backtrace (needs separate E/F predecessor matrices) -- G1.0d
34// - Run-length compression of ops (e.g., "MMMM" -> "4M") -- G1.4c
35// - CIGAR string serialisation as ASCII -- nx_sam_writer (G6)
36// - Soft-clipping of unaligned query ends -- G1.4c
37// - Explicit = / X distinction (we use M for both) -- G1.4c
38//
39// API:
40// smith_waterman_linear_backtrace(
41// a, n, b, m,
42// match_s, mis_s, gap_s,
43// out_ops, max_ops, out_meta
44// ) -> i64 op count (>= 0); -1 on overflow / bad input
45//
46// out_ops : byte array of NX_CIGAR_* codes in alignment order
47// (left-to-right; out_ops[0] is the op at the start
48// of the alignment region)
49// out_meta : writes 5 i64s [score, max_i, max_j, start_i, start_j]
50//
51// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml)
52// intended_use: "SW backtrace -- recovers the alignment
53// operation stream for CIGAR construction;
54// prerequisite for SAM/BAM output and for
55// variant calling (pileup needs op stream)"
56// sil_target: SIL2
57// asil_target: QM
58// dal_target: DAL C
59// iec_62304_class: B
60// evidence: [no_floating_point, deterministic,
61// bit_equal_reproducible,
62// diag_up_left_tie_break_stable,
63// zero_floor_stop_textbook,
64// composes_nx_const_NX_CIGAR_codes,
65// gap_path_hand_traced_KAT,
66// license_tier_INDEPENDENT_REDERIVE]
67// hazard_register: [bug-tape-backtrace-no-predecessor-found-infinite-loop,
68// bug-tape-backtrace-ops-not-reversed,
69// bug-tape-backtrace-start-coords-wrong-by-one,
70// bug-tape-backtrace-i-or-j-off-end-overrun]
71// residual_risk: "Pred selection on ties (diag > up > left)
72// means alignments with equal-score alternative
73// paths return ONE representative; downstream
74// tools that depend on a specific tie-break
75// must agree with this convention or wrap."
76// verdict: NOT_YET_EVALUATED
77
78import "nx_syscalls.nx"
79import "nx_const.nx"
80import "nx_align.nx"
81import "nx_align_affine.nx" // for NX_NEG_INF_I64 sentinel
82
83// State-machine backtrace states (which Gotoh matrix the path is in).
84const NX_AFFINE_STATE_H: i64 = 0
85const NX_AFFINE_STATE_E: i64 = 1
86const NX_AFFINE_STATE_F: i64 = 2
87
88func smith_waterman_linear_backtrace(a: *u8, n: i64,
89 b: *u8, m: i64,
90 match_s: i64,
91 mis_s: i64,
92 gap_s: i64,
93 out_ops: *u8,
94 max_ops: i64,
95 out_meta: *i64) -> i64 {
96 if n <= 0 {
97 out_meta[0] = 0
98 out_meta[1] = 0
99 out_meta[2] = 0
100 out_meta[3] = 0
101 out_meta[4] = 0
102 return 0
103 }
104 if m <= 0 {
105 out_meta[0] = 0
106 out_meta[1] = 0
107 out_meta[2] = 0
108 out_meta[3] = 0
109 out_meta[4] = 0
110 return 0
111 }
112 if max_ops <= 0 { return -1 }
113
114 let cols: i64 = m + 1
115 let rows: i64 = n + 1
116 let bytes: i64 = rows * cols * 8
117 let h: *i64 = sys_mmap(bytes) as *i64
118
119 // ---- DP fill (mirrors smith_waterman_linear from nx_align.nx) ----
120 var best_score: i64 = 0
121 var best_i: i64 = 0
122 var best_j: i64 = 0
123
124 var i: i64 = 1
125 while i <= n {
126 let a_byte: i64 = a[i - 1] & 0xff
127 let row_off: i64 = i * cols
128 let prev_row_off: i64 = (i - 1) * cols
129 var j: i64 = 1
130 while j <= m {
131 let b_byte: i64 = b[j - 1] & 0xff
132
133 var sub_pen: i64 = mis_s
134 if a_byte == b_byte { sub_pen = match_s }
135 let s_diag: i64 = h[prev_row_off + j - 1] + sub_pen
136 let s_up: i64 = h[prev_row_off + j] + gap_s
137 let s_left: i64 = h[row_off + j - 1] + gap_s
138
139 var cell: i64 = s_diag
140 if s_up > cell { cell = s_up }
141 if s_left > cell { cell = s_left }
142 if cell < 0 { cell = 0 }
143 h[row_off + j] = cell
144
145 if cell > best_score {
146 best_score = cell
147 best_i = i
148 best_j = j
149 }
150 j = j + 1
151 }
152 i = i + 1
153 }
154
155 out_meta[0] = best_score
156 out_meta[1] = best_i
157 out_meta[2] = best_j
158
159 // ---- Backtrace ----
160 // Emit ops into a reverse buffer, then reverse into out_ops.
161 let rev_ops: *u8 = sys_mmap(max_ops + 8)
162 var op_count: i64 = 0
163 var ci: i64 = best_i
164 var cj: i64 = best_j
165
166 while ci > 0 {
167 if cj <= 0 { ci = 0 }
168 else {
169 let cur_off: i64 = ci * cols
170 let h_cur: i64 = h[cur_off + cj]
171 if h_cur == 0 {
172 ci = 0 // signal stop
173 } else {
174 let a_byte: i64 = a[ci - 1] & 0xff
175 let b_byte: i64 = b[cj - 1] & 0xff
176 var sub_pen: i64 = mis_s
177 if a_byte == b_byte { sub_pen = match_s }
178
179 let prev_row_off: i64 = (ci - 1) * cols
180 let v_diag: i64 = h[prev_row_off + cj - 1] + sub_pen
181 let v_up: i64 = h[prev_row_off + cj] + gap_s
182 let v_left: i64 = h[cur_off + cj - 1] + gap_s
183
184 if op_count >= max_ops { return -1 }
185
186 // Diag -> up -> left tie-break.
187 if h_cur == v_diag {
188 rev_ops[op_count] = NX_CIGAR_M & 0xff
189 op_count = op_count + 1
190 ci = ci - 1
191 cj = cj - 1
192 } else {
193 if h_cur == v_up {
194 rev_ops[op_count] = NX_CIGAR_I & 0xff
195 op_count = op_count + 1
196 ci = ci - 1
197 } else {
198 if h_cur == v_left {
199 rev_ops[op_count] = NX_CIGAR_D & 0xff
200 op_count = op_count + 1
201 cj = cj - 1
202 } else {
203 // Defensive: no predecessor matches -- stop.
204 ci = 0
205 }
206 }
207 }
208 }
209 }
210 }
211
212 out_meta[3] = ci
213 out_meta[4] = cj
214
215 // Reverse rev_ops into out_ops to get alignment-order stream.
216 var k: i64 = 0
217 while k < op_count {
218 out_ops[k] = rev_ops[op_count - 1 - k]
219 k = k + 1
220 }
221
222 return op_count
223}
224
225// G1.0d -- Affine-gap Smith-Waterman backtrace.
226//
227// Composes Gotoh 1982 3-matrix DP (H / E / F per nx_align_affine)
228// with a state-machine backtrace that tracks which matrix the path
229// is currently in:
230//
231// STATE_H: in H matrix. At cell (i,j), check H[i][j] == 0 (stop),
232// else determine if path came via diagonal (M), via E
233// (switch state), or via F (switch state).
234//
235// STATE_E: in E matrix (gap in a -- left walk, emit D per step).
236// Determine if the previous E cell was opened from H
237// (switch back to STATE_H next iteration) or extended
238// from another E (stay).
239//
240// STATE_F: in F matrix (gap in b -- up walk, emit I per step).
241// Symmetric to STATE_E.
242//
243// Tie-break preferences:
244// - In STATE_H: diag > E > F (favour matches over gap entries)
245// - In STATE_E/F: extend > open (favour staying in the gap so the
246// walk completes the entire extension
247// rather than splitting it).
248//
249// API:
250// smith_waterman_affine_backtrace(
251// a, n, b, m,
252// fields, // [match_s, mismatch_s, gap_open_s, gap_extend_s]
253// out_ops, max_ops,
254// out_meta // [score, max_i, max_j, start_i, start_j]
255// ) -> i64 op_count or -1
256//
257// Memory: 3*(n+1)*(m+1)*8 bytes for H/E/F. Same as smith_waterman_affine.
258func smith_waterman_affine_backtrace(a: *u8, n: i64,
259 b: *u8, m: i64,
260 fields: *i64,
261 out_ops: *u8, max_ops: i64,
262 out_meta: *i64) -> i64 {
263 let match_s: i64 = fields[0]
264 let mis_s: i64 = fields[1]
265 let gap_open: i64 = fields[2]
266 let gap_extend: i64 = fields[3]
267
268 if n <= 0 {
269 out_meta[0] = 0; out_meta[1] = 0; out_meta[2] = 0
270 out_meta[3] = 0; out_meta[4] = 0
271 return 0
272 }
273 if m <= 0 {
274 out_meta[0] = 0; out_meta[1] = 0; out_meta[2] = 0
275 out_meta[3] = 0; out_meta[4] = 0
276 return 0
277 }
278 if max_ops <= 0 { return -1 }
279
280 let cols: i64 = m + 1
281 let rows: i64 = n + 1
282 let bytes: i64 = rows * cols * 8
283
284 let h: *i64 = sys_mmap(bytes) as *i64
285 let e: *i64 = sys_mmap(bytes) as *i64
286 let f: *i64 = sys_mmap(bytes) as *i64
287
288 // sys_mmap returns zeroed; H boundary already 0.
289 // Set E/F boundaries to NEG_INF so "open a gap into nothing" can't
290 // contribute a real score.
291 var j_init: i64 = 0
292 while j_init <= m {
293 f[j_init] = NX_NEG_INF_I64
294 j_init = j_init + 1
295 }
296 var i_init: i64 = 0
297 while i_init <= n {
298 e[i_init * cols] = NX_NEG_INF_I64
299 i_init = i_init + 1
300 }
301
302 let gap_first: i64 = gap_open + gap_extend
303
304 // ---- Forward DP fill (Gotoh) ----
305 var best_score: i64 = 0
306 var best_i: i64 = 0
307 var best_j: i64 = 0
308
309 var i: i64 = 1
310 while i <= n {
311 let a_byte: i64 = a[i - 1] & 0xff
312 let row_off: i64 = i * cols
313 let prev_row_off: i64 = (i - 1) * cols
314
315 var j: i64 = 1
316 while j <= m {
317 let b_byte: i64 = b[j - 1] & 0xff
318
319 let e_open: i64 = h[row_off + j - 1] + gap_first
320 let e_ext: i64 = e[row_off + j - 1] + gap_extend
321 var e_cell: i64 = e_open
322 if e_ext > e_cell { e_cell = e_ext }
323 e[row_off + j] = e_cell
324
325 let f_open: i64 = h[prev_row_off + j] + gap_first
326 let f_ext: i64 = f[prev_row_off + j] + gap_extend
327 var f_cell: i64 = f_open
328 if f_ext > f_cell { f_cell = f_ext }
329 f[row_off + j] = f_cell
330
331 var sub_pen: i64 = mis_s
332 if a_byte == b_byte { sub_pen = match_s }
333 let h_diag: i64 = h[prev_row_off + j - 1] + sub_pen
334
335 var h_cell: i64 = 0
336 if h_diag > h_cell { h_cell = h_diag }
337 if e_cell > h_cell { h_cell = e_cell }
338 if f_cell > h_cell { h_cell = f_cell }
339 h[row_off + j] = h_cell
340
341 if h_cell > best_score {
342 best_score = h_cell
343 best_i = i
344 best_j = j
345 }
346 j = j + 1
347 }
348 i = i + 1
349 }
350
351 out_meta[0] = best_score
352 out_meta[1] = best_i
353 out_meta[2] = best_j
354
355 // ---- Backtrace state machine ----
356 let rev_ops: *u8 = sys_mmap(max_ops + 8)
357 var op_count: i64 = 0
358 var ci: i64 = best_i
359 var cj: i64 = best_j
360 var state: i64 = NX_AFFINE_STATE_H
361
362 var loop_active: i64 = 1
363 while loop_active == 1 {
364 if state == NX_AFFINE_STATE_H {
365 if ci <= 0 { loop_active = 0 }
366 else { if cj <= 0 { loop_active = 0 }
367 else {
368 let cur_off: i64 = ci * cols
369 let h_cur: i64 = h[cur_off + cj]
370 if h_cur == 0 {
371 loop_active = 0
372 } else {
373 let a_byte: i64 = a[ci - 1] & 0xff
374 let b_byte: i64 = b[cj - 1] & 0xff
375 var sub_pen: i64 = mis_s
376 if a_byte == b_byte { sub_pen = match_s }
377 let prev_row_off: i64 = (ci - 1) * cols
378 let h_diag: i64 = h[prev_row_off + cj - 1] + sub_pen
379 let e_cur: i64 = e[cur_off + cj]
380 let f_cur: i64 = f[cur_off + cj]
381
382 if op_count >= max_ops { return -1 }
383
384 // Tie-break: diag > E > F.
385 if h_cur == h_diag {
386 rev_ops[op_count] = NX_CIGAR_M & 0xff
387 op_count = op_count + 1
388 ci = ci - 1
389 cj = cj - 1
390 // stay STATE_H
391 } else {
392 if h_cur == e_cur {
393 // switch to STATE_E -- no op emitted on switch
394 state = NX_AFFINE_STATE_E
395 } else {
396 if h_cur == f_cur {
397 state = NX_AFFINE_STATE_F
398 } else {
399 loop_active = 0 // defensive: no path
400 }
401 }
402 }
403 }
404 } }
405 } else {
406 if state == NX_AFFINE_STATE_E {
407 if cj <= 0 { loop_active = 0 }
408 else {
409 let cur_off: i64 = ci * cols
410 let e_cur: i64 = e[cur_off + cj]
411 let h_prev: i64 = h[cur_off + cj - 1]
412 let e_prev: i64 = e[cur_off + cj - 1]
413 let e_open_cand: i64 = h_prev + gap_first
414 let e_ext_cand: i64 = e_prev + gap_extend
415
416 if op_count >= max_ops { return -1 }
417 rev_ops[op_count] = NX_CIGAR_D & 0xff
418 op_count = op_count + 1
419 cj = cj - 1
420
421 // Decide next state: extend > open tie-break.
422 if e_cur == e_ext_cand {
423 // stay STATE_E (continue extension)
424 } else {
425 if e_cur == e_open_cand {
426 state = NX_AFFINE_STATE_H // gap opened from H
427 } else {
428 loop_active = 0 // defensive
429 }
430 }
431 }
432 } else {
433 // STATE_F
434 if ci <= 0 { loop_active = 0 }
435 else {
436 let cur_off: i64 = ci * cols
437 let f_cur: i64 = f[cur_off + cj]
438 let prev_row_off: i64 = (ci - 1) * cols
439 let h_prev: i64 = h[prev_row_off + cj]
440 let f_prev: i64 = f[prev_row_off + cj]
441 let f_open_cand: i64 = h_prev + gap_first
442 let f_ext_cand: i64 = f_prev + gap_extend
443
444 if op_count >= max_ops { return -1 }
445 rev_ops[op_count] = NX_CIGAR_I & 0xff
446 op_count = op_count + 1
447 ci = ci - 1
448
449 if f_cur == f_ext_cand {
450 // stay STATE_F
451 } else {
452 if f_cur == f_open_cand {
453 state = NX_AFFINE_STATE_H
454 } else {
455 loop_active = 0
456 }
457 }
458 }
459 }
460 }
461 }
462
463 out_meta[3] = ci
464 out_meta[4] = cj
465
466 // Reverse rev_ops into out_ops.
467 var k: i64 = 0
468 while k < op_count {
469 out_ops[k] = rev_ops[op_count - 1 - k]
470 k = k + 1
471 }
472
473 return op_count
474}