nx_align_score.nx source
↩ module page · 153 lines · 6295 B
1// nx_align_score.nx -- score chained alignment via Smith-Waterman.
2//
3// license_tier: INDEPENDENT_REDERIVE
4// genealogy_id: international-research-sources/li-2018-minimap2-seed-and-extend
5//
6// G1.4 of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. The integration
7// primitive that closes the seed-and-extend loop:
8//
9// query+ref packed DNA
10// -> nx_align_minimizer (G1.1)
11// -> nx_align_match (G1.2)
12// -> caller pre-sort by r_pos
13// -> nx_align_chain (G1.3)
14// -> THIS PRIMITIVE (G1.4) -> alignment score + region
15//
16// Given a chain (indices into the original seed arrays), extract the
17// implied alignment region:
18// q_start = seeds_q[chain[0]]
19// q_end = seeds_q[chain[chain_n-1]] + k // last seed spans k bases
20// r_start = seeds_r[chain[0]]
21// r_end = seeds_r[chain[chain_n-1]] + k
22//
23// Then unpack the 2-bit DNA in [q_start, q_end) and [r_start, r_end)
24// into per-byte 2-bit codes and run Smith-Waterman over them. The
25// returned score is the local-alignment score for the implied
26// region; it is >= chain_n * k * match_score under exact match and
27// degrades with gaps + mismatches in the inter-seed regions.
28//
29// Why not score the chain seeds alone:
30// - A chain proves co-linearity of seeds but says nothing about
31// the inter-seed gap quality. Two seeds 3bp apart in q and
32// 50bp apart in r implies a large insertion that should drop
33// the alignment score. Only SW captures that geometry.
34// - Real BWA-MEM / minimap2 use the chain to GUIDE extension
35// (banded SW around chain diagonal) rather than re-aligning
36// the full region; that optimization is G1.5. For G1.4
37// reference correctness we re-align the full region.
38//
39// What G1.4 does NOT do (deferred):
40// - Banded SW restricted to chain diagonal (G1.5 perf path)
41// - Affine gap (G1.0b)
42// - Soft-clipping / chimera detection (G1.4b)
43// - Mapping quality computation (G1.4b)
44// - Multi-chain reporting (primary + supplementary) (G1.4b)
45// - Strand-aware extension (forward + revcomp query) (G1.3c)
46//
47// API:
48// score_chained_alignment(
49// q_bases, q_len, r_bases, r_len,
50// seeds_q, seeds_r,
51// chain_idx, chain_n,
52// k,
53// match_s, mis_s, gap_s,
54// out_region, // [q_start, q_end, r_start, r_end] -- 4 i64s
55// out_sw_end // [max_i, max_j] -- 2 i64s
56// ) -> i64 SW score over the region (>= 0); -1 on bad input
57//
58// Output bundles avoid hitting nxc2 per-call arg limits while
59// preserving deterministic output ordering.
60//
61// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml)
62// intended_use: "Final scoring stage of seed-and-extend
63// read alignment; the integration primitive
64// that composes minimizer + match + chain +
65// Smith-Waterman into one alignment score"
66// sil_target: SIL2
67// asil_target: QM
68// dal_target: DAL C
69// iec_62304_class: B
70// evidence: [no_floating_point, deterministic,
71// bit_equal_reproducible,
72// composes_G0_1+G1_0+G1_3_KAT_chain,
73// end_to_end_pipeline_KAT_acgt_fixture,
74// license_tier_INDEPENDENT_REDERIVE]
75// hazard_register: [bug-tape-chain-endpoint-k-off-by-one,
76// bug-tape-region-out-of-bounds-vs-q_len,
77// bug-tape-unpack-region-byte-walk-off-end,
78// bug-tape-score-lower-than-perfect-match-on-exact]
79// residual_risk: "G1.4 unpacks regions to per-byte 2-bit
80// codes for SW; sufficient for KAT and
81// small-region perf but for chromosome-scale
82// chains needs banded-SW + packed-SW (G1.5)."
83// verdict: NOT_YET_EVALUATED
84
85import "nx_syscalls.nx"
86import "nx_sequence.nx"
87import "nx_align.nx"
88
89// Unpack a packed 2-bit DNA region into per-byte 2-bit codes
90// (one base per output byte, value 0..3). out must have capacity
91// >= count bytes. Returns 0.
92func unpack_region_2bit(bases: *u8, start: i64, count: i64, out: *u8) -> i64 {
93 var i: i64 = 0
94 while i < count {
95 let code: i64 = dna_get_base(bases, start + i)
96 out[i] = code & 0xff
97 i = i + 1
98 }
99 return 0
100}
101
102// Score a chain via Smith-Waterman over its implied region.
103//
104// out_region writes 4 i64s: [q_start, q_end, r_start, r_end]
105// out_sw_end writes 2 i64s: [sw_max_i, sw_max_j]
106func score_chained_alignment(q_bases: *u8, q_len: i64,
107 r_bases: *u8, r_len: i64,
108 seeds_q: *i64, seeds_r: *i64,
109 chain_idx: *i64, chain_n: i64,
110 k: i64,
111 match_s: i64, mis_s: i64, gap_s: i64,
112 out_region: *i64,
113 out_sw_end: *i64) -> i64 {
114 if chain_n <= 0 { return -1 }
115 if k <= 0 { return -1 }
116 if q_len <= 0 { return -1 }
117 if r_len <= 0 { return -1 }
118
119 let first_idx: i64 = chain_idx[0]
120 let last_idx: i64 = chain_idx[chain_n - 1]
121 let q_start: i64 = seeds_q[first_idx]
122 let r_start: i64 = seeds_r[first_idx]
123 let q_end: i64 = seeds_q[last_idx] + k
124 let r_end: i64 = seeds_r[last_idx] + k
125
126 if q_start < 0 { return -1 }
127 if r_start < 0 { return -1 }
128 if q_end > q_len { return -1 }
129 if r_end > r_len { return -1 }
130
131 out_region[0] = q_start
132 out_region[1] = q_end
133 out_region[2] = r_start
134 out_region[3] = r_end
135
136 let q_span: i64 = q_end - q_start
137 let r_span: i64 = r_end - r_start
138
139 let q_buf: *u8 = sys_mmap(q_span + 8)
140 let r_buf: *u8 = sys_mmap(r_span + 8)
141 unpack_region_2bit(q_bases, q_start, q_span, q_buf)
142 unpack_region_2bit(r_bases, r_start, r_span, r_buf)
143
144 let sw_mi: *i64 = sys_mmap(16) as *i64
145 let sw_mj: *i64 = sys_mmap(16) as *i64
146 let score: i64 = smith_waterman_linear(q_buf, q_span,
147 r_buf, r_span,
148 match_s, mis_s, gap_s,
149 sw_mi, sw_mj)
150 out_sw_end[0] = sw_mi[0]
151 out_sw_end[1] = sw_mj[0]
152 return score
153}