gwseq-io 0.2.1

//! Checks against real CRAM files, which run only when those files are there.
//!
//! Everything else in this crate's suite builds what it reads, and that is
//! worth keeping: a test that only runs on one machine is a test that stops
//! running. But a synthetic CRAM exercises the paths *its writer* chose, and
//! for a format with eight block methods, thirty data series and a reference
//! outside the file, those are a small corner of it. A file `samtools` wrote
//! reaches the rest.
//!
//! So these skip themselves when `local/test_data` is empty, and
//! `tests/cram_roundtrip.rs` covers the same ground with a file it builds. The
//! two answer different questions: "does this agree with the specification" and
//! "does this agree with the tool everyone else uses".
//!
//! # The check that is not here
//!
//! The strongest one cannot live in `cargo test`, because it needs `samtools`:
//!
//! ```text
//! samtools view -T ref.fa file.cram             > a.sam
//! cargo run --release --example cramdump -- \
//!     file.cram -T ref.fa                       > b.sam
//! cmp a.sam b.sam
//! ```
//!
//! Run over a 10-million-read CRAM 3.1 — rANS N×16, the name tokeniser, gzip
//! and raw blocks, multi-reference slices, indels, soft clips, and sequences
//! rebuilt from a bgzip-compressed reference — that comparison is byte for
//! byte identical, `MD` and `NM` included. That is what `examples/cramdump.rs`
//! exists for.
//!
//! `local/probe/cram_codec_matrix.sh` runs that comparison over the same file
//! re-encoded every way `samtools` offers — both versions, all four profiles,
//! and the options no profile turns on. It is the **only** check that reaches
//! rANS 4x8, bzip2, the adaptive arithmetic coder or fqzcomp with real data:
//! the test file is a `samtools`-default CRAM 3.1 and uses none of them.
//! Between them they are what `-O cram,version=3.0` and `-O cram,archive`
//! produce, so "the test file passes" says very little about them.

use super::*;
use crate::cram::crai::{CramIndex, IndexEntry};
use crate::source::ByteSource;

fn test_file(name: &str) -> Option<String> {
    let path = format!(
        "{}/../../local/test_data/{name}",
        env!("CARGO_MANIFEST_DIR")
    );
    std::path::Path::new(&path).exists().then_some(path)
}

/// Walk every container and decode every block, which is the container, block
/// and codec layers over bytes this crate did not write.
#[test]
fn every_block_of_a_real_cram_decodes() {
    let Some(path) = test_file("AtTPax7_H3K4me.10M.cram") else {
        return;
    };
    let source = crate::source::open(&path, None, None).expect("opens");
    let head = source.read_at(0, 64).expect("head");
    let definition = container::FileDefinition::parse(&head, &path).expect("file definition");
    assert_eq!((definition.major, definition.minor), (3, 1));

    let mut offset = container::FILE_DEFINITION_SIZE as u64;
    let mut containers = 0usize;
    let mut blocks = 0usize;
    let mut methods = std::collections::BTreeMap::new();
    let mut names: Vec<String> = Vec::new();
    let file_len = source.len().expect("length");

    while offset < file_len {
        let header = container::ContainerHeader::read(source.as_ref(), offset).expect("container");
        if header.is_eof() {
            break;
        }
        let body = source
            .read_exact_at(header.blocks_offset(), header.length as usize)
            .expect("container body");
        let mut at = 0usize;
        for _ in 0..header.n_blocks {
            let Some(rest) = body.get(at..) else { break };
            if rest.is_empty() {
                break;
            }
            let block = container::Block::parse(rest, header.blocks_offset() + at as u64, &path)
                .unwrap_or_else(|e| panic!("container at {offset}, block at {at}: {e}"));
            at += block.total_size;
            blocks += 1;
            *methods.entry(block.method.name()).or_insert(0usize) += 1;

            // The name tokeniser is the one codec whose output has a shape
            // worth checking beyond its length: nul-terminated read names.
            if block.method == container::CompressionMethod::NameTok {
                let decoded: Vec<&[u8]> = block
                    .data
                    .split(|b| *b == 0)
                    .filter(|s| !s.is_empty())
                    .collect();
                assert!(
                    decoded.len() > 9000,
                    "only {} names in a slice of ten thousand",
                    decoded.len()
                );
                for name in &decoded {
                    assert!(
                        name.iter().all(|b| b.is_ascii_graphic()),
                        "a read name with something unprintable in it: {:?}",
                        String::from_utf8_lossy(name)
                    );
                }
                if names.is_empty() {
                    names = decoded
                        .iter()
                        .take(2)
                        .map(|n| String::from_utf8_lossy(n).into_owned())
                        .collect();
                }
            }
        }
        containers += 1;
        offset = header.end_offset();
        // Twenty containers is two hundred thousand alignments and every codec
        // the file uses; walking all of it belongs in a benchmark.
        if containers >= 20 {
            break;
        }
    }
    assert!(containers >= 20, "only {containers} containers walked");
    assert!(blocks > 400, "only {blocks} blocks decoded");
    assert!(!names.is_empty(), "no read names were decoded");
    println!("{containers} containers, {blocks} blocks: {methods:?}");
    println!("read names: {names:?}");
}

/// Open through the public API and read a locus: index build, container and
/// compression headers, slice decode, record reconstruction against a real
/// reference, and the BAM record decoder behind all of it.
#[test]
fn a_real_cram_reads_its_entries() {
    let Some(path) = test_file("AtTPax7_H3K4me.10M.cram") else {
        return;
    };
    let reference = test_file("mm10.fa.gz");
    let reader = CramReader::open(&path, None, reference.as_deref(), 4, None, None).expect("opens");
    assert_eq!(reader.version(), (3, 1));
    assert!(reader.is_indexed(), "{}", reader.index_error());

    let request = crate::bam::EntriesRequest::new(
        crate::genomic::Locs::spans(&["chr1".to_string()], &[3_100_000], &[3_101_000])
            .expect("locs"),
    );
    let entries = reader.read_entries(&request).expect("reads");
    assert_eq!(entries.len(), 1);
    let alignments = &entries[0];
    assert!(!alignments.is_empty(), "no alignments in the locus");

    for entry in alignments {
        assert_eq!(entry.chr(), "chr1");
        assert!(entry.end() > 3_100_000 && entry.start() < 3_101_000);
        // The sequence was rebuilt from the reference, so it must be as long as
        // the cigar says the read is — and made of bases.
        assert_eq!(entry.sequence().len() as i64, entry.query_length());
        assert!(entry.sequence().bytes().all(|b| b"ACGTN=".contains(&b)));
        assert_eq!(entry.qualities().len(), entry.sequence().len());
    }
}

/// With no reference, everything but the sequence still decodes — which is the
/// property the whole reference layer is designed around.
/// An index built from container headers answers a locus with the same slices
/// the `.crai` does.
///
/// `samtools` writes a multi-reference slice wherever a container crosses a
/// chromosome boundary, so a sorted file has one per chromosome. Filed as
/// "could be anywhere", every one of them joined every query: a locus that
/// should read one slice read eighteen, each a ten-thousand-record decode.
/// A `.crai` does not have the problem because `samtools index` lists such a
/// slice once per reference it touches, with the real span — which is what
/// `build` now does too.
#[test]
fn an_index_built_from_containers_answers_a_locus_like_the_crai() {
    let Some(path) = test_file("AtTPax7_H3K4me.10M.cram") else {
        return;
    };
    let Some(crai) = test_file("AtTPax7_H3K4me.10M.cram.crai") else {
        return;
    };
    let source = crate::source::open(&path, None, None).expect("opens");
    let built = CramIndex::build(source.as_ref()).expect("builds");
    let listed = CramIndex::parse(
        &crate::source::open(&crai, None, None)
            .expect("opens")
            .read_to_end(0)
            .expect("reads"),
        &crai,
    )
    .expect("parses");

    // Nothing is left filed as "could be anywhere".
    assert_eq!(
        built.multi_len(),
        0,
        "a built index should place every slice on a reference"
    );

    for (ref_id, start, end) in [
        (0, 3_100_000, 3_101_000),
        (1, 5_000_000, 5_001_000),
        (2, 0, 100_000),
    ] {
        let from_build = built.slices(ref_id, start, end);
        let from_crai = listed.slices(ref_id, start, end);
        assert_eq!(
            from_build.len(),
            from_crai.len(),
            "reference {ref_id}:{start}-{end}: {} built against {} listed",
            from_build.len(),
            from_crai.len()
        );
        let key = |e: &IndexEntry| (e.container_offset, e.landmark);
        assert_eq!(
            from_build.iter().map(key).collect::<Vec<_>>(),
            from_crai.iter().map(key).collect::<Vec<_>>(),
        );
    }
}

#[test]
fn a_real_cram_without_a_reference_reads_everything_but_its_sequences() {
    let Some(path) = test_file("AtTPax7_H3K4me.10M.cram") else {
        return;
    };
    // A hundred kilobases, unfiltered. The thousand-base default this used to
    // ask for held five reads, none of them with a substitution, so the
    // all-N assertion below held for the wrong reason: a substituted base was
    // being invented from the matrix's `N` row and nothing here reached one.
    let locs = || {
        crate::genomic::Locs::spans(&["chr1".to_string()], &[3_100_000], &[3_200_000])
            .expect("locs")
    };
    let request = || crate::bam::EntriesRequest::new(locs()).filter(false);
    // `/dev/null` is a readable path that is not a FASTA, so resolution fails
    // outright rather than falling back to the header's `UR`.
    let without = CramReader::open(&path, None, Some("/dev/null"), 2, None, None).expect("opens");
    assert!(
        !without.reference_error().is_empty(),
        "a reference that cannot be read should say so"
    );
    let bare = without
        .read_entries(&request())
        .expect("reads without a reference");

    let reference = test_file("mm10.fa.gz");
    let with = CramReader::open(&path, None, reference.as_deref(), 2, None, None).expect("opens");
    let full = with
        .read_entries(&request())
        .expect("reads with a reference");

    assert_eq!(bare[0].len(), full[0].len());
    for (bare, full) in bare[0].iter().zip(&full[0]) {
        assert_eq!(bare.read_name(), full.read_name());
        assert_eq!(bare.start(), full.start());
        assert_eq!(bare.end(), full.end());
        assert_eq!(bare.flag(), full.flag());
        assert_eq!(bare.cigar(), full.cigar());
        assert_eq!(bare.qualities(), full.qualities());
        assert_eq!(bare.next_start(), full.next_start());
        assert_eq!(bare.template_length(), full.template_length());
        // The sequence is the one thing that needed the reference. It is still
        // the right length; it is simply not known.
        assert_eq!(bare.sequence().len(), full.sequence().len());
        assert!(
            bare.sequence().bytes().all(|b| b == b'N'),
            "{}: {} without a reference",
            bare.read_name(),
            bare.sequence()
        );
    }
    // And the window really does contain substitutions, so the assertion above
    // is not vacuous: without one, every read is a pure match and the `X`
    // branch is never taken.
    let substituted = full[0]
        .iter()
        .filter(|e| {
            e.tags().is_ok_and(|tags| {
                tags.iter()
                    .any(|(tag, value)| tag == "NM" && format!("{value:?}") != "Int(0)")
            })
        })
        .count();
    assert!(
        substituted > 0,
        "this window holds no mismatching read, so it cannot test substitution"
    );
}