Repository navigation
feat: add header-only handlers for VCF/gVCF, BCF, BAM, CRAM, SAM, FASTQ and FASTA - #127
Conversation
|
@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.
620c0d5 to
9e151b4
Compare
|
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
left a comment
There was a problem hiding this comment.
@renato-umeton @rafiattrach Looks good and independently reproduced. Minor changes we may want to do before merge:
- 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=VCFand 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 withread_prefix_chunksand a per-line limit. VCF should do the same. BCF is not affected, because it caps l_text at 64 MiB. - The docs promise more privacy than the flag gives.
--genomic-sample-idswithholds sample column names and@RGSM tags. File names still appear in every FileObject'scontentUrl, 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. ##reference=goes into the VCF description verbatim. Callers often write it as a local path, for examplefile:///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.- 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.
|
@slobentanzer thx for the careful pass. 1 to 3 fixed here, 4 as a follow up:
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
left a comment
There was a problem hiding this comment.
Looks good; disclaimer @rafiattrach I didn't independently run again, but don't want to hold up the process. :)
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.
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.pyclaims on the##fileformat=VCFmagic (so vCard.vcffiles are refused), reads##and#CHROMheader 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),INFOwith typed sub-fields from##INFO, and, when samples are present,FORMATsub-fields plus one repeatedsamplesfield. Integer tocr:Int64, Float tocr:Float64, Flag tosc:Boolean, String and Character tosc:Text;Numberother than 0 or 1 marks the sub-field repeated. gVCF is detected from##GVCFBlockor##ALT=<ID=NON_REF.handlers/bcf_handler.pysubclasses the VCF handler. BCF is BGZF-wrapped and the compression layer does not strip.bcf, so the handler gunzips in place, reads theBCF\2magic,l_text(capped) and the VCF header text, and emits the same RecordSet. A.bcf.gzarrives 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.pygunzips in place, reads the magic,l_text(capped at 64 MiB) and the header, and stops.handlers/cram_handler.pyreads 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.pyreads@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.pyreads 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.pyreads 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 SMtags in BAM, CRAM and SAM form a sample manifest. By default the output carries counts only;--genomic-sample-idsopts 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.peekalso catchesEOFError,lzma.LZMAErrorandzlib.error(named once assources.UNREADABLE), so a corrupt or truncated compressed input falls throughclaimsto a refusal with a reason instead of escaping dispatch; every handler catches the same set inextractand names the file.read_exactly,read_prefix_chunks,decompress_prefix,plural,MAX_HEADER_BYTES) live inhandlers/utils.py._file_objects_forhonours an optionaldescriptionkey fromextract; a bake of the MIMIC-IV demo is byte-identical before and after..bai,.crai,.csi,.tbi,.fai,.gzi) are reported as unclaimed; nothing opens them.Testing
tracemalloc, full bakes validated under mlcroissant, refusals driven through a bake and read back from the scan report), plusSAMPLESentries 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.pychecks 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 withmainafter the HDF5 handler landed.