code wiki / (root) / nx_align_affine_test.nx

nx_align_affine_test.nx source

↩ module page · 139 lines · 6168 B

1// nx_align_affine_test.nx -- KAT for Gotoh 1982 affine-gap SW. 2// 3// Test strategy: 4// 1. Pure match -- gaps irrelevant -- affine == linear == n * match. 5// 2. Degenerate affine (open=0) reduces to linear with gap=extend. 6// 3. Affine advantage cases -- length-2 / length-3 gaps where 7// affine differentiates favourably vs. equivalent linear. 8// 9// expect_exit: 0 10// 11// license_tier: ORIGINAL 12 13import "nx_syscalls.nx" 14import "nx_align.nx" 15import "nx_align_affine.nx" 16 17func main() -> i64 { 18 let a: *u8 = sys_mmap(32) 19 let b: *u8 = sys_mmap(32) 20 let mi: *i64 = sys_mmap(16) as *i64 21 let mj: *i64 = sys_mmap(16) as *i64 22 23 // ============================================================ 24 // Section A -- identical "ACGT" vs "ACGT": no gaps, score 8. 25 // ============================================================ 26 27 a[0]=0x41; a[1]=0x43; a[2]=0x47; a[3]=0x54 28 b[0]=0x41; b[1]=0x43; b[2]=0x47; b[3]=0x54 29 let s1: i64 = smith_waterman_affine(a, 4, b, 4, 2, -1, -2, -1, mi, mj) 30 if s1 != 8 { return 1 } 31 if mi[0] != 4 { return 2 } 32 if mj[0] != 4 { return 3 } 33 34 // ============================================================ 35 // Section B -- pure match invariant under extreme gap penalties. 36 // open=-100 extend=-100 should still return 8 because no gap is 37 // ever opened along the optimal path. 38 // ============================================================ 39 40 let s2: i64 = smith_waterman_affine(a, 4, b, 4, 2, -1, -100, -100, mi, mj) 41 if s2 != 8 { return 5 } 42 43 // ============================================================ 44 // Section C -- degenerate affine reduces to linear. 45 // open=0 extend=-3 -> length-L gap = L * -3, identical to linear -3. 46 // Replays the gap-in-A test from nx_align_test: 47 // "ACT" vs "AGCT" expected score 4, end (3, 4). 48 // ============================================================ 49 50 a[0]=0x41; a[1]=0x43; a[2]=0x54; a[3]=0 51 b[0]=0x41; b[1]=0x47; b[2]=0x43; b[3]=0x54 52 let s3a: i64 = smith_waterman_affine(a, 3, b, 4, 2, -1, 0, -3, mi, mj) 53 if s3a != 4 { return 10 } 54 if mi[0] != 3 { return 11 } 55 if mj[0] != 4 { return 12 } 56 57 // Cross-check: linear with gap=-3 produces the same on same input. 58 let s3b: i64 = smith_waterman_linear(a, 3, b, 4, 2, -1, -3, mi, mj) 59 if s3b != s3a { return 13 } 60 61 // ============================================================ 62 // Section D -- 1bp gap: affine vs linear give same answer. 63 // Fixture: "AAAAGGGG" (8) vs "AAAATGGGG" (9, insert T). 64 // AAAA-GGGG / AAAATGGGG -> 8 matches + 1 gap of length 1 65 // open=-2 extend=-1 -> length-1 gap = -3 = linear gap=-3 66 // Score = 8*2 - 3 = 13 both ways. 67 // ============================================================ 68 69 a[0]=0x41; a[1]=0x41; a[2]=0x41; a[3]=0x41 70 a[4]=0x47; a[5]=0x47; a[6]=0x47; a[7]=0x47 71 b[0]=0x41; b[1]=0x41; b[2]=0x41; b[3]=0x41 72 b[4]=0x54; b[5]=0x47; b[6]=0x47; b[7]=0x47; b[8]=0x47 73 let s4a: i64 = smith_waterman_affine(a, 8, b, 9, 2, -1, -2, -1, mi, mj) 74 if s4a != 13 { return 20 } 75 let s4b: i64 = smith_waterman_linear(a, 8, b, 9, 2, -1, -3, mi, mj) 76 if s4b != s4a { return 21 } 77 78 // ============================================================ 79 // Section E -- 2bp gap: AFFINE DIFFERENTIATES. 80 // Fixture: "AAAAAGGGGG" (10) vs "AAAAATTGGGGG" (12, insert TT). 81 // AAAAA--GGGGG / AAAAATTGGGGG -> 10 matches + length-2 gap 82 // Affine (open=-2 extend=-1): 10*2 - (2 + 2*1) = 16 83 // Linear (gap=-3): 10*2 - (2*3) = 14 84 // ============================================================ 85 86 a[0]=0x41; a[1]=0x41; a[2]=0x41; a[3]=0x41; a[4]=0x41 87 a[5]=0x47; a[6]=0x47; a[7]=0x47; a[8]=0x47; a[9]=0x47 88 b[0]=0x41; b[1]=0x41; b[2]=0x41; b[3]=0x41; b[4]=0x41 89 b[5]=0x54; b[6]=0x54 90 b[7]=0x47; b[8]=0x47; b[9]=0x47; b[10]=0x47; b[11]=0x47 91 let s5a: i64 = smith_waterman_affine(a, 10, b, 12, 2, -1, -2, -1, mi, mj) 92 if s5a != 16 { return 30 } 93 let s5b: i64 = smith_waterman_linear(a, 10, b, 12, 2, -1, -3, mi, mj) 94 if s5b != 14 { return 31 } 95 96 // ============================================================ 97 // Section F -- 3bp gap: AFFINE prefers gap-traversal, LINEAR gives up. 98 // Fixture: "AAAGGGG" (7) vs "AAATTTGGGG" (10, insert TTT). 99 // Gap-traversing alignment AAA---GGGG / AAATTTGGGG: 100 // Affine (open=-2 extend=-1): 7*2 - (2 + 3*1) = 9 101 // Linear (gap=-3): 7*2 - (3*3) = 5 102 // Local-only alternative -- just match GGGG (4 bases) -> 4*2 = 8 103 // Affine picks gap-traverse (9 > 8); linear picks GGGG-only (8 > 5). 104 // This IS the affine advantage in stark form: affine keeps the 105 // long alignment hypothesis alive; linear discards it. 106 // ============================================================ 107 108 a[0]=0x41; a[1]=0x41; a[2]=0x41 109 a[3]=0x47; a[4]=0x47; a[5]=0x47; a[6]=0x47 110 b[0]=0x41; b[1]=0x41; b[2]=0x41 111 b[3]=0x54; b[4]=0x54; b[5]=0x54 112 b[6]=0x47; b[7]=0x47; b[8]=0x47; b[9]=0x47 113 let s6a: i64 = smith_waterman_affine(a, 7, b, 10, 2, -1, -2, -1, mi, mj) 114 if s6a != 9 { return 40 } 115 let s6b: i64 = smith_waterman_linear(a, 7, b, 10, 2, -1, -3, mi, mj) 116 if s6b != 8 { return 41 } 117 118 // ============================================================ 119 // Section G -- degenerate empty input returns 0, (0, 0). 120 // ============================================================ 121 122 let s7: i64 = smith_waterman_affine(a, 0, b, 5, 2, -1, -2, -1, mi, mj) 123 if s7 != 0 { return 50 } 124 if mi[0] != 0 { return 51 } 125 let s8: i64 = smith_waterman_affine(a, 5, b, 0, 2, -1, -2, -1, mi, mj) 126 if s8 != 0 { return 52 } 127 128 // ============================================================ 129 // Section H -- all-mismatch returns 0 (local-alignment 0-floor). 130 // "AAAA" vs "TTTT" with sane scoring. 131 // ============================================================ 132 133 a[0]=0x41; a[1]=0x41; a[2]=0x41; a[3]=0x41 134 b[0]=0x54; b[1]=0x54; b[2]=0x54; b[3]=0x54 135 let s9: i64 = smith_waterman_affine(a, 4, b, 4, 2, -1, -2, -1, mi, mj) 136 if s9 != 0 { return 60 } 137 138 return 0 139}