Carry aligned reads as CRAM, with --aligned_format cram|bam - #207
Draft
ljwharbers wants to merge 13 commits into
Draft
ljwharbers wants to merge 13 commits into
ljwharbers wants to merge 13 commits into
Conversation
MINIMAP2_ALIGN (patched) writes a reference-compressed CRAM + .crai when bam_format is 'cram'; SAMTOOLS_MERGE gets the reference and the joins follow the cram/crai outputs; LONGPHASE_HAPLOTAG writes CRAM (--cram) and PHASING_HAPLOTYPING joins on bai or crai. ASCAT and CRAMINO_POST now get the FASTA so they can decode CRAM. CLAIRSTO drops its intermediate haplotagged copy of the tumor reads (as large as the input) once done. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
samtools fastq -T '*' carried the basecaller's RG tag through minimap2 -y, next to the @rg added by -R, so every record had two RG tags. BAM keeps both (readers see the first), but CRAM keeps one read group per record: the last, the basecaller's, which has no @rg line in the aligned header. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Probe on the test-profile tasks (CRAM input, reference not reachable through the header, no EBI fallback): - Clair3 and ClairS-TO's Verdict alleleCounter bundle an htslib that cannot decode CRAM 3.1 slices; Clair3 then calls nothing and still exits 0. CRAM 3.0 decodes fine. - ClairS (internal mpileup), Severus and SAVANA's read counting open the reads without the reference; ClairS then writes empty VCFs and exits 0. embed_ref=1 makes every CRAM decodable without one, at ~1-3% size. longphase haplotag's own --cram cannot embed the reference, so with ext.cram it writes into a FIFO that samtools encodes (ext.args2). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The ASCAT in the module's image has no ref.fasta argument, so passing the FASTA fails ascat.prepareHTS outright; the CRAMs embed their reference, so alleleCounter decodes them without one. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
0.5.1 writes tmp_<sample_name>/, so the cleanup missed the 162 GiB haplotagged copy of BL1_ont's tumour reads. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
MINIMAP2_ALIGN, SAMTOOLS_MERGE and LONGPHASE_HAPLOTAG write CRAM 3.0 with the reference embedded only when aligned_format is cram; bam keeps BAM and .bai throughout. The merge/index joins take whichever output the format produced. The uBAM RG-tag fix applies to both formats. Tests: clair_only (the replicate samplesheet) runs --aligned_format bam and keeps its merged-reads MD5 check; the other pipeline tests expect CRAM. Snapshots not yet regenerated. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
All five pipeline tests pass on apptainer (default, consensus, deep_only, union on CRAM; clair_only on --aligned_format bam, snapshot unchanged). The only drift is bamfiles/ (.bam/.bai -> .cram/.crai) and two files whose header names the input file: samtools stats' command line and Severus' breakpoints_double.csv column names. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
|
samtools view --threads lays CRAM containers out differently on every run, so the published CRAM and CRAI bytes never repeat: PR #207 CI got a different set of md5s on each shard and attempt while flagstat, idxstats and stats (the read content) matched. Re-encoding one BAM three times with the module's image confirms it: with 4 threads the .crai changes between identical runs, with 1 thread only the path-derived file ID and @pg differ. Ignore */bamfiles/*.cram{,.crai} in tests/.nftignore (the names stay in stable_name) and drop the 10 md5 entries from the four CRAM snapshots. Verified on Mindwell: default (twice), consensus, deep_only and union pass without --update-snapshot. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The CHANGELOG said ASCAT now gets the FASTA, which 77bcf54 reverted, and quoted the CRAM 3.1 size savings and a 1-3 % embed_ref cost. The shipped 3.0 + embed_ref files are 33 % (ONT), 52 % (older PacBio) and 69 % (Revio) smaller, and the embedded reference costs 1-7 %. It now also says that a resumed run realigns, since -x RG changes MINIMAP2_ALIGN in both formats. The patched LONGPHASE_HAPLOTAG pipes longphase's output through samtools view, which its container has but its conda environment did not. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
clair_only is the only pipeline test with replicates, but it ran --aligned_format bam, so the default CRAM merge (SAMTOOLS_MERGE with the FASTA and version=3.0/embed_ref=1) had no test. It now runs the default. nft-bam's htsjdk decodes the embedded-reference CRAM without a FASTA, and sample4's merged reads MD5 is the one the BAM merge gave (88c8d3cf..., 7272 reads), so that snapshot entry is unchanged. deep_only takes --aligned_format bam instead. Its snapshot md5s are now identical to dev's, and the test drops from 79 to 11 min, since DeepVariant/DeepSomatic make_examples are slow on CRAM only on the thin chr19 test data. consensus and union still run both deep callers on CRAM. Snapshots regenerated with apptainer on Mindwell; the drift is bamfiles/, samtools .stats and Severus breakpoints_double.csv, whose headers name the input file. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Carry the aligned and haplotagged reads as CRAM instead of BAM, behind a new
--aligned_format cram|bam(defaultcram). On long reads this makes the aligned files a third (ONT) to two thirds (Revio) smaller, with identical records and identical calls. Full write-up with all measurements (shared on request): https://claude.ai/artifact/9gUVWBRFuLtv5kSFBgp3aRDraft: not ready for review yet (see "Before review" below).
Why these CRAM settings
Every tool downstream of alignment was re-run on CRAM with its own command and container, with the reference not reachable through the CRAM header. Two settings turned out to be required, because the failures without them are silent:
embed_ref=1). ClairS's internalsamtools mpileup, Severus and SAVANA's read counting open the reads without the FASTA. ClairS then writes empty VCFs and exits 0. It only appears to work on VSC because theUR:path in the header, which points into the minimap2 work dir, happens to be mounted inside the containers. An embedded reference adds 1-7 % to the file and makes each CRAM readable anywhere.Changes
MINIMAP2_ALIGN(nf-core, patched):bam_format = 'cram'makessamtools sortwrite-O cram --reference, with acramemit and a.craiindex.SAMTOOLS_MERGEgets the FASTA. The merge/index joins take whichever ofbam/cramandbai/craithe format produced, so merged replicates no longer depend on.out.bam.LONGPHASE_HAPLOTAG(nf-core, patched): withext.cram, longphase writes into a FIFO thatsamtoolsencodes, because longphase's own--cramcannot embed the reference. No temporary BAM reaches the disk.PHASING_HAPLOTYPINGjoins onbaiorcrai.CRAMINO_POSTgets--reference. ASCAT keeps no FASTA: the module's image has noref.fastaargument, so passing one failsascat.prepareHTS.--aligned_format(schema,nextflow.config,conf/modules.configclosures).bamkeeps today's files and.baiindexes.RGtags fixed: every aligned record carried twoRGtags, the basecaller's (viasamtools fastq -T '*'and-y) and the pipeline's-Rread group. CRAM keeps one per record and kept the basecaller's, which has no@RGline.samtools reset -x RGnow drops the uBAM's copy.CLAIRSTOcleanup: deletestmp_<sample>/phasing_output/phased_bam_outputwhen done, the full haplotagged copy of the tumour reads (157-162 GiB for one ONT sample).docs/usage.md,docs/output.md,CHANGELOG.md.Validation
test_sheet_2, sample4 merged),--aligned_format bamvscram: 185 files identical, merged reads included.read_ids.csv, which also differs between two BAM runs.deep_onlyruns--aligned_format bam; the other four run the CRAM default. That includesclair_only, the only test with replicates, so sample4's two replicates go through the CRAM merge. nft-bam reads the merged CRAM without a FASTA (embedded reference), and its reads MD5 equals the one the BAM merge gave (88c8d3cf…, 7272 reads), so the snapshot entry is unchanged.deep_onlyis also the test DeepVariant/DeepSomatic slow down most on the thin CRAM test data (79 → 11 min); on BAM it reproduces dev's snapshot md5s exactly. Snapshots were regenerated; the only drift isbamfiles/(.bam/.bai↔.cram/.crai) and two files whose header names the input file (samtools.stats, Severusbreakpoints_double.csv). CRAM/CRAI md5s are not snapshotted, because multi-threadedsamtoolslays the containers out differently on every run.Before review
make_examplesare 13-14× slower on CRAM. On real data there is no slowdown: 88 s BAM vs 87 s CRAM on chr20:10-15 Mb, same 111 224 examples.deep_onlyruns on BAM, andconsensus/unionstill run both deep callers on CRAM.-profile condaneedssamtoolsin the patchedlongphase/haplotagenvironment for the FIFO; it has been added, but conda is not in the CI matrix.PR checklist
nf-core pipelines lint).docs/usage.mdis updated.docs/output.mdis updated.CHANGELOG.mdis updated.🤖 Generated with Claude Code