nx_align_chain_stranded.nx source

1// nx_align_chain_stranded.nx -- strand-partitioned co-linear chain DP. 2// 3// license_tier: INDEPENDENT_REDERIVE 4// genealogy_id: international-research-sources/li-2018-minimap2-section-2.2-strand-aware 5// 6// G1.3d of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. Closes the 7// strand-aware-pipeline gap: nx_align_chain (G1.3) is strand-blind 8// and would happily chain forward-strand seeds together with 9// reverse-strand seeds, producing geometrically incoherent 10// alignments. This primitive partitions seeds by strand tag, 11// runs the strand-blind chain DP on each partition, and returns 12// the LONGER of the two chains plus its strand orientation. 13// 14// Why partition not unified-DP: 15// Forward and reverse alignments live in DIFFERENT coordinate 16// systems on the query side (the reverse case uses revcomp-query 17// coordinates). A single DP that mixed them would chain across 18// strand boundaries and produce alignments that don't correspond 19// to any real read placement. Partition + per-strand DP is the 20// simplest correct primitive. 21// 22// Input contract: 23// Caller provides pre-sorted (by r_pos asc) seed-pair arrays 24// that come ALREADY PARTITIONED into a forward block then a 25// reverse block, with the partition boundary count given by 26// n_fwd (so seeds [0..n_fwd) are forward, [n_fwd..n_total) are 27// reverse). This matches the output layout of 28// nx_align_match_stranded after a per-block nx_seed_sort. 29// 30// Tie-break: 31// When forward and reverse produce equal-length chains, the 32// FORWARD chain wins. Deterministic and matches the "default 33// to + strand on tie" convention used by minimap2. 34// 35// API: 36// chain_pick_best_strand( 37// q_arr, r_arr, n_fwd, n_total, 38// out_chain_idx, out_strand, 39// max_out 40// ) -> i64 length of winning chain (>= 0); -1 on bad input 41// 42// - q_arr / r_arr : combined fwd+rev seed arrays, contiguous, 43// each partition pre-sorted by r ascending 44// - n_fwd : count of forward seeds at the head 45// - n_total : total seed count (so rev count = n_total - n_fwd) 46// - out_chain_idx : indices INTO the combined arrays (not into 47// per-partition arrays); for the rev winner 48// these indices are >= n_fwd 49// - out_strand : single i64 written; NX_STRAND_FWD or 50// NX_STRAND_REV, or NX_STRAND_UNKNOWN if both 51// partitions are empty 52// 53// What G1.3d does NOT do (deferred): 54// - Affine-cost chain (minimap2 alpha + beta gap penalty) -- G1.5 55// - Secondary chain reporting (both strands when both meaningful) -- G1.4c 56// - Chain joining across small inversions -- G1.7 57// 58// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml) 59// intended_use: "Best-of-both-strands chain selection; 60// the orchestration primitive that turns a 61// strand-tagged seed stream into a single 62// coherent alignment hypothesis" 63// sil_target: SIL2 64// asil_target: QM 65// dal_target: DAL C 66// iec_62304_class: B 67// evidence: [no_floating_point, deterministic, 68// bit_equal_reproducible, 69// composes_nx_align_chain_KAT, 70// fwd_wins_tie_stable, 71// partition_boundary_explicit_no_silent_mix, 72// license_tier_INDEPENDENT_REDERIVE] 73// hazard_register: [bug-tape-chain-stranded-mixes-strands-silently, 74// bug-tape-rev-indices-not-offset-by-n_fwd, 75// bug-tape-tie-break-reversed-from-spec, 76// bug-tape-empty-partition-junk-out_strand] 77// residual_risk: "Discards the runner-up strand entirely; 78// for reads with similar fwd + rev evidence 79// (chimeras / inverted repeats) downstream 80// consumers should request both via the 81// G1.4c multi-chain primitive." 82// verdict: NOT_YET_EVALUATED 83 84import "nx_syscalls.nx" 85import "nx_const.nx" 86import "nx_align_chain.nx" 87 88// Pick the longer co-linear chain across forward and reverse partitions. 89// Returns winning chain length (>=0); writes indices INTO the combined 90// q_arr/r_arr to out_chain_idx and winning strand to out_strand[0]. 91func chain_pick_best_strand(q_arr: *i64, r_arr: *i64, 92 n_fwd: i64, n_total: i64, 93 out_chain_idx: *i64, 94 out_strand: *i64, 95 max_out: i64) -> i64 { 96 if n_total < 0 { return -1 } 97 if n_fwd < 0 { return -1 } 98 if n_fwd > n_total { return -1 } 99 if max_out <= 0 { return -1 } 100 101 let n_rev: i64 = n_total - n_fwd 102 103 // Scratch buffers for each partition's chain output. Allocate 104 // for the worst case = full partition length. 105 let fwd_buf: *i64 = sys_mmap((n_fwd + 8) * 8) as *i64 106 let rev_buf: *i64 = sys_mmap((n_rev + 8) * 8) as *i64 107 108 // Run chain DP on the forward partition (q_arr+0 .. n_fwd). 109 var fwd_len: i64 = 0 110 if n_fwd > 0 { 111 fwd_len = seed_chain(q_arr, r_arr, n_fwd, fwd_buf, n_fwd + 8) 112 if fwd_len < 0 { return -1 } 113 } 114 115 // Run chain DP on the reverse partition. Pointer arithmetic via 116 // cast-int-add-cast (pattern from nx_arith.nx) to get the offset 117 // sub-array pointer; rev seeds live at q_arr[n_fwd..n_total). 118 var rev_len: i64 = 0 119 if n_rev > 0 { 120 let off_bytes: i64 = n_fwd * 8 121 let rev_q: *i64 = ((q_arr as i64) + off_bytes) as *i64 122 let rev_r: *i64 = ((r_arr as i64) + off_bytes) as *i64 123 rev_len = seed_chain(rev_q, rev_r, n_rev, rev_buf, n_rev + 8) 124 if rev_len < 0 { return -1 } 125 } 126 127 // Both partitions empty -> no alignment. 128 if fwd_len == 0 { 129 if rev_len == 0 { 130 out_strand[0] = NX_STRAND_UNKNOWN 131 return 0 132 } 133 } 134 135 // Tie-break: forward wins on equal length. 136 var pick_strand: i64 = NX_STRAND_FWD 137 var pick_len: i64 = fwd_len 138 if rev_len > fwd_len { 139 pick_strand = NX_STRAND_REV 140 pick_len = rev_len 141 } 142 143 if pick_len > max_out { return -1 } 144 145 // Copy winning chain into out_chain_idx. For reverse winner, 146 // offset indices by n_fwd so they index into the combined array. 147 var i: i64 = 0 148 if pick_strand == NX_STRAND_FWD { 149 while i < pick_len { 150 out_chain_idx[i] = fwd_buf[i] 151 i = i + 1 152 } 153 } else { 154 while i < pick_len { 155 out_chain_idx[i] = rev_buf[i] + n_fwd 156 i = i + 1 157 } 158 } 159 160 out_strand[0] = pick_strand 161 return pick_len 162}