nx_variant_test.nx source
↩ module page · 323 lines · 13424 B
1// nx_variant_test.nx -- KAT for pileup aggregator + first SNV caller.
2//
3// expect_exit: 0
4//
5// license_tier: ORIGINAL
6
7import "nx_syscalls.nx"
8import "nx_const.nx"
9import "nx_pileup.nx"
10import "nx_variant.nx"
11
12func main() -> i64 {
13
14 let counts: *i64 = sys_mmap(64) as *i64
15 let codes: *i64 = sys_mmap(128) as *i64
16 let flags: *i64 = sys_mmap(128) as *i64
17 let alt: *i64 = sys_mmap(16) as *i64
18 let altcnt: *i64 = sys_mmap(16) as *i64
19
20 // ============================================================
21 // Section A -- pileup_aggregate_counts: simple cases.
22 // 4 reads all contributing A (BASE flag) -> counts[A]=4, total=4.
23 // ============================================================
24
25 var i: i64 = 0
26 while i < 4 {
27 codes[i] = NX_DNA_A
28 flags[i] = NX_PILEUP_FLAG_BASE
29 i = i + 1
30 }
31 // Re-zero counts via sys_mmap fresh page.
32 let c1: *i64 = sys_mmap(64) as *i64
33 if pileup_aggregate_counts(codes, flags, 4, c1) != 4 { return 1 }
34 if c1[NX_DNA_A] != 4 { return 2 }
35 if c1[NX_DNA_C] != 0 { return 3 }
36 if c1[NX_PILEUP_COUNT_N] != 0 { return 4 }
37 if c1[NX_PILEUP_COUNT_DEL] != 0 { return 5 }
38
39 // ============================================================
40 // Section B -- mixed contributions across all flag types.
41 // 10 reads: 4A 3G 1N 1DEL 1OUT -> total contributing = 9.
42 // counts[A]=4, counts[G]=3, counts[N]=1, counts[DEL]=1
43 // ============================================================
44
45 let c2: *i64 = sys_mmap(64) as *i64
46 codes[0]=NX_DNA_A; flags[0]=NX_PILEUP_FLAG_BASE
47 codes[1]=NX_DNA_A; flags[1]=NX_PILEUP_FLAG_BASE
48 codes[2]=NX_DNA_A; flags[2]=NX_PILEUP_FLAG_BASE
49 codes[3]=NX_DNA_A; flags[3]=NX_PILEUP_FLAG_BASE
50 codes[4]=NX_DNA_G; flags[4]=NX_PILEUP_FLAG_BASE
51 codes[5]=NX_DNA_G; flags[5]=NX_PILEUP_FLAG_BASE
52 codes[6]=NX_DNA_G; flags[6]=NX_PILEUP_FLAG_BASE
53 codes[7]=NX_DNA_A; flags[7]=NX_PILEUP_FLAG_N // N flag; code ignored
54 codes[8]=NX_DNA_A; flags[8]=NX_PILEUP_FLAG_DEL // DEL flag
55 codes[9]=NX_DNA_A; flags[9]=NX_PILEUP_FLAG_OUT // OUT skipped
56 if pileup_aggregate_counts(codes, flags, 10, c2) != 9 { return 10 }
57 if c2[NX_DNA_A] != 4 { return 11 }
58 if c2[NX_DNA_C] != 0 { return 12 }
59 if c2[NX_DNA_G] != 3 { return 13 }
60 if c2[NX_DNA_T] != 0 { return 14 }
61 if c2[NX_PILEUP_COUNT_N] != 1 { return 15 }
62 if c2[NX_PILEUP_COUNT_DEL] != 1 { return 16 }
63
64 // ============================================================
65 // Section C -- snv_call: all-ref returns NONE.
66 // ref=A, 10A reads, min_cov=5 alt_freq=20% hom_freq=80% -> NONE
67 // ============================================================
68
69 let c3: *i64 = sys_mmap(64) as *i64
70 c3[NX_DNA_A] = 10
71 if snv_call(NX_DNA_A, c3, 5, 20, 80, alt, altcnt) != NX_VARIANT_NONE { return 20 }
72 if altcnt[0] != 0 { return 21 }
73
74 // ============================================================
75 // Section D -- all-alt homozygous variant.
76 // ref=A, 10G reads -> HOM with alt=G, count=10.
77 // alt_freq = 100% which is >= hom_freq_pct=80 -> HOM.
78 // ============================================================
79
80 let c4: *i64 = sys_mmap(64) as *i64
81 c4[NX_DNA_G] = 10
82 if snv_call(NX_DNA_A, c4, 5, 20, 80, alt, altcnt) != NX_VARIANT_HOM { return 30 }
83 if alt[0] != NX_DNA_G { return 31 }
84 if altcnt[0] != 10 { return 32 }
85
86 // ============================================================
87 // Section E -- 50/50 heterozygous variant.
88 // ref=A, 5A + 5G reads -> alt_freq=50%; >= 20 (min) and < 80 (hom) -> HET.
89 // ============================================================
90
91 let c5: *i64 = sys_mmap(64) as *i64
92 c5[NX_DNA_A] = 5
93 c5[NX_DNA_G] = 5
94 if snv_call(NX_DNA_A, c5, 5, 20, 80, alt, altcnt) != NX_VARIANT_HET { return 40 }
95 if alt[0] != NX_DNA_G { return 41 }
96 if altcnt[0] != 5 { return 42 }
97
98 // ============================================================
99 // Section F -- below min_coverage returns NONE.
100 // ref=A, 2A+1G reads (total=3), min_coverage=5 -> NONE.
101 // ============================================================
102
103 let c6: *i64 = sys_mmap(64) as *i64
104 c6[NX_DNA_A] = 2
105 c6[NX_DNA_G] = 1
106 if snv_call(NX_DNA_A, c6, 5, 20, 80, alt, altcnt) != NX_VARIANT_NONE { return 50 }
107
108 // ============================================================
109 // Section G -- alt frequency below threshold returns NONE.
110 // ref=A, 90A+10G (total=100), alt_freq=10%, min=20 -> NONE.
111 // ============================================================
112
113 let c7: *i64 = sys_mmap(64) as *i64
114 c7[NX_DNA_A] = 90
115 c7[NX_DNA_G] = 10
116 if snv_call(NX_DNA_A, c7, 5, 20, 80, alt, altcnt) != NX_VARIANT_NONE { return 60 }
117
118 // Bump min_alt to 10 -> now HET (since 10% >= 10% and < 80%).
119 if snv_call(NX_DNA_A, c7, 5, 10, 80, alt, altcnt) != NX_VARIANT_HET { return 61 }
120 if alt[0] != NX_DNA_G { return 62 }
121 if altcnt[0] != 10 { return 63 }
122
123 // ============================================================
124 // Section H -- HOM boundary: exactly at hom_freq threshold.
125 // ref=A, 2A+8G (total=10), alt_freq=80% -> HOM (>=).
126 // ============================================================
127
128 let c8: *i64 = sys_mmap(64) as *i64
129 c8[NX_DNA_A] = 2
130 c8[NX_DNA_G] = 8
131 if snv_call(NX_DNA_A, c8, 5, 20, 80, alt, altcnt) != NX_VARIANT_HOM { return 70 }
132
133 // Just below: 3A+7G -> 70% -> HET.
134 let c8b: *i64 = sys_mmap(64) as *i64
135 c8b[NX_DNA_A] = 3
136 c8b[NX_DNA_G] = 7
137 if snv_call(NX_DNA_A, c8b, 5, 20, 80, alt, altcnt) != NX_VARIANT_HET { return 71 }
138
139 // ============================================================
140 // Section I -- leftmost-tie-break across alt alleles.
141 // ref=A, 4C + 4G + 4T (total=12). Three alts tied at 4 each.
142 // Leftmost (smallest code) wins: C (code 1).
143 // ============================================================
144
145 let c9: *i64 = sys_mmap(64) as *i64
146 c9[NX_DNA_C] = 4
147 c9[NX_DNA_G] = 4
148 c9[NX_DNA_T] = 4
149 if snv_call(NX_DNA_A, c9, 5, 20, 80, alt, altcnt) != NX_VARIANT_HET { return 80 }
150 if alt[0] != NX_DNA_C { return 81 }
151 if altcnt[0] != 4 { return 82 }
152
153 // ============================================================
154 // Section J -- N/DEL counts excluded from coverage and decisions.
155 // ref=A, 4A + 0 alt + 100N + 100DEL -> total_ACGT=4 < min_cov=5 -> NONE
156 // (the N+DEL counts do NOT contribute to "informative coverage")
157 // ============================================================
158
159 let cA: *i64 = sys_mmap(64) as *i64
160 cA[NX_DNA_A] = 4
161 cA[NX_PILEUP_COUNT_N] = 100
162 cA[NX_PILEUP_COUNT_DEL] = 100
163 if snv_call(NX_DNA_A, cA, 5, 20, 80, alt, altcnt) != NX_VARIANT_NONE { return 90 }
164
165 // ============================================================
166 // Section K -- bad ref code rejection.
167 // ============================================================
168
169 if snv_call(-1, c5, 5, 20, 80, alt, altcnt) != -1 { return 100 }
170 if snv_call(4, c5, 5, 20, 80, alt, altcnt) != -1 { return 101 }
171
172 // ============================================================
173 // Section L -- indel_call (G2.3c).
174 // ============================================================
175
176 let ihist: *i64 = sys_mmap(128) as *i64
177 let elen: *i64 = sys_mmap(16) as *i64
178 let ecnt: *i64 = sys_mmap(16) as *i64
179
180 // Case 1: all no-insertion (ins_hist[0]=10, others 0) -> NONE.
181 ihist[0] = 10
182 if indel_call(ihist, 5, 10, 5, 20, 80, elen, ecnt) != NX_VARIANT_NONE { return 200 }
183 if elen[0] != 0 { return 201 }
184
185 // Case 2: 8 total, 7 with 2bp insertion -> 87% -> HOM at 80% threshold.
186 let ihist2: *i64 = sys_mmap(128) as *i64
187 ihist2[0] = 1
188 ihist2[2] = 7
189 if indel_call(ihist2, 5, 8, 5, 20, 80, elen, ecnt) != NX_VARIANT_HOM { return 210 }
190 if elen[0] != 2 { return 211 }
191 if ecnt[0] != 7 { return 212 }
192
193 // Case 3: 8 total, 6 with 1bp insertion -> 75% -> HET (>=20% min, <80% hom).
194 let ihist3: *i64 = sys_mmap(128) as *i64
195 ihist3[0] = 2
196 ihist3[1] = 6
197 if indel_call(ihist3, 5, 8, 5, 20, 80, elen, ecnt) != NX_VARIANT_HET { return 220 }
198 if elen[0] != 1 { return 221 }
199 if ecnt[0] != 6 { return 222 }
200
201 // Case 4: mixed lengths -- 3@1bp + 5@3bp out of 10 -> alt=3bp, 50% -> HET.
202 let ihist4: *i64 = sys_mmap(128) as *i64
203 ihist4[0] = 2
204 ihist4[1] = 3
205 ihist4[3] = 5
206 if indel_call(ihist4, 5, 10, 5, 20, 80, elen, ecnt) != NX_VARIANT_HET { return 230 }
207 if elen[0] != 3 { return 231 }
208 if ecnt[0] != 5 { return 232 }
209
210 // Case 5: tie between two lengths -- leftmost (smaller L) wins.
211 // 4@2bp + 4@5bp out of 10 -> alt=2bp.
212 let ihist5: *i64 = sys_mmap(128) as *i64
213 ihist5[0] = 2
214 ihist5[2] = 4
215 ihist5[5] = 4
216 if indel_call(ihist5, 5, 10, 5, 20, 80, elen, ecnt) != NX_VARIANT_HET { return 240 }
217 if elen[0] != 2 { return 241 }
218
219 // Case 6: below min_coverage -> NONE.
220 let ihist6: *i64 = sys_mmap(128) as *i64
221 ihist6[0] = 1; ihist6[1] = 2 // total = 3, min_coverage=10
222 if indel_call(ihist6, 5, 3, 10, 20, 80, elen, ecnt) != NX_VARIANT_NONE { return 250 }
223
224 // Case 7: alt-freq below threshold -> NONE.
225 let ihist7: *i64 = sys_mmap(128) as *i64
226 ihist7[0] = 90; ihist7[1] = 10 // total 100, alt 10%
227 if indel_call(ihist7, 5, 100, 10, 20, 80, elen, ecnt) != NX_VARIANT_NONE { return 260 }
228 // Lower threshold to 10% -> HET.
229 if indel_call(ihist7, 5, 100, 10, 10, 80, elen, ecnt) != NX_VARIANT_HET { return 261 }
230 if elen[0] != 1 { return 262 }
231
232 // Case 8: bad input.
233 if indel_call(ihist7, -1, 100, 10, 20, 80, elen, ecnt) != -1 { return 270 }
234 if indel_call(ihist7, 5, -1, 10, 20, 80, elen, ecnt) != -1 { return 271 }
235
236 // ============================================================
237 // Section M -- del_call (G2.3d) reuses indel_call shape.
238 // 8 reads, 7 with 3bp deletion -> 87% -> HOM.
239 // ============================================================
240
241 let dhist: *i64 = sys_mmap(128) as *i64
242 dhist[0] = 1
243 dhist[3] = 7
244 if del_call(dhist, 5, 8, 5, 20, 80, elen, ecnt) != NX_VARIANT_HOM { return 280 }
245 if elen[0] != 3 { return 281 }
246 if ecnt[0] != 7 { return 282 }
247
248 // No deletions -> NONE.
249 let dhist2: *i64 = sys_mmap(128) as *i64
250 dhist2[0] = 10
251 if del_call(dhist2, 5, 10, 5, 20, 80, elen, ecnt) != NX_VARIANT_NONE { return 283 }
252
253 // ============================================================
254 // Section N -- snv_call_multi (G2.4c).
255 // ============================================================
256
257 let m_alts: *i64 = sys_mmap(64) as *i64
258 let m_cnts: *i64 = sys_mmap(64) as *i64
259
260 // Case 1: tri-allelic. ref=A, counts[C]=5, counts[G]=3, counts[T]=2.
261 // Total ACGT = 10. All three non-ref have freq >= 20%. Emit
262 // in descending count order: [C=5, G=3, T=2].
263 let m_c1: *i64 = sys_mmap(64) as *i64
264 m_c1[NX_DNA_C] = 5
265 m_c1[NX_DNA_G] = 3
266 m_c1[NX_DNA_T] = 2
267 if snv_call_multi(NX_DNA_A, m_c1, 20, m_alts, m_cnts, 4) != 3 { return 300 }
268 if m_alts[0] != NX_DNA_C { return 301 }
269 if m_cnts[0] != 5 { return 302 }
270 if m_alts[1] != NX_DNA_G { return 303 }
271 if m_cnts[1] != 3 { return 304 }
272 if m_alts[2] != NX_DNA_T { return 305 }
273 if m_cnts[2] != 2 { return 306 }
274
275 // Case 2: bi-allelic at threshold. 8 C + 2 G out of 10.
276 // C: 80%, G: 20% (meets >=20% threshold). Emit both.
277 let m_c2: *i64 = sys_mmap(64) as *i64
278 m_c2[NX_DNA_C] = 8
279 m_c2[NX_DNA_G] = 2
280 if snv_call_multi(NX_DNA_A, m_c2, 20, m_alts, m_cnts, 4) != 2 { return 310 }
281 if m_alts[0] != NX_DNA_C { return 311 }
282 if m_cnts[0] != 8 { return 312 }
283 if m_alts[1] != NX_DNA_G { return 313 }
284 if m_cnts[1] != 2 { return 314 }
285
286 // Case 3: bi-allelic below threshold. 9 C + 1 G out of 10.
287 // C: 90%, G: 10% (below 20%). Emit only C.
288 let m_c3: *i64 = sys_mmap(64) as *i64
289 m_c3[NX_DNA_C] = 9
290 m_c3[NX_DNA_G] = 1
291 if snv_call_multi(NX_DNA_A, m_c3, 20, m_alts, m_cnts, 4) != 1 { return 320 }
292 if m_alts[0] != NX_DNA_C { return 321 }
293
294 // Case 4: all-ref. Returns 0 alts.
295 let m_c4: *i64 = sys_mmap(64) as *i64
296 m_c4[NX_DNA_A] = 10
297 if snv_call_multi(NX_DNA_A, m_c4, 20, m_alts, m_cnts, 4) != 0 { return 330 }
298
299 // Case 5: empty (no coverage at all).
300 let m_c5: *i64 = sys_mmap(64) as *i64
301 if snv_call_multi(NX_DNA_A, m_c5, 20, m_alts, m_cnts, 4) != 0 { return 340 }
302
303 // Case 6: capacity cap. Tri-allelic but max_alts=2 caps emission.
304 if snv_call_multi(NX_DNA_A, m_c1, 20, m_alts, m_cnts, 2) != 2 { return 350 }
305 if m_alts[0] != NX_DNA_C { return 351 }
306 if m_alts[1] != NX_DNA_G { return 352 }
307
308 // Case 7: tie-break -- ref=A, counts[C]=3, counts[G]=3. Equal.
309 // Leftmost-tie-break (smaller code first) -> [C, G].
310 let m_c7: *i64 = sys_mmap(64) as *i64
311 m_c7[NX_DNA_C] = 3
312 m_c7[NX_DNA_G] = 3
313 if snv_call_multi(NX_DNA_A, m_c7, 20, m_alts, m_cnts, 4) != 2 { return 360 }
314 if m_alts[0] != NX_DNA_C { return 361 }
315 if m_alts[1] != NX_DNA_G { return 362 }
316
317 // Case 8: bad input.
318 if snv_call_multi(-1, m_c1, 20, m_alts, m_cnts, 4) != -1 { return 370 }
319 if snv_call_multi(4, m_c1, 20, m_alts, m_cnts, 4) != -1 { return 371 }
320 if snv_call_multi(NX_DNA_A, m_c1, 20, m_alts, m_cnts, 0) != -1 { return 372 }
321
322 return 0
323}