nx_bam.nx source

1// nx_bam.nx -- BAM binary record serialiser (SAMv1 §4). 2// 3// G6.0d.1: nx_bam_pack_seq -- 4-bit IUPAC nucleotide packing. 4// Each base in the ASCII input is mapped to a 4-bit code per 5// SAMv1 §4.2 Table 1 (= 0, A=1, C=2, M=3, G=4, R=5, S=6, 6// V=7, T=8, W=9, Y=10, H=11, K=12, D=13, B=14, N=15). 7// Bases are packed two-per-byte, HIGH nibble first. 8// Output length: ceil(seq_len / 2) bytes. 9// On odd seq_len, last low nibble is zero per spec. 10// Unrecognised ASCII -> NX_BAM_NT_N (15) -- conservative. 11// 12// G6.0d.2: nx_bam_pack_cigar -- pack CIGAR ops into u32 LE words. 13// Each op packs as (op_len << 4) | op_code; written little-endian. 14// Input: parallel arrays ops[] + lens[] of n_op entries. 15// Output: n_op * 4 bytes. 16// 17// G6.0d.3: nx_bam_write_record -- canonical BAM alignment record. 18// block_size + 32-byte fixed header + read_name + cigar + seq + qual 19// per SAMv1 §4.2. Header fields written little-endian via nx_le. 20// 21// expect_exit: 0 22// 23// license_tier: ORIGINAL 24 25import "nx_syscalls.nx" 26import "nx_const.nx" 27import "nx_le.nx" 28const NX_MAGIC_65535: i64 = 65535 29 30// ASCII -> 4-bit BAM code. Case-insensitive (upper + lower mapped). 31// Anything not recognised -> NX_BAM_NT_N (conservative; preserves 32// downstream usability rather than failing the encode). 33func nx_bam_iupac_to_nibble(c: i64) -> i64 { 34 let ch: i64 = c & 0xff 35 if ch == 0x3D { return NX_BAM_NT_EQ } // = 36 if ch == 0x41 { return NX_BAM_NT_A } // A 37 if ch == 0x61 { return NX_BAM_NT_A } // a 38 if ch == 0x43 { return NX_BAM_NT_C } // C 39 if ch == 0x63 { return NX_BAM_NT_C } // c 40 if ch == 0x4D { return NX_BAM_NT_M } // M 41 if ch == 0x6D { return NX_BAM_NT_M } // m 42 if ch == 0x47 { return NX_BAM_NT_G } // G 43 if ch == 0x67 { return NX_BAM_NT_G } // g 44 if ch == 0x52 { return NX_BAM_NT_R } // R 45 if ch == 0x72 { return NX_BAM_NT_R } // r 46 if ch == 0x53 { return NX_BAM_NT_S } // S 47 if ch == 0x73 { return NX_BAM_NT_S } // s 48 if ch == 0x56 { return NX_BAM_NT_V } // V 49 if ch == 0x76 { return NX_BAM_NT_V } // v 50 if ch == 0x54 { return NX_BAM_NT_T } // T 51 if ch == 0x74 { return NX_BAM_NT_T } // t 52 if ch == 0x55 { return NX_BAM_NT_T } // U -> T (RNA collapsed to DNA) 53 if ch == 0x75 { return NX_BAM_NT_T } // u 54 if ch == 0x57 { return NX_BAM_NT_W } // W 55 if ch == 0x77 { return NX_BAM_NT_W } // w 56 if ch == 0x59 { return NX_BAM_NT_Y } // Y 57 if ch == 0x79 { return NX_BAM_NT_Y } // y 58 if ch == 0x48 { return NX_BAM_NT_H } // H 59 if ch == 0x68 { return NX_BAM_NT_H } // h 60 if ch == 0x4B { return NX_BAM_NT_K } // K 61 if ch == 0x6B { return NX_BAM_NT_K } // k 62 if ch == 0x44 { return NX_BAM_NT_D } // D 63 if ch == 0x64 { return NX_BAM_NT_D } // d 64 if ch == 0x42 { return NX_BAM_NT_B } // B 65 if ch == 0x62 { return NX_BAM_NT_B } // b 66 return NX_BAM_NT_N // N or any unrecognised 67} 68 69// G6.0d.1 -- pack ASCII IUPAC sequence into BAM 4-bit nibbles. 70// Returns number of bytes written, or -1 on error. 71func nx_bam_pack_seq(seq_ascii: *u8, seq_len: i64, 72 out_bytes: *u8, max_out: i64) -> i64 { 73 if seq_len < 0 { return -1 } 74 if max_out < 0 { return -1 } 75 76 let n_bytes: i64 = (seq_len + 1) / 2 77 if n_bytes > max_out { return -1 } 78 79 var i: i64 = 0 80 while i < n_bytes { 81 let hi_idx: i64 = i * 2 82 let lo_idx: i64 = hi_idx + 1 83 let hi: i64 = nx_bam_iupac_to_nibble(seq_ascii[hi_idx] as i64) 84 var lo: i64 = 0 85 if lo_idx < seq_len { 86 lo = nx_bam_iupac_to_nibble(seq_ascii[lo_idx] as i64) 87 } 88 let byte_val: i64 = ((hi & 0xf) << 4) | (lo & 0xf) 89 out_bytes[i] = byte_val as u8 90 i = i + 1 91 } 92 93 return n_bytes 94} 95 96// Field-array indices for nx_bam_write_record's bundled args. 97// Caller passes a 12-entry i64[] in this order (3 length fields 98// folded in to keep the function-arg count at the nxc2 codegen 99// safe budget): 100const NX_BAM_REC_REF_ID: i64 = 0 101const NX_BAM_REC_POS: i64 = 1 // 0-indexed (NOT 1-indexed like SAM) 102const NX_BAM_REC_MAPQ: i64 = 2 103const NX_BAM_REC_BIN: i64 = 3 // UCSC bin (0 for unbinned) 104const NX_BAM_REC_FLAG: i64 = 4 105const NX_BAM_REC_NEXT_REF: i64 = 5 106const NX_BAM_REC_NEXT_POS: i64 = 6 // 0-indexed 107const NX_BAM_REC_TLEN: i64 = 7 // signed 108const NX_BAM_REC_NO_QUAL: i64 = 8 // 0 = qual present, 1 = fill 0xFF 109const NX_BAM_REC_READ_NAME_LEN: i64 = 9 // bytes (NUL appended internally) 110const NX_BAM_REC_N_CIGAR_OP: i64 = 10 111const NX_BAM_REC_L_SEQ: i64 = 11 112 113// G6.0d.2 -- pack CIGAR ops into BAM u32 LE words. 114// Each op encodes as (op_len << 4) | op_code. 115// Returns total bytes written (n_op * 4), or -1 on error. 116func nx_bam_pack_cigar(ops: *i64, lens: *i64, n_op: i64, 117 out_bytes: *u8, max_out: i64) -> i64 { 118 if n_op < 0 { return -1 } 119 if max_out < 0 { return -1 } 120 let n_bytes: i64 = n_op * 4 121 if n_bytes > max_out { return -1 } 122 123 var i: i64 = 0 124 while i < n_op { 125 let op: i64 = ops[i] 126 let lng: i64 = lens[i] 127 if op < 0 { return -1 } 128 if op > 8 { return -1 } // valid CIGAR op codes 0..8 per spec 129 if lng < 0 { return -1 } 130 let packed: i64 = (lng << 4) | (op & 0xf) 131 let off: i64 = i * 4 132 let r: i64 = nx_le_write_u32(out_bytes, off, packed) 133 if r < 0 { return -1 } 134 i = i + 1 135 } 136 137 return n_bytes 138} 139 140// G6.0d.3 -- write a full BAM alignment record per SAMv1 §4.2. 141// 142// Wire layout (every multi-byte field little-endian): 143// [0..4) block_size (i32) -- bytes 4..end (i.e., total - 4) 144// [4..8) refID (i32, -1 = *) 145// [8..12) pos (i32, 0-indexed; -1 = unmapped) 146// [12] l_read_name (u8, includes trailing NUL) 147// [13] mapq (u8) 148// [14..16) bin (u16; 0 acceptable if not binned) 149// [16..18) n_cigar_op (u16) 150// [18..20) flag (u16) 151// [20..24) l_seq (u32) 152// [24..28) next_refID (i32) 153// [28..32) next_pos (i32, 0-indexed) 154// [32..36) tlen (i32, signed) 155// [36..] read_name (l_read_name bytes, NUL-terminated) 156// + cigar (n_cigar_op * 4 bytes) 157// + seq ((l_seq+1)/2 bytes packed) 158// + qual (l_seq bytes; 0xFF * l_seq if NO_QUAL=1) 159// 160// Bundled-args API per the recurring nxc2 codegen quirk (8 args): 161// read_name_ascii -- WITHOUT trailing NUL (we append one) 162// cigar_ops/lens -- arrays of length fields[N_CIGAR_OP] 163// seq_ascii -- IUPAC string, length = fields[L_SEQ] 164// qual_phred -- raw Phred bytes (NOT ASCII +33), length = L_SEQ. 165// Ignored when fields[NX_BAM_REC_NO_QUAL] == 1. 166// fields -- 12-entry i64 array indexed by NX_BAM_REC_* 167// (carries the 3 length counts so the function 168// arg count stays at 8). 169// 170// Returns total bytes written (including the 4-byte block_size 171// prefix), or -1 on validation / capacity failure. 172func nx_bam_write_record(read_name_ascii: *u8, 173 cigar_ops: *i64, cigar_lens: *i64, 174 seq_ascii: *u8, 175 qual_phred: *u8, 176 fields: *i64, 177 out_bytes: *u8, max_out: i64) -> i64 { 178 let read_name_len: i64 = fields[NX_BAM_REC_READ_NAME_LEN] 179 let n_cigar_op: i64 = fields[NX_BAM_REC_N_CIGAR_OP] 180 let l_seq: i64 = fields[NX_BAM_REC_L_SEQ] 181 182 if read_name_len < 1 { return -1 } 183 if n_cigar_op < 0 { return -1 } 184 if l_seq < 0 { return -1 } 185 if max_out < 36 { return -1 } // even an empty record has 36-byte hdr 186 187 let l_read_name: i64 = read_name_len + 1 // BAM spec: includes NUL 188 if l_read_name > 255 { return -1 } // u8 cap 189 190 let cigar_bytes: i64 = n_cigar_op * 4 191 let seq_bytes: i64 = (l_seq + 1) / 2 192 let qual_bytes: i64 = l_seq 193 194 // total bytes including the leading 4-byte block_size field itself. 195 let total: i64 = 4 + 32 + l_read_name + cigar_bytes + seq_bytes + qual_bytes 196 if total > max_out { return -1 } 197 198 // block_size excludes the 4-byte block_size field itself. 199 let block_size: i64 = total - 4 200 201 // Reject obvious out-of-range u16 fields up front. 202 let bin: i64 = fields[NX_BAM_REC_BIN] 203 let flag: i64 = fields[NX_BAM_REC_FLAG] 204 if bin < 0 { return -1 } 205 if bin > NX_MAGIC_65535 { return -1 } 206 if flag < 0 { return -1 } 207 if flag > NX_MAGIC_65535 { return -1 } 208 if n_cigar_op > NX_MAGIC_65535 { return -1 } 209 210 let mapq: i64 = fields[NX_BAM_REC_MAPQ] 211 if mapq < 0 { return -1 } 212 if mapq > 255 { return -1 } 213 214 let no_qual: i64 = fields[NX_BAM_REC_NO_QUAL] 215 216 // ---- fixed-size header ---- 217 if nx_le_write_u32(out_bytes, 0, block_size) < 0 { return -1 } 218 if nx_le_write_u32(out_bytes, 4, fields[NX_BAM_REC_REF_ID]) < 0 { return -1 } 219 if nx_le_write_u32(out_bytes, 8, fields[NX_BAM_REC_POS]) < 0 { return -1 } 220 out_bytes[12] = l_read_name as u8 221 out_bytes[13] = mapq as u8 222 if nx_le_write_u16(out_bytes, 14, bin) < 0 { return -1 } 223 if nx_le_write_u16(out_bytes, 16, n_cigar_op) < 0 { return -1 } 224 if nx_le_write_u16(out_bytes, 18, flag) < 0 { return -1 } 225 if nx_le_write_u32(out_bytes, 20, l_seq) < 0 { return -1 } 226 if nx_le_write_u32(out_bytes, 24, fields[NX_BAM_REC_NEXT_REF]) < 0 { return -1 } 227 if nx_le_write_u32(out_bytes, 28, fields[NX_BAM_REC_NEXT_POS]) < 0 { return -1 } 228 if nx_le_write_u32(out_bytes, 32, fields[NX_BAM_REC_TLEN]) < 0 { return -1 } 229 230 // ---- variable sections ---- 231 var off: i64 = 36 232 233 // read_name (NUL-terminated) 234 var i: i64 = 0 235 while i < read_name_len { 236 out_bytes[off + i] = read_name_ascii[i] 237 i = i + 1 238 } 239 out_bytes[off + read_name_len] = 0 // trailing NUL 240 off = off + l_read_name 241 242 // cigar 243 if n_cigar_op > 0 { 244 let cig_dst: *u8 = ((out_bytes as i64) + off) as *u8 245 let cn: i64 = nx_bam_pack_cigar(cigar_ops, cigar_lens, n_cigar_op, 246 cig_dst, cigar_bytes) 247 if cn != cigar_bytes { return -1 } 248 off = off + cigar_bytes 249 } 250 251 // seq 252 if l_seq > 0 { 253 let seq_dst: *u8 = ((out_bytes as i64) + off) as *u8 254 let sn: i64 = nx_bam_pack_seq(seq_ascii, l_seq, seq_dst, seq_bytes) 255 if sn != seq_bytes { return -1 } 256 off = off + seq_bytes 257 } 258 259 // qual 260 if l_seq > 0 { 261 if no_qual == 1 { 262 var qi: i64 = 0 263 while qi < l_seq { 264 out_bytes[off + qi] = NX_BAM_NO_QUAL_BYTE as u8 265 qi = qi + 1 266 } 267 } else { 268 var qj: i64 = 0 269 while qj < l_seq { 270 out_bytes[off + qj] = qual_phred[qj] 271 qj = qj + 1 272 } 273 } 274 off = off + l_seq 275 } 276 277 return off // total bytes written; equals `total` above. 278}