nx_align_score_test.nx source
↩ module page · 185 lines · 7023 B
1// nx_align_score_test.nx -- end-to-end seed-and-extend pipeline KAT.
2//
3// Full chain composed in one test:
4// 1. pack reference "ACGTACGTACGT" (12 bases) and query "ACGTACGT" (8 bases)
5// 2. minimizer_extract from each (k=4 w=3)
6// 3. minimizer_match to produce seed pairs
7// 4. hand-sort seed pairs by (r_pos, q_pos)
8// 5. seed_chain to find longest co-linear chain
9// 6. score_chained_alignment to run SW over the chain region
10//
11// Hand-traced expected values:
12//
13// Reference minimizers (5): (0x1B, 0) (0x6C, 1) (0x1B, 4) (0x6C, 5) (0x1B, 8)
14// Query minimizers (3): (0x1B, 0) (0x6C, 1) (0x1B, 4)
15//
16// Match pairs (outer q, inner r) = 8:
17// q_idx 0 (v=0x1B,q_pos=0): (0,0) (0,4) (0,8)
18// q_idx 1 (v=0x6C,q_pos=1): (1,1) (1,5)
19// q_idx 2 (v=0x1B,q_pos=4): (4,0) (4,4) (4,8)
20//
21// After sort by (r,q) ascending:
22// index 0: (q=0, r=0)
23// index 1: (q=4, r=0)
24// index 2: (q=1, r=1)
25// index 3: (q=0, r=4)
26// index 4: (q=4, r=4)
27// index 5: (q=1, r=5)
28// index 6: (q=0, r=8)
29// index 7: (q=4, r=8)
30//
31// Chain DP picks chain length 3 at end-index 4 (leftmost tie):
32// indices [0, 2, 4] -> seeds (q=0,r=0), (q=1,r=1), (q=4,r=4)
33//
34// Alignment region:
35// q_start = 0, q_end = 4 + 4 = 8
36// r_start = 0, r_end = 4 + 4 = 8
37// q[0..8) = "ACGTACGT" == r[0..8) -> perfect 8-base match
38// SW score (m=+2 mis=-1 gap=-2) = 8 * 2 = 16, end (8, 8)
39//
40// expect_exit: 0
41//
42// license_tier: ORIGINAL
43
44import "nx_syscalls.nx"
45import "nx_sequence.nx"
46import "nx_align.nx"
47import "nx_align_minimizer.nx"
48import "nx_align_match.nx"
49import "nx_align_chain.nx"
50import "nx_align_score.nx"
51
52func main() -> i64 {
53
54 // ============================================================
55 // Stage 1 -- pack reference + query into 2-bit form.
56 // ============================================================
57
58 let r_bases: *u8 = sys_mmap(8)
59 r_bases[0] = 0x1B // ACGT
60 r_bases[1] = 0x1B // ACGT
61 r_bases[2] = 0x1B // ACGT
62
63 let q_bases: *u8 = sys_mmap(8)
64 q_bases[0] = 0x1B // ACGT
65 q_bases[1] = 0x1B // ACGT
66
67 // ============================================================
68 // Stage 2 -- extract minimizers from each.
69 // ============================================================
70
71 let r_mins_v: *i64 = sys_mmap(128) as *i64
72 let r_mins_p: *i64 = sys_mmap(128) as *i64
73 let r_n: i64 = minimizer_extract(r_bases, 12, 4, 3, r_mins_v, r_mins_p, 16)
74 if r_n != 5 { return 1 }
75 if r_mins_v[0] != 0x1B { return 2 }
76 if r_mins_p[0] != 0 { return 3 }
77 if r_mins_v[4] != 0x1B { return 4 }
78 if r_mins_p[4] != 8 { return 5 }
79
80 let q_mins_v: *i64 = sys_mmap(128) as *i64
81 let q_mins_p: *i64 = sys_mmap(128) as *i64
82 let q_n: i64 = minimizer_extract(q_bases, 8, 4, 3, q_mins_v, q_mins_p, 16)
83 if q_n != 3 { return 10 }
84
85 // ============================================================
86 // Stage 3 -- match query minimizers against reference.
87 // ============================================================
88
89 let mq: *i64 = sys_mmap(256) as *i64
90 let mr: *i64 = sys_mmap(256) as *i64
91 let n_pairs: i64 = minimizer_match(q_mins_v, q_mins_p, q_n,
92 r_mins_v, r_mins_p, r_n,
93 mq, mr, 64)
94 if n_pairs != 8 { return 20 }
95 // Spot-check emission order: outer q first, inner r in r-array order.
96 if mq[0] != 0 { return 21 }
97 if mr[0] != 0 { return 22 }
98 if mq[7] != 4 { return 23 }
99 if mr[7] != 8 { return 24 }
100
101 // ============================================================
102 // Stage 4 -- hand-sort the 8 pairs by (r_pos, q_pos) ascending.
103 // ============================================================
104
105 let sq: *i64 = sys_mmap(256) as *i64
106 let sr: *i64 = sys_mmap(256) as *i64
107 sq[0]=0; sr[0]=0
108 sq[1]=4; sr[1]=0
109 sq[2]=1; sr[2]=1
110 sq[3]=0; sr[3]=4
111 sq[4]=4; sr[4]=4
112 sq[5]=1; sr[5]=5
113 sq[6]=0; sr[6]=8
114 sq[7]=4; sr[7]=8
115
116 // ============================================================
117 // Stage 5 -- chain DP.
118 // ============================================================
119
120 let chain: *i64 = sys_mmap(128) as *i64
121 let chain_n: i64 = seed_chain(sq, sr, 8, chain, 16)
122 if chain_n != 3 { return 30 }
123 if chain[0] != 0 { return 31 }
124 if chain[1] != 2 { return 32 }
125 if chain[2] != 4 { return 33 }
126
127 // ============================================================
128 // Stage 6 -- score the chained alignment via SW.
129 // ============================================================
130
131 let region: *i64 = sys_mmap(64) as *i64
132 let sw_end: *i64 = sys_mmap(32) as *i64
133
134 let score: i64 = score_chained_alignment(q_bases, 8,
135 r_bases, 12,
136 sq, sr,
137 chain, chain_n,
138 4,
139 2, -1, -2,
140 region, sw_end)
141 if score != 16 { return 40 }
142 if region[0] != 0 { return 41 } // q_start
143 if region[1] != 8 { return 42 } // q_end
144 if region[2] != 0 { return 43 } // r_start
145 if region[3] != 8 { return 44 } // r_end
146 if sw_end[0] != 8 { return 45 } // sw_max_i
147 if sw_end[1] != 8 { return 46 } // sw_max_j
148
149 // ============================================================
150 // Stage 7 -- isolated G1.4 single-seed chain (length 1).
151 // chain=[0]: q=0, r=0; region [0,4) on both -> ACGT vs ACGT -> 8.
152 // ============================================================
153
154 let s1q: *i64 = sys_mmap(16) as *i64
155 let s1r: *i64 = sys_mmap(16) as *i64
156 s1q[0] = 0; s1r[0] = 0
157 let ch1: *i64 = sys_mmap(16) as *i64
158 ch1[0] = 0
159 let score1: i64 = score_chained_alignment(q_bases, 8,
160 r_bases, 12,
161 s1q, s1r,
162 ch1, 1,
163 4,
164 2, -1, -2,
165 region, sw_end)
166 if score1 != 8 { return 50 }
167 if region[0] != 0 { return 51 } // q_start
168 if region[1] != 4 { return 52 } // q_end
169
170 // ============================================================
171 // Stage 8 -- bad input rejection.
172 // ============================================================
173
174 if score_chained_alignment(q_bases, 8, r_bases, 12,
175 sq, sr, chain, 0, // chain_n=0
176 4, 2, -1, -2,
177 region, sw_end) != -1 { return 60 }
178
179 if score_chained_alignment(q_bases, 8, r_bases, 12,
180 sq, sr, chain, 3,
181 0, 2, -1, -2, // k=0
182 region, sw_end) != -1 { return 61 }
183
184 return 0
185}