From eca3fc169b3ef6fb12c500442cdc2f3d9779e30e Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Fri, 6 Nov 2026 22:16:15 +0100 Subject: [PATCH 1/8] fix(sv): judge a Delly RDCN on chrX and chrY by the sample's sex and the PARs The class 4-5 postfilter read a Delly RDCN on chrX and chrY as if the sex were unknown, so a male's chrX deletion with RDCN 1 and a female's chrX duplication with RDCN 2 counted as depth-supported, and a duplication inside the PARs needed only 2 copies. Step 05 and ANNOTSV now pass --sex and --par-bed to annotsv_postfilter.py, which applies the rule cnvpytor_to_vcf.py already uses when it adds a Delly call. Bead: pgp-5bja.24 --- bin/annotsv_postfilter.py | 51 ++++++++++++++++++++++++-------- docs/05-annotsv.md | 2 +- modules/local/annotsv/main.nf | 5 ++-- scripts/05-annotsv.sh | 8 +++-- tests/test_annotsv_postfilter.py | 44 +++++++++++++++++++++++++-- 5 files changed, 89 insertions(+), 21 deletions(-) diff --git a/bin/annotsv_postfilter.py b/bin/annotsv_postfilter.py index f09905a7..02f6f355 100755 --- a/bin/annotsv_postfilter.py +++ b/bin/annotsv_postfilter.py @@ -61,9 +61,13 @@ class 4 or 5 without being a real rare deletion: record's own read depth, duphold's DHFFC (a deletion needs < 0.7) or DHBFC (a duplication needs > 1.3), step 15's thresholds; else, for a Delly call over 1,000 bp, - its RDCN (a deletion needs 0 or 1 copies, a duplication - 3 or more, or 2 or more on chrX and chrY, where the sex - is not known here); a CNVpytor call is a read-depth call, + its RDCN, read as bin/cnvpytor_to_vcf.py reads it when + step 05 adds a Delly call: a deletion needs fewer copies + than normal and a duplication more, where normal is 1 on + a male's chrX and chrY outside the PARs (--sex, --par-bed) + and 2 elsewhere; with no --sex either 1 or 2 is normal + there, so a deletion needs 0 or 1 copies and a + duplication 2 or more; a CNVpytor call is a read-depth call, so yes. A Delly call step 05 added (FORMAT GT only) over 1,000 bp is yes: cnvpytor_to_vcf.py set aside the ones whose RDCN does not show the change. "not checked" @@ -94,7 +98,8 @@ class 4 or 5 without being a real rare deletion: other gene is left. Usage: annotsv_postfilter.py --class45 CLASS45.tsv [--genes GENES.bed] - [--rules RULES.tsv] [--consensus CONSENSUS.vcf.gz] --out CHECKED.tsv + [--rules RULES.tsv] [--consensus CONSENSUS.vcf.gz] + [--sex male|female] [--par-bed PAR.bed] --out CHECKED.tsv """ import argparse import csv @@ -102,6 +107,10 @@ class 4 or 5 without being a real rare deletion: import os import sys +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +import cnvpytor_to_vcf # noqa: E402 the Delly RDCN rule step 05 applies when it adds a call +import family_report # noqa: E402 the PAR intervals + EDGE_BP = 5 COMMON_AF = 0.01 CONSENSUS_BP = 1000 @@ -229,8 +238,10 @@ def sv_span(row): return start, max(end, start + 1) -def depth_support(row, caller): - """(yes|no|not checked|., detail) for a deletion or duplication.""" +def depth_support(row, caller, sex="", par=None): + """(yes|no|not checked|., detail) for a deletion or duplication. sex: + male, female or "" (not given); par: family_report.read_par() of the + pseudoautosomal regions, or None.""" t = sv_type(row) if t not in ("DEL", "DUP"): return ".", "" @@ -247,9 +258,17 @@ def depth_support(row, caller): except (KeyError, ValueError): length = 0 if rdcn is not None and rdcn >= 0 and length > DELLY_DEPTH_BP: - sex_chrom = chrom_key(row.get("SV_chrom", "")).upper() in ("X", "Y") - yes = rdcn <= 1 if t == "DEL" else rdcn >= (2 if sex_chrom else 3) - return ("yes" if yes else "no"), f"Delly RDCN {rdcn:g}" + # cnvpytor_to_vcf.py's rule, so a call reads the same in both places: + # the normal copies are 1 on a male's chrX and chrY outside the PARs, + # 2 elsewhere, and either one there when the sex is not given. + chrom = "chr" + chrom_key(row.get("SV_chrom", "")).upper() + try: + start = int(row["SV_start"]) + except (KeyError, ValueError): + start = 0 + no = cnvpytor_to_vcf.delly_depth_contradicts(t, chrom, start, start + length, ["RDCN"], [f"{rdcn:g}"], sex, + family_report.in_par(par or {}, chrom, start)) + return ("no" if no else "yes"), f"Delly RDCN {rdcn:g}" if caller == "delly" and rdcn is None and length > DELLY_DEPTH_BP: # Step 05 adds a Delly call through cnvpytor_to_vcf.py, which keeps # FORMAT GT only, after it has set aside every call over 1,000 bp @@ -325,10 +344,11 @@ def region_rules(row, rules, transcripts): return hits, new_genes, downgrade -def check_rows(rows, transcripts=None, rules=None): +def check_rows(rows, transcripts=None, rules=None, sex="", par=None): """The class 4 and 5 rows with the added columns. transcripts None: the genes file was not given, so the exon edges and gene aliases are not - checked. rules: read_rules() of the rule table, or None.""" + checked. rules: read_rules() of the rule table, or None. sex and par: + as depth_support takes them.""" out = [] for r in rows: o = {k: r.get(k, "") for k in KEEP} @@ -349,7 +369,7 @@ def check_rows(rows, transcripts=None, rules=None): o["Population_source"] = src o["Common"] = "yes" if af is not None and af >= COMMON_AF else "no" o["AnnotSV_class"] = r.get("ACMG_class", "") - o["Depth_support"], depth_detail = depth_support(r, o["Caller"]) + o["Depth_support"], depth_detail = depth_support(r, o["Caller"], sex, par) hits, new_genes, downgrade = region_rules(r, rules or {}, transcripts) if o["Depth_support"] == "no": hits.insert(0, ("no_depth_support", f"read depth does not show the change ({depth_detail})")) @@ -432,8 +452,13 @@ def main(argv=None): ap.add_argument("--genes", help="AnnotSV's genes.RefSeq.sorted.bed") ap.add_argument("--rules", help="assets/sv_rule_regions.tsv") ap.add_argument("--consensus", help="step 22's consensus VCF") + ap.add_argument("--sex", choices=("male", "female", ""), default="", + help="the sample's sex, for the normal copies on chrX and chrY (default: not given)") + ap.add_argument("--par-bed", help="the pseudoautosomal regions (assets/par_grch38.bed)") ap.add_argument("--out", required=True) a = ap.parse_args(argv) + if a.par_bed and not os.path.isfile(a.par_bed): + sys.exit(f"ERROR: no PAR file {a.par_bed}") rows = read_tsv(a.class45) if os.path.getsize(a.class45) else [] if rows and not {"SV_chrom", "SV_start", "SV_end", "SV_type", "ACMG_class"} <= set(rows[0]): sys.exit(f"ERROR: {a.class45} lacks one of SV_chrom, SV_start, SV_end, SV_type, ACMG_class") @@ -447,7 +472,7 @@ def main(argv=None): if not os.path.isfile(a.rules): sys.exit(f"ERROR: no rule table {a.rules}") rules = read_rules(a.rules, {chrom_key(r.get("SV_chrom", "")) for r in rows}) - checked = check_rows(rows, transcripts, rules) + checked = check_rows(rows, transcripts, rules, a.sex, family_report.read_par(a.par_bed)) cols = KEEP + ADDED if a.consensus: for r, s in zip(checked, consensus_support(checked, a.consensus)): diff --git a/docs/05-annotsv.md b/docs/05-annotsv.md index 27da6e35..af11279b 100644 --- a/docs/05-annotsv.md +++ b/docs/05-annotsv.md @@ -117,7 +117,7 @@ AnnotSV gives class 4 or 5 to some calls that are not a rare deletion: a call on Some class 4 and 5 calls are common or unreliable for a reason AnnotSV does not weigh. The table keeps them, with `ACMG_class` 3, AnnotSV's class in `AnnotSV_class` and the reason in `Rule_detail`. The reports list such a call apart from the class 4 and 5 section, with AnnotSV's class, the rule and `Rule_detail`, and leave it out of the section's count. -- **No read-depth support** (`no_depth_support`). A real deletion lowers the read depth and a real duplication raises it. `Depth_support` reads the record's own value: duphold's `DHFFC` below 0.7 for a deletion and `DHBFC` above 1.3 for a duplication, the thresholds of [step 15](15-duphold.md) (Pedersen and Quinlan 2019, PMID 31222198); for a Delly call over 1,000 bp without them, Delly's `RDCN`, 0 or 1 copies for a deletion and 3 or more for a duplication (2 or more on chrX and chrY, since the table does not know the sex). A CNVpytor call is a read-depth call, so it has support. A Delly call step 05 added to AnnotSV's input keeps only `FORMAT/GT`, so its `RDCN` does not reach the table; over 1,000 bp it says `yes`, because `bin/cnvpytor_to_vcf.py` already set aside every such call whose `RDCN` does not show the change. `RDCN` itself is read here only for `SV_VCF` input, where the caller's FORMAT fields survive. `no` moves the call to class 3. A record with no depth value says `not checked` and keeps its class: that happens for a Delly call of 1,000 bp or less, when step 15 did not run (the step then annotates Manta's calls and says so), or for `SV_VCF` input without these fields. Step 15 already drops the Manta calls its depth contradicts, so in a default run this check mostly confirms; it catches an `SV_VCF` file duphold annotated but did not filter. +- **No read-depth support** (`no_depth_support`). A real deletion lowers the read depth and a real duplication raises it. `Depth_support` reads the record's own value: duphold's `DHFFC` below 0.7 for a deletion and `DHBFC` above 1.3 for a duplication, the thresholds of [step 15](15-duphold.md) (Pedersen and Quinlan 2019, PMID 31222198); for a Delly call over 1,000 bp without them, Delly's `RDCN`, judged by the rule `bin/cnvpytor_to_vcf.py` applies when it adds a Delly call: fewer copies than normal for a deletion and more for a duplication, where normal is 1 on a male's chrX and chrY outside the pseudoautosomal regions and 2 everywhere else, PARs included. Step 05 and `ANNOTSV` pass the sample's sex and `assets/par_grch38.bed` for this. With no sex given, 1 and 2 both count as normal on chrX and chrY outside the PARs, so there a deletion needs 0 or 1 copies and a duplication 2 or more. A CNVpytor call is a read-depth call, so it has support. A Delly call step 05 added to AnnotSV's input keeps only `FORMAT/GT`, so its `RDCN` does not reach the table; over 1,000 bp it says `yes`, because `bin/cnvpytor_to_vcf.py` already set aside every such call whose `RDCN` does not show the change. `RDCN` itself is read here only for `SV_VCF` input, where the caller's FORMAT fields survive. `no` moves the call to class 3. A record with no depth value says `not checked` and keeps its class: that happens for a Delly call of 1,000 bp or less, when step 15 did not run (the step then annotates Manta's calls and says so), or for `SV_VCF` input without these fields. Step 15 already drops the Manta calls its depth contradicts, so in a default run this check mostly confirms; it catches an `SV_VCF` file duphold annotated but did not filter. - **Inside a coding tandem repeat** (`tandem_repeat`). An SV whose two breakpoints lie within 10 bp of a tandem repeat that overlaps a MANE Select coding exon changes the repeat's length. Short-read SV callers place such breakpoints loosely and do not size repeats, and AnnotSV can read the change as a frameshift. A length change of a repeat is a question for a repeat genotyper (the repeats of [step 9](09-str-expansions.md)'s catalog are read there), not an SV call. ACMG's PM4 and BP3 (Richards 2015, PMID 25741868) already weigh in-frame changes inside a repeat differently from other length changes. - **The IG and TR loci** (`ig_tr_locus`). An SV more than half inside IGH, IGK, IGL, TRA (with TRD inside it), TRB or TRG. B and T cells cut DNA out of these loci when they rearrange their receptor genes (V(D)J recombination), so DNA from blood carries such deletions in a share of its cells without any inherited change. - **The MHC** (`mhc`). An SV more than half inside the MHC, GRCh38 chr6:28,510,120-33,480,577. Its haplotypes differ by tens of kilobases between people, and gnomAD-SV marks every SV over half inside IGH or the MHC as low confidence (FILTER `IGH_MHC_OVERLAP`, Collins 2020, PMID 32461652). CYP21A2 lies in the MHC; short reads cannot call it here, and [step 35b](35-paralogs.md) reads it from PacBio HiFi reads. diff --git a/modules/local/annotsv/main.nf b/modules/local/annotsv/main.nf index e2fc74c4..b4de6897 100644 --- a/modules/local/annotsv/main.nf +++ b/modules/local/annotsv/main.nf @@ -25,7 +25,8 @@ 5 rows are written to their own table, found by column name, and bin/annotsv_postfilter.py writes each of them with its caller, whether a deletion sits on the exon edges of one transcript of AnnotSV's genes - file, its benign-set frequency, its read-depth support and the class + file, its benign-set frequency, its read-depth support (a Delly RDCN + against the same sex- and PAR-aware normal copies) and the class after the rule table assets/sv_rule_regions.tsv, which moves to class 3 a call inside a coding tandem repeat, the IG and TR loci or the MHC, one the depth contradicts, and one whose gene symbol is a pseudogene's @@ -101,7 +102,7 @@ process ANNOTSV { genes_arg=() if [ -s "\$genes" ]; then genes_arg=(--genes "\$genes"); fi annotsv_postfilter.py --class45 ${meta.id}_sv_class4_5.tsv \${genes_arg[@]+"\${genes_arg[@]}"} \\ - --rules ${sv_rules} --out ${meta.id}_sv_class4_5_checked.tsv + --rules ${sv_rules} --sex '${meta.sex ?: ''}' --par-bed ${par_bed} --out ${meta.id}_sv_class4_5_checked.tsv cat <<-END_VERSIONS > versions.yml "${task.process}": diff --git a/scripts/05-annotsv.sh b/scripts/05-annotsv.sh index bfbcec2e..ec4638a6 100755 --- a/scripts/05-annotsv.sh +++ b/scripts/05-annotsv.sh @@ -174,7 +174,8 @@ fi # caller, whether a deletion sits on the exon edges of one transcript of # AnnotSV's genes file (a processed transcript aligned across an intron), the # benign-set frequency of its type, and whether its read depth shows the -# change. The rule table moves to class 3 a call inside a coding tandem +# change (a Delly RDCN against the same sex- and PAR-aware normal copies as +# above). The rule table moves to class 3 a call inside a coding tandem # repeat, the IG and TR loci or the MHC, one the depth contradicts, and one # whose gene symbol is a pseudogene's alias. How many callers saw it comes # from step 22's consensus, which the reports join when they are written. @@ -187,9 +188,10 @@ else echo "Exon edges: no ${GENES_BED#"$GENOME_DIR"/} in the annotation data, so no deletion is checked against them." fi run_in --cpus 1 --memory 2g -v "${PGP_ROOT}/bin:/pgp-bin:ro" \ - -v "${RULES_TSV}:/pgp/sv_rule_regions.tsv:ro" "${ANNOTSV_IMAGE}" \ + -v "${RULES_TSV}:/pgp/sv_rule_regions.tsv:ro" \ + -v "${PGP_ROOT}/assets/par_grch38.bed:/pgp/par_grch38.bed:ro" "${ANNOTSV_IMAGE}" \ python3 /pgp-bin/annotsv_postfilter.py --class45 "$(cpath "$CLASS45")" ${GENES_ARGS[@]+"${GENES_ARGS[@]}"} \ - --rules /pgp/sv_rule_regions.tsv --out "$(cpath "$CHECKED")" + --rules /pgp/sv_rule_regions.tsv --sex "$SEX" --par-bed /pgp/par_grch38.bed --out "$(cpath "$CHECKED")" echo "=== AnnotSV complete ===" echo "Input: ${INPUT_VCF}" diff --git a/tests/test_annotsv_postfilter.py b/tests/test_annotsv_postfilter.py index 0e84c9b2..0f9b9fe1 100755 --- a/tests/test_annotsv_postfilter.py +++ b/tests/test_annotsv_postfilter.py @@ -32,7 +32,10 @@ Manta deletions at one place, DHFFC 0.95 and 0.5: the first is class 3 (no depth support), the second stays 4; a duplication with DHBFC 1.1 is class 3; Delly RDCN 2 on a deletion is no support, RDCN 1 is; on - chrX a duplication with RDCN 2 counts; a record with no depth value is + chrX a duplication with RDCN 2 counts with no sex given; with --sex + and --par-bed a male's chrX deletion needs RDCN 0, a female's chrX + duplication 3 or more, and inside the PAR a duplication needs 3 for + any sex; a record with no depth value is "not checked" and keeps its class; a CNVpytor call has support; a Delly call with FORMAT GT only (as step 05 adds it) has support over 1,000 bp and is "not checked" at 800 bp; the @@ -157,6 +160,29 @@ def run7(work, name, with_rules=True): return {r["ID"]: r for r in csv.DictReader(f, delimiter="\t")}, p.stderr +SEX_ROWS = "".join([ + row7("XDEL1", "X", 50000000, 50050000, "DEL", 4, "delly", "G20", "GT:RDCN", "1/1:1"), + row7("XDEL0", "X", 51000000, 51050000, "DEL", 4, "delly", "G21", "GT:RDCN", "1/1:0"), + row7("XDUP2", "X", 52000000, 52050000, "DUP", 4, "delly", "G22", "GT:RDCN", "0/1:2"), + row7("PARDUP2", "X", 1000000, 1050000, "DUP", 4, "delly", "G23", "GT:RDCN", "0/1:2"), +]) + + +def run_sex(work, name, sex, par): + """{ID: written row} of the command line on SEX_ROWS with --sex and --par-bed.""" + c45 = os.path.join(work, name + ".tsv") + with open(c45, "w") as f: + f.write("\t".join(HEAD7) + "\n" + SEX_ROWS) + dst = os.path.join(work, name + "_checked.tsv") + p = subprocess.run([sys.executable, os.path.join(REPO, "bin", "annotsv_postfilter.py"), "--class45", c45, + "--sex", sex, "--par-bed", par, "--out", dst], capture_output=True, text=True) + if p.returncode != 0 or not os.path.isfile(dst): + print(f" (the command failed: {p.stderr.strip()[-300:]})") + return {}, p.stderr + with open(dst) as f: + return {r["ID"]: r for r in csv.DictReader(f, delimiter="\t")}, p.stderr + + def rules_section(work): got, err = run7(work, "rules") @@ -193,7 +219,7 @@ def col(i, c): check("Delly RDCN 2 on a deletion: no support; RDCN 1: support", (col("DELLY2", "Depth_support"), col("DELLY1", "Depth_support")) == ("no", "yes"), (got.get("DELLY2"), got.get("DELLY1"))) - check("a chrX duplication with RDCN 2 counts (the sex is not known here)", col("XDUP", "Depth_support") == "yes", + check("a chrX duplication with RDCN 2 counts when the sex is not given", col("XDUP", "Depth_support") == "yes", got.get("XDUP")) check("no depth value: not checked, class kept", (col("NODEPTH", "Depth_support"), col("NODEPTH", "ACMG_class")) == ("not checked", "4"), got.get("NODEPTH")) @@ -203,6 +229,20 @@ def col(i, c): check("the same call of 800 bp had no RDCN gate: not checked, class kept", (col("DELLY_SHORT", "Depth_support"), col("DELLY_SHORT", "ACMG_class")) == ("not checked", "4"), got.get("DELLY_SHORT")) + # Delly RDCN on chrX and chrY by sex and PAR (--sex, --par-bed), the rule + # cnvpytor_to_vcf.py applies when step 05 adds a Delly call. + par = os.path.join(REPO, "assets", "par_grch38.bed") + for sex, want_sex in (("male", {"XDEL1": ("no", "3"), "XDUP2": ("yes", "4"), "PARDUP2": ("no", "3"), + "XDEL0": ("yes", "4")}), + ("female", {"XDEL1": ("yes", "4"), "XDUP2": ("no", "3"), "PARDUP2": ("no", "3"), + "XDEL0": ("yes", "4")}), + ("", {"XDEL1": ("yes", "4"), "XDUP2": ("yes", "4"), "PARDUP2": ("no", "3"), + "XDEL0": ("yes", "4")})): + sx, _ = run_sex(work, "sex_" + (sex or "none"), sex, par) + got_sex = {i: (sx.get(i, {}).get("Depth_support"), sx.get(i, {}).get("ACMG_class")) for i in want_sex} + check(f"sex {sex or 'not given'}: chrX RDCN against 1 copy (male), 2 (female), either (not given); " + f"inside the PAR 2 copies", got_sex == want_sex, got_sex) + moved = sorted(i for i, r in got.items() if r["ACMG_class"] != r["AnnotSV_class"]) want = sorted(["TR_IN", "IGH_IN", "IGH_MOST", "MHC_IN", "ALIAS", "D_LOW", "DUP_LOW", "DELLY2"]) check("exactly the eight rows above are moved to class 3", moved == want, moved) From bcc082e4755e4a3aa4d13d7ec5da03aec41103c5 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Fri, 6 Nov 2026 22:16:19 +0100 Subject: [PATCH 2/8] fix(sv): stop step 05 on an SV VCF with more than one sample The reports read the genotype and phase of the first sample column of AnnotSV's input, so a multi-sample SV_VCF let another person's genotype, or a 0/0, decide whether a deletion counts. sv_svlen_fix.py, the first pass of both step 05 and ANNOTSV, now exits 2 and names the samples. Bead: pgp-5bja.24 --- bin/sv_svlen_fix.py | 12 ++++++++++++ docs/05-annotsv.md | 2 +- tests/test_annotsv_sv_inputs.py | 11 +++++++++++ 3 files changed, 24 insertions(+), 1 deletion(-) diff --git a/bin/sv_svlen_fix.py b/bin/sv_svlen_fix.py index 30d8ac56..7c119463 100755 --- a/bin/sv_svlen_fix.py +++ b/bin/sv_svlen_fix.py @@ -18,6 +18,12 @@ breakend records pass unchanged. A rewritten record keeps its old values in ORIG_SVLEN and ORIG_END. +The input holds one sample. Step 05 and the reports read the genotype and +phase of the first sample column (bin/annotsv_postfilter.py sample_format), +so a VCF with more sample columns would let another person's genotype, or a +0/0, decide whether a deletion counts: such an input stops the script (exit +2) with the sample names. + Usage: sv_svlen_fix.py IN.vcf[.gz]|- > OUT.vcf """ import gzip @@ -82,6 +88,12 @@ def main(argv): for line in fh: if line.startswith("#"): if line.startswith("#CHROM"): + samples = line.rstrip("\n").split("\t")[9:] + if len(samples) > 1: + sys.stderr.write("ERROR: the SV VCF holds {0} samples ({1}); AnnotSV's input must hold one. " + "Pick it first: bcftools view -s \n".format( + len(samples), ", ".join(samples[:5]) + (", ..." if len(samples) > 5 else ""))) + return 2 for h in NEW_INFO: if h not in seen: out.write(h + "\n") diff --git a/docs/05-annotsv.md b/docs/05-annotsv.md index af11279b..406d77c4 100644 --- a/docs/05-annotsv.md +++ b/docs/05-annotsv.md @@ -91,7 +91,7 @@ docker run --rm --user "$(id -u):$(id -g)" \ -bracketedAltMode DEL,DUP,INS ``` -To annotate another SV VCF (for example Sniffles2 output), set `SV_VCF` to its host path; the file must be inside `${GENOME_DIR}`. +To annotate another SV VCF (for example Sniffles2 output), set `SV_VCF` to its host path; the file must be inside `${GENOME_DIR}` and hold one sample. The reports read the genotype and phase of the first sample column, so a VCF with more sample columns stops the step (`bin/sv_svlen_fix.py`) instead of letting another person's genotype decide; pick the sample first with `bcftools view -s `. ## Output - `${GENOME_DIR}/${SAMPLE}/annotsv/${SAMPLE}_sv_annotated.tsv` — main annotated output (one row per SV, with ACMG class) diff --git a/tests/test_annotsv_sv_inputs.py b/tests/test_annotsv_sv_inputs.py index e014e5ce..3f4693bd 100755 --- a/tests/test_annotsv_sv_inputs.py +++ b/tests/test_annotsv_sv_inputs.py @@ -226,6 +226,17 @@ def test_svlen(w): rc, out2, err2 = run("sv_svlen_fix.py", again) check("a second run changes nothing", rc == 0 and out2 == out, err2) + # One sample only: the reports read the GT and PS of the first sample + # column, so a second person's column (here the carrier, after a 0/0) + # must stop the step, not be read silently. + two = os.path.join(w, "two.vcf") + with open(two, "w") as f: + f.write(SV.replace("\ts1\n", "\ts1\ts2\n") + + "chr20\t1300000\tsym\tN\t\t50\tPASS\tEND=1310000;SVTYPE=DEL\tGT\t0/0\t0/1\n") + rc, out3, err3 = run("sv_svlen_fix.py", two) + check("a VCF with two sample columns stops (exit 2) and names them", + rc == 2 and "holds 2 samples (s1, s2)" in err3 and not records(out3), (rc, err3)) + BASE = ( "##fileformat=VCFv4.1\n##contig=\n##contig=\n" From 6f406670d6ed5570532e37b8cf1ae9365c9f59b7 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Fri, 6 Nov 2026 22:23:16 +0100 Subject: [PATCH 3/8] fix(recheck): see image digest changes, skip opt-in steps the run left out, never crash on a bad header An image re-pinned to a new digest at the same tag compared equal, because the check reduced every image to its tag; it now compares the full reference. A pin with no upstream lookup that the run's manifest lacks is now "unknown" instead of "no-newer-release". A pin read only by steps run-all.sh recorded as skipped is left out. A malformed Last-Modified of a PGS file and an offline release with no tag are now errors of that lookup, not a crash. --refresh-clinvar now holds GENOME_DIR/clinvar/.refresh.lock exclusive while setup.sh refreshes; run-all.sh and step 06 hold it shared, so a run that starts during a refresh waits for it. Bead: pgp-5bja.23 Bead: pgp-ax3g --- bin/owngenome_recheck.py | 93 ++++++++++-- docs/recheck.md | 14 +- scripts/06-clinvar-screen.sh | 12 ++ scripts/run-all.sh | 11 ++ tests/test_owngenome_recheck.py | 249 ++++++++++++++++++++++++++++++++ 5 files changed, 360 insertions(+), 19 deletions(-) diff --git a/bin/owngenome_recheck.py b/bin/owngenome_recheck.py index 6dd4c69d..0430e60e 100755 --- a/bin/owngenome_recheck.py +++ b/bin/owngenome_recheck.py @@ -22,6 +22,8 @@ in-pipeline-release a newer pin is in a pipeline release: update the pipeline upstream-only upstream has a newer release that no pipeline release pins no-newer-release the lookup answered and nothing is newer + unknown a pin with no upstream lookup that the run's manifest + does not name: nothing was compared error the lookup failed (never read as up to date); exit 1 The steps an item feeds come from the reads column of assets/steps.tsv, then every step downstream of them through its upstream column, then the reports. @@ -37,7 +39,9 @@ check named, and lists the named script steps (04b, 14, 29, 03a-03d), which skip when their output exists, as not redone. --refresh-clinvar runs setup.sh --refresh clinvar first, and refuses while a run of any sample of -GENOME_DIR is going on, since they share the ClinVar files. Standard library only. +GENOME_DIR is going on, since they share the ClinVar files. It holds +GENOME_DIR/clinvar/.refresh.lock while it refreshes: run-all.sh and step 06 +wait on it. Standard library only. """ import argparse import fcntl @@ -66,8 +70,11 @@ REPORT_STEP = "24" # the HTML report; the text report is written after it # ClinVar's ##fileDate is the release; NCBI's Last-Modified comes a day or two later. CLINVAR_SLACK_DAYS = 2 -CLASSES = ("rerun-ready", "refresh-available", "in-pipeline-release", "upstream-only", "no-newer-release", "error") +CLASSES = ("rerun-ready", "refresh-available", "in-pipeline-release", "upstream-only", "no-newer-release", "unknown", + "error") LOG_HEADER = "utc\tevent\titem\tclass\tdetail\n" +# Held exclusive by --refresh-clinvar, shared by scripts/run-all.sh and scripts/06-clinvar-screen.sh. +CLINVAR_LOCK = "clinvar/.refresh.lock" CLINVAR_FILES = ("clinvar/clinvar_pathogenic_chr.vcf.gz", "clinvar/clinvar_chr.vcf.gz", "clinvar/clinvar.vcf.gz") # The pinned items: (manifest key, label, freshness lookup, how to read its upstream). PINNED = ( @@ -290,7 +297,10 @@ def pgs(self, pgs_id): except self.fr.LOOKUP_ERRORS as e: raise Failed(str(e)) import email.utils - return email.utils.parsedate_to_datetime(lm).strftime("%Y%m%d") + try: # a malformed date is a failed lookup of this score, never a crash of the whole check + return email.utils.parsedate_to_datetime(lm).strftime("%Y%m%d") + except (TypeError, ValueError, AttributeError, OverflowError): + raise Failed(f"cannot read the Last-Modified of {url}: {lm!r}") def release(self): try: @@ -334,7 +344,10 @@ def pgs(self, pgs_id): def release(self): e = self._get("release") - return e["tag"], dict(e.get("pins") or {}) + tag, pins = e.get("tag"), e.get("pins") or {} + if not isinstance(tag, str) or not tag or not isinstance(pins, dict): + raise Failed("the offline answers' release needs a tag and a pins object") + return tag, dict(pins) # --- the check ----------------------------------------------------------------------------- @@ -362,6 +375,10 @@ def classify(item, release_tag, release_pins, release_error, clinvar=False): ran, inst, up = item.get("ran_with"), item.get("installed"), item.get("upstream") if inst and ran and ver(inst) > ver(ran): return "rerun-ready" + # An image re-pinned at the same tag (a new digest): the tags compare equal, the references do not. + ran_ref, inst_ref = item.get("ran_with_ref"), item.get("installed_ref") + if inst and ran and ver(inst) == ver(ran) and ran_ref and inst_ref and ran_ref != inst_ref: + return "rerun-ready" # ClinVar's date is in the manifest whenever step 06 had a file, so a run without one cannot be compared # with the file here. A pin with no run value (HLA data not installed then) compares this checkout's. if clinvar and not ran: @@ -391,6 +408,14 @@ def classify(item, release_tag, release_pins, release_error, clinvar=False): return "upstream-only" +def skipped_steps(sample, sample_dir): + """The steps run-all.sh recorded as skipped in this run (logs/run_status.tsv).""" + try: + return {s["step"] for s in run_state.collect(sample, sample_dir)["steps"] if s["state"] == "skipped"} + except (OSError, ValueError, KeyError): + return set() + + def action(item, sample): c, key = item["class"], item["item"] if c == "rerun-ready": @@ -415,6 +440,9 @@ def action(item, sample): return "no pipeline release pins it yet; nothing to do until one does" if c == "error": return "the lookup failed; check again later" + if c == "unknown": + return ("the run's manifest does not name this pin (a run of an older version): nothing to compare; " + "the next run records it") return "" @@ -446,18 +474,31 @@ def add(item, fn, *args, clinvar=False): ran = pin_value(key, man.get(("data", key)) if key == "hla_database" else version.get(key)) inst = this_hla if key == "hla_database" else pin_value(key, this_pins.get(key)) item = {"item": key, "label": label, "kind": "pin", "ran_with": ran, "installed": inst} + if key.endswith("_IMAGE"): + item["ran_with_ref"], item["installed_ref"] = version.get(key) or None, this_pins.get(key) or None if key == "hla_database" and not this_hla: item["error"] = "could not read this checkout's HLA_DB_RELEASE (scripts/lib/common.sh)" add(item, answers.pinned, key, lookup, pattern) # The other versions.env pins a step reads (steps 50 to 53): no freshness source exists for them, so - # nothing is looked up; a pin here that differs from the run's makes its step rerun-ready. + # nothing is looked up; a pin here that differs from the run's (an image by its full reference, so a + # new digest at the same tag counts) makes its step rerun-ready. A pin the run's manifest lacks is + # unknown, never no-newer-release. A pin read only by steps run-all.sh recorded as skipped (an opt-in + # step the run left out) is left out: a rerun would not run that step either. known = {k for k, *_ in PINNED} + skipped = skipped_steps(sample, sample_dir) for key in dict.fromkeys(r[:-5] if r.endswith("@self") else r for s in steps for r in reads_of(s)): if key in known or key not in this_pins: continue - ran, inst = pin_value(key, version.get(key)), pin_value(key, this_pins.get(key)) + readers = [s["step"] for s in steps if any(key_matches(key, r) for r in reads_of(s))] + if readers and all(n in skipped for n in readers): + continue + ran_ref, inst_ref = version.get(key) or None, this_pins.get(key) or None + ran, inst = pin_value(key, ran_ref), pin_value(key, inst_ref) + cls = "unknown" if not (ran_ref and inst_ref) else "rerun-ready" if ran_ref != inst_ref else "no-newer-release" item = {"item": key, "label": f"{key} (no upstream lookup)", "kind": "pin", "ran_with": ran, - "installed": inst, "class": "rerun-ready" if ran and inst and ran != inst else "no-newer-release"} + "installed": inst, "class": cls} + if key.endswith("_IMAGE"): + item["ran_with_ref"], item["installed_ref"] = ran_ref, inst_ref item["steps"] = affected_steps(steps, [key]) item["action"] = action(item, sample) items.append(item) @@ -520,6 +561,9 @@ def print_check(state, out=None): w(f"Database re-check of {state['sample']} ({state['checked_utc']}):\n") for i in state["items"]: detail = i.get("error") or f"run used {i.get('ran_with') or 'unknown'}, here {i.get('installed') or 'unknown'}, upstream {i.get('upstream') or 'unknown'}" + if not i.get("error") and i.get("ran_with_ref") and i.get("installed_ref") and i["ran_with_ref"] != i["installed_ref"] \ + and i.get("ran_with") == i.get("installed"): + detail += f"; same tag, new image: run used {i['ran_with_ref']}, here {i['installed_ref']}" w(f" {i['label']}: {i['class']} ({detail})\n") if i["class"] not in ("no-newer-release",) and i.get("action"): w(f" {i['action']}\n") @@ -624,6 +668,31 @@ def running_samples(genome_dir): and run_state.collect(n, os.path.join(genome_dir, n))["run"]["state"] == "running"] +def refresh_clinvar(genome_dir): + """setup.sh --refresh clinvar under the ClinVar lock of GENOME_DIR (CLINVAR_LOCK). The ClinVar files + are shared by every sample, so a run of any of them blocks the refresh: run-all.sh and step 06 hold + the lock shared while they read the files, and wait while this holds it. A run with no lock (a bare + nextflow run, an older checkout) is caught by its logs/run_status.tsv. The lock is let go before + --rerun starts run-all.sh, which takes it shared.""" + os.makedirs(os.path.join(genome_dir, "clinvar"), exist_ok=True) + with open(os.path.join(genome_dir, CLINVAR_LOCK), "a") as lock: + try: + fcntl.flock(lock, fcntl.LOCK_EX | fcntl.LOCK_NB) + except OSError: + print(f"ERROR: a run or step 06 holds {os.path.join(genome_dir, CLINVAR_LOCK)}: not refreshing the " + "ClinVar files it reads.", file=sys.stderr) + return 3 + busy = running_samples(genome_dir) + if busy: + print(f"ERROR: a run of {', '.join(busy)} is going on: not refreshing the ClinVar files it reads.", + file=sys.stderr) + return 3 + r = subprocess.run([os.path.join(ROOT, "scripts", "setup.sh"), "--refresh", "clinvar", genome_dir]).returncode + if r != 0: + print("ERROR: setup.sh --refresh clinvar failed: see above.", file=sys.stderr) + return r + + def main(argv=None): ap = argparse.ArgumentParser(description=__doc__.split("\n\n")[0]) ap.add_argument("--sample") @@ -659,16 +728,8 @@ def main(argv=None): print(f"ERROR: no {sample_dir}/run_manifest.tsv: run the pipeline on {a.sample} first.", file=sys.stderr) return 2 if a.refresh_clinvar: - # The ClinVar files are shared by every sample of GENOME_DIR, so a run of any of them blocks the - # refresh. This narrows the window; a run that starts during the refresh is not stopped. - busy = running_samples(a.genome_dir) - if busy: - print(f"ERROR: a run of {', '.join(busy)} is going on: not refreshing the ClinVar files it reads.", - file=sys.stderr) - return 3 - r = subprocess.run([os.path.join(ROOT, "scripts", "setup.sh"), "--refresh", "clinvar", a.genome_dir]).returncode + r = refresh_clinvar(a.genome_dir) if r != 0: - print("ERROR: setup.sh --refresh clinvar failed: see above.", file=sys.stderr) return r fr = load_freshness() this_pins, _ = fr.read_pins(os.path.join(ROOT, "versions.env")) diff --git a/docs/recheck.md b/docs/recheck.md index e3ece402..4eae712a 100644 --- a/docs/recheck.md +++ b/docs/recheck.md @@ -23,6 +23,8 @@ Standard library only, so the plain `python3` on the host is enough. It needs th | IPD-IMGT/HLA | `hla_database` | The `release_version.txt` of the IPD-IMGT/HLA repository | 08 | | Each PGS score | `pgs:`, its `HmPOS_date` | The PGS Catalog's REST API, then the `Last-Modified` date of the harmonised GRCh38 file | 25, then 26 | +The pins of steps 50 to 53 (`MEGANE_IMAGE`, `MOCHA_PANEL_VERSION`, `MOCHA_IMAGE`, `VNTYPER_IMAGE`, `LOCITYPER_PANEL_VERSION`, `LOCITYPER_IMAGE`) have no upstream lookup: each is compared only with this checkout's `versions.env`. An image is compared by its full reference, so a new `@sha256` digest at the same tag counts as a change, for these and for `PHARMCAT_IMAGE`. A pin read only by a step that `run-all.sh` recorded as skipped in this run (an opt-in step it left out) is not listed: a rerun would not run that step either. + The HTML report (24) and the text report are written again after any of them. The ClinVar release inside the PCGR bundle (`cpsr_clinvar_release`) moves with the PCGR bundle, never with step 06's ClinVar file. The lookups are the ones the monthly CI check uses (`scripts/ci/freshness.py`). It also reads the newest release of this pipeline (`--owngenome-repo`, default `GeiserX/OwnGenome`) and that release's `versions.env`, to say whether a newer pin is already in a release. @@ -37,11 +39,12 @@ Each item gets one class: | Class | Meaning | What to do | |---|---|---| -| `rerun-ready` | What is installed here is newer than what the run used: a newer ClinVar file or PGS file on disk, or a pin in this checkout | Rerun (for a pin, run `./scripts/setup.sh` first so its data is there) | -| `refresh-available` | ClinVar only: NCBI's file is more than 2 days newer than the one installed | `./scripts/setup.sh --refresh clinvar $GENOME_DIR`, then rerun; `--rerun --refresh-clinvar` does both, and refuses while any sample of `$GENOME_DIR` is running, since they share the ClinVar files | +| `rerun-ready` | What is installed here is newer than what the run used: a newer ClinVar file or PGS file on disk, or a pin in this checkout (an image re-pinned to a new digest at the same tag too) | Rerun (for a pin, run `./scripts/setup.sh` first so its data is there) | +| `refresh-available` | ClinVar only: NCBI's file is more than 2 days newer than the one installed | `./scripts/setup.sh --refresh clinvar $GENOME_DIR`, then rerun; `--rerun --refresh-clinvar` does both, and refuses while any sample of `$GENOME_DIR` is running, since they share the ClinVar files ([the ClinVar lock](#the-clinvar-lock)) | | `in-pipeline-release` | A pipeline release pins a newer version | Update the pipeline to that release, run `setup.sh`, then rerun | | `upstream-only` | Upstream has a newer release that no pipeline release pins | Nothing yet for a pinned item; for a PGS score, download it again ([step 25](25-prs.md)) | | `no-newer-release` | The lookup answered and nothing is newer | Nothing | +| `unknown` | A pin of steps 50 to 53 that the run's manifest does not name (a run of an older version): nothing was compared | Nothing; the next run records it. It does not make the command exit 1 | | `error` | The lookup failed | Check again later. A failed lookup is never read as "up to date", and the command exits 1 | The 2-day margin for ClinVar is there because NCBI dates the file a day or two after the release in its `##fileDate`. A newer PGS file date means the catalog wrote the harmonised file again; it does not prove the weights changed. @@ -67,6 +70,10 @@ When the run ends, it compares the steps that ran a task again with the steps th A sample last run by an older `run-all.sh` has no `meta command` row: run the same command once by hand. +### The ClinVar lock + +The ClinVar files are shared by every sample of `$GENOME_DIR`. `--refresh-clinvar` holds `$GENOME_DIR/clinvar/.refresh.lock` exclusive while `setup.sh --refresh clinvar` replaces them, and lets go before the rerun starts. `run-all.sh` holds the same lock shared from its start to its end, and `scripts/06-clinvar-screen.sh` while it runs: a run that starts during a refresh waits for it ("Waiting for the ClinVar refresh ..."), and a refresh refuses (exit 3) while a run holds the lock. It also refuses while any sample's `logs/run_status.tsv` says a run is going on, which covers a run of an older checkout. The shell side needs `flock` (part of util-linux on Linux; on macOS `brew install flock`); without it nothing waits. A bare `nextflow run`, and `setup.sh --refresh clinvar` run by hand, take no lock. + ### A bare Nextflow run A bare `nextflow run main.nf` writes no run manifest, so there is nothing to check. To get the command back, `--launch-dir` prints the last command of that launch directory (its `.nextflow/history`) with `-resume` added. It never runs it. @@ -77,9 +84,10 @@ python3 bin/owngenome_recheck.py --launch-dir /path/to/launch/dir ## The files it writes -- `/recheck/state.json`, replaced whole on each check: `checked_utc`, the hosts asked, the pipeline release it read, one entry per item (`item`, `label`, `ran_with`, `installed`, `upstream`, `class`, `steps`, `action`, and `error` or `in_owngenome_release` when they apply), `affected_steps` (the steps a rerun would redo for the `rerun-ready` and `refresh-available` items), and after a rerun a `rerun` block (`started_utc`, `finished_utc`, `rc`, `command`, `predicted_vs_redone`). The last rerun stays until another one replaces it. +- `/recheck/state.json`, replaced whole on each check: `checked_utc`, the hosts asked, the pipeline release it read, one entry per item (`item`, `label`, `ran_with`, `installed`, `upstream`, `class`, `steps`, `action`, and `error` or `in_owngenome_release` when they apply; an image item also has `ran_with_ref` and `installed_ref`, the full references with any digest), `affected_steps` (the steps a rerun would redo for the `rerun-ready` and `refresh-available` items), and after a rerun a `rerun` block (`started_utc`, `finished_utc`, `rc`, `command`, `predicted_vs_redone`). The last rerun stays until another one replaces it. - `/recheck/log.tsv`, only ever appended: one row per item and check, and a row as each rerun starts and ends. - `/recheck/rerun.lock`, an empty file each `--rerun` locks while it runs. +- `$GENOME_DIR/clinvar/.refresh.lock`, an empty file: [the ClinVar lock](#the-clinvar-lock). `run-all.sh` and step 06 make it too when it is missing. ## In the reports and on the live site diff --git a/scripts/06-clinvar-screen.sh b/scripts/06-clinvar-screen.sh index 7023c861..4d938d89 100755 --- a/scripts/06-clinvar-screen.sh +++ b/scripts/06-clinvar-screen.sh @@ -61,6 +61,18 @@ done echo "=== ClinVar Pathogenic Screen: ${SAMPLE} ===" +# bin/owngenome_recheck.py --refresh-clinvar holds this lock while it replaces the ClinVar files every sample +# shares; the screen holds it shared on fd 7 until it ends, so it waits for a refresh going on and a refresh +# refuses while it reads. Without flock (macOS: brew install flock) nothing waits here. +CVLOCK="${GENOME_DIR}/clinvar/.refresh.lock" +if [ -d "${GENOME_DIR}/clinvar" ] && command -v flock >/dev/null 2>&1; then + [ -e "$CVLOCK" ] || : 2>/dev/null >> "$CVLOCK" || true + if [ -r "$CVLOCK" ]; then + exec 7<"$CVLOCK" + flock -s -n 7 || { echo "Waiting for the ClinVar refresh that holds ${CVLOCK} to end ..."; flock -s 7; } + fi +fi + for f in "$VCF" "${VCF}.tbi" "$CLINVAR" "${CLINVAR}.tbi" "$REF" "${REF}.fai" "$SITES"; do if [ ! -f "$f" ]; then echo "ERROR: File not found: ${f}" >&2 diff --git a/scripts/run-all.sh b/scripts/run-all.sh index ab6c3cdf..c5d15d59 100755 --- a/scripts/run-all.sh +++ b/scripts/run-all.sh @@ -62,6 +62,17 @@ export GENOME_DIR=${GENOME_DIR:?Set GENOME_DIR to your data directory} validate_sample "$SAMPLE" G=$GENOME_DIR S="${GENOME_DIR}/${SAMPLE}" LOG_DIR="${GENOME_DIR}/${SAMPLE}/logs" T=${TOOLS:-} C=${EXTRA_CALLERS:-} on() { [[ "${!1:-}" =~ ^(true|1)$ ]]; } +# The ClinVar files are shared by every sample of GENOME_DIR. bin/owngenome_recheck.py --refresh-clinvar holds this +# lock while it replaces them; the run holds it shared on fd 7 until it ends, so a run waits for a refresh going on +# and a refresh refuses while a run goes on. Without flock (macOS: brew install flock) nothing waits here. +CVLOCK="${G}/clinvar/.refresh.lock" +if [ -d "${G}/clinvar" ] && command -v flock >/dev/null 2>&1; then + [ -e "$CVLOCK" ] || : 2>/dev/null >> "$CVLOCK" || true + if [ -r "$CVLOCK" ]; then + exec 7<"$CVLOCK" + flock -s -n 7 || { echo "Waiting for the ClinVar refresh that holds ${CVLOCK} to end ..."; flock -s 7; } + fi +fi java_major=$(java -version 2>&1 | awk -F'"' '/version "/ {split($2, v, "."); print (v[1] == 1 ? v[2] : v[1]); exit}') || true if ! command -v nextflow >/dev/null 2>&1 || [[ ! "${java_major:-}" =~ ^[0-9]+$ ]] || [ "$java_major" -lt 17 ]; then printf '%s\n' "ERROR: run-all.sh starts the Nextflow pipeline and needs Java 17 or later (found: ${java_major:-none}) and Nextflow ($(command -v nextflow || echo 'not on PATH'))." \ diff --git a/tests/test_owngenome_recheck.py b/tests/test_owngenome_recheck.py index bb5b8a50..ea901237 100644 --- a/tests/test_owngenome_recheck.py +++ b/tests/test_owngenome_recheck.py @@ -27,6 +27,18 @@ 12. check() itself names step 53 rerun-ready when the run used an older LOCITYPER_PANEL_VERSION than this checkout pins, and not when they match; its action names setup.sh --locityper-data, which the main setup leaves out. + 13. an image re-pinned at the same tag with a new digest is rerun-ready, for a pin + with no upstream lookup (step 50's) and for a looked-up one (PharmCAT); the + same reference on both sides is no-newer-release; + 14. a pin with no upstream lookup that the run's manifest lacks is unknown, never + no-newer-release, names no step and does not make the command exit 1; + 15. a pin read only by a step run-all.sh recorded as skipped is left out; + 16. a malformed Last-Modified of a PGS file, and an offline release with no tag, + are errors of that lookup, never a crash; + 17. --refresh-clinvar refuses (exit 3) while a run holds the ClinVar lock shared, + and holds it exclusive while setup.sh refreshes; + 18. the ClinVar lock block of run-all.sh and of step 06 waits while a refresh holds + the lock and goes on once it ends (needs flock; it fails in CI without it). Each check is then run against a planted mistake and must fail there. Run: python3 tests/test_owngenome_recheck.py @@ -272,6 +284,8 @@ def all_checks(m): cmp["redone"] == ["06", "07"] and cmp["predicted_not_redone"] == ["13"] and cmp["redone_not_predicted"] == ["07"] and cmp["not_redone_script_steps"] == ["14"], cmp) + new_checks(m, tmp, steps, out) + ld = os.path.join(tmp, "launch") os.makedirs(os.path.join(ld, ".nextflow")) with open(os.path.join(ld, ".nextflow", "history"), "w") as f: @@ -290,6 +304,204 @@ def all_checks(m): return out +def manifest_add(g, rows): + with open(os.path.join(g, "S1", "run_manifest.tsv"), "a") as f: + f.writelines(f"version\t{k}\t{v}\n" for k, v in rows) + + +def skip_step(g, step): + with open(os.path.join(g, "S1", "logs", "run_status.tsv"), "a") as f: + f.write(f"step\t{step}\tskipped (opt-in: add mei to TOOLS)\n") + + +def new_checks(m, tmp, steps, out): + # Image references built from the fixture's own, so no image name is a literal here. + ref = THIS_PINS["PHARMCAT_IMAGE"] + d1, d2 = ref + "@sha256:" + "a" * 64, ref + "@sha256:" + "b" * 64 + g = genome(os.path.join(tmp, "digest")) + manifest_add(g, [("MEGANE_IMAGE", d1)]) + new = by_item(m.check("S1", g, answers(m, g), steps, dict(THIS_PINS, MEGANE_IMAGE=d2), "3.66.0")) + same = by_item(m.check("S1", g, answers(m, g), steps, dict(THIS_PINS, MEGANE_IMAGE=d1), "3.66.0")) + mg, ms = new.get("MEGANE_IMAGE", {}), same.get("MEGANE_IMAGE", {}) + out["a new digest at the same tag of a pin with no lookup is rerun-ready"] = ( + mg.get("class") == "rerun-ready" and "50" in mg.get("steps", []) and mg.get("installed_ref") == d2 + and ms.get("class") == "no-newer-release", (mg, ms)) + + gp = genome(os.path.join(tmp, "digest-pharmcat")) + rows = [x for x in open(os.path.join(gp, "S1", "run_manifest.tsv")) if not x.startswith("version\tPHARMCAT_IMAGE\t")] + open(os.path.join(gp, "S1", "run_manifest.tsv"), "w").writelines(rows + [f"version\tPHARMCAT_IMAGE\t{d1}\n"]) + ph = by_item(m.check("S1", gp, answers(m, gp), steps, dict(THIS_PINS, PHARMCAT_IMAGE=d2), "3.66.0"))["PHARMCAT_IMAGE"] + ph1 = by_item(m.check("S1", gp, answers(m, gp), steps, dict(THIS_PINS, PHARMCAT_IMAGE=d1), "3.66.0"))["PHARMCAT_IMAGE"] + text = io.StringIO() + m.print_check({"sample": "S1", "checked_utc": "x", "items": [ph], "affected_steps": []}, text) + out["a new digest at the same tag of a looked-up image is rerun-ready"] = ( + ph["class"] == "rerun-ready" and ph1["class"] == "no-newer-release" and "same tag, new image" in text.getvalue(), + (ph, ph1, text.getvalue())) + + gu = genome(os.path.join(tmp, "unknown")) + st = m.check("S1", gu, answers(m, gu), steps, dict(THIS_PINS, MEGANE_IMAGE=d1), "3.66.0") + mu = by_item(st).get("MEGANE_IMAGE", {}) + with contextlib.redirect_stdout(io.StringIO()): + ans = write(tmp, "unknown.json", json.load(open(os.path.join(FIX, "offline.json")))) + keep = m.check + m.check = lambda s_, g_, a_, st_, pins, hla: keep(s_, g_, a_, st_, dict(pins, MEGANE_IMAGE=d1), hla) + try: + rc = m.main(["--sample", "S1", "--genome-dir", gu, "--offline", ans]) + finally: + m.check = keep + out["a pin with no lookup the run's manifest lacks is unknown, exit 0"] = ( + mu.get("class") == "unknown" and "50" not in st["affected_steps"] and "50" not in st["rerun_ready_steps"] + and "older version" in mu.get("action", "") and rc == 0, (mu, st["affected_steps"], rc)) + + gs = genome(os.path.join(tmp, "skipped")) + manifest_add(gs, [("MEGANE_IMAGE", d1)]) + skip_step(gs, "50") + sk = m.check("S1", gs, answers(m, gs), steps, dict(THIS_PINS, MEGANE_IMAGE=d2), "3.66.0") + out["a pin read only by a step the run skipped is left out"] = ( + "MEGANE_IMAGE" not in by_item(sk) and "50" not in sk["rerun_ready_steps"], by_item(sk).get("MEGANE_IMAGE")) + + srv, port = serve_pgs("not a date") + try: + fr = m.load_freshness() + fr.RETRIES = 0 + on = m.Online(fr, {}, "o/r", pgs_api=f"http://127.0.0.1:{port}/rest/score", raw=f"http://127.0.0.1:{port}", + out=io.StringIO()) + try: + gm = genome(os.path.join(tmp, "lastmod")) + pg = by_item(m.check("S1", gm, Mixed(on, answers(m, gm)), steps, THIS_PINS, "3.66.0"))["pgs:PGS000001"] + got = (pg["class"] == "error" and "Last-Modified" in pg.get("error", ""), pg) + except Exception as e: # noqa: BLE001 the crash this check is about + got = (False, repr(e)) + finally: + srv.shutdown() + out["a malformed Last-Modified of a PGS file is error, not a crash"] = got + + try: + go = genome(os.path.join(tmp, "notag")) + so = m.check("S1", go, answers(m, go, release={"pins": {}}), steps, THIS_PINS, "3.66.0") + got = ("error" in so["pipeline_release"] and "tag" in so["pipeline_release"]["error"], so["pipeline_release"]) + except Exception as e: # noqa: BLE001 + got = (False, repr(e)) + out["an offline release with no tag is error, not a crash"] = got + + gl = genome(os.path.join(tmp, "cvlock")) + lock = os.path.join(gl, "clinvar", ".refresh.lock") + keep_sp, ran = m.subprocess, [] + m.subprocess = type("SP", (), {"run": staticmethod(lambda *a, **k: ran.append(a) or type("R", (), {"returncode": 0})())}) + e = io.StringIO() + try: + with open(lock, "a") as held: + fcntl.flock(held, fcntl.LOCK_SH | fcntl.LOCK_NB) + with contextlib.redirect_stdout(io.StringIO()), contextlib.redirect_stderr(e): + rc = m.main(["--sample", "S1", "--genome-dir", gl, "--offline", os.path.join(FIX, "offline.json"), + "--rerun", "--refresh-clinvar"]) + finally: + m.subprocess = keep_sp + out["--refresh-clinvar refuses while a run holds the ClinVar lock"] = ( + rc == 3 and ".refresh.lock" in e.getvalue() and not ran, (rc, e.getvalue(), ran)) + + seen = [] + + def setup_run(argv, *a, **k): + if "--refresh" in argv: + with open(lock, "a") as f: + try: + fcntl.flock(f, fcntl.LOCK_SH | fcntl.LOCK_NB) + seen.append("free") + except OSError: + seen.append("held") + return type("R", (), {"returncode": 0})() + m.subprocess = type("SP", (), {"run": staticmethod(setup_run)}) + try: + with contextlib.redirect_stdout(io.StringIO()), contextlib.redirect_stderr(io.StringIO()): + rc = m.main(["--sample", "S1", "--genome-dir", gl, "--offline", os.path.join(FIX, "offline.json"), + "--rerun", "--refresh-clinvar"]) + with open(lock, "a") as f: # let go before the replay's run-all.sh takes it shared + fcntl.flock(f, fcntl.LOCK_SH | fcntl.LOCK_NB) + after = True + except OSError: + after = False + finally: + m.subprocess = keep_sp + out["--refresh-clinvar holds the ClinVar lock while setup.sh refreshes"] = ( + seen == ["held"] and after and rc == 0, (seen, after, rc)) + + +def lock_block(path): + """The ClinVar lock block of a script: from 'CVLOCK=' through its closing 'fi'.""" + x = re.search(r"^CVLOCK=.*?\n fi\nfi\n", open(path).read(), re.S | re.M) + return x.group(0) if x else "" + + +def lock_waits(block): + """block, run in bash while this process holds the lock exclusive: it must wait, then go on + once the lock is let go.""" + import subprocess + import time + tmp = tempfile.mkdtemp(prefix="recheck-cvlock.") + try: + g = os.path.join(tmp, "genome") + os.makedirs(os.path.join(g, "clinvar")) + lock = os.path.join(g, "clinvar", ".refresh.lock") + with open(lock, "a") as held: + fcntl.flock(held, fcntl.LOCK_EX | fcntl.LOCK_NB) + p = subprocess.Popen(["bash", "-c", "set -euo pipefail\n" + block + "echo through\n"], + env=dict(os.environ, G=g, GENOME_DIR=g), stdout=subprocess.PIPE, text=True) + time.sleep(1.5) + waited = p.poll() is None + try: + o, _ = p.communicate(timeout=20) + except subprocess.TimeoutExpired: + p.kill() + o, _ = p.communicate() + return (bool(block) and waited and p.returncode == 0 and "Waiting" in o and o.rstrip().endswith("through"), + (waited, p.returncode, o)) + finally: + shutil.rmtree(tmp, ignore_errors=True) + + +def serve_pgs(last_modified): + """A PGS Catalog score answer, and a HEAD of its file with this Last-Modified.""" + class H(http.server.BaseHTTPRequestHandler): + def do_GET(self): + port = self.server.server_address[1] + body = json.dumps({"id": "PGS000001", "ftp_harmonized_scoring_files": {"GRCh38": { + "positions": f"http://127.0.0.1:{port}/f/PGS000001_hmPOS_GRCh38.txt.gz"}}}).encode() + self.send_response(200) + self.send_header("Content-Type", "application/json") + self.send_header("Content-Length", str(len(body))) + self.end_headers() + self.wfile.write(body) + + def do_HEAD(self): + self.send_response(200) + self.send_header("Last-Modified", last_modified) + self.send_header("Content-Length", "0") + self.end_headers() + + def log_message(self, *a): + pass + srv = http.server.ThreadingHTTPServer(("127.0.0.1", 0), H) + threading.Thread(target=srv.serve_forever, daemon=True).start() + return srv, srv.server_address[1] + + +def from_source(old, new): + """bin/owngenome_recheck.py with old replaced by new (which must be in it), as a module.""" + import types + src = open(TOOL).read() + if src.count(old) != 1: + raise AssertionError(f"planted mistake not in the copy: {old[:60]!r}") + mod = types.ModuleType("owngenome_recheck_planted") + mod.__file__ = TOOL + keep, sys.dont_write_bytecode = sys.dont_write_bytecode, True + try: + exec(compile(src.replace(old, new), TOOL, "exec"), mod.__dict__) + finally: + sys.dont_write_bytecode = keep + return mod + + class Mixed: """The online release and PGS lookups, the offline answers for the rest.""" @@ -428,6 +640,30 @@ def regex_resume(d): return cmd if re.search(r"(^|\s)-resume(\s|$)", cmd) else cmd + " -resume" m.launch_dir_command = regex_resume ms["--launch-dir: an argument holding -resume is not the option"] = m + ms["a new digest at the same tag of a pin with no lookup is rerun-ready"] = from_source( + '"rerun-ready" if ran_ref != inst_ref else', '"rerun-ready" if ran != inst else') + ms["a new digest at the same tag of a looked-up image is rerun-ready"] = from_source( + "and ran_ref != inst_ref:\n return \"rerun-ready\"", "and False:\n return \"rerun-ready\"") + ms["a pin with no lookup the run's manifest lacks is unknown, exit 0"] = from_source( + 'cls = "unknown" if not', 'cls = "no-newer-release" if not') + m = load() + m.skipped_steps = lambda *a: set() + ms["a pin read only by a step the run skipped is left out"] = m + ms["a malformed Last-Modified of a PGS file is error, not a crash"] = from_source( + " try: # a malformed date is a failed lookup of this score, never a crash of the whole check\n" + " return email.utils.parsedate_to_datetime(lm).strftime(\"%Y%m%d\")\n" + " except (TypeError, ValueError, AttributeError, OverflowError):\n", + " return email.utils.parsedate_to_datetime(lm).strftime(\"%Y%m%d\")\n if False:\n") + ms["an offline release with no tag is error, not a crash"] = from_source( + 'tag, pins = e.get("tag"), e.get("pins") or {}\n if not isinstance(tag, str) or not tag or not isinstance(pins, dict):', + 'tag, pins = e["tag"], e.get("pins") or {}\n if not isinstance(pins, dict):') + m = load() + m.fcntl = type("F", (), {"LOCK_EX": 0, "LOCK_NB": 0, "flock": staticmethod(lambda *a: None)}) + ms["--refresh-clinvar refuses while a run holds the ClinVar lock"] = m + ms["--refresh-clinvar holds the ClinVar lock while setup.sh refreshes"] = from_source( + " fcntl.flock(lock, fcntl.LOCK_EX | fcntl.LOCK_NB)\n except OSError:\n" + " print(f\"ERROR: a run or step 06", + " pass\n except OSError:\n print(f\"ERROR: a run or step 06") return ms @@ -595,6 +831,19 @@ def main(): unfiltered = plugin_from_source(src.replace(stale, "if False:"), "owngenome_live_recheck_unfiltered") check("fails on a planted mistake: the live-site plugin leaves out a check older than the run's manifest", not plugin_stale_ok(unfiltered)[0], "passed on the mutant") + import shutil as _sh + for script in ("scripts/run-all.sh", "scripts/06-clinvar-screen.sh"): + block = lock_block(os.path.join(REPO, script)) + if not _sh.which("flock"): + check(f"the ClinVar lock block of {script} waits for a refresh (flock is not installed: " + f"{'FAIL in CI' if os.environ.get('CI') else 'not run here'})", not os.environ.get("CI")) + continue + ok, detail = lock_waits(block) + check(f"the ClinVar lock block of {script} waits for a refresh, then goes on", ok, detail) + planted = re.sub(r"(?m)^ flock -s -n 7 .*$", " :", block) + check(f"planted mistake is in the copy: the ClinVar lock block of {script}", planted != block) + check(f"fails on a planted mistake: the ClinVar lock block of {script} waits for a refresh (no flock call)", + not lock_waits(planted)[0], "passed on the mutant") results = all_checks(load()) for name, (ok, detail) in results.items(): check(name, ok, detail) From ae61a3455e31456cfa697c31e6953594ed930bf0 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Fri, 6 Nov 2026 22:23:31 +0100 Subject: [PATCH 4/8] build(renovate): track the nf-weblog plugin pin of the live site NF_WEBLOG_PLUGIN in scripts/run-all.sh had no Renovate manager, so no bot saw a new nf-weblog release. The rule also moves the copy in docs/live-site.md. Bead: pgp-5bja.23 --- renovate.json | 8 ++++++++ scripts/run-all.sh | 2 +- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/renovate.json b/renovate.json index 4cf6b687..76858499 100644 --- a/renovate.json +++ b/renovate.json @@ -63,6 +63,14 @@ "extractVersionTemplate": "^v(?.+)$", "versioningTemplate": "regex:^(?\\d+)\\.(?\\d+)\\.(?\\d+)$" }, + { + "customType": "regex", + "description": "nf-weblog, the Nextflow plugin of the live site: NF_WEBLOG_PLUGIN in scripts/run-all.sh and the plugins block docs/live-site.md shows, in one PR. Its releases are tagged without a v.", + "managerFilePatterns": ["/^scripts/run-all\\.sh$/", "/^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 c5d15d59..9988346d 100755 --- a/scripts/run-all.sh +++ b/scripts/run-all.sh @@ -306,7 +306,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. +# without it neither downloads the plugin nor warns about an unknown weblog scope. Renovate bumps the pin (renovate.json). NF_WEBLOG_PLUGIN=nf-weblog@1.2.0 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. From feb1f11d474069ba39b02930e8f544b0f7235809 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Fri, 6 Nov 2026 22:28:59 +0100 Subject: [PATCH 5/8] fix(sv): pass a large integral RDCN to the depth rule as digits A Delly RDCN of 1000000 was formatted as 1e+06, which the shared rule reads as no estimate. Adds chrY, PAR deletion and PAR-boundary cases. Bead: pgp-5bja.24 --- bin/annotsv_postfilter.py | 3 ++- tests/test_annotsv_postfilter.py | 29 +++++++++++++++++++++-------- 2 files changed, 23 insertions(+), 9 deletions(-) diff --git a/bin/annotsv_postfilter.py b/bin/annotsv_postfilter.py index 02f6f355..004a021a 100755 --- a/bin/annotsv_postfilter.py +++ b/bin/annotsv_postfilter.py @@ -266,7 +266,8 @@ def depth_support(row, caller, sex="", par=None): start = int(row["SV_start"]) except (KeyError, ValueError): start = 0 - no = cnvpytor_to_vcf.delly_depth_contradicts(t, chrom, start, start + length, ["RDCN"], [f"{rdcn:g}"], sex, + cn = str(int(rdcn)) if rdcn.is_integer() else f"{rdcn:g}" + no = cnvpytor_to_vcf.delly_depth_contradicts(t, chrom, start, start + length, ["RDCN"], [cn], sex, family_report.in_par(par or {}, chrom, start)) return ("no" if no else "yes"), f"Delly RDCN {rdcn:g}" if caller == "delly" and rdcn is None and length > DELLY_DEPTH_BP: diff --git a/tests/test_annotsv_postfilter.py b/tests/test_annotsv_postfilter.py index 0f9b9fe1..4335f122 100755 --- a/tests/test_annotsv_postfilter.py +++ b/tests/test_annotsv_postfilter.py @@ -165,6 +165,14 @@ def run7(work, name, with_rules=True): row7("XDEL0", "X", 51000000, 51050000, "DEL", 4, "delly", "G21", "GT:RDCN", "1/1:0"), row7("XDUP2", "X", 52000000, 52050000, "DUP", 4, "delly", "G22", "GT:RDCN", "0/1:2"), row7("PARDUP2", "X", 1000000, 1050000, "DUP", 4, "delly", "G23", "GT:RDCN", "0/1:2"), + row7("PARDEL1", "X", 1000000, 1050000, "DEL", 4, "delly", "G24", "GT:RDCN", "0/1:1"), + row7("YDEL1", "Y", 20000000, 20050000, "DEL", 4, "delly", "G25", "GT:RDCN", "1/1:1"), + row7("YDUP2", "Y", 21000000, 21050000, "DUP", 4, "delly", "G26", "GT:RDCN", "0/1:2"), + # PAR1 ends at chrX:2,781,479 (assets/par_grch38.bed): the last base in + # it, then the first base past it. + row7("PAREDGE_IN", "X", 2781479, 2831479, "DUP", 4, "delly", "G27", "GT:RDCN", "0/1:2"), + row7("PAREDGE_OUT", "X", 2781480, 2831480, "DUP", 4, "delly", "G28", "GT:RDCN", "0/1:2"), + row7("BIGDUP", "20", 48000000, 48050000, "DUP", 4, "delly", "G29", "GT:RDCN", "0/1:1000000"), ]) @@ -232,16 +240,21 @@ def col(i, c): # Delly RDCN on chrX and chrY by sex and PAR (--sex, --par-bed), the rule # cnvpytor_to_vcf.py applies when step 05 adds a Delly call. par = os.path.join(REPO, "assets", "par_grch38.bed") - for sex, want_sex in (("male", {"XDEL1": ("no", "3"), "XDUP2": ("yes", "4"), "PARDUP2": ("no", "3"), - "XDEL0": ("yes", "4")}), - ("female", {"XDEL1": ("yes", "4"), "XDUP2": ("no", "3"), "PARDUP2": ("no", "3"), - "XDEL0": ("yes", "4")}), - ("", {"XDEL1": ("yes", "4"), "XDUP2": ("yes", "4"), "PARDUP2": ("no", "3"), - "XDEL0": ("yes", "4")})): + yes, no = ("yes", "4"), ("no", "3") + for sex, want_sex in (("male", {"XDEL1": no, "XDUP2": yes, "PARDUP2": no, "XDEL0": yes, "PARDEL1": yes, + "YDEL1": no, "YDUP2": yes, "PAREDGE_IN": no, "PAREDGE_OUT": yes, + "BIGDUP": yes}), + ("female", {"XDEL1": yes, "XDUP2": no, "PARDUP2": no, "XDEL0": yes, "PARDEL1": yes, + "YDEL1": yes, "YDUP2": no, "PAREDGE_IN": no, "PAREDGE_OUT": no, + "BIGDUP": yes}), + ("", {"XDEL1": yes, "XDUP2": yes, "PARDUP2": no, "XDEL0": yes, "PARDEL1": yes, + "YDEL1": yes, "YDUP2": yes, "PAREDGE_IN": no, "PAREDGE_OUT": yes, + "BIGDUP": yes})): sx, _ = run_sex(work, "sex_" + (sex or "none"), sex, par) got_sex = {i: (sx.get(i, {}).get("Depth_support"), sx.get(i, {}).get("ACMG_class")) for i in want_sex} - check(f"sex {sex or 'not given'}: chrX RDCN against 1 copy (male), 2 (female), either (not given); " - f"inside the PAR 2 copies", got_sex == want_sex, got_sex) + check(f"sex {sex or 'not given'}: chrX and chrY RDCN against 1 copy (male), 2 (female), either (not " + f"given); inside the PAR, up to its last base, 2 copies; RDCN 1000000 is a gain", + got_sex == want_sex, got_sex) moved = sorted(i for i, r in got.items() if r["ACMG_class"] != r["AnnotSV_class"]) want = sorted(["TR_IN", "IGH_IN", "IGH_MOST", "MHC_IN", "ALIAS", "D_LOW", "DUP_LOW", "DELLY2"]) From a154ca74ae8e5342c2a1b995c2cb3be12108dca6 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Fri, 6 Nov 2026 23:00:17 +0100 Subject: [PATCH 6/8] fix(sv): let only a PASS, non-reference Manta call hide a Delly or CNVpytor call A Manta DEL or DUP that Manta filtered (MinGQ and the rest) or genotyped 0/0 or ./. went into the overlap index, so it dropped the Delly PASS or CNVpytor deletion it overlapped, and the summary later skips the 0/0 record: the real deletion was lost. Index only base records with FILTER PASS or . and an allele other than 0; a sites-only base still covers. --- bin/cnvpytor_to_vcf.py | 23 ++++++++++++++++++++++- docs/05-annotsv.md | 2 +- tests/test_annotsv_sv_inputs.py | 21 +++++++++++++++++++++ 3 files changed, 44 insertions(+), 2 deletions(-) diff --git a/bin/cnvpytor_to_vcf.py b/bin/cnvpytor_to_vcf.py index 6747c720..45c5f20d 100755 --- a/bin/cnvpytor_to_vcf.py +++ b/bin/cnvpytor_to_vcf.py @@ -25,6 +25,10 @@ on the same chromosome overlaps it reciprocally by at least --min-overlap (default 0.5) of the longer of the two. Delly's calls are added before CNVpytor's, so a deletion both of them call and Manta misses is added once. +Only a base record with FILTER PASS (or .) and an allele other than 0 in its +first sample covers a call: a filtered Manta record (MinGQ and the rest) or +one genotyped 0/0 or ./. stays in the file but hides nothing. A sites-only +base covers with every DEL and DUP. Every record says where it came from in INFO SVCALLER (the base's records get --base-name, default manta). An added record keeps the caller's ID, END, @@ -155,6 +159,23 @@ def first_gt(fields): return values[i] if i < len(values) and values[i] else "./." +def base_call_holds(fields, n_samples): + """True when a base record can stand for a call another caller makes. + + A record Manta filtered (MinGQ, MaxDepth and the rest) or genotyped as + reference (0/0) or no-call (./.) is not a variant this sample carries; + letting it cover a Delly or CNVpytor call would drop that call, and the + summary (collect_summary.py) skips a 0/0 record, so the deletion would be + lost. A sites-only base (no sample column) has no genotype and counts as + a call. + """ + if fields[6] not in ("PASS", "."): + return False + if not n_samples: + return True + return any(a not in ("0", ".", "") for a in first_gt(fields).replace("|", "/").split("/")) + + def as_float(value): try: return float(value) @@ -214,7 +235,7 @@ def main(argv=None): for f in base: info = info_dict(f[7]) kind, sp = sv_type(f, info), span(f, info) - if kind in ("DEL", "DUP") and sp: + if kind in ("DEL", "DUP") and sp and base_call_holds(f, n_samples): known.add(f[0], kind, sp[0], sp[1]) f[7] = ("" if f[7] in ("", ".") else f[7] + ";") + "SVCALLER=" + args.base_name diff --git a/docs/05-annotsv.md b/docs/05-annotsv.md index 406d77c4..3041eb38 100644 --- a/docs/05-annotsv.md +++ b/docs/05-annotsv.md @@ -1,7 +1,7 @@ # Step 5: Structural Variant Annotation (AnnotSV) ## What This Does -Classifies every structural variant from Manta that passes duphold's depth filter (step 15; when that file is missing or older than Manta's calls, the script annotates Manta's calls and says so), plus the Delly and CNVpytor deletions and duplications that no Manta call covers, using ACMG guidelines (class 1-5), adding gene overlap, population frequency, and clinical significance. +Classifies every structural variant from Manta that passes duphold's depth filter (step 15; when that file is missing or older than Manta's calls, the script annotates Manta's calls and says so), plus the Delly and CNVpytor deletions and duplications that no Manta call covers (a Manta record that is filtered or genotyped 0/0 covers nothing), using ACMG guidelines (class 1-5), adding gene overlap, population frequency, and clinical significance. ## Why Raw Manta output contains thousands of SVs with no clinical interpretation. AnnotSV tells you which ones matter by cross-referencing known pathogenic SVs, gene databases, and population data. diff --git a/tests/test_annotsv_sv_inputs.py b/tests/test_annotsv_sv_inputs.py index 3f4693bd..ba4ee45c 100755 --- a/tests/test_annotsv_sv_inputs.py +++ b/tests/test_annotsv_sv_inputs.py @@ -306,6 +306,27 @@ def test_extra_callers(w): "delly: 2 PASS DEL/DUP calls, 1 added" in err and "cnvpytor: 3 DEL/DUP calls with e-val < 0.01 and q0 < 0.5, 2 added" in err, err) + # A Manta record that is filtered or genotyped 0/0 is no call: it must not + # hide the Delly PASS deletion it overlaps. + for label, gt, filt in (("a 0/0", "0/0", "PASS"), ("a MinGQ", "0/1", "MinGQ"), ("a ./.", "./.", "PASS")): + b = os.path.join(w, "base_nocall.vcf") + with open(b, "w") as f: + f.write(BASE + "chr20\t2000100\tMantaDEL:9\tN\t\t500\t{0}\tEND=2049900;SVTYPE=DEL;SVLEN=-49800\tGT\t{1}\n" + .format(filt, gt)) + rc, out, err = run("cnvpytor_to_vcf.py", "--base", b, "--delly", delly) + ids = [r[2] for r in records(out)] + check("{0} Manta deletion stays and does not hide Delly's DEL00000002".format(label), + rc == 0 and "MantaDEL:9" in ids and "DEL00000002" in ids, (ids, err)) + with open(b, "w") as f: + f.write(BASE + "chr20\t2000100\tMantaDEL:9\tN\t\t500\tPASS\tEND=2049900;SVTYPE=DEL;SVLEN=-49800\tGT\t0|1\n") + rc, out, err = run("cnvpytor_to_vcf.py", "--base", b, "--delly", delly) + check("control: a PASS 0|1 Manta deletion still covers Delly's DEL00000002", + rc == 0 and "DEL00000002" not in [r[2] for r in records(out)], err) + sites = "".join("\t".join(l.split("\t")[:8]) + "\n" for l in (BASE + "chr20\t2000100\tMantaDEL:9\tN\t\t500\tPASS\tEND=2049900;SVTYPE=DEL\tGT\t0/0\n").splitlines()) + rc, out, err = run("cnvpytor_to_vcf.py", "--base", "-", "--delly", delly, stdin=sites) + check("a sites-only base deletion covers Delly's DEL00000002 (no genotype to read)", + rc == 0 and "DEL00000002" not in [r[2] for r in records(out)], (records(out), err)) + rc, out, err = run("cnvpytor_to_vcf.py", "--base", "-", "--base-name", "sv_vcf", stdin=BASE) check("no other caller: stdin in, the base out with SVCALLER=sv_vcf", rc == 0 and "none given" in err and [info(r).get("SVCALLER") for r in records(out)] == ["sv_vcf"] * 2, err) From acf47c99de0f426f770b65fd4389dc38722aeb6d Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Fri, 6 Nov 2026 23:03:56 +0100 Subject: [PATCH 7/8] fix(recheck): leave out the pins of opt-in steps a default run never planned run-all.sh writes no status row for an opt-in step 50-53 left out of TOOLS, so the check saw it as 'not planned', not 'skipped', and still called its changed pin rerun-ready. Count both states, and leave nothing out when the status file records no step at all. --- bin/owngenome_recheck.py | 19 +++++++++++++------ docs/recheck.md | 2 +- tests/test_owngenome_recheck.py | 30 ++++++++++++++++++++++++------ 3 files changed, 38 insertions(+), 13 deletions(-) diff --git a/bin/owngenome_recheck.py b/bin/owngenome_recheck.py index 0430e60e..5e940be7 100755 --- a/bin/owngenome_recheck.py +++ b/bin/owngenome_recheck.py @@ -408,10 +408,17 @@ def classify(item, release_tag, release_pins, release_error, clinvar=False): return "upstream-only" -def skipped_steps(sample, sample_dir): - """The steps run-all.sh recorded as skipped in this run (logs/run_status.tsv).""" +def steps_not_run(sample, sample_dir): + """The steps this run did not run: the ones run-all.sh recorded as skipped (no BAM, missing + data), and the ones it wrote no row for, which is how it leaves out an opt-in step not in TOOLS + (run_state's "not planned"). Empty when logs/run_status.tsv records no step at all (a run of an + older version), so nothing is left out on a guess.""" try: - return {s["step"] for s in run_state.collect(sample, sample_dir)["steps"] if s["state"] == "skipped"} + status = run_state.read_run_status(os.path.join(sample_dir, "logs", "run_status.tsv")) + if not (status or {}).get("steps"): + return set() + return {s["step"] for s in run_state.collect(sample, sample_dir)["steps"] + if s["state"] in ("skipped", "not planned")} except (OSError, ValueError, KeyError): return set() @@ -482,10 +489,10 @@ def add(item, fn, *args, clinvar=False): # The other versions.env pins a step reads (steps 50 to 53): no freshness source exists for them, so # nothing is looked up; a pin here that differs from the run's (an image by its full reference, so a # new digest at the same tag counts) makes its step rerun-ready. A pin the run's manifest lacks is - # unknown, never no-newer-release. A pin read only by steps run-all.sh recorded as skipped (an opt-in - # step the run left out) is left out: a rerun would not run that step either. + # unknown, never no-newer-release. A pin read only by steps this run did not run (an opt-in step not + # in TOOLS, or one skipped for missing data) is left out: a rerun would not run that step either. known = {k for k, *_ in PINNED} - skipped = skipped_steps(sample, sample_dir) + skipped = steps_not_run(sample, sample_dir) for key in dict.fromkeys(r[:-5] if r.endswith("@self") else r for s in steps for r in reads_of(s)): if key in known or key not in this_pins: continue diff --git a/docs/recheck.md b/docs/recheck.md index 4eae712a..221c6fc0 100644 --- a/docs/recheck.md +++ b/docs/recheck.md @@ -23,7 +23,7 @@ Standard library only, so the plain `python3` on the host is enough. It needs th | IPD-IMGT/HLA | `hla_database` | The `release_version.txt` of the IPD-IMGT/HLA repository | 08 | | Each PGS score | `pgs:`, its `HmPOS_date` | The PGS Catalog's REST API, then the `Last-Modified` date of the harmonised GRCh38 file | 25, then 26 | -The pins of steps 50 to 53 (`MEGANE_IMAGE`, `MOCHA_PANEL_VERSION`, `MOCHA_IMAGE`, `VNTYPER_IMAGE`, `LOCITYPER_PANEL_VERSION`, `LOCITYPER_IMAGE`) have no upstream lookup: each is compared only with this checkout's `versions.env`. An image is compared by its full reference, so a new `@sha256` digest at the same tag counts as a change, for these and for `PHARMCAT_IMAGE`. A pin read only by a step that `run-all.sh` recorded as skipped in this run (an opt-in step it left out) is not listed: a rerun would not run that step either. +The pins of steps 50 to 53 (`MEGANE_IMAGE`, `MOCHA_PANEL_VERSION`, `MOCHA_IMAGE`, `VNTYPER_IMAGE`, `LOCITYPER_PANEL_VERSION`, `LOCITYPER_IMAGE`) have no upstream lookup: each is compared only with this checkout's `versions.env`. An image is compared by its full reference, so a new `@sha256` digest at the same tag counts as a change, for these and for `PHARMCAT_IMAGE`. A pin read only by steps this run did not run (an opt-in step not in `TOOLS`, or one skipped for missing data) is not listed: a rerun would not run that step either. The HTML report (24) and the text report are written again after any of them. The ClinVar release inside the PCGR bundle (`cpsr_clinvar_release`) moves with the PCGR bundle, never with step 06's ClinVar file. diff --git a/tests/test_owngenome_recheck.py b/tests/test_owngenome_recheck.py index ea901237..bb98fcb7 100644 --- a/tests/test_owngenome_recheck.py +++ b/tests/test_owngenome_recheck.py @@ -32,7 +32,9 @@ same reference on both sides is no-newer-release; 14. a pin with no upstream lookup that the run's manifest lacks is unknown, never no-newer-release, names no step and does not make the command exit 1; - 15. a pin read only by a step run-all.sh recorded as skipped is left out; + 15. a pin read only by a step the run did not run is left out: an opt-in step not in + TOOLS (run-all.sh writes no row for it) and one recorded as skipped (no BAM); + a status file with no step rows at all leaves every pin in; 16. a malformed Last-Modified of a PGS file, and an offline release with no tag, are errors of that lookup, never a crash; 17. --refresh-clinvar refuses (exit 3) while a run holds the ClinVar lock shared, @@ -309,9 +311,10 @@ def manifest_add(g, rows): f.writelines(f"version\t{k}\t{v}\n" for k, v in rows) -def skip_step(g, step): +def step_rows(g, rows): + """run_status.tsv step rows as run-all.sh writes them; a step it left out gets no row.""" with open(os.path.join(g, "S1", "logs", "run_status.tsv"), "a") as f: - f.write(f"step\t{step}\tskipped (opt-in: add mei to TOOLS)\n") + f.writelines(f"step\t{k}\t{v}\n" for k, v in rows) def new_checks(m, tmp, steps, out): @@ -353,12 +356,25 @@ def new_checks(m, tmp, steps, out): mu.get("class") == "unknown" and "50" not in st["affected_steps"] and "50" not in st["rerun_ready_steps"] and "older version" in mu.get("action", "") and rc == 0, (mu, st["affected_steps"], rc)) + # A default run (no TOOLS): run-all.sh writes rows for the steps it ran and none for step 50. + gd = genome(os.path.join(tmp, "not-in-tools")) + manifest_add(gd, [("MEGANE_IMAGE", d1)]) + step_rows(gd, [("06", "planned"), ("06", "ok")]) + nd = m.check("S1", gd, answers(m, gd), steps, dict(THIS_PINS, MEGANE_IMAGE=d2), "3.66.0") + out["a pin read only by an opt-in step not in TOOLS is left out"] = ( + "MEGANE_IMAGE" not in by_item(nd) and "50" not in nd["rerun_ready_steps"], by_item(nd).get("MEGANE_IMAGE")) + gs = genome(os.path.join(tmp, "skipped")) manifest_add(gs, [("MEGANE_IMAGE", d1)]) - skip_step(gs, "50") + step_rows(gs, [("06", "ok"), ("50", "skipped (no BAM)")]) sk = m.check("S1", gs, answers(m, gs), steps, dict(THIS_PINS, MEGANE_IMAGE=d2), "3.66.0") + gr = genome(os.path.join(tmp, "ran")) + manifest_add(gr, [("MEGANE_IMAGE", d1)]) + step_rows(gr, [("06", "ok"), ("50", "planned"), ("50", "ok")]) + rn = by_item(m.check("S1", gr, answers(m, gr), steps, dict(THIS_PINS, MEGANE_IMAGE=d2), "3.66.0")) out["a pin read only by a step the run skipped is left out"] = ( - "MEGANE_IMAGE" not in by_item(sk) and "50" not in sk["rerun_ready_steps"], by_item(sk).get("MEGANE_IMAGE")) + "MEGANE_IMAGE" not in by_item(sk) and "50" not in sk["rerun_ready_steps"] + and rn.get("MEGANE_IMAGE", {}).get("class") == "rerun-ready", (by_item(sk).get("MEGANE_IMAGE"), rn.get("MEGANE_IMAGE"))) srv, port = serve_pgs("not a date") try: @@ -646,8 +662,10 @@ def regex_resume(d): "and ran_ref != inst_ref:\n return \"rerun-ready\"", "and False:\n return \"rerun-ready\"") ms["a pin with no lookup the run's manifest lacks is unknown, exit 0"] = from_source( 'cls = "unknown" if not', 'cls = "no-newer-release" if not') + ms["a pin read only by an opt-in step not in TOOLS is left out"] = from_source( + 'if s["state"] in ("skipped", "not planned")}', 'if s["state"] in ("skipped",)}') m = load() - m.skipped_steps = lambda *a: set() + m.steps_not_run = lambda *a: set() ms["a pin read only by a step the run skipped is left out"] = m ms["a malformed Last-Modified of a PGS file is error, not a crash"] = from_source( " try: # a malformed date is a failed lookup of this score, never a crash of the whole check\n" From 9a3b0e32d474794e4c6c0dab24350bf20ccc71a6 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sat, 7 Nov 2026 02:06:10 +0100 Subject: [PATCH 8/8] fix(sv): keep a Delly 0/0 or ./. call out of AnnotSV's input A Delly PASS deletion genotyped 0/0 or ./. was added to the AnnotSV input and then covered the CNVpytor deletion at the same place, which the summary lost when it skipped the 0/0 record. Delly calls now need a non-reference allele too, the same rule the Manta base follows. --- bin/cnvpytor_to_vcf.py | 16 +++++++++------- docs/05-annotsv.md | 2 +- tests/test_annotsv_sv_inputs.py | 12 ++++++++++++ 3 files changed, 22 insertions(+), 8 deletions(-) diff --git a/bin/cnvpytor_to_vcf.py b/bin/cnvpytor_to_vcf.py index 45c5f20d..873ca582 100755 --- a/bin/cnvpytor_to_vcf.py +++ b/bin/cnvpytor_to_vcf.py @@ -28,7 +28,8 @@ Only a base record with FILTER PASS (or .) and an allele other than 0 in its first sample covers a call: a filtered Manta record (MinGQ and the rest) or one genotyped 0/0 or ./. stays in the file but hides nothing. A sites-only -base covers with every DEL and DUP. +base covers with every DEL and DUP. A Delly PASS call genotyped 0/0 or ./. is +not added either, so it covers no CNVpytor call. Every record says where it came from in INFO SVCALLER (the base's records get --base-name, default manta). An added record keeps the caller's ID, END, @@ -159,15 +160,16 @@ def first_gt(fields): return values[i] if i < len(values) and values[i] else "./." -def base_call_holds(fields, n_samples): - """True when a base record can stand for a call another caller makes. +def call_holds(fields, n_samples): + """True when a record is a call this sample carries. A record Manta filtered (MinGQ, MaxDepth and the rest) or genotyped as reference (0/0) or no-call (./.) is not a variant this sample carries; letting it cover a Delly or CNVpytor call would drop that call, and the summary (collect_summary.py) skips a 0/0 record, so the deletion would be - lost. A sites-only base (no sample column) has no genotype and counts as - a call. + lost. The same holds for a Delly PASS record genotyped 0/0 or ./.: it is + not added, and so does not cover a CNVpytor call. A sites-only record (no + sample column) has no genotype and counts as a call. """ if fields[6] not in ("PASS", "."): return False @@ -235,7 +237,7 @@ def main(argv=None): for f in base: info = info_dict(f[7]) kind, sp = sv_type(f, info), span(f, info) - if kind in ("DEL", "DUP") and sp and base_call_holds(f, n_samples): + if kind in ("DEL", "DUP") and sp and call_holds(f, n_samples): known.add(f[0], kind, sp[0], sp[1]) f[7] = ("" if f[7] in ("", ".") else f[7] + ";") + "SVCALLER=" + args.base_name @@ -270,7 +272,7 @@ def take(caller, records, keep): no_depth = [0] def delly_keep(f, info): - if f[6] != "PASS": + if f[6] != "PASS" or not call_holds(f, 1 if len(f) > 9 else 0): return False sp = span(f, info) fmt = f[8].split(":") if len(f) > 9 else [] diff --git a/docs/05-annotsv.md b/docs/05-annotsv.md index 3041eb38..d2bd8aa9 100644 --- a/docs/05-annotsv.md +++ b/docs/05-annotsv.md @@ -45,7 +45,7 @@ Step 05 and the Nextflow `ANNOTSV` stop with a pointer to `--annotsv-data` when AnnotSV classifies only the records of the one VCF it is given. The step builds that VCF, `annotsv/_annotsv_input.vcf`, in two passes inside the AnnotSV image: - **The net length change** (`bin/sv_svlen_fix.py`). Manta writes a deletion that also inserts a few bases with its REF and ALT sequences, `END` at the last deleted base and `SVLEN` as the deleted length. AnnotSV measures the overlapped coding sequence over `POS-1..END`, so REF 137 bp with ALT 2 bp read as 137 deleted coding bases: `Frameshift` "yes" and criterion 2E-1 (+0.90, class 4), although the change is -135 bp and keeps the frame. For every record whose REF and ALT are both sequence, the script sets `SVLEN` to `len(ALT) - len(REF)` and, for a net deletion, `END` to `POS - 1 + len(REF) - len(ALT)`: one base less than AnnotSV's own default `END` (`POS + len(REF) - len(ALT)`), so that its `POS-1..END` interval, which counts the anchor base, is as long as the net change. The old values stay in `ORIG_SVLEN` and `ORIG_END`. Symbolic records (``) and breakends pass unchanged. -- **The other callers' calls** (`bin/cnvpytor_to_vcf.py`). A large deletion that CNVpytor (step 18) or Delly (step 19) calls, but Manta calls as something else or misses, got no class before. The script adds Delly's PASS deletions and duplications that read depth does not contradict (below), then CNVpytor's with e-val (`pytorP1`) below 0.01 and a fraction of mapping-quality-0 reads (`pytorQ0`) below 0.5, each only when no record already in the file has the same type on the same chromosome with a reciprocal overlap of at least half the longer call. A deletion both of them call is added once, as Delly's. CNVpytor's copy-neutral LOH calls are left out. Each file is used when it exists, and the step says when one is missing. +- **The other callers' calls** (`bin/cnvpytor_to_vcf.py`). A large deletion that CNVpytor (step 18) or Delly (step 19) calls, but Manta calls as something else or misses, got no class before. The script adds Delly's PASS deletions and duplications with a non-reference genotype (not 0/0 or ./.) that read depth does not contradict (below), then CNVpytor's with e-val (`pytorP1`) below 0.01 and a fraction of mapping-quality-0 reads (`pytorQ0`) below 0.5, each only when no record already in the file has the same type on the same chromosome with a reciprocal overlap of at least half the longer call. A deletion both of them call is added once, as Delly's. CNVpytor's copy-neutral LOH calls are left out. Each file is used when it exists, and the step says when one is missing. - **Read depth for long Delly calls.** Short-read Delly can mark PASS a long deletion or duplication, even one of many megabases, whose read depth shows no change, and AnnotSV then classes it as the loss or gain of every gene it spans. So a Delly call longer than 1,000 bp joins the input only when Delly's own read-depth copy number (`FORMAT/RDCN`) agrees: below a normal region's copies for a deletion, above them for a duplication. The normal is 1 on a male's chrX or chrY outside the pseudoautosomal regions (`assets/par_grch38.bed`) and 2 everywhere else. With no sex given, a stretch of chrX or chrY outside them counts 1 and 2 as normal, so there a deletion needs RDCN 0 or 1 and a duplication 2 or more: give the sex (the second argument, or the samplesheet's `sex` column) for the strict check. A long call with no estimate (`.` or a negative RDCN, no reads in the flanks) is set aside too. A call of 1,000 bp or less rests on its split reads and read pairs alone, and a Delly VCF without the RDCN field is read as before. The couple screen of [step 39](39-couple.md) uses the same rule, from the same code. The step's log line `Callers beside manta: ...`, repeated in the input VCF's `##svcallers_added` header line, counts the Delly calls set aside. Every record names its caller in `INFO/SVCALLER` (`manta`, `delly` or `cnvpytor`), which AnnotSV copies into its `INFO` column. With `SV_VCF` set, that file alone is annotated, after the length pass, with `SVCALLER=sv_vcf`. diff --git a/tests/test_annotsv_sv_inputs.py b/tests/test_annotsv_sv_inputs.py index ba4ee45c..a153490d 100755 --- a/tests/test_annotsv_sv_inputs.py +++ b/tests/test_annotsv_sv_inputs.py @@ -327,6 +327,18 @@ def test_extra_callers(w): check("a sites-only base deletion covers Delly's DEL00000002 (no genotype to read)", rc == 0 and "DEL00000002" not in [r[2] for r in records(out)], (records(out), err)) + # A Delly PASS deletion genotyped 0/0 or ./. is no call either: it is not + # added, and CNVpytor's deletion at the same place is added instead. + for gt in ("0/0", "./."): + with gzip.open(delly, "wt") as f: + f.write(DELLY.replace("\t1/1:40\n", "\t{0}:40\n".format(gt))) + rc, out, err = run("cnvpytor_to_vcf.py", "--base", base, "--delly", delly, "--cnvpytor", cnvp) + ids = [r[2] for r in records(out)] + check("a {0} Delly deletion is not added and does not hide CNVpytor_del0".format(gt), + rc == 0 and "DEL00000002" not in ids and "CNVpytor_del0" in ids, (ids, err)) + with gzip.open(delly, "wt") as f: + f.write(DELLY) + rc, out, err = run("cnvpytor_to_vcf.py", "--base", "-", "--base-name", "sv_vcf", stdin=BASE) check("no other caller: stdin in, the base out with SVCALLER=sv_vcf", rc == 0 and "none given" in err and [info(r).get("SVCALLER") for r in records(out)] == ["sv_vcf"] * 2, err)