From 3a13913a4654b4dc7564bafe55399374ab5507dc Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 17:14:35 +0100 Subject: [PATCH 01/20] fix(report): read a header-only Y tree path table as an empty path y_tree_place.py writes the header alone when no branch is derived; the report marked the whole Y tree section unreadable. Require the Branch column in the header, not a row. Bead: pgp-5bja.25 --- bin/collect_summary.py | 11 ++++++++--- tests/test_collect_summary.py | 12 ++++++++++-- 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/bin/collect_summary.py b/bin/collect_summary.py index 0caf4f4b..6dba8cea 100755 --- a/bin/collect_summary.py +++ b/bin/collect_summary.py @@ -2638,9 +2638,14 @@ def sec_y_tree(d, s): "explore": haplogroup_explore("y", hg)} lvl = first_existing(d, [f"y_haplogroup/{s}_y_tree_path.tsv"]) if lvl: - prows = read_tsv(lvl) - if not prows or "Branch" not in prows[0]: - raise Unreadable(f"{os.path.basename(lvl)} has no Branch column or no rows") + # bin/y_tree_place.py writes the header alone when no branch is derived: + # an empty path, not an unreadable file. + with open_text(lvl) as f: + rd = csv.DictReader(f, delimiter="\t") + prows = list(rd) + header = rd.fieldnames or [] + if "Branch" not in header: + raise Unreadable(f"{os.path.basename(lvl)} has no Branch column") out["path"] = [{k: (row.get(k) or "") for k in ("Level", "Branch", "SNPs", "On_tables", "Derived", "Ancestral", "No_call", "Ancestral_SNPs")} for row in prows] diff --git a/tests/test_collect_summary.py b/tests/test_collect_summary.py index 501aa928..fb800a62 100644 --- a/tests/test_collect_summary.py +++ b/tests/test_collect_summary.py @@ -2416,11 +2416,19 @@ def main(): "Path_no_call\tAncestral_in_terminal\tSNPs_called\tMin_depth\tTree\tPositions\t" "Outvoted_ancestral\nS\tNA\t0\t0\t0\t0\t.\t12\t10\tYFull YTree 14.06.0\t" "Yleaf 4.1.6 tables\t.\n") + # bin/y_tree_place.py writes the path table's header alone when no branch is derived. + put(f"{h}/y_haplogroup/S_y_tree_path.tsv", "Level\tBranch\tSNPs\tOn_tables\tDerived\tAncestral\t" + "No_call\tAncestral_SNPs\n") summ_na, txt_na, _ = render(h) check("Y tree: Hg NA reads as no derived call, and nothing ancestral is printed", - summ_na["sections"]["y_tree"]["data"].get("haplogroup") == "no branch has a derived call" + summ_na["sections"]["y_tree"]["state"] == "ok" + and summ_na["sections"]["y_tree"]["data"].get("haplogroup") == "no branch has a derived call" + and summ_na["sections"]["y_tree"]["data"].get("path") == [] and "Ancestral inside the terminal branch" not in txt_na, - summ_na["sections"]["y_tree"]["data"]) + summ_na["sections"]["y_tree"]) + put(f"{h}/y_haplogroup/S_y_tree_path.tsv", "Level\tSNPs\n") + check("Y tree: a path table with no Branch column is unreadable", + render(h)[0]["sections"]["y_tree"]["state"] == "unreadable") os.remove(f"{h}/y_haplogroup/S_y_tree.txt") os.remove(f"{h}/y_haplogroup/S_y_tree_path.tsv") From 5c2776d13b36d9b2f2b44938dbc5e4b0c4b85173 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 17:16:24 +0100 Subject: [PATCH 02/20] fix(report): read no copy from a phased GT with an uncalled allele phase_haps gave {1} for '.|1', so a deletion '.|1' and a hit '0|1' in one phase set were called in cis and dropped as no second hit. A phased GT with a missing allele now gives no phase. Every e2e case that runs bin/cnvpytor_to_vcf.py now lists bin/family_report.py, which it imports. Bead: pgp-5bja.24 --- bin/collect_summary.py | 9 +++++++-- scripts/ci/e2e-case-paths.tsv | 18 +++++++++--------- tests/test_collect_summary.py | 10 ++++++++++ 3 files changed, 26 insertions(+), 11 deletions(-) diff --git a/bin/collect_summary.py b/bin/collect_summary.py index 6dba8cea..6df0cd61 100755 --- a/bin/collect_summary.py +++ b/bin/collect_summary.py @@ -3901,10 +3901,15 @@ def over(chrom, a, b): def phase_haps(gt): """The copies (haplotype indexes) a phased GT puts the alt allele on, or - None when the GT is not phased: '0|1' -> {1}, '1|0' -> {0}, '0/1' -> None.""" + None when the GT is not phased or a copy is not called: '0|1' -> {1}, + '1|0' -> {0}, '0/1' -> None, '.|1' -> None (the missing copy could carry + the other allele, so no cis or trans can be read from it).""" if not gt or "|" not in gt: return None - return frozenset(i for i, a in enumerate(gt.split("|")) if a not in ("0", ".")) + alleles = gt.split("|") + if "." in alleles: + return None + return frozenset(i for i, a in enumerate(alleles) if a != "0") def pair_phase(dl, v): diff --git a/scripts/ci/e2e-case-paths.tsv b/scripts/ci/e2e-case-paths.tsv index 5dc19dc9..7021aaec 100644 --- a/scripts/ci/e2e-case-paths.tsv +++ b/scripts/ci/e2e-case-paths.tsv @@ -48,9 +48,9 @@ alignment-markdup-gvcf-sex-3-homref assets/pgs_scores.tsv bin/collect_summary.py alignment-markdup-gvcf-sex-4-bwamem2 scripts/02a-alignment-bwamem2.sh scripts/ci/fake-docker/lib.sh ancestry-fine-grained-1 assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/mito_confirmed.tsv assets/pgs_scores.tsv assets/stub/sample.vcf.gz bin/ancestry_fine_report.py bin/collect_summary.py bin/render_report.py bin/results_index.py modules/local/ancestry_fine/* modules/local/html_report/* modules/local/prs/* modules/local/vcf_precheck/* scripts/24-html-report.sh scripts/25-prs.sh scripts/26-ancestry.sh scripts/26b-ancestry-fine.sh scripts/chip-to-vcf.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh scripts/validate-setup.sh ancient-similarity-1 assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/mito_confirmed.tsv assets/stub/sample.vcf.gz bin/ancient_report.py bin/collect_summary.py bin/render_report.py bin/results_index.py modules/local/ancient/* modules/local/html_report/* modules/local/vcf_precheck/* scripts/24-html-report.sh scripts/44-ancient.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh scripts/validate-setup.sh -annotsv-delly-depth-1-chain assets/par_grch38.bed assets/sv_rule_regions.tsv bin/annotsv_postfilter.py bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/sv_svlen_fix.py scripts/05-annotsv.sh scripts/15-duphold.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/ci/fake-docker/lib.sh tests/fixtures/annotsv/class45.sh -annotsv-report-section-1-chain assets/par_grch38.bed assets/sv_rule_regions.tsv bin/annotsv_postfilter.py bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/sv_svlen_fix.py scripts/05-annotsv.sh scripts/15-duphold.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/22-survivor-merge.sh scripts/24-html-report.sh scripts/ci/fake-docker/lib.sh scripts/generate-report.sh tests/fixtures/annotsv/class45.sh -annotsv-sv-inputs-1-manta assets/par_grch38.bed assets/sv_rule_regions.tsv bin/annotsv_postfilter.py bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/manta_inv_filter.py bin/sv_svlen_fix.py scripts/04-manta.sh scripts/05-annotsv.sh scripts/15-duphold.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/ci/fake-docker/lib.sh +annotsv-delly-depth-1-chain assets/par_grch38.bed assets/sv_rule_regions.tsv bin/annotsv_postfilter.py bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/family_report.py bin/sv_svlen_fix.py scripts/05-annotsv.sh scripts/15-duphold.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/ci/fake-docker/lib.sh tests/fixtures/annotsv/class45.sh +annotsv-report-section-1-chain assets/par_grch38.bed assets/sv_rule_regions.tsv bin/annotsv_postfilter.py bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/family_report.py bin/sv_svlen_fix.py scripts/05-annotsv.sh scripts/15-duphold.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/22-survivor-merge.sh scripts/24-html-report.sh scripts/ci/fake-docker/lib.sh scripts/generate-report.sh tests/fixtures/annotsv/class45.sh +annotsv-sv-inputs-1-manta assets/par_grch38.bed assets/sv_rule_regions.tsv bin/annotsv_postfilter.py bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/family_report.py bin/manta_inv_filter.py bin/sv_svlen_fix.py scripts/04-manta.sh scripts/05-annotsv.sh scripts/15-duphold.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/ci/fake-docker/lib.sh archaic-ancestry-1 assets/pgs_scores.tsv bin/archaic_report.py bin/collect_summary.py modules/local/archaic/* modules/local/prs/* modules/local/vcf_precheck/* scripts/25-prs.sh scripts/26-ancestry.sh scripts/43-archaic.sh scripts/ci/fake-docker/lib.sh scripts/generate-report.sh scripts/setup.sh scripts/validate-setup.sh bash-step-08-hla - bash-step-09-str assets/str_extra_loci.json assets/stub/db/catalog.json bin/rfc1_motif.py bin/strchive_convert.py scripts/09-expansion-hunter.sh scripts/09b-stranger.sh scripts/ci/fake-docker/lib.sh @@ -100,18 +100,18 @@ hla-second-caller-2-nextflow assets/acmg_sf.tsv assets/callability_genes.tsv ass ibd-segments-1 assets/families.example.tsv assets/par_grch38.bed assets/pedigree.example.tsv assets/stub/pedigree.tsv bin/family_report.py bin/family_segments.py bin/ibd_report.py scripts/38-family.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh live-server-1-weblog assets/live/index.html assets/live/live.css assets/live/login.html assets/viewer/index.html bin/owngenome_live.py bin/render_report.py bin/run_state.py modules/local/roh/* modules/local/vcf_precheck/* scripts/11-roh-analysis.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh lpa-1 assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/mito_confirmed.tsv assets/par_grch38.bed assets/stub/sample.vcf.gz bin/cn_loci.tsv bin/collect_summary.py bin/copy_number.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/render_report.py bin/results_index.py modules/local/deepvariant/* modules/local/html_report/* modules/local/lpa/* modules/local/vcf_precheck/* scripts/24-html-report.sh scripts/46-lpa.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh -mei-1 assets/mei_sites.tsv assets/par_grch38.bed assets/sv_rule_regions.tsv bin/annotsv_postfilter.py bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/mei_report.py bin/sv_svlen_fix.py modules/local/align_minimap2/* modules/local/annotsv/* modules/local/deepvariant/* modules/local/fastp/* modules/local/mei/* modules/local/vcf_precheck/* scripts/01b-fastp-qc.sh scripts/05-annotsv.sh scripts/15-duphold.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/50-mei.sh scripts/ci/check-images.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh scripts/validate-setup.sh tests/fixtures/annotsv/class45.sh +mei-1 assets/mei_sites.tsv assets/par_grch38.bed assets/sv_rule_regions.tsv bin/annotsv_postfilter.py bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/family_report.py bin/mei_report.py bin/sv_svlen_fix.py modules/local/align_minimap2/* modules/local/annotsv/* modules/local/deepvariant/* modules/local/fastp/* modules/local/mei/* modules/local/vcf_precheck/* scripts/01b-fastp-qc.sh scripts/05-annotsv.sh scripts/15-duphold.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/50-mei.sh scripts/ci/check-images.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh scripts/validate-setup.sh tests/fixtures/annotsv/class45.sh mito-extended-1-deletions assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/construct_variants.tsv assets/deep_intronic_sites.tsv assets/mito_confirmed.tsv assets/stub/db/clinvar.vcf.gz assets/stub/db/clinvar.vcf.gz.tbi assets/stub/sample.bam assets/stub/sample.bam.bai assets/stub/sample.vcf.gz bin/cdna_scan.py bin/clinvar_hits.awk bin/clinvar_join_adjacent.awk bin/clip_check.py bin/collect_summary.py bin/render_report.py bin/results_index.py modules/local/clinvar_screen/* modules/local/html_report/* modules/local/mito_variants/* modules/local/vcf_precheck/* scripts/06-clinvar-screen.sh scripts/16b-mosdepth.sh scripts/20-mtoolbox.sh scripts/24-html-report.sh scripts/ci/build-mito-confirmed.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh mocha-1 assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/mito_confirmed.tsv assets/stub/sample.vcf.gz bin/collect_summary.py bin/mocha_report.py bin/render_report.py bin/results_index.py modules/local/html_report/* modules/local/mocha/* modules/local/vcf_precheck/* scripts/24-html-report.sh scripts/51-mocha.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh scripts/validate-setup.sh mosaic-loy-lox-1 assets/callability_genes.tsv bin/collect_summary.py bin/mosaic_sex_loss.py scripts/16b-mosdepth.sh scripts/33-sample-qc.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh nextflow-from-fastq-1-bash assets/callability_genes.tsv assets/construct_variants.tsv assets/deep_intronic_sites.tsv assets/pgs_scores.tsv assets/pgx_genes.tsv assets/stub/db/hla.dat assets/stub/db/kir.dat bin/cdna_scan.py bin/clinvar_hits.awk bin/clinvar_join_adjacent.awk bin/clip_check.py bin/cn_loci.tsv bin/collect_summary.py bin/copy_number.py bin/cyp2d6_depth_check.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/pgx_mito_merge.py bin/pgx_outside_calls.py bin/pgx_refblock_fill.py scripts/01b-fastp-qc.sh scripts/06-clinvar-screen.sh scripts/07-pharmacogenomics.sh scripts/08-hla-typing.sh scripts/11-roh-analysis.sh scripts/16b-mosdepth.sh scripts/20-mtoolbox.sh scripts/25-prs.sh scripts/32-pypgx.sh scripts/36-pgx-consensus.sh scripts/ci/fake-docker/lib.sh scripts/ci/parity-diff.sh scripts/validate-setup.sh nextflow-from-fastq-2-nextflow assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/construct_variants.tsv assets/deep_intronic_sites.tsv assets/live/index.html assets/mito_confirmed.tsv assets/par_grch38.bed assets/pgs_scores.tsv assets/pgx_contraindications.tsv assets/pgx_genes.tsv assets/pgx_inhibitors.tsv assets/stub/db/hla.dat assets/stub/db/kir.dat assets/stub/sample.bam assets/stub/sample.bam.bai assets/stub/sample.vcf.gz assets/viewer/index.html bin/cdna_scan.py bin/clinvar_hits.awk bin/clinvar_join_adjacent.awk bin/clip_check.py bin/cn_loci.tsv bin/collect_summary.py bin/copy_number.py bin/cyp2d6_depth_check.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/pgx_mito_merge.py bin/pgx_outside_calls.py bin/pgx_parse.py bin/pgx_refblock_fill.py bin/render_report.py bin/results_index.py modules/local/align_minimap2/* modules/local/clinvar_screen/* modules/local/cpic_lookup/* modules/local/cyrius/* modules/local/deepvariant/* modules/local/fastp/* modules/local/hla_typing/* modules/local/html_report/* modules/local/indexcov/* modules/local/mito_haplogroup/* modules/local/mito_variants/* modules/local/mosdepth/* modules/local/multiqc/* modules/local/pharmcat/* modules/local/prs/* modules/local/roh/* modules/local/t1k_build/* modules/local/telomere_hunter/* modules/local/vcf_precheck/* modules/local/vcfanno/* scripts/01b-fastp-qc.sh scripts/06-clinvar-screen.sh scripts/07-pharmacogenomics.sh scripts/08-hla-typing.sh scripts/10-telomere-hunter.sh scripts/11-roh-analysis.sh scripts/12-mito-haplogroup.sh scripts/16-indexcov.sh scripts/16b-mosdepth.sh scripts/20-mtoolbox.sh scripts/21-cyrius.sh scripts/24-html-report.sh scripts/25-prs.sh scripts/26-ancestry.sh scripts/27-cpic-lookup.sh scripts/28-multiqc.sh scripts/30-vcfanno.sh scripts/32-pypgx.sh scripts/36-pgx-consensus.sh scripts/ci/check-images.sh scripts/ci/fake-docker/lib.sh scripts/cyrius-constraints.txt scripts/validate-setup.sh nextflow-from-fastq-3-sexcheck assets/callability_genes.tsv assets/live/index.html assets/viewer/index.html bin/collect_summary.py modules/local/indexcov/* modules/local/mosdepth/* modules/local/vcf_precheck/* scripts/16-indexcov.sh scripts/16b-mosdepth.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh -nextflow-hardening-1-manta assets/par_grch38.bed assets/sv_rule_regions.tsv bin/annotsv_postfilter.py bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/manta_inv_filter.py bin/sv_svlen_fix.py scripts/04-manta.sh scripts/05-annotsv.sh scripts/15-duphold.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/ci/fake-docker/lib.sh -nextflow-hardening-2-nf assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/construct_variants.tsv assets/deep_intronic_sites.tsv assets/live/index.html assets/mito_confirmed.tsv assets/par_grch38.bed assets/pgx_genes.tsv assets/stub/sample.bam assets/stub/sample.bam.bai assets/stub/sample.vcf.gz assets/sv_rule_regions.tsv assets/viewer/index.html bin/annotsv_postfilter.py bin/cdna_scan.py bin/clinvar_hits.awk bin/clinvar_join_adjacent.awk bin/clip_check.py bin/cn_loci.tsv bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/copy_number.py bin/cyp2d6_depth_check.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/manta_inv_filter.py bin/pgx_mito_merge.py bin/pgx_outside_calls.py bin/pgx_refblock_fill.py bin/render_report.py bin/results_index.py bin/sv_svlen_fix.py bin/y_tree_place.py bin/yleaf_run.py modules/local/clinvar_screen/* modules/local/delly/* modules/local/html_report/* modules/local/indexcov/* modules/local/manta/* modules/local/mito_haplogroup/* modules/local/mosdepth/* modules/local/pgx_consensus/* modules/local/pharmcat/* modules/local/pypgx/* modules/local/survivor_merge/* modules/local/telomere_hunter/* modules/local/vcf_precheck/* modules/local/vcfanno/* modules/local/y_haplogroup/* scripts/04-manta.sh scripts/05-annotsv.sh scripts/06-clinvar-screen.sh scripts/07-pharmacogenomics.sh scripts/10-telomere-hunter.sh scripts/12-mito-haplogroup.sh scripts/15-duphold.sh scripts/16-indexcov.sh scripts/16b-mosdepth.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/20-mtoolbox.sh scripts/21-cyrius.sh scripts/22-survivor-merge.sh scripts/24-html-report.sh scripts/30-vcfanno.sh scripts/32-pypgx.sh scripts/36-pgx-consensus.sh scripts/37-y-haplogroup.sh scripts/ci/fake-docker/lib.sh scripts/ci/parity-diff.sh scripts/cyrius-constraints.txt scripts/setup.sh scripts/validate-setup.sh -nextflow-report-parity-1-sections assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/construct_variants.tsv assets/deep_intronic_sites.tsv assets/mito_confirmed.tsv assets/par_grch38.bed assets/pgx_genes.tsv assets/str_extra_loci.json assets/stub/db/catalog.json assets/stub/sample.bam assets/stub/sample.bam.bai assets/stub/sample.vcf.gz assets/sv_rule_regions.tsv bin/annotsv_postfilter.py bin/cdna_scan.py bin/clinvar_hits.awk bin/clinvar_join_adjacent.awk bin/clip_check.py bin/cn_loci.tsv bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/comphet_phase.py bin/constraint_join.awk bin/copy_number.py bin/csq_rank.awk bin/cyp2d6_depth_check.py bin/export_consumer.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/pgx_outside_calls.py bin/render_report.py bin/results_index.py bin/rfc1_motif.py bin/strchive_convert.py bin/sv_svlen_fix.py modules/local/annotsv/* modules/local/clinvar_screen/* modules/local/expansion_hunter/* modules/local/export_consumer/* modules/local/html_report/* modules/local/mito_variants/* modules/local/pypgx/* modules/local/slivar/* modules/local/stranger/* modules/local/vcf_precheck/* modules/local/vep/* scripts/05-annotsv.sh scripts/06-clinvar-screen.sh scripts/09-expansion-hunter.sh scripts/09b-stranger.sh scripts/13-vep-annotation.sh scripts/15-duphold.sh scripts/16b-mosdepth.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/20-mtoolbox.sh scripts/24-html-report.sh scripts/31-slivar.sh scripts/32-pypgx.sh scripts/36-pgx-consensus.sh scripts/40-export-consumer.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh tests/fixtures/consumer_export/sites.tsv.gz +nextflow-hardening-1-manta assets/par_grch38.bed assets/sv_rule_regions.tsv bin/annotsv_postfilter.py bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/family_report.py bin/manta_inv_filter.py bin/sv_svlen_fix.py scripts/04-manta.sh scripts/05-annotsv.sh scripts/15-duphold.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/ci/fake-docker/lib.sh +nextflow-hardening-2-nf assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/construct_variants.tsv assets/deep_intronic_sites.tsv assets/live/index.html assets/mito_confirmed.tsv assets/par_grch38.bed assets/pgx_genes.tsv assets/stub/sample.bam assets/stub/sample.bam.bai assets/stub/sample.vcf.gz assets/sv_rule_regions.tsv assets/viewer/index.html bin/annotsv_postfilter.py bin/cdna_scan.py bin/clinvar_hits.awk bin/clinvar_join_adjacent.awk bin/clip_check.py bin/cn_loci.tsv bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/copy_number.py bin/cyp2d6_depth_check.py bin/family_report.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/manta_inv_filter.py bin/pgx_mito_merge.py bin/pgx_outside_calls.py bin/pgx_refblock_fill.py bin/render_report.py bin/results_index.py bin/sv_svlen_fix.py bin/y_tree_place.py bin/yleaf_run.py modules/local/clinvar_screen/* modules/local/delly/* modules/local/html_report/* modules/local/indexcov/* modules/local/manta/* modules/local/mito_haplogroup/* modules/local/mosdepth/* modules/local/pgx_consensus/* modules/local/pharmcat/* modules/local/pypgx/* modules/local/survivor_merge/* modules/local/telomere_hunter/* modules/local/vcf_precheck/* modules/local/vcfanno/* modules/local/y_haplogroup/* scripts/04-manta.sh scripts/05-annotsv.sh scripts/06-clinvar-screen.sh scripts/07-pharmacogenomics.sh scripts/10-telomere-hunter.sh scripts/12-mito-haplogroup.sh scripts/15-duphold.sh scripts/16-indexcov.sh scripts/16b-mosdepth.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/20-mtoolbox.sh scripts/21-cyrius.sh scripts/22-survivor-merge.sh scripts/24-html-report.sh scripts/30-vcfanno.sh scripts/32-pypgx.sh scripts/36-pgx-consensus.sh scripts/37-y-haplogroup.sh scripts/ci/fake-docker/lib.sh scripts/ci/parity-diff.sh scripts/cyrius-constraints.txt scripts/setup.sh scripts/validate-setup.sh +nextflow-report-parity-1-sections assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/construct_variants.tsv assets/deep_intronic_sites.tsv assets/mito_confirmed.tsv assets/par_grch38.bed assets/pgx_genes.tsv assets/str_extra_loci.json assets/stub/db/catalog.json assets/stub/sample.bam assets/stub/sample.bam.bai assets/stub/sample.vcf.gz assets/sv_rule_regions.tsv bin/annotsv_postfilter.py bin/cdna_scan.py bin/clinvar_hits.awk bin/clinvar_join_adjacent.awk bin/clip_check.py bin/cn_loci.tsv bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/comphet_phase.py bin/constraint_join.awk bin/copy_number.py bin/csq_rank.awk bin/cyp2d6_depth_check.py bin/export_consumer.py bin/family_report.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/pgx_outside_calls.py bin/render_report.py bin/results_index.py bin/rfc1_motif.py bin/strchive_convert.py bin/sv_svlen_fix.py modules/local/annotsv/* modules/local/clinvar_screen/* modules/local/expansion_hunter/* modules/local/export_consumer/* modules/local/html_report/* modules/local/mito_variants/* modules/local/pypgx/* modules/local/slivar/* modules/local/stranger/* modules/local/vcf_precheck/* modules/local/vep/* scripts/05-annotsv.sh scripts/06-clinvar-screen.sh scripts/09-expansion-hunter.sh scripts/09b-stranger.sh scripts/13-vep-annotation.sh scripts/15-duphold.sh scripts/16b-mosdepth.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/20-mtoolbox.sh scripts/24-html-report.sh scripts/31-slivar.sh scripts/32-pypgx.sh scripts/36-pgx-consensus.sh scripts/40-export-consumer.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh tests/fixtures/consumer_export/sites.tsv.gz nf-wiring-1-report-tables assets/acmg_sf.tsv assets/test_samplesheet.csv assets/stub/* bin/clip_check.py bin/collect_summary.py bin/cpsr_sample_id bin/render_report.py bin/results_index.py conf/test_all.config conf/test.config modules/local/align_minimap2/* modules/local/clinical_filter/* modules/local/clinvar_screen/* modules/local/cpsr/* modules/local/html_report/* modules/local/slivar/* modules/local/vep/* scripts/24-html-report.sh tests/fixtures/consumer_export/sites.tsv.gz -nf-wiring-2-sv-str-tools bin/cnvpytor_to_vcf.py bin/manta_inv_filter.py bin/strchive_convert.py bin/sv_svlen_fix.py modules/local/cnvpytor/* modules/local/duphold/* modules/local/expansion_hunter/* modules/local/indexcov/* modules/local/manta/* modules/local/vcf_precheck/* +nf-wiring-2-sv-str-tools bin/cnvpytor_to_vcf.py bin/family_report.py bin/manta_inv_filter.py bin/strchive_convert.py bin/sv_svlen_fix.py modules/local/cnvpytor/* modules/local/duphold/* modules/local/expansion_hunter/* modules/local/indexcov/* modules/local/manta/* modules/local/vcf_precheck/* output-contract-1-hg002 assets/* bin/* modules/* scripts/* paralog-catalog-1-gba assets/paralog_loci.tsv assets/parascopy_models_v1.7.txt bin/cn_loci.tsv bin/copy_number.py bin/hba_copy_number.py bin/lpa_kiv2.py modules/local/parascopy/* modules/local/vcf_precheck/* scripts/35-paralogs.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh scripts/validate-setup.sh paralog-checks-and-doc-fixes-1-hba assets/par_grch38.bed assets/paralog_loci.tsv bin/cn_loci.tsv bin/collect_summary.py bin/copy_number.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/paraphase_report.py modules/local/align_minimap2/* modules/local/deepvariant/* modules/local/fastp/* modules/local/paraphase/* modules/local/parascopy/* modules/local/vcf_precheck/* scripts/01b-fastp-qc.sh scripts/35-paralogs.sh scripts/35b-paraphase.sh scripts/ci/check-images.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh scripts/validate-setup.sh @@ -139,7 +139,7 @@ qc-1-sample-qc bin/collect_summary.py bin/mosaic_sex_loss.py scripts/28-multiqc. qc-2-nextflow-sample-qc assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/families.example.tsv assets/live/index.html assets/mito_confirmed.tsv assets/par_grch38.bed assets/pedigree.example.tsv assets/pgx_contraindications.tsv assets/pgx_genes.tsv assets/pgx_inhibitors.tsv assets/stub/db/hla.dat assets/stub/db/kir.dat assets/stub/pedigree.tsv assets/stub/sample.vcf.gz assets/viewer/index.html bin/cn_loci.tsv bin/collect_summary.py bin/copy_number.py bin/cyp2d6_depth_check.py bin/family_report.py bin/family_segments.py bin/hba_copy_number.py bin/ibd_report.py bin/lpa_kiv2.py bin/mosaic_sex_loss.py bin/pgx_mito_merge.py bin/pgx_outside_calls.py bin/pgx_parse.py bin/pgx_refblock_fill.py bin/render_report.py bin/results_index.py modules/local/align_minimap2/* modules/local/cpic_lookup/* modules/local/deepvariant/* modules/local/fastp/* modules/local/hla_typing/* modules/local/html_report/* modules/local/indexcov/* modules/local/mito_haplogroup/* modules/local/mito_variants/* modules/local/mosdepth/* modules/local/pharmcat/* modules/local/pypgx/* modules/local/roh/* modules/local/somalier/* modules/local/t1k_build/* modules/local/telomere_hunter/* modules/local/vcf_precheck/* modules/local/verifybamid2/* scripts/01b-fastp-qc.sh scripts/07-pharmacogenomics.sh scripts/08-hla-typing.sh scripts/10-telomere-hunter.sh scripts/11-roh-analysis.sh scripts/12-mito-haplogroup.sh scripts/16-indexcov.sh scripts/16b-mosdepth.sh scripts/20-mtoolbox.sh scripts/24-html-report.sh scripts/27-cpic-lookup.sh scripts/32-pypgx.sh scripts/33-sample-qc.sh scripts/36-pgx-consensus.sh scripts/38-family.sh scripts/ci/check-images.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh qc-3-cram assets/callability_genes.tsv bin/collect_summary.py modules/local/cram_archive/* modules/local/mosdepth/* modules/local/vcf_precheck/* scripts/16b-mosdepth.sh scripts/34-cram-archive.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh reference-noalt scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh -report-followups-1-sections assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/genotype_rules.tsv assets/mito_confirmed.tsv assets/par_grch38.bed assets/str_extra_loci.json assets/stub/db/catalog.json assets/stub/sample.vcf.gz assets/sv_rule_regions.tsv assets/traits.tsv bin/annotsv_postfilter.py bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/export_consumer.py bin/manta_inv_filter.py bin/render_report.py bin/results_index.py bin/rfc1_motif.py bin/strchive_convert.py bin/sv_svlen_fix.py bin/traits_report.py modules/local/annotsv/* modules/local/expansion_hunter/* modules/local/export_consumer/* modules/local/html_report/* modules/local/manta/* modules/local/traits/* modules/local/vcf_precheck/* scripts/04-manta.sh scripts/05-annotsv.sh scripts/09-expansion-hunter.sh scripts/15-duphold.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/24-html-report.sh scripts/40-export-consumer.sh scripts/42-traits.sh scripts/ci/fake-docker/lib.sh scripts/generate-report.sh scripts/validate-setup.sh tests/fixtures/consumer_export/sites.tsv.gz +report-followups-1-sections assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/genotype_rules.tsv assets/mito_confirmed.tsv assets/par_grch38.bed assets/str_extra_loci.json assets/stub/db/catalog.json assets/stub/sample.vcf.gz assets/sv_rule_regions.tsv assets/traits.tsv bin/annotsv_postfilter.py bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/export_consumer.py bin/family_report.py bin/manta_inv_filter.py bin/render_report.py bin/results_index.py bin/rfc1_motif.py bin/strchive_convert.py bin/sv_svlen_fix.py bin/traits_report.py modules/local/annotsv/* modules/local/expansion_hunter/* modules/local/export_consumer/* modules/local/html_report/* modules/local/manta/* modules/local/traits/* modules/local/vcf_precheck/* scripts/04-manta.sh scripts/05-annotsv.sh scripts/09-expansion-hunter.sh scripts/15-duphold.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/24-html-report.sh scripts/40-export-consumer.sh scripts/42-traits.sh scripts/ci/fake-docker/lib.sh scripts/generate-report.sh scripts/validate-setup.sh tests/fixtures/consumer_export/sites.tsv.gz report-sections-mito-roh-telomere-cnv-1-mito assets/callability_genes.tsv bin/collect_summary.py scripts/10-telomere-hunter.sh scripts/16b-mosdepth.sh scripts/20-mtoolbox.sh scripts/ci/fake-docker/lib.sh scripts/generate-report.sh scripts/validate-setup.sh reports-and-pgx-correctness-clinical assets/construct_variants.tsv bin/clip_check.py bin/comphet_phase.py bin/constraint_join.awk bin/csq_rank.awk scripts/23-clinical-filter.sh scripts/31-slivar.sh scripts/ci/fake-docker/lib.sh reports-and-pgx-correctness-pgx assets/pgx_contraindications.tsv assets/pgx_genes.tsv assets/pgx_inhibitors.tsv bin/cn_loci.tsv bin/collect_summary.py bin/copy_number.py bin/cyp2d6_depth_check.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/pgx_outside_calls.py bin/pgx_parse.py scripts/27-cpic-lookup.sh scripts/32-pypgx.sh scripts/36-pgx-consensus.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh diff --git a/tests/test_collect_summary.py b/tests/test_collect_summary.py index fb800a62..4079a649 100644 --- a/tests/test_collect_summary.py +++ b/tests/test_collect_summary.py @@ -1215,6 +1215,16 @@ def run(name, gt): _, t = run("homref", "0/0") check("sv phase: a deletion genotyped 0/0 pairs with nothing", t == [], t) + # A phased GT with one copy not called says nothing about cis or trans. + check("sv phase: '.|1' gives no copy, '0|1' gives copy 1", + (collect_summary.phase_haps(".|1"), collect_summary.phase_haps("0|1")) == (None, frozenset({1})), + (collect_summary.phase_haps(".|1"), collect_summary.phase_haps("0|1"))) + dl = {"inside": False, "chrom": "chr1", "start": 100, "end": 200, "gt": ".|1", "ps": "90"} + v = {"chrom": "chr1", "pos": 300, "gt": "0|1", "ps": "90", "copies": 1, "ploidy": 2} + verdict, why = collect_summary.pair_phase(dl, v) + check("sv phase: a deletion '.|1' and a hit '0|1' in one phase set are not called in cis", + verdict == "possibly biallelic" and "in cis" not in why, (verdict, why)) + def report_followups_sections(work, schema): """26. step 42's genotype rules, step 9b's motif table of the other From 9f19f150e8c333f523466ae0648c46bf73aaceb9 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 17:24:26 +0100 Subject: [PATCH 03/20] fix(live): pin the nf-weblog plugin in versions.env with a Renovate rule The plugin pin sat in run-all.sh, the one image or plugin version kept outside versions.env, and no Renovate rule saw it. It is now NF_WEBLOG_PLUGIN in versions.env; a custom manager bumps it there and in the example config of docs/live-site.md in one PR; the e2e case reads it from versions.env. Bead: pgp-5bja.23 --- docs/live-site.md | 2 +- renovate.json | 14 ++++++++++++++ scripts/run-all.sh | 3 +-- tests/e2e/live-server-1-weblog.sh | 6 +++--- versions.env | 4 ++++ 5 files changed, 23 insertions(+), 6 deletions(-) diff --git a/docs/live-site.md b/docs/live-site.md index 26018b78..fe03b28e 100644 --- a/docs/live-site.md +++ b/docs/live-site.md @@ -29,7 +29,7 @@ OWNGENOME_LIVE_URL=http://127.0.0.1:8788 ./scripts/run-all.sh your_name male # `./scripts/owngenome-live.sh your_name` does the same and prints the address of that run's page. `hash-password` writes `GENOME_DIR/.owngenome-live/password` with mode 600. The file holds a PBKDF2-SHA256 hash, never the password. Open http://127.0.0.1:8787, log in, and pick the run. -With `OWNGENOME_LIVE_URL` set, `run-all.sh` writes `/nextflow/live.config` and passes it to Nextflow with `-c`. That config loads the nf-weblog plugin (pinned in `run-all.sh`; Nextflow downloads it on the first live run) and points it at `OWNGENOME_LIVE_URL/ingest/`, so every task's progress is posted. The numbered script steps that run after the pipeline (GRIDSS, imputation prep, the somatic caller, the extra callers) post when they start and how they end (the benchmark, `benchmark-variants.sh`, has no step number and posts nothing), and `run-all.sh` posts once more as it ends, with its exit code. A site that is down or slow costs one warning; the run never waits for it and never fails for it. `OWNGENOME_LIVE_URL` must be a plain `http://host:port` URL. +With `OWNGENOME_LIVE_URL` set, `run-all.sh` writes `/nextflow/live.config` and passes it to Nextflow with `-c`. That config loads the nf-weblog plugin (pinned as `NF_WEBLOG_PLUGIN` in `versions.env`; Nextflow downloads it on the first live run) and points it at `OWNGENOME_LIVE_URL/ingest/`, so every task's progress is posted. The numbered script steps that run after the pipeline (GRIDSS, imputation prep, the somatic caller, the extra callers) post when they start and how they end (the benchmark, `benchmark-variants.sh`, has no step number and posts nothing), and `run-all.sh` posts once more as it ends, with its exit code. A site that is down or slow costs one warning; the run never waits for it and never fails for it. `OWNGENOME_LIVE_URL` must be a plain `http://host:port` URL. A bare `nextflow run main.nf` can report too. Write a config like this one, with the sample id at the end of the URL (the page counts only the tasks tagged with it), and pass it with `-c`: diff --git a/renovate.json b/renovate.json index 4cf6b687..a91850c8 100644 --- a/renovate.json +++ b/renovate.json @@ -63,6 +63,20 @@ "extractVersionTemplate": "^v(?.+)$", "versioningTemplate": "regex:^(?\\d+)\\.(?\\d+)\\.(?\\d+)$" }, + { + "customType": "regex", + "description": "nf-weblog, the Nextflow plugin the live site's config loads: NF_WEBLOG_PLUGIN in versions.env and the example config in docs/live-site.md, in one PR. The e2e case live-server-1-weblog reads the pin from versions.env, so it holds no copy to match.", + "managerFilePatterns": [ + "/^versions\\.env$/", + "/^docs/live-site\\.md$/" + ], + "matchStrings": [ + "NF_WEBLOG_PLUGIN=\"nf-weblog@(?\\d+\\.\\d+\\.\\d+)\"", + "id 'nf-weblog@(?\\d+\\.\\d+\\.\\d+)'" + ], + "depNameTemplate": "nextflow-io/nf-weblog", + "datasourceTemplate": "github-releases" + }, { "customType": "regex", "description": "actionlint, downloaded by version in the lint workflow.", diff --git a/scripts/run-all.sh b/scripts/run-all.sh index ec353a40..5ec9d955 100755 --- a/scripts/run-all.sh +++ b/scripts/run-all.sh @@ -296,8 +296,7 @@ else echo " benchmark-variants skipped (only one caller VCF: set EXTRA_CALLERS, export NXF_VER=${NXF_VER:-$NEXTFLOW_VERSION} # The live site (docs/live-site.md). Nextflow's weblog is the nf-weblog plugin, which Nextflow 26.10 loads only when # a config names it, so the block lives in a config of this run's own, passed only when OWNGENOME_LIVE_URL is set: a run -# without it neither downloads the plugin nor warns about an unknown weblog scope. Bump the pin with its release. -NF_WEBLOG_PLUGIN=nf-weblog@1.2.0 +# without it neither downloads the plugin nor warns about an unknown weblog scope. The pin is NF_WEBLOG_PLUGIN in versions.env. if [ -n "${OWNGENOME_LIVE_URL:-}" ]; then # The ingest listener speaks plain HTTP and serves /ingest/ at its root, so https or a path would lose every post. [[ "$OWNGENOME_LIVE_URL" =~ ^http://[A-Za-z0-9.:_-]+/?$ ]] \ diff --git a/tests/e2e/live-server-1-weblog.sh b/tests/e2e/live-server-1-weblog.sh index c98b4e40..75035997 100644 --- a/tests/e2e/live-server-1-weblog.sh +++ b/tests/e2e/live-server-1-weblog.sh @@ -1,15 +1,15 @@ #!/usr/bin/env bash # The live site (bin/owngenome_live.py) on a real Nextflow run: the pinned Nextflow # loads the nf-weblog plugin from the config run-all.sh writes when -# OWNGENOME_LIVE_URL is set (the same text, the pin read from run-all.sh), posts each +# OWNGENOME_LIVE_URL is set (the same text, the pin read from versions.env), posts each # task to the site's ingest port, and the site stores the messages and shows # step 11 (ROH) done in state.json. A control reads state.json before the run: # step 11 is 'not planned' then, so the state check can fail. . "$(dirname "$0")/lib.sh" command -v nextflow >/dev/null || { fail "nextflow is not on PATH"; finish; } -PLUGIN=$(sed -n 's/^NF_WEBLOG_PLUGIN=\([^ ]*\)$/\1/p' "${REPO}/scripts/run-all.sh") -check "run-all.sh pins nf-weblog with a version" has '^nf-weblog@[0-9]+\.[0-9]+\.[0-9]+$' "$PLUGIN" +PLUGIN=$(sed -n 's/^NF_WEBLOG_PLUGIN="\([^"]*\)"$/\1/p' "${REPO}/versions.env") +check "versions.env pins nf-weblog with a version" has '^nf-weblog@[0-9]+\.[0-9]+\.[0-9]+$' "$PLUGIN" G="$GENOME_DIR" ROOT="${CASE_TMP}/root" STATE="${CASE_TMP}/state" diff --git a/versions.env b/versions.env index fb7d3876..b05b4bff 100644 --- a/versions.env +++ b/versions.env @@ -224,6 +224,10 @@ GRIDSS_BLACKLIST_COMMIT="bd7daff8fc32f349752ff713e1fd776361a57d7e" # times larger, has been read and the two steps adapted. # Nextflow release the pipeline is validated with. NEXTFLOW_VERSION="26.10.5" +# The nf-weblog plugin (github.com/nextflow-io/nf-weblog) the live site's +# config loads: run-all.sh writes it into /nextflow/live.config when +# OWNGENOME_LIVE_URL is set, and Nextflow downloads it on the first live run. +NF_WEBLOG_PLUGIN="nf-weblog@1.2.0" # Release of pgsc_calc (github.com/PGScatalog/pgsc_calc) steps 25 and 26 run; # the PGSC_*_IMAGE lines above belong to it. setup.sh pulls it with Nextflow. PGSC_CALC_VERSION="v2.3.0" From c7a136a435a9d03abe181f075089cbdc851728ae Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 17:24:54 +0100 Subject: [PATCH 04/20] fix(report): read the CFTR and UGT1A1 repeats without overclaiming - The setup hint is now the command that runs: STR_EXTRA_LOCI=yes ./scripts/setup.sh --strchive-data . - No 5T copy no longer says the tract changes nothing: the poly-T still sets how R117H reads, so the line reads R117H and names its copy. - A TA length with no named allele gets no function claim and no PharmCAT comparison. - The readings use the traits and CPIC section states: a stale table is read and said to be so, a failed or unreadable one is not read. Bead: pgp-5bja.28 --- bin/collect_summary.py | 103 ++++++++++++++++++++++++---------- docs/09-str-expansions.md | 6 +- tests/test_collect_summary.py | 47 ++++++++++++++-- 3 files changed, 118 insertions(+), 38 deletions(-) diff --git a/bin/collect_summary.py b/bin/collect_summary.py index 6df0cd61..953f7e0d 100755 --- a/bin/collect_summary.py +++ b/bin/collect_summary.py @@ -2099,6 +2099,10 @@ def str_not_graded(d, s): # The repeats STRchive leaves out (assets/str_extra_loci.json), read from # their counts when step 09 sized them (STR_EXTRA_LOCI=yes; docs/09-str-expansions.md). +STR_EXTRA_SETUP = "STR_EXTRA_LOCI=yes ./scripts/setup.sh --strchive-data " +# What a section state other than ok says about a table another step wrote. +OTHER_STEP_STATE = {"stale": "from an earlier run", "failed": "its step failed in the last run", + "unreadable": "unreadable"} # CFTR intron 9 (legacy IVS8): ExpansionHunter writes the TG and the poly-T # tract as two records, so which TG count sits on the copy of which T count # is known only when one of the two is the same on both copies. R117H @@ -2139,13 +2143,25 @@ def trait_genotype(d, s, rsid): return (row.get("status") or "", row.get("genotype") or "") if row else None -def cftr_reading(by, d, s): +def r117h_genotype(d, s, state): + """(genotype or None, why not read, note) of R117H in step 42's traits + table, given the traits section's state: a table of a failed or + unreadable step is not read, one from an earlier run is read and said so.""" + if state in ("failed", "unreadable"): + return None, f"step 42's traits table is {OTHER_STEP_STATE[state]}", "" + r = trait_genotype(d, s, CFTR_R117H) + if not r or r[0] != "called": + return None, "no step 42 traits row for rs78655421 was called", "" + return r[1], "", (" (R117H from step 42's traits table of an earlier run)" if state == "stale" else "") + + +def cftr_reading(by, d, s, traits_state="ok"): tg, t = by.get("CFTR_TG"), by.get("CFTR_T") out = {"locus": "CFTR", "repeat": "intron 9 (legacy IVS8) TG and poly-T tract", "source": CFTR_SOURCE} if not tg or not t: return dict(out, status="not_sized", genotype=".", phase=".", - reading="not in step 09's output: STR_EXTRA_LOCI=yes setup.sh --strchive-data adds the tract " - "to the catalog (docs/09-str-expansions.md).") + reading=f"not in step 09's output: {STR_EXTRA_SETUP} adds the tract to the catalog " + "(docs/09-str-expansions.md).") a, b = str_counts(tg), str_counts(t) if not a or not b: return dict(out, status="not_sized", genotype=".", phase=".", @@ -2161,10 +2177,21 @@ def cftr_reading(by, d, s): other = sorted({n for n in b if n not in (5, 7, 9)}) odd = (f" A poly-T length of {', '.join(f'{n}T' for n in other)} is none of the usual 5T, 7T and 9T." if other else "") + r117, why, note = r117h_genotype(d, s, traits_state) if not five: - return dict(out, five_t=False, phase="not_needed", - reading=f"No 5T copy ({b[0]}T/{b[1]}T), so the tract does not change how a CFTR finding is " - f"read.{odd}") + # No 5T risk of its own, but the poly-T still sets how R117H reads: + # with no 5T copy, an R117H copy sits on a 7T or 9T (or other) copy. + line = (f"No 5T copy ({b[0]}T/{b[1]}T): the tract adds no 5T allele of its own. The poly-T length still " + "matters beside R117H, which is milder and of variable penetrance on a 7T or 9T copy and CF-causing " + "with 5T in cis.") + if r117 is None: + line += f" R117H was not read ({why})." + elif r117.split("/").count("A"): + copies = f"{b[0]}T" if b[0] == b[1] else f"{b[0]}T or {b[1]}T" + line += f" R117H ({r117}) is read: with no 5T copy it sits on a {copies} copy.{note}" + else: + line += f" No R117H copy.{note}" + return dict(out, five_t=False, phase="not_needed", r117h=r117 or ".", reading=line + odd) if len(pairs) == 1: # One TG count on both copies, or 5T on both copies with their two TG counts. tgs = [a[0]] if a[0] == a[1] else a @@ -2174,36 +2201,34 @@ def cftr_reading(by, d, s): line = (f"5T: {five} copy{'s' if five > 1 else ''} with the 5T poly-T tract, {tg_text}. 5T lowers the splicing " "of exon 10 and is a variable-penetrance allele: with TG12 or TG13 on the same copy it more often " "causes CBAVD or another CFTR-related disorder when the other copy carries a CF-causing variant.") - r = trait_genotype(d, s, CFTR_R117H) - if not r or r[0] != "called": + if r117 is None: phase = "not_read" - r117 = (" R117H was not read (no step 42 traits row for rs78655421 was called), so whether it sits on the " - "5T copy is not known.") + text = f" R117H was not read ({why}), so whether it sits on the 5T copy is not known." else: - n = r[1].split("/").count("A") + n = r117.split("/").count("A") if n == 0: phase = "not_needed" - r117 = " No R117H copy: the 5T is read on its own." + text = " No R117H copy: the 5T is read on its own." elif n == 2 or five == 2: phase = "in_cis" - r117 = (f" R117H ({r[1]}) and 5T share a copy: R117H with 5T in cis is a CF-causing allele of " + text = (f" R117H ({r117}) and 5T share a copy: R117H with 5T in cis is a CF-causing allele of " "variable penetrance.") else: phase = "unknown" - r117 = (" One R117H copy and one 5T copy: phase unknown. R117H lies 17.6 kb from the tract, farther " + text = (" One R117H copy and one 5T copy: phase unknown. R117H lies 17.6 kb from the tract, farther " "than a read pair spans, so whether they share a copy (in cis, a CF-causing allele of variable " "penetrance) or not (in trans) is not read; testing a parent tells them apart.") - return dict(out, five_t=True, phase=phase, r117h=r[1] if r and r[0] == "called" else ".", - reading=line + r117 + odd) + text += note + return dict(out, five_t=True, phase=phase, r117h=r117 or ".", reading=line + text + odd) -def ugt1a1_reading(by, d, s): +def ugt1a1_reading(by, d, s, cpic_state="ok"): item = by.get("UGT1A1") out = {"locus": "UGT1A1", "repeat": "promoter TA repeat", "source": UGT1A1_SOURCE} if not item: return dict(out, status="not_sized", genotype=".", pharmcat=".", agreement="not_compared", - reading="not in step 09's output: STR_EXTRA_LOCI=yes setup.sh --strchive-data adds the repeat " - "to the catalog (docs/09-str-expansions.md). PharmCAT (step 7) still reads UGT1A1 *28.") + reading=f"not in step 09's output: {STR_EXTRA_SETUP} adds the repeat to the catalog " + "(docs/09-str-expansions.md). PharmCAT (step 7) still reads UGT1A1 *28.") c = str_counts(item) if not c: return dict(out, status="not_sized", genotype=".", pharmcat=".", agreement="not_compared", @@ -2213,15 +2238,20 @@ def ugt1a1_reading(by, d, s): genotype = "/".join(names) dec = sum(n in UGT1A1_DECREASED for n in names) ta = "/".join(f"TA{n - 1}" for n in c) - reading = (f"{genotype} ({ta}, {c[0]}/{c[1]} TA units as ExpansionHunter counts them): " - + ("no decreased-function TA allele." if not dec else - f"{dec} decreased-function allele{'s' if dec > 1 else ''} (*28 and *37 are decreased function).") - + ("" if all(n in UGT1A1_ALLELES for n in c) else " A TA length outside TA5 to TA8 is not a named allele here.")) out.update(status="read", genotype=genotype, ta_units=c) + reading = f"{genotype} ({ta}, {c[0]}/{c[1]} TA units as ExpansionHunter counts them): " + if not all(n in UGT1A1_ALLELES for n in c): + # A length with no named allele has no function to read, and so no + # count of decreased-function copies to set beside PharmCAT's. + return dict(out, pharmcat=".", agreement="not_compared", + reading=reading + "a TA length outside TA5 to TA8 is not a named allele here, so its function " + "is not read and it is not compared with PharmCAT (step 7).") + reading += ("no decreased-function TA allele." if not dec else + f"{dec} decreased-function allele{'s' if dec > 1 else ''} (*28 and *37 are decreased function).") # PharmCAT's UGT1A1 call (step 27's table, with the reference-block fill # of bin/pgx_refblock_fill.py): its count of *28 and *37 copies. pc, agreement = ".", "not_compared" - p = first_existing(d, [f"cpic/{s}_phenotypes.tsv"]) + p = first_existing(d, [f"cpic/{s}_phenotypes.tsv"]) if cpic_state not in ("failed", "unreadable") else None row = None if p: try: @@ -2235,20 +2265,26 @@ def ugt1a1_reading(by, d, s): # A combination allele reads "*80+*28" or "[*80 + *28]". n = sum(any(x in re.split(r"[+ \[\]]", a) for x in UGT1A1_DECREASED) for a in alleles) agreement = "agrees" if n == dec else "differs" + old = " (step 27's table is from an earlier run)" if cpic_state == "stale" and agreement != "not_compared" else "" if agreement == "agrees": - reading += f" PharmCAT (step 7) called {pc}: the same number of *28 and *37 copies." + reading += f" PharmCAT (step 7) called {pc}: the same number of *28 and *37 copies{old}." elif agreement == "differs": - reading += (f" PharmCAT (step 7) called {pc}, a different number of *28 and *37 copies: read the reads at " - "chr2:233760234 before using either.") + reading += (f" PharmCAT (step 7) called {pc}, a different number of *28 and *37 copies{old}: read the reads " + "at chr2:233760234 before using either.") + elif cpic_state in ("failed", "unreadable"): + reading += f" Not compared with PharmCAT (step 27's table is {OTHER_STEP_STATE[cpic_state]})." else: reading += " Not compared with PharmCAT (no unambiguous UGT1A1 call from step 27)." return dict(out, pharmcat=pc, agreement=agreement, reading=reading) -def str_readings(by, d, s): +def str_readings(by, d, s, states=None): """What the counts of the extra repeats mean, one entry per locus, and the - 5-HTTLPR line that no step can test.""" - return [cftr_reading(by, d, s), ugt1a1_reading(by, d, s), + 5-HTTLPR line that no step can test. states: the traits and cpic section + states, once collect() knows them (a stale, failed or unreadable table is + said so, not read as current).""" + states = states or {} + return [cftr_reading(by, d, s, states.get("traits", "ok")), ugt1a1_reading(by, d, s, states.get("cpic", "ok")), {"locus": "SLC6A4", "repeat": "5-HTTLPR promoter length polymorphism", "status": "not_tested", "genotype": ".", "reading": HTTLPR_READING, "source": "Bousman et al. 2023, PMID 37032427"}] @@ -4795,6 +4831,13 @@ def collect(sample, sample_dir, declared_sex=None): sections[key] = sec cyp2d6 = cyp2d6_block(sample, sample_dir, sections) + # The repeat readings lean on the traits (R117H) and CPIC (UGT1A1) + # tables, whose sections come after the expansions one: read them again + # with those sections' states. + ex = sections["expansions"] + if ex["state"] in ("ok", "stale"): + ex["data"]["readings"] = str_readings({l["locus"]: l for l in ex["data"]["loci"]}, sample_dir, sample, + {k: sections[k]["state"] for k in ("traits", "cpic")}) if sections["locityper"]["state"] in ("ok", "stale"): primary = locityper_primary(sections) for x in sections["locityper"]["data"]["loci"]: diff --git a/docs/09-str-expansions.md b/docs/09-str-expansions.md index 4e2cce94..2261c0a1 100644 --- a/docs/09-str-expansions.md +++ b/docs/09-str-expansions.md @@ -49,9 +49,9 @@ Each entry has a `NoThresholds` reason and no Stranger thresholds: step 9b gives The reports read the counts (`bin/collect_summary.py`, under the repeat expansions section, `readings` in `summary.json`): -- **CFTR**: the TG and poly-T counts of both copies, such as `TG11-7T/TG11-7T`. ExpansionHunter writes the two tracts as two records, so when both the TG and the poly-T counts differ between the copies, which TG sits with which poly-T is not known: the line says `TG11/TG12 5T/7T` and names both pairings. A 5T copy gets the 5T line: 5T is a variable-penetrance allele, more often disease-causing with TG12 or TG13 on the same copy, and it matters when the other copy carries a CF-causing variant. The line then reads R117H (rs78655421), a row of step 42's traits table: no R117H copy; R117H and 5T in cis when one of them is on both copies; otherwise phase unknown, because R117H lies 17.6 kb from the tract, farther than a read pair spans. No 5T copy gives no 5T line. -- **UGT1A1**: the TA count of each copy as its allele (`*36`, `*1`, `*28`, `*37`), and the number of decreased-function copies (`*28`, `*37`) compared with PharmCAT's UGT1A1 call of step 27, which reads `*28` from the small-variant calls with the reference-block fill of `bin/pgx_refblock_fill.py`. The line says whether the two agree; when they differ, read the reads before using either. PharmCAT's call is still the one step 27's drug guidance uses. -- Without the opt-in, both lines say the repeat was not sized and how to add it. +- **CFTR**: the TG and poly-T counts of both copies, such as `TG11-7T/TG11-7T`. ExpansionHunter writes the two tracts as two records, so when both the TG and the poly-T counts differ between the copies, which TG sits with which poly-T is not known: the line says `TG11/TG12 5T/7T` and names both pairings. A 5T copy gets the 5T line: 5T is a variable-penetrance allele, more often disease-causing with TG12 or TG13 on the same copy, and it matters when the other copy carries a CF-causing variant. The line then reads R117H (rs78655421), a row of step 42's traits table: no R117H copy; R117H and 5T in cis when one of them is on both copies; otherwise phase unknown, because R117H lies 17.6 kb from the tract, farther than a read pair spans. No 5T copy gives no 5T line, but the poly-T still matters beside R117H (milder on a 7T or 9T copy, CF-causing with 5T in cis), so the line still reads R117H and says which poly-T copy it sits on. A traits table from an earlier run is read and said to be so; one of a failed or unreadable step 42 is not read. +- **UGT1A1**: the TA count of each copy as its allele (`*36`, `*1`, `*28`, `*37`), and the number of decreased-function copies (`*28`, `*37`) compared with PharmCAT's UGT1A1 call of step 27, which reads `*28` from the small-variant calls with the reference-block fill of `bin/pgx_refblock_fill.py`. The line says whether the two agree; when they differ, read the reads before using either. A TA length outside TA5 to TA8 has no named allele, so the line makes no function claim and does not compare it. A step 27 table from an earlier run is compared and said to be so; one of a failed or unreadable step is not compared. PharmCAT's call is still the one step 27's drug guidance uses. +- Without the opt-in, both lines say the repeat was not sized and give the command that adds it: `STR_EXTRA_LOCI=yes ./scripts/setup.sh --strchive-data `. `tests/test_str_catalog.py` checks each region against the GRCh38 sequence, `tests/test_collect_summary.py` checks the readings on made-up counts (a TG11-7T/TG11-7T sample is the red control of the 5T line), the e2e case `genotype-rule-rows-1-str-extra` sizes both tracts on made-up reads with ExpansionHunter 5.0.0, and the e2e case `str-extra-repeats-1-readings` runs step 09 and both reports on made-up reads and checks the lines. diff --git a/tests/test_collect_summary.py b/tests/test_collect_summary.py index 4079a649..5f0f6e4a 100644 --- a/tests/test_collect_summary.py +++ b/tests/test_collect_summary.py @@ -1328,7 +1328,8 @@ def str_extra_readings(work, schema): import html as htmlmod fmt = "GT:REPCN" - def made_up(name, tg, t, ugt, r117h=None, ugt_pc=None, stranger=False): + def made_up(name, tg, t, ugt, r117h=None, ugt_pc=None, stranger=False, traits_text=None, side_age=0, + run_status=None): x = os.path.join(work, "str-extra", name) eh = (f"chr7\t117548606\t.\tT\t\t.\tPASS\tEND=117548628;REPID=CFTR;VARID=CFTR_TG\t{fmt}\t0/1:{tg}\n" f"chr7\t117548628\t.\tT\t\t.\tPASS\tEND=117548635;REPID=CFTR;VARID=CFTR_T\t{fmt}\t0/1:{t}\n" @@ -1342,12 +1343,15 @@ def made_up(name, tg, t, ugt, r117h=None, ugt_pc=None, stranger=False): put(f"{x}/expansion_hunter/{name}_eh_catalog.json", json.dumps([ {"LocusId": "CFTR", "NoThresholds": "the CFTR tract"}, {"LocusId": "UGT1A1", "NoThresholds": "the UGT1A1 repeat"}])) - if r117h: - put(f"{x}/traits/{name}_traits.tsv", "\t".join(collect_summary.traits_report.OUT_COLUMNS) + "\n" - + f"rs78655421\tCFTR\tchr7\t117530975\tCFTR R117H\t{r117h}\tcalled\tr\tno\tn\t32404922\n") + if r117h or traits_text: + put(f"{x}/traits/{name}_traits.tsv", traits_text or "\t".join(collect_summary.traits_report.OUT_COLUMNS) + + "\n" + f"rs78655421\tCFTR\tchr7\t117530975\tCFTR R117H\t{r117h}\tcalled\tr\tno\tn\t32404922\n", + age_days=side_age) if ugt_pc: put(f"{x}/cpic/{name}_phenotypes.tsv", "Gene\tDiplotype\tPhenotype\tStatus\n" - f"UGT1A1\t{ugt_pc}\tsome function\tnon-normal\n") + f"UGT1A1\t{ugt_pc}\tsome function\tnon-normal\n", age_days=side_age) + if run_status: + put(f"{x}/logs/run_status.tsv", run_status) summ, txt, html_ = render(x, name) ed = summ["sections"]["expansions"]["data"] return summ, ed, {r["locus"]: r for r in ed.get("readings") or []}, txt, htmlmod.unescape(html_) @@ -1413,6 +1417,39 @@ def made_up(name, tg, t, ugt, r117h=None, ugt_pc=None, stranger=False): [rd[k]["status"] for k in ("CFTR", "UGT1A1", "SLC6A4")] == ["not_sized", "not_sized", "not_tested"] and "STR_EXTRA_LOCI=yes" in rd["CFTR"]["reading"] and "CFTR intron 9 (legacy IVS8) TG and poly-T tract: " "not sized" in txt, rd) + check("extra repeats: the setup hint is a command that runs, with the script path and the genome folder", + all("STR_EXTRA_LOCI=yes ./scripts/setup.sh --strchive-data " in rd[k]["reading"] + for k in ("CFTR", "UGT1A1")), rd) + + # No 5T copy: the poly-T still matters beside R117H, so the line must not say it changes nothing. + _, _, rd, _, _ = made_up("EXH", "11/11", "7/9", "7/7", r117h="G/A") + c = rd["CFTR"] + check("extra repeats: no 5T copy with one R117H copy says R117H sits on the 7T or 9T copy", + c.get("five_t") is False and "R117H (G/A) is read" in c["reading"] and "7T or 9T copy" in c["reading"] + and "does not change how a CFTR finding is read" not in c["reading"] and c.get("r117h") == "G/A", c) + _, _, rd, _, _ = made_up("EXI", "11/11", "7/7", "7/7") + check("extra repeats: no 5T copy and no traits table says R117H was not read", + "R117H was not read" in rd["CFTR"]["reading"], rd["CFTR"]) + + # A TA length with no named allele: no function claim and no PharmCAT comparison. + _, _, rd, _, _ = made_up("EXJ", None, None, "7/11", ugt_pc="*1/*1") + u = rd["UGT1A1"] + check("extra repeats: UGT1A1 TA6/TA10 makes no function claim and is not compared with PharmCAT", + u.get("genotype") == "*1/TA10" and u.get("agreement") == "not_compared" + and "decreased-function" not in u["reading"] and "not compared with PharmCAT" in u["reading"], u) + + # The traits and CPIC tables are read with their sections' states. + bad_traits = "rsid\tgenotype\tstatus\nrs78655421\tG/A\tcalled\n" + _, _, rd, _, _ = made_up("EXK", "12/12", "5/7", "7/8", traits_text=bad_traits) + check("extra repeats: an unreadable traits table gives no R117H, and says why", + rd["CFTR"].get("phase") == "not_read" and "traits table is unreadable" in rd["CFTR"]["reading"], rd["CFTR"]) + status = f"meta\tstarted_epoch\t{time.time() - DAY}\nstep\t09\tok\nstep\t27\tskipped\nstep\t42\tskipped\n" + _, _, rd, _, _ = made_up("EXL", "12/12", "5/7", "7/8", r117h="G/A", ugt_pc="*1/*28", side_age=3, + run_status=status) + check("extra repeats: R117H and PharmCAT from tables of an earlier run are said to be so", + "R117H from step 42's traits table of an earlier run" in rd["CFTR"]["reading"] + and rd["UGT1A1"]["agreement"] == "agrees" and "step 27's table is from an earlier run" in rd["UGT1A1"]["reading"], + rd) SMN_COLS = ("chrom\tstart\tend\tlocus\tagCN_filter\tagCN\tagCN_qual\tpsCN_filter\tpsCN\tpsCN_qual\t" From 70c074227012ff5ca03828ff4c566beabf870179 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 17:25:16 +0100 Subject: [PATCH 05/20] docs: say what the report search, the Nextflow report and the tests really cover - 24-html-report: a section's locus name (repeat, Locityper, copy-number locus) is indexed too, not only gene symbols; the Nextflow report lacks only slivar's compound-het counts. - testing: nf-wiring-2-sv-str-tools runs the Nextflow duphold, CNVpytor and ExpansionHunter modules with their real images. - 53-locityper: the image test repeats the Locityper commands of the database build and the step; it does not run setup.sh or step 53. Bead: pgp-5bja.27 Bead: pgp-5bja.28 --- docs/24-html-report.md | 4 ++-- docs/53-locityper.md | 2 +- docs/testing.md | 6 +++--- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/24-html-report.md b/docs/24-html-report.md index 5efddee2..31ab1503 100644 --- a/docs/24-html-report.md +++ b/docs/24-html-report.md @@ -68,7 +68,7 @@ After them, "What this report leaves out" lists the steps not run and what no st - **Links to the step's own files**, relative to the report, when the summary says they exist: the PharmCAT HTML report, the CPIC recommendations per gene, the CPSR HTML report (or the `cpsr/` folder that holds it) and the Exomiser ranking table. The report sits in the sample folder beside each step's folder, in the bash layout and in a Nextflow `--outdir`. MultiQC and the couple report are not linked: the summary does not record them, and the couple report sits in the family folder. - **Nothing is fetched.** A Content-Security-Policy (`default-src 'none'; img-src data:; style-src 'unsafe-inline'`) lets the browser load nothing but the file itself and its inline pictures, and there is no JavaScript. The only addresses off this computer are the haplogroup links (FamilyTreeDNA, YFull, the docs) and the project link in the footer, each marked "opens a website"; the FamilyTreeDNA and YFull ones add that the address carries the haplogroup. The genome viewer links go to `127.0.0.1`, this computer. - **A Viewer column on every findings table.** The ClinVar hits, the secondary-findings genes, the mitochondrial variants, the known deep intronic sites, the cDNA scan, the AnnotSV class 4 and 5 SVs and the calls moved to class 3, slivar's pairs in trans, the mosaic alterations and Exomiser's ranked variants each end with a Viewer column: one link per row, to that row's position (a range for an SV, a region or a pair). The link opens the reads there in the [genome viewer](viewer.md) once `scripts/owngenome-view.sh ` runs; a note under the table says so. The column is new, so every other cell keeps its text; it replaces the single line of links that sat under three of these tables. A row with no position has an empty Viewer cell. On the [live site](live-site.md) the served report points the same links at the site's own viewer of the run. -- **Genes and rsIDs in these results.** The closing list is built by `bin/results_index.py`. It holds every gene symbol and rsID that a section in state complete or stale names in `summary.json`, by the same rule `scripts/owngenome-lookup.sh` uses for "in the reports" (`names_of` in `bin/genome_lookup.py`). It adds the genes of the clinical filter (step 23) and variant prioritization (step 31), whose cards only count them: from their summary tables on the bash path, from the VCFs on the Nextflow path. The two give the same genes when a variant's most severe consequence is on the transcript the report ranks first (MANE, else canonical); the step 31 table keeps the gene of the most severe consequence, so where an overlapping gene's transcript carries the worse one, the two forms name different genes. A ClinVar hit's rsIDs come from step 13's VEP output (`Existing_variation`), because the hits file's ID column holds ClinVar's own id; without step 13 a hit has none. The Nextflow `HTML_REPORT` process does not get step 13's VEP output, so its index has no ClinVar rsIDs; step 24, which `run-all.sh` runs over the output folder afterwards, adds them. A failed, missing or unreadable section adds nothing. The browser's find (Ctrl+F) searches the list; there is still no script. A position is not in it: ask `scripts/owngenome-lookup.sh` for a genotype ([genotype lookup](lookup.md)). +- **Genes and rsIDs in these results.** The closing list is built by `bin/results_index.py`. It holds every gene symbol and rsID that a section in state complete or stale names in `summary.json`, by the same rule `scripts/owngenome-lookup.sh` uses for "in the reports" (`names_of` in `bin/genome_lookup.py`). A section's locus name counts as a term too, when it holds an upper-case letter: a repeat expansion (HTT, FMR1, CFTR_TG), a Locityper locus (HLA-DRB1, FCGR) and a copy-number locus are indexed by their `locus` field, so a term is not always a gene symbol. It adds the genes of the clinical filter (step 23) and variant prioritization (step 31), whose cards only count them: from their summary tables on the bash path, from the VCFs on the Nextflow path. The two give the same genes when a variant's most severe consequence is on the transcript the report ranks first (MANE, else canonical); the step 31 table keeps the gene of the most severe consequence, so where an overlapping gene's transcript carries the worse one, the two forms name different genes. A ClinVar hit's rsIDs come from step 13's VEP output (`Existing_variation`), because the hits file's ID column holds ClinVar's own id; without step 13 a hit has none. The Nextflow `HTML_REPORT` process does not get step 13's VEP output, so its index has no ClinVar rsIDs; step 24, which `run-all.sh` runs over the output folder afterwards, adds them. A failed, missing or unreadable section adds nothing. The browser's find (Ctrl+F) searches the list; there is still no script. A position is not in it: ask `scripts/owngenome-lookup.sh` for a genotype ([genotype lookup](lookup.md)). - **Dark on screen, light on paper.** - **Print one section.** Open a section from the contents (or its "link to this section"), then print: only that section prints, with the header (sample, generation date, pipeline commit, ClinVar file date, ClinVar in CPSR's bundle, HLA database), the list of what is not assessed and the disclaimer. A chapter prints the same way from its title in the contents. Without an anchor in the address the whole report prints, as it does in a browser without CSS `:has()`. @@ -177,4 +177,4 @@ start ${GENOME_DIR}/${SAMPLE}/${SAMPLE}_report.html - Steps that were not run show "Not run" and are listed under Steps Not Run — this is expected. A step recorded as failed shows "Failed" instead - The report contains health findings: the ClinVar table lists your pathogenic and likely pathogenic variants, and other sections summarise pharmacogenomics, CPSR and the other steps. It holds no raw reads or full VCF, but share it only as you would share a medical record - Re-run this script anytime to update the report after running additional steps -- The Nextflow `HTML_REPORT` module runs the same `bin/render_report.py` on the outputs of its run's QC, ClinVar, PharmCAT, CPIC, CPSR, clinical filter, slivar, ROH and mitochondrial haplogroup steps, so both reports print the same numbers for the same file. Its clinical filter and slivar cards show counts only. For every section on a Nextflow run, run this script with `GENOME_DIR` set to the `--outdir` +- The Nextflow `HTML_REPORT` module runs the same `bin/render_report.py` on the outputs of its run's QC, ClinVar, PharmCAT, CPIC, CPSR, clinical filter, slivar, ROH and mitochondrial haplogroup steps, so both reports print the same numbers for the same file. Its slivar card has no compound-het variant and gene counts: the Nextflow slivar module writes no `_compound_hets.tsv`. For every section on a Nextflow run, run this script with `GENOME_DIR` set to the `--outdir` diff --git a/docs/53-locityper.md b/docs/53-locityper.md index 5b4f415d..9c89183c 100644 --- a/docs/53-locityper.md +++ b/docs/53-locityper.md @@ -88,7 +88,7 @@ Both reports carry a section "Locityper Second Opinion at Hard Loci" in the medi ## How it was checked -The image test (`tests/smoke/locityper.sh`, the `LOCITYPER_IMAGE` row of `tests/smoke/commands.tsv`) runs the database build and the step on the fixture's GIAB HG002 reads at HLA-DRB1, HLA-DQB1 and CYP2D6. The panel holds HG002's own two assembled haplotypes, which are the truth. Both called haplotypes at each locus must match them at QV 20 or more (one difference per 100 bases or fewer), aligned with `locityper align`. That QV measures how close each called haplotype is to HG002's true one; it is not the `quality` column, Locityper's own genotype quality. +The image test (`tests/smoke/locityper.sh`, the `LOCITYPER_IMAGE` row of `tests/smoke/commands.tsv`) repeats, inside the image, the Locityper commands that `setup.sh --locityper-data` and step 53 run (the database build, `preproc` and `genotype`) on the fixture's GIAB HG002 reads at HLA-DRB1, HLA-DQB1 and CYP2D6. It does not run `setup.sh` or `scripts/53-locityper.sh` themselves. The panel holds HG002's own two assembled haplotypes, which are the truth. Both called haplotypes at each locus must match them at QV 20 or more (one difference per 100 bases or fewer), aligned with `locityper align`. That QV measures how close each called haplotype is to HG002's true one; it is not the `quality` column, Locityper's own genotype quality. The red control takes out of the panel HG002's haplotypes and every other haplotype within QV 20 of them, rebuilds the database and genotypes again; the same check must then fail at every locus. Leaving HG002 out with `--leave-out` is not enough: other people in the panel carry near copies of HG002's haplotypes, and the call stays close. diff --git a/docs/testing.md b/docs/testing.md index fdabfd79..43b14519 100644 --- a/docs/testing.md +++ b/docs/testing.md @@ -168,15 +168,15 @@ This is every step script that no e2e case runs, plus the full-size index build | 02b (long-read alignment) | long reads; the e2e fixture has none | the minimap2 image's smoke test, which maps synthetic HiFi reads with `-x map-hifi`; nothing runs the script | | 03e (Clair3) | long reads; the e2e fixture has none | the image's smoke test, with a recall check on synthetic HiFi reads and a check that the two model directories the script names exist; the fake-docker suite runs the script | | 04c (Sniffles2) | long reads; the e2e fixture has none | the image's smoke test on the synthetic HiFi alignment; `tests/fake-docker/bash-step-partial-outputs.sh` runs the script | -| 15 (duphold) | no case yet | the image's smoke test on a planted deletion; `tests/fake-docker/bash-step-sv.sh` runs the script; the Nextflow module is stub-run | -| 18 (CNVpytor) | no case yet; it needs its own resource files | the image's smoke test; the SV merge case reads a planted CNVpytor VCF, not one CNVpytor wrote | +| 15 (duphold) | the bash step: no case yet | case `nf-wiring-2-sv-str-tools` runs the Nextflow `DUPHOLD` and `DUPHOLD_FILTER` with the real image on the fixture's Manta calls; the image's smoke test on a planted deletion; `tests/fake-docker/bash-step-sv.sh` runs the script | +| 18 (CNVpytor) | the bash step: no case yet | case `nf-wiring-2-sv-str-tools` runs the Nextflow `CNVPYTOR` and `CNVPYTOR_VCF` with the real image and the pinned GC and mask files; the image's smoke test; the SV merge case reads a planted CNVpytor VCF, not one CNVpytor wrote | | 29 (somatic Mutect2) | no case yet | `tests/fake-docker/bash-step-somatic-mito.sh` checks the Mutect2 command line | | 03c (Strelka2) | no case yet | the image's smoke test, with a recall check against the truth slice | | 03d (Octopus) | no case yet | the image's smoke test, with a recall check against the truth slice; case `bash-step-bash42` runs the script under bash 4.2 against the fake docker | | 04a (TIDDIT) | no case yet | the image's smoke test; settle-doubts question 4 runs the script with only a BWA-MEM2 index; case `bash-step-bash42` and `tests/fake-docker/bash-step-sv.sh` run it against the fake docker | | chip-to-vcf (array data to VCF) on a vendor download | case `consumer-export-1` runs it only on the 23andMe and AncestryDNA files step 40 wrote | `tests/test_chip_to_vcf.sh` in the unit-tests job runs the script with the real bcftools and Picard images on two ten-row vendor files | -The Nextflow ExpansionHunter module is only stub-run (`-profile test_all -stub` in `nextflow.yml`); the bash step 09 runs in case `bash-step-09-str`. +Case `nf-wiring-2-sv-str-tools` runs the Nextflow `EXPANSION_HUNTER` module with the real image on the fixture's BAM (the bundled catalog cut to the fixture's contigs); the bash step 09 runs in case `bash-step-09-str`. ## The settle-doubts job From 7a4c99e1fe449835e619c89a83001c10232e2997 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 17:27:45 +0100 Subject: [PATCH 06/20] fix(schema): name the cram,crai and gvcf,gvcf_index samplesheet columns in --input main.nf reads both pairs, but the schema's --input description named only the FASTQ, BAM and VCF columns. Bead: pgp-5bja.29 --- nextflow_schema.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nextflow_schema.json b/nextflow_schema.json index 6795fd58..d22fecf3 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -15,7 +15,7 @@ "type": "string", "format": "file-path", "mimetype": "text/csv", - "description": "Path to comma-separated samplesheet: sample, then fastq_1,fastq_2 (aligned and called), or bam,bam_index (called), or vcf,vcf_index (a VCF with one sample) with an optional bam,bam_index; sex (male or female) is required on every row that is called and when expansion_hunter is enabled.", + "description": "Path to comma-separated samplesheet: sample, then fastq_1,fastq_2 (aligned and called), or bam,bam_index (called), or vcf,vcf_index (a VCF with one sample) with an optional bam,bam_index; sex (male or female) is required on every row that is called and when expansion_hunter is enabled. A row can also give cram,crai instead of the BAM, and gvcf,gvcf_index beside the VCF.", "fa_icon": "fas fa-file-csv" }, "outdir": { From 7d310a41b5c6cd73d8c9504e14bc2c36b413c673 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 17:28:00 +0100 Subject: [PATCH 07/20] test(docs): fail when a schema parameter has no row in docs/parameters.md The page is written by hand from nextflow_schema.json, so a new parameter could ship without a row. The test also fails on a row for a parameter the schema lacks (main.nf's --pedigree is allowed by name), and runs a dropped row and a planted row as controls. Bead: pgp-5bja.30 --- tests/test_parameters_doc.py | 76 ++++++++++++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100644 tests/test_parameters_doc.py diff --git a/tests/test_parameters_doc.py b/tests/test_parameters_doc.py new file mode 100644 index 00000000..a04274e7 --- /dev/null +++ b/tests/test_parameters_doc.py @@ -0,0 +1,76 @@ +#!/usr/bin/env python3 +"""docs/parameters.md lists every parameter of nextflow_schema.json. + +Checks: + 1. every parameter of the schema (each group's properties, and any + top-level ones) has a row "| `--` |" in docs/parameters.md; + 2. every row of the page names a schema parameter, or a parameter main.nf + reads that the schema does not hold yet (DOC_ONLY, each checked against + main.nf), so a removed or renamed parameter cannot keep its row; + 3. controls: the page with one row taken out fails check 1 and names that + parameter; the page with a planted row fails check 2. +""" +import json +import os +import re +import sys + +ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) +SCHEMA = os.path.join(ROOT, "nextflow_schema.json") +PAGE = os.path.join(ROOT, "docs", "parameters.md") +MAIN = os.path.join(ROOT, "main.nf") +# Read by main.nf, not (yet) in the schema; the page writes their rows by hand. +DOC_ONLY = {"pedigree"} +ROW = re.compile(r"^\| `--([A-Za-z0-9_]+)` \|", re.M) + +failures = 0 + + +def check(desc, ok, detail=""): + global failures + print(f"[{'PASS' if ok else 'FAIL'}] {desc}" + ("" if ok else f" -- {detail}")) + if not ok: + failures += 1 + + +def schema_params(schema): + names = set(schema.get("properties", {})) + for group in schema.get("definitions", schema.get("$defs", {})).values(): + names |= set(group.get("properties", {})) + return names + + +def problems(params, page_text, main_text): + rows = ROW.findall(page_text) + missing = sorted(params - set(rows)) + extra = sorted(r for r in set(rows) - params + if r not in DOC_ONLY or not re.search(rf"\bparams\.{r}\b|'{r}'", main_text)) + return missing, extra + + +with open(SCHEMA, encoding="utf-8") as f: + params = schema_params(json.load(f)) +with open(PAGE, encoding="utf-8") as f: + page = f.read() +with open(MAIN, encoding="utf-8") as f: + main_nf = f.read() + +check("the schema has parameters to compare", len(params) > 20, len(params)) +missing, extra = problems(params, page, main_nf) +check("every schema parameter has a row in docs/parameters.md", not missing, f"no row for: {', '.join(missing)}") +check("every row of docs/parameters.md is a schema parameter or one main.nf reads", not extra, + f"not in the schema: {', '.join(extra)}") +check("each parameter has one row", len(ROW.findall(page)) == len(set(ROW.findall(page))), + sorted(r for r in set(ROW.findall(page)) if ROW.findall(page).count(r) > 1)) + +# Controls: the same comparison must catch a dropped row and a planted one. +victim = sorted(params)[0] +dropped = re.sub(rf"^\| `--{victim}` \|.*\n", "", page, flags=re.M) +check(f"control: the page without the --{victim} row fails, naming it", + problems(params, dropped, main_nf)[0] == [victim], problems(params, dropped, main_nf)) +planted = page + "| `--no_such_parameter` | none | planted |\n" +check("control: a planted row for a parameter the schema lacks fails", + problems(params, planted, main_nf)[1] == ["no_such_parameter"], problems(params, planted, main_nf)) + +print(f"\nRESULT: {'ALL PASS' if not failures else f'{failures} FAILED'}") +sys.exit(1 if failures else 0) From d85e9d5759bf05d3d6aeec7dfb54cd3a3b8d2067 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 18:07:56 +0100 Subject: [PATCH 08/20] fix(report): say a homozygous R117H with no 5T copy sits on both copies --- bin/collect_summary.py | 3 +++ tests/test_collect_summary.py | 9 +++++++++ 2 files changed, 12 insertions(+) diff --git a/bin/collect_summary.py b/bin/collect_summary.py index 953f7e0d..25f2e285 100755 --- a/bin/collect_summary.py +++ b/bin/collect_summary.py @@ -2186,6 +2186,9 @@ def cftr_reading(by, d, s, traits_state="ok"): "with 5T in cis.") if r117 is None: line += f" R117H was not read ({why})." + elif r117.split("/").count("A") == 2: + copies = f"both {b[0]}T copies" if b[0] == b[1] else f"both copies ({b[0]}T and {b[1]}T)" + line += f" R117H ({r117}) is read: it sits on {copies}.{note}" elif r117.split("/").count("A"): copies = f"{b[0]}T" if b[0] == b[1] else f"{b[0]}T or {b[1]}T" line += f" R117H ({r117}) is read: with no 5T copy it sits on a {copies} copy.{note}" diff --git a/tests/test_collect_summary.py b/tests/test_collect_summary.py index 5f0f6e4a..5ef81ccc 100644 --- a/tests/test_collect_summary.py +++ b/tests/test_collect_summary.py @@ -1427,6 +1427,15 @@ def made_up(name, tg, t, ugt, r117h=None, ugt_pc=None, stranger=False, traits_te check("extra repeats: no 5T copy with one R117H copy says R117H sits on the 7T or 9T copy", c.get("five_t") is False and "R117H (G/A) is read" in c["reading"] and "7T or 9T copy" in c["reading"] and "does not change how a CFTR finding is read" not in c["reading"] and c.get("r117h") == "G/A", c) + # Homozygous R117H with no 5T copy sits on both copies, not on one or the other. + _, _, rd, _, _ = made_up("EXK", "11/11", "7/9", "7/7", r117h="A/A") + c = rd["CFTR"] + check("extra repeats: no 5T copy with two R117H copies says R117H sits on both copies (7T and 9T)", + c.get("five_t") is False and "R117H (A/A) is read: it sits on both copies (7T and 9T)" in c["reading"] + and "it sits on a 7T or 9T copy" not in c["reading"], c) + _, _, rd, _, _ = made_up("EXL", "11/11", "7/7", "7/7", r117h="A/A") + check("extra repeats: two R117H copies on two 7T copies says both 7T copies", + "R117H (A/A) is read: it sits on both 7T copies" in rd["CFTR"]["reading"], rd["CFTR"]) _, _, rd, _, _ = made_up("EXI", "11/11", "7/7", "7/7") check("extra repeats: no 5T copy and no traits table says R117H was not read", "R117H was not read" in rd["CFTR"]["reading"], rd["CFTR"]) From 6d8c717987198d936debef029da04083a4d199b3 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 18:33:17 +0100 Subject: [PATCH 09/20] docs(report): say a chrX, chrY or chrM expansion position can be a search term --- docs/24-html-report.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/24-html-report.md b/docs/24-html-report.md index 31ab1503..ee465b71 100644 --- a/docs/24-html-report.md +++ b/docs/24-html-report.md @@ -68,7 +68,7 @@ After them, "What this report leaves out" lists the steps not run and what no st - **Links to the step's own files**, relative to the report, when the summary says they exist: the PharmCAT HTML report, the CPIC recommendations per gene, the CPSR HTML report (or the `cpsr/` folder that holds it) and the Exomiser ranking table. The report sits in the sample folder beside each step's folder, in the bash layout and in a Nextflow `--outdir`. MultiQC and the couple report are not linked: the summary does not record them, and the couple report sits in the family folder. - **Nothing is fetched.** A Content-Security-Policy (`default-src 'none'; img-src data:; style-src 'unsafe-inline'`) lets the browser load nothing but the file itself and its inline pictures, and there is no JavaScript. The only addresses off this computer are the haplogroup links (FamilyTreeDNA, YFull, the docs) and the project link in the footer, each marked "opens a website"; the FamilyTreeDNA and YFull ones add that the address carries the haplogroup. The genome viewer links go to `127.0.0.1`, this computer. - **A Viewer column on every findings table.** The ClinVar hits, the secondary-findings genes, the mitochondrial variants, the known deep intronic sites, the cDNA scan, the AnnotSV class 4 and 5 SVs and the calls moved to class 3, slivar's pairs in trans, the mosaic alterations and Exomiser's ranked variants each end with a Viewer column: one link per row, to that row's position (a range for an SV, a region or a pair). The link opens the reads there in the [genome viewer](viewer.md) once `scripts/owngenome-view.sh ` runs; a note under the table says so. The column is new, so every other cell keeps its text; it replaces the single line of links that sat under three of these tables. A row with no position has an empty Viewer cell. On the [live site](live-site.md) the served report points the same links at the site's own viewer of the run. -- **Genes and rsIDs in these results.** The closing list is built by `bin/results_index.py`. It holds every gene symbol and rsID that a section in state complete or stale names in `summary.json`, by the same rule `scripts/owngenome-lookup.sh` uses for "in the reports" (`names_of` in `bin/genome_lookup.py`). A section's locus name counts as a term too, when it holds an upper-case letter: a repeat expansion (HTT, FMR1, CFTR_TG), a Locityper locus (HLA-DRB1, FCGR) and a copy-number locus are indexed by their `locus` field, so a term is not always a gene symbol. It adds the genes of the clinical filter (step 23) and variant prioritization (step 31), whose cards only count them: from their summary tables on the bash path, from the VCFs on the Nextflow path. The two give the same genes when a variant's most severe consequence is on the transcript the report ranks first (MANE, else canonical); the step 31 table keeps the gene of the most severe consequence, so where an overlapping gene's transcript carries the worse one, the two forms name different genes. A ClinVar hit's rsIDs come from step 13's VEP output (`Existing_variation`), because the hits file's ID column holds ClinVar's own id; without step 13 a hit has none. The Nextflow `HTML_REPORT` process does not get step 13's VEP output, so its index has no ClinVar rsIDs; step 24, which `run-all.sh` runs over the output folder afterwards, adds them. A failed, missing or unreadable section adds nothing. The browser's find (Ctrl+F) searches the list; there is still no script. A position is not in it: ask `scripts/owngenome-lookup.sh` for a genotype ([genotype lookup](lookup.md)). +- **Genes and rsIDs in these results.** The closing list is built by `bin/results_index.py`. It holds every gene symbol and rsID that a section in state complete or stale names in `summary.json`, by the same rule `scripts/owngenome-lookup.sh` uses for "in the reports" (`names_of` in `bin/genome_lookup.py`). A section's locus name counts as a term too, when it holds an upper-case letter: a repeat expansion (HTT, FMR1, CFTR_TG), a Locityper locus (HLA-DRB1, FCGR) and a copy-number locus are indexed by their `locus` field, so a term is not always a gene symbol. An expansion record with no REPID or VARID takes its position as its locus (`chrX:147912050`), and that counts as a term only when it holds an upper-case letter (a chrX, chrY or chrM position). It adds the genes of the clinical filter (step 23) and variant prioritization (step 31), whose cards only count them: from their summary tables on the bash path, from the VCFs on the Nextflow path. The two give the same genes when a variant's most severe consequence is on the transcript the report ranks first (MANE, else canonical); the step 31 table keeps the gene of the most severe consequence, so where an overlapping gene's transcript carries the worse one, the two forms name different genes. A ClinVar hit's rsIDs come from step 13's VEP output (`Existing_variation`), because the hits file's ID column holds ClinVar's own id; without step 13 a hit has none. The Nextflow `HTML_REPORT` process does not get step 13's VEP output, so its index has no ClinVar rsIDs; step 24, which `run-all.sh` runs over the output folder afterwards, adds them. A failed, missing or unreadable section adds nothing. The browser's find (Ctrl+F) searches the list; there is still no script. Apart from those locus names a position is not in it, and the search does not answer a position: ask `scripts/owngenome-lookup.sh` for a genotype ([genotype lookup](lookup.md)). - **Dark on screen, light on paper.** - **Print one section.** Open a section from the contents (or its "link to this section"), then print: only that section prints, with the header (sample, generation date, pipeline commit, ClinVar file date, ClinVar in CPSR's bundle, HLA database), the list of what is not assessed and the disclaimer. A chapter prints the same way from its title in the contents. Without an anchor in the address the whole report prints, as it does in a browser without CSS `:has()`. From b8073df36c719e5e9f32168af99bee9203063479 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 19:39:10 +0100 Subject: [PATCH 10/20] feat(traits): read copy-aware markers from the reads at every GRCh38 copy of the site The LPA KIV-2 4925G>A and SMN1 g.27134T>G sites sit in sequence the genome holds more than once, so a gVCF genotype cannot read them. MARKERS pins each one's dbSNP b157 record and the GRCh38 context of every copy; the reads at all copies are added up and read as a fraction of copies. The SMN1 reading sits beside step 35's psCN and always says it does not prove 2+0. Bead: pgp-pfw6.21 --- bin/traits_report.py | 268 ++++++++++++++++++++++++++++++++++++ tests/test_traits_report.py | 117 +++++++++++++++- 2 files changed, 384 insertions(+), 1 deletion(-) diff --git a/bin/traits_report.py b/bin/traits_report.py index 6b56e039..e23c5099 100755 --- a/bin/traits_report.py +++ b/bin/traits_report.py @@ -36,6 +36,27 @@ and exits 1 on the first list of problems; tests/test_traits_table.py and tests/test_traits_report.py run the same checks. Standard library only (it runs in the plain python image and in the WhatsHap image). + +Copy-aware markers (MARKERS): two well-sourced markers of the rule catalog's +list sit in sequence the genome holds more than once, so a gVCF genotype at +one position cannot read them and they are not rule lines. Each is read from +the BAM instead, by the step that counts its copies: the reads at every +GRCh38 copy of the site are added up, and the share that carries the marker's +base is a fraction of copies, not a diploid genotype. + LPA_KIV2_4925 the LPA KIV-2 4925G>A splice variant (Coassin et al. 2017, + PMID 28444229), read by step 46 beside the KIV-2 count + (bin/lpa_kiv2.py null-allele) + SMN1_27134 SMN1 g.27134T>G (rs143838139; Luo et al. 2014, PMID + 23788250), read by step 35 beside the SMN1 psCN + + traits_report.py --marker-regions NAME + print the marker's copy sites as samtools regions (chrom:pos-pos) + traits_report.py --marker SMN1_27134 --pileup P [--copy-table SMN.tsv] --marker-out OUT.tsv + read `samtools mpileup` at those sites (no reference, every MAPQ) and + write OUT.tsv (metricvalue). SMN.tsv is step 35's SMN table: the + SMN1 psCN at exon 7 and the agCN at the marker site come from it. Exit + 0 whatever the status (an unreadable pileup is 'not assessed'); 2 when + OUT.tsv cannot be written. """ import argparse import itertools @@ -504,6 +525,228 @@ def report(rows, recs, rules=None, phased=None, lines=None): return out +# --- Copy-aware markers -------------------------------------------------------------- + +# Each marker: the dbSNP record of one copy (rsID, chrom, POS, REF, ALT of +# dbSNP b157's GRCh38 VCF, read from NCBI's Variation API on 2026-11-07), every +# GRCh38 copy of the site (1-based), the plus-strand base pair read there, and +# the GRCh38 sequence of 10 bases on each side of every copy (UCSC hg38 and +# Ensembl GRCh38 agree, 2026-11-07): the same 21 bases at each copy, which is +# why the reads cannot be placed on one copy. tests/test_traits_report.py +# checks every copy against the context and the record against the marker, and +# tests/e2e/genotype-rules-2-smn1-marker.sh compares the records with dbSNP. +# min_depth fewer reads over all copies together: too little depth +# min_alt reads with the marker's base for 'present' ... +# min_fraction ... and a share of the reads at least this +# max_absent this many reads with the marker's base or fewer: 'not seen' +# (one or two reads can be sequencing errors); between the two +# the marker is 'uncertain' +MARKERS = { + # KIV-2 4925G>A: the last base of the second KIV-2 exon (its donor splice + # site), G>A on LPA's minus strand, so C>T on the plus strand. GRCh38 holds + # it at the end of each of the six KIV-2 units (exons 5, 7, ..., 15 of + # ENST00000316300.10; NM_005577.4:c.733G>A at the first, p.Ala245Thr) and + # at the near-identical last exon of KIV-1 (exon 3, c.391): seven copies, + # the block bin/lpa_kiv2.py counts as N + 2 copies (N KIV-2 repeats and + # one KIV-1 per chromosome). dbSNP names the first KIV-2 copy only. + "LPA_KIV2_4925": { + "name": "LPA KIV-2 4925G>A", "chrom": "chr6", "ref": "C", "alt": "T", + "dbsnp": ("rs1484919544", "chr6", 160640667, "C", "T"), + "copies": (160646214, 160640667, 160635123, 160629574, 160624030, 160618484, 160612932), + "context": "GACTCCTTACCTTGTTCGGAA", + "min_depth": 100, "min_alt": 5, "min_fraction": 0.008, "max_absent": 2, + }, + # SMN1 g.27134T>G is NM_000344.4:c.*3+80T>G in intron 7 (plus strand). + # SMN2 holds the same 21 bases 875,420 bases upstream, as at exon 7. + "SMN1_27134": { + "name": "SMN1 g.27134T>G", "chrom": "chr5", "ref": "T", "alt": "G", + "dbsnp": ("rs143838139", "chr5", 70952074, "T", "G"), + "copies": (70952074, 70076654), + "context": "TTTTTGAACATTTAAAAAGTT", + "min_depth": 10, "min_alt": 3, "min_fraction": 0.05, "max_absent": 1, + }, +} +MARKER_STATUSES = ("present", "not seen", "uncertain", "too little depth", "not assessed") +# SMN1 exon 7 (c.840, chr5:70,951,946): where step 35's report reads the SMN1 psCN. +SMN1_EXON7 = ("chr5", 70951946) +SMN1_NOT_PROOF = "This marker does not prove 2+0, and its absence does not rule 2+0 out." + + +def marker_errors(name, m=None, dbsnp=None): + """Problems of one marker (empty when it is consistent): the dbSNP record + is one of its copies with its ref and alt, on the plus strand and not + swapped; every copy sits in the middle of the same context, whose middle + base is the ref.""" + m = MARKERS[name] if m is None else m + rs, chrom, pos, ref, alt = m["dbsnp"] if dbsnp is None else dbsnp + comp = {"A": "T", "C": "G", "G": "C", "T": "A"} + errs = [] + if chrom != m["chrom"] or pos not in m["copies"]: + errs.append(f"{name}: dbSNP places {rs} at {chrom}:{pos}, not at one of the marker's copies") + if (m["ref"], m["alt"]) == (alt, ref): + errs.append(f"{name}: ref and alt are swapped (marker {m['ref']}>{m['alt']}, dbSNP {ref}>{alt})") + elif (m["ref"], m["alt"]) == (comp.get(ref), comp.get(alt)): + errs.append(f"{name}: the alleles are the minus-strand bases (marker {m['ref']}>{m['alt']}, " + f"dbSNP {ref}>{alt}): write the plus-strand bases") + elif (m["ref"], m["alt"]) != (ref, alt): + errs.append(f"{name}: marker {m['ref']}>{m['alt']}, dbSNP {ref}>{alt}") + ctx = m["context"] + if len(ctx) % 2 != 1 or ctx[len(ctx) // 2] != m["ref"]: + errs.append(f"{name}: the middle base of the context is not the ref {m['ref']}") + if len(set(m["copies"])) != len(m["copies"]): + errs.append(f"{name}: a copy is listed twice") + return errs + + +def marker_regions(name): + """The copy sites as samtools regions, chrom:pos-pos, in the marker's order.""" + m = MARKERS[name] + return [f"{m['chrom']}:{p}-{p}" for p in m["copies"]] + + +def read_marker_pileup(path, name): + """(ref reads, marker reads, other reads) over every copy of the marker, + from `samtools mpileup` lines (chrom, pos, ref, depth, bases, ...) at its + sites, without a reference: a read that carries the marker's base counts + wherever the aligner put it. Lines at other positions are left out; a + copy with no line has no read. Raises ValueError on a line that is not + an mpileup line.""" + import copy_number # bin/: its mpileup base counter, shared with the CYP21A2 sites + m = MARKERS[name] + want = set(m["copies"]) + n = dict.fromkeys("ACGT", 0) + with open(path, errors="replace") as f: + for line in f: + if not line.strip(): + continue + c = line.rstrip("\n").split("\t") + if len(c) < 5 or not c[1].isdigit(): + raise ValueError(f"{path}: expected samtools mpileup lines (chrom, pos, ref, depth, bases, ...); " + f"got {line.strip()[:80]!r}") + if c[0] != m["chrom"] or int(c[1]) not in want: + continue + for k, v in copy_number.count_bases(c[4], c[2]).items(): + n[k] += v + ref, alt = n[m["ref"]], n[m["alt"]] + return ref, alt, sum(n.values()) - ref - alt + + +def assess_marker(name, counts): + """{status, alt_reads, depth, fraction} of a marker from read_marker_pileup's + counts (None: the pileup could not be read). depth counts the reads with + the ref or the marker's base.""" + m = MARKERS[name] + if counts is None: + return {"status": "not assessed", "alt_reads": ".", "depth": ".", "fraction": "."} + ref, alt, _ = counts + depth = ref + alt + out = {"alt_reads": str(alt), "depth": str(depth), "fraction": f"{alt / depth:.4f}" if depth else "."} + if depth < m["min_depth"]: + out["status"] = "too little depth" + elif alt >= m["min_alt"] and alt / depth >= m["min_fraction"]: + out["status"] = "present" + elif alt <= m["max_absent"]: + out["status"] = "not seen" + else: + out["status"] = "uncertain" + return out + + +def smn_copy_numbers(path): + """(SMN1 psCN at exon 7 or None, its text, agCN at the marker site or None) + from step 35's SMN table. A copy number counts only with filter PASS and + quality 20 or more, the threshold of step 35's reading.""" + def regions(r): + out = [(r["chrom"], int(r["start"]) + 1, int(r["end"]))] + for h in (r.get("homologous_regions") or "").split(","): + mm = re.match(r"^([^:]+):(\d+)-(\d+)", h.strip()) + if mm: + out.append((mm.group(1), int(mm.group(2)), int(mm.group(3)))) + return out + + def value(r, kind, idx): + vals = (r.get(kind) or "").split(",") + quals = (r.get(kind + "_qual") or "").split(",") + v = vals[idx] if idx < len(vals) else "." + q = quals[idx] if idx < len(quals) else "." + try: + ok = (r.get(kind + "_filter") == "PASS" and float(q) >= 20 and v.isdigit()) + except ValueError: + ok = False + return ((int(v) if ok else None), + f"{v} (quality {q}, filter {r.get(kind + '_filter') or '.'})" + ("" if ok else ", not resolved")) + + sm = MARKERS["SMN1_27134"] + ps, ps_text, ag = None, "not read (no SMN table from step 35)", None + if not path or not os.path.exists(path): + return ps, ps_text, ag + ps_text = "not read (no row of the SMN table covers SMN1 exon 7)" + with open(path, errors="replace") as f: + header = f.readline().rstrip("\n").split("\t") + for line in f: + r = dict(zip(header, line.rstrip("\n").split("\t"))) + try: + regs = regions(r) + except (KeyError, ValueError): + continue + for i, (c, s, e) in enumerate(regs): + if c == SMN1_EXON7[0] and s <= SMN1_EXON7[1] <= e: + ps, ps_text = value(r, "psCN", i) + # agCN is one value per row: the row whose region, or a + # homologous region, holds a copy of the marker site + if c == sm["chrom"] and any(s <= p <= e for p in sm["copies"]): + ag = value(r, "agCN", 0)[0] + return ps, ps_text, ag + + +def smn1_marker(counts, smn_table=None): + """[(metric, value)] of step 35's SMN1 marker table: the marker's reading + beside the SMN1 psCN. The message always says the marker does not prove + 2+0 (SMN1_NOT_PROOF).""" + m = MARKERS["SMN1_27134"] + a = assess_marker("SMN1_27134", counts) + ps, ps_text, ag = smn_copy_numbers(smn_table) + copies = "." + if a["status"] in ("present", "uncertain") and ag is not None and a["fraction"] != ".": + copies = f"{float(a['fraction']) * ag:.1f}" + st = a["status"] + if st == "not assessed": + head = f"{m['name']} ({m['dbsnp'][0]}): not assessed, the reads at the marker could not be read." + else: + head = (f"{m['name']} ({m['dbsnp'][0]}): {st}, {a['alt_reads']} of {a['depth']} reads carry G at the " + "SMN1 and SMN2 copies of the site together" + + (f" (about {copies} of the {ag} SMN1 and SMN2 copies)" if copies != "." else "") + + (f", under the {m['min_depth']} reads a reading needs" if st == "too little depth" else "") + + ". The reads cannot tell whether a G sits on SMN1 or SMN2.") + if ps == 2: + body = {"present": " With two SMN1 copies, this marker goes with a higher chance that both sit on one " + "chromosome (2+0, a silent carrier): in Luo 2014 it raised the residual carrier risk of a " + "two-copy result (in African Americans from about 1 in 121 to 1 in 34); the figures " + "depend on ancestry.", + "not seen": " With two SMN1 copies and no marker, the chance of 2+0 is lower than for a two-copy " + "result alone (in Luo 2014, in African Americans, about 1 in 396 against 1 in 121), " + "not zero; the figures depend on ancestry."}.get( + st, " With two SMN1 copies the marker would refine the chance of 2+0; here it gives no reading.") + elif ps is not None: + body = (f" SMN1 has {ps} cop{'y' if ps == 1 else 'ies'} here; the marker refines the carrier risk of a " + "two-copy result only.") + else: + body = f" SMN1 copy number: {ps_text}; the marker refines the carrier risk of a two-copy result only." + return [("marker", m["name"]), ("rsid", m["dbsnp"][0]), ("status", st), ("alt_reads", a["alt_reads"]), + ("depth", a["depth"]), ("fraction", a["fraction"]), ("copies", copies), + ("smn1_pscn", ps_text), ("agcn", "." if ag is None else str(ag)), + ("message", head + body + " " + SMN1_NOT_PROOF)] + + +def write_metric_table(path, rows): + tmp = path + ".tmp" + with open(tmp, "w", encoding="utf-8") as f: + f.write("metric\tvalue\n") + for k, v in rows: + f.write(f"{k}\t{v}\n") + os.replace(tmp, path) + + def write_tsv(path, out, columns=OUT_COLUMNS): tmp = path + ".tmp" with open(tmp, "w", encoding="utf-8") as f: @@ -524,7 +767,32 @@ def main(argv=None): ap.add_argument("--rules-out", help="the sample's genotype rules (TSV)") ap.add_argument("--phase-sites", metavar="VCF", help="write the sites WhatsHap must phase and exit") ap.add_argument("--check-table", metavar="TABLE", help="check the table's schema (and --rules') and exit") + ap.add_argument("--marker-regions", metavar="NAME", choices=sorted(MARKERS), + help="print a copy-aware marker's sites as samtools regions and exit") + ap.add_argument("--marker", choices=["SMN1_27134"], help="read the SMN1 marker from --pileup") + ap.add_argument("--pileup", help="samtools mpileup at the --marker-regions sites") + ap.add_argument("--copy-table", help="step 35's SMN table (the SMN1 psCN and the agCN)") + ap.add_argument("--marker-out", help="the marker's table (metricvalue)") a = ap.parse_args(argv) + if a.marker_regions: + print("\n".join(marker_regions(a.marker_regions))) + return 0 + if a.marker: + if not (a.pileup and a.marker_out): + ap.error("--marker needs --pileup and --marker-out") + try: + counts = read_marker_pileup(a.pileup, a.marker) + except (OSError, ValueError) as e: + print(f"WARNING: {e}; the marker is not assessed", file=sys.stderr) + counts = None + rows = smn1_marker(counts, a.copy_table) + try: + write_metric_table(a.marker_out, rows) + except OSError as e: + print(f"ERROR: cannot write: {e}", file=sys.stderr) + return 2 + print(dict(rows)["message"]) + return 0 try: if a.check_table: rows = read_table(a.check_table) diff --git a/tests/test_traits_report.py b/tests/test_traits_report.py index be38c68d..ca368b62 100644 --- a/tests/test_traits_report.py +++ b/tests/test_traits_report.py @@ -46,7 +46,20 @@ site and the mistake, on planted copies: ref and alt swapped, the minus-strand bases, an rsID the dbSNP table lacks, a wrong position; and on a rules file with no site, and on a dbSNP table of another build than - versions.env pins. + versions.env pins; + 11. the copy-aware markers (MARKERS: LPA KIV-2 4925G>A, SMN1 g.27134T>G): + each pinned dbSNP record is one of the marker's copies with its ref and + alt, and a planted swap, a minus-strand pair and a position off the + copies are each caught; the context's middle base is the ref; + --marker-regions prints one region per copy; on mpileup lines made up + for each case at the SMN1 and SMN2 copies: a quarter of the reads with G + is 'present' (copies = fraction x agCN), T only (the reference) 'not + seen', 2 reads with G 'uncertain', under 10 reads 'too little depth', an + unreadable pileup 'not assessed' (exit 0); the SMN1 psCN and the agCN + come from step 35's table (homologous region included, quality under 20 + not resolved); every message, for every status and every psCN, says the + marker does not prove 2+0 (control: a message without the sentence fails + the same test). Run: python3 tests/test_traits_report.py """ @@ -255,6 +268,107 @@ def line(rid, gts): shutil.rmtree(work) +SMN_HEADER = "chrom\tstart\tend\tlocus\tagCN_filter\tagCN\tagCN_qual\tpsCN_filter\tpsCN\tpsCN_qual\thomologous_regions\n" + + +def says_not_proof(message): + """The SMN1 marker's message must say the marker does not prove 2+0.""" + return "does not prove 2+0" in message + + +def marker_checks(): + """11. The copy-aware markers.""" + for name, m in tr.MARKERS.items(): + check(f"marker {name}: the pinned dbSNP record agrees with the marker", tr.marker_errors(name) == [], + tr.marker_errors(name)) + rs, chrom, pos, ref, alt = m["dbsnp"] + comp = {"A": "T", "C": "G", "G": "C", "T": "A"} + e = tr.marker_errors(name, dict(m, ref=alt, alt=ref)) + check(f"marker {name} control: ref and alt swapped is caught", any("swapped" in x for x in e), e) + e = tr.marker_errors(name, dict(m, ref=comp[ref], alt=comp[alt], + context=m["context"][:10] + comp[ref] + m["context"][11:])) + check(f"marker {name} control: the minus-strand bases are caught", any("minus-strand" in x for x in e), e) + e = tr.marker_errors(name, dbsnp=(rs, chrom, pos + 1, ref, alt)) + check(f"marker {name} control: a dbSNP position off the copies is caught", + any("not at one of the marker's copies" in x for x in e), e) + e = tr.marker_errors(name, dict(m, context=m["context"][:10] + alt + m["context"][11:])) + check(f"marker {name} control: a context whose middle base is not the ref is caught", + any("middle base" in x for x in e), e) + res = run("--marker-regions", name) + check(f"marker {name}: --marker-regions prints one region per copy", + res.returncode == 0 and res.stdout.split() == [f"{chrom}:{p}-{p}" for p in m["copies"]], res) + check("LPA_KIV2_4925: seven copies (six KIV-2 units and KIV-1)", len(tr.MARKERS["LPA_KIV2_4925"]["copies"]) == 7) + + sm = tr.MARKERS["SMN1_27134"] + work = tempfile.mkdtemp() + try: + def pileup(t_reads, g_reads, extra=""): + p = os.path.join(work, "p.tsv") + write(p, "".join(f"chr5\t{pos}\tN\t{t_reads + g_reads}\t^I{'T' * t_reads}{'g' * g_reads}$\t.\n" + for pos in sm["copies"]) + extra) + return p + + def smn(pscn="2,2", pqual="30,30", pfilt="PASS", agcn="4", homolog="chr5:70049000-70078000"): + p = os.path.join(work, "smn.tsv") + write(p, SMN_HEADER + f"chr5\t70925000\t70953000\tSMN1\tPASS\t{agcn}\t50\t{pfilt}\t{pscn}\t{pqual}\t{homolog}\n") + return p + + def table(p, smn_path): + out = os.path.join(work, "o.tsv") + res = run("--marker", "SMN1_27134", "--pileup", p, "--copy-table", smn_path, "--marker-out", out) + return res, dict(x.split("\t", 1) for x in open(out).read().splitlines()[1:]) if os.path.exists(out) else {} + + res, t = table(pileup(15, 5), smn()) + check("SMN1 marker: a quarter of the reads with G is present, exit 0", + res.returncode == 0 and t.get("status") == "present" and t.get("alt_reads") == "10" + and t.get("depth") == "40" and t.get("fraction") == "0.2500", (res, t)) + check("SMN1 marker: copies = fraction x agCN (0.25 x 4)", t.get("copies") == "1.0" and t.get("agcn") == "4", t) + check("SMN1 marker: the SMN1 psCN at exon 7 is read", t.get("smn1_pscn") == "2 (quality 30, filter PASS)", t) + check("SMN1 marker: present with two copies says it raises the chance of 2+0 and does not prove it", + "higher chance" in t.get("message", "") and says_not_proof(t.get("message", "")), t) + res, t = table(pileup(20, 0), smn()) + check("SMN1 marker: reads with T only (the reference) is not seen, copies '.'", + t.get("status") == "not seen" and t.get("copies") == "." and "lower" in t.get("message", ""), t) + res, t = table(pileup(20, 0, extra="chr5\t70952074\tN\t2\tGG\t.\n"), smn()) + check("SMN1 marker: 2 reads with G is uncertain", t.get("status") == "uncertain", t) + res, t = table(pileup(3, 1), smn()) + check("SMN1 marker: under 10 reads is too little depth", t.get("status") == "too little depth", t) + res, t = table(pileup(15, 5, extra="chr5\t70952075\tN\t30\t" + "G" * 30 + "\t.\n"), smn()) + check("SMN1 marker: reads at another position are not counted", t.get("alt_reads") == "10", t) + bad = os.path.join(work, "bad.tsv") + write(bad, "nonsense\n") + res, t = table(bad, smn()) + check("SMN1 marker: an unreadable pileup is not assessed, exit 0", + res.returncode == 0 and t.get("status") == "not assessed", (res, t)) + res, t = table(pileup(15, 5), smn(pscn="1,3")) + check("SMN1 marker: another SMN1 copy number says the marker refines a two-copy result only", + "SMN1 has 1 copy here" in t.get("message", ""), t) + res, t = table(pileup(15, 5), smn(pqual="12,30")) + check("SMN1 marker: a psCN of quality under 20 is not resolved", + t.get("smn1_pscn") == "2 (quality 12, filter PASS), not resolved" + and "SMN1 copy number: 2 (quality 12, filter PASS), not resolved;" in t.get("message", ""), t) + # the SMN table with SMN2's region first: exon 7 of SMN1 is the homologous region + p = os.path.join(work, "smn2.tsv") + write(p, SMN_HEADER + "chr5\t70049000\t70078000\tSMN1\tPASS\t4\t50\tPASS\t2,2\t30,40\tchr5:70925000-70953000\n") + res, t = table(pileup(15, 5), p) + check("SMN1 marker: exon 7 found through the homologous region (its own psCN and quality)", + t.get("smn1_pscn") == "2 (quality 40, filter PASS)" and t.get("agcn") == "4", t) + res, t = table(pileup(15, 5), os.path.join(work, "missing.tsv")) + check("SMN1 marker: no SMN table says so", "no SMN table" in t.get("smn1_pscn", ""), t) + + # Every status with every psCN: the message says the marker does not prove 2+0. + msgs = [] + for counts in [(30, 10, 0), (40, 0, 0), (38, 2, 0), (3, 1, 0), None]: + for path in (smn(), smn(pscn="1,3"), smn(pqual="5,5"), os.path.join(work, "missing.tsv")): + msgs.append(dict(tr.smn1_marker(counts, path))["message"]) + check(f"SMN1 marker: every message ({len(msgs)}) says it does not prove 2+0", + all(says_not_proof(x) for x in msgs), [x for x in msgs if not says_not_proof(x)]) + check("SMN1 marker control: a message without the sentence fails that test", + not says_not_proof(msgs[0].replace(tr.SMN1_NOT_PROOF, ""))) + finally: + shutil.rmtree(work) + + def main(): rules = tr.read_rules(RULES) dbsnp_checks(rules) @@ -480,6 +594,7 @@ def rule(rid, *lines, phased=None): shutil.rmtree(work) multi_snp_checks(rules) + marker_checks() print("\nRESULT:", "ALL PASS" if FAILS == 0 else f"{FAILS} FAILED") return 1 if FAILS else 0 From ff45aae04df7833c7c7f1f785182d353e4efd915 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 19:39:17 +0100 Subject: [PATCH 11/20] feat(lpa): print the KIV-2 4925G>A null allele beside the KIV-2 count Step 46 and the LPA processes pile up the reads at the seven GRCh38 copies of the 4925G>A site and add the kiv2_4925_* rows to the table, with one sentence after the message, so both reports print it under the count. Bead: pgp-pfw6.21 --- bin/lpa_kiv2.py | 104 +++++++++++++++++++++++++++++++++++++- modules/local/lpa/main.nf | 39 ++++++++------ scripts/46-lpa.sh | 46 ++++++++++++----- tests/test_lpa_kiv2.py | 89 +++++++++++++++++++++++++++++++- 4 files changed, 247 insertions(+), 31 deletions(-) diff --git a/bin/lpa_kiv2.py b/bin/lpa_kiv2.py index 7e1405db..36ef6144 100755 --- a/bin/lpa_kiv2.py +++ b/bin/lpa_kiv2.py @@ -26,6 +26,14 @@ over the same regions, and write LPA.tsv (metricvalue). Exit 0 whatever the status; 2 on unreadable input. Two empty tables mean the BAM has no chr6: status 'not assessed'. + lpa_kiv2.py null-allele --pileup PILEUP.tsv --table LPA.tsv + read `samtools mpileup` (no reference, every MAPQ) at the seven copies + of the KIV-2 4925G>A site (`traits_report.py --marker-regions + LPA_KIV2_4925`) and add its reading to LPA.tsv, written by `estimate` + or `copy_number.py call --locus LPA_KIV2`: the kiv2_4925_* rows, and + one sentence after the message, so both reports print it beside the + count. Exit 0 whatever the reading (an unreadable pileup reads 'not + assessed'); 2 when LPA.tsv cannot be read or written. The regions are exons of the MANE Select transcript ENST00000316300.10, numbered from its first exon: @@ -51,11 +59,24 @@ of an aligner or a library. status is 'estimated', 'too little depth' or 'not assessed'; the message says -why. Research estimate, not a diagnostic test: the actionable test is a blood +why. + +The KIV-2 4925G>A splice variant (Coassin et al. 2017, PMID 28444229) sits at +the donor splice site of the second KIV-2 exon: carriers had a much lower +Lp(a) than non-carriers with the same isoform size, most of all with small +isoforms, and it is usually in one or two of a person's KIV-2 repeats. GRCh38 +holds the site in each of its six units and in KIV-1, so null-allele adds up +the reads at all seven (bin/traits_report.py MARKERS) and reads the share +with the variant as a fraction of copies: kiv2_4925_status is 'present', +'not seen', 'uncertain', 'too little depth' or 'not assessed', and with a +count, kiv2_4925_copies is the fraction x (N + 2), the copies that carry it. +A short-read fraction cannot say which chromosome or which repeat carries +it, so it cannot say whether it sits on the small or the large isoform. Research estimate, not a diagnostic test: the actionable test is a blood Lp(a) measurement. Standard library only: it runs in the Python image and under tests/test_lpa_kiv2.py. """ import argparse +import os import sys import copy_number as cn @@ -141,6 +162,77 @@ def assess(m, q1): "mapq0_share", "message", "advice") +NULL_ORDER = ("kiv2_4925_status", "kiv2_4925_alt_reads", "kiv2_4925_depth", "kiv2_4925_fraction", + "kiv2_4925_copies", "kiv2_4925_message") +NULL_LIMIT = ("A short-read fraction cannot say which chromosome or which repeat carries it, so it cannot say " + "whether it sits on the small or the large isoform.") + + +def null_allele(counts, kiv2_copies="."): + """[(metric, value)] of the KIV-2 4925G>A reading from + traits_report.read_marker_pileup's counts (None: unreadable), with the + count N of the table ('.' when there is none).""" + import traits_report as tr # bin/: the copy-aware marker catalog + m = tr.MARKERS["LPA_KIV2_4925"] + a = tr.assess_marker("LPA_KIV2_4925", counts) + st = a["status"] + try: + total = float(kiv2_copies) + 2 + except ValueError: + total = None + copies = "." + if st in ("present", "uncertain") and total is not None: + copies = f"{float(a['fraction']) * total:.1f}" + name = f"{m['name']} (a splice variant that lowers the Lp(a) made from its chromosome, Coassin 2017)" + if st == "not assessed": + msg = f"{name}: not assessed, the reads at its seven GRCh38 copies could not be read." + else: + msg = (f"{name}: {st}, {a['alt_reads']} of {a['depth']} reads at its seven GRCh38 copies carry it" + + (f" (about {copies} of the {total:.0f} copies counted)" if copies != "." else "") + + (f", under the {m['min_depth']} reads a reading needs" if st == "too little depth" else "") + + ".") + if st == "present": + msg += (" Carriers had a much lower Lp(a) than non-carriers with the same number of repeats, most " + "of all with few repeats. " + NULL_LIMIT) + elif st == "not seen": + msg += (" Not seen in the reads is not proof of absence: one carrying repeat among many gives few " + "reads, and a short-read fraction never says which chromosome or repeat holds a variant.") + else: + msg += " " + NULL_LIMIT + return [("kiv2_4925_status", st), ("kiv2_4925_alt_reads", a["alt_reads"]), ("kiv2_4925_depth", a["depth"]), + ("kiv2_4925_fraction", a["fraction"]), ("kiv2_4925_copies", copies), ("kiv2_4925_message", msg)] + + +def add_null_allele(table_path, pileup_path): + """Add the 4925G>A rows to step 46's table, and the sentence after its message.""" + import traits_report as tr + with open(table_path, errors="replace") as f: + lines = [ln.rstrip("\n").split("\t", 1) for ln in f if ln.strip()] + if not lines or lines[0] != ["metric", "value"] or any(len(c) != 2 for c in lines): + raise ValueError(f"{table_path}: not a metricvalue table") + rows = [(k, v) for k, v in lines[1:] if not k.startswith("kiv2_4925_")] + kv = dict(rows) + try: + counts = tr.read_marker_pileup(pileup_path, "LPA_KIV2_4925") + except (OSError, ValueError) as e: + print(f"WARNING: {e}; KIV-2 4925G>A is not assessed", file=sys.stderr) + counts = None + new = null_allele(counts, kv.get("kiv2_copies", ".")) + msg = dict(new)["kiv2_4925_message"] + def joined(v): + # a sentence of an earlier run goes first, so a second run adds it once + v = v.split(tr.MARKERS["LPA_KIV2_4925"]["name"] + " (", 1)[0].strip() + return f"{v}{'' if not v or v.endswith('.') else '.'} {msg}".strip() + out = [(k, joined(v) if k == "message" else v) for k, v in rows] + new + tmp = table_path + ".tmp" + with open(tmp, "w") as f: + f.write("metric\tvalue\n") + for k, v in out: + f.write(f"{k}\t{v}\n") + os.replace(tmp, table_path) + return msg + + def main(argv=None): ap = argparse.ArgumentParser(prog="lpa_kiv2.py", description=__doc__.split("\n", 1)[0]) sub = ap.add_subparsers(dest="cmd", required=True) @@ -149,7 +241,17 @@ def main(argv=None): c.add_argument("--bedcov", required=True) c.add_argument("--bedcov-mapq1", required=True, help="samtools bedcov -j -Q 1 over the same regions") c.add_argument("--out", required=True) + n = sub.add_parser("null-allele", help="add the KIV-2 4925G>A reading to the table") + n.add_argument("--pileup", required=True, help="samtools mpileup at traits_report.py --marker-regions LPA_KIV2_4925") + n.add_argument("--table", required=True, help="the table estimate wrote") args = ap.parse_args(argv) + if args.cmd == "null-allele": + try: + print("LPA " + add_null_allele(args.table, args.pileup)) + except (OSError, ValueError) as e: + print(f"ERROR: {e}", file=sys.stderr) + return 2 + return 0 if args.cmd == "bed": print("\n".join(bed_lines())) return 0 diff --git a/modules/local/lpa/main.nf b/modules/local/lpa/main.nf index 3cd83c7e..24345b77 100644 --- a/modules/local/lpa/main.nf +++ b/modules/local/lpa/main.nf @@ -4,11 +4,15 @@ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Opt-in ('lpa' in --tools). Three processes, as in scripts/46-lpa.sh: LPA_REGIONS bin/copy_number.py bed --locus LPA_KIV2, once: the KIV-2 - and single-copy LPA exons (the LPA_KIV2 row of bin/cn_loci.tsv) - LPA_DEPTH samtools bedcov -j over them, all MAPQ and MAPQ >= 1, for - each sample with a BAM; empty tables when the BAM has no chr6 + and single-copy LPA exons (the LPA_KIV2 row of bin/cn_loci.tsv), + and the seven GRCh38 copies of the KIV-2 4925G>A site + (bin/traits_report.py --marker-regions LPA_KIV2_4925) + LPA_DEPTH samtools bedcov -j over the exons, all MAPQ and MAPQ >= 1, + and samtools mpileup at the 4925G>A sites, for each sample + with a BAM; empty tables when the BAM has no chr6 LPA bin/copy_number.py call --locus LPA_KIV2: the count, its - 95% interval and the message (laid out by bin/lpa_kiv2.py) + 95% interval and the message (laid out by bin/lpa_kiv2.py), + then bin/lpa_kiv2.py null-allele: the 4925G>A reading beside it Equivalent to: scripts/46-lpa.sh ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -18,8 +22,8 @@ process LPA_REGIONS { label 'process_single' output: - path "lpa_regions.bed", emit: bed - path "versions.yml", emit: versions + path "lpa_*", emit: bed // lpa_regions.bed and lpa_sites.txt + path "versions.yml", emit: versions when: task.ext.when == null || task.ext.when @@ -27,6 +31,7 @@ process LPA_REGIONS { script: """ copy_number.py bed --locus LPA_KIV2 > lpa_regions.bed + traits_report.py --marker-regions LPA_KIV2_4925 > lpa_sites.txt cat <<-END_VERSIONS > versions.yml "${task.process}": @@ -36,7 +41,7 @@ process LPA_REGIONS { stub: """ - touch lpa_regions.bed + touch lpa_regions.bed lpa_sites.txt cat <<-END_VERSIONS > versions.yml "${task.process}": @@ -53,25 +58,28 @@ process LPA_DEPTH { input: tuple val(meta), path(bam), path(bai) - path(regions) // LPA_REGIONS' BED + path(regions) // LPA_REGIONS' lpa_regions.bed and lpa_sites.txt, read by name output: - tuple val(meta), path("${meta.id}_lpa_kiv2_depth.tsv"), path("${meta.id}_lpa_kiv2_depth_q1.tsv"), emit: depth - path "versions.yml", emit: versions + tuple val(meta), path("${meta.id}_lpa_kiv2_depth.tsv"), path("${meta.id}_lpa_kiv2_depth_q1.tsv"), + path("${meta.id}_lpa_kiv2_4925_pileup.tsv"), emit: depth + path "versions.yml", emit: versions when: task.ext.when == null || task.ext.when script: - // bedcov stops on a contig the BAM has no @SQ line for: without chr6 both + // bedcov stops on a contig the BAM has no @SQ line for: without chr6 the // tables stay empty. The same commands as scripts/46-lpa.sh. """ if samtools view -H ${bam} | awk -F '\\t' '\$1 == "@SQ" {for (i = 2; i <= NF; i++) if (\$i == "SN:chr6") f = 1} END {exit !f}'; then - samtools bedcov -j ${regions} ${bam} > ${meta.id}_lpa_kiv2_depth.tsv - samtools bedcov -j -Q 1 ${regions} ${bam} > ${meta.id}_lpa_kiv2_depth_q1.tsv + samtools bedcov -j lpa_regions.bed ${bam} > ${meta.id}_lpa_kiv2_depth.tsv + samtools bedcov -j -Q 1 lpa_regions.bed ${bam} > ${meta.id}_lpa_kiv2_depth_q1.tsv + while read -r r; do samtools mpileup -A -B -q 0 -Q 13 -d 100000 -r "\$r" ${bam}; done < lpa_sites.txt > ${meta.id}_lpa_kiv2_4925_pileup.tsv else : > ${meta.id}_lpa_kiv2_depth.tsv : > ${meta.id}_lpa_kiv2_depth_q1.tsv + : > ${meta.id}_lpa_kiv2_4925_pileup.tsv fi cat <<-END_VERSIONS > versions.yml @@ -82,7 +90,7 @@ process LPA_DEPTH { stub: """ - touch ${meta.id}_lpa_kiv2_depth.tsv ${meta.id}_lpa_kiv2_depth_q1.tsv + touch ${meta.id}_lpa_kiv2_depth.tsv ${meta.id}_lpa_kiv2_depth_q1.tsv ${meta.id}_lpa_kiv2_4925_pileup.tsv cat <<-END_VERSIONS > versions.yml "${task.process}": @@ -98,7 +106,7 @@ process LPA { publishDir { "${params.outdir}/${meta.id}/lpa" }, mode: params.publish_dir_mode input: - tuple val(meta), path(depth), path(depth_q1) + tuple val(meta), path(depth), path(depth_q1), path(pileup) output: tuple val(meta), path("${meta.id}_lpa_kiv2.tsv"), emit: table @@ -110,6 +118,7 @@ process LPA { script: """ copy_number.py call --locus LPA_KIV2 --bedcov ${depth} --bedcov-mapq1 ${depth_q1} --table ${meta.id}_lpa_kiv2.tsv + lpa_kiv2.py null-allele --pileup ${pileup} --table ${meta.id}_lpa_kiv2.tsv cat <<-END_VERSIONS > versions.yml "${task.process}": diff --git a/scripts/46-lpa.sh b/scripts/46-lpa.sh index b0e42802..73ec9fe4 100755 --- a/scripts/46-lpa.sh +++ b/scripts/46-lpa.sh @@ -8,7 +8,12 @@ # the KIV-2 exons against the depth over single-copy LPA exons # (bin/copy_number.py on the LPA_KIV2 row of bin/cn_loci.tsv, laid out by # bin/lpa_kiv2.py; the same code as the LPA Nextflow processes), with a 95% -# interval. The two common risk variants, rs10455872 and rs3798220, are rows +# interval. Beside the count it reads the KIV-2 4925G>A splice variant +# (Coassin et al. 2017), which lowers the Lp(a) made from the chromosome that +# carries it: samtools mpileup at the site's seven GRCh38 copies (the six +# KIV-2 units and KIV-1), read as a fraction of copies by bin/lpa_kiv2.py +# null-allele (the kiv2_4925_* rows, and one sentence after the message). +# The two common risk variants, rs10455872 and rs3798220, are rows # of step 42's traits panel; the report's Lp(a) section reads both. A research # estimate, not a diagnostic test: the actionable test is a blood Lp(a) # measurement. @@ -21,6 +26,8 @@ # its 95% interval, the depths, message # lpa/_lpa_kiv2_depth.tsv samtools bedcov -j over the regions (all MAPQ) # lpa/_lpa_kiv2_depth_q1.tsv the same, reads of MAPQ 1 or more +# lpa/_lpa_kiv2_4925_pileup.tsv samtools mpileup at the seven +# copies of the 4925G>A site (every MAPQ) set -euo pipefail SAMPLE=${1:?Usage: $0 } @@ -36,39 +43,50 @@ OUTDIR="${S}/lpa" TABLE="${OUTDIR}/${SAMPLE}_lpa_kiv2.tsv" DEPTH="${OUTDIR}/${SAMPLE}_lpa_kiv2_depth.tsv" DEPTH_Q1="${OUTDIR}/${SAMPLE}_lpa_kiv2_depth_q1.tsv" +PILEUP="${OUTDIR}/${SAMPLE}_lpa_kiv2_4925_pileup.tsv" BED="${OUTDIR}/regions.bed" +SITES="${OUTDIR}/sites.txt" echo "=== Lp(a): LPA KIV-2 repeat count from read depth: ${SAMPLE} ===" [ -f "$BAM" ] || { echo "ERROR: File not found: ${BAM}" >&2; exit 1; } [ -f "${BAM}.bai" ] || { echo "ERROR: File not found: ${BAM}.bai" >&2; exit 1; } echo "BAM: ${BAM}" mkdir -p "$OUTDIR" -rm -f "$TABLE" "$DEPTH" "$DEPTH_Q1" "$BED" +rm -f "$TABLE" "$DEPTH" "$DEPTH_Q1" "$PILEUP" "$BED" "$SITES" -echo "[1/3] The regions (bin/copy_number.py bed --locus LPA_KIV2)..." -# shellcheck disable=SC2016 # $1 belongs to the inner shell +echo "[1/4] The regions (bin/copy_number.py bed --locus LPA_KIV2) and the 4925G>A sites..." +# shellcheck disable=SC2016 # $1 and $2 belong to the inner shell run_in -v "${PGP_ROOT}/bin:/pgp-bin:ro" --cpus 1 --memory 1g "$PYTHON_IMAGE" \ - sh -c 'python3 /pgp-bin/copy_number.py bed --locus LPA_KIV2 > "$1"' _ "$(cpath "$BED")" + sh -c 'python3 /pgp-bin/copy_number.py bed --locus LPA_KIV2 > "$1" && + python3 /pgp-bin/traits_report.py --marker-regions LPA_KIV2_4925 > "$2"' _ "$(cpath "$BED")" "$(cpath "$SITES")" -echo "[2/3] Read depth over the KIV-2 and single-copy LPA exons (samtools bedcov -j, all MAPQ and MAPQ >= 1)..." -# bedcov stops on a contig the BAM has no @SQ line for: without chr6 both -# tables stay empty and the count is not assessed. The LPA_DEPTH process runs -# the same commands. -# shellcheck disable=SC2016 # $1 to $5 belong to the inner bash, the program to awk +echo "[2/4] Read depth over the KIV-2 and single-copy LPA exons (samtools bedcov -j, all MAPQ and MAPQ >= 1)" +echo " and the reads at the 4925G>A sites (samtools mpileup, every MAPQ)..." +# bedcov stops on a contig the BAM has no @SQ line for: without chr6 the +# tables stay empty and the count is not assessed. mpileup runs once per site +# (-r), every MAPQ, duplicates left out, bases of quality 13 or more. The +# LPA_DEPTH process runs the same commands. +# shellcheck disable=SC2016 # $1 to $7 belong to the inner bash, the program to awk run_in --cpus 1 --memory 2g "$SAMTOOLS_IMAGE" bash -euo pipefail -c ' if samtools view -H "$2" | awk -F "\t" "$5"; then samtools bedcov -j "$1" "$2" > "$3" samtools bedcov -j -Q 1 "$1" "$2" > "$4" - else : > "$3"; : > "$4"; fi' \ + while read -r r; do samtools mpileup -A -B -q 0 -Q 13 -d 100000 -r "$r" "$2"; done < "$7" > "$6" + else : > "$3"; : > "$4"; : > "$6"; fi' \ _ "$(cpath "$BED")" "$(cpath "$BAM")" "$(cpath "$DEPTH")" "$(cpath "$DEPTH_Q1")" \ - '$1 == "@SQ" {for (i = 2; i <= NF; i++) if ($i == "SN:chr6") f = 1} END {exit !f}' -rm -f "$BED" + '$1 == "@SQ" {for (i = 2; i <= NF; i++) if ($i == "SN:chr6") f = 1} END {exit !f}' \ + "$(cpath "$PILEUP")" "$(cpath "$SITES")" +rm -f "$BED" "$SITES" -echo "[3/3] The count and its interval (bin/copy_number.py call --locus LPA_KIV2)..." +echo "[3/4] The count and its interval (bin/copy_number.py call --locus LPA_KIV2)..." run_in -v "${PGP_ROOT}/bin:/pgp-bin:ro" --cpus 1 --memory 1g "$PYTHON_IMAGE" \ python3 /pgp-bin/copy_number.py call --locus LPA_KIV2 --bedcov "$(cpath "$DEPTH")" \ --bedcov-mapq1 "$(cpath "$DEPTH_Q1")" --table "$(cpath "$TABLE")" +echo "[4/4] KIV-2 4925G>A beside the count (bin/lpa_kiv2.py null-allele)..." +run_in -v "${PGP_ROOT}/bin:/pgp-bin:ro" --cpus 1 --memory 1g "$PYTHON_IMAGE" \ + python3 /pgp-bin/lpa_kiv2.py null-allele --pileup "$(cpath "$PILEUP")" --table "$(cpath "$TABLE")" + [ -s "$TABLE" ] || { echo "ERROR: no table written (${TABLE})" >&2; exit 1; } echo "" echo "Table: ${TABLE}" diff --git a/tests/test_lpa_kiv2.py b/tests/test_lpa_kiv2.py index be06694d..bcc764d0 100644 --- a/tests/test_lpa_kiv2.py +++ b/tests/test_lpa_kiv2.py @@ -28,12 +28,22 @@ reasons and the advice; a typical count with no risk allele is not flagged; without step 42's table the variants are 'not read'; 'too little depth' is shown as such; no card and no text section without - step 46's table (only the line in the list of steps that did not run). + step 46's table (only the line in the list of steps that did not run); + 8. KIV-2 4925G>A (null-allele) on mpileup lines made up for each case, at + the seven GRCh38 copies of the site: 5% of the reads with T is + 'present' with copies = fraction x (N + 2); reads with C only (the + reference) is 'not seen'; 3 reads with T is 'uncertain'; under 100 reads + is 'too little depth'; an unreadable pileup is 'not assessed', exit 0; + reads at other positions and read starts, ends and indels in the bases + column are not counted; the rows go after the table and the sentence + after the message, once, and the text report prints it beside the count; + a table that is not metricvalue exits 2. Run: python3 tests/test_lpa_kiv2.py """ import math import os +import shutil import subprocess import sys import tempfile @@ -241,6 +251,83 @@ def sample(t, name, n_copies, flank=30.0, genotypes=None): check("report: no card and no text section without step 46", "

Lp(a)" not in render_report.html_report(nothing) and "## Lp(a)" not in render_report.text_report(nothing)) + # 8. KIV-2 4925G>A beside the count + sites = traits_report.MARKERS["LPA_KIV2_4925"]["copies"] + + def pileup(path, per_site, extra=""): + """mpileup lines at the seven copies; per_site is (C reads, T reads) at each.""" + with open(path, "w") as f: + for p in sites: + c, t_ = per_site + bases = "^]" + "C" * c + "t" * t_ + "$" + f.write(f"chr6\t{p}\tN\t{c + t_}\t{bases}\t{'I' * (c + t_)}\n") + f.write(extra) + + def table_after(t, name, per_site, n_copies=38, extra=""): + d = os.path.join(t, name) + os.makedirs(d) + a, b, o, p = (os.path.join(d, x) for x in ("a.tsv", "b.tsv", "S_lpa_kiv2.tsv", "p.tsv")) + bedcov(a, depths(n_copies)) + bedcov(b, q1(depths(n_copies))) + lpa.main(["estimate", "--bedcov", a, "--bedcov-mapq1", b, "--out", o]) + pileup(p, per_site, extra) + rc = subprocess.run([sys.executable, SCRIPT, "null-allele", "--pileup", p, "--table", o], + capture_output=True, text=True) + return rc, dict(x.split("\t", 1) for x in open(o).read().splitlines()[1:]), o, p + + with tempfile.TemporaryDirectory() as t: + rc, tab, o, p = table_after(t, "present", (95, 5)) + check("4925G>A: 5% of the reads with T at the seven copies is present, exit 0", + rc.returncode == 0 and tab.get("kiv2_4925_status") == "present" + and tab.get("kiv2_4925_alt_reads") == "35" and tab.get("kiv2_4925_depth") == "700", (rc, tab)) + check("4925G>A: copies = fraction x (N + 2) = 0.05 x 40", + tab.get("kiv2_4925_fraction") == "0.0500" and tab.get("kiv2_4925_copies") == "2.0", tab) + check("4925G>A: the sentence follows the count's message, once", + tab["message"].startswith("about 38 KIV-2 repeats") and tab["message"].count("4925G>A") == 1 + and tab["message"].endswith(tab["kiv2_4925_message"]), tab["message"]) + check("4925G>A: the message says what a short-read fraction cannot say", + "cannot say which chromosome or which repeat" in tab["kiv2_4925_message"], tab) + before = open(o).read() + subprocess.run([sys.executable, SCRIPT, "null-allele", "--pileup", p, "--table", o], capture_output=True) + check("4925G>A: a second run leaves the table as it was (the rows and the sentence once)", + open(o).read() == before, open(o).read()) + rc, tab, *_ = table_after(t, "ref", (100, 0)) + check("4925G>A: reads with the reference base only is not seen", + tab.get("kiv2_4925_status") == "not seen" and tab.get("kiv2_4925_copies") == "." + and "not proof of absence" in tab["kiv2_4925_message"], tab) + rc, tab, *_ = table_after(t, "unc", (100, 0), extra="chr6\t160640667\tN\t3\tTTT\tIII\n") + check("4925G>A: 3 reads with T is uncertain", tab.get("kiv2_4925_status") == "uncertain" + and tab.get("kiv2_4925_alt_reads") == "3", tab) + rc, tab, *_ = table_after(t, "thin", (10, 1)) + check("4925G>A: under 100 reads is too little depth", tab.get("kiv2_4925_status") == "too little depth", tab) + rc, tab, *_ = table_after(t, "other", (100, 0), + extra="chr6\t160640668\tN\t50\t" + "T" * 50 + "\t.\nchr7\t160640667\tN\t50\t" + + "T" * 50 + "\t.\nchr6\t160635123\tN\t2\tC+2TTC-1T\tII\n") + check("4925G>A: reads at other positions and indel bases are not counted", + tab.get("kiv2_4925_status") == "not seen" and tab.get("kiv2_4925_alt_reads") == "0" + and tab.get("kiv2_4925_depth") == "702", tab) + rc, tab, o, p = table_after(t, "bad", (95, 5)) + with open(p, "w") as f: + f.write("not a pileup\n") + rc = subprocess.run([sys.executable, SCRIPT, "null-allele", "--pileup", p, "--table", o], + capture_output=True, text=True) + tab = dict(x.split("\t", 1) for x in open(o).read().splitlines()[1:]) + check("4925G>A: an unreadable pileup is not assessed, exit 0", + rc.returncode == 0 and tab.get("kiv2_4925_status") == "not assessed", (rc, tab)) + with open(o, "w") as f: + f.write("a\tb\tc\n") + rc = subprocess.run([sys.executable, SCRIPT, "null-allele", "--pileup", p, "--table", o], + capture_output=True, text=True) + check("4925G>A: a table that is not metricvalue exits 2", rc.returncode == 2, rc) + d = os.path.join(t, "rep", "S") + os.makedirs(os.path.join(d, "lpa")) + _, _, o, _ = table_after(t, "rep2", (95, 5)) + shutil.copy(o, os.path.join(d, "lpa", "S_lpa_kiv2.tsv")) + txt = render_report.text_report(collect_summary.collect("S", d)) + check("4925G>A: the text report prints it beside the count", + "KIV-2 repeats (both chromosomes): 38.0" in txt and "4925G>A" in txt + and txt.index("4925G>A") > txt.index("KIV-2 repeats (both chromosomes)"), txt[-2000:]) + print("\nRESULT:", "ALL PASS" if FAILS == 0 else f"{FAILS} FAILED") return 1 if FAILS else 0 From c925be1778e34b61171eff3a7354d1639f17e27f Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 19:39:23 +0100 Subject: [PATCH 12/20] feat(paralogs): read the SMN1 g.27134T>G silent-carrier marker beside the SMN1 psCN Step 35 and the PARASCOPY module write _smn1_marker.tsv from the reads at the marker's SMN1 and SMN2 copies and the SMN table. Bead: pgp-pfw6.21 --- modules/local/parascopy/main.nf | 24 ++++++++++++++++++ scripts/35-paralogs.sh | 45 ++++++++++++++++++++++++++++++--- 2 files changed, 66 insertions(+), 3 deletions(-) diff --git a/modules/local/parascopy/main.nf b/modules/local/parascopy/main.nf index 8d39f0ab..ce481860 100644 --- a/modules/local/parascopy/main.nf +++ b/modules/local/parascopy/main.nf @@ -44,6 +44,12 @@ 30 kb chimeric deletion and a possible conversion show (_cyp21a2_psv.tsv, docs/cn-loci.md). + Beside the SMN1 psCN, SMN1 g.27134T>G (rs143838139, Luo et al. 2014), a + marker of SMN1 duplication alleles: samtools mpileup at its SMN1 and SMN2 + copies, read as a fraction of copies by bin/traits_report.py --marker + SMN1_27134 with the SMN table (_smn1_marker.tsv). It never proves a + silent 2+0 carrier, and its message says so. + The interval rows of that catalog (ClinGen's recurrent CNV regions, the carrier deletion alleles and the Y chromosome AZF regions) take one more bedcov pass over all their segments, one line per row in @@ -83,6 +89,7 @@ process PARASCOPY { tuple val(meta), path("${meta.id}_smn_copy_number.tsv"), emit: copy_number, optional: true // the SMN1 row's tuple val(meta), path("${meta.id}_parascopy"), emit: results tuple val(meta), path("${meta.id}_hba_copy_number.tsv"), emit: hba + tuple val(meta), path("${meta.id}_smn1_marker.tsv"), emit: smn1_marker // SMN1 g.27134T>G beside the psCN tuple val(meta), path("${meta.id}_*_copy_number.tsv"), emit: tables // every catalog row's table, and HBA's tuple val(meta), path("${meta.id}_cn_dosage.tsv"), emit: dosage // HP, LCE3BC, CCL3L1 from depth tuple val(meta), path("${meta.id}_cyp21a2_psv.tsv"), emit: cyp21a2 // CYP21A2 from depth and PSVs @@ -149,6 +156,22 @@ process PARASCOPY { rm -rf depth # end of the paralog loci + # SMN1 g.27134T>G beside the SMN1 psCN (as scripts/35-paralogs.sh): mpileup at its SMN1 and SMN2 + # copies, then the reading with the SMN table; not assessed when it fails + if traits_report.py --marker-regions SMN1_27134 > smn1_marker_sites.txt \ + && while read -r r; do samtools mpileup -A -B -q 0 -Q 13 -d 100000 -r "\$r" ${bam}; done \ + < smn1_marker_sites.txt > smn1_marker_pileup.tsv \ + && traits_report.py --marker SMN1_27134 --pileup smn1_marker_pileup.tsv \ + --copy-table ${meta.id}_smn_copy_number.tsv --marker-out ${meta.id}_smn1_marker.tsv.tmp; then + mv ${meta.id}_smn1_marker.tsv.tmp ${meta.id}_smn1_marker.tsv + else + rm -f ${meta.id}_smn1_marker.tsv.tmp + echo "WARNING: the SMN1 g.27134T>G check failed; ${meta.id}_smn1_marker.tsv says not assessed." >&2 + printf 'metric\\tvalue\\nmarker\\tSMN1 g.27134T>G\\nrsid\\trs143838139\\nstatus\\tnot assessed\\nalt_reads\\t.\\ndepth\\t.\\nfraction\\t.\\ncopies\\t.\\nsmn1_pscn\\t.\\nagcn\\t.\\nmessage\\t%s\\n' \ + 'SMN1 g.27134T>G (rs143838139): not assessed, the check failed; see the log of step 35. This marker does not prove 2+0, and its absence does not rule 2+0 out.' > ${meta.id}_smn1_marker.tsv + fi + rm -f smn1_marker_sites.txt smn1_marker_pileup.tsv + # HBA1/HBA2 from depth (as scripts/35-paralogs.sh): all MAPQ and MAPQ 1 or more, then the verdict if copy_number.py bed --locus HBA > hba_regions.bed \\ && samtools bedcov hba_regions.bed ${bam} > hba_bedcov.tsv \\ @@ -231,6 +254,7 @@ process PARASCOPY { printf 'metric\\tvalue\\nstatus\\tnot assessed\\nalpha_globin_genes\\t.\\nestimate\\t.\\nhba2_depth\\t0\\nhba1_depth\\t0\\nflank_left_depth\\t0\\nflank_right_depth\\t0\\nmapq0_share\\t0\\nmessage\\tstub\\n' \\ > ${meta.id}_hba_copy_number.tsv printf 'locus\\tstatus\\tcopies\\testimate\\treason\\treading\\n' > ${meta.id}_cn_dosage.tsv + printf 'metric\\tvalue\\nstatus\\tnot assessed\\nmessage\\tstub\\n' > ${meta.id}_smn1_marker.tsv printf 'metric\\tvalue\\nstatus\\tnot assessed\\nreason\\tstub\\n' > ${meta.id}_cyp21a2_psv.tsv printf 'locus\\tname\\tchrom\\tstart\\tend\\tstatus\\tcopies\\texpected\\tchange\\testimate\\treason\\n' > ${meta.id}_cn_intervals.tsv diff --git a/scripts/35-paralogs.sh b/scripts/35-paralogs.sh index 350a8571..e984661b 100755 --- a/scripts/35-paralogs.sh +++ b/scripts/35-paralogs.sh @@ -36,6 +36,14 @@ # chimeric deletion and a possible conversion show; ${SAMPLE}_cyp21a2_psv.tsv # (docs/cn-loci.md). # +# Beside the SMN1 psCN the step reads SMN1 g.27134T>G (rs143838139; Luo et +# al. 2014), a marker of SMN1 duplication alleles: samtools mpileup at its +# SMN1 and SMN2 copies (the reads cannot tell them apart there), read as a +# fraction of copies by bin/traits_report.py --marker SMN1_27134. With two +# SMN1 copies it raises or lowers the chance of a silent 2+0 carrier; it +# never proves 2+0, and the table's message says so; +# ${SAMPLE}_smn1_marker.tsv. +# # The interval rows of that catalog (ClinGen's recurrent CNV regions such as # 22q11.2 and 1q21.1, the carrier deletion alleles and the Y chromosome AZF # regions) come from one more bedcov pass over all their segments, one line @@ -58,6 +66,7 @@ # # Output: ${GENOME_DIR}/${SAMPLE}/paralogs/${SAMPLE}__copy_number.tsv # (one per catalog row; the SMN1 row's is ${SAMPLE}_smn_copy_number.tsv) +# ${GENOME_DIR}/${SAMPLE}/paralogs/${SAMPLE}_smn1_marker.tsv # ${GENOME_DIR}/${SAMPLE}/paralogs/${SAMPLE}_hba_copy_number.tsv # ${GENOME_DIR}/${SAMPLE}/paralogs/${SAMPLE}_cn_dosage.tsv # ${GENOME_DIR}/${SAMPLE}/paralogs/${SAMPLE}_cyp21a2_psv.tsv @@ -158,9 +167,9 @@ mkdir -p "$OUTDIR" WORK="${OUTDIR}/${SAMPLE}_parascopy.part" rm -rf "$WORK" "${OUTDIR}/depth.part" # The depth tables of an earlier run go first: a step that stops after the -# new SMN tables are in place must not leave an old HBA, dosage, CYP21A2 or -# interval table that reads as this run's result. -rm -f "${OUTDIR}/${SAMPLE}_hba_copy_number.tsv" "${OUTDIR}/${SAMPLE}_cn_dosage.tsv" \ +# new SMN tables are in place must not leave an old SMN1 marker, HBA, dosage, +# CYP21A2 or interval table that reads as this run's result. +rm -f "${OUTDIR}/${SAMPLE}_smn1_marker.tsv" "${OUTDIR}/${SAMPLE}_hba_copy_number.tsv" "${OUTDIR}/${SAMPLE}_cn_dosage.tsv" \ "${OUTDIR}/${SAMPLE}_cyp21a2_psv.tsv" "${OUTDIR}/${SAMPLE}_cn_intervals.tsv" mkdir -p "$WORK" # Background depth once, then the copy number of each catalog row @@ -234,6 +243,34 @@ for t in "${TABLES[@]}"; do mv "${WORK}/${t}" "${OUTDIR}/${t}"; done rm -rf "${OUTDIR}/depth.part" "${OUTDIR}/${SAMPLE}_parascopy" mv "$WORK" "${OUTDIR}/${SAMPLE}_parascopy" +# SMN1 g.27134T>G beside the SMN1 psCN (as the PARASCOPY module): samtools +# mpileup at its SMN1 and SMN2 copies, one -r per site, every MAPQ, +# duplicates left out, bases of quality 13 or more; then +# bin/traits_report.py --marker SMN1_27134 with the SMN table (the SMN1 psCN +# at exon 7 and the agCN). A failure writes a table that says not assessed, +# and the step goes on. +SMN1M_OUT="${OUTDIR}/${SAMPLE}_smn1_marker.tsv" +SMN1M_WORK="${OUTDIR}/smn1_marker.part" +rm -rf "$SMN1M_OUT" "$SMN1M_WORK" +mkdir -p "$SMN1M_WORK" +# shellcheck disable=SC2016 # $1 to $3 belong to the inner bash +if run_in -v "${PGP_ROOT}/bin:/pgp-bin:ro" "${PARASCOPY_IMAGE}" \ + python3 /pgp-bin/traits_report.py --marker-regions SMN1_27134 > "${SMN1M_WORK}/sites.txt" \ + && run_in --cpus 1 --memory 2g "${PARASCOPY_IMAGE}" bash -euo pipefail -c ' + while read -r r; do samtools mpileup -A -B -q 0 -Q 13 -d 100000 -r "$r" "$1"; done < "$2" > "$3"' \ + _ "$(cpath "$BAM")" "$(cpath "${SMN1M_WORK}/sites.txt")" "$(cpath "${SMN1M_WORK}/pileup.tsv")" \ + && run_in -v "${PGP_ROOT}/bin:/pgp-bin:ro" "${PARASCOPY_IMAGE}" \ + python3 /pgp-bin/traits_report.py --marker SMN1_27134 --pileup "$(cpath "${SMN1M_WORK}/pileup.tsv")" \ + --copy-table "$(cpath "${OUTDIR}/${SAMPLE}_smn_copy_number.tsv")" --marker-out "$(cpath "${SMN1M_OUT}.tmp")"; then + mv "${SMN1M_OUT}.tmp" "$SMN1M_OUT" +else + rm -f "${SMN1M_OUT}.tmp" + echo "WARNING: the SMN1 g.27134T>G check failed; ${SMN1M_OUT##*/} says not assessed." >&2 + printf 'metric\tvalue\nmarker\tSMN1 g.27134T>G\nrsid\trs143838139\nstatus\tnot assessed\nalt_reads\t.\ndepth\t.\nfraction\t.\ncopies\t.\nsmn1_pscn\t.\nagcn\t.\nmessage\t%s\n' \ + 'SMN1 g.27134T>G (rs143838139): not assessed, the check failed; see the log of step 35. This marker does not prove 2+0, and its absence does not rule 2+0 out.' > "$SMN1M_OUT" +fi +rm -rf "$SMN1M_WORK" + # HBA1/HBA2 from depth (as the PARASCOPY module): samtools bedcov over the # regions of the HBA row of bin/cn_loci.tsv, all MAPQ and MAPQ 1 or more (a # BAM with no MAPQ 0 read over the genes left out the reads that fit both), @@ -365,6 +402,8 @@ for t in "${TABLES[@]}"; do echo "Copy number: ${OUTDIR}/${t}" column -t -s $'\t' "${OUTDIR}/${t}" 2>/dev/null || cat "${OUTDIR}/${t}" done +echo "SMN1 marker: ${SMN1M_OUT}" +awk -F'\t' '$1 == "message" {print " " $2}' "$SMN1M_OUT" echo "" echo "A row that says not assessed names the reason. SMN1 sits at chr5:70.92-70.95 Mb and SMN2 at chr5:70.05-70.08 Mb (GRCh38). psCN lists the copy" echo "number of the row's region first, then of each region in homologous_regions. A quality under 20" From 84a6072ed3f5db83bd7f3e3c86f7e9ef880e22a3 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 19:39:26 +0100 Subject: [PATCH 13/20] test(e2e): plant and clear the LPA and SMN1 markers in HG002 reads; docs and contract rows Bead: pgp-pfw6.21 --- assets/output_contract.tsv | 2 + docs/35-paralogs.md | 18 +++ docs/42-traits.md | 5 +- docs/46-lpa.md | 36 ++++-- scripts/ci/e2e-case-paths.tsv | 2 + scripts/ci/e2e-case-times.tsv | 2 + tests/e2e/genotype-rules-2-lpa-null.sh | 134 ++++++++++++++++++++ tests/e2e/genotype-rules-2-smn1-marker.sh | 147 ++++++++++++++++++++++ 8 files changed, 336 insertions(+), 10 deletions(-) create mode 100755 tests/e2e/genotype-rules-2-lpa-null.sh create mode 100755 tests/e2e/genotype-rules-2-smn1-marker.sh diff --git a/assets/output_contract.tsv b/assets/output_contract.tsv index 5477c732..7fec14eb 100644 --- a/assets/output_contract.tsv +++ b/assets/output_contract.tsv @@ -136,6 +136,7 @@ file {sample}/kir/database_release.txt 08 KIR_TYPING nextflow text experimental file {sample}/kir/{sample}_kir_genotype.tsv 08 KIR_TYPING nextflow tsv experimental - - file {sample}/logs/run_status.tsv - - nextflow tsv stable - - file {sample}/lpa/{sample}_lpa_kiv2.tsv 46 LPA nextflow tsv experimental - - +file {sample}/lpa/{sample}_lpa_kiv2_4925_pileup.tsv 46 LPA_DEPTH both tsv experimental - - file {sample}/lpa/{sample}_lpa_kiv2_depth.tsv 46 LPA_DEPTH nextflow tsv experimental - - file {sample}/lpa/{sample}_lpa_kiv2_depth_q1.tsv 46 LPA_DEPTH nextflow tsv experimental - - file {sample}/manta/results/variants/candidateSV.vcf.gz 04 MANTA nextflow vcf stable - - @@ -181,6 +182,7 @@ file {sample}/paralogs/{sample}_paraphase_copy_number.tsv 35b PARAPHASE_REPORT n file {sample}/paralogs/{sample}_paraphase_variants.tsv 35b PARAPHASE_REPORT nextflow tsv experimental - - file {sample}/paralogs/{sample}_parascopy 35 PARASCOPY nextflow dir experimental - - file {sample}/paralogs/{sample}_smn_copy_number.tsv 35 PARASCOPY nextflow tsv experimental - - +file {sample}/paralogs/{sample}_smn1_marker.tsv 35 PARASCOPY both tsv experimental - metric,value file {sample}/pgx_consensus/{sample}_outside_calls.tsv 36 PGX_CONSENSUS nextflow tsv stable - - file {sample}/pgx_consensus/{sample}_pgx_consensus.tsv 36 PGX_CONSENSUS nextflow tsv stable - - file {sample}/pharmcat/{sample}.match.json 07 PHARMCAT nextflow json stable - - diff --git a/docs/35-paralogs.md b/docs/35-paralogs.md index 5262c197..50dbb113 100644 --- a/docs/35-paralogs.md +++ b/docs/35-paralogs.md @@ -66,6 +66,7 @@ With Nextflow: `--tools ...,parascopy --parascopy_data ${GENOME_DIR}/reference/p 5. The gene-dosage rows (HP, LCE3BC, CCL3L1) the same way, into `${SAMPLE}_cn_dosage.tsv` ([the gene-dosage rows](cn-loci.md#the-gene-dosage-rows)). 6. `samtools bedcov` over the CYP21A2 row's segments (one pass, every read), `samtools mpileup` without a reference at its PSV sites, then `bin/copy_number.py call --locus CYP21A2 --pileup ...` writes `${SAMPLE}_cyp21a2_psv.tsv`. A failure prints a warning and writes a table that says not assessed (`check_failed`). 7. The interval rows into `${SAMPLE}_cn_intervals.tsv` ([the interval rows](cn-loci.md#the-interval-rows)). +8. SMN1 g.27134T>G beside the SMN1 psCN (before the HBA count): `samtools mpileup` without a reference at its SMN1 and SMN2 positions, then `bin/traits_report.py --marker SMN1_27134` with the SMN table writes `${SAMPLE}_smn1_marker.tsv` ([the SMN1 silent-carrier marker](#the-smn1-silent-carrier-marker)). A failure prints a warning and writes a table that says not assessed. ## Output @@ -75,11 +76,28 @@ All output is written to `${GENOME_DIR}/${SAMPLE}/paralogs/`. |---|---| | `${SAMPLE}_smn_copy_number.tsv` | One row per region: `chrom`, `start`, `end`, `locus`, `agCN_filter`, `agCN`, `agCN_qual`, `psCN_filter`, `psCN`, `psCN_qual`, `homologous_regions`. A locus that was not assessed has one row over the catalog region (`chrom`, `start`, `end`), with `locus` (the model name), `agCN_filter` and `psCN_filter` (both `not assessed: `) and `.` in every other column | | `${SAMPLE}_
_copy_number.tsv` | The same, for every other catalog row (`gba1`, `pms2_pms2cl`, `c4` and the rest of the table above) | +| `${SAMPLE}_smn1_marker.tsv` | `metric`, `value`: `marker`, `rsid`, `status` (`present`, `not seen`, `uncertain`, `too little depth`, `not assessed`), `alt_reads`, `depth`, `fraction`, `copies` (the fraction times agCN), `smn1_pscn` (the SMN1 psCN at exon 7 with its quality and filter), `agcn`, `message`. The step's log prints the message after the SMN table | | `${SAMPLE}_hba_copy_number.tsv` | `metric`, `value`: `status` (`called` or `not assessed`), `alpha_globin_genes`, `estimate`, `hba2_depth`, `hba1_depth`, `flank_left_depth`, `flank_right_depth`, `mapq0_share` (the part of the depth over the two genes from reads of mapping quality 0), `message` | | `${SAMPLE}_cn_intervals.tsv` | One line per interval row of `bin/cn_loci.tsv` (ClinGen's recurrent CNV regions, the carrier deletion alleles, the Y chromosome AZF regions): `locus`, `name`, `chrom`, `start`, `end`, `status`, `copies`, `expected`, `change` (`deletion`, `duplication`, `no change`, or `.` when not called), `estimate`, `reason` ([the interval rows](cn-loci.md#the-interval-rows)). Both reports show it, and the couple screen (step 39) reads it for the carrier deletion alleles and AZF | | `${SAMPLE}_cyp21a2_psv.tsv` | `metric`, `value`: `status`, `cyp21a2_copies` (the least CYP21A2 copies over the low sites, else the median site's), `module_total` (CYP21A2 plus CYP21A1P copies), `expected_total` (4), `estimate`, `finding` (`none`, `deletion`, `conversion`, `pseudogene_deletion`), `conversion_sites`, one `psv_` row per site (gene, pseudogene and other reads, copies, ratio to the other sites), `target_depth`, `control_depth`, `reason`, `reading` | | `${SAMPLE}_parascopy/
/` | Parascopy's own output for each catalog row (`smn/` for SMN1): `res.samples.bed.gz`, `res.paralog.bed.gz`, `psvs.vcf.gz` and the rest ([described upstream](https://github.com/tprodanov/parascopy/blob/main/docs/cn_output.md)) | +## The SMN1 silent-carrier marker + +Two SMN1 copies is the usual result, but some people with two copies have both on one chromosome and none on the other (2+0): a silent carrier, whom copy number alone cannot find. SMN1 g.27134T>G (rs143838139, `NM_000344.4:c.*3+80T>G`, chr5:70,952,074 in GRCh38) is a marker of SMN1 duplication alleles. In Luo et al. 2014 (Genet Med, PMID 23788250) a two-copy result with the marker raised the residual carrier risk (in African Americans from about 1 in 121 to 1 in 34) and one without it lowered it (to about 1 in 396); the figures depend on ancestry, and in Ashkenazi Jewish people a negative result lowered it from about 1 in 345 to 1 in 580. + +SMN2 holds the same 21 bases 875,420 bases upstream, so the reads cannot tell an SMN1 copy of the site from an SMN2 copy. The step adds up the reads at both (every mapping quality, duplicates left out, bases of quality 13 or more) and reads the share with G as a fraction of copies; with a resolved agCN, the copies that carry it are that fraction times agCN. + +| `status` | When | +|---|---| +| `present` | 3 reads or more with G, and at least 5% of the reads | +| `not seen` | 1 read with G or none | +| `uncertain` | between the two | +| `too little depth` | under 10 reads over the two positions | +| `not assessed` | the pileup could not be read, or the check failed | + +The message reads the marker beside the SMN1 psCN at exon 7 (filter `PASS` and quality 20 or more, the rule of the report's SMN1 line): with two SMN1 copies it says whether the marker raises or lowers the chance of 2+0; with another count, or a psCN that is not resolved, it says the marker refines a two-copy result only. Every message ends with "This marker does not prove 2+0, and its absence does not rule 2+0 out." The reports do not read this table yet; it sits beside the SMN table in the same folder. + ## Runtime On the CI fixture with the one-row catalog the step took 8 seconds: one thread, the background depth over 4,000 windows of 100 bp on the chr20 slice, then the SMN1/SMN2 locus on the chr5 slice (case `pgx-4-parascopy`, e2e run 39544031036). The e2e case `paralog-rows-1` times the whole 37-row catalog on HG002 reads over every catalog locus on the fixture's chromosomes and fails past 30 minutes. A 30x genome has not been timed: there the background depth reads Parascopy's windows across the whole genome, which the fixture's slices do not have. diff --git a/docs/42-traits.md b/docs/42-traits.md index 2fa79a21..f9f20113 100644 --- a/docs/42-traits.md +++ b/docs/42-traits.md @@ -113,7 +113,10 @@ The six rules after TAS2R38 need no phase. Their sites lie close together on one - HFE: the C282Y homozygote is also a secondary finding (`assets/acmg_sf.tsv`, [step 23](23-clinical-filter.md)); the rule gives the reading of every genotype. - GC: the protein levels are those of Powe et al. (2013), whose monoclonal assay was questioned in later letters. A protein type, not a vitamin D deficiency. -Two more markers of the planned list are not rules here: the LPA KIV-2 4925G>A null allele (Coassin et al. 2017) lies inside the KIV-2 repeat, whose copies a gVCF genotype at one position cannot tell apart, and the SMN1 g.27134T>G marker (Luo et al. 2014) lies in SMN1, whose reads cannot be told from SMN2's at that site. Both need a copy-aware reader beside [step 46](46-lpa.md) and [step 35](35-paralogs.md)'s SMN1 copy number, and are not read today. +Two more markers of the planned list are not rules here, because each sits in sequence the genome holds more than once and a gVCF genotype at one position cannot read it. They are read from the BAM by the step that counts the copies, as a fraction of copies rather than a genotype (`MARKERS` in `bin/traits_report.py`, which pins each one's dbSNP b157 record and the GRCh38 sequence around every copy): + +- **LPA KIV-2 4925G>A** (Coassin et al. 2017, PMID 28444229), a splice variant that lowers the Lp(a) made from the chromosome that carries it. GRCh38 holds the site in each of its six KIV-2 units and in KIV-1. [Step 46](46-lpa.md) reads it beside the KIV-2 count. +- **SMN1 g.27134T>G** (rs143838139; Luo et al. 2014, PMID 23788250), a marker of SMN1 duplication alleles that SMN2 holds too. [Step 35](35-paralogs.md) reads it beside the SMN1 copy number. It refines the chance of a silent 2+0 carrier for a two-copy result, and its line always says that it does not prove 2+0. Kim et al. (2003) found five TAS2R38 haplotypes worldwide. The rule reads the two common ones. A copy with AAV, AAI, PVI or a rarer one is `other_haplotype`, so the GRCh38 reference itself, which carries AAI (Ala49, Ala262, Ile296), is not read when every site matches it. The readings follow Bufe et al. (2005): PAV/PAV rated PTC and PROP the most bitter, PAV/AVI recognised them as bitter at nearly the same concentrations with a lower and more variable intensity, and AVI/AVI recognised PTC only at far higher concentrations. diff --git a/docs/46-lpa.md b/docs/46-lpa.md index e45d7e9c..09816977 100644 --- a/docs/46-lpa.md +++ b/docs/46-lpa.md @@ -2,7 +2,7 @@ ## What This Does -Estimates how many kringle IV type 2 (KIV-2) repeats the LPA gene carries on both chromosomes together, from read depth, with a 95% interval. The report's Lp(a) section puts that count beside the two common LPA risk variants that [step 42](42-traits.md) reads, rs10455872 and rs3798220, and flags a genetically high Lp(a) when either points that way. Opt-in. +Estimates how many kringle IV type 2 (KIV-2) repeats the LPA gene carries on both chromosomes together, from read depth, with a 95% interval. Beside the count it reads the KIV-2 4925G>A splice variant, which lowers the Lp(a) made from the chromosome that carries it. The report's Lp(a) section puts that count beside the two common LPA risk variants that [step 42](42-traits.md) reads, rs10455872 and rs3798220, and flags a genetically high Lp(a) when either points that way. Opt-in. It is a research estimate, not a diagnostic test. The test that counts is a blood Lp(a) measurement. @@ -15,16 +15,18 @@ Two parts of LPA predict the level: - **The KIV-2 repeat count.** The KIV-2 repeat is a 5.5 kb unit, repeated a variable number of times. Fewer repeats give a smaller apo(a) and, on average, a higher Lp(a). In a Danish population study the total over both chromosomes ranged from 6 to 99, and mean Lp(a) was 56 mg/dL in the quarter with the fewest repeats against 15 mg/dL in the quarter with the most; the count explained 21% to 27% of the variation in Lp(a) level (Kamstrup et al., JAMA 2009, PMID 19509380). - **Two common variants.** rs10455872 and rs3798220 were each associated with a higher Lp(a) level, a smaller apo(a) and coronary disease, with odds ratios of 1.70 and 1.92; one risk allele at either went with an odds ratio of 1.51, two or more with 2.57 (Clarke et al., NEJM 2009, PMID 20032323). The study was of people of European descent. -A VCF cannot show the repeat count: the six KIV-2 units in GRCh38 are near-identical, so short reads from any copy pile up on them. That pile-up is what this step measures. +- **The KIV-2 4925G>A splice variant.** It sits at the donor splice site of the second KIV-2 exon and is usually in one or two of a person's KIV-2 repeats. In a European population 22% carried it, most of them with isoforms of 19 to 25 kringle IV repeats; carriers had a lower Lp(a) than non-carriers across the whole isoform range, and the variant explained 20.6% of the variation in Lp(a) among people with small isoforms (Coassin et al., Eur Heart J 2017, PMID 28444229). + +A VCF cannot show the repeat count or the 4925G>A variant: the six KIV-2 units in GRCh38 are near-identical, so short reads from any copy pile up on them. That pile-up is what this step measures. ## Tool -- samtools `bedcov -j` for the read depth over the exons +- samtools `bedcov -j` for the read depth over the exons, `mpileup` for the reads at the 4925G>A sites - `bin/copy_number.py` (standard library Python), the copy-number engine, on the LPA_KIV2 row of `bin/cn_loci.tsv` for the count, its interval and the reading; `bin/lpa_kiv2.py` lays out the table. See [Copy number from depth](cn-loci.md) ## Docker Image -- `SAMTOOLS_IMAGE` for the depth, `PYTHON_IMAGE` for the count +- `SAMTOOLS_IMAGE` for the depth and the pileup, `PYTHON_IMAGE` for the count and the 4925G>A reading Pinned in `versions.env`; [Image versions](versions.md) lists the current tags. The step adds no image and no reference data. @@ -39,7 +41,7 @@ Pinned in `versions.env`; [Image versions](versions.md) lists the current tags. ./scripts/46-lpa.sh your_name ``` -With `run-all.sh`: `TOOLS=...,lpa`. It is opt-in and needs no data, so a run lists step 46 only when `TOOLS` names it, and it is skipped for a sample with no BAM. In the Nextflow pipeline: `lpa` in `--tools` (processes `LPA_REGIONS`, `LPA_DEPTH` and `LPA`, for each sample with a BAM). The two risk variants come from step 42 (`traits`, on by default). +With `run-all.sh`: `TOOLS=...,lpa`. It is opt-in and needs no data, so a run lists step 46 only when `TOOLS` names it, and it is skipped for a sample with no BAM. In the Nextflow pipeline: `lpa` in `--tools` (processes `LPA_REGIONS`, `LPA_DEPTH` and `LPA`, for each sample with a BAM, with the same commands). The two risk variants come from step 42 (`traits`, on by default). ## How the count is made @@ -79,23 +81,39 @@ The count is not made when: | `too little depth` | the single-copy exons have under 10x | | `not assessed` | reads of MAPQ 0 make up under 20% of the depth over the KIV-2 block (the aligner left out the reads that fit several copies), the KIV-2 depth is too low for any repeat count, or the BAM has no chr6 | +## How the 4925G>A variant is read + +The variant is the last base of the second KIV-2 exon, G>A on LPA's minus strand, so C>T on the plus strand. GRCh38 holds that base at the end of each of its six KIV-2 units (exons 5, 7, ..., 15; `NM_005577.4:c.733G>A` at the first, the one dbSNP names as rs1484919544) and at the near-identical last exon of KIV-1 (exon 3). The 21 bases around it are the same at all seven, so a read from any copy can land on any of them. `samtools mpileup` without a reference reads every read at the seven positions, every mapping quality, duplicates left out, bases of quality 13 or more, and `bin/lpa_kiv2.py null-allele` adds them up. The share of reads with T is a fraction of the copies, not a genotype. With a count, the copies that carry it are that fraction times N + 2 (the KIV-2 repeats plus one KIV-1 per chromosome). + +| `kiv2_4925_status` | When | +|---|---| +| `present` | 5 reads or more with T, and at least 0.8% of the reads (half a copy among 60) | +| `not seen` | 2 reads with T or fewer (one or two can be sequencing errors) | +| `uncertain` | between the two | +| `too little depth` | under 100 reads over the seven positions together | +| `not assessed` | the pileup could not be read | + +What a short-read fraction cannot say: which chromosome or which repeat carries the variant, so whether it sits on the small or the large isoform. A variant in one repeat among many gives few reads, so `not seen` is not proof of absence at low depth. + ## Output | File | Contents | |---|---| -| `lpa/${SAMPLE}_lpa_kiv2.tsv` | `metric`, `value`: `status`, `kiv2_copies`, `ci95_low`, `ci95_high`, `class`, `kiv2_depth`, `flank_depth`, `mapq0_share`, `message`, `advice` | +| `lpa/${SAMPLE}_lpa_kiv2.tsv` | `metric`, `value`: `status`, `kiv2_copies`, `ci95_low`, `ci95_high`, `class`, `kiv2_depth`, `flank_depth`, `mapq0_share`, `message`, `advice`, then the 4925G>A rows `kiv2_4925_status`, `kiv2_4925_alt_reads`, `kiv2_4925_depth`, `kiv2_4925_fraction`, `kiv2_4925_copies`, `kiv2_4925_message`. The `message` ends with the 4925G>A sentence | | `lpa/${SAMPLE}_lpa_kiv2_depth.tsv` | `samtools bedcov -j` over the regions, every mapping quality | | `lpa/${SAMPLE}_lpa_kiv2_depth_q1.tsv` | the same, reads of MAPQ 1 or more | +| `lpa/${SAMPLE}_lpa_kiv2_4925_pileup.tsv` | `samtools mpileup` at the seven 4925G>A positions | -Both reports have an Lp(a) section: the count with its interval and class, the two risk variants from step 42 with their risk allele counts (rs10455872 G, rs3798220 C), a flag when there is at least one risk allele or the count is low, and the advice to measure Lp(a) in blood. Without step 42's table the two variants are shown as not read. +Both reports have an Lp(a) section: the count with its interval and class, the 4925G>A sentence under it (it is the end of the table's `message`), the two risk variants from step 42 with their risk allele counts (rs10455872 G, rs3798220 C), a flag when there is at least one risk allele or the count is low, and the advice to measure Lp(a) in blood. Without step 42's table the two variants are shown as not read. ## Runtime -Seconds: `bedcov` reads only the 29 exons through the BAM index. +Seconds: `bedcov` reads only the 29 exons and `mpileup` the seven 4925G>A positions, through the BAM index. ## Notes - Research estimate from read depth, not a diagnostic test, and not a measurement of Lp(a). Lp(a) also depends on variants this step does not read. - The count is the total over both chromosomes. Which chromosome carries the shorter allele matters for the level, and depth cannot tell. - Exome data can work in principle, since the regions are exons, but the step has only been run on whole-genome reads. -- On the e2e fixture, a slice of HG002 with no LPA, the step reports the count as not assessed or as too little depth. The e2e case `lpa-1` runs it on HG002's real reads over LPA and checks the arithmetic on a copy with half the KIV-2 reads. +- The 4925G>A reading and the count are separate: a sample with too little depth for the count can still have a 4925G>A reading, then without the copies that carry it. +- On the e2e fixture, a slice of HG002 with no LPA, the step reports the count as not assessed or as too little depth. The e2e case `lpa-1` runs it on HG002's real reads over LPA and checks the arithmetic on a copy with half the KIV-2 reads. The case `genotype-rules-2-lpa-null` plants T in a share of HG002's reads at the seven 4925G>A positions (`present`) and sets them all to C (`not seen`). diff --git a/scripts/ci/e2e-case-paths.tsv b/scripts/ci/e2e-case-paths.tsv index 1254dfcb..77ed2a0b 100644 --- a/scripts/ci/e2e-case-paths.tsv +++ b/scripts/ci/e2e-case-paths.tsv @@ -94,6 +94,8 @@ genome-lookup-1 scripts/ci/fake-docker/lib.sh scripts/owngenome-lookup.sh genotype-rule-rows-1-str-extra assets/str_extra_loci.json assets/stub/db/catalog.json scripts/09-expansion-hunter.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh genotype-rules-1-phase assets/genotype_rules.tsv assets/traits.tsv bin/traits_report.py modules/local/traits/* modules/local/vcf_precheck/* scripts/42-traits.sh scripts/ci/check-rule-alleles.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh genotype-rules-2-g2 assets/genotype_rules.tsv bin/collect_summary.py bin/genotype_rules_indels.py bin/traits_report.py +genotype-rules-2-lpa-null bin/cn_loci.tsv bin/copy_number.py bin/lpa_kiv2.py bin/traits_report.py scripts/46-lpa.sh +genotype-rules-2-smn1-marker bin/copy_number.py bin/traits_report.py modules/local/parascopy/* scripts/35-paralogs.sh haplogroup-explore-1-files assets/callability_genes.tsv assets/live/index.html assets/mito_confirmed.tsv assets/viewer/index.html bin/collect_summary.py bin/y_tree_place.py bin/yleaf_run.py modules/local/indexcov/* modules/local/mito_haplogroup/* modules/local/mito_variants/* modules/local/vcf_precheck/* modules/local/y_haplogroup/* scripts/12-mito-haplogroup.sh scripts/16-indexcov.sh scripts/16b-mosdepth.sh scripts/20-mtoolbox.sh scripts/37-y-haplogroup.sh scripts/ci/fake-docker/lib.sh scripts/generate-report.sh scripts/setup.sh scripts/validate-setup.sh hba-short-read-accuracy-1-nf-report assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/mito_confirmed.tsv assets/par_grch38.bed assets/stub/sample.vcf.gz bin/cn_loci.tsv bin/collect_summary.py bin/copy_number.py bin/export_consumer.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/render_report.py bin/results_index.py modules/local/export_consumer/* modules/local/html_report/* modules/local/vcf_precheck/* scripts/24-html-report.sh scripts/40-export-consumer.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh tests/fixtures/consumer_export/sites.tsv.gz hla-second-caller-1-bash assets/pgx_genes.tsv assets/stub/db/hla.dat assets/stub/db/kir.dat bin/cn_loci.tsv bin/copy_number.py bin/cyp2d6_depth_check.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/pgx_outside_calls.py bin/render_report.py scripts/08-hla-typing.sh scripts/24-html-report.sh scripts/32-pypgx.sh scripts/36-pgx-consensus.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh diff --git a/scripts/ci/e2e-case-times.tsv b/scripts/ci/e2e-case-times.tsv index 6bf12c02..1af463f8 100644 --- a/scripts/ci/e2e-case-times.tsv +++ b/scripts/ci/e2e-case-times.tsv @@ -93,6 +93,8 @@ genome-lookup-1 shard 7 - genotype-rule-rows-1-str-extra shard 8 - genotype-rules-1-phase shard 37 - genotype-rules-2-g2 shard 20 - +genotype-rules-2-lpa-null shard 45 - +genotype-rules-2-smn1-marker shard 30 - haplogroup-explore-1-files shard 59 - hba-short-read-accuracy-1-nf-report shard 111 - hla-second-caller-1-bash shard 482 - diff --git a/tests/e2e/genotype-rules-2-lpa-null.sh b/tests/e2e/genotype-rules-2-lpa-null.sh new file mode 100755 index 00000000..8a46aab8 --- /dev/null +++ b/tests/e2e/genotype-rules-2-lpa-null.sh @@ -0,0 +1,134 @@ +#!/usr/bin/env bash +# Step 46's KIV-2 4925G>A reading (bin/lpa_kiv2.py null-allele, MARKERS in +# bin/traits_report.py) on real reads: +# - dbSNP orientation, online: the record pinned for the marker is what +# dbSNP's GRCh38 VCF of the build versions.env pins holds for rs1484919544 +# (same position and REF, the pinned ALT among dbSNP's); control: the same +# compare with REF and ALT swapped fails; +# - HG002's Illumina reads over LPA (chr6:160.53-160.67 Mb, the 35x BWA-MEM +# GRCh38 no-ALT BAM of DeepVariant's case study, as case lpa-1 uses), +# carried onto the fixture's header. Three samples, each run through +# scripts/46-lpa.sh: +# unmodified HG002's own reading, whatever it is, goes to the notes; its +# reads at the seven copies are counted +# planted every 8th read with C at one of the seven copies (a read +# of one M block) gets T there: 'present', with the copies +# that carry it about fraction x (N + 2) +# reference every read with T at a copy gets C, and reads that start +# up to 500 bases before a copy with any other CIGAR are +# dropped: 'not seen' +# and the 4925G>A sentence follows the count's message in each table. +# The Nextflow side runs the same commands; case lpa-1 compares the two +# tables, 4925G>A rows included, byte for byte. +. "$(dirname "$0")/lib.sh" + +G="$GENOME_DIR" T="$CASE_TMP" +SRC="https://storage.googleapis.com/deepvariant/case-study-testdata/HG002.novaseq.pcr-free.35x.dedup.grch38_no_alt.bam" +BASE="${G}/${SAMPLE}/aligned/${SAMPLE}_sorted.bam" +dk() { docker run --rm -u "$(id -u):$(id -g)" -e HOME=/tmp -v "${E2E_WORK}:${E2E_WORK}" -w "$T" "$1" "${@:2}"; } +retry() { local i; for i in 1 2 3 4; do "$@" && return 0; [ "$i" -lt 4 ] && { echo " try ${i} of 4 failed; again in $((i * 15)) s"; sleep $((i * 15)); }; done; return 1; } +metric() { awk -F'\t' -v k="$2" '$1 == k {print $2}' "${G}/$1/lpa/$1_lpa_kiv2.tsv" 2>/dev/null; } +sample() { + rm -rf "${G:?}/$1"; mkdir -p "${G}/$1/aligned" + cp "$2" "${G}/$1/aligned/$1_sorted.bam"; cp "$2.bai" "${G}/$1/aligned/$1_sorted.bam.bai" +} +SITES=$(python3 -c "import sys; sys.path.insert(0, '${REPO}/bin'); import traits_report as t; \ +print(','.join(map(str, t.MARKERS['LPA_KIV2_4925']['copies'])))") +check_eq "seven copies of the 4925G>A site" "$(tr ',' '\n' <<< "$SITES" | wc -l | tr -d ' ')" 7 + +# --- dbSNP orientation, online -------------------------------------------------------- +URL="https://ftp.ncbi.nih.gov/snp/archive/${CONSUMER_DBSNP_BUILD:?}/VCF/GCF_000001405.40.gz" +PINNED=$(python3 -c "import sys; sys.path.insert(0, '${REPO}/bin'); import traits_report as t; \ +d = t.MARKERS['LPA_KIV2_4925']['dbsnp']; print(d[0], d[2], d[3], d[4], sep='\t')") +same_record() { + awk -F'\t' 'NR == 1 { id = $1; pos = $2; ref = $3; alt = $4; next } + $1 == id && $2 == pos && $3 == ref && index("," $4 ",", "," alt ",") { ok = 1 } + END { exit ok ? 0 : 1 }' <(echo "$1") <(echo "$2") +} +LIVE=$(retry docker run --rm "${BCFTOOLS_HTTPS_IMAGE:?}" bcftools query -r NC_000006.12:160640667-160640667 \ + -i 'ID="rs1484919544"' -f '%ID\t%POS\t%REF\t%ALT\n' "$URL" 2>"${T}/dbsnp.err") +echo "pinned: ${PINNED}; dbSNP ${CONSUMER_DBSNP_BUILD}: ${LIVE:-}" +check "the pinned 4925G>A record is dbSNP's" same_record "$PINNED" "$LIVE" +SWAPPED=$(awk -F'\t' -v OFS='\t' '{ t = $3; $3 = $4; $4 = t; print }' <<< "$PINNED") +check "control: the record with REF and ALT swapped does not match dbSNP" \ + test "$(same_record "$SWAPPED" "$LIVE"; echo $?)" -ne 0 + +# --- HG002's reads over LPA, on the fixture's header --------------------------------- +fetch_reads() { + rm -f "${T}/win.bam" "${T}/$(basename "$SRC").bai" + docker run --rm -u "$(id -u):$(id -g)" -e HOME=/tmp -v /etc/ssl/certs:/etc/ssl/certs:ro \ + -e CURL_CA_BUNDLE=/etc/ssl/certs/ca-certificates.crt -v "${E2E_WORK}:${E2E_WORK}" -w "$T" "$SAMTOOLS_IMAGE" \ + samtools view -b -o "${T}/win.bam" "$SRC" chr6:160530000-160670000 +} +check "case 20's BAM is there" test -s "${BASE}.bai" +check "HG002's reads over LPA download" retry fetch_reads +{ dk "$SAMTOOLS_IMAGE" samtools view -H "$BASE" | grep -v -e '^@RG' -e '^@PG' + dk "$SAMTOOLS_IMAGE" samtools view -H "${T}/win.bam" | grep '^@RG' + dk "$SAMTOOLS_IMAGE" samtools view "${T}/win.bam" | awk -F'\t' '$7 == "=" || $7 == "*" || $7 == "chr6"' +} > "${T}/orig.sam" +rm -f "${T}/win.bam" + +# plant MODE: the SAM with the reads at the seven copies changed (see the header). +plant() { + awk -F'\t' -v OFS='\t' -v sites="$SITES" -v mode="$1" ' + BEGIN { n = split(sites, s, ",") } + /^@/ { print; next } + $3 == "chr6" { + for (i = 1; i <= n; i++) { + if (s[i] < $4 || s[i] > $4 + 500) continue + if ($6 !~ /^[0-9]+M$/) { if (mode == "reference") { drop = 1; break }; continue } + if (s[i] >= $4 + length($10)) continue + o = s[i] - $4 + 1; b = substr($10, o, 1) + if (mode == "planted" && b == "C" && (k++ % 8) == 0) { $10 = substr($10, 1, o - 1) "T" substr($10, o + 1); f++ } + if (mode == "reference" && b == "T") { $10 = substr($10, 1, o - 1) "C" substr($10, o + 1); f++ } + } + if (drop) { drop = 0; d++; next } + } + { print } + END { print (f + 0) " bases changed, " (d + 0) " reads dropped" > "/dev/stderr" }' "${T}/orig.sam" +} +plant planted > "${T}/planted.sam" 2> "${T}/planted.log" +plant reference > "${T}/reference.sam" 2> "${T}/reference.log" +echo "planted: $(cat "${T}/planted.log"); reference: $(cat "${T}/reference.log")" +check_ge "planted: bases changed to T" "$(awk '{print $1}' "${T}/planted.log")" 10 + +H="${SAMPLE}nul" +for v in unmodified planted reference; do + src="${T}/${v}.sam"; [ "$v" = unmodified ] && src="${T}/orig.sam" + dk "$SAMTOOLS_IMAGE" samtools sort -o "${T}/${v}.bam" "$src" + dk "$SAMTOOLS_IMAGE" samtools index "${T}/${v}.bam" + sample "${H}${v:0:1}" "${T}/${v}.bam" + run_step 46-lpa.sh "${H}${v:0:1}"; check_step_exit 46-lpa.sh + cat "${G}/${H}${v:0:1}/lpa/${H}${v:0:1}_lpa_kiv2.tsv" 2>/dev/null +done +U="${H}u" P="${H}p" R="${H}r" +check_ge "HG002 reads over LPA" "$(dk "$SAMTOOLS_IMAGE" samtools view -c "${T}/unmodified.bam" 2>/dev/null || echo 0)" 20000 + +check_eq "unmodified: a count" "$(metric "$U" status)" estimated +check "unmodified: a 4925G>A reading of the seven copies" has '^(present|not seen|uncertain)$' "$(metric "$U" kiv2_4925_status)" +check_ge "unmodified: 100 reads or more at the seven copies" "$(metric "$U" kiv2_4925_depth)" 100 +check "unmodified: the pileup is published, one line per copy with reads" \ + test "$(wc -l < "${G}/${U}/lpa/${U}_lpa_kiv2_4925_pileup.tsv" | tr -d ' ')" -eq 7 +check "unmodified: the message ends with the 4925G>A sentence" \ + has "^about [0-9]+ KIV-2 repeats.*LPA KIV-2 4925G>A \(.*\): (present|not seen|uncertain)" "$(metric "$U" message)" + +check_eq "planted: present" "$(metric "$P" kiv2_4925_status)" present +check "planted: the fraction ($(metric "$P" kiv2_4925_fraction)) is about an 8th of the C reads (0.08 to 0.25)" \ + awk -v f="$(metric "$P" kiv2_4925_fraction)" 'BEGIN {exit !(f >= 0.08 && f <= 0.25)}' +check "planted: copies that carry it = fraction x (N + 2)" \ + awk -v f="$(metric "$P" kiv2_4925_fraction)" -v n="$(metric "$P" kiv2_copies)" -v c="$(metric "$P" kiv2_4925_copies)" \ + 'BEGIN {d = f * (n + 2) - c; exit !(c > 0 && d < 0.06 && d > -0.06)}' +check "planted: the sentence names what a short-read fraction cannot say" \ + has "cannot say which chromosome or which repeat" "$(metric "$P" message)" +check_eq "planted: the count is the unmodified one (only bases changed)" "$(metric "$P" kiv2_copies)" "$(metric "$U" kiv2_copies)" + +check_eq "reference: not seen" "$(metric "$R" kiv2_4925_status)" "not seen" +check_eq "reference: no read with T" "$(metric "$R" kiv2_4925_alt_reads)" 0 +check "control: the planted and the reference tables differ" bash -c '! cmp -s "$1" "$2"' _ \ + "${G}/${P}/lpa/${P}_lpa_kiv2.tsv" "${G}/${R}/lpa/${R}_lpa_kiv2.tsv" +echo "genotype-rules-2-lpa-null: HG002 KIV-2 4925G>A $(metric "$U" kiv2_4925_status)," \ + "$(metric "$U" kiv2_4925_alt_reads) of $(metric "$U" kiv2_4925_depth) reads (fraction $(metric "$U" kiv2_4925_fraction))," \ + "KIV-2 count $(metric "$U" kiv2_copies); planted $(metric "$P" kiv2_4925_fraction), copies $(metric "$P" kiv2_4925_copies)" >> "$E2E_NOTES" + +rm -rf "${G:?}/${U}" "${G:?}/${P}" "${G:?}/${R}" "${T:?}"/*.sam "${T:?}"/*.bam "${T:?}"/*.bai +finish diff --git a/tests/e2e/genotype-rules-2-smn1-marker.sh b/tests/e2e/genotype-rules-2-smn1-marker.sh new file mode 100755 index 00000000..5cdbaf03 --- /dev/null +++ b/tests/e2e/genotype-rules-2-smn1-marker.sh @@ -0,0 +1,147 @@ +#!/usr/bin/env bash +# Step 35's SMN1 g.27134T>G reading (bin/traits_report.py --marker +# SMN1_27134) on real reads: +# - dbSNP orientation, online: the record pinned for the marker is what +# dbSNP's GRCh38 VCF of the build versions.env pins holds for rs143838139; +# control: the same compare with REF and ALT swapped fails; +# - HG002's Illumina reads over the marker's SMN1 and SMN2 copies (the 35x +# BWA-MEM GRCh38 no-ALT BAM of DeepVariant's case study), three BAMs: +# unmodified HG002's own reading, whatever it is, goes to the notes +# planted every 4th read with T at a copy (a read of one M block) +# gets G there: 'present' +# reference every read with G at a copy gets T, and reads that start up +# to 500 bases before a copy with any other CIGAR are dropped: +# 'not seen' +# Running Parascopy on these windows would measure nothing, so the SMN +# table beside them is a stand-in written here: two SMN1 copies at exon 7 +# (PASS, quality 30) and agCN 4. Step 35's marker block runs on each (cut +# from scripts/35-paralogs.sh, with scripts/lib/common.sh's run_in), and +# the PARASCOPY module's block, cut from its main.nf, writes the same table +# in its image, byte for byte, on the planted BAM. +# - every message says the marker does not prove 2+0. +. "$(dirname "$0")/lib.sh" + +SRC="https://storage.googleapis.com/deepvariant/case-study-testdata/HG002.novaseq.pcr-free.35x.dedup.grch38_no_alt.bam" +P=HG002smn +W="${GENOME_DIR}/${P}" +rm -rf "$W"; mkdir -p "${W}/aligned" "${W}/paralogs" +retry() { local i; for i in 1 2 3 4; do "$@" && return 0; [ "$i" -lt 4 ] && { echo " try ${i} of 4 failed; again in $((i * 15)) s"; sleep $((i * 15)); }; done; return 1; } +col() { awk -F'\t' -v k="$2" '$1 == k {print $2}' "$1" 2>/dev/null; } +SITES=$(python3 -c "import sys; sys.path.insert(0, '${REPO}/bin'); import traits_report as t; \ +print(','.join(map(str, t.MARKERS['SMN1_27134']['copies'])))") +check_eq "two copies of the marker site (SMN1, SMN2)" "$SITES" "70952074,70076654" + +# --- dbSNP orientation, online -------------------------------------------------------- +URL="https://ftp.ncbi.nih.gov/snp/archive/${CONSUMER_DBSNP_BUILD:?}/VCF/GCF_000001405.40.gz" +PINNED=$(python3 -c "import sys; sys.path.insert(0, '${REPO}/bin'); import traits_report as t; \ +d = t.MARKERS['SMN1_27134']['dbsnp']; print(d[0], d[2], d[3], d[4], sep='\t')") +same_record() { + awk -F'\t' 'NR == 1 { id = $1; pos = $2; ref = $3; alt = $4; next } + $1 == id && $2 == pos && $3 == ref && index("," $4 ",", "," alt ",") { ok = 1 } + END { exit ok ? 0 : 1 }' <(echo "$1") <(echo "$2") +} +LIVE=$(retry docker run --rm "${BCFTOOLS_HTTPS_IMAGE:?}" bcftools query -r NC_000005.10:70952074-70952074 \ + -i 'ID="rs143838139"' -f '%ID\t%POS\t%REF\t%ALT\n' "$URL" 2>"${CASE_TMP}/dbsnp.err") +echo "pinned: ${PINNED}; dbSNP ${CONSUMER_DBSNP_BUILD}: ${LIVE:-}" +check "the pinned SMN1 g.27134T>G record is dbSNP's" same_record "$PINNED" "$LIVE" +SWAPPED=$(awk -F'\t' -v OFS='\t' '{ t = $3; $3 = $4; $4 = t; print }' <<< "$PINNED") +check "control: the record with REF and ALT swapped does not match dbSNP" \ + test "$(same_record "$SWAPPED" "$LIVE"; echo $?)" -ne 0 + +# --- HG002's reads over both copies --------------------------------------------------- +fetch_reads() { + rm -f "${W}/win.bam" "${W}/$(basename "$SRC").bai" + docker run --rm -u "$(id -u):$(id -g)" -e HOME=/tmp -v /etc/ssl/certs:/etc/ssl/certs:ro \ + -e CURL_CA_BUNDLE=/etc/ssl/certs/ca-certificates.crt -v "${GENOME_DIR}:/genome" -w "/genome/${P}" "$SAMTOOLS_IMAGE" \ + samtools view -b -o "/genome/${P}/win.bam" "$SRC" chr5:70075600-70077700 chr5:70951000-70953100 +} +check "HG002's reads over the SMN1 and SMN2 copies download" retry fetch_reads +rm -f "${W}/$(basename "$SRC").bai" +sam view -h "${P}/win.bam" > "${W}/orig.sam" +rm -f "${W}/win.bam" + +plant() { + awk -F'\t' -v OFS='\t' -v sites="$SITES" -v mode="$1" ' + BEGIN { n = split(sites, s, ",") } + /^@/ { print; next } + $3 == "chr5" { + for (i = 1; i <= n; i++) { + if (s[i] < $4 || s[i] > $4 + 500) continue + if ($6 !~ /^[0-9]+M$/) { if (mode == "reference") { drop = 1; break }; continue } + if (s[i] >= $4 + length($10)) continue + o = s[i] - $4 + 1; b = substr($10, o, 1) + if (mode == "planted" && b == "T" && (k++ % 4) == 0) { $10 = substr($10, 1, o - 1) "G" substr($10, o + 1); f++ } + if (mode == "reference" && b == "G") { $10 = substr($10, 1, o - 1) "T" substr($10, o + 1); f++ } + } + if (drop) { drop = 0; d++; next } + } + { print } + END { print (f + 0) " bases changed, " (d + 0) " reads dropped" > "/dev/stderr" }' "${W}/orig.sam" +} +plant planted > "${W}/planted.sam" 2> "${CASE_TMP}/planted.log" +plant reference > "${W}/reference.sam" 2> "${CASE_TMP}/reference.log" +cp "${W}/orig.sam" "${W}/unmodified.sam" +echo "planted: $(cat "${CASE_TMP}/planted.log"); reference: $(cat "${CASE_TMP}/reference.log")" +check_ge "planted: bases changed to G" "$(awk '{print $1}' "${CASE_TMP}/planted.log")" 5 +for v in unmodified planted reference; do + sam sort -o "${P}/${v}.bam" "${P}/${v}.sam" && sam index "${P}/${v}.bam" +done +rm -f "${W}"/*.sam +check_ge "HG002 reads over both copies" "$(sam view -c "${P}/unmodified.bam" 2>/dev/null || echo 0)" 200 + +# The stand-in SMN table: two SMN1 copies at exon 7, agCN 4. +printf 'chrom\tstart\tend\tlocus\tagCN_filter\tagCN\tagCN_qual\tpsCN_filter\tpsCN\tpsCN_qual\thomologous_regions\n' \ + > "${W}/paralogs/${P}_smn_copy_number.tsv" +printf 'chr5\t70925000\t70953000\tSMN1\tPASS\t4\t50\tPASS\t2,2\t30,30\tchr5:70049000-70078000\n' \ + >> "${W}/paralogs/${P}_smn_copy_number.tsv" + +# --- step 35's marker block, with common.sh's run_in ---------------------------------- +sed -n '/^SMN1M_OUT=/,/^rm -rf "\$SMN1M_WORK"$/p' "${REPO}/scripts/35-paralogs.sh" > "${CASE_TMP}/step_body.sh" +check_ge "step 35's marker block was cut from the script" "$(wc -l < "${CASE_TMP}/step_body.sh" | tr -d ' ')" 15 +step_block() { # step_block NAME: the block on ${P}/.bam, its table copied to CASE_TMP/.tsv + { echo 'set -euo pipefail' + printf 'SAMPLE=%q GENOME_DIR=%q\n' "$P" "$GENOME_DIR" + printf '. %q\n' "${REPO}/scripts/lib/common.sh" + printf 'OUTDIR=%q BAM=%q\n' "${W}/paralogs" "${W}/$1.bam" + cat "${CASE_TMP}/step_body.sh"; } > "${CASE_TMP}/block.sh" + echo "+ step 35's marker block ($1)" + bash "${CASE_TMP}/block.sh" 2>&1 | tail -n 5 + cp "${W}/paralogs/${P}_smn1_marker.tsv" "${CASE_TMP}/$1.tsv" 2>/dev/null || : > "${CASE_TMP}/$1.tsv" + cat "${CASE_TMP}/$1.tsv" +} +for v in unmodified planted reference; do step_block "$v"; done +U="${CASE_TMP}/unmodified.tsv" L="${CASE_TMP}/planted.tsv" R="${CASE_TMP}/reference.tsv" + +check "unmodified: a reading" has '^(present|not seen|uncertain)$' "$(col "$U" status)" +check_ge "unmodified: 10 reads or more at the two copies" "$(col "$U" depth)" 10 +check_eq "the SMN1 psCN comes from the SMN table" "$(col "$U" smn1_pscn)" "2 (quality 30, filter PASS)" +check_eq "planted: present" "$(col "$L" status)" present +check "planted: the fraction ($(col "$L" fraction)) is about a quarter of the T reads (0.15 to 0.45)" \ + awk -v f="$(col "$L" fraction)" 'BEGIN {exit !(f >= 0.15 && f <= 0.45)}' +check "planted: copies = fraction x agCN" \ + awk -v f="$(col "$L" fraction)" -v c="$(col "$L" copies)" 'BEGIN {d = f * 4 - c; exit !(c > 0 && d < 0.06 && d > -0.06)}' +check "planted: with two SMN1 copies it raises the chance of 2+0" has "higher chance" "$(col "$L" message)" +check_eq "reference: not seen" "$(col "$R" status)" "not seen" +check_eq "reference: no read with G" "$(col "$R" alt_reads)" 0 +for v in U L R; do + check "${v}: the message says the marker does not prove 2+0" has "does not prove 2\+0" "$(col "${!v}" message)" +done +check "control: the planted and the reference tables differ" bash -c '! cmp -s "$1" "$2"' _ "$L" "$R" +echo "genotype-rules-2-smn1-marker: HG002 SMN1 g.27134T>G $(col "$U" status), $(col "$U" alt_reads) of" \ + "$(col "$U" depth) reads (fraction $(col "$U" fraction)); planted $(col "$L" fraction)" >> "$E2E_NOTES" + +# --- the PARASCOPY module's marker block, in its image -------------------------------- +MOD="${REPO}/modules/local/parascopy/main.nf" +sed -n '/# SMN1 g.27134T>G beside the SMN1 psCN (as scripts\/35-paralogs.sh)/,/rm -f smn1_marker_sites.txt smn1_marker_pileup.tsv/p' "$MOD" \ + | sed -e 's/\\\\/\\/g' -e 's/\\\$/$/g' > "${CASE_TMP}/module_body.sh" +check_ge "the module's marker block was cut from main.nf" "$(wc -l < "${CASE_TMP}/module_body.sh" | tr -d ' ')" 12 +M="${W}/module" +rm -rf "$M"; mkdir -p "$M" +cp "${W}/paralogs/${P}_smn_copy_number.tsv" "$M/" +sed -e "s|\${bam}|/genome/${P}/planted.bam|g" -e "s|\${meta.id}|${P}|g" "${CASE_TMP}/module_body.sh" > "${M}/block.sh" +docker run --rm --user "$(id -u):$(id -g)" -v "${GENOME_DIR}:/genome" -v "${REPO}/bin:/pgp-bin:ro" -w "/genome/${P}/module" \ + "$PARASCOPY_IMAGE" sh -c 'PATH="/pgp-bin:$PATH"; sh -e block.sh' 2>&1 | tail -n 3 +check "the module's block writes step 35's table, byte for byte" cmp -s "$L" "${M}/${P}_smn1_marker.tsv" + +rm -rf "$W" +finish From daa0f55bb1f68a3c8f87bc8f4ffbd6b86add64db Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 19:39:59 +0100 Subject: [PATCH 14/20] test(ci): plant a flipped LPA marker and drop the 2+0 sentence, to watch the new checks go red Reverted in the next commit. Bead: pgp-pfw6.21 --- bin/traits_report.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/bin/traits_report.py b/bin/traits_report.py index e23c5099..20331a2d 100755 --- a/bin/traits_report.py +++ b/bin/traits_report.py @@ -550,7 +550,7 @@ def report(rows, recs, rules=None, phased=None, lines=None): # the block bin/lpa_kiv2.py counts as N + 2 copies (N KIV-2 repeats and # one KIV-1 per chromosome). dbSNP names the first KIV-2 copy only. "LPA_KIV2_4925": { - "name": "LPA KIV-2 4925G>A", "chrom": "chr6", "ref": "C", "alt": "T", + "name": "LPA KIV-2 4925G>A", "chrom": "chr6", "ref": "T", "alt": "C", "dbsnp": ("rs1484919544", "chr6", 160640667, "C", "T"), "copies": (160646214, 160640667, 160635123, 160629574, 160624030, 160618484, 160612932), "context": "GACTCCTTACCTTGTTCGGAA", @@ -735,7 +735,7 @@ def smn1_marker(counts, smn_table=None): return [("marker", m["name"]), ("rsid", m["dbsnp"][0]), ("status", st), ("alt_reads", a["alt_reads"]), ("depth", a["depth"]), ("fraction", a["fraction"]), ("copies", copies), ("smn1_pscn", ps_text), ("agcn", "." if ag is None else str(ag)), - ("message", head + body + " " + SMN1_NOT_PROOF)] + ("message", head + body)] def write_metric_table(path, rows): From 2ef9c14552e711ee6dd8a523d099a2d8847b5b8c Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 19:45:10 +0100 Subject: [PATCH 15/20] Revert the planted flipped LPA marker and the dropped 2+0 sentence The checks went red on the planted commit; this restores the marker and the sentence. Bead: pgp-pfw6.21 --- bin/traits_report.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/bin/traits_report.py b/bin/traits_report.py index 20331a2d..e23c5099 100755 --- a/bin/traits_report.py +++ b/bin/traits_report.py @@ -550,7 +550,7 @@ def report(rows, recs, rules=None, phased=None, lines=None): # the block bin/lpa_kiv2.py counts as N + 2 copies (N KIV-2 repeats and # one KIV-1 per chromosome). dbSNP names the first KIV-2 copy only. "LPA_KIV2_4925": { - "name": "LPA KIV-2 4925G>A", "chrom": "chr6", "ref": "T", "alt": "C", + "name": "LPA KIV-2 4925G>A", "chrom": "chr6", "ref": "C", "alt": "T", "dbsnp": ("rs1484919544", "chr6", 160640667, "C", "T"), "copies": (160646214, 160640667, 160635123, 160629574, 160624030, 160618484, 160612932), "context": "GACTCCTTACCTTGTTCGGAA", @@ -735,7 +735,7 @@ def smn1_marker(counts, smn_table=None): return [("marker", m["name"]), ("rsid", m["dbsnp"][0]), ("status", st), ("alt_reads", a["alt_reads"]), ("depth", a["depth"]), ("fraction", a["fraction"]), ("copies", copies), ("smn1_pscn", ps_text), ("agcn", "." if ag is None else str(ag)), - ("message", head + body)] + ("message", head + body + " " + SMN1_NOT_PROOF)] def write_metric_table(path, rows): From 655087bc7e6f68a2dda815719ffecc3363221b89 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 19:51:31 +0100 Subject: [PATCH 16/20] fix(lpa): add the 4925G>A rows through traits_report.py, so step 46 calls no table layout script The fake-docker case cn-engine-calls holds that step 46 calls copy_number.py and never lpa_kiv2.py. traits_report.py --marker LPA_KIV2_4925 reads the marker and hands the rows to lpa_kiv2.add_null_allele. Bead: pgp-pfw6.21 --- bin/traits_report.py | 19 ++++++++++++++++--- docs/46-lpa.md | 2 +- modules/local/lpa/main.nf | 5 +++-- scripts/46-lpa.sh | 9 +++++---- tests/test_lpa_kiv2.py | 11 +++++++++++ 5 files changed, 36 insertions(+), 10 deletions(-) diff --git a/bin/traits_report.py b/bin/traits_report.py index e23c5099..e9f496ba 100755 --- a/bin/traits_report.py +++ b/bin/traits_report.py @@ -45,12 +45,16 @@ base is a fraction of copies, not a diploid genotype. LPA_KIV2_4925 the LPA KIV-2 4925G>A splice variant (Coassin et al. 2017, PMID 28444229), read by step 46 beside the KIV-2 count - (bin/lpa_kiv2.py null-allele) + (laid out by bin/lpa_kiv2.py) SMN1_27134 SMN1 g.27134T>G (rs143838139; Luo et al. 2014, PMID 23788250), read by step 35 beside the SMN1 psCN traits_report.py --marker-regions NAME print the marker's copy sites as samtools regions (chrom:pos-pos) + traits_report.py --marker LPA_KIV2_4925 --pileup P --marker-out LPA.tsv + add the 4925G>A reading to step 46's table (bin/lpa_kiv2.py + add_null_allele: the kiv2_4925_* rows and a sentence after the + message); exit 2 when LPA.tsv cannot be read or written traits_report.py --marker SMN1_27134 --pileup P [--copy-table SMN.tsv] --marker-out OUT.tsv read `samtools mpileup` at those sites (no reference, every MAPQ) and write OUT.tsv (metricvalue). SMN.tsv is step 35's SMN table: the @@ -769,10 +773,11 @@ def main(argv=None): ap.add_argument("--check-table", metavar="TABLE", help="check the table's schema (and --rules') and exit") ap.add_argument("--marker-regions", metavar="NAME", choices=sorted(MARKERS), help="print a copy-aware marker's sites as samtools regions and exit") - ap.add_argument("--marker", choices=["SMN1_27134"], help="read the SMN1 marker from --pileup") + ap.add_argument("--marker", choices=sorted(MARKERS), help="read a copy-aware marker from --pileup") ap.add_argument("--pileup", help="samtools mpileup at the --marker-regions sites") ap.add_argument("--copy-table", help="step 35's SMN table (the SMN1 psCN and the agCN)") - ap.add_argument("--marker-out", help="the marker's table (metricvalue)") + ap.add_argument("--marker-out", help="SMN1_27134: the marker's table (metricvalue); LPA_KIV2_4925: " + "step 46's table, which gets the kiv2_4925_* rows") a = ap.parse_args(argv) if a.marker_regions: print("\n".join(marker_regions(a.marker_regions))) @@ -780,6 +785,14 @@ def main(argv=None): if a.marker: if not (a.pileup and a.marker_out): ap.error("--marker needs --pileup and --marker-out") + if a.marker == "LPA_KIV2_4925": + import lpa_kiv2 # bin/: step 46's table layout + try: + print("LPA " + lpa_kiv2.add_null_allele(a.marker_out, a.pileup)) + except (OSError, ValueError) as e: + print(f"ERROR: {e}", file=sys.stderr) + return 2 + return 0 try: counts = read_marker_pileup(a.pileup, a.marker) except (OSError, ValueError) as e: diff --git a/docs/46-lpa.md b/docs/46-lpa.md index 09816977..a63abd80 100644 --- a/docs/46-lpa.md +++ b/docs/46-lpa.md @@ -83,7 +83,7 @@ The count is not made when: ## How the 4925G>A variant is read -The variant is the last base of the second KIV-2 exon, G>A on LPA's minus strand, so C>T on the plus strand. GRCh38 holds that base at the end of each of its six KIV-2 units (exons 5, 7, ..., 15; `NM_005577.4:c.733G>A` at the first, the one dbSNP names as rs1484919544) and at the near-identical last exon of KIV-1 (exon 3). The 21 bases around it are the same at all seven, so a read from any copy can land on any of them. `samtools mpileup` without a reference reads every read at the seven positions, every mapping quality, duplicates left out, bases of quality 13 or more, and `bin/lpa_kiv2.py null-allele` adds them up. The share of reads with T is a fraction of the copies, not a genotype. With a count, the copies that carry it are that fraction times N + 2 (the KIV-2 repeats plus one KIV-1 per chromosome). +The variant is the last base of the second KIV-2 exon, G>A on LPA's minus strand, so C>T on the plus strand. GRCh38 holds that base at the end of each of its six KIV-2 units (exons 5, 7, ..., 15; `NM_005577.4:c.733G>A` at the first, the one dbSNP names as rs1484919544) and at the near-identical last exon of KIV-1 (exon 3). The 21 bases around it are the same at all seven, so a read from any copy can land on any of them. `samtools mpileup` without a reference reads every read at the seven positions, every mapping quality, duplicates left out, bases of quality 13 or more, and `bin/traits_report.py --marker LPA_KIV2_4925` adds them up (`bin/lpa_kiv2.py` lays out the rows). The share of reads with T is a fraction of the copies, not a genotype. With a count, the copies that carry it are that fraction times N + 2 (the KIV-2 repeats plus one KIV-1 per chromosome). | `kiv2_4925_status` | When | |---|---| diff --git a/modules/local/lpa/main.nf b/modules/local/lpa/main.nf index 24345b77..aa2be971 100644 --- a/modules/local/lpa/main.nf +++ b/modules/local/lpa/main.nf @@ -12,7 +12,8 @@ with a BAM; empty tables when the BAM has no chr6 LPA bin/copy_number.py call --locus LPA_KIV2: the count, its 95% interval and the message (laid out by bin/lpa_kiv2.py), - then bin/lpa_kiv2.py null-allele: the 4925G>A reading beside it + then bin/traits_report.py --marker LPA_KIV2_4925: the 4925G>A + reading beside it Equivalent to: scripts/46-lpa.sh ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -118,7 +119,7 @@ process LPA { script: """ copy_number.py call --locus LPA_KIV2 --bedcov ${depth} --bedcov-mapq1 ${depth_q1} --table ${meta.id}_lpa_kiv2.tsv - lpa_kiv2.py null-allele --pileup ${pileup} --table ${meta.id}_lpa_kiv2.tsv + traits_report.py --marker LPA_KIV2_4925 --pileup ${pileup} --marker-out ${meta.id}_lpa_kiv2.tsv cat <<-END_VERSIONS > versions.yml "${task.process}": diff --git a/scripts/46-lpa.sh b/scripts/46-lpa.sh index 73ec9fe4..9b9a4e11 100755 --- a/scripts/46-lpa.sh +++ b/scripts/46-lpa.sh @@ -11,8 +11,9 @@ # interval. Beside the count it reads the KIV-2 4925G>A splice variant # (Coassin et al. 2017), which lowers the Lp(a) made from the chromosome that # carries it: samtools mpileup at the site's seven GRCh38 copies (the six -# KIV-2 units and KIV-1), read as a fraction of copies by bin/lpa_kiv2.py -# null-allele (the kiv2_4925_* rows, and one sentence after the message). +# KIV-2 units and KIV-1), read as a fraction of copies by +# bin/traits_report.py --marker LPA_KIV2_4925 (the kiv2_4925_* rows, laid out +# by bin/lpa_kiv2.py, and one sentence after the message). # The two common risk variants, rs10455872 and rs3798220, are rows # of step 42's traits panel; the report's Lp(a) section reads both. A research # estimate, not a diagnostic test: the actionable test is a blood Lp(a) @@ -83,9 +84,9 @@ run_in -v "${PGP_ROOT}/bin:/pgp-bin:ro" --cpus 1 --memory 1g "$PYTHON_IMAGE" \ python3 /pgp-bin/copy_number.py call --locus LPA_KIV2 --bedcov "$(cpath "$DEPTH")" \ --bedcov-mapq1 "$(cpath "$DEPTH_Q1")" --table "$(cpath "$TABLE")" -echo "[4/4] KIV-2 4925G>A beside the count (bin/lpa_kiv2.py null-allele)..." +echo "[4/4] KIV-2 4925G>A beside the count (bin/traits_report.py --marker LPA_KIV2_4925)..." run_in -v "${PGP_ROOT}/bin:/pgp-bin:ro" --cpus 1 --memory 1g "$PYTHON_IMAGE" \ - python3 /pgp-bin/lpa_kiv2.py null-allele --pileup "$(cpath "$PILEUP")" --table "$(cpath "$TABLE")" + python3 /pgp-bin/traits_report.py --marker LPA_KIV2_4925 --pileup "$(cpath "$PILEUP")" --marker-out "$(cpath "$TABLE")" [ -s "$TABLE" ] || { echo "ERROR: no table written (${TABLE})" >&2; exit 1; } echo "" diff --git a/tests/test_lpa_kiv2.py b/tests/test_lpa_kiv2.py index bcc764d0..08eda715 100644 --- a/tests/test_lpa_kiv2.py +++ b/tests/test_lpa_kiv2.py @@ -291,6 +291,17 @@ def table_after(t, name, per_site, n_copies=38, extra=""): subprocess.run([sys.executable, SCRIPT, "null-allele", "--pileup", p, "--table", o], capture_output=True) check("4925G>A: a second run leaves the table as it was (the rows and the sentence once)", open(o).read() == before, open(o).read()) + rc, _, o2, p2 = table_after(t, "viatr", (100, 0)) + os.makedirs(os.path.join(t, "viatr2")) + o3 = os.path.join(t, "viatr2", "S_lpa_kiv2.tsv") + a, b = os.path.join(t, "viatr2", "a.tsv"), os.path.join(t, "viatr2", "b.tsv") + bedcov(a, depths(38)) + bedcov(b, q1(depths(38))) + lpa.main(["estimate", "--bedcov", a, "--bedcov-mapq1", b, "--out", o3]) + res = subprocess.run([sys.executable, os.path.join(REPO, "bin", "traits_report.py"), "--marker", "LPA_KIV2_4925", + "--pileup", p2, "--marker-out", o3], capture_output=True, text=True) + check("4925G>A: traits_report.py --marker LPA_KIV2_4925 (what step 46 and the LPA process run) writes " + "the same table as null-allele", res.returncode == 0 and open(o3).read() == open(o2).read(), res) rc, tab, *_ = table_after(t, "ref", (100, 0)) check("4925G>A: reads with the reference base only is not seen", tab.get("kiv2_4925_status") == "not seen" and tab.get("kiv2_4925_copies") == "." From 46149a4b5483eb5fa949207b478de22e94606846 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 22:00:07 +0100 Subject: [PATCH 17/20] fix(lpa): keep step 46's table as it was and write the 4925G>A reading to its own file Appending rows and a sentence to lpa_kiv2.tsv changed a table that 1.0.0 keeps stable. The reading now goes to lpa/_lpa_kiv2_4925.tsv on both paths, with the count read from step 46's table, and the e2e case compares the module's pileup and table lines with step 46's byte for byte. Bead: pgp-pfw6.21 --- assets/output_contract.tsv | 1 + bin/lpa_kiv2.py | 76 +++++++++----------- bin/traits_report.py | 18 ++--- docs/46-lpa.md | 9 +-- modules/local/lpa/main.nf | 11 +-- scripts/46-lpa.sh | 13 ++-- tests/e2e/genotype-rules-2-lpa-null.sh | 64 +++++++++++------ tests/test_lpa_kiv2.py | 97 ++++++++++++-------------- 8 files changed, 153 insertions(+), 136 deletions(-) diff --git a/assets/output_contract.tsv b/assets/output_contract.tsv index 7fec14eb..5a751524 100644 --- a/assets/output_contract.tsv +++ b/assets/output_contract.tsv @@ -136,6 +136,7 @@ file {sample}/kir/database_release.txt 08 KIR_TYPING nextflow text experimental file {sample}/kir/{sample}_kir_genotype.tsv 08 KIR_TYPING nextflow tsv experimental - - file {sample}/logs/run_status.tsv - - nextflow tsv stable - - file {sample}/lpa/{sample}_lpa_kiv2.tsv 46 LPA nextflow tsv experimental - - +file {sample}/lpa/{sample}_lpa_kiv2_4925.tsv 46 LPA both tsv experimental - metric,value file {sample}/lpa/{sample}_lpa_kiv2_4925_pileup.tsv 46 LPA_DEPTH both tsv experimental - - file {sample}/lpa/{sample}_lpa_kiv2_depth.tsv 46 LPA_DEPTH nextflow tsv experimental - - file {sample}/lpa/{sample}_lpa_kiv2_depth_q1.tsv 46 LPA_DEPTH nextflow tsv experimental - - diff --git a/bin/lpa_kiv2.py b/bin/lpa_kiv2.py index 36ef6144..2b0d3e60 100755 --- a/bin/lpa_kiv2.py +++ b/bin/lpa_kiv2.py @@ -26,14 +26,14 @@ over the same regions, and write LPA.tsv (metricvalue). Exit 0 whatever the status; 2 on unreadable input. Two empty tables mean the BAM has no chr6: status 'not assessed'. - lpa_kiv2.py null-allele --pileup PILEUP.tsv --table LPA.tsv + lpa_kiv2.py null-allele --pileup PILEUP.tsv --table LPA.tsv --out NULL.tsv read `samtools mpileup` (no reference, every MAPQ) at the seven copies of the KIV-2 4925G>A site (`traits_report.py --marker-regions - LPA_KIV2_4925`) and add its reading to LPA.tsv, written by `estimate` - or `copy_number.py call --locus LPA_KIV2`: the kiv2_4925_* rows, and - one sentence after the message, so both reports print it beside the - count. Exit 0 whatever the reading (an unreadable pileup reads 'not - assessed'); 2 when LPA.tsv cannot be read or written. + LPA_KIV2_4925`) and write its reading to NULL.tsv (metricvalue), + with the count N of LPA.tsv (written by `estimate` or `copy_number.py + call --locus LPA_KIV2`, left as it is). Exit 0 whatever the reading (an + unreadable pileup reads 'not assessed'); 2 when NULL.tsv cannot be + written. The regions are exons of the MANE Select transcript ENST00000316300.10, numbered from its first exon: @@ -67,9 +67,9 @@ isoforms, and it is usually in one or two of a person's KIV-2 repeats. GRCh38 holds the site in each of its six units and in KIV-1, so null-allele adds up the reads at all seven (bin/traits_report.py MARKERS) and reads the share -with the variant as a fraction of copies: kiv2_4925_status is 'present', -'not seen', 'uncertain', 'too little depth' or 'not assessed', and with a -count, kiv2_4925_copies is the fraction x (N + 2), the copies that carry it. +with the variant as a fraction of copies: status is 'present', 'not seen', +'uncertain', 'too little depth' or 'not assessed', and with a count, copies +is the fraction x (N + 2), the copies that carry it. A short-read fraction cannot say which chromosome or which repeat carries it, so it cannot say whether it sits on the small or the large isoform. Research estimate, not a diagnostic test: the actionable test is a blood Lp(a) measurement. Standard library only: it runs in the Python image and @@ -162,16 +162,14 @@ def assess(m, q1): "mapq0_share", "message", "advice") -NULL_ORDER = ("kiv2_4925_status", "kiv2_4925_alt_reads", "kiv2_4925_depth", "kiv2_4925_fraction", - "kiv2_4925_copies", "kiv2_4925_message") NULL_LIMIT = ("A short-read fraction cannot say which chromosome or which repeat carries it, so it cannot say " "whether it sits on the small or the large isoform.") def null_allele(counts, kiv2_copies="."): - """[(metric, value)] of the KIV-2 4925G>A reading from + """[(metric, value)] of the KIV-2 4925G>A table from traits_report.read_marker_pileup's counts (None: unreadable), with the - count N of the table ('.' when there is none).""" + count N of step 46's table ('.' when there is none).""" import traits_report as tr # bin/: the copy-aware marker catalog m = tr.MARKERS["LPA_KIV2_4925"] a = tr.assess_marker("LPA_KIV2_4925", counts) @@ -192,45 +190,36 @@ def null_allele(counts, kiv2_copies="."): + (f", under the {m['min_depth']} reads a reading needs" if st == "too little depth" else "") + ".") if st == "present": - msg += (" Carriers had a much lower Lp(a) than non-carriers with the same number of repeats, most " - "of all with few repeats. " + NULL_LIMIT) + msg += (" Carriers had a lower Lp(a) than non-carriers with the same number of repeats, most of " + "all with few repeats. " + NULL_LIMIT) elif st == "not seen": msg += (" Not seen in the reads is not proof of absence: one carrying repeat among many gives few " "reads, and a short-read fraction never says which chromosome or repeat holds a variant.") else: msg += " " + NULL_LIMIT - return [("kiv2_4925_status", st), ("kiv2_4925_alt_reads", a["alt_reads"]), ("kiv2_4925_depth", a["depth"]), - ("kiv2_4925_fraction", a["fraction"]), ("kiv2_4925_copies", copies), ("kiv2_4925_message", msg)] + return [("marker", m["name"]), ("rsid", m["dbsnp"][0] + " (first KIV-2 copy)"), ("status", st), + ("alt_reads", a["alt_reads"]), ("depth", a["depth"]), ("fraction", a["fraction"]), ("copies", copies), + ("kiv2_copies", kiv2_copies), ("message", msg)] -def add_null_allele(table_path, pileup_path): - """Add the 4925G>A rows to step 46's table, and the sentence after its message.""" +def write_null_allele(table_path, pileup_path, out_path): + """Write the KIV-2 4925G>A table beside step 46's table (read for the + count N; a table that is missing or unreadable gives no N). Returns the message.""" import traits_report as tr - with open(table_path, errors="replace") as f: - lines = [ln.rstrip("\n").split("\t", 1) for ln in f if ln.strip()] - if not lines or lines[0] != ["metric", "value"] or any(len(c) != 2 for c in lines): - raise ValueError(f"{table_path}: not a metricvalue table") - rows = [(k, v) for k, v in lines[1:] if not k.startswith("kiv2_4925_")] - kv = dict(rows) + n = "." + try: + with open(table_path, errors="replace") as f: + n = dict(ln.rstrip("\n").split("\t", 1) for ln in f if "\t" in ln).get("kiv2_copies", ".") + except OSError: + pass try: counts = tr.read_marker_pileup(pileup_path, "LPA_KIV2_4925") except (OSError, ValueError) as e: print(f"WARNING: {e}; KIV-2 4925G>A is not assessed", file=sys.stderr) counts = None - new = null_allele(counts, kv.get("kiv2_copies", ".")) - msg = dict(new)["kiv2_4925_message"] - def joined(v): - # a sentence of an earlier run goes first, so a second run adds it once - v = v.split(tr.MARKERS["LPA_KIV2_4925"]["name"] + " (", 1)[0].strip() - return f"{v}{'' if not v or v.endswith('.') else '.'} {msg}".strip() - out = [(k, joined(v) if k == "message" else v) for k, v in rows] + new - tmp = table_path + ".tmp" - with open(tmp, "w") as f: - f.write("metric\tvalue\n") - for k, v in out: - f.write(f"{k}\t{v}\n") - os.replace(tmp, table_path) - return msg + rows = null_allele(counts, n) + tr.write_metric_table(out_path, rows) + return dict(rows)["message"] def main(argv=None): @@ -241,14 +230,15 @@ def main(argv=None): c.add_argument("--bedcov", required=True) c.add_argument("--bedcov-mapq1", required=True, help="samtools bedcov -j -Q 1 over the same regions") c.add_argument("--out", required=True) - n = sub.add_parser("null-allele", help="add the KIV-2 4925G>A reading to the table") + n = sub.add_parser("null-allele", help="the KIV-2 4925G>A table beside the count") n.add_argument("--pileup", required=True, help="samtools mpileup at traits_report.py --marker-regions LPA_KIV2_4925") - n.add_argument("--table", required=True, help="the table estimate wrote") + n.add_argument("--table", required=True, help="the table estimate wrote (for the count)") + n.add_argument("--out", required=True, help="the 4925G>A table") args = ap.parse_args(argv) if args.cmd == "null-allele": try: - print("LPA " + add_null_allele(args.table, args.pileup)) - except (OSError, ValueError) as e: + print("LPA " + write_null_allele(args.table, args.pileup, args.out)) + except OSError as e: print(f"ERROR: {e}", file=sys.stderr) return 2 return 0 diff --git a/bin/traits_report.py b/bin/traits_report.py index e9f496ba..62a21310 100755 --- a/bin/traits_report.py +++ b/bin/traits_report.py @@ -51,10 +51,10 @@ traits_report.py --marker-regions NAME print the marker's copy sites as samtools regions (chrom:pos-pos) - traits_report.py --marker LPA_KIV2_4925 --pileup P --marker-out LPA.tsv - add the 4925G>A reading to step 46's table (bin/lpa_kiv2.py - add_null_allele: the kiv2_4925_* rows and a sentence after the - message); exit 2 when LPA.tsv cannot be read or written + traits_report.py --marker LPA_KIV2_4925 --pileup P --copy-table LPA.tsv --marker-out OUT.tsv + the 4925G>A table beside step 46's table LPA.tsv, whose count gives + the copies that carry it (laid out by bin/lpa_kiv2.py); exit 2 when + OUT.tsv cannot be written traits_report.py --marker SMN1_27134 --pileup P [--copy-table SMN.tsv] --marker-out OUT.tsv read `samtools mpileup` at those sites (no reference, every MAPQ) and write OUT.tsv (metricvalue). SMN.tsv is step 35's SMN table: the @@ -775,9 +775,9 @@ def main(argv=None): help="print a copy-aware marker's sites as samtools regions and exit") ap.add_argument("--marker", choices=sorted(MARKERS), help="read a copy-aware marker from --pileup") ap.add_argument("--pileup", help="samtools mpileup at the --marker-regions sites") - ap.add_argument("--copy-table", help="step 35's SMN table (the SMN1 psCN and the agCN)") - ap.add_argument("--marker-out", help="SMN1_27134: the marker's table (metricvalue); LPA_KIV2_4925: " - "step 46's table, which gets the kiv2_4925_* rows") + ap.add_argument("--copy-table", help="SMN1_27134: step 35's SMN table (the SMN1 psCN and the agCN); " + "LPA_KIV2_4925: step 46's table (the KIV-2 count)") + ap.add_argument("--marker-out", help="the marker's table (metricvalue)") a = ap.parse_args(argv) if a.marker_regions: print("\n".join(marker_regions(a.marker_regions))) @@ -788,8 +788,8 @@ def main(argv=None): if a.marker == "LPA_KIV2_4925": import lpa_kiv2 # bin/: step 46's table layout try: - print("LPA " + lpa_kiv2.add_null_allele(a.marker_out, a.pileup)) - except (OSError, ValueError) as e: + print("LPA " + lpa_kiv2.write_null_allele(a.copy_table or "", a.pileup, a.marker_out)) + except OSError as e: print(f"ERROR: {e}", file=sys.stderr) return 2 return 0 diff --git a/docs/46-lpa.md b/docs/46-lpa.md index a63abd80..75dd202e 100644 --- a/docs/46-lpa.md +++ b/docs/46-lpa.md @@ -83,9 +83,9 @@ The count is not made when: ## How the 4925G>A variant is read -The variant is the last base of the second KIV-2 exon, G>A on LPA's minus strand, so C>T on the plus strand. GRCh38 holds that base at the end of each of its six KIV-2 units (exons 5, 7, ..., 15; `NM_005577.4:c.733G>A` at the first, the one dbSNP names as rs1484919544) and at the near-identical last exon of KIV-1 (exon 3). The 21 bases around it are the same at all seven, so a read from any copy can land on any of them. `samtools mpileup` without a reference reads every read at the seven positions, every mapping quality, duplicates left out, bases of quality 13 or more, and `bin/traits_report.py --marker LPA_KIV2_4925` adds them up (`bin/lpa_kiv2.py` lays out the rows). The share of reads with T is a fraction of the copies, not a genotype. With a count, the copies that carry it are that fraction times N + 2 (the KIV-2 repeats plus one KIV-1 per chromosome). +The variant is the last base of the second KIV-2 exon, G>A on LPA's minus strand, so C>T on the plus strand. GRCh38 holds that base at the end of each of its six KIV-2 units (exons 5, 7, ..., 15; `NM_005577.4:c.733G>A` at the first, the one dbSNP names as rs1484919544) and at the near-identical last exon of KIV-1 (exon 3). The 21 bases around it are the same at all seven, so a read from any copy can land on any of them. `samtools mpileup` without a reference reads every read at the seven positions, every mapping quality, duplicates left out, bases of quality 13 or more, and `bin/traits_report.py --marker LPA_KIV2_4925` adds them up into its own table (`bin/lpa_kiv2.py` lays it out). The share of reads with T is a fraction of the copies, not a genotype. With a count, the copies that carry it are that fraction times N + 2 (the KIV-2 repeats plus one KIV-1 per chromosome). -| `kiv2_4925_status` | When | +| `status` | When | |---|---| | `present` | 5 reads or more with T, and at least 0.8% of the reads (half a copy among 60) | | `not seen` | 2 reads with T or fewer (one or two can be sequencing errors) | @@ -99,12 +99,13 @@ What a short-read fraction cannot say: which chromosome or which repeat carries | File | Contents | |---|---| -| `lpa/${SAMPLE}_lpa_kiv2.tsv` | `metric`, `value`: `status`, `kiv2_copies`, `ci95_low`, `ci95_high`, `class`, `kiv2_depth`, `flank_depth`, `mapq0_share`, `message`, `advice`, then the 4925G>A rows `kiv2_4925_status`, `kiv2_4925_alt_reads`, `kiv2_4925_depth`, `kiv2_4925_fraction`, `kiv2_4925_copies`, `kiv2_4925_message`. The `message` ends with the 4925G>A sentence | +| `lpa/${SAMPLE}_lpa_kiv2.tsv` | `metric`, `value`: `status`, `kiv2_copies`, `ci95_low`, `ci95_high`, `class`, `kiv2_depth`, `flank_depth`, `mapq0_share`, `message`, `advice` | +| `lpa/${SAMPLE}_lpa_kiv2_4925.tsv` | `metric`, `value`: `marker`, `rsid`, `status`, `alt_reads`, `depth`, `fraction`, `copies` (the copies that carry it, `.` without a count or a reading), `kiv2_copies`, `message` | | `lpa/${SAMPLE}_lpa_kiv2_depth.tsv` | `samtools bedcov -j` over the regions, every mapping quality | | `lpa/${SAMPLE}_lpa_kiv2_depth_q1.tsv` | the same, reads of MAPQ 1 or more | | `lpa/${SAMPLE}_lpa_kiv2_4925_pileup.tsv` | `samtools mpileup` at the seven 4925G>A positions | -Both reports have an Lp(a) section: the count with its interval and class, the 4925G>A sentence under it (it is the end of the table's `message`), the two risk variants from step 42 with their risk allele counts (rs10455872 G, rs3798220 C), a flag when there is at least one risk allele or the count is low, and the advice to measure Lp(a) in blood. Without step 42's table the two variants are shown as not read. +Both reports have an Lp(a) section: the count with its interval and class, the two risk variants from step 42 with their risk allele counts (rs10455872 G, rs3798220 C), a flag when there is at least one risk allele or the count is low, and the advice to measure Lp(a) in blood. Without step 42's table the two variants are shown as not read. The reports do not read the 4925G>A table yet; step 46 prints its message beside the count. ## Runtime diff --git a/modules/local/lpa/main.nf b/modules/local/lpa/main.nf index aa2be971..b8a92928 100644 --- a/modules/local/lpa/main.nf +++ b/modules/local/lpa/main.nf @@ -13,7 +13,7 @@ LPA bin/copy_number.py call --locus LPA_KIV2: the count, its 95% interval and the message (laid out by bin/lpa_kiv2.py), then bin/traits_report.py --marker LPA_KIV2_4925: the 4925G>A - reading beside it + reading beside it, _lpa_kiv2_4925.tsv Equivalent to: scripts/46-lpa.sh ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -110,8 +110,9 @@ process LPA { tuple val(meta), path(depth), path(depth_q1), path(pileup) output: - tuple val(meta), path("${meta.id}_lpa_kiv2.tsv"), emit: table - path "versions.yml", emit: versions + tuple val(meta), path("${meta.id}_lpa_kiv2.tsv"), emit: table + tuple val(meta), path("${meta.id}_lpa_kiv2_4925.tsv"), emit: null_allele // KIV-2 4925G>A beside the count + path "versions.yml", emit: versions when: task.ext.when == null || task.ext.when @@ -119,7 +120,8 @@ process LPA { script: """ copy_number.py call --locus LPA_KIV2 --bedcov ${depth} --bedcov-mapq1 ${depth_q1} --table ${meta.id}_lpa_kiv2.tsv - traits_report.py --marker LPA_KIV2_4925 --pileup ${pileup} --marker-out ${meta.id}_lpa_kiv2.tsv + traits_report.py --marker LPA_KIV2_4925 --pileup ${pileup} --copy-table ${meta.id}_lpa_kiv2.tsv \\ + --marker-out ${meta.id}_lpa_kiv2_4925.tsv cat <<-END_VERSIONS > versions.yml "${task.process}": @@ -130,6 +132,7 @@ process LPA { stub: """ printf 'metric\\tvalue\\nstatus\\tnot assessed\\nmessage\\tstub\\n' > ${meta.id}_lpa_kiv2.tsv + printf 'metric\\tvalue\\nstatus\\tnot assessed\\nmessage\\tstub\\n' > ${meta.id}_lpa_kiv2_4925.tsv cat <<-END_VERSIONS > versions.yml "${task.process}": diff --git a/scripts/46-lpa.sh b/scripts/46-lpa.sh index 9b9a4e11..e0914bab 100755 --- a/scripts/46-lpa.sh +++ b/scripts/46-lpa.sh @@ -12,8 +12,8 @@ # (Coassin et al. 2017), which lowers the Lp(a) made from the chromosome that # carries it: samtools mpileup at the site's seven GRCh38 copies (the six # KIV-2 units and KIV-1), read as a fraction of copies by -# bin/traits_report.py --marker LPA_KIV2_4925 (the kiv2_4925_* rows, laid out -# by bin/lpa_kiv2.py, and one sentence after the message). +# bin/traits_report.py --marker LPA_KIV2_4925 into its own table, laid out by +# bin/lpa_kiv2.py. # The two common risk variants, rs10455872 and rs3798220, are rows # of step 42's traits panel; the report's Lp(a) section reads both. A research # estimate, not a diagnostic test: the actionable test is a blood Lp(a) @@ -27,6 +27,8 @@ # its 95% interval, the depths, message # lpa/_lpa_kiv2_depth.tsv samtools bedcov -j over the regions (all MAPQ) # lpa/_lpa_kiv2_depth_q1.tsv the same, reads of MAPQ 1 or more +# lpa/_lpa_kiv2_4925.tsv metricvalue: the 4925G>A reading +# (status, reads, fraction, copies, message) # lpa/_lpa_kiv2_4925_pileup.tsv samtools mpileup at the seven # copies of the 4925G>A site (every MAPQ) set -euo pipefail @@ -45,6 +47,7 @@ TABLE="${OUTDIR}/${SAMPLE}_lpa_kiv2.tsv" DEPTH="${OUTDIR}/${SAMPLE}_lpa_kiv2_depth.tsv" DEPTH_Q1="${OUTDIR}/${SAMPLE}_lpa_kiv2_depth_q1.tsv" PILEUP="${OUTDIR}/${SAMPLE}_lpa_kiv2_4925_pileup.tsv" +NULL="${OUTDIR}/${SAMPLE}_lpa_kiv2_4925.tsv" BED="${OUTDIR}/regions.bed" SITES="${OUTDIR}/sites.txt" @@ -53,7 +56,7 @@ echo "=== Lp(a): LPA KIV-2 repeat count from read depth: ${SAMPLE} ===" [ -f "${BAM}.bai" ] || { echo "ERROR: File not found: ${BAM}.bai" >&2; exit 1; } echo "BAM: ${BAM}" mkdir -p "$OUTDIR" -rm -f "$TABLE" "$DEPTH" "$DEPTH_Q1" "$PILEUP" "$BED" "$SITES" +rm -f "$TABLE" "$DEPTH" "$DEPTH_Q1" "$PILEUP" "$NULL" "$BED" "$SITES" echo "[1/4] The regions (bin/copy_number.py bed --locus LPA_KIV2) and the 4925G>A sites..." # shellcheck disable=SC2016 # $1 and $2 belong to the inner shell @@ -86,8 +89,10 @@ run_in -v "${PGP_ROOT}/bin:/pgp-bin:ro" --cpus 1 --memory 1g "$PYTHON_IMAGE" \ echo "[4/4] KIV-2 4925G>A beside the count (bin/traits_report.py --marker LPA_KIV2_4925)..." run_in -v "${PGP_ROOT}/bin:/pgp-bin:ro" --cpus 1 --memory 1g "$PYTHON_IMAGE" \ - python3 /pgp-bin/traits_report.py --marker LPA_KIV2_4925 --pileup "$(cpath "$PILEUP")" --marker-out "$(cpath "$TABLE")" + python3 /pgp-bin/traits_report.py --marker LPA_KIV2_4925 --pileup "$(cpath "$PILEUP")" \ + --copy-table "$(cpath "$TABLE")" --marker-out "$(cpath "$NULL")" [ -s "$TABLE" ] || { echo "ERROR: no table written (${TABLE})" >&2; exit 1; } echo "" echo "Table: ${TABLE}" +echo "KIV-2 4925G>A: ${NULL}" diff --git a/tests/e2e/genotype-rules-2-lpa-null.sh b/tests/e2e/genotype-rules-2-lpa-null.sh index 8a46aab8..7614ed17 100755 --- a/tests/e2e/genotype-rules-2-lpa-null.sh +++ b/tests/e2e/genotype-rules-2-lpa-null.sh @@ -17,9 +17,10 @@ # reference every read with T at a copy gets C, and reads that start # up to 500 bases before a copy with any other CIGAR are # dropped: 'not seen' -# and the 4925G>A sentence follows the count's message in each table. -# The Nextflow side runs the same commands; case lpa-1 compares the two -# tables, 4925G>A rows included, byte for byte. +# The reading goes to lpa/_lpa_kiv2_4925.tsv; step 46's own table +# keeps its rows. +# - the LPA module's pileup and 4925G>A lines, cut from its main.nf, write the +# same pileup and table in their images, byte for byte, on the planted BAM. . "$(dirname "$0")/lib.sh" G="$GENOME_DIR" T="$CASE_TMP" @@ -28,6 +29,7 @@ BASE="${G}/${SAMPLE}/aligned/${SAMPLE}_sorted.bam" dk() { docker run --rm -u "$(id -u):$(id -g)" -e HOME=/tmp -v "${E2E_WORK}:${E2E_WORK}" -w "$T" "$1" "${@:2}"; } retry() { local i; for i in 1 2 3 4; do "$@" && return 0; [ "$i" -lt 4 ] && { echo " try ${i} of 4 failed; again in $((i * 15)) s"; sleep $((i * 15)); }; done; return 1; } metric() { awk -F'\t' -v k="$2" '$1 == k {print $2}' "${G}/$1/lpa/$1_lpa_kiv2.tsv" 2>/dev/null; } +null() { awk -F'\t' -v k="$2" '$1 == k {print $2}' "${G}/$1/lpa/$1_lpa_kiv2_4925.tsv" 2>/dev/null; } sample() { rm -rf "${G:?}/$1"; mkdir -p "${G}/$1/aligned" cp "$2" "${G}/$1/aligned/$1_sorted.bam"; cp "$2.bai" "${G}/$1/aligned/$1_sorted.bam.bai" @@ -99,36 +101,58 @@ for v in unmodified planted reference; do dk "$SAMTOOLS_IMAGE" samtools index "${T}/${v}.bam" sample "${H}${v:0:1}" "${T}/${v}.bam" run_step 46-lpa.sh "${H}${v:0:1}"; check_step_exit 46-lpa.sh - cat "${G}/${H}${v:0:1}/lpa/${H}${v:0:1}_lpa_kiv2.tsv" 2>/dev/null + cat "${G}/${H}${v:0:1}/lpa/${H}${v:0:1}_lpa_kiv2.tsv" "${G}/${H}${v:0:1}/lpa/${H}${v:0:1}_lpa_kiv2_4925.tsv" 2>/dev/null done U="${H}u" P="${H}p" R="${H}r" check_ge "HG002 reads over LPA" "$(dk "$SAMTOOLS_IMAGE" samtools view -c "${T}/unmodified.bam" 2>/dev/null || echo 0)" 20000 check_eq "unmodified: a count" "$(metric "$U" status)" estimated -check "unmodified: a 4925G>A reading of the seven copies" has '^(present|not seen|uncertain)$' "$(metric "$U" kiv2_4925_status)" -check_ge "unmodified: 100 reads or more at the seven copies" "$(metric "$U" kiv2_4925_depth)" 100 +check "unmodified: a 4925G>A reading of the seven copies" has '^(present|not seen|uncertain)$' "$(null "$U" status)" +check_ge "unmodified: 100 reads or more at the seven copies" "$(null "$U" depth)" 100 check "unmodified: the pileup is published, one line per copy with reads" \ test "$(wc -l < "${G}/${U}/lpa/${U}_lpa_kiv2_4925_pileup.tsv" | tr -d ' ')" -eq 7 -check "unmodified: the message ends with the 4925G>A sentence" \ - has "^about [0-9]+ KIV-2 repeats.*LPA KIV-2 4925G>A \(.*\): (present|not seen|uncertain)" "$(metric "$U" message)" +check "unmodified: the count's message is the count's alone" bash -c '! grep -q 4925G "$1"' _ "${G}/${U}/lpa/${U}_lpa_kiv2.tsv" +check_eq "unmodified: the 4925G>A table carries the count" "$(null "$U" kiv2_copies)" "$(metric "$U" kiv2_copies)" -check_eq "planted: present" "$(metric "$P" kiv2_4925_status)" present -check "planted: the fraction ($(metric "$P" kiv2_4925_fraction)) is about an 8th of the C reads (0.08 to 0.25)" \ - awk -v f="$(metric "$P" kiv2_4925_fraction)" 'BEGIN {exit !(f >= 0.08 && f <= 0.25)}' +check_eq "planted: present" "$(null "$P" status)" present +check "planted: the fraction ($(null "$P" fraction)) is about an 8th of the C reads (0.08 to 0.25)" \ + awk -v f="$(null "$P" fraction)" 'BEGIN {exit !(f >= 0.08 && f <= 0.25)}' check "planted: copies that carry it = fraction x (N + 2)" \ - awk -v f="$(metric "$P" kiv2_4925_fraction)" -v n="$(metric "$P" kiv2_copies)" -v c="$(metric "$P" kiv2_4925_copies)" \ + awk -v f="$(null "$P" fraction)" -v n="$(null "$P" kiv2_copies)" -v c="$(null "$P" copies)" \ 'BEGIN {d = f * (n + 2) - c; exit !(c > 0 && d < 0.06 && d > -0.06)}' check "planted: the sentence names what a short-read fraction cannot say" \ - has "cannot say which chromosome or which repeat" "$(metric "$P" message)" + has "cannot say which chromosome or which repeat" "$(null "$P" message)" check_eq "planted: the count is the unmodified one (only bases changed)" "$(metric "$P" kiv2_copies)" "$(metric "$U" kiv2_copies)" -check_eq "reference: not seen" "$(metric "$R" kiv2_4925_status)" "not seen" -check_eq "reference: no read with T" "$(metric "$R" kiv2_4925_alt_reads)" 0 +check_eq "reference: not seen" "$(null "$R" status)" "not seen" +check_eq "reference: no read with T" "$(null "$R" alt_reads)" 0 check "control: the planted and the reference tables differ" bash -c '! cmp -s "$1" "$2"' _ \ - "${G}/${P}/lpa/${P}_lpa_kiv2.tsv" "${G}/${R}/lpa/${R}_lpa_kiv2.tsv" -echo "genotype-rules-2-lpa-null: HG002 KIV-2 4925G>A $(metric "$U" kiv2_4925_status)," \ - "$(metric "$U" kiv2_4925_alt_reads) of $(metric "$U" kiv2_4925_depth) reads (fraction $(metric "$U" kiv2_4925_fraction))," \ - "KIV-2 count $(metric "$U" kiv2_copies); planted $(metric "$P" kiv2_4925_fraction), copies $(metric "$P" kiv2_4925_copies)" >> "$E2E_NOTES" + "${G}/${P}/lpa/${P}_lpa_kiv2_4925.tsv" "${G}/${R}/lpa/${R}_lpa_kiv2_4925.tsv" +echo "genotype-rules-2-lpa-null: HG002 KIV-2 4925G>A $(null "$U" status)," \ + "$(null "$U" alt_reads) of $(null "$U" depth) reads (fraction $(null "$U" fraction))," \ + "KIV-2 count $(metric "$U" kiv2_copies); planted $(null "$P" fraction), copies $(null "$P" copies)" >> "$E2E_NOTES" -rm -rf "${G:?}/${U}" "${G:?}/${P}" "${G:?}/${R}" "${T:?}"/*.sam "${T:?}"/*.bam "${T:?}"/*.bai +# --- the LPA module's pileup and 4925G>A lines, in their images ---------------------- +MOD="${REPO}/modules/local/lpa/main.nf" +M="${T}/module" +rm -rf "$M"; mkdir -p "$M" +unesc() { sed -e 's/\\\\/\\/g' -e 's/\\\$/$/g' -e "s|\${bam}|${G}/${P}/aligned/${P}_sorted.bam|g" -e "s|\${meta.id}|${P}|g" \ + -e "s|\${pileup}|${P}_lpa_kiv2_4925_pileup.tsv|g"; } +grep -E '^ *traits_report\.py --marker-regions LPA_KIV2_4925 > lpa_sites\.txt$' "$MOD" | unesc > "${M}/sites.sh" +grep -E '^ *while read -r r; do samtools mpileup ' "$MOD" | unesc > "${M}/pileup.sh" +sed -n '/^ *traits_report\.py --marker LPA_KIV2_4925 --pileup/,/--marker-out/p' "$MOD" | unesc > "${M}/null.sh" +check_eq "the module's sites, pileup and 4925G>A lines were cut from main.nf" \ + "$(cat "${M}/sites.sh" "${M}/pileup.sh" "${M}/null.sh" | wc -l | tr -d ' ')" 4 +cat "${M}/sites.sh" "${M}/pileup.sh" "${M}/null.sh" +cp "${G}/${P}/lpa/${P}_lpa_kiv2.tsv" "$M/" +mod() { docker run --rm -u "$(id -u):$(id -g)" -e HOME=/tmp -v "${E2E_WORK}:${E2E_WORK}" -v "${REPO}/bin:/pgp-bin:ro" \ + -w "$M" "$1" sh -c 'PATH="/pgp-bin:$PATH"; set -e; . "./$1"' _ "$2"; } +mod "$PYTHON_IMAGE" sites.sh && mod "$SAMTOOLS_IMAGE" pileup.sh && mod "$PYTHON_IMAGE" null.sh 2>&1 | tail -n 3 +check "the module writes step 46's pileup, byte for byte" \ + cmp -s "${G}/${P}/lpa/${P}_lpa_kiv2_4925_pileup.tsv" "${M}/${P}_lpa_kiv2_4925_pileup.tsv" +check "the module writes step 46's 4925G>A table, byte for byte" \ + cmp -s "${G}/${P}/lpa/${P}_lpa_kiv2_4925.tsv" "${M}/${P}_lpa_kiv2_4925.tsv" +check "the module leaves step 46's table byte for byte" cmp -s "${G}/${P}/lpa/${P}_lpa_kiv2.tsv" "${M}/${P}_lpa_kiv2.tsv" + +rm -rf "${G:?}/${U}" "${G:?}/${P}" "${G:?}/${R}" "${M:?}" "${T:?}"/*.sam "${T:?}"/*.bam "${T:?}"/*.bai finish diff --git a/tests/test_lpa_kiv2.py b/tests/test_lpa_kiv2.py index 08eda715..c4533561 100644 --- a/tests/test_lpa_kiv2.py +++ b/tests/test_lpa_kiv2.py @@ -35,15 +35,14 @@ reference) is 'not seen'; 3 reads with T is 'uncertain'; under 100 reads is 'too little depth'; an unreadable pileup is 'not assessed', exit 0; reads at other positions and read starts, ends and indels in the bases - column are not counted; the rows go after the table and the sentence - after the message, once, and the text report prints it beside the count; - a table that is not metricvalue exits 2. + column are not counted; the reading goes to its own table and step 46's + table is left byte for byte; without a count the copies are '.'; an + output that cannot be written exits 2. Run: python3 tests/test_lpa_kiv2.py """ import math import os -import shutil import subprocess import sys import tempfile @@ -264,80 +263,74 @@ def pileup(path, per_site, extra=""): f.write(extra) def table_after(t, name, per_site, n_copies=38, extra=""): + """Step 46's table, a pileup, then null-allele; returns (rc, the 4925G>A table, LPA table, pileup, out).""" d = os.path.join(t, name) os.makedirs(d) - a, b, o, p = (os.path.join(d, x) for x in ("a.tsv", "b.tsv", "S_lpa_kiv2.tsv", "p.tsv")) + a, b, o, p, n = (os.path.join(d, x) for x in ("a.tsv", "b.tsv", "S_lpa_kiv2.tsv", "p.tsv", "S_lpa_kiv2_4925.tsv")) bedcov(a, depths(n_copies)) bedcov(b, q1(depths(n_copies))) lpa.main(["estimate", "--bedcov", a, "--bedcov-mapq1", b, "--out", o]) pileup(p, per_site, extra) - rc = subprocess.run([sys.executable, SCRIPT, "null-allele", "--pileup", p, "--table", o], + rc = subprocess.run([sys.executable, SCRIPT, "null-allele", "--pileup", p, "--table", o, "--out", n], capture_output=True, text=True) - return rc, dict(x.split("\t", 1) for x in open(o).read().splitlines()[1:]), o, p + tab = dict(x.split("\t", 1) for x in open(n).read().splitlines()[1:]) if os.path.exists(n) else {} + return rc, tab, o, p, n with tempfile.TemporaryDirectory() as t: - rc, tab, o, p = table_after(t, "present", (95, 5)) + rc, tab, o, p, n = table_after(t, "present", (95, 5)) check("4925G>A: 5% of the reads with T at the seven copies is present, exit 0", - rc.returncode == 0 and tab.get("kiv2_4925_status") == "present" - and tab.get("kiv2_4925_alt_reads") == "35" and tab.get("kiv2_4925_depth") == "700", (rc, tab)) + rc.returncode == 0 and tab.get("status") == "present" + and tab.get("alt_reads") == "35" and tab.get("depth") == "700", (rc, tab)) check("4925G>A: copies = fraction x (N + 2) = 0.05 x 40", - tab.get("kiv2_4925_fraction") == "0.0500" and tab.get("kiv2_4925_copies") == "2.0", tab) - check("4925G>A: the sentence follows the count's message, once", - tab["message"].startswith("about 38 KIV-2 repeats") and tab["message"].count("4925G>A") == 1 - and tab["message"].endswith(tab["kiv2_4925_message"]), tab["message"]) + tab.get("fraction") == "0.0500" and tab.get("copies") == "2.0" and tab.get("kiv2_copies") == "38.0", tab) check("4925G>A: the message says what a short-read fraction cannot say", - "cannot say which chromosome or which repeat" in tab["kiv2_4925_message"], tab) - before = open(o).read() - subprocess.run([sys.executable, SCRIPT, "null-allele", "--pileup", p, "--table", o], capture_output=True) - check("4925G>A: a second run leaves the table as it was (the rows and the sentence once)", - open(o).read() == before, open(o).read()) - rc, _, o2, p2 = table_after(t, "viatr", (100, 0)) - os.makedirs(os.path.join(t, "viatr2")) - o3 = os.path.join(t, "viatr2", "S_lpa_kiv2.tsv") - a, b = os.path.join(t, "viatr2", "a.tsv"), os.path.join(t, "viatr2", "b.tsv") - bedcov(a, depths(38)) - bedcov(b, q1(depths(38))) - lpa.main(["estimate", "--bedcov", a, "--bedcov-mapq1", b, "--out", o3]) + "cannot say which chromosome or which repeat" in tab.get("message", ""), tab) + check("4925G>A: the first rows name the marker and its dbSNP record", + open(n).read().splitlines()[:3] == ["metric\tvalue", "marker\tLPA KIV-2 4925G>A", + "rsid\trs1484919544 (first KIV-2 copy)"], open(n).read()) + d0 = os.path.join(t, "plain") + os.makedirs(d0) + a0, b0, o0 = (os.path.join(d0, x) for x in ("a.tsv", "b.tsv", "S_lpa_kiv2.tsv")) + bedcov(a0, depths(38)) + bedcov(b0, q1(depths(38))) + lpa.main(["estimate", "--bedcov", a0, "--bedcov-mapq1", b0, "--out", o0]) + check("4925G>A: step 46's table is left byte for byte", open(o).read() == open(o0).read(), open(o).read()) + n3 = os.path.join(t, "viatr.tsv") res = subprocess.run([sys.executable, os.path.join(REPO, "bin", "traits_report.py"), "--marker", "LPA_KIV2_4925", - "--pileup", p2, "--marker-out", o3], capture_output=True, text=True) + "--pileup", p, "--copy-table", o, "--marker-out", n3], capture_output=True, text=True) check("4925G>A: traits_report.py --marker LPA_KIV2_4925 (what step 46 and the LPA process run) writes " - "the same table as null-allele", res.returncode == 0 and open(o3).read() == open(o2).read(), res) + "the same table as null-allele", res.returncode == 0 and open(n3).read() == open(n).read(), res) + res = subprocess.run([sys.executable, SCRIPT, "null-allele", "--pileup", p, "--table", + os.path.join(t, "none.tsv"), "--out", n3], capture_output=True, text=True) + tab = dict(x.split("\t", 1) for x in open(n3).read().splitlines()[1:]) + check("4925G>A: without step 46's table the reading stands and the copies are '.'", + res.returncode == 0 and tab.get("status") == "present" and tab.get("copies") == "." + and tab.get("kiv2_copies") == ".", (res, tab)) rc, tab, *_ = table_after(t, "ref", (100, 0)) check("4925G>A: reads with the reference base only is not seen", - tab.get("kiv2_4925_status") == "not seen" and tab.get("kiv2_4925_copies") == "." - and "not proof of absence" in tab["kiv2_4925_message"], tab) + tab.get("status") == "not seen" and tab.get("copies") == "." + and "not proof of absence" in tab.get("message", ""), tab) rc, tab, *_ = table_after(t, "unc", (100, 0), extra="chr6\t160640667\tN\t3\tTTT\tIII\n") - check("4925G>A: 3 reads with T is uncertain", tab.get("kiv2_4925_status") == "uncertain" - and tab.get("kiv2_4925_alt_reads") == "3", tab) + check("4925G>A: 3 reads with T is uncertain", tab.get("status") == "uncertain" + and tab.get("alt_reads") == "3", tab) rc, tab, *_ = table_after(t, "thin", (10, 1)) - check("4925G>A: under 100 reads is too little depth", tab.get("kiv2_4925_status") == "too little depth", tab) + check("4925G>A: under 100 reads is too little depth", tab.get("status") == "too little depth", tab) rc, tab, *_ = table_after(t, "other", (100, 0), extra="chr6\t160640668\tN\t50\t" + "T" * 50 + "\t.\nchr7\t160640667\tN\t50\t" + "T" * 50 + "\t.\nchr6\t160635123\tN\t2\tC+2TTC-1T\tII\n") check("4925G>A: reads at other positions and indel bases are not counted", - tab.get("kiv2_4925_status") == "not seen" and tab.get("kiv2_4925_alt_reads") == "0" - and tab.get("kiv2_4925_depth") == "702", tab) - rc, tab, o, p = table_after(t, "bad", (95, 5)) + tab.get("status") == "not seen" and tab.get("alt_reads") == "0" and tab.get("depth") == "702", tab) + rc, tab, o, p, n = table_after(t, "bad", (95, 5)) with open(p, "w") as f: f.write("not a pileup\n") - rc = subprocess.run([sys.executable, SCRIPT, "null-allele", "--pileup", p, "--table", o], + rc = subprocess.run([sys.executable, SCRIPT, "null-allele", "--pileup", p, "--table", o, "--out", n], capture_output=True, text=True) - tab = dict(x.split("\t", 1) for x in open(o).read().splitlines()[1:]) + tab = dict(x.split("\t", 1) for x in open(n).read().splitlines()[1:]) check("4925G>A: an unreadable pileup is not assessed, exit 0", - rc.returncode == 0 and tab.get("kiv2_4925_status") == "not assessed", (rc, tab)) - with open(o, "w") as f: - f.write("a\tb\tc\n") - rc = subprocess.run([sys.executable, SCRIPT, "null-allele", "--pileup", p, "--table", o], - capture_output=True, text=True) - check("4925G>A: a table that is not metricvalue exits 2", rc.returncode == 2, rc) - d = os.path.join(t, "rep", "S") - os.makedirs(os.path.join(d, "lpa")) - _, _, o, _ = table_after(t, "rep2", (95, 5)) - shutil.copy(o, os.path.join(d, "lpa", "S_lpa_kiv2.tsv")) - txt = render_report.text_report(collect_summary.collect("S", d)) - check("4925G>A: the text report prints it beside the count", - "KIV-2 repeats (both chromosomes): 38.0" in txt and "4925G>A" in txt - and txt.index("4925G>A") > txt.index("KIV-2 repeats (both chromosomes)"), txt[-2000:]) + rc.returncode == 0 and tab.get("status") == "not assessed", (rc, tab)) + rc = subprocess.run([sys.executable, SCRIPT, "null-allele", "--pileup", p, "--table", o, + "--out", os.path.join(t, "no", "such", "dir.tsv")], capture_output=True, text=True) + check("4925G>A: an output that cannot be written exits 2", rc.returncode == 2, rc) print("\nRESULT:", "ALL PASS" if FAILS == 0 else f"{FAILS} FAILED") return 1 if FAILS else 0 From 07c92920ccc1d7148d45dede033a9392aaccb4c8 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 22:05:33 +0100 Subject: [PATCH 18/20] fix(traits): read the markers as not assessed when their reads cannot be trusted A failed mpileup at an SMN1 or SMN2 copy could hand the reader a partial pileup, and a failed LPA pileup stopped the KIV-2 count. An aligner that drops the reads of MAPQ 0 leaves out the reads a carrying repeat gives, so the 4925G>A reading is now not assessed there. An SMN table that cannot be read no longer discards the SMN1 marker's reading. Bead: pgp-pfw6.21 --- bin/lpa_kiv2.py | 27 +++++++++++++++++++------- bin/traits_report.py | 6 +++++- docs/46-lpa.md | 2 +- modules/local/lpa/main.nf | 7 +++++-- modules/local/parascopy/main.nf | 4 ++-- scripts/46-lpa.sh | 9 ++++++--- tests/e2e/genotype-rules-2-lpa-null.sh | 4 ++-- tests/test_lpa_kiv2.py | 14 ++++++++++++- tests/test_traits_report.py | 4 ++++ 9 files changed, 58 insertions(+), 19 deletions(-) diff --git a/bin/lpa_kiv2.py b/bin/lpa_kiv2.py index 2b0d3e60..ff8183e6 100755 --- a/bin/lpa_kiv2.py +++ b/bin/lpa_kiv2.py @@ -166,13 +166,15 @@ def assess(m, q1): "whether it sits on the small or the large isoform.") -def null_allele(counts, kiv2_copies="."): +def null_allele(counts, kiv2_copies=".", why_not=None): """[(metric, value)] of the KIV-2 4925G>A table from traits_report.read_marker_pileup's counts (None: unreadable), with the - count N of step 46's table ('.' when there is none).""" + count N of step 46's table ('.' when there is none). why_not: the reason + the reads cannot be trusted (the aligner left out the reads of MAPQ 0), + which makes it 'not assessed' whatever the counts.""" import traits_report as tr # bin/: the copy-aware marker catalog m = tr.MARKERS["LPA_KIV2_4925"] - a = tr.assess_marker("LPA_KIV2_4925", counts) + a = tr.assess_marker("LPA_KIV2_4925", None if why_not else counts) st = a["status"] try: total = float(kiv2_copies) + 2 @@ -183,7 +185,7 @@ def null_allele(counts, kiv2_copies="."): copies = f"{float(a['fraction']) * total:.1f}" name = f"{m['name']} (a splice variant that lowers the Lp(a) made from its chromosome, Coassin 2017)" if st == "not assessed": - msg = f"{name}: not assessed, the reads at its seven GRCh38 copies could not be read." + msg = f"{name}: not assessed, " + (why_not or "the reads at its seven GRCh38 copies could not be read.") else: msg = (f"{name}: {st}, {a['alt_reads']} of {a['depth']} reads at its seven GRCh38 copies carry it" + (f" (about {copies} of the {total:.0f} copies counted)" if copies != "." else "") @@ -206,18 +208,29 @@ def write_null_allele(table_path, pileup_path, out_path): """Write the KIV-2 4925G>A table beside step 46's table (read for the count N; a table that is missing or unreadable gives no N). Returns the message.""" import traits_report as tr - n = "." + t = {} try: with open(table_path, errors="replace") as f: - n = dict(ln.rstrip("\n").split("\t", 1) for ln in f if "\t" in ln).get("kiv2_copies", ".") + t = dict(ln.rstrip("\n").split("\t", 1) for ln in f if "\t" in ln) except OSError: pass + n = t.get("kiv2_copies", ".") + why_not = None + try: + share = float(t.get("mapq0_share", ".")) + except ValueError: + share = None + if share is not None and share < row()["_mapq0_min"]: + # The reads that fit several copies are the ones a carrying repeat + # gives; without them the share at the seven copies is biased. + why_not = (f"reads of mapping quality 0 make up {100 * share:.1f}% of the KIV-2 depth, so the aligner " + "left out the reads that fit several copies, the ones a carrying repeat gives.") try: counts = tr.read_marker_pileup(pileup_path, "LPA_KIV2_4925") except (OSError, ValueError) as e: print(f"WARNING: {e}; KIV-2 4925G>A is not assessed", file=sys.stderr) counts = None - rows = null_allele(counts, n) + rows = null_allele(counts, n, why_not) tr.write_metric_table(out_path, rows) return dict(rows)["message"] diff --git a/bin/traits_report.py b/bin/traits_report.py index 62a21310..cd4fbac9 100755 --- a/bin/traits_report.py +++ b/bin/traits_report.py @@ -709,7 +709,11 @@ def smn1_marker(counts, smn_table=None): 2+0 (SMN1_NOT_PROOF).""" m = MARKERS["SMN1_27134"] a = assess_marker("SMN1_27134", counts) - ps, ps_text, ag = smn_copy_numbers(smn_table) + try: + ps, ps_text, ag = smn_copy_numbers(smn_table) + except OSError as e: + print(f"WARNING: {e}; the SMN copy numbers are not read", file=sys.stderr) + ps, ps_text, ag = None, "not read (the SMN table could not be read)", None copies = "." if a["status"] in ("present", "uncertain") and ag is not None and a["fraction"] != ".": copies = f"{float(a['fraction']) * ag:.1f}" diff --git a/docs/46-lpa.md b/docs/46-lpa.md index 75dd202e..3c0c31f2 100644 --- a/docs/46-lpa.md +++ b/docs/46-lpa.md @@ -91,7 +91,7 @@ The variant is the last base of the second KIV-2 exon, G>A on LPA's minus strand | `not seen` | 2 reads with T or fewer (one or two can be sequencing errors) | | `uncertain` | between the two | | `too little depth` | under 100 reads over the seven positions together | -| `not assessed` | the pileup could not be read | +| `not assessed` | the pileup could not be read (also when `samtools mpileup` failed; the count still runs), or the count found that the aligner left out the reads of MAPQ 0, which are the reads a carrying repeat gives | What a short-read fraction cannot say: which chromosome or which repeat carries the variant, so whether it sits on the small or the large isoform. A variant in one repeat among many gives few reads, so `not seen` is not proof of absence at low depth. diff --git a/modules/local/lpa/main.nf b/modules/local/lpa/main.nf index b8a92928..ac3229c5 100644 --- a/modules/local/lpa/main.nf +++ b/modules/local/lpa/main.nf @@ -71,12 +71,15 @@ process LPA_DEPTH { script: // bedcov stops on a contig the BAM has no @SQ line for: without chr6 the - // tables stay empty. The same commands as scripts/46-lpa.sh. + // tables stay empty. A failed mpileup leaves the count to run: the pileup + // then holds one line that is not mpileup output, so 4925G>A reads 'not + // assessed'. The same commands as scripts/46-lpa.sh. """ if samtools view -H ${bam} | awk -F '\\t' '\$1 == "@SQ" {for (i = 2; i <= NF; i++) if (\$i == "SN:chr6") f = 1} END {exit !f}'; then samtools bedcov -j lpa_regions.bed ${bam} > ${meta.id}_lpa_kiv2_depth.tsv samtools bedcov -j -Q 1 lpa_regions.bed ${bam} > ${meta.id}_lpa_kiv2_depth_q1.tsv - while read -r r; do samtools mpileup -A -B -q 0 -Q 13 -d 100000 -r "\$r" ${bam}; done < lpa_sites.txt > ${meta.id}_lpa_kiv2_4925_pileup.tsv + ( while read -r r; do samtools mpileup -A -B -q 0 -Q 13 -d 100000 -r "\$r" ${bam} || exit 1; done < lpa_sites.txt ) > ${meta.id}_lpa_kiv2_4925_pileup.tsv \\ + || echo '# samtools mpileup failed: 4925G>A not assessed' > ${meta.id}_lpa_kiv2_4925_pileup.tsv else : > ${meta.id}_lpa_kiv2_depth.tsv : > ${meta.id}_lpa_kiv2_depth_q1.tsv diff --git a/modules/local/parascopy/main.nf b/modules/local/parascopy/main.nf index ce481860..6263776b 100644 --- a/modules/local/parascopy/main.nf +++ b/modules/local/parascopy/main.nf @@ -159,8 +159,8 @@ process PARASCOPY { # SMN1 g.27134T>G beside the SMN1 psCN (as scripts/35-paralogs.sh): mpileup at its SMN1 and SMN2 # copies, then the reading with the SMN table; not assessed when it fails if traits_report.py --marker-regions SMN1_27134 > smn1_marker_sites.txt \ - && while read -r r; do samtools mpileup -A -B -q 0 -Q 13 -d 100000 -r "\$r" ${bam}; done \ - < smn1_marker_sites.txt > smn1_marker_pileup.tsv \ + && ( while read -r r; do samtools mpileup -A -B -q 0 -Q 13 -d 100000 -r "\$r" ${bam} || exit 1; done \ + < smn1_marker_sites.txt ) > smn1_marker_pileup.tsv \ && traits_report.py --marker SMN1_27134 --pileup smn1_marker_pileup.tsv \ --copy-table ${meta.id}_smn_copy_number.tsv --marker-out ${meta.id}_smn1_marker.tsv.tmp; then mv ${meta.id}_smn1_marker.tsv.tmp ${meta.id}_smn1_marker.tsv diff --git a/scripts/46-lpa.sh b/scripts/46-lpa.sh index e0914bab..7bbfe417 100755 --- a/scripts/46-lpa.sh +++ b/scripts/46-lpa.sh @@ -68,14 +68,17 @@ echo "[2/4] Read depth over the KIV-2 and single-copy LPA exons (samtools bedcov echo " and the reads at the 4925G>A sites (samtools mpileup, every MAPQ)..." # bedcov stops on a contig the BAM has no @SQ line for: without chr6 the # tables stay empty and the count is not assessed. mpileup runs once per site -# (-r), every MAPQ, duplicates left out, bases of quality 13 or more. The -# LPA_DEPTH process runs the same commands. +# (-r), every MAPQ, duplicates left out, bases of quality 13 or more. A +# failed mpileup leaves the count to run: the pileup then holds one line that +# is not mpileup output, so 4925G>A reads 'not assessed'. The LPA_DEPTH +# process runs the same commands. # shellcheck disable=SC2016 # $1 to $7 belong to the inner bash, the program to awk run_in --cpus 1 --memory 2g "$SAMTOOLS_IMAGE" bash -euo pipefail -c ' if samtools view -H "$2" | awk -F "\t" "$5"; then samtools bedcov -j "$1" "$2" > "$3" samtools bedcov -j -Q 1 "$1" "$2" > "$4" - while read -r r; do samtools mpileup -A -B -q 0 -Q 13 -d 100000 -r "$r" "$2"; done < "$7" > "$6" + ( while read -r r; do samtools mpileup -A -B -q 0 -Q 13 -d 100000 -r "$r" "$2" || exit 1; done < "$7" ) > "$6" \ + || echo "# samtools mpileup failed: 4925G>A not assessed" > "$6" else : > "$3"; : > "$4"; : > "$6"; fi' \ _ "$(cpath "$BED")" "$(cpath "$BAM")" "$(cpath "$DEPTH")" "$(cpath "$DEPTH_Q1")" \ '$1 == "@SQ" {for (i = 2; i <= NF; i++) if ($i == "SN:chr6") f = 1} END {exit !f}' \ diff --git a/tests/e2e/genotype-rules-2-lpa-null.sh b/tests/e2e/genotype-rules-2-lpa-null.sh index 7614ed17..104fbed5 100755 --- a/tests/e2e/genotype-rules-2-lpa-null.sh +++ b/tests/e2e/genotype-rules-2-lpa-null.sh @@ -139,10 +139,10 @@ rm -rf "$M"; mkdir -p "$M" unesc() { sed -e 's/\\\\/\\/g' -e 's/\\\$/$/g' -e "s|\${bam}|${G}/${P}/aligned/${P}_sorted.bam|g" -e "s|\${meta.id}|${P}|g" \ -e "s|\${pileup}|${P}_lpa_kiv2_4925_pileup.tsv|g"; } grep -E '^ *traits_report\.py --marker-regions LPA_KIV2_4925 > lpa_sites\.txt$' "$MOD" | unesc > "${M}/sites.sh" -grep -E '^ *while read -r r; do samtools mpileup ' "$MOD" | unesc > "${M}/pileup.sh" +sed -n '/^ *( while read -r r; do samtools mpileup /,/|| echo /p' "$MOD" | unesc > "${M}/pileup.sh" sed -n '/^ *traits_report\.py --marker LPA_KIV2_4925 --pileup/,/--marker-out/p' "$MOD" | unesc > "${M}/null.sh" check_eq "the module's sites, pileup and 4925G>A lines were cut from main.nf" \ - "$(cat "${M}/sites.sh" "${M}/pileup.sh" "${M}/null.sh" | wc -l | tr -d ' ')" 4 + "$(cat "${M}/sites.sh" "${M}/pileup.sh" "${M}/null.sh" | wc -l | tr -d ' ')" 5 cat "${M}/sites.sh" "${M}/pileup.sh" "${M}/null.sh" cp "${G}/${P}/lpa/${P}_lpa_kiv2.tsv" "$M/" mod() { docker run --rm -u "$(id -u):$(id -g)" -e HOME=/tmp -v "${E2E_WORK}:${E2E_WORK}" -v "${REPO}/bin:/pgp-bin:ro" \ diff --git a/tests/test_lpa_kiv2.py b/tests/test_lpa_kiv2.py index c4533561..f17a37a2 100644 --- a/tests/test_lpa_kiv2.py +++ b/tests/test_lpa_kiv2.py @@ -35,7 +35,8 @@ reference) is 'not seen'; 3 reads with T is 'uncertain'; under 100 reads is 'too little depth'; an unreadable pileup is 'not assessed', exit 0; reads at other positions and read starts, ends and indels in the bases - column are not counted; the reading goes to its own table and step 46's + column are not counted; when step 46's table shows the aligner left out + the reads of MAPQ 0, 'not assessed'; the reading goes to its own table and step 46's table is left byte for byte; without a count the copies are '.'; an output that cannot be written exits 2. @@ -306,6 +307,17 @@ def table_after(t, name, per_site, n_copies=38, extra=""): check("4925G>A: without step 46's table the reading stands and the copies are '.'", res.returncode == 0 and tab.get("status") == "present" and tab.get("copies") == "." and tab.get("kiv2_copies") == ".", (res, tab)) + for share, want in (("0.050", "not assessed"), ("0.900", "present")): + lt = os.path.join(t, f"share{share}.tsv") + with open(lt, "w") as f: + f.write(f"metric\tvalue\nstatus\tnot assessed\nkiv2_copies\t.\nmapq0_share\t{share}\n") + res = subprocess.run([sys.executable, SCRIPT, "null-allele", "--pileup", p, "--table", lt, "--out", n3], + capture_output=True, text=True) + tab = dict(x.split("\t", 1) for x in open(n3).read().splitlines()[1:]) + check(f"4925G>A: a MAPQ 0 share of {share} in step 46's table gives {want} on the same reads", + res.returncode == 0 and tab.get("status") == want + and (want == "present" or "left out the reads that fit several copies" in tab.get("message", "")), + (res, tab)) rc, tab, *_ = table_after(t, "ref", (100, 0)) check("4925G>A: reads with the reference base only is not seen", tab.get("status") == "not seen" and tab.get("copies") == "." diff --git a/tests/test_traits_report.py b/tests/test_traits_report.py index ca368b62..375926cb 100644 --- a/tests/test_traits_report.py +++ b/tests/test_traits_report.py @@ -355,6 +355,10 @@ def table(p, smn_path): t.get("smn1_pscn") == "2 (quality 40, filter PASS)" and t.get("agcn") == "4", t) res, t = table(pileup(15, 5), os.path.join(work, "missing.tsv")) check("SMN1 marker: no SMN table says so", "no SMN table" in t.get("smn1_pscn", ""), t) + res, t = table(pileup(15, 5), work) + check("SMN1 marker: an SMN table that cannot be read keeps the marker's reading, exit 0", + res.returncode == 0 and t.get("status") == "present" and t.get("agcn") == "." + and "could not be read" in t.get("smn1_pscn", ""), (res, t)) # Every status with every psCN: the message says the marker does not prove 2+0. msgs = [] From 9a5229cc44a30307f9ec6112c74d19084b998093 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sun, 8 Nov 2026 01:06:55 +0100 Subject: [PATCH 19/20] test(e2e): select the LPA null-allele case when the LPA module changes The case checks the LPA module's command lines against step 46, so a change to modules/local/lpa alone must run it. --- scripts/ci/e2e-case-paths.tsv | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/ci/e2e-case-paths.tsv b/scripts/ci/e2e-case-paths.tsv index 77ed2a0b..0d3c471c 100644 --- a/scripts/ci/e2e-case-paths.tsv +++ b/scripts/ci/e2e-case-paths.tsv @@ -94,7 +94,7 @@ genome-lookup-1 scripts/ci/fake-docker/lib.sh scripts/owngenome-lookup.sh genotype-rule-rows-1-str-extra assets/str_extra_loci.json assets/stub/db/catalog.json scripts/09-expansion-hunter.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh genotype-rules-1-phase assets/genotype_rules.tsv assets/traits.tsv bin/traits_report.py modules/local/traits/* modules/local/vcf_precheck/* scripts/42-traits.sh scripts/ci/check-rule-alleles.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh genotype-rules-2-g2 assets/genotype_rules.tsv bin/collect_summary.py bin/genotype_rules_indels.py bin/traits_report.py -genotype-rules-2-lpa-null bin/cn_loci.tsv bin/copy_number.py bin/lpa_kiv2.py bin/traits_report.py scripts/46-lpa.sh +genotype-rules-2-lpa-null bin/cn_loci.tsv bin/copy_number.py bin/lpa_kiv2.py bin/traits_report.py modules/local/lpa/* scripts/46-lpa.sh genotype-rules-2-smn1-marker bin/copy_number.py bin/traits_report.py modules/local/parascopy/* scripts/35-paralogs.sh haplogroup-explore-1-files assets/callability_genes.tsv assets/live/index.html assets/mito_confirmed.tsv assets/viewer/index.html bin/collect_summary.py bin/y_tree_place.py bin/yleaf_run.py modules/local/indexcov/* modules/local/mito_haplogroup/* modules/local/mito_variants/* modules/local/vcf_precheck/* modules/local/y_haplogroup/* scripts/12-mito-haplogroup.sh scripts/16-indexcov.sh scripts/16b-mosdepth.sh scripts/20-mtoolbox.sh scripts/37-y-haplogroup.sh scripts/ci/fake-docker/lib.sh scripts/generate-report.sh scripts/setup.sh scripts/validate-setup.sh hba-short-read-accuracy-1-nf-report assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/mito_confirmed.tsv assets/par_grch38.bed assets/stub/sample.vcf.gz bin/cn_loci.tsv bin/collect_summary.py bin/copy_number.py bin/export_consumer.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/render_report.py bin/results_index.py modules/local/export_consumer/* modules/local/html_report/* modules/local/vcf_precheck/* scripts/24-html-report.sh scripts/40-export-consumer.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh tests/fixtures/consumer_export/sites.tsv.gz From c5c465c396f710786c604595b7108d5bfc66a4a8 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sun, 8 Nov 2026 02:46:52 +0100 Subject: [PATCH 20/20] docs(export): say why the round trip leaves chrM out now that chip-to-vcf moves 23andMe MT sites --- docs/40-export-consumer.md | 2 +- tests/e2e/consumer-export-1.sh | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/docs/40-export-consumer.md b/docs/40-export-consumer.md index c2d1bbc5..c281cf47 100644 --- a/docs/40-export-consumer.md +++ b/docs/40-export-consumer.md @@ -81,4 +81,4 @@ On the CI fixture (the gVCF of the fixture's 15 regions, 7,623 sites of the CI t - The AncestryDNA file has the 23andMe v5 sites, not the sites of an AncestryDNA chip; the two chips share part of their sites. A service that matches on the sites it expects uses the ones the files have in common. - The 23andMe file has no internal `i` sites and none of the dropped sites, so it has fewer rows than a file from 23andMe. - In a duplicated region the list's GRCh37 position and dbSNP's GRCh38 position can sit on two different copies. The file then carries the list's position with the genotype read at dbSNP's. On the CI table this is one site of 7,338 off chrM (in SMN1/SMN2), found by the round trip below. -- `scripts/chip-to-vcf.sh` reads both files back. The e2e case `consumer-export-1` does that on the CI fixture, lifts the result to GRCh38 and checks it against the fixture's calls; chrM is left out there, because the hg19 chain maps hg19's chrM, not the rCRS positions the file uses. +- `scripts/chip-to-vcf.sh` reads both files back. The e2e case `consumer-export-1` does that on the CI fixture, lifts the result to GRCh38 and checks it against the fixture's calls; chrM is left out of both files there. `chip-to-vcf.sh` treats the MT positions of a 23andMe-layout file as 23andMe numbers them, on hg19's chrM, and moves each one to the rCRS through the chain's chrM block. This file already writes rCRS positions, so its MT sites would land one or two bases off. The AncestryDNA file's MT positions are read as rCRS positions and keep their place, but the case checks neither. diff --git a/tests/e2e/consumer-export-1.sh b/tests/e2e/consumer-export-1.sh index b858f3ea..238d3aa9 100755 --- a/tests/e2e/consumer-export-1.sh +++ b/tests/e2e/consumer-export-1.sh @@ -15,8 +15,10 @@ # table's GRCh38 position and has the alleles of case 21's calls there # (hom-ref where case 21 called nothing), but for at most 0.1% of the # sites, where the list's GRCh37 position and dbSNP's GRCh38 position -# fall on two copies of a duplicated region. chrM is left out: the hg19 -# chain maps hg19's chrM, not the rCRS positions the file uses; +# fall on two copies of a duplicated region. chrM is left out of both +# files: chip-to-vcf.sh moves a 23andMe-layout file's MT positions from +# hg19's chrM numbering to the rCRS, and step 40 already writes rCRS +# positions, so those sites would land one or two bases off; # 5. red controls: the table shifted by one GRCh38 position stops the # export (reference bases differ); shifted by one GRCh37 position, the # export runs and the round trip fails;