nx_align_affine.nx source
↩ module page · 199 lines · 7640 B
1// nx_align_affine.nx -- Smith-Waterman local alignment, AFFINE gap.
2//
3// license_tier: INDEPENDENT_REDERIVE
4// genealogy_id: international-research-sources/gotoh-1982
5//
6// G1.0b of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. Three-matrix DP
7// (Gotoh 1982) layering an affine gap penalty on top of the linear-
8// gap reference impl in nx_align.nx (G1.0).
9//
10// Why affine: a single 5bp deletion is biologically ONE event, not
11// 5 independent events. Linear-gap penalises it 5*gap; affine
12// penalises it open + 5*extend with extend << open, which matches
13// the empirical indel-length distribution and the model that
14// BWA-MEM / minimap2 / Bowtie2 all use by default.
15//
16// Penalty model:
17// Cost of a length-L gap = gap_open_score + L * gap_extend_score
18// Both scores are NEGATIVE (added to running scores).
19// BWA-MEM defaults: match=+1, mis=-4, gap_open=-6, gap_extend=-1.
20// Test convention here: match=+2, mis=-1, gap_open=-2, gap_extend=-1
21// so length-1 gap = -3 (= linear gap=-3), length-2 gap = -4
22// (vs linear -6) -- this differentiation is the whole point.
23//
24// Recurrence (Gotoh 1982):
25// E[i][j] = max( H[i ][j-1] + (gap_open + gap_extend),
26// E[i ][j-1] + gap_extend ) // gap in a
27// F[i][j] = max( H[i-1][j ] + (gap_open + gap_extend),
28// F[i-1][j ] + gap_extend ) // gap in b
29// H[i][j] = max( 0,
30// H[i-1][j-1] + s(a[i-1], b[j-1]),
31// E[i][j],
32// F[i][j] )
33// Boundary: H[0][*] = H[*][0] = 0; E[*][0] = F[0][*] = NEG_INF
34// max_score = max over all H[i][j]
35//
36// Negative-infinity sentinel:
37// NEG_INF = -1000000000 (i.e., -1e9). All real scores fit in
38// roughly ±(n * match_score) << 1e9, so saturating-add semantics
39// are avoided. Bounded n * m * |max_score_per_cell| << 1e9.
40//
41// What G1.0b does NOT do (deferred):
42// - Backtrace (alignment string output) -- G1.4b
43// - Banded SW around chain diagonal -- G1.5
44// - Striped/SIMD vectorisation -- G1.5
45// - Two-piece affine (long-gap penalty) -- G1.7
46// - Position-specific scoring matrices -- G2.1 (variant calling)
47//
48// API:
49// smith_waterman_affine(
50// a, n, b, m,
51// match_score, mismatch_score,
52// gap_open_score, gap_extend_score,
53// out_max_i, out_max_j) -> i64 max score (>= 0)
54//
55// Memory:
56// Internal sys_mmap of 3 * (n+1) * (m+1) * 8 bytes (H, E, F tables).
57// Triple the linear-gap variant. For n=m=200 -> ~1 MB.
58//
59// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml)
60// intended_use: "Affine-gap local alignment -- the realistic
61// gap model used by BWA-MEM / minimap2 /
62// Bowtie2. Replaces or composes with
63// smith_waterman_linear (G1.0) when biology
64// of long single-event indels matters."
65// sil_target: SIL2
66// asil_target: QM
67// dal_target: DAL C
68// iec_62304_class: B
69// evidence: [no_floating_point, deterministic,
70// bit_equal_reproducible,
71// gotoh_1982_three_matrix_recurrence,
72// degenerate_open_zero_equals_linear_KAT,
73// two_bp_gap_outscores_linear_KAT,
74// license_tier_INDEPENDENT_REDERIVE]
75// hazard_register: [bug-tape-affine-e-f-init-not-neg-inf,
76// bug-tape-affine-open-cost-double-counted,
77// bug-tape-affine-h-max-misses-e-or-f,
78// bug-tape-affine-neg-inf-arithmetic-overflow,
79// bug-tape-affine-zero-floor-applied-to-e-f]
80// residual_risk: "NEG_INF = -1e9 is safe for n*m up to ~1e8;
81// substrate must enforce n*m bound or upgrade
82// sentinel before use on whole-genome scale."
83// verdict: NOT_YET_EVALUATED
84
85import "nx_syscalls.nx"
86import "nx_align.nx"
87import "nx_const.nx"
88
89// Negative-infinity sentinel for E/F boundary cells is
90// NX_NEG_INF_I64 from nx_const. NishiLang's `const` requires
91// integer literals so we reference NX_NEG_INF_I64 directly
92// throughout this module rather than re-aliasing. See
93// docs/NISHI_CONSTANT_REGISTRY.md.
94
95// 4-way max on i64. Inlines nx_max3 + one more compare.
96func nx_max4(a: i64, b: i64, c: i64, d: i64) -> i64 {
97 var m: i64 = a
98 if b > m { m = b }
99 if c > m { m = c }
100 if d > m { m = d }
101 return m
102}
103
104// Smith-Waterman local alignment with affine gap penalty (Gotoh 1982).
105// gap_open_score + gap_extend_score are both expected NEGATIVE.
106func smith_waterman_affine(a: *u8, n: i64,
107 b: *u8, m: i64,
108 match_score: i64,
109 mismatch_score: i64,
110 gap_open_score: i64,
111 gap_extend_score: i64,
112 out_max_i: *i64,
113 out_max_j: *i64) -> i64 {
114 if n <= 0 {
115 out_max_i[0] = 0
116 out_max_j[0] = 0
117 return 0
118 }
119 if m <= 0 {
120 out_max_i[0] = 0
121 out_max_j[0] = 0
122 return 0
123 }
124
125 let cols: i64 = m + 1
126 let rows: i64 = n + 1
127 let bytes: i64 = rows * cols * 8
128
129 let h: *i64 = sys_mmap(bytes) as *i64
130 let e: *i64 = sys_mmap(bytes) as *i64
131 let f: *i64 = sys_mmap(bytes) as *i64
132
133 // sys_mmap returns zeroed pages, so H[0][*] = H[*][0] = 0 already.
134 // E and F need NEG_INF on their boundaries to prevent the
135 // "open a gap into nothing" path from contributing a real score.
136 var j_init: i64 = 0
137 while j_init <= m {
138 f[j_init] = NX_NEG_INF_I64 // F[0][j]
139 j_init = j_init + 1
140 }
141 var i_init: i64 = 0
142 while i_init <= n {
143 e[i_init * cols] = NX_NEG_INF_I64 // E[i][0]
144 i_init = i_init + 1
145 }
146
147 // gap_first = gap_open + gap_extend = total cost of a length-1 gap.
148 let gap_first: i64 = gap_open_score + gap_extend_score
149
150 var best_score: i64 = 0
151 var best_i: i64 = 0
152 var best_j: i64 = 0
153
154 var i: i64 = 1
155 while i <= n {
156 let a_byte: i64 = a[i - 1] & 0xff
157 let row_off: i64 = i * cols
158 let prev_row_off: i64 = (i - 1) * cols
159
160 var j: i64 = 1
161 while j <= m {
162 let b_byte: i64 = b[j - 1] & 0xff
163
164 // E[i][j]: gap in a (consume b[j-1], not a[i-1]).
165 let e_open: i64 = h[row_off + j - 1] + gap_first
166 let e_ext: i64 = e[row_off + j - 1] + gap_extend_score
167 var e_cell: i64 = e_open
168 if e_ext > e_cell { e_cell = e_ext }
169 e[row_off + j] = e_cell
170
171 // F[i][j]: gap in b (consume a[i-1], not b[j-1]).
172 let f_open: i64 = h[prev_row_off + j] + gap_first
173 let f_ext: i64 = f[prev_row_off + j] + gap_extend_score
174 var f_cell: i64 = f_open
175 if f_ext > f_cell { f_cell = f_ext }
176 f[row_off + j] = f_cell
177
178 // H[i][j]: max(0, diag, E, F).
179 var sub_pen: i64 = mismatch_score
180 if a_byte == b_byte { sub_pen = match_score }
181 let h_diag: i64 = h[prev_row_off + j - 1] + sub_pen
182
183 var h_cell: i64 = nx_max4(0, h_diag, e_cell, f_cell)
184 h[row_off + j] = h_cell
185
186 if h_cell > best_score {
187 best_score = h_cell
188 best_i = i
189 best_j = j
190 }
191 j = j + 1
192 }
193 i = i + 1
194 }
195
196 out_max_i[0] = best_i
197 out_max_j[0] = best_j
198 return best_score
199}