nx_alignment.nx source

1// nx_alignment.nx -- AlignmentRecord struct + mapq + revcomp coord remap. 2// 3// license_tier: INDEPENDENT_REDERIVE 4// genealogy_id: international-research-sources/li-2018-minimap2-mapq + sam-spec-v1 5// 6// G1.4b of NISHI_GENOMICS_SUBSTRATE_ROADMAP.md. The output structure 7// that downstream consumers (variant callers, BAM writers, alignment 8// browsers, audit tooling) consume. Closes the "what does an 9// aligner return?" gap by giving every read-vs-reference alignment 10// a uniform queryable record. 11// 12// Field semantics (SAM-spec aligned but Nishi-native): 13// query_id : caller-assigned numeric id for the query/read 14// ref_id : caller-assigned numeric id for the reference sequence 15// q_start : 0-indexed start position in the ORIGINAL query 16// (NOT the rc-query; rc alignments are remapped via 17// nx_remap_revcomp_coords before record construction) 18// q_end : 0-indexed end (exclusive) in original query 19// r_start : 0-indexed start in reference 20// r_end : 0-indexed end (exclusive) in reference 21// strand : NX_STRAND_FWD (+1) or NX_STRAND_REV (-1) from nx_const 22// score : SW alignment score (linear or affine, caller's choice) 23// mapq : mapping quality, 0..NX_MAPQ_MAX clamped 24// n_seeds : number of seeds in the chain that anchored this 25// alignment (passes through from nx_align_chain output) 26// 27// Mapping-quality model (G1.4b reference impl): 28// mapq = clamp(0, 60, floor(60 * (primary - secondary) / max(1, primary))) 29// - primary = score of the best alignment for this query 30// - secondary = score of the best ALTERNATIVE alignment, or 0 if unique 31// - secondary == 0 -> mapq = 60 (uniquely placed) 32// - secondary == primary -> mapq = 0 (cannot distinguish primary from a tie) 33// - secondary == primary/2 -> mapq = 30 (moderate confidence) 34// - secondary > primary -> mapq = 0 (negative clamps to 0) 35// - primary <= 0 -> mapq = 0 (no real alignment) 36// This is the SIMPLEST coherent integer-arithmetic mapq that 37// matches the high-end BWA-MEM 60-cap convention without 38// pretending to model minimap2's full probabilistic formula 39// (which uses chain count + log + sub-best heuristics, all of 40// which require nx_fixed_point primitives -- G3+). 41// 42// Reverse-strand coordinate remapping: 43// When a chain is found against the REVCOMP query, the chain's 44// (q_start, q_end) are in rc-query coordinates. To present a 45// record in original-query coordinates the substrate flips: 46// orig_q_start = q_len - rc_q_end 47// orig_q_end = q_len - rc_q_start 48// (rc-coordinate i corresponds to original-coordinate (q_len-1-i) 49// for a single base; intervals flip and re-bound.) 50// 51// What G1.4b does NOT do (deferred): 52// - CIGAR string output -- requires SW backtrace (G1.0c) and 53// soft-clipping logic; both queued as G1.4b.2 54// - SAM/BAM serialisation -- nx_sam_writer (G6) 55// - Probabilistic mapq from full minimap2 formula (needs log + 56// float-equivalent fixed-point math) -- G1.5 57// - Secondary alignment reporting (multi-hit reads) -- G1.4c 58// - Paired-end fragment-pair grouping -- G1.6 59// 60// API: 61// nx_alignment_new(...10 args bundled as fields[10] array...) -> *AlignmentRecord 62// nx_alignment_compute_mapq(primary_score, secondary_score) -> i64 0..60 63// nx_remap_revcomp_coords(q_len, rc_start, rc_end, 64// out_orig_start, out_orig_end) -> i64 0 65// 66// nx_safety_envelope: (schema: nishi-library/seeds/safety-critical-standards.toml) 67// intended_use: "Output record for read-to-reference 68// alignment; the structure downstream 69// variant-calling / BAM-writing / alignment- 70// audit primitives consume" 71// sil_target: SIL2 72// asil_target: QM 73// dal_target: DAL C 74// iec_62304_class: B 75// evidence: [no_floating_point, deterministic, 76// bit_equal_reproducible, 77// integer_arithmetic_mapq_clamped_0_60, 78// revcomp_coord_remap_involution_KAT, 79// composes_nx_const_strand_constants, 80// license_tier_INDEPENDENT_REDERIVE] 81// hazard_register: [bug-tape-mapq-negative-not-clamped, 82// bug-tape-mapq-divide-by-zero-primary, 83// bug-tape-revcomp-remap-off-by-one, 84// bug-tape-revcomp-remap-not-involutive, 85// bug-tape-record-strand-tag-out-of-domain] 86// residual_risk: "G1.4b uses a simplified mapq; comparison 87// with BWA-MEM mapq scores on identical 88// alignments will diverge. Documented in 89// the formula comment so consumers can pick 90// the heuristic that matches their tooling." 91// verdict: NOT_YET_EVALUATED 92 93import "nx_syscalls.nx" 94import "nx_const.nx" 95 96// Mapping-quality upper bound (BWA-MEM / minimap2 / Bowtie2 / SAM-spec convention). 97// Not in nx_const because "60 is the mapq ceiling" is an aligner-convention 98// choice (Bowtie1 used 255), not a universal mathematical constant. 99const NX_MAPQ_MAX: i64 = 60 100const NX_MAPQ_MIN: i64 = 0 101 102struct AlignmentRecord { 103 query_id: i64, 104 ref_id: i64, 105 q_start: i64, 106 q_end: i64, 107 r_start: i64, 108 r_end: i64, 109 strand: i64, 110 score: i64, 111 mapq: i64, 112 n_seeds: i64, 113} 114 115const NX_ALIGNMENT_RECORD_BYTES: i64 = 80 // 10 i64 fields 116 117// Build a fresh AlignmentRecord from a packed fields[] i64 array. 118// fields layout matches struct field order: 119// [0]=query_id [1]=ref_id [2]=q_start [3]=q_end [4]=r_start 120// [5]=r_end [6]=strand [7]=score [8]=mapq [9]=n_seeds 121// Caller is responsible for valid strand value (NX_STRAND_FWD or 122// NX_STRAND_REV). Returns pointer to the newly mmapped record. 123func nx_alignment_new(fields: *i64) -> *AlignmentRecord { 124 let rec: *AlignmentRecord = sys_mmap(NX_ALIGNMENT_RECORD_BYTES) as *AlignmentRecord 125 rec.query_id = fields[0] 126 rec.ref_id = fields[1] 127 rec.q_start = fields[2] 128 rec.q_end = fields[3] 129 rec.r_start = fields[4] 130 rec.r_end = fields[5] 131 rec.strand = fields[6] 132 rec.score = fields[7] 133 rec.mapq = fields[8] 134 rec.n_seeds = fields[9] 135 return rec 136} 137 138// Integer mapping quality clamped to [0, NX_MAPQ_MAX]. 139// mapq = clamp(0, 60, 60 * (primary - secondary) / max(1, primary)) 140// Primary <= 0 short-circuits to 0 (no real alignment). 141// Secondary >= primary -> 0 (cannot distinguish from a tie / better hit). 142func nx_alignment_compute_mapq(primary_score: i64, secondary_score: i64) -> i64 { 143 if primary_score <= 0 { return NX_MAPQ_MIN } 144 if secondary_score >= primary_score { return NX_MAPQ_MIN } 145 146 let diff: i64 = primary_score - secondary_score 147 let raw: i64 = (NX_MAPQ_MAX * diff) / primary_score 148 149 if raw < NX_MAPQ_MIN { return NX_MAPQ_MIN } 150 if raw > NX_MAPQ_MAX { return NX_MAPQ_MAX } 151 return raw 152} 153 154// Remap an interval from revcomp-query coordinates back to original-query 155// coordinates. rc-coord i corresponds to original-coord (q_len-1-i) for 156// a single base; for the half-open interval [rc_start, rc_end), the 157// remapped half-open interval is [q_len - rc_end, q_len - rc_start). 158// 159// Returns 0 on success; -1 if input range is out of bounds. 160// 161// Involution property: remap(remap(x)) == x. Verified in KAT. 162func nx_remap_revcomp_coords(q_len: i64, 163 rc_start: i64, rc_end: i64, 164 out_orig_start: *i64, 165 out_orig_end: *i64) -> i64 { 166 if q_len <= 0 { return -1 } 167 if rc_start < 0 { return -1 } 168 if rc_end <= rc_start { return -1 } 169 if rc_end > q_len { return -1 } 170 171 out_orig_start[0] = q_len - rc_end 172 out_orig_end[0] = q_len - rc_start 173 return 0 174}