nx_align_match_stranded_test.nx source
↩ module page · 173 lines · 6190 B
1// nx_align_match_stranded_test.nx -- KAT for strand-aware seed matcher.
2//
3// End-to-end fixture demonstrating that canonical k-mer minimizers
4// have STRAND-INVARIANT VALUES (same set for fwd and rc query) but
5// the strand-aware matcher correctly tags pairs by which orientation
6// the match came from.
7//
8// Query : "AACGT" (5 bases)
9// Query rc : "ACGTT" (5 bases)
10// Reference : "AACGTACGTT" (10 bases) -- contains AACGT at pos 0 + ACGTT at pos 5
11//
12// Hand-traced canonical k-mers (k=4):
13// Fwd query:
14// pos 0: AACG raw 0x06 rc CGTT 0x6F canon 0x06
15// pos 1: ACGT palindrome 0x1B canon 0x1B
16// Rc query:
17// pos 0: ACGT palindrome 0x1B canon 0x1B
18// pos 1: CGTT raw 0x6F rc AACG 0x06 canon 0x06
19//
20// w=1 minimizers (every k-mer is its own window with w=1):
21// Fwd query: (0x06, 0), (0x1B, 1) count 2
22// Rc query: (0x1B, 0), (0x06, 1) count 2
23// (same value set, different positions)
24//
25// Reference k-mers (k=4):
26// pos 0: AACG canon 0x06
27// pos 1: ACGT canon 0x1B
28// pos 2: CGTA raw 0x6C rc TACG 0xC6 canon 0x6C
29// pos 3: GTAC palindrome 0xB1
30// pos 4: TACG raw 0xC6 rc CGTA 0x6C canon 0x6C
31// pos 5: ACGT canon 0x1B
32// pos 6: CGTT raw 0x6F rc AACG 0x06 canon 0x06
33//
34// w=1 ref minimizers (7):
35// (0x06, 0), (0x1B, 1), (0x6C, 2), (0xB1, 3), (0x6C, 4), (0x1B, 5), (0x06, 6)
36//
37// Strand-aware match output:
38// Forward block (tag = +1):
39// q=0 (val 0x06): hits ref 0x06 at positions 0, 6 -> (0,0)+, (0,6)+
40// q=1 (val 0x1B): hits ref 0x1B at positions 1, 5 -> (1,1)+, (1,5)+
41// Reverse block (tag = -1):
42// q=0 (val 0x1B): hits ref 0x1B at positions 1, 5 -> (0,1)-, (0,5)-
43// q=1 (val 0x06): hits ref 0x06 at positions 0, 6 -> (1,0)-, (1,6)-
44// Total: 8 pairs.
45//
46// expect_exit: 0
47//
48// license_tier: ORIGINAL
49
50import "nx_syscalls.nx"
51import "nx_const.nx"
52import "nx_sequence.nx"
53import "nx_align_minimizer.nx"
54import "nx_align_match_stranded.nx"
55
56func main() -> i64 {
57
58 // ============================================================
59 // Stage 1 -- pack reference + forward query + revcomp query.
60 // ============================================================
61
62 // Reference "AACGTACGTT" (10 bases).
63 // byte 0: AACG = 00 00 01 10 = 0x06
64 // byte 1: TACG = 11 00 01 10 = 0xC6
65 // byte 2: TT = 11 11 (low 4 bits padding zero) = 0xF0
66 let r_bases: *u8 = sys_mmap(8)
67 r_bases[0] = 0x06
68 r_bases[1] = 0xC6
69 r_bases[2] = 0xF0
70
71 // Forward query "AACGT" (5 bases).
72 // byte 0: AACG = 0x06
73 // byte 1: T = 11 (high 2 bits) padding = 0xC0
74 let q_fwd_bases: *u8 = sys_mmap(4)
75 q_fwd_bases[0] = 0x06
76 q_fwd_bases[1] = 0xC0
77
78 // Revcomp query "ACGTT" (5 bases).
79 // byte 0: ACGT = 0x1B
80 // byte 1: T = 11 (high 2 bits) padding = 0xC0
81 // Verify by calling dna_revcomp on q_fwd_bases:
82 let q_rc_bases: *u8 = sys_mmap(4)
83 dna_revcomp(q_fwd_bases, 5, q_rc_bases)
84 if (q_rc_bases[0] & 0xff) != 0x1B { return 1 }
85 if (q_rc_bases[1] & 0xff) != 0xC0 { return 2 }
86
87 // ============================================================
88 // Stage 2 -- extract minimizers from all three sequences (k=4 w=1).
89 // ============================================================
90
91 let r_v: *i64 = sys_mmap(128) as *i64
92 let r_p: *i64 = sys_mmap(128) as *i64
93 let r_n: i64 = minimizer_extract(r_bases, 10, 4, 1, r_v, r_p, 16)
94 if r_n != 7 { return 10 }
95
96 let q_fwd_v: *i64 = sys_mmap(64) as *i64
97 let q_fwd_p: *i64 = sys_mmap(64) as *i64
98 let q_fwd_n: i64 = minimizer_extract(q_fwd_bases, 5, 4, 1, q_fwd_v, q_fwd_p, 16)
99 if q_fwd_n != 2 { return 11 }
100 if q_fwd_v[0] != 0x06 { return 12 }
101 if q_fwd_p[0] != 0 { return 13 }
102 if q_fwd_v[1] != 0x1B { return 14 }
103 if q_fwd_p[1] != 1 { return 15 }
104
105 let q_rc_v: *i64 = sys_mmap(64) as *i64
106 let q_rc_p: *i64 = sys_mmap(64) as *i64
107 let q_rc_n: i64 = minimizer_extract(q_rc_bases, 5, 4, 1, q_rc_v, q_rc_p, 16)
108 if q_rc_n != 2 { return 16 }
109 if q_rc_v[0] != 0x1B { return 17 }
110 if q_rc_p[0] != 0 { return 18 }
111 if q_rc_v[1] != 0x06 { return 19 }
112 if q_rc_p[1] != 1 { return 20 }
113
114 // ============================================================
115 // Stage 3 -- strand-aware match.
116 // ============================================================
117
118 let oq: *i64 = sys_mmap(256) as *i64
119 let or: *i64 = sys_mmap(256) as *i64
120 let os: *i64 = sys_mmap(256) as *i64
121
122 let total: i64 = minimizer_match_stranded(
123 q_fwd_v, q_fwd_p, q_fwd_n,
124 q_rc_v, q_rc_p, q_rc_n,
125 r_v, r_p, r_n,
126 oq, or, os, 32)
127 if total != 8 { return 30 }
128
129 // Forward pairs in order (outer q, inner r):
130 // q=0 (0x06) hits ref 0x06 at r_pos 0 and 6
131 // q=1 (0x1B) hits ref 0x1B at r_pos 1 and 5
132 if oq[0] != 0 { return 40 }
133 if or[0] != 0 { return 41 }
134 if os[0] != NX_STRAND_FWD { return 42 }
135 if oq[1] != 0 { return 43 }
136 if or[1] != 6 { return 44 }
137 if os[1] != NX_STRAND_FWD { return 45 }
138 if oq[2] != 1 { return 46 }
139 if or[2] != 1 { return 47 }
140 if os[2] != NX_STRAND_FWD { return 48 }
141 if oq[3] != 1 { return 49 }
142 if or[3] != 5 { return 50 }
143 if os[3] != NX_STRAND_FWD { return 51 }
144
145 // Reverse pairs in order:
146 // q=0 (0x1B) hits ref 0x1B at r_pos 1 and 5
147 // q=1 (0x06) hits ref 0x06 at r_pos 0 and 6
148 if oq[4] != 0 { return 60 }
149 if or[4] != 1 { return 61 }
150 if os[4] != NX_STRAND_REV { return 62 }
151 if oq[5] != 0 { return 63 }
152 if or[5] != 5 { return 64 }
153 if os[5] != NX_STRAND_REV { return 65 }
154 if oq[6] != 1 { return 66 }
155 if or[6] != 0 { return 67 }
156 if os[6] != NX_STRAND_REV { return 68 }
157 if oq[7] != 1 { return 69 }
158 if or[7] != 6 { return 70 }
159 if os[7] != NX_STRAND_REV { return 71 }
160
161 // ============================================================
162 // Stage 4 -- empty-input cases.
163 // ============================================================
164
165 let t2: i64 = minimizer_match_stranded(
166 q_fwd_v, q_fwd_p, 0,
167 q_rc_v, q_rc_p, 0,
168 r_v, r_p, r_n,
169 oq, or, os, 32)
170 if t2 != 0 { return 80 }
171
172 return 0
173}