Skip to content

seqair 0.3.0

Choose a tag to compare

@killercup killercup released this 22 Sep 14:48
· 327 commits to main since this release

Breaking

Every genomic interval in the API is now one closed span, core::range::RangeInclusive<Pos0>,
instead of two loose positions (see r[interval.span_type]). The reference for a region is fetched
with the very value that queried it, so there is no + 1 at any boundary for a caller to get wrong,
and a span ending on Pos0::MAX names its last base where a half-open end could not.

  • fetch_into(tid, start, end, store) → fetch_into(tid, span, store), on Readers,
    IndexedReader and every format reader; same for fetch_into_customized,
    estimate_region_bytes, IndexedBamReader::query and PileupEngine::new. These were already
    inclusive on both ends: (start..=end).into() is the whole migration.
  • IndexedFastaReader::fetch_seq / fetch_seq_into and Readers::fetch_base_seq were half-open
    and are now closed.
    fetch_seq(name, p(0), p(4)) becomes fetch_seq(name, (p(0)..=p(3)).into()).
    fetch_seq_into_u64, the side door for reaching the last representable base, is gone — the
    closed span reaches it. FastaError::RegionOutOfBounds reports last (inclusive) instead of
    end, and fires for start > last or last >= seq_len.
  • The segment target (resolver, start, end) is (resolver, span).
  • Segment::end() → Segment::last(): it returned the last covered position under a name that
    reads as one past it. Segment::core_range() → core_span(), now a core::range::RangeInclusive
    (.start / .last; Copy, 8 bytes, not an iterator). New Segment::span().
  • Pos0::max_value() → Pos0::MAX.
  • MSRV 1.92.0 → 1.98.1. It was already effectively 1.98.1 via seqair-types; core::range needs
    1.96.
  • seqair-types is depended on with default-features = false and a forwarding serde feature was
    added. seqair serializes nothing itself, so its users no longer compile serde for nothing; enable
    seqair/serde to serialize seqair-types values.

Fixed

  • A reversed query returned the reads spanning it. fetch_into with start > end handed the
    overlap test an interval that names no positions, and the test — pos <= end && end_pos >= start
    — is satisfied by exactly the records covering the whole gap. Every reader now treats a span with
    last < start as what it is, empty, and returns nothing. Found by a property test written for the
    new span API; the behaviour predates it.

  • A query on a truncated BAM never returned. Not slowly — never. IndexedBamReader::fetch_into
    on a file truncated anywhere past its header, with its index left intact, spun forever with no
    panic, no allocation and no error, which is why 26 fuzz targets never saw it: libFuzzer could
    only ever have reported it as a timeout. RegionBuf::read_record reported two different facts as
    the same BgzfError::UnexpectedEof — "the window ran out mid-chunk, refill" and "the file ended
    here" — and the query loop read every one as the first, looping back on the assumption that the
    chunk step would happen on the next turn. At the end of the planned ranges the cursor stops
    advancing, so it never did. read_record now returns Ok(None) when the ranges are exhausted at
    a record boundary and keeps UnexpectedEof for running out partway through a record, and the
    chunk step is explicit, so every turn of the loop consumes either a record or a chunk. A
    truncated file now yields the records before the cut, as htslib does with one.

  • CRAM: a multi-reference container could decode reads as N with no error. Such a container
    holds one index entry per slice per reference, and the reference window was taken from the first
    entry whose container offset matched. When a later slice reached further along the reference, the
    window stopped short and the tail of those reads came back as N — silently, since running past
    the fetched reference is only a warning. It is the union of every matching entry now. htslib turns
    multi-reference slices on by itself once a container would hold few records per reference, so an
    ordinary file with short contigs reaches this; two slices one base apart are enough.

  • CRAM: embed_ref=2 files could not be opened at all. The slice-header MD5 was checked against
    the FASTA whether or not the slice carried its own reference. Under embed_ref=2 htslib embeds a
    consensus computed from the reads and digests that, so it matches the external reference only
    where the reads happen to agree with it — every such file with low coverage or a real difference
    failed with ReferenceMd5Mismatch. A slice with an embedded reference is now exempt.

  • CRAM: a no-reference file mis-decoded any read with an insertion. The Q and q features
    carry quality and nothing else, but both were decoded as an anchoring reference match, adding a
    base and an M operation. In no-reference mode htslib emits one Q per inserted base next to
    the I feature, so those reads came back one base too long and were refused with
    QualLenMismatch. It needs no option to reach: htslib drops into no-reference mode by itself when
    embed_ref meets multi_seq_per_slice.

  • A CSI could silently drop records from a region query. A bin's loffset was written as
    that bin's own first chunk offset. htslib derives it from the linear index instead — the first
    record at or after the start of the bin's leftmost leaf window — and the two differ whenever a
    record in another bin begins earlier inside that window, which a record straddling a window
    boundary does routinely. The bin's own chunk is then the larger of the two, and that is the
    dangerous direction: a reader takes min_off from loffset and discards every chunk ending
    before it, so tabix and bcftools view -r both returned short, with exit status 0. The record
    was never lost from the file and a query for its exact position still found it; only a query
    whose range started earlier missed it. Smallest case: a record at 0-based 16384 with a second
    straddling the next window, queried from position 1.

  • The VCF/BCF writer could not index a contig over 512 Mbp. Its index was built with
    min_shift=14, depth=5 whatever the header said, so bins ran out at 2^29 — a record above that
    was written to the file, pushed to the index, and then unreachable, with bcftools view -r
    returning nothing while the same query over bcftools index -c's CSI returned it. Nothing
    errored. The depth is now derived from the longest reference, as htslib's
    hts_adjust_csi_settings does. This is the case CSI exists for; TBI and BAI cannot express it.
    Indexes that already fitted the depth-5 scheme are byte-identical, since the search floors at 5.

  • OwnedBamRecord::to_bam_bytes accepted a record whose quality no longer matched its sequence.
    set_seq validates the new sequence against the CIGAR and deliberately leaves qual alone, so
    set_seq with no following set_qual was a reachable state where the two disagreed. It
    serialized: l_seq governs how many quality bytes a reader consumes, so the record did not
    bounce, it was misread from the sequence onwards. It now raises the same SeqQualLengthMismatch
    the builder and set_qual raise. An empty qual remains legal at any length.

  • Non-finite floats in VCF text now match htslib's spelling. write_float_g sent NaN and the
    infinities through Display, writing NaN where C's %g — and so htslib — writes nan. VCF 4.3
    §1.3 admits INF/INFINITY/NAN case-insensitively as Float values and §6.3.3 gives quiet NaN
    first-class status, distinct from the missing sentinel, so these are values rather than errors and
    are deliberately not written as .: that would turn a value into a missing value and put the
    text output at odds with the BCF output, which writes the caller's bits through unchanged.

  • The SAM reader rejected every negative element of a signed B array. B:c, B:s and
    B:i are the signed subtypes, but all three were converted through the unsigned type of the
    same width, so XX:B:c,-1 failed the whole record with InvalidAuxValue. Each subtype is now
    converted through the type it names. An unrecognised subtype is also an error rather than a
    B tag whose element count promises more bytes than follow it.

  • VCF float text now matches bcftools exactly. write_float_g claimed to write C's %g
    with six significant digits — the format htslib emits — but only ever produced the fixed form,
    so a value outside 1e-4 .. 1e6 came out as 0.0000610352 or 1234567 where bcftools writes
    6.10352e-05 and 1.23457e+06. It now switches forms on the value's exponent after rounding
    to six significant digits (so 0.0001 stays 0.0001 and does not become 1e-04), writes the
    exponent C's way with a sign and at least two digits, and writes negative zero as -0. Values
    in 1e-4 .. 1e6 are unchanged, which is every quality score and most everything else.

  • A float below about 1e-25 could not be written to VCF at all. write_float_g sized its
    decimal places from the value's magnitude and formatted into a 32-byte buffer, so six significant
    digits of 5.169879e-26 overflowed it and the write failed with FormattedFloatLongerThan32Chars,
    taking the whole record with it. Those magnitudes do not occur in QUAL, but they do in a p-value
    carried in INFO. The writer now falls back to scientific notation there, which is what %g and
    htslib emit. Values that already fit are byte-identical to before.

Added

  • IndexBuilder::csi builds a CSI whose depth covers a given longest-reference length, and
    IndexBuilder::csi_depth_for exposes that calculation on its own. IndexBuilder::BAI_DEPTH names
    the 5 that BAI and TBI are fixed at and that CSI now treats as a floor.
  • seqair::bam::DecodeError is re-exported. RecordStore::set_alignment, push_raw and
    push_fields are public and return it, but it was only reachable through a private module, so
    a caller outside the crate could not match on the error.

Changed

  • RegionBuf::read_record returns Result<Option<&[u8]>, BgzfError>. Ok(None) means the planned
    ranges are exhausted at a record boundary — there is no next record; Err(UnexpectedEof) now
    means only that a record was cut in half. Callers that treated the old single error as "retry"
    must handle Ok(None) as an end, not a hiccup; see the truncation fix above for why.

Internal

  • Property-based tests moved from proptest to hegel (hegeltest). Case
    counts live in hegel.toml at the workspace root: 256 locally, 1000 on CI, and a thorough
    profile at 10000. No public API change; proptest is gone from the dev-dependencies.
  • BgzfWriter's virtual offsets, the two CustomizeRecordStore filter hooks, and the
    Pos/QPos type walls all gained property tests. The BGZF ones were checked against the
    offset bug fixed in the previous release: reverting it makes them fail on a two-write case.
  • Six cross-implementation comparisons that ran on fixed fixtures now run on generated inputs:
    the pileup against htslib's bam_plp_auto (and now comparing per-alignment qpos, not only
    depth); the write → index → query round-trip against both samtools view and the generated
    records; CRAM against BAM decodes of the same reads; multi-sample FORMAT arrays against
    bcftools; SAM through the reader and write_store_record and back out through samtools; and
    records built with OwnedBamRecord::builder through BamWriter and back.
  • Generated coverage now reaches the corners that were fixture-only, and found the six bugs listed
    above: the CRAM writer-option matrix (version × embed_ref × seqs_per_slice ×
    slices_per_container, with multi-reference slices and span-0 index entries asserted per case
    rather than counted across the run), OwnedBamRecord under sequences of mutations refereed by
    samtools, reader-level filtering compared field by field against the unfiltered fetch, extras
    through the pileup, INFO types and missing-value integer arrays against bcftools, CSI output
    including contigs past 2^29, mate fields refereed by samtools fixmate, and which typed error a
    corrupt or truncated file actually surfaces.
  • Spec traceability: 13 dangling rule ids resolved to zero — six were renames, seven rules the code
    was annotated against but nobody had written — one stale reference fixed, and 28 rules in the
    bcf_writer.* and csi.* namespaces annotated against round-trips that already pinned them.
    Five rules were deliberately left unverified rather than annotated on tests that would not fail if
    they broke. Three rules whose prose had drifted from the code (record_store.customize.trait,
    perf.reuse_alignment_vec, unified.segment_options) were retold to match it.
  • The six bare compile_fail blocks in the VCF writer each gained a passing twin that uses every
    name the rejection uses, so a rename fails the guard rather than silently making the rejection
    free; each was mutated into its legal form to confirm it rejects for the reason it claims.
  • The nightly workflow (renamed fuzz.yml → nightly.yml) grew from one job to three, because
    there are three axes and they do not substitute for one another: fuzz as before, thorough at
    the 50000-case profile, and deep-oracles with HEGEL_TEST_CASES=500. The last is the only
    thing that reaches the 83 tests which pin their own case count — the subprocess round-trips
    against samtools, bcftools, htslib and noodles — and it found the CSI loffset bug above on its
    first outing. Both numbers are measured, not guessed, and the curves are recorded in hegel.toml
    and the workflow comments.