nx_align.nx source

1// nx_align.nx -- Smith-Waterman local alignment, linear-gap. 2// 3// license_tier: INDEPENDENT_REDERIVE 4// genealogy_id: international-research-sources/smith-waterman-1981 5// 6// G1.0 of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. Reference-impl 7// Smith-Waterman with linear gap penalty. This is the EXTENSION 8// half of BWA-MEM-class seed-and-extend; future modules: 9// 10// nx_align_minimizer.nx (G1.1) -- minimizer-based seeding 11// nx_align_chain.nx (G1.2) -- compatible-seed chaining 12// nx_align_affine.nx (G1.0b) -- affine gap (Gotoh 1982) 13// nx_align_bwa_mem.nx (G1.3) -- full BWA-MEM-class pipeline 14// 15// Recurrence (Smith & Waterman 1981): 16// H[0][j] = H[i][0] = 0 17// H[i][j] = max( 0, 18// H[i-1][j-1] + s(a[i-1], b[j-1]), // diagonal 19// H[i-1][j] + gap, // up (gap in b) 20// H[i][j-1] + gap ) // left (gap in a) 21// max_score = max over all H[i][j] 22// 23// Why local + linear-gap for G1.0: 24// - Local (the 0-floor) is what read aligners need: a read may 25// match only part of the reference, with mismatch / clipping 26// at the ends. Global alignment (Needleman-Wunsch) forces 27// end-to-end match which is wrong for read alignment. 28// - Linear gap is one parameter; affine gap (Gotoh 1982) is three 29// parameters and three DP matrices. Affine is more realistic 30// (a single 5bp indel is usually one event, not 5 events) but 31// adds complexity. G1.0b layers affine on top of this. 32// 33// Tie-breaking: 34// On equal scores, the FIRST max found wins (strict > update). 35// Means: when multiple equally-good local alignments exist, the 36// one closest to the top-left of the DP table is reported. 37// Deterministic across hosts; matches the audit discipline. 38// 39// API: 40// nx_max3(a, b, c) -> i64 41// nx_score_match(a_byte, b_byte, match, mis) -> i64 42// smith_waterman_linear( 43// a, n, b, m, 44// match_score, mismatch_score, gap_score, 45// out_max_i, out_max_j) -> i64 max score 46// 47// Memory: 48// Internal sys_mmap of (n+1) * (m+1) * 8 bytes for the DP table. 49// For typical short reads (n=100-300) against window (m=300-1000) 50// this is 80-300 KB. Long-read or whole-chromosome alignment 51// requires the banded / striped variants (G1.2 / G1.3). 52// 53// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml) 54// intended_use: "Local alignment extension half of read-to- 55// reference alignment + CRISPR off-target 56// scoring + primer-target similarity" 57// sil_target: SIL2 58// asil_target: QM 59// dal_target: DAL C 60// iec_62304_class: B 61// evidence: [no_floating_point, deterministic, 62// bit_equal_reproducible, 63// smith_waterman_1981_textbook_recurrence, 64// published_kat_vectors_verified, 65// strict_greater_tie_break_stable, 66// license_tier_INDEPENDENT_REDERIVE] 67// hazard_register: [bug-tape-sw-zero-floor-missed-negative, 68// bug-tape-sw-tie-break-nondeterministic, 69// bug-tape-sw-table-mn-overflow-page, 70// bug-tape-sw-max-tracking-wrong-end-position] 71// residual_risk: "G1.0 is O(mn) memory + time -- fine for 72// single-seed extension (n ~ 50-300, m ~ 50-1000) 73// but quadratic blow-up beyond. Banded / striped / 74// SIMD-vectorised variants live in G1.2-G1.3." 75// verdict: NOT_YET_EVALUATED 76 77import "nx_syscalls.nx" 78 79// 3-way max on i64. 80func nx_max3(a: i64, b: i64, c: i64) -> i64 { 81 var m: i64 = a 82 if b > m { m = b } 83 if c > m { m = c } 84 return m 85} 86 87// Single-base score: match_score if equal, else mismatch_score. 88// Bytes masked to 0..255 before compare; works for ASCII bases, 89// packed 2-bit codes one per byte, or any byte alphabet. 90func nx_score_match(a_byte: i64, b_byte: i64, 91 match_score: i64, mismatch_score: i64) -> i64 { 92 if (a_byte & 0xff) == (b_byte & 0xff) { return match_score } 93 return mismatch_score 94} 95 96// Smith-Waterman local alignment with linear gap penalty. 97// Returns the maximum H-cell value (>= 0); writes the (i, j) of that 98// max cell into out_max_i, out_max_j (both 1-indexed into the DP 99// table, so the alignment ends at a[out_max_i-1] and b[out_max_j-1] 100// in the original sequences). 101// 102// On the degenerate cases n==0 or m==0, returns 0 with (0, 0). 103func smith_waterman_linear(a: *u8, n: i64, 104 b: *u8, m: i64, 105 match_score: i64, 106 mismatch_score: i64, 107 gap_score: i64, 108 out_max_i: *i64, 109 out_max_j: *i64) -> i64 { 110 if n <= 0 { 111 out_max_i[0] = 0 112 out_max_j[0] = 0 113 return 0 114 } 115 if m <= 0 { 116 out_max_i[0] = 0 117 out_max_j[0] = 0 118 return 0 119 } 120 121 let cols: i64 = m + 1 122 let rows: i64 = n + 1 123 let bytes: i64 = rows * cols * 8 124 let h: *i64 = sys_mmap(bytes) as *i64 125 126 // sys_mmap returns zeroed pages, so H[0][j] = H[i][0] = 0 holds 127 // automatically. We do not re-write the boundary. 128 129 var best_score: i64 = 0 130 var best_i: i64 = 0 131 var best_j: i64 = 0 132 133 var i: i64 = 1 134 while i <= n { 135 let a_byte: i64 = a[i - 1] & 0xff 136 let row_off: i64 = i * cols 137 let prev_row_off: i64 = (i - 1) * cols 138 var j: i64 = 1 139 while j <= m { 140 let b_byte: i64 = b[j - 1] & 0xff 141 142 let s_diag: i64 = h[prev_row_off + j - 1] + 143 nx_score_match(a_byte, b_byte, match_score, mismatch_score) 144 let s_up: i64 = h[prev_row_off + j] + gap_score 145 let s_left: i64 = h[row_off + j - 1] + gap_score 146 147 var cell: i64 = nx_max3(s_diag, s_up, s_left) 148 if cell < 0 { cell = 0 } // local-alignment 0-floor 149 h[row_off + j] = cell 150 151 if cell > best_score { 152 best_score = cell 153 best_i = i 154 best_j = j 155 } 156 j = j + 1 157 } 158 i = i + 1 159 } 160 161 out_max_i[0] = best_i 162 out_max_j[0] = best_j 163 return best_score 164}