rustalign_io/

sam.rs

//! SAM format output

use rustalign_common::Nuc;
use std::io::Write;

/// SAM output configuration
#[derive(Debug, Clone)]
pub struct SamConfig {
    /// Print AS:i: alignment score
    pub print_as: bool,

    /// Print XS:i: secondary alignment score
    pub print_xs: bool,

    /// Print YS:i: mate alignment score
    pub print_ys: bool,

    /// Print X0:i: number of best alignments
    pub print_x0: bool,

    /// Print X1:i: number of suboptimal alignments
    pub print_x1: bool,

    /// Print XN:i: number of Ns in reference
    pub print_xn: bool,

    /// Print XM:i: number of mismatches
    pub print_xm: bool,

    /// Print XO:i: number of gap opens
    pub print_xo: bool,

    /// Print XG:i: number of gap extensions
    pub print_xg: bool,

    /// Print NM:i: edit distance
    pub print_nm: bool,

    /// Print MD:Z: mismatch string
    pub print_md: bool,

    /// Print YF:Z: mate alignment flags
    pub print_yf: bool,

    /// Print YD:i: mate distance
    pub print_yd: bool,

    /// Print YT:Z: mate type
    pub print_yt: bool,

    /// Print ZS:Z: strand
    pub print_zs: bool,

    /// Print NH:i: hit count
    pub print_nh: bool,

    /// Print HI:i: hit index
    pub print_hi: bool,

    /// Print nM:i: genome mismatches
    pub print_nm_genome: bool,

    /// Print MD:n: mismatch as number
    pub print_md_number: bool,

    /// Print jM:i: junction mismatch
    pub print_jm: bool,

    /// Print jI:i: junction insertions
    pub print_ji: bool,

    /// Print full ref coords
    pub full_ref: bool,
}

impl Default for SamConfig {
    fn default() -> Self {
        Self {
            print_as: true,
            print_xs: true,
            print_ys: false,
            print_x0: true,
            print_x1: false,
            print_xn: false,
            print_xm: true,
            print_xo: true,
            print_xg: false,
            print_nm: true,
            print_md: true,
            print_yf: false,
            print_yd: false,
            print_yt: true,
            print_zs: false,
            print_nh: true,
            print_hi: false,
            print_nm_genome: false,
            print_md_number: false,
            print_jm: false,
            print_ji: false,
            full_ref: false,
        }
    }
}

/// A SAM record
#[derive(Debug, Clone)]
pub struct SamRecord {
    /// Query name
    pub qname: String,

    /// Bitwise flag
    pub flag: u16,

    /// Reference name
    pub rname: String,

    /// Position (1-based)
    pub pos: i32,

    /// Mapping quality
    pub mapq: u8,

    /// CIGAR string
    pub cigar: String,

    /// Reference name of mate
    pub rnext: String,

    /// Position of mate
    pub pnext: i32,

    /// Template length
    pub tlen: i32,

    /// Sequence
    pub seq: Vec<Nuc>,

    /// Quality scores
    pub qual: Vec<u8>,

    /// Optional fields
    pub opt: Vec<SamOpt>,
}

/// Optional SAM field
#[derive(Debug, Clone)]
pub enum SamOpt {
    /// Integer field
    Int(String, i64),
    /// String field
    String(String, String),
    /// Float field
    Float(String, f64),
}

impl SamRecord {
    /// Create a new SAM record
    pub fn new(qname: String) -> Self {
        Self {
            qname,
            flag: 0,
            rname: "*".to_string(),
            pos: 0,
            mapq: 255,
            cigar: "*".to_string(),
            rnext: "*".to_string(),
            pnext: 0,
            tlen: 0,
            seq: vec![Nuc::N],
            qual: vec![255],
            opt: Vec::new(),
        }
    }

    /// Add an optional field
    pub fn add_opt(&mut self, opt: SamOpt) {
        self.opt.push(opt);
    }
}

impl std::fmt::Display for SamRecord {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        // Required fields
        write!(
            f,
            "{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t",
            self.qname,
            self.flag,
            self.rname,
            self.pos,
            self.mapq,
            self.cigar,
            self.rnext,
            self.pnext,
            self.tlen,
            self.seq
                .iter()
                .map(|n| n.to_ascii() as char)
                .collect::<String>()
        )?;

        // Quality
        let qual_str: String = self.qual.iter().map(|&q| (q + 33) as char).collect();
        write!(f, "{}", qual_str)?;

        // Optional fields
        for opt in &self.opt {
            match opt {
                SamOpt::Int(tag, val) => write!(f, "\t{}:i:{}", tag, val)?,
                SamOpt::String(tag, val) => write!(f, "\t{}:Z:{}", tag, val)?,
                SamOpt::Float(tag, val) => write!(f, "\t{}:f:{}", tag, val)?,
            }
        }

        Ok(())
    }
}

/// SAM file writer
pub struct SamWriter<W: Write> {
    writer: W,
    #[allow(dead_code)]
    config: SamConfig,
}

impl<W: Write> SamWriter<W> {
    /// Create a new SAM writer
    pub fn new(writer: W, config: SamConfig) -> Self {
        Self { writer, config }
    }

    /// Write the SAM header
    ///
    /// # Arguments
    /// * `ref_names` - Reference sequence names
    /// * `ref_lengths` - Reference sequence lengths (must be same length as ref_names)
    pub fn write_header(
        &mut self,
        ref_names: &[String],
        ref_lengths: &[u32],
    ) -> std::io::Result<()> {
        writeln!(self.writer, "@HD\tVN:1.0\tSO:unsorted")?;
        for (name, len) in ref_names.iter().zip(ref_lengths.iter()) {
            writeln!(self.writer, "@SQ\tSN:{}\tLN:{}", name, len)?;
        }
        writeln!(self.writer, "@PG\tID:rustalign\tPN:rustalign\tVN:0.1.0")?;
        Ok(())
    }

    /// Write a SAM record
    pub fn write(&mut self, record: &SamRecord) -> std::io::Result<()> {
        writeln!(self.writer, "{}", record)?;
        Ok(())
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_sam_config_default() {
        let config = SamConfig::default();
        assert!(config.print_as);
        assert!(config.print_md);
        assert!(!config.print_ys);
    }

    #[test]
    fn test_sam_record_new() {
        let record = SamRecord::new("read1".to_string());
        assert_eq!(record.qname, "read1");
        assert_eq!(record.flag, 0);
    }
}