Skip to content

feat: add header-only handlers for VCF/gVCF, BCF, BAM, CRAM, SAM, FASTQ and FASTA - #127

Merged
rafiattrach merged 46 commits into
MIT-LCP:mainfrom
renato-umeton:feat/genomic-handlers
Sep 30, 2026
Merged

rafiattrach merged 46 commits into
MIT-LCP:mainfrom
renato-umeton:feat/genomic-handlers

Conversation

@renato-umeton

@renato-umeton renato-umeton commented Sep 9, 2026 •

Copy link
Copy Markdown
Collaborator

Closes #125.

Adds header-only handlers for the genomic formats a sequencing delivery is made of, with no new dependencies: VCF/gVCF, BCF, BAM, CRAM, SAM, FASTQ and FASTA. Together with the existing TSV handler this covers a complete St. Jude Cloud Genomics Platform delivery (BAM or CRAM, gVCF, somatic VCF, CNV, feature counts) and the reads and reference that sit next to it. Every handler reads a header or a first record and stops; no alignment, variant or sequence record is ever read, and the read is bounded by the header, not the file.

Variant callsets: VCF, gVCF and BCF

  • handlers/vcf_handler.py claims on the ##fileformat=VCF magic (so vCard .vcf files are refused), reads ## and #CHROM header lines only, and refuses a header that does not declare the eight mandatory columns. One RecordSet per file: CHROM, POS, ID, REF, ALT (repeated), QUAL, FILTER (repeated), INFO with typed sub-fields from ##INFO, and, when samples are present, FORMAT sub-fields plus one repeated samples field. Integer to cr:Int64, Float to cr:Float64, Flag to sc:Boolean, String and Character to sc:Text; Number other than 0 or 1 marks the sub-field repeated. gVCF is detected from ##GVCFBlock or ##ALT=<ID=NON_REF.
  • handlers/bcf_handler.py subclasses the VCF handler. BCF is BGZF-wrapped and the compression layer does not strip .bcf, so the handler gunzips in place, reads the BCF\2 magic, l_text (capped) and the VCF header text, and emits the same RecordSet. A .bcf.gz arrives with one layer already removed and is read the same way. BCF1 is claimed and refused with a reason.

Alignments: BAM, CRAM and SAM

All three carry the same SAM text header, parsed once in handlers/sam_header.py. Each emits a described FileObject (sort order, reference count and assembly, read groups with platform and centre, program chain) and no RecordSet, since aligned reads are not records of a dataset schema.

  • handlers/bam_handler.py gunzips in place, reads the magic, l_text (capped at 64 MiB) and the header, and stops.
  • handlers/cram_handler.py reads the 26-byte file definition, walks the first container header (ITF8/LTF8 integers, record counter width per major version, CRC on version 3) and the FILE_HEADER block. Versions 2.x and 3.x are read; 1.x and 4.x are refused with a reason. Raw, gzip (including bare zlib, as htslib reads it), bzip2 and lzma header blocks are decoded incrementally against the size the block declares, so a block that lies about its size is refused before it is expanded; rANS-coded header blocks are refused with a reason. The reference the CRAM was encoded against is not needed, because only the header is read.
  • handlers/sam_handler.py reads @ lines in bounded chunks and stops at the first alignment line. Claimed on extension and header line together, since FASTQ also opens with @.

Reads and references: FASTQ and FASTA

  • handlers/fastq_handler.py reads the first record only and reports its read length. Read names are never emitted: Illumina names carry instrument, run and flowcell identifiers.
  • handlers/fasta_handler.py reads the first line only and reports the format. Record names are never emitted: a per-sample assembly names its sample there.

Governance default

VCF and BCF sample columns and @RG SM tags in BAM, CRAM and SAM form a sample manifest. By default the output carries counts only; --genomic-sample-ids opts in to listing identifiers. Metadata-only exports of controlled-access data therefore do not leak the manifest. FASTQ read names and FASTA record names have no opt-in and are never emitted.

Supporting changes

  • FileSource.peek also catches EOFError, lzma.LZMAError and zlib.error (named once as sources.UNREADABLE), so a corrupt or truncated compressed input falls through claims to a refusal with a reason instead of escaping dispatch; every handler catches the same set in extract and names the file.
  • Shared helpers (read_exactly, read_prefix_chunks, decompress_prefix, plural, MAX_HEADER_BYTES) live in handlers/utils.py.
  • _file_objects_for honours an optional description key from extract; a bake of the MIMIC-IV demo is byte-identical before and after.
  • Index files (.bai, .crai, .csi, .tbi, .fai, .gzi) are reported as unclaimed; nothing opens them.
  • Formats table, CLI reference and the supported-formats guide regenerated and extended with one section per format.
  • The README's supported-formats table lists all seven genomic formats, and its key-features list names the sample-identifier opt-in.

Testing

  • One test module per handler (magic and extension claims, refusals with reasons, type mapping, sample opt-in, wrapped inputs, bounded reads asserted with a counting stream, memory-bounded CRAM decode asserted with tracemalloc, full bakes validated under mlcroissant, refusals driven through a bake and read back from the scan report), plus SAMPLES entries that drive the contract sweep and the compression matrix.
  • tests/data/input/genomics_htslib/ holds tiny synthetic files written by pysam 0.24.1 (htslib 1.24): SAM, BAM, CRAM 2.1, 3.0 and 3.1, VCF, BCF, FASTQ, FASTA, their index files and a CRAM 1.0 file definition. tests/test_genomics_real_files.py checks the handlers against them as an independent reading of the specifications: all five alignment containers describe the same header, VCF and BCF yield structurally identical RecordSets, index files are unclaimed, CRAM 1.0 is refused with its reason, and no sample identifier appears in the document without the opt-in.
  • uv run pytest -v: 1323 passed. uv run pre-commit run --all-files: all hooks pass. docs/generate.py: no diff. Merged with main after the HDF5 handler landed.

@renato-umeton renato-umeton changed the title feat: add VCF/gVCF and BAM header handlers feat: add header-only handlers for VCF/gVCF, BCF, BAM, CRAM, SAM, FASTQ and FASTA Sep 10, 2026
renato-umeton added a commit to renato-umeton/croissant-maker that referenced this pull request Sep 11, 2026
@rafiattrach

Copy link
Copy Markdown
Collaborator

@slobentanzer your review here would be greatly appreciated 🙏

Genomics is the one gap in the handler set that a St. Jude or CCDI delivery
falls into on its first file, and a VCF header is the cheapest structure in the
tree to recover: the producer has already written the schema down. The header
declares its reference, its contigs, and every INFO and FORMAT key with a type
and a cardinality, so a record set built from it is traceable to header bytes
alone. The read stops at the #CHROM line and no record is ever parsed, which is
the same header-only commitment DICOM, NIfTI and WFDB make.

The claim is the ##fileformat=VCF declaration rather than the extension,
because .vcf is also the vCard extension. That makes this the first handler
that reads bytes before looking at a name, so it meets every file in the
dataset: a stream it cannot decompress is one it does not claim, rather than an
exception out of dispatch for the handlers behind it.

INFO and FORMAT are per-record key-value bags, not columns, so their
declarations become sub-fields of the column that carries them. Number other
than 0 or 1 marks a sub-field repeated, which covers A, R, G and the unbounded
form in one rule.

Sample column names are a cohort manifest, so the count is emitted and the
names are not. --genomic-sample-ids opts in, taking the same route through the
generator that --count-csv-rows takes.
Every other handler has a section stating what it reads and what it refuses;
the VCF one has to state two things a reader will otherwise assume wrong: that
the claim is the magic and not the extension, and that sample identifiers are
withheld unless asked for.
BAM is the primary delivered format of a genomics platform, and its SAM header
is the cheapest provenance in the tree to recover: which assembly the reads
were placed against, which platform and centre produced them, and which
programs touched them, in order. The read stops after the reference count and
no alignment record is ever touched.

A .bam is a compressed container, not a compressed file, so the input layer is
left alone: registering .bam as a wrapper would make every BAM look like some
other format that happened to arrive gzipped. The handler decompresses what it
is given instead, and accepts either spelling the pipeline can hand over, since
a file wrapped a second time arrives with one layer already off.

No RecordSet: aligned reads are records of a genome, not of a dataset schema,
and naming columns no consumer can read through Croissant would be a promise
nobody can keep. The properties go on the FileObject instead, through an
optional description key the generator now honours. That is the smallest hook
that lets a handler describe a file it builds no nodes for, and the contract
sweep now counts it as a third way of describing one.

@rg SM tags are a cohort manifest in the same way the VCF sample columns are,
so they take the same --genomic-sample-ids opt-in and are absent, from the
metadata and from the description alike, without it.
Two things a reader will otherwise assume wrong: that a .bam is a format rather
than a compressed file, and that a handler always produces a record set. Both
are stated where the other handlers state theirs.
House style for this work: the punctuation, not the prose. No behaviour change.
Two ways a malformed file got more than it should out of these handlers.

l_text is a signed 32-bit integer the BAM chooses, and it was trusted: a file
declaring 2147483647 pulled its whole body through the reader, 4 MiB of it in
the test that now guards this, which is exactly the header-only commitment the
handler exists to keep. The length is refused against a 64 MiB cap before a
byte of it is read, so the same file now costs ten bytes.

The VCF record set was built from the fixed columns the format defines rather
than from the ones the file declared, so a #CHROM line spelled with spaces, or
truncated to three columns, still produced a full eight-field record set
describing a schema that header never stated. Extraction now refuses when the
column line does not open with the eight mandatory names, which is what makes
the guide's sentence about declared columns true.

The contract sweep counted a FileObject description as describing a batch if
any one file carried it; a handler that described the first file and dropped
the rest would have passed. It has to hold for every file of the batch.
FileSource.peek promises b"" for a file that cannot be read, and caught only
OSError. A wrapper that will not open does not raise one in general: lzma
raises LZMAError and a corrupt deflate body raises zlib.error, neither of which
derives from it. Every caller of peek is a handler sniffing magic while the
registry is still deciding who owns the file, so one of those escaping ends
dispatch for every handler behind it, and the file gets a traceback instead of
a reason. A corrupt .dcm.xz did exactly that before the VCF handler existed.

Fixed where the bytes are read rather than in each handler, which is what lets
the two handlers added here drop the blanket except they were carrying. The
prefix the BAM handler decompresses itself is still its own to guard, now by
the types a refused member actually raises.
The SAM text header is the same whether it sits in a BAM, a CRAM or a SAM
file, and the VCF header is the same in a BCF. Move the SAM parser and the
alignment description into handlers/sam_header.py, and let the VCF header be
read from any iterable of lines, so the next containers reuse both instead
of copying them.
SAM is the text spelling of what BAM holds in binary, and it carries the
same header. The handler reads the @ lines and stops at the first line
that is not one, so the cost of describing a file is the size of its
header rather than the size of the file; a BAM has a declared length in
front of the header, and a SAM has nothing but that stop.

The claim needs both the .sam extension and a real header line: a FASTQ
opens with @ as well, and a .sam carrying only alignment records states
no sort order, assembly or read group, so it is reported as unclaimed
rather than described as something it does not say it is.

No RecordSet, for the reason BAM emits none. @rg SM is withheld under
the existing --genomic-sample-ids opt-in.
What the header gives, why the claim needs both the extension and a
header line, and why no record set follows. The formats table is
regenerated from the handler's class attributes.
FASTA files carry no structure a dataset schema can hold: a description
line names a record and everything under it is bases. The handler reads
that first line to recognise the format, discards it, and emits a
described FileObject with no record set, the way BAM does.

The claim needs the extension and the leading '>' together. Either alone
is wrong: '>' is one character that ordinary text also opens with, and
the extension would claim any file a user named .fa.

Record names are withheld. A per-sample assembly names its sample on the
description line, so the name is treated as the sample identifiers
elsewhere in the tree are, and there is no opt-in for it.
What is read (one description line), why the claim needs the extension
and the leading byte together, and what is deliberately left out: record
names, comment text, record counts and sequence lengths. Regenerates the
formats table from the new handler's class attributes.
BCF is a VCF whose records are packed into a binary encoding behind a
magic and a length; the header it declares is the same text. So the
handler subclasses the VCF one and overrides only the seam that finds
that text: it decompresses the container itself, as the BAM handler
does, reads the magic, the declared length and the header bytes, and
hands them to the VCF header reader. Nothing below the header is
decoded, and one callset describes the same way in either container.

Both minor versions of BCF 2 are read, since they differ in the record
encoding and not in the header. BCF1 carries no VCF header text at all
and is refused with that reason rather than half-described, as is a
declared header length above the 64 MiB cap, which is checked before a
byte of it is pulled.
The section says what the container adds over a VCF and what it does
not: the magic and length the handler reads, the two claim spellings,
why BCF1 is reported instead of described, and that everything after
the header text is the VCF handler's, sample-id opt-in included. The
formats table is regenerated from the registry.
A CRAM carries the same SAM text header a BAM does, in the first block of
its first container. Reading it means walking the container header field by
field, because every field is written in one of CRAM's two variable-width
integer encodings and the block behind them cannot be seeked to.

The header block is decoded from raw, gzip, bzip2 or LZMA. Major versions
other than 2 and 3, a block coded with rANS, a first block that is not the
file header, and a declared size larger than any real header are each
refused with a reason rather than guessed at. CRC32 values are read past
rather than verified: the header text is what is described, and a mismatch
is a decoder's corruption report.

No record set, for the reason BAM emits none. The read-group sample tags
are withheld under the same --genomic-sample-ids opt-in.

read_exactly moves to sam_header, where both containers reach it.
What is read, what is refused and why, and that the reference the file was
encoded against is never needed because only the header is read. The
formats table is regenerated from the handler's class attributes.
read_exactly, plural and the compressed-prefix peek were spread across
sam_header, bam_handler and bcf_handler, with FASTQ reaching into the SAM
header module for a string helper that has nothing to do with SAM. They are
cross-handler helpers, so they live where the other cross-handler helpers do,
leaving sam_header.py to parse SAM headers.

The 64 MiB header cap is one constant now as well: BAM, BCF and CRAM each
declared their own copy of it for the same reason, that a container stating its
own header length may not be believed about it.

Behaviour and every refusal message are unchanged.
The compressed size and the declared raw size were both checked against the
64 MiB cap, but the decode itself was not bounded and the declared raw size was
then discarded. A 39 KB CRAM whose file header block is an LZMA stream of
256 MiB of NULs declaring a raw size of 100 was expanded in full, and then
described as a header of no reference sequences and no read groups, because the
text length read out of the NULs was zero.

The block is now decoded incrementally to one byte past what it declares, which
is enough to see that it holds more without expanding it to find out how much,
and a block holding anything other than its declared size is refused. A stream
that ends before its end-of-member marker is refused too, rather than being
described from the part of it that decoded.

Method 1 is inflated with a window argument that accepts either a gzip or a
bare zlib header, which is what htslib does; a writer emitting the zlib
spelling was previously refused as "not a gzipped file".
The supported-formats table stopped at HDF5 while the branch adds seven
handlers; a reader of the README would not know a BAM or a VCF is
described. Add one row per format in the wording the guide uses, and name
the sample-identifier opt-in among the key features, since it changes what
a metadata-only export of controlled-access data carries.
The sample-identifier bullet said read-group SM tags were counted; the
alignment handlers count read groups, not samples, and withhold the tags
outright. The BCF row put the header inside the BGZF layer, which carries
none. The image and WFDB rows are brought in line with the handlers'
declared extensions.
The previous commit carried a lock rewrite from the local uv, which
downgrades the lock revision; the lock is upstream's to change.
This branch adds no dependency, so its lock should match main's byte for
byte. An earlier commit here pinned the croissant-baker version line back
to 0.4.0, which was upstream's value when the branch started; main records
0.5.0 today, and the rebase carried the older line forward.
@renato-umeton

Copy link
Copy Markdown
Collaborator Author

Rebased on main at 0.6.0. Conflicts were in the utils imports, one docstring, the test helpers, the registry and the two docs tables, resolved by keeping both sides. No handler or test module changed beyond that. Suite green at 1455.

@slobentanzer slobentanzer left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@renato-umeton @rafiattrach Looks good and independently reproduced. Minor changes we may want to do before merge:

  1. The VCF header read has no size limit (vcf_handler.py, _read_header). It loops for raw in stream, so a file that starts with ##fileformat=VCF and has no newline gets read into memory in one go. Tested with a 300 MB file: memory use peaked at 1.4 GB before the handler refused the file. A truncated or corrupt multi-GB VCF could run the process out of memory. The SAM, FASTQ and FASTA handlers already guard against this with read_prefix_chunks and a per-line limit. VCF should do the same. BCF is not affected, because it caps l_text at 64 MiB.
  2. The docs promise more privacy than the flag gives. --genomic-sample-ids withholds sample column names and @RG SM tags. File names still appear in every FileObject's contentUrl, and sequencing deliveries often name files by sample (SJ001234_D1.bam). The "Sample identifiers" section and the README bullet could say that file names are emitted either way.
  3. ##reference= goes into the VCF description verbatim. Callers often write it as a local path, for example file:///gpfs/.../GRCh38.fa, which exposes the producer's filesystem layout. Minor, but I would either keep only the file name part, or leave it as is and mention the behaviour in the docs.
  4. Modularity: there is some duplicated code in BAM, CRAM, SAM, BCF, but could also be cleaned up later.

Otherwise looks good and ready to merge!

A VCF header states no length, and the handler read it a line at a time,
so a file opening with ##fileformat=VCF and holding no line ending was
read whole as one line: 300 MB peaked at 1.4 GB before the refusal.

Read it a chunk at a time through read_prefix_chunks, as SAM, FASTQ and
FASTA do, and refuse a line above 32 MiB or a header above 64 MiB with
the file named. The line cap sits above the SAM one because the #CHROM
line of a biobank cohort names every sample. The same 300 MB file now
peaks at 169 MB. BCF is unchanged: it states its header length and is
already capped at 64 MiB.
##reference went into the record set description verbatim, and callers
often write it as a local path such as file:///gpfs/.../GRCh38.fa, which
publishes the producer's filesystem layout.

A path, a file:// URI or a bucket URI now keeps its file name alone. An
http, https or ftp address is a public location and is kept, less any
login and any query or fragment, where a signed link carries its
credential. A build name like GRCh38 is kept as declared. BCF reads the
same header code, so it gets the same treatment.
--genomic-sample-ids gates the sample columns and @rg SM tags, but the
file name is in every FileObject name and contentUrl and in the record
set ids built from it, and deliveries often name files by sample. Say
so in the Sample identifiers section and the README bullet.
reference_name read parts.hostname and parts.port, and urlsplit and
.port raise ValueError on an unclosed IPv6 bracket, a port out of range
or a port that is no number. A VCF or BCF the previous revision
described was then refused. Keep the host text as written from netloc,
leave .port alone, and read a URL urlsplit refuses as a path. This also
keeps the brackets on an IPv6 host.
One enclosing <...> is now taken off, and a structured <ID=...,URL=...>
keeps its id while its URL or Path value is read as a bare reference.
Only a value shaped like a path is cut to its file name, so GRCh38/hg38
stays whole, and a path ending in a directory, such as gs://bkt/, is
left out, since every component of it is the producer's layout.
…ording

An http address is an assumption about where the reference is
published, and internal hosts exist, so say so. State the #CHROM line
size as several to tens of MiB, say the refused read is the cap plus at
most one chunk, and say name holds the base name and contentUrl the
path, with directory names showing too. Rewrap the long line.
A structured <ID=...,URL=...> reference kept every key besides URL and
Path verbatim, so a Description, Source or File key carrying a path
leaked the layout the location was cut to hide. State the id and the
cleaned URL or Path, and drop every other key.
A relative path with no . in its last part, such as jdoe/refs/hg38, and
an scp target such as host:/gpfs/refs/hg38 were kept whole. A value
with :/ or three or more slash-separated parts now counts as a path,
while GRCh38/hg38 stays whole. An empty reference, <> or <ID=,URL=>
now states nothing, as the docstring said.
The structured form read as key=value pairs in the description. State
it as the id followed by the location in parentheses, GRCh38
(https://h.org/GRCh38.fa), or as whichever of the two it has, and drop
anything trailing the closing bracket.
… url

An id written as a path carries the same layout the location is cut to hide, and a published URL tells a reader more than a bare file name.
@renato-umeton

Copy link
Copy Markdown
Collaborator Author

@slobentanzer thx for the careful pass. 1 to 3 fixed here, 4 as a follow up:

  1. vcf header now read in chunks via read_prefix_chunks like sam/fastq/fasta, refused at a 32 MiB line or 64 MiB header w/ the file named. line cap is above sam's 1 MiB bc the #CHROM line of a biobank cohort names every sample. your 300 MB no-newline case now costs ~48 MB over the import baseline vs ~1.3 GB before. bcf untouched 756017d (+ docstring nit dcecdcc, crlf + chunk split tests 08ae020)
  2. sample ids section + README bullet now say file and dir names show either way (name, contentUrl, record set ids) and suggest rename or --exclude 63d856a
  3. ##reference keeps only the file name when it looks like a path, file://, scp or bucket uri (file:///gpfs/.../GRCh38.fa -> GRCh38.fa, jdoe/refs/hg38 -> hg38), a bare dir or empty value is left out, build names like GRCh38 or GRCh38/hg38 as is. http/https/ftp kept as likely public, minus login + query/fragment (signed links). <ID=..,URL=..> keeps only the id + cleaned URL (else Path), id cut the same way, reads as GRCh38 (https://.../GRCh38.fa), any other key (Description, Source..) dropped since it can hold a path. bcf gets the same 61b7275 4d57c45 7b58a87 2be2197 8825769 9d9dc33
    also fixed a regression from that: a malformed url (bad port, unclosed ipv6 bracket) made the file get refused, now its still described 47232d1
    docs c3d59d9 dd675de
  4. opened refactor: dedupe the shared bits of the bam/cram/sam/bcf (and vcf) handlers #160 for 4

still merges clean w/ main and the lock passes uv lock --check on the merged tree. thx

Resolve tests/helpers.py by keeping both import sets and exporting bake_validated.

@slobentanzer slobentanzer left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good; disclaimer @rafiattrach I didn't independently run again, but don't want to hold up the process. :)

@rafiattrach
rafiattrach merged commit a57f2db into MIT-LCP:main Sep 30, 2026
3 checks passed
renato-umeton added a commit to renato-umeton/croissant-maker that referenced this pull request Oct 1, 2026
Brings in MIT-LCP#127 (VCF, BCF, BAM, CRAM, SAM, FASTQ and FASTA handlers),
MIT-LCP#154 (stable summary, file set and record set order) and the 0.7.0
release.

Conflicts:
- README.md: main's WFDB and Images rows, then this branch's
  whole-slide row and its DICOM row with the whole-slide fields.
- handlers/utils.py: union of the typing imports (Iterable from this
  branch, BinaryIO and Iterator from main).
- tests/helpers.py: union of the imports (lzma, struct, zlib, Fraction)
  and of the __all__ names (the whole-slide samples and builders, plus
  VCF_HEADER_TEXT, bake_validated, bcf_payload, cram_payload, cut_gzip).

docs/_generated is regenerated with docs/generate.py and matches the
merged file. The whole-slide golden is unchanged: the golden test
passes in both discovery orders.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

genomics support

3 participants