nx_align_minimizer.nx source
↩ module page · 154 lines · 6008 B
1// nx_align_minimizer.nx -- canonical k-mer minimizer extraction.
2//
3// license_tier: INDEPENDENT_REDERIVE
4// genealogy_id: international-research-sources/roberts-hunt-2004-schleimer-2003
5//
6// G1.1 of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. Composes the existing
7// nx_sequence k-mer + canonical-k-mer primitives.
8//
9// Roberts/Schleimer minimizer scheme:
10// - Slide a window of w consecutive k-mers along the sequence.
11// - For each window, the minimizer is the leftmost canonical k-mer
12// with the smallest value.
13// - As the window slides, emit a new (value, position) pair only
14// when it differs from the previous window's pair. This
15// deduplicates the common case of one minimizer persisting
16// across many overlapping windows.
17//
18// Why this scheme:
19// - Any two sequences sharing a substring of length k + w - 1 share
20// at least one minimizer (the "window guarantee"). This makes
21// minimizers a correct seeding primitive for read alignment.
22// - Sub-samples ~2/(w+1) of the k-mers, vastly reducing index size
23// vs. storing every k-mer (minimap2: ~0.2x reference size with
24// k=15 w=10 vs ~6x with all-kmers).
25// - Composable: same scheme for indexing a reference + for
26// extracting query seeds; matching is just (value, refpos) lookup.
27//
28// Tie-break:
29// When multiple k-mers in a window have the same min value, the
30// LEFTMOST one wins. This matches minimap2's convention and is
31// stable across hosts.
32//
33// Naive complexity:
34// O((n - k - w + 2) * w) -- one full scan per window. Fine for
35// the G1.1 reference impl + KAT. The deque-based O(n) variant
36// lives in nx_align_minimizer_fast.nx (G1.5).
37//
38// API:
39// minimizer_count_windows(n, k, w) -> i64
40// minimizer_extract(bases, n, k, w,
41// out_vals, out_pos, max_out) -> i64 count emitted
42//
43// bases: packed 2-bit DNA per nx_sequence encoding
44// n: base count
45// k: 1..32
46// w: >= 1
47// Pre: n >= k + w - 1
48//
49// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml)
50// intended_use: "Read alignment seeding -- composes with
51// nx_align Smith-Waterman extension to form
52// seed-and-extend pipeline (G1.2 + G1.3)"
53// sil_target: SIL2
54// asil_target: QM
55// dal_target: DAL C
56// iec_62304_class: B
57// evidence: [no_floating_point, deterministic,
58// bit_equal_reproducible,
59// roberts_2004_window_guarantee,
60// canonical_kmer_strand_invariant,
61// leftmost_tie_break_stable,
62// license_tier_INDEPENDENT_REDERIVE]
63// hazard_register: [bug-tape-minimizer-tie-break-rightmost-drift,
64// bug-tape-minimizer-window-off-by-one,
65// bug-tape-minimizer-dedup-misses-true-change,
66// bug-tape-minimizer-window-guarantee-violated]
67// residual_risk: "G1.1 naive O((n-k-w+2)*w); for reference
68// index over a 3 Gbp chromosome use the
69// deque-based O(n) impl in G1.5. Window
70// guarantee depends on canonical k-mer
71// correctness; nx_sequence k-mer KAT (G0.1)
72// verifies that."
73// verdict: NOT_YET_EVALUATED
74
75import "nx_syscalls.nx"
76import "nx_sequence.nx"
77
78// Number of windows of w consecutive k-mers in a sequence of n bases.
79// Returns 0 if the sequence is too short to fit even one window.
80func minimizer_count_windows(n: i64, k: i64, w: i64) -> i64 {
81 let num_kmers: i64 = n - k + 1
82 let num_windows: i64 = num_kmers - w + 1
83 if num_windows < 0 { return 0 }
84 return num_windows
85}
86
87// Extract minimizers. Returns count written to out_vals + out_pos,
88// or -1 on bad parameters / capacity exceeded.
89//
90// out_vals[i] = canonical k-mer value of the i-th emitted minimizer
91// out_pos[i] = base position (0-indexed) where that k-mer starts
92//
93// Deduplication: emit (val, pos) for window 0 unconditionally; for
94// each subsequent window, emit only if its (val, pos) differs from
95// the previous window's (val, pos).
96func minimizer_extract(bases: *u8, n: i64, k: i64, w: i64,
97 out_vals: *i64, out_pos: *i64,
98 max_out: i64) -> i64 {
99 if k <= 0 { return -1 }
100 if k > 32 { return -1 }
101 if w <= 0 { return -1 }
102 if max_out <= 0 { return -1 }
103
104 let num_windows: i64 = minimizer_count_windows(n, k, w)
105 if num_windows <= 0 { return 0 }
106
107 var emitted: i64 = 0
108 var prev_val: i64 = -1
109 var prev_pos: i64 = -1
110
111 var win: i64 = 0
112 while win < num_windows {
113 // Scan k-mers at positions win, win+1, ..., win+w-1.
114 // Find the leftmost minimum canonical value.
115 var min_val: i64 = -1
116 var min_pos: i64 = -1
117 var j: i64 = 0
118 while j < w {
119 let kmer_pos: i64 = win + j
120 let raw_kmer: i64 = dna_kmer_at(bases, kmer_pos, k)
121 let canon: i64 = dna_kmer_canonical(raw_kmer, k)
122 if j == 0 {
123 min_val = canon
124 min_pos = kmer_pos
125 } else {
126 if canon < min_val {
127 min_val = canon
128 min_pos = kmer_pos
129 }
130 }
131 j = j + 1
132 }
133
134 // Emit if differs from previous window's minimizer.
135 var emit: i64 = 0
136 if win == 0 {
137 emit = 1
138 } else {
139 if min_val != prev_val { emit = 1 }
140 if min_pos != prev_pos { emit = 1 }
141 }
142 if emit == 1 {
143 if emitted >= max_out { return -1 }
144 out_vals[emitted] = min_val
145 out_pos[emitted] = min_pos
146 emitted = emitted + 1
147 prev_val = min_val
148 prev_pos = min_pos
149 }
150 win = win + 1
151 }
152
153 return emitted
154}