Skip to content

Severus fix #223

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
robert-a-forsyth merged 2 commits into from
Oct 14, 2026
Merged

Severus fix #223

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c29be20
Severus fix: run Severus on the haplotagged BAMs
robert-a-forsyth Oct 14, 2026
b6ed884
Add PR number to CHANGELOG
robert-a-forsyth Oct 14, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0

### `Fixed`

- [#223](https://github.com/IntGenomicsLab/lrsomatic/pull/223) - Fixed the Severus input channel (@robert-a-forsyth).
- [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - `PREPARE_AA_DATA_REPO` no longer reads a `versions` output `UNTAR` does not have, which crashed every default GRCh38 run with CoRAL; overriding `--aa_data_repo_url` drops the genome's default MD5 (@robert-a-forsyth).
- [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - `-profile test_full` starts again: it skips CoRAL and signatures and carries the chr19 SAVANA workarounds from `test` (@robert-a-forsyth).
- [#205](https://github.com/IntGenomicsLab/lrsomatic/pull/205) - `CLAIRSTO_VERDICT_TAG` no longer crashes when ASCAT finds no solution and writes an empty segments table (@robert-a-forsyth).
Expand Down
99 changes: 37 additions & 62 deletions workflows/lrsomatic.nf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -867,31 +867,47 @@ workflow LRSOMATIC {
)
}

// SEVERUS input: tumor-only and paired samples with their phased germline VCF; [] normal BAM/BAI = tumor-only mode
branched_minimap.tumor_only
.map{ meta, bam, bai ->
def new_meta = meta.subMap('id',
'paired_data',
'platform',
'sex',
'fiber',
'clair3_model',
'clairS_model',
'clairSTO_model',
'kinetics')
return[new_meta, bam, bai]
// SEVERUS input: tumor-only and paired samples with their phased germline VCF; [] normal BAM/BAI = tumor-only mode.
// Severus, Wakhan and SAVANA read the HP (haplotype) tag per read, so build it from
// PHASING_HAPLOTYPING's haplotagged BAMs rather than the unphased minimap2 BAMs.
def sv_meta_keys = ['id', 'paired_data', 'platform', 'sex', 'fiber',
'clair3_model', 'clairS_model', 'clairSTO_model', 'kinetics']

PHASING_HAPLOTYPING.out.tumor_normal_hapbams_ch

@claude claude Bot Oct 14, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Snapshot is not name-only for Severus. The PR description says "snapshots list Severus outputs by name only, so no snapshot changes expected". That's wrong. tests/.nftignore only ignores Severus's read_ids.csv, severus.log and the *.vcf.gz/.tbi files. breakpoints_double.csv, read_qual.txt and {all_SVs,somatic_SVs}/breakpoint_clusters*.tsv are md5-checked in tests/default.nf.test.snap (lines 612–634), and default.nf.test is tagged small.

Now that Severus gets HP-tagged reads, the per-haplotype support in breakpoints_double.csv should change. One of two things will happen:

  • the md5s drift and the small CI test fails until the snapshot is regenerated, or
  • they don't drift, which suggests the test data doesn't show the fix working.

Please regenerate tests/default.nf.test.snap and explain the Severus md5 changes in the PR. Also regenerate the extended snapshots (clair_only, consensus, deep_only, union, chm13) that include Severus outputs.

.branch { meta, _bam, _bai ->
tumor_only: !meta.paired_data
paired_tumor: meta.paired_data && meta.type == 'tumor'
paired_normal: meta.paired_data && meta.type == 'normal'
}
.map{meta, tumor_bam, tumor_bai->
.set { branched_hapbams }

branched_hapbams.tumor_only
.map { meta, bam, bai ->
def normal_bam = []
def normal_bai = []
return [meta, tumor_bam, tumor_bai, normal_bam, normal_bai]
return [meta.subMap(sv_meta_keys), bam, bai, normal_bam, normal_bai]
}
// Mix with paired T/N input (which already has normal BAM/BAI from somatic_smallvar_input)
.mix(somatic_smallvar_input)
.set { tumoronly_hapbams }

branched_hapbams.paired_tumor
.map { meta, bam, bai -> return [meta.subMap(sv_meta_keys), bam, bai] }
.set { paired_tumor_hapbams }

branched_hapbams.paired_normal
.map { meta, bam, bai -> return [meta.subMap(sv_meta_keys), bam, bai] }
.set { paired_normal_hapbams }

paired_tumor_hapbams
.join(paired_normal_hapbams)
.set { paired_hapbams }
// paired_hapbams: [meta, tumor_bam, tumor_bai, normal_bam, normal_bai]

tumoronly_hapbams
.mix(paired_hapbams)
// Attach phased germline VCF (used by SEVERUS for phased SV calling)
.join(PHASING_HAPLOTYPING.out.phased_germline_vcf)
.set{severus_input}
Comment on lines 907 to 909

@claude claude Bot Oct 14, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Missing test coverage for the fix, and Wakhan's input change is untested.

  1. No test checks the behaviour this PR fixes. That bug was severus.log reporting H1 / H2 coverage 0.0 / 0.0 and every hVAF being x,0,0. default.nf.test only checks that severus_somatic.vcf.gz exists, and severus.log is in .nftignore. If someone switches severus_input back to the unphased BAMs, CI would still pass. Please add an assertion to the small default test, e.g. that sample1/variants/severus/severus.log reports non-zero H1/H2 coverage, or that some Severus record has a non-zero haplotype VAF.
  2. Wakhan now gets the haplotagged BAMs too, but conf/test.config sets skip_wakhan = true. The only small test that enables Wakhan (sv_annotation_compose.nf.test) is a preview run that only checks process names. So no CI test runs Wakhan on the new input. Per REVIEW.md, please run the Slurm test array (or a real --skip_wakhan false run) and report the result before merging.

// severus_input: [meta, tumor_bam, tumor_bai, normal_bam, normal_bai, phased_germline_vcf, phased_germline_tbi]
// severus_input: [meta, tumor_hapbam, tumor_bai, normal_hapbam, normal_bai, phased_germline_vcf, phased_germline_tbi]
// normal_bam/bai are empty lists [] for tumor-only samples

//
Expand Down Expand Up @@ -1046,55 +1062,14 @@ workflow LRSOMATIC {
savana_allele_counts = channel.empty()

if (!params.skip_savana) {
// SAVANA reads the HP (haplotype) tag per read and its README recommends phased BAMs,
// so build its input from PHASING_HAPLOTYPING's haplotagged BAMs rather than the
// unphased ones severus_input carries.
def savana_meta_keys = ['id', 'paired_data', 'platform', 'sex', 'fiber',
'clair3_model', 'clairS_model', 'clairSTO_model', 'kinetics']

PHASING_HAPLOTYPING.out.tumor_normal_hapbams_ch
.branch { meta, _bam, _bai ->
tumor_only: !meta.paired_data
paired_tumor: meta.paired_data && meta.type == 'tumor'
paired_normal: meta.paired_data && meta.type == 'normal'
}
.set { branched_savana_hapbams }

branched_savana_hapbams.tumor_only
.map { meta, bam, bai ->
def normal_bam = []
def normal_bai = []
return [meta.subMap(savana_meta_keys), bam, bai, normal_bam, normal_bai]
}
.set { savana_tumoronly_hapbams }

branched_savana_hapbams.paired_tumor
.map { meta, bam, bai -> return [meta.subMap(savana_meta_keys), bam, bai] }
.set { savana_paired_tumor_hapbams }

branched_savana_hapbams.paired_normal
.map { meta, bam, bai -> return [meta.subMap(savana_meta_keys), bam, bai] }
.set { savana_paired_normal_hapbams }

savana_paired_tumor_hapbams
.join(savana_paired_normal_hapbams)
.set { savana_paired_hapbams }
// savana_paired_hapbams: [meta, tumor_bam, tumor_bai, normal_bam, normal_bai]

savana_tumoronly_hapbams
.mix(savana_paired_hapbams)
.join(PHASING_HAPLOTYPING.out.phased_germline_vcf)
.set { savana_input }
// savana_input: [meta, tumor_bam, tumor_bai, normal_bam, normal_bai, phased_germline_vcf, phased_germline_tbi]
// normal_bam/bai are [] for tumor-only samples

savana_input
// SAVANA shares severus_input: haplotagged BAMs with the phased germline VCF
severus_input
.branch { meta, _tumor_bam, _tumor_bai, normal_bam, _normal_bai, _phased_vcf, _phased_tbi ->
tumor_only: !normal_bam
paired: normal_bam
}
.set { branched_savana_input }
// branched_savana_input.tumor_only / .paired: same 7-tuple shape as savana_input
// branched_savana_input.tumor_only / .paired: same 7-tuple shape as severus_input

branched_savana_input.tumor_only
.map { meta, tumor_bam, tumor_bai, _normal_bam, _normal_bai, phased_vcf, phased_tbi ->
Expand Down
Loading