Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
51 changes: 51 additions & 0 deletions crates/fgumi-raw-bam/src/cigar.rs
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,34 @@ pub fn get_cigar_ops(bam: &[u8]) -> Vec<u32> {
cigar_bytes.chunks_exact(4).map(|c| u32::from_le_bytes([c[0], c[1], c[2], c[3]])).collect()
}

/// CIGAR operation type to SAM character mapping.
const CIGAR_OP_CHARS: [u8; 9] = [b'M', b'I', b'D', b'N', b'S', b'H', b'P', b'=', b'X'];

/// Formats CIGAR operations from a raw BAM record as a SAM-style CIGAR string.
///
/// Returns an empty string if the record has no CIGAR operations, or `"*"` is not used
/// (callers should check `n_cigar_op == 0` separately if they need `"*"`).
///
/// # Panics
///
/// Panics if a CIGAR operation has an invalid type (>= 9).
#[must_use]
pub fn cigar_to_string_from_raw(bam: &[u8]) -> String {
let ops = get_cigar_ops(bam);
if ops.is_empty() {
return String::new();
}
let mut result = String::with_capacity(ops.len() * 4);
for &op in &ops {
let op_len = op >> 4;
let op_type = (op & 0xF) as usize;
assert!(op_type < CIGAR_OP_CHARS.len(), "invalid CIGAR op type: {op_type}");
result.push_str(&op_len.to_string());
result.push(CIGAR_OP_CHARS[op_type] as char);
}
result
}

/// Calculate reference-consuming length from CIGAR operations.
///
/// This is the sum of M/D/N/=/X operations, which represents how many
Expand Down Expand Up @@ -1934,4 +1962,27 @@ mod tests {
rec.truncate(36);
assert_eq!(unclipped_5prime_raw(&rec, 100, false), 100);
}

// ========================================================================
// cigar_to_string_from_raw tests
// ========================================================================

#[test]
fn test_cigar_to_string_simple() {
let rec = make_bam_bytes(0, 0, 0, b"r1", &[encode_op(0, 100)], 100, -1, -1, &[]);
assert_eq!(cigar_to_string_from_raw(&rec), "100M");
}

#[test]
fn test_cigar_to_string_complex() {
let ops = [encode_op(4, 5), encode_op(0, 50), encode_op(1, 3), encode_op(0, 42)];
let rec = make_bam_bytes(0, 0, 0, b"r1", &ops, 100, -1, -1, &[]);
assert_eq!(cigar_to_string_from_raw(&rec), "5S50M3I42M");
}

#[test]
fn test_cigar_to_string_no_cigar() {
let rec = make_bam_bytes(0, 0, 0, b"r1", &[], 4, -1, -1, &[]);
assert_eq!(cigar_to_string_from_raw(&rec), "");
}
}
88 changes: 88 additions & 0 deletions crates/fgumi-raw-bam/src/fields.rs
Original file line number Diff line number Diff line change
Expand Up @@ -193,6 +193,42 @@ pub fn set_flags(bam: &mut [u8], new_flags: u16) {
bam[14..16].copy_from_slice(&new_flags.to_le_bytes());
}

/// Set the mapping quality in a BAM record.
#[inline]
pub fn set_mapq(bam: &mut [u8], new_mapq: u8) {
bam[9] = new_mapq;
}

/// Set the mate reference sequence ID in a BAM record.
#[inline]
pub fn set_mate_ref_id(bam: &mut [u8], val: i32) {
bam[20..24].copy_from_slice(&val.to_le_bytes());
}

/// Set the mate 0-based position in a BAM record.
#[inline]
pub fn set_mate_pos(bam: &mut [u8], val: i32) {
bam[24..28].copy_from_slice(&val.to_le_bytes());
}

/// Set the template length (TLEN) in a BAM record.
#[inline]
pub fn set_template_length(bam: &mut [u8], val: i32) {
bam[28..32].copy_from_slice(&val.to_le_bytes());
}

/// Set the reference sequence ID in a BAM record.
#[inline]
pub fn set_ref_id(bam: &mut [u8], val: i32) {
bam[0..4].copy_from_slice(&val.to_le_bytes());
}

/// Set the 0-based leftmost position in a BAM record.
#[inline]
pub fn set_pos(bam: &mut [u8], val: i32) {
bam[4..8].copy_from_slice(&val.to_le_bytes());
}

/// Calculate the offset to auxiliary data in a BAM record.
///
/// `aux_offset = 32 + l_read_name + n_cigar_op*4 + (l_seq+1)/2 + l_seq`
Expand Down Expand Up @@ -841,4 +877,56 @@ mod tests {
assert!(fields.flags.paired());
assert_eq!(fields.name, b"rd");
}

// ========================================================================
// Field setter tests
// ========================================================================

#[test]
fn test_set_mapq_roundtrip() {
let mut rec = make_bam_bytes(0, 0, 0, b"r1", &[], 4, -1, -1, &[]);
assert_eq!(mapq(&rec), 0);
set_mapq(&mut rec, 42);
assert_eq!(mapq(&rec), 42);
}

#[test]
fn test_set_mate_ref_id_roundtrip() {
let mut rec = make_bam_bytes(0, 0, 0, b"r1", &[], 4, -1, -1, &[]);
assert_eq!(mate_ref_id(&rec), -1);
set_mate_ref_id(&mut rec, 5);
assert_eq!(mate_ref_id(&rec), 5);
}

#[test]
fn test_set_mate_pos_roundtrip() {
let mut rec = make_bam_bytes(0, 0, 0, b"r1", &[], 4, -1, -1, &[]);
assert_eq!(mate_pos(&rec), -1);
set_mate_pos(&mut rec, 12345);
assert_eq!(mate_pos(&rec), 12345);
}

#[test]
fn test_set_template_length_roundtrip() {
let mut rec = make_bam_bytes(0, 0, 0, b"r1", &[], 4, -1, -1, &[]);
assert_eq!(template_length(&rec), 0);
set_template_length(&mut rec, -300);
assert_eq!(template_length(&rec), -300);
}

#[test]
fn test_set_ref_id_roundtrip() {
let mut rec = make_bam_bytes(0, 0, 0, b"r1", &[], 4, -1, -1, &[]);
assert_eq!(ref_id(&rec), 0);
set_ref_id(&mut rec, 3);
assert_eq!(ref_id(&rec), 3);
}

#[test]
fn test_set_pos_roundtrip() {
let mut rec = make_bam_bytes(0, 0, 0, b"r1", &[], 4, -1, -1, &[]);
assert_eq!(pos(&rec), 0);
set_pos(&mut rec, 9999);
assert_eq!(pos(&rec), 9999);
}
}
Loading
Loading