gwseq-io 0.2.1

//! Print a CRAM's alignments as SAM, for comparing against another reader.
//!
//! ```text
//! cargo run --release -p gwseq-io --example cramdump -- FILE.cram [-T ref.fa] [region] [limit]
//! ```
//!
//! The point of it is `diff`. CRAM has five entropy coders and thirty data
//! series, and a decoder that is subtly wrong produces *plausible* records
//! rather than an error — so the check that matters is whether this reader and
//! an independent one agree, line for line, on a real file:
//!
//! ```text
//! samtools view -T ref.fa file.cram chr1:3,000,000-3,100,000 > a.sam
//! cargo run --release -p gwseq-io --example cramdump -- file.cram -T ref.fa chr1:3000000-3100000 > b.sam
//! diff a.sam b.sam
//! ```
//!
//! Tags are printed in the order the record stores them, which is the order
//! `samtools view` prints them in too.

use gwseq_io::bam::{EntriesRequest, TagValue};
use gwseq_io::cram::CramReader;
use gwseq_io::genomic::Locs;

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let args: Vec<String> = std::env::args().skip(1).collect();
    if args.is_empty() {
        eprintln!("usage: cramdump FILE.cram [-T reference.fa] [region] [limit]");
        std::process::exit(2);
    }

    let mut path = None;
    let mut reference = None;
    let mut region = None;
    let mut limit = usize::MAX;
    let mut rest = args.iter();
    while let Some(arg) = rest.next() {
        match arg.as_str() {
            "-T" => reference = rest.next().cloned(),
            other if path.is_none() => path = Some(other.to_string()),
            other if region.is_none() => region = Some(other.to_string()),
            other => limit = other.parse().unwrap_or(usize::MAX),
        }
    }
    let path = path.expect("a file");

    let reader = CramReader::open(&path, None, reference.as_deref(), -1, None, None)?;
    if !reader.reference_error().is_empty() {
        eprintln!("no reference: {}", reader.reference_error());
    }

    // No flag filtering: this is meant to line up with `samtools view`, which
    // prints every record it is given.
    let whole_file = region.is_none();
    let entries = match &region {
        Some(region) => {
            let (chr, start, end) = parse_region(region, &reader)?;
            let request = EntriesRequest::new(Locs::spans(&[chr], &[start], &[end])?).filter(false);
            reader
                .read_entries(&request)?
                .into_iter()
                .flatten()
                .collect()
        }
        None => Vec::new(),
    };

    let out = std::io::stdout();
    let mut out = std::io::BufWriter::new(out.lock());
    use std::io::Write as _;

    // Streamed, not collected. `read_all_entries` holds every record of the
    // file at once — ten million `BamRecord`s, each pinning the slice buffer
    // it borrows from — which peaked at 8.3 GB against `samtools view`'s
    // 334 MB. Since this example is also the timing comparison the README
    // quotes, measuring the wrong shape would be worse than not measuring.
    let mut streamed: Option<_> = None;
    if whole_file {
        let request = EntriesRequest::new(Locs::spans(&[], &[], &[])?).filter(false);
        streamed = Some(reader.iter_all_entries(&request, 1 << 22)?);
    }

    let mut written = 0usize;
    let emit = |entry: &gwseq_io::bam::BamRecord,
                out: &mut std::io::BufWriter<std::io::StdoutLock<'_>>|
     -> Result<(), Box<dyn std::error::Error>> {
        write!(
            out,
            "{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}",
            entry.read_name(),
            entry.flag(),
            entry.chr(),
            entry.start() + 1,
            entry.mapping_quality(),
            if entry.cigar().is_empty() {
                "*"
            } else {
                entry.cigar()
            },
            mate_name(entry.chr(), entry.next_chr()),
            entry.next_start() + 1,
            entry.template_length(),
            entry.sequence(),
            entry.qualities(),
        )?;
        for (tag, value) in entry.tags()? {
            write!(out, "\t{tag}:{}", format_tag(value))?;
        }
        writeln!(out)?;
        Ok(())
    };

    match streamed {
        Some(windows) => {
            for batch in windows {
                let batch = batch?;
                for entry in batch.iter() {
                    if written >= limit {
                        return Ok(());
                    }
                    emit(entry, &mut out)?;
                    written += 1;
                }
            }
        }
        None => {
            for entry in entries.iter().take(limit) {
                emit(entry, &mut out)?;
            }
        }
    }
    Ok(())
}

/// SAM writes `=` for a mate on the same reference as its read.
fn mate_name<'a>(chr: &str, next: &'a str) -> &'a str {
    if next == chr {
        "="
    } else {
        next
    }
}

fn format_tag(value: &TagValue) -> String {
    match value {
        TagValue::Char(c) => format!("A:{c}"),
        TagValue::Int(v) => format!("i:{v}"),
        TagValue::Float(v) => format!("f:{v}"),
        TagValue::Str(s) => format!("Z:{s}"),
        TagValue::IntArray(values) => format!(
            "B:i{}",
            values.iter().map(|v| format!(",{v}")).collect::<String>()
        ),
        TagValue::FloatArray(values) => format!(
            "B:f{}",
            values.iter().map(|v| format!(",{v}")).collect::<String>()
        ),
    }
}

/// `chr1`, `chr1:1000`, `chr1:1000-2000`, with commas allowed as samtools
/// allows them.
fn parse_region(
    region: &str,
    reader: &CramReader,
) -> Result<(String, i64, i64), Box<dyn std::error::Error>> {
    let (chr, span) = match region.split_once(':') {
        Some((chr, span)) => (chr, Some(span)),
        None => (region, None),
    };
    let size = reader.chr_sizes().resolve(chr)?.size;
    let Some(span) = span else {
        return Ok((chr.to_string(), 0, size));
    };
    let span = span.replace(',', "");
    let (start, end) = match span.split_once('-') {
        Some((start, end)) => (start.parse::<i64>()?, end.parse::<i64>()?),
        None => (span.parse::<i64>()?, size),
    };
    // Regions are 1-based and inclusive, as samtools writes them.
    Ok((chr.to_string(), (start - 1).max(0), end.min(size)))
}