Repository navigation
Conversation
`dnascope` now accepts a population VCF and transfers its annotations onto the raw DNAscope VCF before DNAModelApply, the way `dnascope-longread`, `dnascope-hybrid` and the pangenome pipelines already do. The existing `TransferApplyStage` call gains a `TransferSpec`; with `--gvcf` the transfer runs on the gVCF, as in `dnascope-pangenome`. bcftools >= 1.22 is required only when `--pop_vcf` is supplied. A new `validate_bundle` reads the bundle's `bundle_info.json` and follows the `dnascope-longread` rules: a bundle with a `SentieonVcfID` requires a `--pop_vcf` whose header carries the same ID (the hidden `--skip_pop_vcf_id_check` skips the comparison), and a bundle without one rejects `--pop_vcf`. The reference index and the pop VCF contigs are read, and the genome sharded, only when `--pop_vcf` is supplied. Because `validate()` now always reads the model bundle, the test fixtures that stood in an empty file for the bundle write an empty ar archive instead -- a bundle without a `bundle_info.json`. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
`validate_bundle` now reads the bundle's `platform` from its `bundle_info.json`, and when it is Ultima (in any case) the DNAscope driver call gets `--read_filter UltimaReadFilter`, as in `dnascope-pangenome`. The small-variant and BND passes share that call, so SV calling is filtered too; duplicate marking, metrics and DNAModelApply are not. A bundle without a `platform`, or without a `bundle_info.json`, adds no filter. The `--pop_vcf` tests move to `test_dnascope_bundle.py`, which now covers both values read from the bundle. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
`indel2cnv.py` and `combine_cnv.py` move to the combine_sv_cnv 2.0.1 release. `indel2cnv.py` now finds an indel's repeat period from a k-mer gap histogram and no longer uses `find_peaks`, so `_indel2cnv_utils.py` and its test are removed. `combine_cnv.py` follows CNVscope's ploidy, marks CNVscope gains as well as losses `SVdup` when a converted SV call describes the same event, filters long converted gains that CNVscope's raw segmentation does not support (`NoDepth`) and, for single-end CNVscope models, short CNVscope losses that no converted loss overlaps (`NoSV`). The `combine-cnv` job passes the new arguments. `--preset` is `SE` for Ultima bundles, which carry a single-end CNVscope model, and `PE` otherwise: `BasePangenome.cnv_combine_preset()` returns `PE` and `SentieonPangenome` overrides it. `--sex` and `--par` are the arguments CNVscope was called with, carried on `CNVResult.sex` / `par`, so both tools use the same ploidy; neither is passed when the sample sex is unknown. `--raw` is the CNVscope output before CNVModelApply, which `combine-cnv` already depends on through `cnv-model-apply`. `cmd_pyexec_combine_cnv` takes `raw_vcf`, `preset`, `sex` and `par` as keyword-only arguments. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
`combine_cnv.py` 2.0.1 required an `HMMCN` header on `--raw`, which the CNVscope output from current bundles does not have, so every `combine-cnv` job failed. 2.1.0 accepts that output: the `NoDepth` check then takes its gain evidence from the CNVModelApply PASS calls and logs a warning, and the `##CombineCNVPreset` header names the evidence used. `--raw` and `--cnv` are checked for CNVscope and CNVModelApply provenance instead, and a warning is logged when no sex is known. Output with HMMCN-bearing CNVscope output is unchanged. `indel2cnv.py` changes only its version constant. The command line is the same as in 2.0.1, so the pipeline is unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…reset `dnascope-pangenome` and `hybrid-pangenome` no longer require `dnascope.model` in the model bundle under `--skip_small_variants`. The model is only passed to the small-variant DNAscope call, and `--segdup_caller`, which reads it in `hybrid-pangenome`, is already rejected with `--skip_small_variants`. `hybrid-pangenome` reads the bundle's `shortReadPlatform` into `sr_platform` and, for Ultima short reads, passes `combine_cnv.py` the single-end `SE` preset, as `dnascope-pangenome` does for an Ultima `platform`. A bundle without the key keeps `PE`. The new bundle tests also cover `bwa.model` being required only with FASTQ input. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
KMC fails when run with more than 128 threads. `cmd_kmc` and
`_kmc_stdin_cmd` now pass `-t{min(threads, KMC_MAX_THREADS)}`, so the
`kmc` job in `dnascope-pangenome` and the `kmc` and `extract-kmc` jobs in
`hybrid-pangenome` use at most 128 threads. Only the KMC thread count is
capped: `samtools fasta` and the ReadWriter driver in the same jobs, and
the jobs' scheduler thread counts, still follow `--cores`.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
`GVCFtyper` and `GVCFtyperStage` take an optional `dbsnp`, emitted as `--dbsnp` only when it is set, and `dnascope` and `dnascope-pangenome` pass their `--dbsnp`. With `--gvcf`, the VCF that GVCFtyper genotypes from the gVCF now carries dbSNP IDs, as the VCF-mode output does. Without `--dbsnp` the commands are unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
With `--gvcf`, `dnascope-hybrid` ran both calling passes in gVCF mode but sent their records through the VCF post-processing into `--output_vcf`, where the `bcftools view -a` in `final-norm` removed `<NON_REF>`; no `.g.vcf.gz` was written. The pipeline now follows the gVCF-then-GVCFtyper pattern of `dnascope` and `dnascope-pangenome`: - The two passes are merged as gVCFs. DNAscope does not clip reference blocks to `--interval`, so `bcftools view -T ^bed` kept blocks that run into the realigned regions and dropped blocks that run out of them. The new `scripts/gvcf_clip.py` keeps pass 1 outside the regions (`subset-calls`, `--exclude`) and pass 2 inside them (`subset-pass2-calls`), clipping reference blocks at the boundaries and keeping variant records by their start position, before the same `concat-calls`. - Annotation, the `--pop_vcf` transfer and DNAModelApply run on the gVCF, and DNAModelApply writes `<output>.g.vcf.gz`. - GVCFtyper, over `--bed` and with `--dbsnp`, genotypes it. GVCFtyper copies the hybrid FORMAT/LAD, LPL, SAD and SPL arrays from the gVCF without subsetting them to the alleles it keeps, so until it does the new `drop-allele-fields` job removes them with `bcftools annotate -x`; `<output>.g.vcf.gz` keeps them. `final-norm` then writes `--output_vcf`. - Under the hidden `--skip_model_apply`, `drop-allele-fields` writes `--output_vcf`, since that path skips `final-norm` without `--gvcf` too. `filter_norm` always excludes hom-ref records now and loses its `exclude_homref` argument. Without `--gvcf` the DAG is unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
`hybrid_stage2.bed`, the `all_bed` output of HybridStage2, covers the whole genome, since the realigned short reads map everywhere, so with `--bed` the second DNAscope pass and the merge added records outside `--bed`. A new `intersect-stage2-bed` job writes its intersection with `--bed`, sorted and merged (`cmd_bedtools_intersect`), and the second pass and the merge steps use it, in both VCF and gVCF mode. Stage 3 and `<output>_sr_realigned.cram` keep the whole stage-2 BED. An empty intersection needs no special case, since the driver treats an empty interval file as zero-length. Without `--bed` the DAG is unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
… ploidies when the sex is unknown `combine_cnv.py` 2.1.2 reads raw CNVscope output with and without HMM states the same way, and it stops when it cannot tell the sample's sex. CNVscope run without `--sex` writes no `##SampleSex`, so a sample whose sex could not be estimated would now fail at the combine step. When the sex is unknown, CNVscope and `combine_cnv.py` now get the same `--blocks` BED file in place of `--sex/--par`: chrX at two copies and chrY at one copy outside the PAR intervals. A warning reports that sex estimation failed and that CNV calling proceeds with assumed sex chromosome ploidies. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
`bcftools merge` sets QUAL to the highest QUAL of the merged records and, under its default `--filter-logic +`, combines their FILTERs, so the transfer's output took QUAL and FILTER from the population VCF wherever its QUAL was higher or it carried a filter, including at sites where the population allele is then trimmed away. Each shard now merges two `--no-index` process substitutions, each a `bcftools annotate` restricted to the shard BED with `--regions-overlap pos` (`-R` cannot be combined with `--no-index`): `-x FILTER` on the input VCF and `-x QUAL,FILTER` on the population VCF. The merge ignores a missing QUAL, so the output keeps the input VCF's QUAL, and FILTER is empty. Contigs the population VCF lacks skip the merge; their subset (`cmd_bcftools_subset_remove`, replacing `cmd_bcftools_view_regions`) drops FILTER too. Population IDs still carry over. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
`validate_bundle` read only the `platform` key of `bundle_info.json`, but `SentieonUltima1.1.bundle` names its platform `shortReadPlatform`, as the hybrid and pangenome bundles do, so the DNAscope call got no `UltimaReadFilter` with that bundle. The pipeline now reads `platform`, then `shortReadPlatform`; `platform` wins when a bundle has both. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
`gvcf_clip.py` keeps or drops variant records by their start position. A dropped record can span positions on the kept side, such as a pass-1 deletion that starts inside a realigned region and ends past it, and when no kept record covered those positions the spliced gVCF had no record for them: 54 gaps of 1 to 118 bp at pass-2 region ends in a whole-genome `dnascope-hybrid --gvcf` run. The script now covers such positions with a reference block that claims no genotype confidence (`0/0`, GQ 0, flat PL, DP and MIN_DP from the dropped record), written only with the FORMAT fields the header declares. A fill block waits for the records after it, which split it where they overlap, so the splice still adds no overlapping blocks, and it is cut at shard boundaries, since the records past a boundary belong to the next shard. The include side, for pass 2, is treated the same way. An output of `-` writes uncompressed VCF to stdout. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
With `--gvcf`, `hybrid_select.py` read the gVCF-mode pass-1 output. DNAscope's gVCF mode emits every candidate site, including `0/0` records that the script's long-read/short-read rules pass, and leaves `STR` off many sites, so a whole-genome run selected 57% more sequence than without `--gvcf` (2.28 Gb vs 1.46 Gb) and realigned and recalled different regions. Reading 564 M gVCF records also took the script 1,721 s rather than 76 s. The pass-1 driver run now carries a second, VCF-mode DNAscope algo with otherwise identical arguments under `--gvcf`, writing `initial.vcf.gz` next to `initial.g.vcf.gz`, and `hybrid-select` reads `initial.vcf.gz` in both modes. The selected and realigned regions are then those of a run without `--gvcf`; the gVCF still comes from `initial.g.vcf.gz`. The extra algo costs about the VCF-mode pass-1 CPU again. `rm-tmp5` removes `initial.vcf.gz`. Without `--gvcf` the DAG is unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
In a whole-genome `dnascope-hybrid --gvcf` run, `gvcf_clip.py` on pass 1 took 8,552 s with about 3 of 32 cores busy, and `hybrid_anno.py` 5,590 s. Both shard their work with vcflib, whose parent process then re-reads, compresses and indexes every output record of the 564 M-record gVCF on its own. Under `--gvcf`, both scripts now write uncompressed VCF to stdout, which `sentieon util vcfconvert` compresses and indexes, as `gvcf_combine.py` already does: `cmd_pyexec_gvcf_clip` always streams, and `cmd_pyexec_hybrid_anno` takes `stream`. `concat-calls` also gets `--threads` under `--gvcf`. Without `--gvcf` the commands are unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The shard merge reads its two inputs with `bcftools merge --no-index`,
which refuses inputs whose headers list different contigs, or the same
contigs in a different order:
[E::bcf_sr_add_hreader] Different number of sequences in the header,
refusing to stream multiple unindexed files
The input VCF lists the reference's contigs and a population VCF
generally lists others, so the transfer failed on such references in
every pipeline that takes `--pop_vcf`. Restricting each input to the
shard's region does not help, since the check is on the headers.
Each input stream now passes through `awk` (`KEEP_CONTIG_AWK`), keeping
only the `##contig` line of the shard's contig, so the two headers
agree; `cmd_bcftools_merge_trim` takes the shard's `contig`. A merged
shard VCF lists that contig alone, and the concat merges the shard
headers in shard order, so the transfer output lists the contigs in
reference order. Contigs that only the population VCF lists no longer
appear in the output header.
New run tests cover a population VCF with other contigs in another
order, and a three-shard transfer whose concatenated header comes back
in reference order.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The README pulled `sentieon/sentieon_cli`, but the images are published as `sentieon/sentieon-cli`, the name `docker.yml` pushes. The Dockerfile's `org.opencontainers.image.base.name` label now matches its `debian:13.7-slim` base image. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Update the CNV scripts to combine_sv_cnv 2.3.0 (combine_sv_cnv 6428c81). CNVModelApply in the shared `CNVscopeStage` now runs with `--emit_conf 0`, so the CNV VCF also holds the segments the model rejects, with FILTER `MLrejected`; the PASS records do not change. This covers `dnascope-hybrid`, `sentieon-pangenome` and `hybrid-pangenome`. `combine_cnv.py` reads the rejected segments from `--cnv`, so the pangenome combine step no longer passes `--raw`. `--emit_conf` is a hidden CNVModelApply option in 202503.04. The combine rules now come from the model bundle. When the bundle has a `combine_cnv.json` member, the combine step gets `--preset <bundle>/combine_cnv.json`; without it, `combine_cnv.py` uses its defaults, the paired-end values. The rule that gave Ultima short reads the SE preset is removed, and hybrid-pangenome no longer keeps the bundle's `shortReadPlatform`. An Ultima bundle needs the member to keep the single-end rules. `combine_cnv.py` 2.3.0 reads `--preset` from a JSON file or a bundle member and makes `--raw` optional for `--emit_conf 0` input; `indel2cnv.py` changes only its version. Tests cover the preset with and without the member on each platform, `--emit_conf 0` in both CNV stages, and loading a bundle preset with the packaged script. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
DNAscope already splits bwa into one job per NUMA node (halved until each piece has at most 64 CPUs) on large hosts. dnascope-pangenome ran bwa and minimap2 as single jobs across every node. With FASTQ input it now uses the same rule: - BwaExtractStage builds one bwa-extract job per lane: `fqidx extract -F j/N` feeds `bwa mem -p`, and bwa and `pgutil extract` are pinned to the lane's CPUs with taskset. Each lane writes sample-bwa-<j>.bam and sample-bwa-<j>.fq.gz. - One mm2-lift job per lane reads that lane's extracted reads, pinned to the same CPUs, and waits only for its own bwa lane. - Both dedup stages read every lane's output. Job names are unchanged, so job metrics still group the lanes as one step. - A hidden --no_split_alignment keeps the single jobs. BAM/CRAM input and hosts with 32 or fewer cores are unchanged. `fqidx extract` rejects -K below 10,000,000, while the pangenome bwa uses -K 2,000,000, so split runs use -K 10,000,000 for both. Unsplit runs keep 2,000,000. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
indel2cnv sorts the candidate variants by length, longest first, so the expensive ones are dispatched early. It then grouped them into batches of 5 and called starmap with its default chunksize, which is about len(batches) / (4 * workers). The first chunk therefore held the longest variants of the whole run and went to a single worker, which finished long after the others went idle. Each task is now one variant, with chunksize=1. Results are collected in the same order as before, so the output does not change. On HG002 35x (HPRC v2.0), indel2cnv took 420 -> 68 s on 64 vCPUs and 156 -> 49 s on 192 vCPUs, and the CNV VCF body was identical in both cases. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
vg-haplotypes only needs the k-mer counts, so it is meant to overlap bwa on hosts with enough memory. It did not: the job reserved every core, so the scheduler held it until bwa released its threads, whatever its dependencies said. When the host has room for vg haplotypes (100 GB) plus each bwa lane (32 GB), the job now reserves no threads, as kmc already does, and runs next to bwa. Otherwise it depends on every bwa lane, as before. This replaces the fixed 70 GB check. On HG002 35x with 4 lanes on a 192-vCPU, 384 GB host, the pipeline took 30:10 -> 26:24, with indel2cnv's change included. The output VCF body was unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Review follow-up: - Job names are unique again: split runs name their lanes bwa-extract-<j> and mm2-lift-<j>, while unsplit runs keep bwa-extract and mm2-lift. - dnascope-pangenome aligns with -K 10,000,000 whether or not it splits, so its results no longer depend on how many NUMA lanes the host has. `fqidx extract` needs at least that, which is why split runs already used it. BwaExtractStage takes the -K as an option, so hybrid-pangenome keeps its current 2,000,000. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
dnascope-pangenome: NUMA alignment lanes, faster indel2cnv, vg haplotypes alongside bwa
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
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.
No description provided.