nx_align_match_stranded.nx source
↩ module page · 144 lines · 6342 B
1// nx_align_match_stranded.nx -- strand-aware seed matcher.
2//
3// license_tier: INDEPENDENT_REDERIVE
4// genealogy_id: international-research-sources/li-2018-minimap2-section-2.1
5//
6// G1.3c of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. Real-world reads
7// arrive from either strand of the reference; a forward-only matcher
8// misses 50% of true alignments. This primitive matches BOTH the
9// forward-query minimizers and the reverse-complement-query
10// minimizers against the reference, tagging each emitted pair with
11// strand orientation.
12//
13// Strand semantics:
14// q_fwd_* -- minimizers from the query as given
15// q_rc_* -- minimizers from the reverse-complement of the query
16// (the caller pre-extracts these via dna_revcomp +
17// minimizer_extract)
18// Match emits forward-strand pairs first (tag = NX_STRAND_FWD = +1),
19// then reverse-strand pairs (tag = NX_STRAND_REV = -1).
20// Within each strand, emission follows minimizer_match's outer-q /
21// inner-r order so chain DP sees a deterministic stream.
22//
23// Why canonical-k-mer minimizers are NOT enough on their own:
24// Canonical k-mers are strand-invariant by construction (the
25// smaller of forward and revcomp wins), so the SET of canonical
26// minimizer VALUES is identical between a query and its revcomp.
27// But the POSITIONS differ: a query base at position i corresponds
28// to revcomp-query base at position (n - 1 - i). Strand info
29// lives in this position-mapping, not in the value. Therefore
30// we run TWO matches and tag the output stream.
31//
32// q_pos semantics in the output:
33// Forward pairs: q_pos refers to the original query's coordinate
34// (where the k-mer starts in the query as given).
35// Reverse pairs: q_pos refers to the REVCOMP-QUERY coordinate.
36// Downstream chain DP and SW-extension are responsible for
37// un-mapping reverse-strand q_pos back to original-query
38// coordinates when needed (see G1.4b nx_alignment_record).
39//
40// API:
41// minimizer_match_stranded(
42// q_fwd_v, q_fwd_p, q_fwd_n,
43// q_rc_v, q_rc_p, q_rc_n,
44// r_v, r_p, r_n,
45// out_q, out_r, out_strand,
46// max_out
47// ) -> i64 total pair count (>= 0); -1 on overflow / bad input
48//
49// Composes:
50// - nx_sequence::dna_revcomp -- caller pre-computes rc bases
51// - nx_align_minimizer -- caller runs twice (fwd + rc)
52// - nx_align_match -- this primitive runs twice internally
53// - nx_const::NX_STRAND_* -- canonical strand tag values
54//
55// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml)
56// intended_use: "Strand-aware seed match -- the real-world
57// version that handles reads from either
58// reference strand; required for any aligner
59// used on actual sequencer output"
60// sil_target: SIL2
61// asil_target: QM
62// dal_target: DAL C
63// iec_62304_class: B
64// evidence: [no_floating_point, deterministic,
65// bit_equal_reproducible,
66// composes_minimizer_match_KAT,
67// canonical_kmer_strand_invariance_documented,
68// strand_tag_via_nx_const_canonical,
69// license_tier_INDEPENDENT_REDERIVE]
70// hazard_register: [bug-tape-strand-tag-sign-swapped,
71// bug-tape-revcomp-coordinate-not-remapped,
72// bug-tape-fwd-and-rc-double-count-palindromes,
73// bug-tape-strand-stream-interleave-instead-of-concat]
74// residual_risk: "Palindromic k-mers (canonical == revcomp)
75// produce identical pairs on both strands;
76// downstream chain DP must dedupe or accept
77// the redundancy. Not addressed in G1.3c
78// because downstream tolerance is the more
79// composable solution."
80// verdict: NOT_YET_EVALUATED
81
82import "nx_syscalls.nx"
83import "nx_const.nx"
84import "nx_align_match.nx"
85
86// Strand-aware match. Internally runs minimizer_match twice and
87// concatenates outputs (fwd block first, rc block second), tagging
88// each emitted pair with its originating strand via out_strand[i].
89func minimizer_match_stranded(
90 q_fwd_v: *i64, q_fwd_p: *i64, q_fwd_n: i64,
91 q_rc_v: *i64, q_rc_p: *i64, q_rc_n: i64,
92 r_v: *i64, r_p: *i64, r_n: i64,
93 out_q: *i64, out_r: *i64, out_strand: *i64,
94 max_out: i64) -> i64 {
95 if q_fwd_n < 0 { return -1 }
96 if q_rc_n < 0 { return -1 }
97 if r_n < 0 { return -1 }
98 if max_out <= 0 { return -1 }
99
100 // ---- Forward block ----
101 let n_fwd: i64 = minimizer_match(q_fwd_v, q_fwd_p, q_fwd_n,
102 r_v, r_p, r_n,
103 out_q, out_r, max_out)
104 if n_fwd < 0 { return -1 }
105
106 // Tag forward pairs.
107 var i: i64 = 0
108 while i < n_fwd {
109 out_strand[i] = NX_STRAND_FWD
110 i = i + 1
111 }
112
113 // ---- Reverse block: append after forward, with offset pointers ----
114 let remaining: i64 = max_out - n_fwd
115 if remaining <= 0 {
116 // Forward filled the buffer; can only succeed if rc would emit zero.
117 // Probe rc with capacity 1 sentinel via a tiny stack array...
118 // Simpler: refuse if rc could emit anything. Tight bound check
119 // delegated to caller via documented residual_risk.
120 if q_rc_n == 0 { return n_fwd }
121 if r_n == 0 { return n_fwd }
122 return -1
123 }
124
125 // Pointer arithmetic via cast-int-add-cast (pattern from nx_arith.nx).
126 let offset_bytes: i64 = n_fwd * 8
127 let rc_out_q: *i64 = ((out_q as i64) + offset_bytes) as *i64
128 let rc_out_r: *i64 = ((out_r as i64) + offset_bytes) as *i64
129 let rc_out_strand: *i64 = ((out_strand as i64) + offset_bytes) as *i64
130
131 let n_rc: i64 = minimizer_match(q_rc_v, q_rc_p, q_rc_n,
132 r_v, r_p, r_n,
133 rc_out_q, rc_out_r, remaining)
134 if n_rc < 0 { return -1 }
135
136 // Tag reverse pairs.
137 var j: i64 = 0
138 while j < n_rc {
139 rc_out_strand[j] = NX_STRAND_REV
140 j = j + 1
141 }
142
143 return n_fwd + n_rc
144}