Repository navigation
Carry aligned reads as CRAM, with --aligned_format cram|bam - #207
ljwharbers wants to merge 24 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>
|
@claude review |
![claude[bot]](https://avatars.githubusercontent.com/in/1236702?s=60&v=4){kind=small}
| tuple val(meta), path("*.paf") , optional: true, emit: paf | ||
| tuple val(meta), path("*.bam") , optional: true, emit: bam | ||
| tuple val(meta), path("*.bam.${bam_index_extension}"), optional: true, emit: index | ||
| tuple val(meta), path("*.cram") , optional: true, emit: cram |
There was a problem hiding this comment.
The patch adds a cram emit, widens the index glob to *.{bam,cram}.${bam_index_extension}, and makes bam_format accept 'cram'/'bam'. But meta.yml is not changed (the patch says "'meta.yml' is unchanged"). It still lists bam_format as a boolean ("output should be in BAM format"), has no cram output, and documents index as *.bam.${bam_index_extension}. REVIEW.md asks for meta.yml to be updated with new inputs/outputs, and nf-core modules lint compares the emits against meta.yml. Please add the cram output, the new index pattern and the string values of bam_format, then regenerate minimap2-align.diff.
| ext.args2 = { params.aligned_format == 'cram' ? '--output-fmt-option version=3.0 --output-fmt-option embed_ref=1' : '' } | ||
| // Drop the uBAM's own RG tag: -T '*' would carry it through -y next to the @RG from -R, and a record | ||
| // with two RG tags keeps only the last one in CRAM -- the basecaller's, which has no @RG line here | ||
| ext.args3 = "-x RG" |
There was a problem hiding this comment.
No test exercises the duplicate-RG fix. The CHANGELOG says -x RG "alters MINIMAP2_ALIGN in both formats". Yet deep_only (--aligned_format bam) keeps every bamfiles/*.bam md5 byte-identical to dev (only the timestamp moved in deep_only.nf.test.snap). So the test uBAMs carry no RG aux tag, and this line is a no-op in the test suite. A regression here (e.g. a samtools reset version without -x, or the uBAM RG coming back through -T '*') would not be caught. Please consider a test uBAM with a basecaller RG tag, or an assertion that every aligned record's RG matches an @RG header line.
| loci_files, | ||
| [], | ||
| [], | ||
| [], // no fasta: the image's ASCAT has no ref.fasta argument; the CRAMs embed their reference |
There was a problem hiding this comment.
ASCAT now reads CRAM with no FASTA and relies on the embedded reference. No nf-test covers this: conf/test.config sets skip_ascat, skip_wakhan and skip_modkit for every pipeline test. So ASCAT (alleleCounter), Wakhan and modkit on CRAM are only validated by the one full-sample run in the PR description (BL1_ont, tumour-only, CHM13). ASCAT on a paired tumour/normal CRAM and GRCh38 is never run. Please run these steps on CRAM (e.g. the Slurm test array or a paired real sample) before merging. If alleleCounter ever falls back to the UR: path, it fails silently, as the PR says ClairS does.
| SAMTOOLS_MERGE( | ||
| ch_aligned_split.multiple, | ||
| [[],[],[],[]] | ||
| ch_fasta.join(ch_fai).map { meta, fasta, fai -> [meta, fasta, fai, []] }.first() |
There was a problem hiding this comment.
PR CI (small tag) never runs this path. Only tests/default.nf.test is tagged small, and its samplesheet has no replicates. The CRAM merge (FASTA passed to SAMTOOLS_MERGE, the out.bam.mix(out.cram) / bai.mix(crai) joins) is only covered by clair_only (extended). The whole --aligned_format bam path, including the true→'bam' change to MINIMAP2_ALIGN and the ext.cram == false branch of LONGPHASE_HAPLOTAG, is only covered by deep_only (extended). The PR says these passed on Mindwell. Still, a later change could break either path with CI staying green, so it may be worth tagging one of them small, or at least rerunning the extended/Slurm array before merge.
|
Automated review: 4 inline findings.
A human reviewer needs to run the extended tests / Slurm array on this branch, plus ASCAT/Wakhan/modkit on a paired CRAM sample. This review read the diff only and did not run any tests. |
…clair_only
minimap2/align's meta.yml now describes the cram emit, the {bam,cram}
index pattern and the string values of bam_format. nf-core modules lint
still reports correct_meta_outputs for this module: for a patched module it
compares the reverse-patched (upstream) meta.yml with the patched main.nf,
so it failed the same way before this change.
clair_only checks every output CRAM: each record's RG must be declared in
the file's header (CRAM keeps one RG per record, so an undeclared second
one would replace the pipeline's), and sample4's merged header must keep
both replicates' read groups. The test uBAMs carry no record-level RG tag,
so this guards merge and haplotag but cannot yet catch a regression of
MINIMAP2_ALIGN's -x RG; that needs an RG-tagged test input.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
clair_only now reads test_sheet_2_rgtagged.csv (IntGenomicsLab/test-datasets#5, must be merged first): test_sheet_2 with sample3's uBAM carrying dorado's RG on every record, as real dorado output does. None of the existing test uBAMs has a record-level RG tag, so the read-group assertion could not see a MINIMAP2_ALIGN that keeps it. Checked on Mindwell from dd3127f: with -x RG the test passes and the snapshot is byte-identical to the one from test_sheet_2; with ext.args3 = "" the read-group assertion fails, with all 3122 sample3 records carrying the undeclared dorado RG, while no other sample or snapshot entry changes. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
On the MT_SCCNA C8 PacBio pair, longphase 2.0.1 modcall called 28 % fewer methylation sites from the CRAM run than from the BAM run (normal 64333 -> 46927, tumour 26972 -> 18886) with identical reads. On chr20 of the normal: BAM 2261 sites twice, the CRAM decoded back to BAM 2261 twice, the CRAM itself 2019/1871 and a slice of the pipeline's CRAM 1761/1697, with MM/ML byte-identical. 2.0.1 walks a sam_itr_querys iterator with sam_itr_multi_next; longphase 2.0.2 replaces that with sam_itr_next and fixes an uninitialised field that made modcall non-deterministic. With 2.0.2 BAM and CRAM give the same 2261 sites, same md5, on repeated runs. longphase phase (still 2.0.1, same iterator pattern) is unaffected: its chr20 output is byte-identical for BAM and CRAM with and without a mod file, so only the local LONGPHASE_MODCALL module changes. It needs only longphase, so it now uses the plain 2.0.2 biocontainer. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
All six pipeline tests pass on Mindwell (apptainer) with 382ba38. The only change in every snapshot is LONGPHASE_MODCALL's longphase version (2.0.1 -> 2.0.2) plus the timestamp; no output md5 moves. So the chr19 test data is too small to trigger 2.0.1's CRAM misread, which is why the tests never showed it. clair_only was run on a local copy of test_sheet_2_rgtagged.csv, as IntGenomicsLab/test-datasets#5 is not merged yet. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
headered_H2009_ONT_UL_chr19_downsampled.bam (+ .bai) itself now carries RG:Z:<its header's dorado @rg ID> on every record, and the separate _rgtagged copy and test_sheet_2_rgtagged.csv are gone. Every lrsomatic pipeline test that reads this file (sample1 tumour and sample3 in samplesheet_lr-somatic.csv, sample1/3/5 in test_sheet_2.csv) now gets an RG-tagged uBAM, as real dorado output is. Merge only once IntGenomicsLab/lrsomatic#207 is in dev: without its `samtools reset -x RG`, the aligned records carry two RG tags and dev's BAM snapshots change. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
| 'https://community-cr-prod.seqera.io/docker/registry/v2/blobs/sha256/83/83fce1d397cf71705cc096fc0e0e52f7013bdd471ef68ee53ae765688e5c439c/data': | ||
| 'community.wave.seqera.io/library/longphase_samtools:8c61296cae7a5fc0' }" | ||
| 'https://depot.galaxyproject.org/singularity/longphase:2.0.2--h4e109e1_0': | ||
| 'quay.io/biocontainers/longphase:2.0.2--h4e109e1_0' }" |
There was a problem hiding this comment.
This breaks the container convention in REVIEW.md. Biocontainers should use the depot.galaxyproject.org / biocontainers/ pair, not quay.io/biocontainers/..., because the registry comes from docker.registry:
| 'quay.io/biocontainers/longphase:2.0.2--h4e109e1_0' }" | |
| 'biocontainers/longphase:2.0.2--h4e109e1_0' }" |
Two smaller problems:
environment.ymlstill pinsbioconda::htslib=1.23.1next tolongphase=2.0.2. This plain biocontainer ships whatever htslib theh4e109e1_0build resolved, so-profile condaand the containers may run different htslib versions. That matters because this bump exists to fix an htslib iterator problem on CRAM. Please either drop the pin or confirm the image ships 1.23.1.- The comment above stops mid-sentence ("so the plain biocontainer").
| input = "https://raw.githubusercontent.com/IntGenomicsLab/test-datasets/refs/heads/main/samplesheets/test_sheet_2.csv" | ||
| // test_sheet_2 with sample3's uBAM RG-tagged on every record, as dorado output is, so the | ||
| // "Read groups" assertion below catches a MINIMAP2_ALIGN that keeps the uBAM's RG | ||
| input = "https://raw.githubusercontent.com/IntGenomicsLab/test-datasets/refs/heads/main/samplesheets/test_sheet_2_rgtagged.csv" |
There was a problem hiding this comment.
PR CI runs only tests tagged small (.github/workflows/nf-test.yml), and that is just default.nf.test. This test is extended, so CI never runs the parts this PR adds here:
- the CRAM replicate merge (
SAMTOOLS_MERGEwith the new FASTA input, thenSAMTOOLS_INDEX_MERGEjoined on.crai); - the read-group assertion that guards the new
samtools reset -x RG; - the
sample4_merged_readscheck.
Also, test_sheet_2_rgtagged.csv doesn't exist on test-datasets main until IntGenomicsLab/test-datasets#5 merges. Until then this test fails at input.
Please either tag a test that covers the merge and RG paths as small, or run the Slurm/extended test array (clair_only in particular) after test-datasets#5 is merged and before this lands.
|
Automated review: 3 findings, posted as inline comments.
A human reviewer needs to run, on the final commit and after test-datasets#5 is merged:
|
| input = "https://raw.githubusercontent.com/IntGenomicsLab/test-datasets/refs/heads/main/samplesheets/test_sheet_2.csv" | ||
| // test_sheet_2 with sample3's uBAM RG-tagged on every record, as dorado output is, so the | ||
| // "Read groups" assertion below catches a MINIMAP2_ALIGN that keeps the uBAM's RG | ||
| input = "https://raw.githubusercontent.com/IntGenomicsLab/test-datasets/refs/heads/main/samplesheets/test_sheet_2_rgtagged.csv" |
There was a problem hiding this comment.
This points at samplesheets/test_sheet_2_rgtagged.csv, but IntGenomicsLab/test-datasets#5 doesn't add that file any more. Its body says: "An earlier version of this PR added the tagged file as a copy plus a test_sheet_2_rgtagged.csv; replaced by tagging the original in place". So the PR now only modifies data/aligned/headered_H2009_ONT_UL_chr19_downsampled.bam(+.bai). As written, clair_only will fail to fetch its samplesheet, whichever order the two PRs merge in.
The merge order also contradicts itself. This PR's "Before review" says to merge test-datasets#5 first. test-datasets#5 says to merge it only right after #207, because tagging in place changes dev's aligned BAMs and would break every open PR's default test.
Suggest switching back to test_sheet_2.csv (sample1/3/5 use the file that gets tagged in place) and fixing the merge order in the description. Until #5 is merged, the read-group assertion below runs on untagged data, so it can't catch a missing -x RG.
| somatic_var_keep = 'deepsomatic' | ||
| // The one test on BAM (the others run the CRAM default); DeepVariant/DeepSomatic make_examples are | ||
| // 13-14x slower on CRAM on the thin chr19 test data only, so this also keeps the test short | ||
| aligned_format = 'bam' |
There was a problem hiding this comment.
PR CI runs only small-tagged tests, and the only one is default.nf.test, on the CRAM default. This test and clair_only are tagged extended, so CI exercises none of these:
--aligned_format bam(only set here). That includes thebam_format = 'bam'string replacing the oldtrueinMINIMAP2_ALIGN, the''ext.args2 branches and the non-FIFO haplotag path.- The CRAM replicate path:
SAMTOOLS_MERGEwith the reference,SAMTOOLS_INDEX_MERGE.craijoins, and sample4 inclair_only. - The new read-group assertion.
Per REVIEW.md, please run the extended/Slurm test array on the final commit (after the test-datasets change lands), or tag one of these small.
| def cram_out = task.ext.cram ?: false | ||
| def cram_start = cram_out ? """ | ||
| mkfifo ${prefix}.bam | ||
| samtools view --threads ${task.cpus} -C -T ${fasta} ${args2} -o ${prefix}.cram ${prefix}.bam & |
There was a problem hiding this comment.
The FIFO encoder runs at the same time as longphase haplotag --threads $task.cpus, and both get task.cpus threads. That's about 2× the process_medium allocation, and CRAM encoding with embed_ref=1 is CPU-heavy. On schedulers that enforce the core count this means contention or the job getting killed. Consider splitting the allocation, e.g. giving samtools --threads ${Math.max(1, task.cpus.intdiv(2))} and longphase the rest.
|
Automated review: 3 findings (inline).
To check by running tests: run the extended/Slurm array on the final commit, after the test-datasets change lands. Confirm that |
Workflow-level concurrency is evaluated before the job's if, so every comment on a PR joined claude-review-<N> and cancelled the running review. Claude's own summary comment did this at the end of every review, so each run on IntGenomicsLab#202 and IntGenomicsLab#207 ended as cancelled (red) even though its comments were posted. A human comment mid-review would kill it outright. Comments that are not '@claude review' now get a per-run group. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Workflow-level concurrency is evaluated before the job's if, so every comment on a PR joined claude-review-<N> and cancelled the running review. Claude's own summary comment did this at the end of every review, so each run on IntGenomicsLab#202 and IntGenomicsLab#207 ended as cancelled (red) even though its comments were posted. A human comment mid-review would kill it outright. Comments that are not '@claude review' now get a per-run group. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…in place test_sheet_2_rgtagged.csv was dropped from test-datasets#5, so the URL 404s. The RG assertion holds either way and starts biting once #5 is merged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Conflicts with dev's module updates (#212): - SAMTOOLS_MERGE keeps the reference tuple (needed for CRAM) and takes dev's new index_format argument (''). - samtools merge/index now emit one `.index` channel for .bai/.crai/.csi, so the bai.mix(crai) joins become .join(...out.index). - longphase/haplotag moves to dev's git_sha and keeps the CRAM FIFO patch; both it and the minimap2/align patch were regenerated with `nf-core modules patch` (haplotag diff byte-identical). - tests/.nftignore and CHANGELOG.md keep both sides. - Pipeline snapshots take dev's version here; they are regenerated in the next commit. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
All six tests pass on Mindwell (apptainer, Slurm array 11667471). Drift vs dev's snapshots is the CRAM switch only: bamfiles/ lists .cram/.crai instead of .bam/.bai, the samtools .stats and Severus breakpoints_double.csv md5s move (their headers name the input file), and longphase reports 2.0.2 (modcall). deep_only runs --aligned_format bam, so only its version line moves. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
|
@claude review |
1 similar comment
|
@claude review |
| somatic_var_keep = 'deepsomatic' | ||
| // The one test on BAM (the others run the CRAM default); DeepVariant/DeepSomatic make_examples are | ||
| // 13-14x slower on CRAM on the thin chr19 test data only, so this also keeps the test short | ||
| aligned_format = 'bam' |
There was a problem hiding this comment.
Test coverage: the --aligned_format bam path is only exercised by an extended test.
deep_only.nf.test is tagged extended, and PR CI only runs small tests (.github/workflows/nf-test.yml). That leaves CI with the CRAM single-replicate path only (default.nf.test). It never checks the BAM branch through the patched MINIMAP2_ALIGN (string bam_format, out.bam.mix(out.cram)), SAMTOOLS_MERGE (now with --reference), or haplotag with ext.cram = false.
Per REVIEW.md, please either add a small case (or a stub test) with aligned_format = 'bam', or confirm that the Slurm test array has been run on the final head commit.
| def ids = reader.getFileHeader().getReadGroups()*.getId() as Set | ||
| def undeclared = reader.iterator().findAll { !(it.getAttribute('RG') in ids) }.size() | ||
| reader.close() | ||
| assert undeclared == 0 : "${name}.cram: ${undeclared} records carry an RG missing from the header ${ids}" |
There was a problem hiding this comment.
This read-group check can pass without testing anything, and CI never runs it.
- The test still reads
test_sheet_2.csv(line 13). The PR description says it readstest_sheet_2_rgtagged.csv, but commit "Read test_sheet_2 again…" moved it back and relies on Tag the H2009 ONT test uBAM with its read group (merge right after lrsomatic#207) test-datasets#5 re-tagging that uBAM in place. Until Add mosdepth module #5 is merged, the H2009 uBAM has no doradoRG. Soundeclared == 0holds even with-x RGremoved fromMINIMAP2_ALIGN. The comment above ("carries dorado's RG on every record") is only true after Add mosdepth module #5. clair_onlyis taggedextended, so PR CI (onlysmalltests) never runs this assertion or the CRAMSAMTOOLS_MERGEpath (sample4) at all. The only safeguard for the-x RGfix and the CRAM merge is therefore a manual Slurm run.
Please merge test-datasets#5 first. Then either run the Slurm array on the final commit, or move a minimal RG / merge check into a small test.
|
Automated review: 2 findings, both about test coverage. I found no code bugs in the diff.
For a human reviewer: I could not run the pipeline or nf-test. Please run the Slurm test array on the final head commit after test-datasets#5 is merged, covering |
Brings in DMR (#204), CoRAL/ecDNA (#205), the FIBERSEQ subworkflow, the REJOIN norm before LongPhase, the Severus whitelist image and Severus on haplotagged reads, among others. Conflicts: - CHANGELOG.md: kept both sides' "Changed" entries (#207 first, then dev's #204/#205 entries). - tests/*.nf.test.snap: resolved per entry, not per side. Timestamps take dev's position after "meta". bamfiles/ md5s stay off the CRAM tests (CRAMs are nftignored); dev's .bam/.bai md5s for default and clair_only sample5 are dropped. samtools .stats keep the branch values (the header names the CRAM). WhatsHap stats for clair_only sample4/sample5 take dev's values (REJOIN changes the phasing). Severus breakpoints_double.csv in default changed on both sides (dev: Severus now reads haplotagged reads; branch: the header names the CRAM), so the branch values are kept as placeholders and are stale until the snapshots are regenerated. CRAM check of the dev code: no change needed. CoRAL (seed, reconstruct, plot) reads the tumour minimap2 CRAM through pysam.AlignmentFile(path, "rb"), which opens an embed_ref CRAM with no reference reachable (checked with pysam 0.23.3). Severus, Wakhan and SAVANA now all read the haplotagged CRAM+CRAI from PHASING_HAPLOTYPING; the patched Severus fork adds no new read access. DMR, FIBERSEQ, REJOIN and the tumour-only DS_VERDICT changes do not touch aligned reads. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Slurm array 11750729 on 955c79e (apptainer, --update-snapshot): all six pipeline tests pass. The only output drift is Severus breakpoints_double.csv, now that Severus reads the haplotagged reads (dev) from CRAM (this branch). chm13 needed no update. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
|
@claude review |
| }, | ||
|
|
||
| // ── sample4: merged BAM content ────────────────────────────── | ||
| // ── Read groups ────────────────────────────────────────────── |
There was a problem hiding this comment.
PR CI never runs this test. clair_only is tagged extended, and .github/workflows/nf-test.yml only runs tests tagged small on pull requests. So these checks never run in CI:
- this read-group assertion, which the PR description calls the only check that catches
MINIMAP2_ALIGNwithoutsamtools reset -x RG; - the CRAM
SAMTOOLS_MERGEpath (sample4's replicates) and its.first()-wrapped FASTA/FAI join.
default is the only small test that actually runs tasks (sv_annotation_compose is -preview only). Its samplesheet has no replicates, and it has no RG assertion.
A future regression in either area would pass PR CI. Either move the RG and merged-reads assertions into a small test, or tag this one small.
There was a problem hiding this comment.
I moved the read-group assertion into the default test, which is small. Its H2009 uBAM puts dorado's RG on every sample1-tumour and sample3 record, so the check is live there. In a negative control with ext.args3 = "" it is the only assertion that fails (sample1_tumor.cram: 3122 records carry an RG missing from the header), and the snapshot stays the same. The replicate merge now runs on every PR through the new small -stub test on the extended samplesheet: it checks that SAMTOOLS_MERGE runs for sample4 in both formats and checks the published files. The merged reads' content (getReadsMD5) stays in clair_only, because a real replicate run would add about 12 minutes to every PR. Fixed in 4df8743.
| somatic_var_keep = 'deepsomatic' | ||
| // The one test on BAM (the others run the CRAM default); DeepVariant/DeepSomatic make_examples are | ||
| // 13-14x slower on CRAM on the thin chr19 test data only, so this also keeps the test short | ||
| aligned_format = 'bam' |
There was a problem hiding this comment.
This is the only test of --aligned_format bam, and it is tagged extended, so PR CI never runs it.
The BAM branches added here get no CI coverage:
bam_format == 'cram' ? … : 'bam'in the patchedMINIMAP2_ALIGNMINIMAP2_ALIGN.out.bam.mix(...cram)- the
ext.cramfalse path inLONGPHASE_HAPLOTAG SAMTOOLS_MERGEwith a reference on BAM input
Consider a small test (a -stub run is enough) that sets aligned_format = 'bam' and checks that bamfiles/*.bam and *.bam.bai are published.
There was a problem hiding this comment.
I added tests/aligned_format.nf.test (small). It runs -stub twice on the extended samplesheet, once on the CRAM default and once with aligned_format = 'bam'. Each run asserts that bamfiles/ holds exactly .bam+.bam.bai (or .cram+.crai) for all seven alignments, and that SAMTOOLS_MERGE ran for sample4. That covers MINIMAP2_ALIGN's output format, the bam/cram mix, the merge with its reference, and haplotag's ext.cram false path, at about 1.5 minutes per case. The pipeline needed two stub fixes to run under -stub: METAEXTRACT reported no model, and the patched LONGPHASE_PHASE stub indexed VCFs it never wrote. Fixed in 4df8743.
| "breakpoint_clusters.tsv:md5,d36a70de292ee130ef30da4a58bced18", | ||
| "breakpoint_clusters_list.tsv:md5,0c0ce62e329f8de492487e8414c30a50", | ||
| "breakpoints_double.csv:md5,003e4e80ec4233315762233e5f1ec0d8", | ||
| "breakpoints_double.csv:md5,3dcda68bb247950df3ce73812a087863", |
There was a problem hiding this comment.
Unexplained snapshot drift on the BAM path. deep_only runs --aligned_format bam, and the PR description says it "reproduces dev's snapshot md5s exactly". Yet all three Severus breakpoints_double.csv md5s change here.
The PR explains this drift as "the header names the input file". That explanation only holds for the CRAM runs (.bam → .cram). In BAM mode the haplotagged file names are unchanged, and the samtools .stats md5s in this snapshot do not move.
So something else changed Severus's output on BAM. The likely cause is the new samtools reset -x RG in MINIMAP2_ALIGN, or the longphase 2.0.2 modcall changing the phasing that Severus reads through HP tags. Please name the cause in the PR before accepting these md5s.
There was a problem hiding this comment.
This drift comes from dev, not from this PR. #223 moved Severus onto the haplotagged reads, whose HP tags depend on the caller, and dev re-recorded only the default snapshot (b7e7382), so the extended snapshots went stale. deep_only on unmodified dev 2b1d510 (Slurm 11748329, committed in #224 as 7c6e9fe) gives exactly these three md5s. Apart from the timestamp, this snapshot differs from #224's only in the two tumour BAM md5s (samtools reset -x RG) and longphase 2.0.2's version line. The "header names the input file" explanation only covers the CRAM tests (Severus' columns are <sample>_tumor.cram_0_…); I've corrected that in the PR description.
| loci_files, | ||
| [], | ||
| [], | ||
| [], // no fasta: the image's ASCAT has no ref.fasta argument; the CRAMs embed their reference |
There was a problem hiding this comment.
CoRAL now gets CRAM by default, and nothing has checked it.
ch_ecdna_tumor_bamis built from this ASCAT input, soCORAL_RECONSTRUCT --lr-bamnow receives the tumour CRAM.skip_coral = falseby default.conf/test.configsetsskip_coral(andskip_ascat,skip_wakhan,skip_modkit), so no nf-test runs it.CORAL_RECONSTRUCThaserrorStrategy ... 'ignore'. If CoRAL cannot read CRAM, the run still succeeds and the ecDNA results are silently empty.
The PR description lists ASCAT, Wakhan and modkit as checked on the real C8/BL1 samples, but not CoRAL. It also says that whichever of #205/#207 merges second must check CoRAL on CRAM. #205 is already in this branch, so that check falls to this PR.
Please run CoRAL (seed/reconstruct/cycle) on a CRAM before merging. If it has not been checked, force BAM for that path.
There was a problem hiding this comment.
I checked it on real data: BL1_ont (CHM13, ONT) with its real ASCAT chr21 gain, CRAM 3.0 with embed_ref=1 from this branch against its BAM. I ran the module commands with the default args in a contained container with no reachable reference. seed, reconstruct, cycle_all and plot_all all succeed. The seeds, graph, cycles, summary, plots and CoRAL's full-file chimeric-read pickle (610,280 reads) are byte-identical to the BAM run, and CRAM was faster (2040 s vs 2965 s for the read pass). One caveat: this amplicon has no discordant edges, so no SV edge was exercised. CoRAL keeps the CRAM; details are in the PR description.
|
Automated review: 4 findings (inline).
For a human reviewer to check by running (not verified here):
|
…tests PR CI runs only `small` tests, and its one real run (default) is CRAM without replicates, so the BAM path, the replicate merge and the read group check in clair_only never ran on a PR (review on #207). - tests/aligned_format.nf.test (small): two -stub runs on the extended samplesheet, CRAM and --aligned_format bam. Each asserts that SAMTOOLS_MERGE ran for sample4 only and that bamfiles/ holds exactly the expected .cram/.crai or .bam/.bai for all seven alignments. 98 s + 78 s on Mindwell. - tests/default.nf.test: the read-group assertion from clair_only. Its H2009 uBAM (sample1 tumour, sample3) carries dorado's RG on every record. With MINIMAP2_ALIGN ext.args3 = "" it is the only one of the seven assertions that fails ("sample1_tumor.cram: 3122 records carry an RG missing from the header"); with the fix the snapshot does not change. The pipeline could not run with -stub before. METAEXTRACT's stub reported no basecalling model, so no Clair3 model resolved and the paired germline arm never started. The patched LONGPHASE_PHASE stub indexed VCFs it never wrote. Both stubs now write their outputs, and the whole pipeline runs with no failed task. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
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.LONGPHASE_MODCALL(local) runs longphase 2.0.2 instead of 2.0.1. 2.0.1 walks its region iterator withsam_itr_multi_next; on CRAM that returns fewer reads, differently every run, so modcall called about 28 % fewer methylation sites on CRAM than on BAM, and the phasing and every haplotype-resolved output moved with it. 2.0.2 swaps insam_itr_nextand fixes an uninitialised field that made modcall non-deterministic.longphase phase(still 2.0.1, same pattern) reads CRAM correctly: its chr20 output is byte-identical for BAM and CRAM, with and without a mod file. This also helpsdev: 2.0.1 loses about 7 % of modcall sites on BAM too (see Validation).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.--aligned_format bamvscram, same commit:BPS=differs). Severusbreakpoint_clusters*.tsvlist tumour/normal in a different order per cluster, andread_ids.csvdiffers between BAM runs too.chm13); the bump changes only the longphase version in the snapshots.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 the md5s ofdeep_onlyrun on unmodified dev 2b1d510 (Slurm 11748329, recorded in Template update for nf-core/tools v4.2.1 #224), apart from the tumour BAMs that-x RGchanges. Snapshots were regenerated after the dev merge. The CRAM-vs-BAM drift isbamfiles/(.bam/.bai↔.cram/.crai) and two files whose header names the input file (samtools.stats, Severusbreakpoints_double.csv, whose columns are<sample>_tumor.cram_0_…). The Severus md5s also moved in all five tests for a reason unrelated to this PR: Severus fix #223 moved Severus onto the haplotagged reads, whose HP tags depend on the caller, and dev re-recorded onlydefault(b7e7382), so dev's extended snapshots were stale. CRAM/CRAI md5s are not snapshotted, because multi-threadedsamtoolslays the containers out differently on every run.Before review
RGon every record, andclair_onlyreads the plaintest_sheet_2.csv. Its read-group assertion is the only check that catches aMINIMAP2_ALIGNwithoutsamtools reset -x RG. Re-run on Mindwell against the tagged data (7e5aa48): all six pipeline tests pass with the committed snapshots, and the assertion now has records to check (with-x RGremoved, all 3122 sample3 records carried the undeclared dorado RG and only this assertion failed).defaulttest fails until this PR merges. dev has no-x RG, so its aligned BAMs now carry the tag andsample1_tumor.bam/sample3_tumor.bam(+.bai) md5s no longer match dev's snapshot (checked on dev 99d00fc; nothing else drifts). Other PRs against dev hit the same failure until they merge dev after this one.WAKHAN. The module gets the FASTA without its.fai, and Wakhan's parallelsamtools mpileup -fstart about 16 s before the index is written. In all four C8 runs (both formats, before and after the longphase bump) one or two chromosomes were piled up with referenceN(which ones varies), and once a contig was missing from the half-written index altogether, so the task failed and was retried ([E::fai_read] Could not understand FASTA index, reproduced outside the pipeline). The final copy-number calls were not affected here. Fix: stagech_faiintoWAKHAN; separate PR againstdev.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.embed_ref=1from this branch vs its BAM, module command lines with default args in a contained container with no reachable reference. seed, reconstruct,cycle_allandplot_allall succeed. Seeds, graph, cycles, summary, plots and CoRAL's chimeric-read pickle (610 280 reads) are byte-identical; the CRAM read pass was faster (2040 s vs 2965 s). That amplicon has no discordant edges, so no SV edge was exercised.small):tests/aligned_format.nf.testruns-stubon the replicate samplesheet in both formats and checks the publishedbamfiles/and thatSAMTOOLS_MERGEran for sample4.defaultnow carries the read-group assertion (the H2009 uBAM has dorado's RG on every record; with-x RGremoved it is the only assertion that fails). This needed two stub fixes: METAEXTRACT's stub reported no basecalling model, and the patched LONGPHASE_PHASE stub indexed VCFs it never wrote. The merged reads' content check stays inclair_only(extended).-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