From 5ba2f3698919798263815171715c954adf060858 Mon Sep 17 00:00:00 2001
From: GeiserX <9169332+GeiserX@users.noreply.github.com>
Date: Fri, 6 Nov 2026 22:22:59 +0100
Subject: [PATCH 1/7] fix(cn): show step 35's CYP21A2, gene-dosage and interval
 tables in both reports

The step 35 card said CYP21A2 was not assessed while step 35 calls it, and
the dosage and interval tables reached no report. The summary now carries
the dosage table, both reports print the CYP21A2 call, one line per dosage
locus and the interval counts with each changed row, and a missing table
says it is not in the run's output. A header without the copy-number
columns is unreadable, also with no row under it.

Bead: pgp-5bja.21
Bead: pgp-7fmg
---
 bin/collect_summary.py        |  32 ++++++++--
 bin/render_report.py          | 100 +++++++++++++++++++++++++++++--
 tests/test_collect_summary.py | 108 +++++++++++++++++++++++++++++++++-
 3 files changed, 228 insertions(+), 12 deletions(-)

diff --git a/bin/collect_summary.py b/bin/collect_summary.py
index aae233df..9534329b 100755
--- a/bin/collect_summary.py
+++ b/bin/collect_summary.py
@@ -783,6 +783,13 @@ def smn_reading(smn1):
 PARASCOPY_COLUMNS = {"chrom", "start", "end", "agCN", "agCN_filter", "agCN_qual", "psCN", "psCN_qual", "psCN_filter"}
 
 
+def tsv_columns(path):
+    """The header of a tab-separated table as a set ({} for an empty file), so
+    a table with a header and no row is checked for its columns too."""
+    with open_text(path) as f:
+        return set(next(csv.reader(f, delimiter="\t"), []))
+
+
 def parascopy_tables(d, s):
     """{table: path} of step 35's catalog tables (<s>_<table>_copy_number.tsv)."""
     out = {}
@@ -800,7 +807,7 @@ def parascopy_locus(name, path, row):
     out = {"table": name, "genes": (row or {}).get("genes") or name, "tag": (row or {}).get("tag") or ".",
            "reading": (row or {}).get("reading") or "."}
     rows = read_tsv(path)
-    if rows and not PARASCOPY_COLUMNS <= set(rows[0]):
+    if not PARASCOPY_COLUMNS <= tsv_columns(path):
         return dict(out, status="unreadable", cn=None,
                     text=f"unreadable: {os.path.basename(path)} lacks the copy-number columns")
     na = next((r["agCN_filter"] for r in rows if (r.get("agCN_filter") or "").startswith("not assessed")), None)
@@ -846,13 +853,14 @@ def sec_parascopy(d, s):
     hp = os.path.join(d, "paralogs", f"{s}_hba_copy_number.tsv")
     ip = os.path.join(d, "paralogs", f"{s}_cn_intervals.tsv")
     cp = os.path.join(d, "paralogs", f"{s}_cyp21a2_psv.tsv")
-    p = (tables.get("smn") or next(iter(tables.values()), None) or (hp if os.path.isfile(hp) else None)
-         or (ip if os.path.isfile(ip) else None) or (cp if os.path.isfile(cp) else None))
+    dp = os.path.join(d, "paralogs", f"{s}_cn_dosage.tsv")
+    p = (tables.get("smn") or next(iter(tables.values()), None)
+         or next((x for x in (hp, ip, cp, dp) if os.path.isfile(x)), None))
     if not p:
         return None, {}
     if "smn" in tables:
         rows = read_tsv(tables["smn"])
-        if rows and not PARASCOPY_COLUMNS <= set(rows[0]):
+        if not PARASCOPY_COLUMNS <= tsv_columns(tables["smn"]):
             raise Unreadable(f"{os.path.basename(tables['smn'])} lacks the chrom, start, end, agCN, agCN_filter, "
                              "agCN_qual, psCN, psCN_qual or psCN_filter column")
         smn1 = parascopy_cn_at(rows, *paraphase_report.SMN1_EXON7)
@@ -882,6 +890,9 @@ def sec_parascopy(d, s):
                # CYP21A2 from depth and PSV sites (bin/copy_number.py, the
                # CYP21A2 row): copies, module total, finding; None without it
                "cyp21a2": read_hba(cp) if os.path.isfile(cp) else None,
+               # the gene-dosage rows (HP, LCE3BC, CCL3L1): total copies from
+               # depth, one entry per locus; None without step 35's table
+               "dosage": read_cn_dosage(dp) if os.path.isfile(dp) else None,
                "reading": reading,
                # every other catalog table, in the catalog's order (a table the
                # catalog does not name, from a PARALOG_LOCI of the user's, last)
@@ -894,7 +905,7 @@ def read_cn_intervals(path):
     """Step 35's interval table (bin/copy_number.py intervals): the rows with
     a deletion or a duplication, and a count of the rows by status."""
     rows = read_tsv(path)
-    if rows and not {"locus", "name", "status", "copies", "expected", "change"} <= set(rows[0]):
+    if not {"locus", "name", "chrom", "start", "end", "status", "copies", "expected", "change"} <= tsv_columns(path):
         raise Unreadable(f"{os.path.basename(path)} lacks the locus, name, status, copies, expected or change column")
     counts = {}
     for r in rows:
@@ -904,6 +915,17 @@ def read_cn_intervals(path):
                         for r in rows if r["change"] in ("deletion", "duplication")]}
 
 
+DOSAGE_COLUMNS = ("locus", "status", "copies", "estimate", "reason", "reading")
+
+
+def read_cn_dosage(path):
+    """Step 35's gene-dosage table (HP, LCE3BC, CCL3L1): one entry per locus,
+    in the table's order."""
+    if not set(DOSAGE_COLUMNS) <= tsv_columns(path):
+        raise Unreadable(f"{os.path.basename(path)} lacks the {', '.join(DOSAGE_COLUMNS[:-1])} or reading column")
+    return [{k: r.get(k) or "." for k in DOSAGE_COLUMNS} for r in read_tsv(path)]
+
+
 def read_hba(path):
     """Step 35's HBA table (metric, value) as a dict, and its CYP21A2 table
     the same way; a table without a status row reads as not assessed, with
diff --git a/bin/render_report.py b/bin/render_report.py
index 1688f9c2..8cc7b84c 100755
--- a/bin/render_report.py
+++ b/bin/render_report.py
@@ -1110,6 +1110,8 @@ def head(key, title=None):
             w(f"  {line}")
         for line in parascopy_locus_lines(S["parascopy"]["data"]):
             w(f"  {line}")
+        for line in cn_interval_lines(S["parascopy"]["data"].get("cn_intervals")):
+            w(f"  {line}")
         w(f"  {PARASCOPY_NOTE}")
         w("")
     at("paralogs")
@@ -1670,10 +1672,12 @@ def paralog_lines(d):
     return out
 
 
-PARASCOPY_NOTE = ("Parascopy reads short reads at SMN1 and SMN2 and at the other loci of its catalog; CYP21A2 has "
-                  "no Parascopy model and is not assessed. A quality under 20 or a filter other than PASS leaves a copy unresolved. HBA1 and "
-                  "HBA2 have no Parascopy model either: step 35 counts the alpha-globin genes together from read "
-                  "depth against flanking sequence, which cannot say which gene or which chromosome lost a copy. "
+PARASCOPY_NOTE = ("Parascopy reads short reads at SMN1 and SMN2 and at the other loci of its catalog. A quality "
+                  "under 20 or a filter other than PASS leaves a copy unresolved. HBA1 and HBA2, CYP21A2 and the "
+                  "gene-dosage loci have no Parascopy model: step 35 counts them from read depth against flanking "
+                  "sequence (CYP21A2 also from six sites where it differs from its pseudogene CYP21A1P), which "
+                  "cannot say which chromosome lost or gained a copy. The copy-number intervals (ClinGen's "
+                  "recurrent CNV regions, carrier deletion alleles, the AZF regions) are read from depth too. "
                   "A research estimate, not a diagnostic test: confirm a result with a clinical laboratory (MLPA "
                   "or qPCR; gap-PCR for alpha-thalassaemia).")
 
@@ -1685,7 +1689,92 @@ def parascopy_lines(d):
             f"SMN2 copies: {smn2.get('text', '.')}",
             f"SMN1 and SMN2 together (agCN) at exon 7: {d.get('agcn', '.')} (quality {d.get('agcn_qual', '.')}, "
             f"filter {d.get('agcn_filter', '.')})",
-            f"Reading: {d.get('reading', '.')}"] + ([hba_line(d["hba"])] if d.get("hba") else [])
+            f"Reading: {d.get('reading', '.')}"] + ([hba_line(d["hba"])] if d.get("hba") else []) \
+        + [cyp21a2_line(d.get("cyp21a2"))] + dosage_lines(d.get("dosage"))
+
+
+def cyp21a2_line(c):
+    """Step 35's CYP21A2 call from depth and PSV sites, one line; not
+    assessed only when the table says so or is not there."""
+    label = "CYP21A2 (RCCX module, from depth and six PSV sites)"
+    if not c:
+        return f"{label}: not assessed (step 35's <id>_cyp21a2_psv.tsv is not in this run's output)"
+    if c.get("status") != "called":
+        why = c.get("message") or c.get("reason") or "."
+        return f"{label}: not assessed ({why.removeprefix('not assessed: ')})"
+    value = (f"{c.get('cyp21a2_copies', '.')} CYP21A2 copies, module total {c.get('module_total', '.')} of "
+             f"{c.get('expected_total', '.')}, finding {c.get('finding', '.')}")
+    if c.get("conversion_sites", ".") != ".":
+        value += f" at {c['conversion_sites']}"
+    if c.get("reading", ".") != ".":
+        value += f". Reading: {c['reading']}"
+    return f"{label}: {value}"
+
+
+DOSAGE_NAMES = {"HP": "HP (haptoglobin)", "LCE3BC": "LCE3B/LCE3C", "CCL3L1": "CCL3L1/CCL3L3"}
+
+
+def dosage_lines(rows):
+    """Step 35's gene-dosage loci (HP, LCE3BC, CCL3L1), one line each."""
+    if rows is None:
+        return ["Gene-dosage loci (HP, LCE3B/LCE3C, CCL3L1): not assessed (step 35's <id>_cn_dosage.tsv is not "
+                "in this run's output)"]
+    if not rows:
+        return ["Gene-dosage loci (HP, LCE3B/LCE3C, CCL3L1): not assessed (step 35's dosage table has no row)"]
+    out = []
+    for r in rows:
+        label = f"{DOSAGE_NAMES.get(r['locus'], r['locus'])} copies (from depth)"
+        if r["status"] == "called":
+            value = f"{r['copies']} (estimate {r['estimate']})"
+            if r["reading"] != ".":
+                value += f". Reading: {r['reading']}"
+        else:
+            value = f"{r['status']} ({r['reason']})"
+        out.append(f"{label}: {value}")
+    return out
+
+
+CN_INTERVALS_LABEL = "Copy-number intervals (ClinGen recurrent regions, carrier deletions, AZF; from depth)"
+
+
+def cn_interval_summary(iv):
+    """The interval table in one line: rows by status and how many changed."""
+    if iv is None:
+        return f"{CN_INTERVALS_LABEL}: not assessed (step 35's <id>_cn_intervals.tsv is not in this run's output)"
+    counts = ", ".join(f"{n} {k}" for k, n in sorted((iv.get("by_status") or {}).items())) or "no row"
+    n = len(iv.get("changes") or [])
+    return (f"{CN_INTERVALS_LABEL}: {iv.get('rows', 0)} rows ({counts}); "
+            f"{n} with a deletion or a duplication" if n else
+            f"{CN_INTERVALS_LABEL}: {iv.get('rows', 0)} rows ({counts}); no deletion or duplication")
+
+
+def cn_interval_place(r):
+    """chrom:start-end (1-based) of an interval row of the table (BED)."""
+    b = r.get("start", ".")
+    return f"{r.get('chrom', '.')}:{int(b) + 1 if str(b).isdigit() else b}-{r.get('end', '.')}"
+
+
+def cn_interval_lines(iv):
+    """The interval table for the text report: the summary, then one line per
+    row with a deletion or a duplication."""
+    return [cn_interval_summary(iv)] + [
+        f"  {r['change'].capitalize()}: {r['name']} ({r['locus']}, {cn_interval_place(r)}), "
+        f"{r['copies']} of {r['expected']} copies" for r in (iv or {}).get("changes") or []]
+
+
+def cn_intervals_html(iv):
+    """The interval table for the HTML card: the summary, then a table of the
+    rows with a deletion or a duplication."""
+    label, value = cn_interval_summary(iv).split(": ", 1)
+    out = [stat(label, E(value))]
+    changes = (iv or {}).get("changes") or []
+    if changes:
+        rows = "".join(f"<tr><td>{E(r['name'])}</td><td>{E(r['locus'])}</td><td>{E(cn_interval_place(r))}</td>"
+                       f"<td>{E(r['copies'])} of {E(r['expected'])}</td><td>{E(r['change'])}</td></tr>"
+                       for r in changes)
+        out.append("    <table><tr><th>Region</th><th>Row</th><th>Place</th><th>Copies</th><th>Change</th></tr>"
+                   + rows + "</table>")
+    return out
 
 
 PARASCOPY_LOCI_NOTE = ("The other loci of step 35's catalog (assets/paralog_loci.tsv): agCN is the copies of the "
@@ -2685,6 +2774,7 @@ def a(key, lines):
                [stat("Status", done_badge(S["parascopy"]))]
                + [stat(*(E(x) for x in line.split(": ", 1))) for line in parascopy_lines(S["parascopy"]["data"])]
                + parascopy_loci_html(S["parascopy"]["data"])
+               + cn_intervals_html(S["parascopy"]["data"].get("cn_intervals"))
                + [small_note(PARASCOPY_NOTE)]))
     unchecked = paralogs_unchecked(S)
     if unchecked:
diff --git a/tests/test_collect_summary.py b/tests/test_collect_summary.py
index 94d90f13..a5c90df5 100644
--- a/tests/test_collect_summary.py
+++ b/tests/test_collect_summary.py
@@ -1409,6 +1409,109 @@ def made_up(name, tg, t, ugt, r117h=None, ugt_pc=None, stranger=False):
             "homologous_regions\n")
 
 
+IV_COLS = "locus\tname\tchrom\tstart\tend\tstatus\tcopies\texpected\tchange\testimate\treason\n"
+CYP_DELETION = ("metric\tvalue\nstatus\tcalled\ncyp21a2_copies\t1\nmodule_total\t3\nexpected_total\t4\n"
+                "estimate\t3.02\nfinding\tdeletion\nconversion_sites\t.\ntarget_depth\t45.30\n"
+                "control_depth\t30.10\nreason\t.\nreading\tone RCCX module fewer and one CYP21A2 copy fewer at the "
+                "differing sites: the pattern of the 30 kb CYP21A1P/CYP21A2 chimeric deletion\n")
+
+
+def cn_report_section(work, schema):
+    """Step 35's tables that no Parascopy model makes reach both reports: the
+    CYP21A2 call, the gene-dosage rows and the interval rows; and a table with
+    a header but without the copy-number columns reads as unreadable."""
+    d = os.path.join(work, "cn_report", "S")
+    put(f"{d}/paralogs/S_smn_copy_number.tsv", SMN_COLS
+        + "chr5\t70950920\t70953320\tSMN1\tPASS\t4\t60\tPASS\t2,2\t33,30\tchr5:70075496-70077900:+\n")
+
+    def card(html):
+        return html.split("<h2>SMN1/SMN2 Copy Number (Parascopy)</h2>", 1)[-1].split("<h2>", 1)[0]
+
+    # No CYP21A2, dosage or interval table: each says it is not in the output, never "has no model".
+    summ, txt, html = render(d)
+    check("cn report: without the tables each line says the table is not in this run's output",
+          "  CYP21A2 (RCCX module, from depth and six PSV sites): not assessed (step 35's <id>_cyp21a2_psv.tsv "
+          "is not in this run's output)\n" in txt
+          and "  Gene-dosage loci (HP, LCE3B/LCE3C, CCL3L1): not assessed (step 35's <id>_cn_dosage.tsv" in txt
+          and "  Copy-number intervals (ClinGen recurrent regions, carrier deletions, AZF; from depth): not assessed "
+              "(step 35's <id>_cn_intervals.tsv" in txt
+          and "is not assessed" not in txt and "CYP21A2 has no Parascopy model" not in txt, txt)
+
+    put(f"{d}/paralogs/S_cyp21a2_psv.tsv", CYP_DELETION)
+    put(f"{d}/paralogs/S_cn_dosage.tsv", "locus\tstatus\tcopies\testimate\treason\treading\n"
+        "HP\tcalled\t3\t3.04\t.\tthree copies of the duplicated HP exons: Hp2-1\n"
+        "LCE3BC\tcalled\t1\t0.98\t.\tLCE3B and LCE3C deleted on one chromosome (LCE3C_LCE3B-del on one)\n"
+        "CCL3L1\tnot assessed\t.\t.\tcheck_failed\tthe depth check failed; see the log of step 35\n")
+    put(f"{d}/paralogs/S_cn_intervals.tsv", IV_COLS
+        + "22q11.2_LCR_A-LCR_D\t22q11.2 recurrent (DGS/VCFS) region (proximal, A-D) (includes TBX1)\tchr22\t"
+          "18924717\t21111383\tcalled\t1\t2\tdeletion\t1.02\t.\n"
+        + "CFTR_dele2_3\tCFTRdele2,3\tchr7\t117479000\t117481000\tcalled\t2\t2\tno change\t1.97\t.\n"
+        + "AZFa\tAZFa\tchrY\t12904108\t13696000\tnot assessed\t.\t.\t.\t.\tno_y\n")
+    summ, txt, html = render(d)
+    pd = (summ["sections"].get("parascopy") or {}).get("data") or {}
+    cs = card(html)
+    cyp = ("CYP21A2 (RCCX module, from depth and six PSV sites): 1 CYP21A2 copies, module total 3 of 4, "
+           "finding deletion. Reading: one RCCX module fewer")
+    check("cn report: a CYP21A2 deletion is in the summary and in both reports",
+          (pd.get("cyp21a2") or {}).get("finding") == "deletion" and f"  {cyp}" in txt
+          and ">CYP21A2 (RCCX module, from depth and six PSV sites)<" in cs
+          and "1 CYP21A2 copies, module total 3 of 4, finding deletion" in cs, (pd.get("cyp21a2"), txt))
+    check("cn report: the dosage rows are in the summary, one entry per locus in the table's order",
+          [(x["locus"], x["status"], x["copies"]) for x in pd.get("dosage") or []]
+          == [("HP", "called", "3"), ("LCE3BC", "called", "1"), ("CCL3L1", "not assessed", ".")], pd.get("dosage"))
+    check("cn report: both reports give one line per dosage locus, a failed one with its reason",
+          "  HP (haptoglobin) copies (from depth): 3 (estimate 3.04). Reading: three copies of the duplicated HP "
+          "exons: Hp2-1\n" in txt
+          and "  LCE3B/LCE3C copies (from depth): 1 (estimate 0.98). Reading: LCE3B and LCE3C deleted" in txt
+          and "  CCL3L1/CCL3L3 copies (from depth): not assessed (check_failed)\n" in txt
+          and ">HP (haptoglobin) copies (from depth)<" in cs and ">CCL3L1/CCL3L3 copies (from depth)<" in cs, txt)
+    check("cn report: the interval table's counts and its one deletion are in the text report",
+          "  Copy-number intervals (ClinGen recurrent regions, carrier deletions, AZF; from depth): 3 rows "
+          "(2 called, 1 not assessed); 1 with a deletion or a duplication\n" in txt
+          and "    Deletion: 22q11.2 recurrent (DGS/VCFS) region (proximal, A-D) (includes TBX1) "
+              "(22q11.2_LCR_A-LCR_D, chr22:18924718-21111383), 1 of 2 copies\n" in txt
+          and "CFTRdele2,3" not in txt, txt)
+    check("cn report: the HTML card has the counts and a row for the deletion",
+          "3 rows (2 called, 1 not assessed); 1 with a deletion or a duplication" in cs
+          and "<td>22q11.2_LCR_A-LCR_D</td><td>chr22:18924718-21111383</td><td>1 of 2</td><td>deletion</td>" in cs,
+          cs)
+    check("cn report: the summary validates", not validate.validate(schema, json.loads(json.dumps(summ))),
+          validate.validate(schema, json.loads(json.dumps(summ))))
+    put(f"{d}/paralogs/S_cyp21a2_psv.tsv", "metric\tvalue\nstatus\tnot assessed\ncyp21a2_copies\t.\n"
+        "reason\tpsv_low_depth\n")
+    _, txt, _ = render(d)
+    check("cn report: a CYP21A2 table that says not assessed gives its reason",
+          "  CYP21A2 (RCCX module, from depth and six PSV sites): not assessed (psv_low_depth)\n" in txt, txt)
+    put(f"{d}/paralogs/S_cn_intervals.tsv", IV_COLS
+        + "CFTR_dele2_3\tCFTRdele2,3\tchr7\t117479000\t117481000\tcalled\t2\t2\tno change\t1.97\t.\n")
+    _, txt, _ = render(d)
+    check("cn report: an interval table without a change says so",
+          "from depth): 1 rows (1 called); no deletion or duplication\n" in txt, txt)
+
+    # A header without the copy-number columns is unreadable, also with no row under it.
+    put(f"{d}/paralogs/S_cn_dosage.tsv", "locus\tstatus\n")
+    ps = collect_summary.collect("S", d)["sections"].get("parascopy") or {}
+    check("cn report: a dosage table without the copies column is unreadable", ps.get("state") == "unreadable", ps)
+    os.remove(f"{d}/paralogs/S_cn_dosage.tsv")
+    put(f"{d}/paralogs/S_cn_intervals.tsv", "locus\tname\tstatus\n")
+    ps = collect_summary.collect("S", d)["sections"].get("parascopy") or {}
+    check("cn report: a header-only interval table without the copy columns is unreadable",
+          ps.get("state") == "unreadable", ps)
+    os.remove(f"{d}/paralogs/S_cn_intervals.tsv")
+    put(f"{d}/paralogs/S_gba1_copy_number.tsv", "chrom\tstart\tend\tlocus\n")
+    pd = (collect_summary.collect("S", d)["sections"].get("parascopy") or {}).get("data") or {}
+    check("cn report: a header-only catalog table without the copy-number columns is unreadable, not not assessed",
+          {x["table"]: x["status"] for x in pd.get("loci", [])}.get("gba1") == "unreadable", pd.get("loci"))
+    put(f"{d}/paralogs/S_gba1_copy_number.tsv", SMN_COLS)
+    pd = (collect_summary.collect("S", d)["sections"].get("parascopy") or {}).get("data") or {}
+    check("cn report: a header-only catalog table with every column stays not assessed (no region row)",
+          {x["table"]: x["status"] for x in pd.get("loci", [])}.get("gba1") == "not assessed", pd.get("loci"))
+    put(f"{d}/paralogs/S_smn_copy_number.tsv", "chrom\tstart\tend\tagCN\n")
+    ps = collect_summary.collect("S", d)["sections"].get("parascopy") or {}
+    check("cn report: a header-only SMN table without the psCN columns is unreadable",
+          ps.get("state") == "unreadable", ps)
+
+
 def parascopy_section(work, schema):
     """Step 35's own section, and the not-checked line without steps 35 and 35b."""
     d = os.path.join(work, "parascopy", "S")
@@ -1436,8 +1539,8 @@ def parascopy_section(work, schema):
     card_s = html.split("<h2>SMN1/SMN2 Copy Number (Parascopy)</h2>", 1)[-1].split("<h2>", 1)[0]
     check("parascopy: the HTML card has both genes, the reading and the short-read note",
           ">SMN1 copies<" in card_s and "2 (quality 33, filter PASS)" in card_s and ">SMN2 copies<" in card_s
-          and "not a carrier by copy number" in card_s and "CYP21A2 has no Parascopy model" in card_s
-          and "counts the alpha-globin genes together from read depth" in card_s
+          and "not a carrier by copy number" in card_s and "CYP21A2 and the gene-dosage loci have no Parascopy model" in card_s
+          and "step 35 counts them from read depth against flanking sequence" in card_s
           and "Not checked" not in html, card_s)
     check("parascopy: HBA1/HBA2 listed as not assessed without step 35's HBA table or step 35b",
           any(x.startswith("HBA1/HBA2 copy number") for x in summ["not_assessed"])
@@ -2762,6 +2865,7 @@ def crow(gene, result, risk, fa, fb):
         report_followups_sections(work, schema)
         str_extra_readings(work, schema)
         parascopy_section(work, schema)
+        cn_report_section(work, schema)
         report_sections_mito_roh_telomere_cnv(work, schema)
         roh_upd_karyotype(work, schema)
         hla_optitype_kir_section(work, schema)

From 3003028437263c152f9dd3b06a1d0e9df9860b55 Mon Sep 17 00:00:00 2001
From: GeiserX <9169332+GeiserX@users.noreply.github.com>
Date: Fri, 6 Nov 2026 22:23:13 +0100
Subject: [PATCH 2/7] fix(couple): pass the CYP21A2 table in step 39 and name
 read-depth alleles under the shared-risk heading

Bead: pgp-q4jt
Bead: pgp-5bja.21
---
 bin/couple_report.py        | 48 +++++++++++++++++++++++++++++++++----
 scripts/39-couple.sh        |  7 ++++--
 tests/test_couple_report.py | 37 ++++++++++++++++++++++++++++
 3 files changed, 86 insertions(+), 6 deletions(-)

diff --git a/bin/couple_report.py b/bin/couple_report.py
index 5d6fc7d9..8b830b32 100755
--- a/bin/couple_report.py
+++ b/bin/couple_report.py
@@ -451,6 +451,11 @@ def cyp21a2_finding(path, genes):
 CN_ALLELES = (("GJB6_D13S1830", "del(GJB6-D13S1830)"), ("CFTR_dele2_3", "CFTRdele2,3"),
               ("HBB_exon3_del", "HBB exon 3 deletion"), ("HBB_619bp_del", "HBB 619 bp deletion"),
               ("HBD_HBB_Sicilian_del", "Sicilian (delta-beta)0 deletion"))
+# The genes each allele of CN_ALLELES is a change of, as the carrier list
+# names them (del(GJB6-D13S1830) is a DFNB1 allele beside GJB2 changes).
+CN_ALLELE_GENES = (("GJB2/GJB6 (DFNB1)", ("GJB2", "GJB6"), ("GJB6_D13S1830",)),
+                   ("CFTR", ("CFTR",), ("CFTR_dele2_3",)),
+                   ("HBB", ("HBB",), ("HBB_exon3_del", "HBB_619bp_del", "HBD_HBB_Sicilian_del")))
 # Named as deletions, so present reads as "the deletion is there" (as for CN_ALLELES).
 CN_AZF = (("AZFa", "AZFa deletion"), ("AZFb", "AZFb deletion"),
           ("AZFc_DAZ", "AZFc deletion (DAZ copies; gr/gr and b2/b3 leave two of four)"))
@@ -514,6 +519,32 @@ def cnv_note(row, dels):
     return ""
 
 
+def cn_shared(a, b, intervals, found):
+    """The genes where both partners carry a change once the read-depth
+    alleles count: one has a present allele of CN_ALLELES and the other a
+    present allele or a finding in the same gene. The gene table leaves
+    these alleles out, so without this a couple who both carry CFTRdele2,3
+    would read as sharing nothing. intervals: {partner: {locus: row}} of the
+    partners with an interval table."""
+    names = dict(CN_ALLELES)
+    out = []
+    for label, genes, loci in CN_ALLELE_GENES:
+        have, by_depth = {}, False
+        for s in (a, b):
+            iv = intervals.get(s) or {}
+            alleles = ["%s present (read depth, step 35: %s)" % (names[x], cn_state(iv.get(x))[1])
+                       for x in loci if cn_state(iv.get(x))[0] == "present"]
+            hits = [f["text"] for f in found.get(s, []) if f["gene"] in genes]
+            by_depth = by_depth or bool(alleles)
+            if alleles or hits:
+                have[s] = "; ".join(alleles + hits)
+        if len(have) == 2 and by_depth:
+            out.append("**%s**: %s: %s; %s: %s. Both partners may carry a change in %s: the gene table does not "
+                       "count read-depth alleles, so a clinical test has to confirm them before the chance for a "
+                       "child can be given." % (label, a, have[a], b, have[b], label))
+    return out
+
+
 def cn_lines(s, sex, intervals, dels):
     """The report lines of a partner's carrier deletion alleles and AZF regions."""
     out = []
@@ -648,7 +679,7 @@ def write_table(path, rows):
            "can say what a negative screen leaves out for your ancestry.")
 
 
-def write_markdown(path, a, b, sexes, rows, assessed, gaps, ngenes, common, cn=None):
+def write_markdown(path, a, b, sexes, rows, assessed, gaps, ngenes, common, cn=None, shared_cn=None):
     def sx(s):
         return sexes.get(s) or "sex unknown"
     out = ["# Couple carrier screen: %s and %s" % (a, b), "",
@@ -671,7 +702,12 @@ def section(title, keep, none):
         out.append("")
 
     section("Shared risk: both partners carry a change in the same gene", "shared",
+            "No gene of the list where both partners carry a change in the gene table; read depth (step 35) "
+            "adds a possible shared risk, below." if shared_cn else
             "No gene of the list where both partners carry a change.")
+    if shared_cn:
+        out.extend(["Possible shared risk from read depth (step 35's carrier deletion alleles, not in the gene "
+                    "table):", ""] + ["- " + x for x in shared_cn] + [""])
     section("X-linked findings", "x_linked", "No change in an X-linked gene of the list.")
     section("Genes where only one partner carries a change", "one_partner", "None.")
     if cn:
@@ -751,10 +787,11 @@ def cmd_screen(a):
             print("%s and %s: not screened, %s %s no ClinVar screen (step 06's <sample>_clinvar_hits.vcf)"
                   % (x, y, " and ".join(lack), "has" if len(lack) == 1 else "have"))
             continue
-        found, assessed, gaps, cn = {}, {}, [], []
+        found, assessed, gaps, cn, ivs = {}, {}, [], [], {}
         for s in (x, y):
             if s in cni:
-                cn += cn_lines(s, sexes[s], read_cn_intervals(cni[s]), cnvpytor_deletions(cnv[s]) if s in cnv else None)
+                ivs[s] = read_cn_intervals(cni[s])
+                cn += cn_lines(s, sexes[s], ivs[s], cnvpytor_deletions(cnv[s]) if s in cnv else None)
             elif cni:
                 gaps.append("Carrier deletion alleles of %s: no interval table (step 35)." % s)
             aside = []
@@ -812,13 +849,16 @@ def cmd_screen(a):
         write_table(base + ".tsv", rows)
         for s in (x, y):
             write_table(os.path.join(a.out_dir, "%s_couple_carrier.tsv" % s), rows)
-        write_markdown(base + ".md", x, y, sexes, rows, assessed, gaps, len(genes), common, cn)
+        shared_cn = cn_shared(x, y, ivs, found)
+        write_markdown(base + ".md", x, y, sexes, rows, assessed, gaps, len(genes), common, cn, shared_cn)
         n = {k: sum(1 for r in rows if r["result"] == k) for k in ("shared", "x_linked", "one_partner")}
         print("%s and %s: %d shared gene%s, %d X-linked, %d with one partner only (%s.tsv)" % (
             x, y, n["shared"], "" if n["shared"] == 1 else "s", n["x_linked"], n["one_partner"], base))
         for r in rows:
             if r["result"] == "shared":
                 print("  shared: %s (%s): %s" % (r["gene"], r["condition"], r["offspring_risk"]))
+        for line in shared_cn:
+            print("  possibly shared, from read depth: %s" % line.split(":", 1)[0].strip("*"))
         done += 1
     if not done:
         die(["no couple of the pedigree has a ClinVar screen for both partners: run step 06 for each"], 1)
diff --git a/scripts/39-couple.sh b/scripts/39-couple.sh
index 8e2cfc50..60a010b3 100755
--- a/scripts/39-couple.sh
+++ b/scripts/39-couple.sh
@@ -12,7 +12,9 @@
 #        ExpansionHunter's (step 09, <sample>/expansion_hunter/<sample>_eh.vcf),
 #        Parascopy's SMN1 table (step 35, <sample>/paralogs/<sample>_smn_copy_number.tsv),
 #        step 35's interval table (<sample>/paralogs/<sample>_cn_intervals.tsv:
-#        the carrier deletion alleles and a male partner's AZF regions) and
+#        the carrier deletion alleles and a male partner's AZF regions), step
+#        35's CYP21A2 table (<sample>/paralogs/<sample>_cyp21a2_psv.tsv: the 30 kb
+#        chimeric deletion and pseudogene conversions) and
 #        CNVpytor's calls, its cross-check (step 18, <sample>/cnvpytor/<sample>_cnvs.vcf.gz),
 #        GENCODE's exons (setup.sh: reference/gencode.v<release>.basic.exons.bed)
 #        and AnnotSV's benign loss set (setup.sh --annotsv-data, with gnomAD-SV's
@@ -119,7 +121,8 @@ for s in ${NAMED[@]+"${NAMED[@]}"}; do
   [ -d "${GENOME_DIR}/${s}" ] || continue
   for kind in "clinvar:clinvar/${s}_clinvar_hits.vcf:06" "delly:delly/${s}_sv.vcf.gz:19" \
               "eh:expansion_hunter/${s}_eh.vcf:09" "smn:paralogs/${s}_smn_copy_number.tsv:35" \
-              "cn-intervals:paralogs/${s}_cn_intervals.tsv:35" "cnvpytor:cnvpytor/${s}_cnvs.vcf.gz:18"; do
+              "cn-intervals:paralogs/${s}_cn_intervals.tsv:35" "cyp21a2:paralogs/${s}_cyp21a2_psv.tsv:35" \
+              "cnvpytor:cnvpytor/${s}_cnvs.vcf.gz:18"; do
     IFS=: read -r k rel step <<< "$kind"
     if [ -f "${GENOME_DIR}/${s}/${rel}" ]; then
       ARGS+=("--${k}" "${s}=/genome/${s}/${rel}")
diff --git a/tests/test_couple_report.py b/tests/test_couple_report.py
index cb380e09..e6f0ca90 100644
--- a/tests/test_couple_report.py
+++ b/tests/test_couple_report.py
@@ -446,6 +446,43 @@ def ivrow(locus, chrom, s, e, status="called", copies="2", expected="2", change=
         check("cn alleles: without --cn-intervals no allele section is written",
               rc == 0 and "Carrier deletion alleles" not in open(f"{t}/nocn/couple_A_B.md").read())
 
+        # 11b. a read-depth allele in both partners, or one beside the other's
+        # finding in the same gene, is a possible shared risk the headline names
+        c_hits = put(f"{t}/C.vcf", HEAD + hit("UGT1A1"))
+        d_hits = put(f"{t}/D.vcf", HEAD + hit("HBB"))
+        c_iv = put(f"{t}/C_cn_intervals.tsv", ivh + ivrow("CFTR_dele2_3", "chr7", 117498312, 117519392, copies="1",
+                                                          change="deletion")
+                   + ivrow("HBB_exon3_del", "chr11", 5225254, 5225875, copies="1", change="deletion"))
+        d_iv = put(f"{t}/D_cn_intervals.tsv", ivh + ivrow("CFTR_dele2_3", "chr7", 117498312, 117519392, copies="1",
+                                                          change="deletion"))
+        ped2 = put(f"{t}/ped2.tsv", "sample\tsex\tfather\tmother\tpartner\nC\tfemale\t.\t.\tD\nD\tmale\t.\t.\tC\n")
+
+        def screen2(civ, div, out):
+            rc, log = run("screen", "--pedigree", ped2, "--genes", GENES, "--clinvar", f"C={c_hits}", "--clinvar",
+                          f"D={d_hits}", "--cn-intervals", f"C={civ}", "--cn-intervals", f"D={div}", "--out-dir", out)
+            p = f"{out}/couple_C_D.md"
+            return rc, log, open(p).read() if os.path.exists(p) else "", rows(f"{out}/couple_C_D.tsv") if os.path.exists(
+                f"{out}/couple_C_D.tsv") else []
+        rc, log, md, tab = screen2(c_iv, d_iv, f"{t}/cn2")
+        head = md.split("## Shared risk", 1)[-1].split("## X-linked", 1)[0]
+        check("cn shared: no gene-table row is shared, and the headline says read depth adds a possible shared risk",
+              rc == 0 and not any(r["result"] == "shared" for r in tab)
+              and "No gene of the list where both partners carry a change in the gene table; read depth (step 35) "
+                  "adds a possible shared risk, below." in head, (rc, log, head))
+        check("cn shared: CFTRdele2,3 present in both partners is named under the headline",
+              "- **CFTR**: C: CFTRdele2,3 present (read depth, step 35: 1 of 2 copies); D: CFTRdele2,3 present "
+              "(read depth, step 35: 1 of 2 copies). Both partners may carry a change in CFTR" in head, head)
+        check("cn shared: an HBB deletion of one partner beside the other's ClinVar HBB hit is named too",
+              "- **HBB**: C: HBB exon 3 deletion present (read depth, step 35: 1 of 2 copies); D: " in head
+              and "possibly shared, from read depth: CFTR" in log and "possibly shared, from read depth: HBB" in log,
+              (head, log))
+        none_iv = put(f"{t}/N_cn_intervals.tsv", ivh + ivrow("CFTR_dele2_3", "chr7", 117498312, 117519392))
+        rc, log, md, _ = screen2(none_iv, d_iv, f"{t}/cn3")
+        head = md.split("## Shared risk", 1)[-1].split("## X-linked", 1)[0]
+        check("cn shared, red: with C's CFTRdele2,3 absent and no HBB deletion the plain negative line stays",
+              rc == 0 and "No gene of the list where both partners carry a change.\n" in head
+              and "Possible shared risk" not in head and "possibly shared" not in log, head)
+
         # 10. genes
         gtf = put(f"{t}/g.gtf", "##description: evidence-based annotation of the human genome (GRCh38), version 99\n" + "".join(
             f'{c}\tHAVANA\tgene\t{s}\t{e}\t.\t+\t.\tgene_id "X"; gene_name "{g}";\n'

From 954395023a2a182f3e0478521746ea1e0a76b996 Mon Sep 17 00:00:00 2001
From: GeiserX <9169332+GeiserX@users.noreply.github.com>
Date: Fri, 6 Nov 2026 22:23:26 +0100
Subject: [PATCH 3/7] test(cn): e2e case for the step 35 tables in both reports
 and step 39, and docs to match

Bead: pgp-5bja.21
---
 docs/35-paralogs.md             |  14 +++-
 docs/39-couple.md               |   3 +-
 docs/cn-loci.md                 |   4 +-
 scripts/ci/e2e-case-paths.tsv   |   1 +
 scripts/ci/e2e-case-times.tsv   |   1 +
 tests/e2e/cn-report-1-tables.sh | 121 ++++++++++++++++++++++++++++++++
 6 files changed, 139 insertions(+), 5 deletions(-)
 create mode 100644 tests/e2e/cn-report-1-tables.sh

diff --git a/docs/35-paralogs.md b/docs/35-paralogs.md
index eca43180..6f1ad846 100644
--- a/docs/35-paralogs.md
+++ b/docs/35-paralogs.md
@@ -76,7 +76,7 @@ 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: <reason>`) and `.` in every other column |
 | `${SAMPLE}_<table>_copy_number.tsv` | The same, for every other catalog row (`gba1`, `pms2_pms2cl`, `c4` and the rest of the table above) |
 | `${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)). The couple screen (step 39) reads it for the carrier deletion alleles and AZF |
+| `${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_<site>` row per site (gene, pseudogene and other reads, copies, ratio to the other sites), `target_depth`, `control_depth`, `reason`, `reading` |
 | `${SAMPLE}_parascopy/<table>/` | 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)) |
 
@@ -225,12 +225,22 @@ The ratio of a site is its estimate over the median of the sites, so one copy co
 - `tests/test_copy_number.py` section 10 runs the CYP21A2 row on made-up bedcov and mpileup tables (`tests/fixtures/cn/CYP21A2_*`): unmodified, 2 copies and no finding; a planted chimeric deletion, 1 copy with the module total 3; a planted In2G conversion on one copy, a ratio of 0.5 at In2G and a possible conversion there. The same depths with every site even give no finding. The CYP21A2 block of step 35, cut out, writes those tables and a not-assessed one when mpileup fails, and the couple screen reads them.
 - The e2e case `cyp21a2-short-reads-1-hg002` runs the CYP21A2 block of step 35 on HG002's Illumina reads over the row, unmodified and with a planted In2G conversion and a planted chimeric deletion made from the same reads, and the block of the PARASCOPY module on the unmodified BAM, which must write the same table.
 - The e2e case `hba-short-read-accuracy-1-nf-report` runs the Nextflow `HTML_REPORT` with and without an HBA table and checks the HBA line and the not-assessed list, then stub-runs every tool and checks that `main.nf` passes the table to each `HTML_REPORT` task that got the SMN table, against a copy of the pipeline with that join removed.
+- `tests/test_collect_summary.py` (`cn_report_section`) reads made-up CYP21A2, dosage and interval tables into the summary and both reports: a CYP21A2 deletion, HP 3, LCE3BC 1, CCL3L1 not assessed, a 22q11.2 deletion with its counts; without the tables each line says the table is not in the output; a header without the copy-number columns is unreadable for the dosage, interval, catalog and SMN tables, also with no row.
+- The e2e case `cn-report-1-tables` runs steps 24, `generate-report.sh` and 39 in their images on made-up tables: both reports show the CYP21A2 deletion, the dosage loci and the 22q11.2 deletion; step 39 passes each partner's CYP21A2 table, and CFTRdele2,3 present in both partners is named under the shared-risk heading. Red controls: the CYP21A2 table removed, for the report and for one partner.
 
 ## In the Reports
 
 Both reports (step 24 and the Nextflow `HTML_REPORT`) read `${SAMPLE}_smn_copy_number.tsv` into a section of its own, **SMN1/SMN2 Copy Number (Parascopy, short reads)**. It shows the copy number of SMN1 and of SMN2 at exon 7 (the row of the table, or its homologous region, that holds chr5:70,951,946 for SMN1 and chr5:70,076,526 for SMN2), with quality and filter; the two together (agCN) there; and a one-line reading: no SMN1 copy is the pattern of SMA, one copy that of a carrier, two copies is not a carrier by copy number with a 2+0 carrier not ruled out. A copy with a quality under 20 or a filter other than `PASS` is shown as not resolved and gets no reading, a table with no row over SMN1 exon 7 says that instead of a reading, and a table that says not assessed gives that reason instead. With step 35b's table too, the Paraphase section also sets the two tools' SMN1 side by side. When `${SAMPLE}_hba_copy_number.tsv` is there, the section adds one line, `HBA1 and HBA2 together (alpha-globin genes, from depth)`, with the count, its reading and the estimate, or the reason it was not assessed. Then one line per other catalog table, in the catalog's order (a table of a row the catalog does not name, from a `PARALOG_LOCI` of your own, last): the genes and their copies together (agCN). One number when every region of the locus profile has filter `PASS` and quality 20 or more and they agree; else each region with its value, quality and filter. The catalog's reading follows a call. A table that says not assessed gives the reason. In the HTML card these are a table: genes, copies together, reading. The HBA line and the other loci do not depend on the SMN table: a catalog without the SMN1 row still shows them, and SMN1 then says it has no table. The Nextflow `PARASCOPY` task writes the same tables to the same folder, and `HTML_REPORT` reads them the same way (`main.nf` passes it every catalog table, and both catalogs).
 
-Neither report shows `${SAMPLE}_cyp21a2_psv.tsv` yet: the run summary carries it (the Parascopy section's `cyp21a2`), and the couple screen reads it with `couple_report.py screen --cyp21a2 ID=<id>_cyp21a2_psv.tsv`, where a deletion or a possible conversion is a CYP21A2 carrier finding. Step 39 and the Nextflow couple screen do not pass it yet.
+The same section shows the tables no Parascopy model makes, each read on its own:
+
+- `${SAMPLE}_cyp21a2_psv.tsv` as one line, `CYP21A2 (RCCX module, from depth and six PSV sites)`: the CYP21A2 copies, the module total of 4, the finding (with the sites of a possible conversion) and its reading, or the reason it was not assessed. The run summary carries it as the Parascopy section's `cyp21a2`.
+- `${SAMPLE}_cn_dosage.tsv` as one line per locus (HP, LCE3B/LCE3C, CCL3L1/CCL3L3): the copies, the estimate and the reading, or the status and its reason. The run summary carries it as `dosage`, one entry per locus.
+- `${SAMPLE}_cn_intervals.tsv` as one line of counts (rows by status, how many show a deletion or a duplication), then each row with a change: its region name, catalog row, place (1-based) and copies of expected. In the HTML card the changed rows are a table. The run summary carries it as `cn_intervals`.
+
+A table that is not in the run's output reads `not assessed (step 35's <id>_... is not in this run's output)`, so a missing table never reads as a normal result. A table with a header but without its columns makes the section unreadable; a catalog table without the copy-number columns makes that locus unreadable, also with no row under the header.
+
+The couple screen reads the CYP21A2 table with `couple_report.py screen --cyp21a2 ID=<id>_cyp21a2_psv.tsv`, where a deletion or a possible conversion is a CYP21A2 carrier finding. Step 39 passes it when a partner's folder holds it. The Nextflow `PARASCOPY` task writes the CYP21A2 and dosage tables, but `main.nf` does not pass them to `HTML_REPORT` or the couple screen yet, so a Nextflow run's reports say they are not in the output.
 
 When neither step 35 nor step 35b left a result, both reports print `Paralog copy number: not checked (neither step 35 nor step 35b ran)`, naming a step that failed instead, so a missing check never reads as a normal result. Unless step 35's depth count made a call, or step 35b called both HBA1 and HBA2 (a run limited to other regions with `--paraphase_regions` has neither, and a `no call` row does not count), the "Not Assessed" list also names HBA1/HBA2 copy number.
 
diff --git a/docs/39-couple.md b/docs/39-couple.md
index 76624c59..cd2d2b95 100644
--- a/docs/39-couple.md
+++ b/docs/39-couple.md
@@ -25,6 +25,7 @@ Before this step a couple had to compare two ClinVar tables by hand ([Multi-Samp
 | ClinVar pathogenic or likely pathogenic hit | step 06, `clinvar/<sample>_clinvar_hits.vcf` (required) | the genotype carries the alternate allele and ClinVar's `GENEINFO` names a gene of the carrier list |
 | Deletion | step 19 (Delly), `delly/<sample>_sv.vcf.gz` | `SVTYPE=DEL`, `FILTER` PASS and a per-sample `FT` that is not a failed filter, the genotype carries it, it is 50 bp or longer, it removes at least one exon of a listed gene (GENCODE exons), when it is longer than 1,000 bp read depth shows the loss, and it does not match a common deletion (see below) |
 | SMN1 copy number | step 35 (Parascopy), `paralogs/<sample>_smn_copy_number.tsv` | SMN1 exon 7 at 1 copy (a carrier) or 0, with filter PASS and quality 20 or more, by step 35's rule |
+| CYP21A2 deletion or conversion | step 35 (depth and six PSV sites), `paralogs/<sample>_cyp21a2_psv.tsv` | the table's `finding` is `deletion` or `conversion` (on both copies when a site has no CYP21A2 copy left); unconfirmed, short reads. Step 39 passes it; the Nextflow couple screen does not yet |
 | FMR1 and FXN repeats | step 09 (ExpansionHunter), `expansion_hunter/<sample>_eh.vcf` | an FMR1 allele of 55 CGG or more (a premutation; more than 200 is a full mutation), an FXN allele of 66 GAA or more |
 
 Short-read Delly can mark PASS a long deletion, even one of many megabases, that read depth contradicts, and such a call removes an exon of hundreds of listed genes at once. So a deletion longer than 1,000 bp counts only when Delly's own read-depth copy number (`FORMAT/RDCN`, the depth inside the deletion against its flanks) is below a normal region's copies: 1 on a male's chrX or chrY outside the pseudoautosomal regions, 2 everywhere else. When the pedigree gives no sex, every region is judged against 2. Delly infers the sex itself and scales RDCN to it, so a long artefact on an unsexed male's chrX or chrY reads 1 and still counts: give the sex in the pedigree and in the samplesheet. A heterozygous loss reads 1, a homozygous or hemizygous loss 0. A long deletion with RDCN at the normal count or more, or with no estimate (`.` or a negative value, no reads in the flanks), is set aside, and the table of what was assessed gives the count. 1,000 bp is Delly's own line between small and large SVs in its read counts (`indelsize` in Delly's source); a deletion that size or smaller rests on its split reads and read pairs alone. A deletion VCF without the RDCN field (not written by Delly 2.x) is read as before, without the depth check.
@@ -38,7 +39,7 @@ The table of what was assessed counts, for each partner, the deletions over a li
 
 Only the ClinVar screen is required. Each of the others is used when its file is there, and the report lists every one that was not.
 
-`couple_report.py screen` also takes step 35's interval table (`--cn-intervals ID=paralogs/<sample>_cn_intervals.tsv`, [the interval rows](cn-loci.md#the-interval-rows)). The report then lists, for each partner, del(GJB6-D13S1830), CFTRdele2,3, the HBB exon 3 deletion, the HBB 619 bp deletion and the Sicilian (delta-beta)0 deletion as present, absent or not assessable from read depth, and a male partner's AZFa, AZFb and AZFc deletions (Y chromosome microdeletions, a cause of male infertility) the same way, where present means the deletion is there; a female partner's AZF line says not applicable. The three beta-globin rows overlap (the Sicilian deletion covers both HBB rows, and the 619 bp and exon 3 rows share 487 bp), so one deletion can read as present on two or three of them; the report says so, and only a clinical test names the allele. With `--cnvpytor ID=cnvpytor/<sample>_cnvs.vcf.gz` (step 18), a CNVpytor call that disagrees with the depth call (a deletion over half the allele's span or more where depth sees none, or none where depth sees one) is printed beside it. A missing call is not printed for an allele shorter than 2 kb (the two HBB rows): step 18 counts in 1 kb bins and cannot see a deletion that short. These lines are a list, not carrier findings: they do not enter the gene table. Step 39 passes both when a partner's folder holds them (`paralogs/<sample>_cn_intervals.tsv` of step 35, `cnvpytor/<sample>_cnvs.vcf.gz` of step 18), and the Nextflow couple screen passes them when `--tools` has `parascopy` and `cnvpytor`.
+`couple_report.py screen` also takes step 35's interval table (`--cn-intervals ID=paralogs/<sample>_cn_intervals.tsv`, [the interval rows](cn-loci.md#the-interval-rows)). The report then lists, for each partner, del(GJB6-D13S1830), CFTRdele2,3, the HBB exon 3 deletion, the HBB 619 bp deletion and the Sicilian (delta-beta)0 deletion as present, absent or not assessable from read depth, and a male partner's AZFa, AZFb and AZFc deletions (Y chromosome microdeletions, a cause of male infertility) the same way, where present means the deletion is there; a female partner's AZF line says not applicable. The three beta-globin rows overlap (the Sicilian deletion covers both HBB rows, and the 619 bp and exon 3 rows share 487 bp), so one deletion can read as present on two or three of them; the report says so, and only a clinical test names the allele. With `--cnvpytor ID=cnvpytor/<sample>_cnvs.vcf.gz` (step 18), a CNVpytor call that disagrees with the depth call (a deletion over half the allele's span or more where depth sees none, or none where depth sees one) is printed beside it. A missing call is not printed for an allele shorter than 2 kb (the two HBB rows): step 18 counts in 1 kb bins and cannot see a deletion that short. These lines are a list, not carrier findings: they do not enter the gene table. When both partners carry something in one gene once these alleles count (CFTRdele2,3 present in both, or an HBB deletion in one and an HBB finding in the other; del(GJB6-D13S1830) counts with GJB2 and GJB6), the shared-risk heading says read depth adds a possible shared risk and names the gene and what each partner carries, so the heading never reads as nothing shared. The couple table and each partner's copy stay as they were. Step 39 passes both when a partner's folder holds them (`paralogs/<sample>_cn_intervals.tsv` of step 35, `cnvpytor/<sample>_cnvs.vcf.gz` of step 18), and the Nextflow couple screen passes them when `--tools` has `parascopy` and `cnvpytor`.
 
 ### The carrier list
 
diff --git a/docs/cn-loci.md b/docs/cn-loci.md
index ca64e612..1d86b9ad 100644
--- a/docs/cn-loci.md
+++ b/docs/cn-loci.md
@@ -72,7 +72,7 @@ Step 35 (and the PARASCOPY module) runs these three rows one after the other and
 - **CCL3L1**: the no-ALT reference holds one copy of CCL3L1/CCL3L3 (CCL3L3), so every read of either lands there and the row counts both together. The count goes back to Gonzalez 2005 (PMID 15637236), but the disease links reported for it are disputed, so no reading is given. The test fails when a CCL3L1 reading names a disease.
 - **DEFB4** has no row. The beta-defensin repeat that carries it (DEFB4B to DEFB103B) lies in the paralog copies of SPAG11A's Parascopy model, so Parascopy owns it under the rule below. The test plants a DEFB4 row there and checks that the owner check catches it.
 
-The report does not show `<sample>_cn_dosage.tsv` yet; the table and the step log are where to read it.
+Both reports show `<sample>_cn_dosage.tsv` in step 35's section, one line per locus ([step 35](35-paralogs.md#in-the-reports)). `main.nf` does not pass it to the Nextflow `HTML_REPORT` yet, so there the report says the table is not in the output; step 24 reads it from the sample folder.
 
 ### The CYP21A2 row
 
@@ -98,7 +98,7 @@ A site is low when its estimate is more than 0.6 copies under the larger of the
 
 A row on chrX or chrY is read against the sample's own count there: its controls sit on the same chromosome, so a male's one X or Y gives 1 copy, and `expected` says so. The sex is the declared one: step 35's second argument, and in Nextflow the samplesheet's sex column. Step 16's inferred sex is not used, so both give the same table. A chrY row on a female is `not applicable`; without a sex, chrX and chrY rows are `not assessed` (`sex_unknown`). Every row needs 8x of control depth (4x on chrY), flanks within 1.5x of each other, and an estimate within 0.3 of a whole number; a control with no read makes that row `not assessed` (`broken_control`) and the others still run. A failed bedcov pass writes every row as `not assessed` (`check_failed`).
 
-The couple screen (step 39) reads the table with `couple_report.py screen --cn-intervals ID=<id>_cn_intervals.tsv`: each carrier allele is present, absent or not assessable for each partner, and a male partner's AZF deletions are listed the same way (present means the deletion is there). With `--cnvpytor ID=<id>_cnvs.vcf.gz` (step 18), a CNVpytor call that disagrees with the depth call is printed beside it. Step 39 and the Nextflow couple screen pass both when the partner has them.
+The couple screen (step 39) reads the table with `couple_report.py screen --cn-intervals ID=<id>_cn_intervals.tsv`: each carrier allele is present, absent or not assessable for each partner, and a male partner's AZF deletions are listed the same way (present means the deletion is there). With `--cnvpytor ID=<id>_cnvs.vcf.gz` (step 18), a CNVpytor call that disagrees with the depth call is printed beside it. Step 39 and the Nextflow couple screen pass both when the partner has them. Both reports show the table in step 35's section: the counts by status and each row with a deletion or a duplication.
 
 To change a row, change the script (its hand rows or a pinned download) and run `scripts/ci/build-cn-intervals.py`; `--check` fails when the catalog is not what the downloads give.
 
diff --git a/scripts/ci/e2e-case-paths.tsv b/scripts/ci/e2e-case-paths.tsv
index 3bf941d0..270c2d13 100644
--- a/scripts/ci/e2e-case-paths.tsv
+++ b/scripts/ci/e2e-case-paths.tsv
@@ -74,6 +74,7 @@ clinvar-cpsr-quality-1	assets/construct_variants.tsv assets/deep_intronic_sites.
 cn-engine-1-parity	assets/bloodtype_markers.tsv bin/bloodtype_report.py bin/cn_loci.tsv bin/copy_number.py bin/hba_copy_number.py bin/lpa_kiv2.py scripts/41-bloodtype.sh scripts/46-lpa.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh
 cn-intervals-1-catalog	bin/cn_loci.tsv scripts/ci/build-cn-intervals.py scripts/ci/fake-docker/lib.sh
 cn-intervals-2-hg002	assets/par_grch38.bed assets/paralog_loci.tsv bin/cn_loci.tsv bin/copy_number.py bin/hba_copy_number.py bin/lpa_kiv2.py modules/local/deepvariant/* modules/local/parascopy/* modules/local/vcf_precheck/* scripts/35-paralogs.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh
+cn-report-1-tables	assets/carrier_genes.tsv assets/par_grch38.bed bin/collect_summary.py bin/couple_report.py bin/family_report.py bin/paraphase_report.py bin/render_report.py bin/results_index.py scripts/24-html-report.sh scripts/39-couple.sh scripts/ci/fake-docker/lib.sh scripts/generate-report.sh
 comphet-phasing-1-reads	assets/construct_variants.tsv assets/stub/sample.bam assets/stub/sample.bam.bai bin/clip_check.py bin/comphet_phase.py bin/constraint_join.awk bin/csq_rank.awk modules/local/slivar/* modules/local/vcf_precheck/* modules/local/vep/* scripts/13-vep-annotation.sh scripts/31-slivar.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh tests/fixtures/comphet/make_fixture.py
 construct-scan-1-both	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/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/vcf_precheck/* scripts/06-clinvar-screen.sh scripts/24-html-report.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh
 consumer-export-1	assets/chip_genes.tsv assets/par_grch38.bed bin/chip_report.py bin/export_consumer.py modules/local/chip/* modules/local/export_consumer/* modules/local/vcf_precheck/* scripts/40-export-consumer.sh scripts/49-chip.sh scripts/chip-to-vcf.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh tests/fixtures/consumer_export/sites.tsv.gz
diff --git a/scripts/ci/e2e-case-times.tsv b/scripts/ci/e2e-case-times.tsv
index 96e74ca1..a73b4f5c 100644
--- a/scripts/ci/e2e-case-times.tsv
+++ b/scripts/ci/e2e-case-times.tsv
@@ -73,6 +73,7 @@ clinvar-cpsr-quality-1	shard	5	clinvar-clip-check-1
 cn-engine-1-parity	shard	33	-
 cn-intervals-1-catalog	shard	0	-
 cn-intervals-2-hg002	shard	19	-
+cn-report-1-tables	shard	60	-
 comphet-phasing-1-reads	shard	16	clinical-filter-mane-and-quality-1-tiers
 construct-scan-1-both	shard	42	-
 consumer-export-1	shard	59	-
diff --git a/tests/e2e/cn-report-1-tables.sh b/tests/e2e/cn-report-1-tables.sh
new file mode 100644
index 00000000..42d98ac3
--- /dev/null
+++ b/tests/e2e/cn-report-1-tables.sh
@@ -0,0 +1,121 @@
+#!/usr/bin/env bash
+# Step 35's tables that no Parascopy model makes reach both reports and the
+# couple screen, with the real images of steps 24, generate-report.sh and
+# 39, on made-up tables in a GENOME_DIR of their own:
+#   1. steps 24 and generate-report.sh on a sample with an SMN table, a
+#      CYP21A2 table (a chimeric deletion: 1 copy, module total 3 of 4), the
+#      gene-dosage table (HP 3, LCE3BC 1, CCL3L1 not assessed) and an interval
+#      table with one 22q11.2 deletion: both reports show each of them, and
+#      neither says CYP21A2 is not assessed. Red control: with the CYP21A2
+#      table gone the text report says it is not in the output, with no
+#      finding.
+#   2. step 39 on a made-up couple, CNR_A (female) and CNR_B (male), each
+#      with an empty ClinVar screen, an interval table with CFTRdele2,3
+#      present and a CYP21A2 table (CNR_A the deletion, CNR_B none): step 39
+#      passes both CYP21A2 tables, the couple table has CNR_A's CYP21A2
+#      deletion, and the shared-risk heading names CFTRdele2,3 from read depth.
+#      Red control: with CNR_A's CYP21A2 table gone there is no CYP21A2 row
+#      and the report names the missing table.
+. "$(dirname "$0")/lib.sh"
+
+G="${CASE_TMP}/genome"
+T="${SAMPLE}cnr"
+rm -rf "$G"
+mkdir -p "${G}/${T}/paralogs"
+printf 'run\tscript\te2e\n' > "${G}/${T}/run_manifest.tsv"
+P="${G}/${T}/paralogs"
+printf 'chrom\tstart\tend\tlocus\tagCN_filter\tagCN\tagCN_qual\tpsCN_filter\tpsCN\tpsCN_qual\thomologous_regions\n%s\n' \
+  "chr5	70950920	70953320	SMN1	PASS	4	60	PASS	2,2	33,30	chr5:70075496-70077900:+" > "${P}/${T}_smn_copy_number.tsv"
+cyp_del() {
+  printf 'metric\tvalue\nstatus\tcalled\ncyp21a2_copies\t1\nmodule_total\t3\nexpected_total\t4\nestimate\t3.02\n'
+  printf 'finding\tdeletion\nconversion_sites\t.\ntarget_depth\t45.30\ncontrol_depth\t30.10\nreason\t.\n'
+  printf 'reading\tthe pattern of the 30 kb CYP21A1P/CYP21A2 chimeric deletion\n'
+}
+cyp_none() {
+  printf 'metric\tvalue\nstatus\tcalled\ncyp21a2_copies\t2\nmodule_total\t4\nexpected_total\t4\nestimate\t4.01\n'
+  printf 'finding\tnone\nconversion_sites\t.\ntarget_depth\t60.20\ncontrol_depth\t30.10\nreason\t.\nreading\t.\n'
+}
+cyp_del > "${P}/${T}_cyp21a2_psv.tsv"
+{
+  printf 'locus\tstatus\tcopies\testimate\treason\treading\n'
+  printf 'HP\tcalled\t3\t3.04\t.\tthree copies of the duplicated HP exons: Hp2-1\n'
+  printf 'LCE3BC\tcalled\t1\t0.98\t.\tLCE3B and LCE3C deleted on one chromosome\n'
+  printf 'CCL3L1\tnot assessed\t.\t.\tcheck_failed\tthe depth check failed; see the log of step 35\n'
+} > "${P}/${T}_cn_dosage.tsv"
+ivh() { printf 'locus\tname\tchrom\tstart\tend\tstatus\tcopies\texpected\tchange\testimate\treason\n'; }
+{
+  ivh
+  printf '22q11.2_LCR_A-LCR_D\t22q11.2 recurrent region (A-D)\tchr22\t18924717\t21111383\tcalled\t1\t2\tdeletion\t1.02\t.\n'
+  printf 'CFTR_dele2_3\tCFTRdele2,3\tchr7\t117498312\t117519392\tcalled\t2\t2\tno change\t1.97\t.\n'
+} > "${P}/${T}_cn_intervals.tsv"
+
+# --- 1. both reports ----------------------------------------------------------------
+GENOME_DIR="$G" run_step 24-html-report.sh "$T"
+check_step_exit 24-html-report.sh
+HTML=$(cat "${G}/${T}/${T}_report.html" 2>/dev/null)
+check "HTML: the CYP21A2 deletion" \
+  has '>CYP21A2 \(RCCX module, from depth and six PSV sites\)<.*1 CYP21A2 copies, module total 3 of 4, finding deletion' "$HTML"
+check "HTML: the HP count" has '>HP \(haptoglobin\) copies \(from depth\)<.*3 \(estimate 3.04\)' "$HTML"
+check "HTML: the 22q11.2 deletion row" \
+  has '<td>22q11.2_LCR_A-LCR_D</td><td>chr22:18924718-21111383</td><td>1 of 2</td><td>deletion</td>' "$HTML"
+check "HTML: no claim that CYP21A2 is not assessed" lacks 'CYP21A2 has no Parascopy model and is not assessed' "$HTML"
+
+GENOME_DIR="$G" run_step generate-report.sh "$T"
+check_step_exit generate-report.sh
+TXT=$(cat "${G}/${T}/${T}_report.txt" 2>/dev/null)
+sed -n '/SMN1\/SMN2 Copy Number (Parascopy/,/^$/p' <<< "$TXT"
+check "text: the CYP21A2 deletion" \
+  has '^  CYP21A2 \(RCCX module, from depth and six PSV sites\): 1 CYP21A2 copies, module total 3 of 4, finding deletion' "$TXT"
+check "text: the three dosage loci" \
+  bash -c "grep -q '^  HP (haptoglobin) copies (from depth): 3' <<< \"\$1\" && grep -q '^  LCE3B/LCE3C copies (from depth): 1' <<< \"\$1\" \
+    && grep -q '^  CCL3L1/CCL3L3 copies (from depth): not assessed (check_failed)' <<< \"\$1\"" _ "$TXT"
+check "text: the interval counts and the deletion" \
+  bash -c "grep -q 'from depth): 2 rows (2 called); 1 with a deletion or a duplication' <<< \"\$1\" \
+    && grep -q '^    Deletion: 22q11.2 recurrent region (A-D) (22q11.2_LCR_A-LCR_D, chr22:18924718-21111383), 1 of 2 copies' <<< \"\$1\"" _ "$TXT"
+
+mv "${P}/${T}_cyp21a2_psv.tsv" "${CASE_TMP}/cyp.tsv"
+GENOME_DIR="$G" run_step generate-report.sh "$T"
+TXT0=$(cat "${G}/${T}/${T}_report.txt" 2>/dev/null)
+check "red: without the CYP21A2 table the text report says it is not in the output" \
+  has "CYP21A2 \(RCCX module, from depth and six PSV sites\): not assessed \(step 35's <id>_cyp21a2_psv.tsv is not in this run's output\)" "$TXT0"
+check "red: and gives no CYP21A2 finding" lacks 'finding deletion' "$TXT0"
+
+# --- 2. step 39 -------------------------------------------------------------------
+C="${G}/couple-case"
+A=CNR_A; B=CNR_B
+for s in "$A" "$B"; do
+  mkdir -p "${C}/${s}/clinvar" "${C}/${s}/paralogs"
+  printf '##fileformat=VCFv4.2\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t%s\n' "$s" \
+    > "${C}/${s}/clinvar/${s}_clinvar_hits.vcf"
+  {
+    ivh
+    printf 'CFTR_dele2_3\tCFTRdele2,3\tchr7\t117498312\t117519392\tcalled\t1\t2\tdeletion\t1.01\t.\n'
+  } > "${C}/${s}/paralogs/${s}_cn_intervals.tsv"
+done
+cyp_del > "${C}/${A}/paralogs/${A}_cyp21a2_psv.tsv"
+cyp_none > "${C}/${B}/paralogs/${B}_cyp21a2_psv.tsv"
+printf 'sample\tsex\tfather\tmother\tpartner\n%s\tfemale\t.\t.\t%s\n%s\tmale\t.\t.\t%s\n' "$A" "$B" "$B" "$A" > "${C}/pedigree.tsv"
+
+GENOME_DIR="$C" run_step 39-couple.sh couple "${C}/pedigree.tsv"
+check_step_exit 39-couple.sh
+check "step 39 passes both CYP21A2 tables" \
+  bash -c "grep -q '${A}: paralogs/${A}_cyp21a2_psv.tsv (step 35)' '$STEP_LOG' && grep -q '${B}: paralogs/${B}_cyp21a2_psv.tsv (step 35)' '$STEP_LOG'"
+TSV=$(cat "${C}/couple/couple_${A}_${B}.tsv" 2>/dev/null)
+MD=$(cat "${C}/couple/couple_${A}_${B}.md" 2>/dev/null)
+check "step 39: CNR_A's CYP21A2 deletion is a row of the couple table" \
+  has "^CYP21A2	.*one_partner	.*CYP21A2 deletion from depth and PSV sites \(CYP21A2 1 copies, RCCX module total 3 of 4\)" "$TSV"
+check "step 39: the shared-risk heading names read depth" \
+  has "No gene of the list where both partners carry a change in the gene table; read depth \(step 35\) adds a possible shared risk, below\." "$MD"
+check "step 39: CFTRdele2,3 present in both is named" \
+  has "^- \*\*CFTR\*\*: ${A}: CFTRdele2,3 present \(read depth, step 35: 1 of 2 copies\); ${B}: CFTRdele2,3 present" "$MD"
+
+mv "${C}/${A}/paralogs/${A}_cyp21a2_psv.tsv" "${CASE_TMP}/a_cyp.tsv"
+GENOME_DIR="$C" run_step 39-couple.sh couple "${C}/pedigree.tsv"
+check_step_exit 39-couple.sh
+check "red: without CNR_A's CYP21A2 table no CYP21A2 row" lacks "^CYP21A2	" "$(cat "${C}/couple/couple_${A}_${B}.tsv" 2>/dev/null)"
+check "red: and the report names the missing table" \
+  has "CYP21A2 deletion and conversions of ${A}: step 35's CYP21A2 table was not passed to the couple screen \(--cyp21a2\)\." \
+  "$(cat "${C}/couple/couple_${A}_${B}.md" 2>/dev/null)"
+docker run --rm -v "${CASE_TMP}:/w" "$PYTHON_IMAGE" rm -rf /w/genome
+
+finish

From d7856cc4974f68467d6c26730d2cb7b77fb3bbe6 Mon Sep 17 00:00:00 2001
From: GeiserX <9169332+GeiserX@users.noreply.github.com>
Date: Fri, 6 Nov 2026 22:24:16 +0100
Subject: [PATCH 4/7] test(cn): make the cn-report e2e case executable

Bead: pgp-5bja.21
---
 tests/e2e/cn-report-1-tables.sh | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 mode change 100644 => 100755 tests/e2e/cn-report-1-tables.sh

diff --git a/tests/e2e/cn-report-1-tables.sh b/tests/e2e/cn-report-1-tables.sh
old mode 100644
new mode 100755

From ab1e98d48b0be839872e01ce7b9acfeba118a47c Mon Sep 17 00:00:00 2001
From: GeiserX <9169332+GeiserX@users.noreply.github.com>
Date: Fri, 6 Nov 2026 23:36:45 +0100
Subject: [PATCH 5/7] fix(cn): an interval table with no row called says none
 called, not no deletion

A failed bedcov pass leaves every interval row not assessed, and the
summary line then read as a negative result. It now ends with 'none
called', and a called table names how many rows the negative covers.
The docs name which header-only tables are unreadable: a CYP21A2 table
without a status row reads as not assessed.
---
 bin/render_report.py          |  9 ++++++---
 docs/35-paralogs.md           |  4 ++--
 tests/test_collect_summary.py | 10 +++++++++-
 3 files changed, 17 insertions(+), 6 deletions(-)

diff --git a/bin/render_report.py b/bin/render_report.py
index 8cc7b84c..c70f962b 100755
--- a/bin/render_report.py
+++ b/bin/render_report.py
@@ -1743,9 +1743,12 @@ def cn_interval_summary(iv):
         return f"{CN_INTERVALS_LABEL}: not assessed (step 35's <id>_cn_intervals.tsv is not in this run's output)"
     counts = ", ".join(f"{n} {k}" for k, n in sorted((iv.get("by_status") or {}).items())) or "no row"
     n = len(iv.get("changes") or [])
-    return (f"{CN_INTERVALS_LABEL}: {iv.get('rows', 0)} rows ({counts}); "
-            f"{n} with a deletion or a duplication" if n else
-            f"{CN_INTERVALS_LABEL}: {iv.get('rows', 0)} rows ({counts}); no deletion or duplication")
+    k = (iv.get("by_status") or {}).get("called", 0)
+    # Rows that were not called are no negative result: with none called the
+    # line says so instead of "no deletion or duplication".
+    end = (f"{n} with a deletion or a duplication" if n else
+           f"no deletion or duplication among the {k} called rows" if k else "none called")
+    return f"{CN_INTERVALS_LABEL}: {iv.get('rows', 0)} rows ({counts}); {end}"
 
 
 def cn_interval_place(r):
diff --git a/docs/35-paralogs.md b/docs/35-paralogs.md
index 6f1ad846..5262c197 100644
--- a/docs/35-paralogs.md
+++ b/docs/35-paralogs.md
@@ -236,9 +236,9 @@ The same section shows the tables no Parascopy model makes, each read on its own
 
 - `${SAMPLE}_cyp21a2_psv.tsv` as one line, `CYP21A2 (RCCX module, from depth and six PSV sites)`: the CYP21A2 copies, the module total of 4, the finding (with the sites of a possible conversion) and its reading, or the reason it was not assessed. The run summary carries it as the Parascopy section's `cyp21a2`.
 - `${SAMPLE}_cn_dosage.tsv` as one line per locus (HP, LCE3B/LCE3C, CCL3L1/CCL3L3): the copies, the estimate and the reading, or the status and its reason. The run summary carries it as `dosage`, one entry per locus.
-- `${SAMPLE}_cn_intervals.tsv` as one line of counts (rows by status, how many show a deletion or a duplication), then each row with a change: its region name, catalog row, place (1-based) and copies of expected. In the HTML card the changed rows are a table. The run summary carries it as `cn_intervals`.
+- `${SAMPLE}_cn_intervals.tsv` as one line of counts (rows by status, how many show a deletion or a duplication, or `none called` when no row was called), then each row with a change: its region name, catalog row, place (1-based) and copies of expected. In the HTML card the changed rows are a table. The run summary carries it as `cn_intervals`.
 
-A table that is not in the run's output reads `not assessed (step 35's <id>_... is not in this run's output)`, so a missing table never reads as a normal result. A table with a header but without its columns makes the section unreadable; a catalog table without the copy-number columns makes that locus unreadable, also with no row under the header.
+A table that is not in the run's output reads `not assessed (step 35's <id>_... is not in this run's output)`, so a missing table never reads as a normal result. An SMN, interval or dosage table with a header but without its columns makes the section unreadable, and a CYP21A2 table without a status row reads as not assessed; a catalog table without the copy-number columns makes that locus unreadable, also with no row under the header.
 
 The couple screen reads the CYP21A2 table with `couple_report.py screen --cyp21a2 ID=<id>_cyp21a2_psv.tsv`, where a deletion or a possible conversion is a CYP21A2 carrier finding. Step 39 passes it when a partner's folder holds it. The Nextflow `PARASCOPY` task writes the CYP21A2 and dosage tables, but `main.nf` does not pass them to `HTML_REPORT` or the couple screen yet, so a Nextflow run's reports say they are not in the output.
 
diff --git a/tests/test_collect_summary.py b/tests/test_collect_summary.py
index a5c90df5..d8cb7b8f 100644
--- a/tests/test_collect_summary.py
+++ b/tests/test_collect_summary.py
@@ -1486,7 +1486,15 @@ def card(html):
         + "CFTR_dele2_3\tCFTRdele2,3\tchr7\t117479000\t117481000\tcalled\t2\t2\tno change\t1.97\t.\n")
     _, txt, _ = render(d)
     check("cn report: an interval table without a change says so",
-          "from depth): 1 rows (1 called); no deletion or duplication\n" in txt, txt)
+          "from depth): 1 rows (1 called); no deletion or duplication among the 1 called rows\n" in txt, txt)
+    put(f"{d}/paralogs/S_cn_intervals.tsv", IV_COLS
+        + "CFTR_dele2_3\tCFTRdele2,3\tchr7\t117479000\t117481000\tnot assessed\t.\t2\t.\t.\tcheck_failed\n"
+        + "AZFc\tAZFc\tchrY\t24000000\t26000000\tnot assessed\t.\t1\t.\t.\tcheck_failed\n")
+    _, txt, html = render(d)
+    check("cn report: an interval table with no row called says none called, never no deletion",
+          "from depth): 2 rows (2 not assessed); none called\n" in txt
+          and "2 rows (2 not assessed); none called" in card(html)
+          and "no deletion or duplication" not in txt, txt)
 
     # A header without the copy-number columns is unreadable, also with no row under it.
     put(f"{d}/paralogs/S_cn_dosage.tsv", "locus\tstatus\n")

From 6723b87a6f588c089bc12c1bbe7f6ad72b567e5d Mon Sep 17 00:00:00 2001
From: GeiserX <9169332+GeiserX@users.noreply.github.com>
Date: Sat, 7 Nov 2026 00:59:18 +0100
Subject: [PATCH 6/7] fix(cn): reject cut-short dosage and interval rows, keep
 the couple depth risk to the carrier list

A dosage or interval row cut short of a column read as a row of '.'
values, or crashed both reports on a None status; such a table is now
unreadable. Step 35 removes its earlier HBA, dosage, CYP21A2 and interval
tables before the new SMN tables land, so a stopped step leaves none that
reads as this run's. The read-depth shared risk of the couple screen now
skips a gene the carrier list leaves out, and the sample reports say their
counts are of the gene table only. Tests for a DFNB1 pair and a folder
with only the dosage table.
---
 bin/collect_summary.py        | 16 ++++++++++++++--
 bin/couple_report.py          |  9 ++++++---
 bin/render_report.py          |  5 +++--
 docs/39-couple.md             |  2 +-
 scripts/35-paralogs.sh        |  5 +++++
 tests/test_collect_summary.py | 18 ++++++++++++++++++
 tests/test_couple_report.py   | 20 ++++++++++++++++++++
 7 files changed, 67 insertions(+), 8 deletions(-)

diff --git a/bin/collect_summary.py b/bin/collect_summary.py
index 9534329b..0caf4f4b 100755
--- a/bin/collect_summary.py
+++ b/bin/collect_summary.py
@@ -905,8 +905,10 @@ def read_cn_intervals(path):
     """Step 35's interval table (bin/copy_number.py intervals): the rows with
     a deletion or a duplication, and a count of the rows by status."""
     rows = read_tsv(path)
-    if not {"locus", "name", "chrom", "start", "end", "status", "copies", "expected", "change"} <= tsv_columns(path):
+    cols = ("locus", "name", "chrom", "start", "end", "status", "copies", "expected", "change")
+    if not set(cols) <= tsv_columns(path):
         raise Unreadable(f"{os.path.basename(path)} lacks the locus, name, status, copies, expected or change column")
+    short_row(path, rows, cols)
     counts = {}
     for r in rows:
         counts[r["status"]] = counts.get(r["status"], 0) + 1
@@ -923,7 +925,17 @@ def read_cn_dosage(path):
     in the table's order."""
     if not set(DOSAGE_COLUMNS) <= tsv_columns(path):
         raise Unreadable(f"{os.path.basename(path)} lacks the {', '.join(DOSAGE_COLUMNS[:-1])} or reading column")
-    return [{k: r.get(k) or "." for k in DOSAGE_COLUMNS} for r in read_tsv(path)]
+    rows = read_tsv(path)
+    short_row(path, rows, DOSAGE_COLUMNS)
+    return [{k: r.get(k) or "." for k in DOSAGE_COLUMNS} for r in rows]
+
+
+def short_row(path, rows, cols):
+    """A row cut short of a column (csv gives None there) makes the table
+    unreadable, so it never reads as a row of '.' values."""
+    bad = next((i for i, r in enumerate(rows, 2) if any(r.get(k) is None for k in cols)), None)
+    if bad:
+        raise Unreadable(f"{os.path.basename(path)} line {bad} lacks a field")
 
 
 def read_hba(path):
diff --git a/bin/couple_report.py b/bin/couple_report.py
index 8b830b32..1096de13 100755
--- a/bin/couple_report.py
+++ b/bin/couple_report.py
@@ -519,16 +519,19 @@ def cnv_note(row, dels):
     return ""
 
 
-def cn_shared(a, b, intervals, found):
+def cn_shared(a, b, intervals, found, carrier):
     """The genes where both partners carry a change once the read-depth
     alleles count: one has a present allele of CN_ALLELES and the other a
     present allele or a finding in the same gene. The gene table leaves
     these alleles out, so without this a couple who both carry CFTRdele2,3
     would read as sharing nothing. intervals: {partner: {locus: row}} of the
-    partners with an interval table."""
+    partners with an interval table; carrier: the genes of the carrier list,
+    so a gene the list leaves out is no shared risk here either."""
     names = dict(CN_ALLELES)
     out = []
     for label, genes, loci in CN_ALLELE_GENES:
+        if not any(g in carrier for g in genes):
+            continue
         have, by_depth = {}, False
         for s in (a, b):
             iv = intervals.get(s) or {}
@@ -849,7 +852,7 @@ def cmd_screen(a):
         write_table(base + ".tsv", rows)
         for s in (x, y):
             write_table(os.path.join(a.out_dir, "%s_couple_carrier.tsv" % s), rows)
-        shared_cn = cn_shared(x, y, ivs, found)
+        shared_cn = cn_shared(x, y, ivs, found, genes)
         write_markdown(base + ".md", x, y, sexes, rows, assessed, gaps, len(genes), common, cn, shared_cn)
         n = {k: sum(1 for r in rows if r["result"] == k) for k in ("shared", "x_linked", "one_partner")}
         print("%s and %s: %d shared gene%s, %d X-linked, %d with one partner only (%s.tsv)" % (
diff --git a/bin/render_report.py b/bin/render_report.py
index c70f962b..585006ac 100755
--- a/bin/render_report.py
+++ b/bin/render_report.py
@@ -85,8 +85,9 @@ def ancient_lines(d):
 
 
 COUPLE_NOTE = ("A screen, not a diagnosis: the couple report (couple_<a>_<b>.md in the family folder) lists the "
-               "genes, variant types and repeats it does not see. See a genetic counsellor before acting on a "
-               "result, or on its absence (docs/39-couple.md).")
+               "genes, variant types and repeats it does not see. These counts are of its gene table: a carrier "
+               "deletion allele from read depth (step 35), and a possible shared risk it adds, is in that report "
+               "only. See a genetic counsellor before acting on a result, or on its absence (docs/39-couple.md).")
 PRS_ANCESTRY_NOTE = ("Most PRS scores were developed in European-ancestry populations and predict less well for "
                      "other ancestries, even with an ancestry-adjusted percentile.")
 
diff --git a/docs/39-couple.md b/docs/39-couple.md
index cd2d2b95..6402c026 100644
--- a/docs/39-couple.md
+++ b/docs/39-couple.md
@@ -39,7 +39,7 @@ The table of what was assessed counts, for each partner, the deletions over a li
 
 Only the ClinVar screen is required. Each of the others is used when its file is there, and the report lists every one that was not.
 
-`couple_report.py screen` also takes step 35's interval table (`--cn-intervals ID=paralogs/<sample>_cn_intervals.tsv`, [the interval rows](cn-loci.md#the-interval-rows)). The report then lists, for each partner, del(GJB6-D13S1830), CFTRdele2,3, the HBB exon 3 deletion, the HBB 619 bp deletion and the Sicilian (delta-beta)0 deletion as present, absent or not assessable from read depth, and a male partner's AZFa, AZFb and AZFc deletions (Y chromosome microdeletions, a cause of male infertility) the same way, where present means the deletion is there; a female partner's AZF line says not applicable. The three beta-globin rows overlap (the Sicilian deletion covers both HBB rows, and the 619 bp and exon 3 rows share 487 bp), so one deletion can read as present on two or three of them; the report says so, and only a clinical test names the allele. With `--cnvpytor ID=cnvpytor/<sample>_cnvs.vcf.gz` (step 18), a CNVpytor call that disagrees with the depth call (a deletion over half the allele's span or more where depth sees none, or none where depth sees one) is printed beside it. A missing call is not printed for an allele shorter than 2 kb (the two HBB rows): step 18 counts in 1 kb bins and cannot see a deletion that short. These lines are a list, not carrier findings: they do not enter the gene table. When both partners carry something in one gene once these alleles count (CFTRdele2,3 present in both, or an HBB deletion in one and an HBB finding in the other; del(GJB6-D13S1830) counts with GJB2 and GJB6), the shared-risk heading says read depth adds a possible shared risk and names the gene and what each partner carries, so the heading never reads as nothing shared. The couple table and each partner's copy stay as they were. Step 39 passes both when a partner's folder holds them (`paralogs/<sample>_cn_intervals.tsv` of step 35, `cnvpytor/<sample>_cnvs.vcf.gz` of step 18), and the Nextflow couple screen passes them when `--tools` has `parascopy` and `cnvpytor`.
+`couple_report.py screen` also takes step 35's interval table (`--cn-intervals ID=paralogs/<sample>_cn_intervals.tsv`, [the interval rows](cn-loci.md#the-interval-rows)). The report then lists, for each partner, del(GJB6-D13S1830), CFTRdele2,3, the HBB exon 3 deletion, the HBB 619 bp deletion and the Sicilian (delta-beta)0 deletion as present, absent or not assessable from read depth, and a male partner's AZFa, AZFb and AZFc deletions (Y chromosome microdeletions, a cause of male infertility) the same way, where present means the deletion is there; a female partner's AZF line says not applicable. The three beta-globin rows overlap (the Sicilian deletion covers both HBB rows, and the 619 bp and exon 3 rows share 487 bp), so one deletion can read as present on two or three of them; the report says so, and only a clinical test names the allele. With `--cnvpytor ID=cnvpytor/<sample>_cnvs.vcf.gz` (step 18), a CNVpytor call that disagrees with the depth call (a deletion over half the allele's span or more where depth sees none, or none where depth sees one) is printed beside it. A missing call is not printed for an allele shorter than 2 kb (the two HBB rows): step 18 counts in 1 kb bins and cannot see a deletion that short. These lines are a list, not carrier findings: they do not enter the gene table. When both partners carry something in one gene once these alleles count (CFTRdele2,3 present in both, or an HBB deletion in one and an HBB finding in the other; del(GJB6-D13S1830) counts with GJB2 and GJB6), the shared-risk heading says read depth adds a possible shared risk and names the gene and what each partner carries, so the heading never reads as nothing shared. A gene the carrier list (`--genes`) leaves out gives no such line. The couple table and each partner's copy stay as they were, so the sample reports count the gene table only and say the read-depth alleles are in the couple report. Step 39 passes both when a partner's folder holds them (`paralogs/<sample>_cn_intervals.tsv` of step 35, `cnvpytor/<sample>_cnvs.vcf.gz` of step 18), and the Nextflow couple screen passes them when `--tools` has `parascopy` and `cnvpytor`.
 
 ### The carrier list
 
diff --git a/scripts/35-paralogs.sh b/scripts/35-paralogs.sh
index 31e8b90e..350a8571 100755
--- a/scripts/35-paralogs.sh
+++ b/scripts/35-paralogs.sh
@@ -157,6 +157,11 @@ echo "Models: ${POP} (${MODELS}), loci: ${LOCI}"
 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" \
+  "${OUTDIR}/${SAMPLE}_cyp21a2_psv.tsv" "${OUTDIR}/${SAMPLE}_cn_intervals.tsv"
 mkdir -p "$WORK"
 # Background depth once, then the copy number of each catalog row
 run_in --cpus "${THREADS}" --memory 8g "${PARASCOPY_IMAGE}" \
diff --git a/tests/test_collect_summary.py b/tests/test_collect_summary.py
index d8cb7b8f..501aa928 100644
--- a/tests/test_collect_summary.py
+++ b/tests/test_collect_summary.py
@@ -1420,6 +1420,14 @@ def cn_report_section(work, schema):
     """Step 35's tables that no Parascopy model makes reach both reports: the
     CYP21A2 call, the gene-dosage rows and the interval rows; and a table with
     a header but without the copy-number columns reads as unreadable."""
+    only = os.path.join(work, "cn_dosage_only", "S")
+    put(f"{only}/paralogs/S_cn_dosage.tsv", "locus\tstatus\tcopies\testimate\treason\treading\n"
+        "HP\tcalled\t3\t3.04\t.\tthree copies\n")
+    summ, txt, html = render(only)
+    dosage = ((summ["sections"].get("parascopy") or {}).get("data") or {}).get("dosage") or []
+    check("cn report: a folder with only the dosage table still shows it in the summary and both reports",
+          [x["copies"] for x in dosage] == ["3"] and "HP (haptoglobin) copies (from depth): 3 (estimate 3.04)" in txt
+          and ">HP (haptoglobin) copies (from depth)<" in html, (dosage, txt))
     d = os.path.join(work, "cn_report", "S")
     put(f"{d}/paralogs/S_smn_copy_number.tsv", SMN_COLS
         + "chr5\t70950920\t70953320\tSMN1\tPASS\t4\t60\tPASS\t2,2\t33,30\tchr5:70075496-70077900:+\n")
@@ -1496,6 +1504,16 @@ def card(html):
           and "2 rows (2 not assessed); none called" in card(html)
           and "no deletion or duplication" not in txt, txt)
 
+    # A row cut short of a column is unreadable, never a row of '.' values.
+    put(f"{d}/paralogs/S_cn_dosage.tsv", "locus\tstatus\tcopies\testimate\treason\treading\nHP\tcalled\n")
+    ps = collect_summary.collect("S", d)["sections"].get("parascopy") or {}
+    check("cn report: a dosage row cut short of its columns is unreadable", ps.get("state") == "unreadable", ps)
+    os.remove(f"{d}/paralogs/S_cn_dosage.tsv")
+    put(f"{d}/paralogs/S_cn_intervals.tsv", IV_COLS
+        + "CFTR_dele2_3\tCFTRdele2,3\tchr7\t117479000\t117481000\tcalled\t2\t2\tno change\t1.97\t.\n"
+        + "AZFc\tAZFc\n")
+    ps = collect_summary.collect("S", d)["sections"].get("parascopy") or {}
+    check("cn report: an interval row cut short before its status is unreadable", ps.get("state") == "unreadable", ps)
     # A header without the copy-number columns is unreadable, also with no row under it.
     put(f"{d}/paralogs/S_cn_dosage.tsv", "locus\tstatus\n")
     ps = collect_summary.collect("S", d)["sections"].get("parascopy") or {}
diff --git a/tests/test_couple_report.py b/tests/test_couple_report.py
index e6f0ca90..60bc5b4b 100644
--- a/tests/test_couple_report.py
+++ b/tests/test_couple_report.py
@@ -482,6 +482,26 @@ def screen2(civ, div, out):
         check("cn shared, red: with C's CFTRdele2,3 absent and no HBB deletion the plain negative line stays",
               rc == 0 and "No gene of the list where both partners carry a change.\n" in head
               and "Possible shared risk" not in head and "possibly shared" not in log, head)
+        # DFNB1: one partner's del(GJB6-D13S1830) beside the other's GJB2 hit
+        g_hits = put(f"{t}/G.vcf", HEAD + hit("GJB2"))
+        g_iv = put(f"{t}/G_cn_intervals.tsv", ivh + ivrow("GJB6_D13S1830", "chr13", 20223037, 20531806, copies="1",
+                                                          change="deletion"))
+        rc, log = run("screen", "--pedigree", ped2, "--genes", GENES, "--clinvar", f"C={c_hits}", "--clinvar",
+                      f"D={g_hits}", "--cn-intervals", f"C={g_iv}", "--out-dir", f"{t}/cn4")
+        md = open(f"{t}/cn4/couple_C_D.md").read() if os.path.exists(f"{t}/cn4/couple_C_D.md") else ""
+        head = md.split("## Shared risk", 1)[-1].split("## X-linked", 1)[0]
+        check("cn shared: C's del(GJB6-D13S1830) beside D's GJB2 hit is a possible DFNB1 shared risk naming both",
+              rc == 0 and "- **GJB2/GJB6 (DFNB1)**: C: del(GJB6-D13S1830) present (read depth, step 35: 1 of 2 "
+              "copies); D: " in head and "GJB2" in head.split("; D: ", 1)[-1], (rc, log, head))
+        # A carrier list without CFTR: CFTRdele2,3 in both is no shared risk; HBB still is
+        no_cftr = put(f"{t}/no_cftr.tsv", "".join(ln for ln in open(GENES) if not ln.startswith("CFTR\t")))
+        rc, log = run("screen", "--pedigree", ped2, "--genes", no_cftr, "--clinvar", f"C={c_hits}", "--clinvar",
+                      f"D={d_hits}", "--cn-intervals", f"C={c_iv}", "--cn-intervals", f"D={d_iv}", "--out-dir", f"{t}/cn5")
+        md = open(f"{t}/cn5/couple_C_D.md").read() if os.path.exists(f"{t}/cn5/couple_C_D.md") else ""
+        head = md.split("## Shared risk", 1)[-1].split("## X-linked", 1)[0]
+        check("cn shared: a carrier list without CFTR gives no CFTR shared risk, HBB stays",
+              rc == 0 and "**CFTR**" not in head and "- **HBB**: C: HBB exon 3 deletion present" in head
+              and "possibly shared, from read depth: CFTR" not in log, (rc, log, head))
 
         # 10. genes
         gtf = put(f"{t}/g.gtf", "##description: evidence-based annotation of the human genome (GRCh38), version 99\n" + "".join(

From 4c19ad12654607e5b125666244968410f47364f7 Mon Sep 17 00:00:00 2001
From: GeiserX <sergio@geiser.cloud>
Date: Sat, 7 Nov 2026 10:43:51 +0100
Subject: [PATCH 7/7] fix(search): find repeat, Locityper and CFTR/UGT1A1
 readings, and keep Locityper settings on a rerun (#321)

## Problem

The results search and `owngenome-lookup.sh` could not find a repeat expansion, a Locityper locus or the CFTR, UGT1A1 and SLC6A4 readings by gene. Those sections name their gene under `locus`, and the search only read `gene`, `genes`, `Gene`, `SYMBOL` and `symbol`. A flagged HTT expansion was not in the index.

Three Locityper gaps sat beside it:

- `setup.sh` stamps the sha256 of the reference's `.fai` and FASTA into the database's `DONE`, but step 53 and the Nextflow `LOCITYPER` checked only the `.fai`. A FASTA with the same contigs and other bases passed.
- An install without `DONE` failed in the module with "built from another reference", which is the wrong reason.
- A result with no `total_reads` or `unexpl_reads` read as 0, so the unexplained-read rule passed on "0 of 0".

And `run-all.sh` read `LOCITYPER_DATA` and `LOCITYPER_BG_REGION` but left them out of `SWITCHES`. The recorded command dropped them, so `owngenome_recheck.py --rerun` replayed step 53 against the default database.

## Fix

- `bin/genome_lookup.py`: `locus` joins `GENE_KEYS`, behind the same upper-case filter. The index and the lookup share `names_of()`, so both find HTT, FMR1, UGT1A1, FCGR and HLA-DRB1 now. Regions sit under other keys, so no region becomes a term.
- `scripts/53-locityper.sh` and `modules/local/locityper/main.nf` check `fai=<sha256> fasta=<sha256> ` in `DONE`. The module first checks that `DONE` exists and says the install is unfinished when it does not. Hashing a GRCh38 FASTA adds about 10 to 20 seconds per run.
- `bin/locityper_report.py`: a result with no read counts fails the filter with the reason, and its count columns read `.`.
- `scripts/run-all.sh`: both switches are in `SWITCHES` and in the header. Two `switch` rows join `assets/output_contract.tsv`.

Output changes: the search index and the report's closing gene list gain terms from these sections. In `<sample>_locityper.tsv`, `total_reads` and `unexplained_reads` read `.` instead of `0` when Locityper gave no count, and such a call now fails the filter. Columns and file names are unchanged.

Decisions:

- Missing read counts fail the filter but keep the call, instead of marking the locus "not called". Locityper did produce a genotype. Failing the filter keeps it visible and stops it from passing as clean.
- The package's file list names `bin/results_search*.py`, which matches no file. The search keys live in `bin/genome_lookup.py` (and its test in `tests/test_results_index.py`), and no open lane edits them, so this PR changes them. `bin/locityper_report.py` is changed for the same reason: its test is in the list.
- Left for another PR, because the files belong to lanes in flight: an early `DONE` check for `--locityper_data` in `main.nf`, next to the other panel checks, and the four STR reading fixes in `bin/collect_summary.py` (the setup hint path, the no-5T wording, unmapped TA counts, and the stale state of traits and CPIC).

## How it was checked

- `tests/test_results_index.py`: a summary with HTT and FMR1 in the expansions, a UGT1A1 reading, and FCGR and HLA-DRB1 in Locityper. All five are terms of their sections, no region is a term, and the lookup finds HTT and FCGR. Red control in the test: without the `locus` key none of them is a term. On the pre-fix `GENE_KEYS` the test fails 2 checks ("HTT: None, FMR1: None ...", "none of the reports mentions HTT").
- `tests/test_locityper_report.py`, new cases 2, 8 and 9: missing read counts fail; the recorded command and `meta switches` keep both Locityper switches (control: a run without them records none); step 53 and the module's guard lines, run in bash, stop on a missing `DONE` and on another FASTA and pass on this one (control: a `DONE` naming only the `.fai` is refused). On the pre-fix code the test fails 8 checks.
- `scripts/ci/check-output-contract.py --static`: 0 problems; without the two rows it reports 2 (`switch LOCITYPER_DATA ... has no row`).
- CI on 6eb53e8c: Guard green, its unit-tests job ran every new check above (https://github.com/GeiserX/OwnGenome/actions/runs/42298606493). Lint, Nextflow, Personal data, Container Test and Docs green too. The first push failed the personal-data scan, which read the test's repeat counts as genotypes; the fixture no longer carries any.
- e2e: this PR adds no e2e case. The fork's e2e run on the final head is the gate.

Bead: pgp-5bja.28


<!-- This is an auto-generated comment: release notes by coderabbit.ai -->

<!-- coderabbit-ai-generated:v1 -->

## Summary by CodeRabbit

* **Bug Fixes**
  * Gene names listed under `locus` are now included in report matching and genome lookup.
  * Reports show `.` for unavailable or invalid read counts, and calls without valid counts fail the read-count filter.
  * Locityper stops before running if its database is missing the required completion marker or its FASTA reference checks do not match.

* **Documentation**
  * Clarified database requirements and how missing read counts affect filtering.
  * Documented the Locityper environment settings used by `run-all`.

<!-- end of auto-generated comment: release notes by coderabbit.ai -->

E2E evidence: https://github.com/geiresx/OwnGenome/actions/runs/42340604784
---
 assets/output_contract.tsv      |   2 +
 bin/genome_lookup.py            |   7 +-
 bin/locityper_report.py         |  34 +++++++--
 docs/53-locityper.md            |   3 +-
 modules/local/locityper/main.nf |   8 ++-
 scripts/53-locityper.sh         |   7 +-
 scripts/run-all.sh              |   7 +-
 tests/test_locityper_report.py  | 121 +++++++++++++++++++++++++++++++-
 tests/test_results_index.py     |  33 +++++++++
 9 files changed, 203 insertions(+), 19 deletions(-)

diff --git a/assets/output_contract.tsv b/assets/output_contract.tsv
index cca7450b..4eb1563f 100644
--- a/assets/output_contract.tsv
+++ b/assets/output_contract.tsv
@@ -419,6 +419,8 @@ switch	EXTRA_CALLERS	-	-	-	-	stable	-	-
 switch	BENCHMARK	-	-	-	-	stable	-	-
 switch	FAMILIES	-	-	-	-	experimental	-	-
 switch	OWNGENOME_LIVE_URL	-	-	-	-	experimental	-	-
+switch	LOCITYPER_DATA	-	-	-	-	experimental	-	-
+switch	LOCITYPER_BG_REGION	-	-	-	-	experimental	-	-
 switch	NXF_VER	-	-	-	-	stable	-	-
 switch	NXF_WORK	-	-	-	-	stable	-	-
 switch	NXF_HOME	-	-	-	-	stable	-	-
diff --git a/bin/genome_lookup.py b/bin/genome_lookup.py
index e8624085..9eee8529 100755
--- a/bin/genome_lookup.py
+++ b/bin/genome_lookup.py
@@ -54,8 +54,11 @@
     ("gnomADe_AF", "gnomAD exomes AF"),
     ("MAX_AF", "highest population AF"),
 ]
-# Keys of a summary.json entry that name a gene.
-GENE_KEYS = ("gene", "genes", "Gene", "SYMBOL", "symbol")
+# Keys of a summary.json entry that name a gene. "locus" is how the repeat
+# expansions, the CFTR/UGT1A1/SLC6A4 readings and the Locityper rows name
+# theirs (HTT, CFTR, HLA-DRB1); their regions sit under other keys (region,
+# chrom), so a locus value is a name.
+GENE_KEYS = ("gene", "genes", "Gene", "SYMBOL", "symbol", "locus")
 SYMBOLIC = {"<*>", "<NON_REF>", "*"}
 IMPACT_RANK = {"HIGH": 0, "MODERATE": 1, "LOW": 2, "MODIFIER": 3}
 
diff --git a/bin/locityper_report.py b/bin/locityper_report.py
index e9f3ff36..96e467f9 100755
--- a/bin/locityper_report.py
+++ b/bin/locityper_report.py
@@ -16,8 +16,10 @@
                          which of their two copies), "." when not called
   quality                Locityper's genotype quality, Phred: the chance that
                          another genotype (often a very similar one) is right
-  total_reads            reads Locityper used at the locus
+  total_reads            reads Locityper used at the locus ("." when its
+                         result gives no count)
   unexplained_reads      reads that fit some haplotype but not the two picked
+                         ("." when its result gives no count)
   weighted_distance      the distance from the call to the other likely
                          genotypes, weighted by their probability (edit
                          distance in bp with a database augmented by
@@ -27,7 +29,8 @@
                          for their HLA and KIR benchmark (Prodanov et al.
                          2025, Zenodo 15879519 README): no warning, a weighted
                          distance of 30 or less, and fewer than 1000 or fewer
-                         than 20% unexplained reads. "." when not called
+                         than 20% unexplained reads; fail when Locityper gave
+                         no read counts. "." when not called
   reason                 why a locus was not called, or why the filter failed
 
 A locus Locityper did not finish is "not called", never a genotype. Exit code
@@ -82,6 +85,10 @@ def fmt(x, digits):
     return f"{v:.{digits}f}"
 
 
+def dot(n):
+    return "." if n is None else str(n)
+
+
 def warnings_of(res):
     w = res.get("warnings")
     if w in (None, "", "*", []):
@@ -91,6 +98,20 @@ def warnings_of(res):
     return [str(x) for x in w]
 
 
+def count(res, key):
+    """A read count of res.json.gz as an int, None when it is absent or not a whole number."""
+    v = res.get(key)
+    if isinstance(v, bool):
+        return None
+    if isinstance(v, int):
+        return v if v >= 0 else None
+    if isinstance(v, float) and v.is_integer() and v >= 0:
+        return int(v)
+    if isinstance(v, str) and v.strip().isdigit():
+        return int(v)
+    return None
+
+
 def judge(res, warns):
     """pass or fail by the authors' filter, with the reason it failed."""
     why = []
@@ -105,9 +126,10 @@ def judge(res, warns):
         why.append("no weighted distance (the database has no haplotype alignments: setup.sh runs locityper augment)")
     elif wd > MAX_WEIGHT_DIST:
         why.append(f"weighted distance {wd:.2f} is above {MAX_WEIGHT_DIST:g}")
-    total = int(res.get("total_reads") or 0)
-    unexpl = int(res.get("unexpl_reads") or 0)
-    if not (unexpl < MAX_UNEXPLAINED or (total > 0 and unexpl < MAX_UNEXPLAINED_SHARE * total)):
+    total, unexpl = count(res, "total_reads"), count(res, "unexpl_reads")
+    if total is None or unexpl is None:
+        why.append("Locityper gave no read counts (total_reads, unexpl_reads), so the unexplained-read rule cannot be checked")
+    elif not (unexpl < MAX_UNEXPLAINED or (total > 0 and unexpl < MAX_UNEXPLAINED_SHARE * total)):
         why.append(f"{unexpl} of {total} reads unexplained (1000 or more, and 20% or more)")
     return ("fail" if why else "pass"), "; ".join(why)
 
@@ -137,7 +159,7 @@ def row_for(name, region, gt_dir, exit_code):
     return {
         "locus": name, "region": region, "status": "called", "genotype": ",".join(haps),
         "haplotype_1": haps[0], "haplotype_2": haps[1], "quality": fmt(res.get("quality"), 1),
-        "total_reads": str(int(res.get("total_reads") or 0)), "unexplained_reads": str(int(res.get("unexpl_reads") or 0)),
+        "total_reads": dot(count(res, "total_reads")), "unexplained_reads": dot(count(res, "unexpl_reads")),
         "weighted_distance": fmt(res.get("weight_dist"), 3), "warnings": ",".join(warns) or ".",
         "filter": filt, "reason": why or ".",
     }
diff --git a/docs/53-locityper.md b/docs/53-locityper.md
index 5972909c..5b4f415d 100644
--- a/docs/53-locityper.md
+++ b/docs/53-locityper.md
@@ -51,6 +51,7 @@ We left C4, AMY1 and LPA out. They are not in the authors' database, and on the
 ## Input
 
 - `${SAMPLE}/aligned/${SAMPLE}_sorted.bam` with its `.bai`, aligned to `REF_FASTA`, the reference the database was built from. Locityper reads the BAM's contigs by their order in the reference.
+- The database `setup.sh --locityper-data` built, with its `DONE` (written last). `DONE` names the sha256 of the `.fai` and of the FASTA it was built from. Step 53 and the Nextflow `LOCITYPER` stop before Locityper runs when `DONE` is missing or names another `.fai` or FASTA: a FASTA with the same contigs and other bases is another reference. Hashing a GRCh38 FASTA takes about 10 to 20 seconds.
 
 ## Command
 
@@ -66,7 +67,7 @@ With `run-all.sh`: `TOOLS=...,locityper`. It is opt-in, so a run lists step 53 o
 
 `bin/locityper_report.py` reads `gt/loci/<locus>/res.json.gz` and writes one row per locus. A locus Locityper did not finish reads "not called", with Locityper's exit code when it stopped.
 
-The filter is the one the authors used for their HLA and KIR benchmark (Zenodo 15879519, README): a call passes with no Locityper warning, a weighted distance of 30 or less, and fewer than 1000 or fewer than 20% unexplained reads. The weighted distance is the edit distance from the call to the other likely genotypes, weighted by their probability: a high value means the reads fit very different genotypes about as well.
+The filter is the one the authors used for their HLA and KIR benchmark (Zenodo 15879519, README): a call passes with no Locityper warning, a weighted distance of 30 or less, and fewer than 1000 or fewer than 20% unexplained reads. A call whose result gives no read counts fails it, and its `total_reads` and `unexplained_reads` read `.`. The weighted distance is the edit distance from the call to the other likely genotypes, weighted by their probability: a high value means the reads fit very different genotypes about as well.
 
 Locityper names haplotypes of panel people, such as `HG00733.1`, not alleles. The two callers are shown side by side and not compared: an allele name from T1K and a panel haplotype are not the same kind of answer.
 
diff --git a/modules/local/locityper/main.nf b/modules/local/locityper/main.nf
index 9a767298..39b1d384 100644
--- a/modules/local/locityper/main.nf
+++ b/modules/local/locityper/main.nf
@@ -39,8 +39,12 @@ process LOCITYPER {
     script:
     def bg = params.locityper_bg_region ? "-b ${params.locityper_bg_region}" : ''
     """
-    # The database's k-mer counts and haplotypes come from one reference: DONE names its .fai.
-    grep -q "fai=\$(sha256sum ${fai} | cut -d' ' -f1) " ${data}/DONE || {
+    # setup.sh --locityper-data writes DONE last: without it the install stopped.
+    [ -s ${data}/DONE ] || {
+        echo "ERROR: Locityper's database ${data} is not finished (no DONE); run scripts/setup.sh --locityper-data <genome_dir> again" >&2; exit 1; }
+    # The database's k-mer counts and haplotypes come from one reference: DONE names
+    # the digests of its .fai and its FASTA (a FASTA with the same contigs and other bases is another reference).
+    grep -q "fai=\$(sha256sum ${fai} | cut -d' ' -f1) fasta=\$(sha256sum ${fasta} | cut -d' ' -f1) " ${data}/DONE || {
         echo "ERROR: Locityper's database ${data} was built from another reference than ${fasta} (setup.sh --locityper-data)" >&2; exit 1; }
     mkdir -p work
     rc=0
diff --git a/scripts/53-locityper.sh b/scripts/53-locityper.sh
index f1d38b1c..dc32f91c 100755
--- a/scripts/53-locityper.sh
+++ b/scripts/53-locityper.sh
@@ -51,9 +51,10 @@ done
 for f in "${DATA}/DONE" "${DATA}/counts.jf" "${DATA}/loci.bed"; do
   [ -s "$f" ] || { echo "ERROR: Locityper's database not found: ${f} (run: ./scripts/setup.sh --locityper-data ${GENOME_DIR})" >&2; exit 1; }
 done
-# The database's k-mer counts and haplotypes come from one reference: DONE names its .fai.
-FAI_SHA=$(_digest sha256 "${REF_FASTA}.fai")
-grep -q "fai=${FAI_SHA} " "${DATA}/DONE" || {
+# The database's k-mer counts and haplotypes come from one reference: DONE names the digests
+# of its .fai and its FASTA (a FASTA with the same contigs and other bases is another reference).
+FAI_SHA=$(_digest sha256 "${REF_FASTA}.fai") FASTA_SHA=$(_digest sha256 "$REF_FASTA")
+grep -q "fai=${FAI_SHA} fasta=${FASTA_SHA} " "${DATA}/DONE" || {
   echo "ERROR: Locityper's database ${DATA} was built from another reference than ${REF_FASTA} (run: ./scripts/setup.sh --locityper-data ${GENOME_DIR})" >&2; exit 1; }
 echo "BAM: ${BAM}"
 echo "Database: ${DATA} ($(grep -vc '^#' "${DATA}/loci.bed") loci)"
diff --git a/scripts/run-all.sh b/scripts/run-all.sh
index ab6c3cdf..ec353a40 100755
--- a/scripts/run-all.sh
+++ b/scripts/run-all.sh
@@ -22,8 +22,9 @@
 #   REF_FASTA, EH_CATALOG (a file, or bundled), MANTA_CALL_REGIONS, PARASCOPY_POPULATION, PARAPHASE_REGIONS, ANCESTRY_PANEL (--ancestry_ref; none: raw
 #   scores), ANCESTRY_FINE_DIR (--ancestry_fine_panel, step 26b,
 #   always 100 bootstrap replicates; none: skip), ARCHAIC_DIR and ARCHAIC_POP (--archaic_data and --archaic_pop,
-#   step 43; ARCHAIC_DIR=none: skip) and PGSC_CALC_DIR (--pgsc_calc) as for
-#   the single steps; the panels and pgsc_calc are passed when installed.
+#   step 43; ARCHAIC_DIR=none: skip), PGSC_CALC_DIR (--pgsc_calc), LOCITYPER_DATA (--locityper_data, step 53's
+#   database) and LOCITYPER_BG_REGION (--locityper_bg_region) as for the single steps; the panels, the databases and
+#   pgsc_calc are passed when installed.
 # Run as scripts after the pipeline: GRIDSS=true (04b), IMPUTATION=true (14),
 #   SOMATIC=true (29), EXTRA_CALLERS=gatk,freebayes,strelka2,octopus (03a-03d), BENCHMARK=true (needs
 #   EXTRA_CALLERS or a second caller VCF); then the HTML report (24) and the text report, also after a failed pipeline.
@@ -51,7 +52,7 @@ shift 2; SCRIPT_DIR=$(cd "$(dirname "$0")" && pwd) USER_THREADS=${THREADS:-}
 SWITCHES=(GENOME_DIR SKIP_VALIDATION THREADS SKIP_TRIM INTERVALS ALIGN_DIR TOOLS KIR HLA_QC HPO REF_FASTA EH_CATALOG
   MANTA_CALL_REGIONS PARASCOPY_POPULATION PARAPHASE_REGIONS ANCESTRY_PANEL ANCESTRY_FINE_DIR ARCHAIC_DIR ARCHAIC_POP
   PGSC_CALC_DIR ANCIENT_PANEL_DIR VIRAL_REF MEGANE_KMERS MOCHA_DATA GRIDSS IMPUTATION SOMATIC EXTRA_CALLERS BENCHMARK
-  FAMILIES OWNGENOME_LIVE_URL NXF_VER NXF_WORK NXF_HOME NXF_OFFLINE)
+  FAMILIES OWNGENOME_LIVE_URL LOCITYPER_DATA LOCITYPER_BG_REGION NXF_VER NXF_WORK NXF_HOME NXF_OFFLINE)
 RERUN_CMD="cd $(printf '%q' "$PWD") && "; for v in "${SWITCHES[@]}"; do [ -z "${!v+x}" ] || RERUN_CMD+="${v}=$(printf '%q' "${!v}") "; done
 RERUN_CMD+=$(printf '%q ' "${SCRIPT_DIR}/run-all.sh" "$SAMPLE" "$SEX" "$@"); RERUN_CMD=${RERUN_CMD% }
 # With -bg nextflow returns at once, and the steps would be recorded and the reports written before the run ends.
diff --git a/tests/test_locityper_report.py b/tests/test_locityper_report.py
index 7ec725c1..f2bf5932 100755
--- a/tests/test_locityper_report.py
+++ b/tests/test_locityper_report.py
@@ -27,7 +27,17 @@
      installed lists it as running and passes --tools locityper,
      --locityper_data and --locityper_bg_region; without the database it is
      skipped, "data not installed". Control: the opted-in checks fail on the
-     default run's output.
+     default run's output. A result without read counts fails the filter
+     and its counts read '.', never 0 of 0 (case 2);
+  8. LOCITYPER_DATA and LOCITYPER_BG_REGION are run-all.sh switches: the
+     command logs/run_status.tsv keeps for owngenome_recheck.py --rerun names
+     both, and 'meta switches' lists them. Control: a run without
+     LOCITYPER_DATA records none;
+  9. step 53 and the LOCITYPER module's guard (its script lines before
+     Locityper, run in bash) stop on a missing DONE and on a DONE that names
+     this .fai with another FASTA, and pass on this .fai and FASTA. Control:
+     a DONE that names only the .fai, which the old guard accepted, is
+     refused.
 
 Run: python3 tests/test_locityper_report.py
 """
@@ -117,6 +127,12 @@ def row(result=None, success=True, raw=None, exit_code=0, truncate=False, damage
     r = row(res(wd=None))
     out.append(("2. no weighted distance fails, pointing at locityper augment",
                 r["filter"] == "fail" and "augment" in r["reason"] and r["weighted_distance"] == ".", r))
+    r = row({k: v for k, v in res().items() if k not in ("total_reads", "unexpl_reads")})
+    out.append(("2. no read counts fails, saying so, and the counts read '.' (not 0 of 0)",
+                r["filter"] == "fail" and "no read counts" in r["reason"] and r["total_reads"] == "."
+                and r["unexplained_reads"] == ".", r))
+    r = row(res(total=None, unexpl=12))
+    out.append(("2. a null total_reads fails too", r["filter"] == "fail" and "no read counts" in r["reason"], r))
     r = row(res(), success=False, exit_code=101)
     out.append(("3. no success file: not called, with Locityper's exit code",
                 r["status"] == "not called" and "exit code 101" in r["reason"] and r["genotype"] == "."
@@ -143,7 +159,7 @@ def row(result=None, success=True, raw=None, exit_code=0, truncate=False, damage
     return out
 
 
-def run_all(w, tools, done=True):
+def run_all(w, tools, done=True, extra=None):
     """scripts/run-all.sh on a placeholder BAM and VCF under the fake docker
     and nextflow; returns its output and the arguments it gave nextflow."""
     import shutil
@@ -164,6 +180,7 @@ def run_all(w, tools, done=True):
            "LOCITYPER_BG_REGION": "chr20:10000001-10500000"}
     if tools:
         env["TOOLS"] = tools
+    env.update(extra or {})
     p = subprocess.run(["bash", os.path.join(REPO, "scripts", "run-all.sh"), "s1", "male"], env=env, stdin=subprocess.DEVNULL,
                        capture_output=True, text=True)
     log = open(os.path.join(w, "docker.log")).read() if os.path.exists(os.path.join(w, "docker.log")) else ""
@@ -203,6 +220,103 @@ def plan_checks(w):
     failed = [d for d, ok, _ in opted_in(out_off, nf_off) if not ok]
     check("7. control: the opted-in checks fail on the default run's output", len(failed) == 2, failed)
 
+    # 8. LOCITYPER_DATA and LOCITYPER_BG_REGION are run-all.sh switches: the
+    # command logs/run_status.tsv keeps for owngenome_recheck.py --rerun names
+    # them, and 'meta switches' lists them, so a replay keeps the database.
+    def meta(w, key):
+        path = os.path.join(w, "g", "s1", "logs", "run_status.tsv")
+        for line in (open(path) if os.path.exists(path) else []):
+            c = line.rstrip("\n").split("\t")
+            if c[:2] == ["meta", key]:
+                return c[2]
+        return ""
+    alt = os.path.join(w, "elsewhere", "locityper")  # outside the folder run_all clears
+    os.makedirs(alt, exist_ok=True)
+    with open(os.path.join(alt, "DONE"), "w") as fh:
+        fh.write("ok\n")
+    out, nf = run_all(os.path.join(w, "alt"), "pharmcat,locityper", done=False, extra={"LOCITYPER_DATA": alt})
+    cmd, sw = meta(os.path.join(w, "alt"), "command"), meta(os.path.join(w, "alt"), "switches").split()
+    check("8. LOCITYPER_DATA reaches Nextflow", f"--locityper_data {alt}" in nf, nf)
+    check("8. the recorded command keeps LOCITYPER_DATA and LOCITYPER_BG_REGION",
+          f"LOCITYPER_DATA={alt} " in cmd and "LOCITYPER_BG_REGION=chr20:10000001-10500000 " in cmd, cmd)
+    check("8. 'meta switches' lists both", {"LOCITYPER_DATA", "LOCITYPER_BG_REGION"} <= set(sw), sw)
+    cmd_off = meta(os.path.join(w, "default"), "command")
+    check("8. control: a run without LOCITYPER_DATA records none, so the check above can fail",
+          bool(cmd_off) and "LOCITYPER_DATA=" not in cmd_off, cmd_off)
+
+
+def done_guards(w):
+    """9. scripts/53-locityper.sh and the LOCITYPER module stop before
+    Locityper when DONE is missing or names another .fai or FASTA: DONE names
+    both digests (setup.sh), and a FASTA with the same contigs and other bases
+    is another reference."""
+    import hashlib
+    import re
+    import shutil
+    g = os.path.join(w, "g")
+    ref = os.path.join(g, "reference", "GRCh38_no_alt_analysis_set.fasta")
+    data = os.path.join(g, "reference", "locityper-" + lr_panel_version())
+    for f in ("s1/aligned/s1_sorted.bam", "s1/aligned/s1_sorted.bam.bai"):
+        os.makedirs(os.path.dirname(os.path.join(g, f)), exist_ok=True)
+        open(os.path.join(g, f), "w").write("placeholder\n")
+    os.makedirs(data)
+    open(ref, "w").write(">chr1\nACGT\n")
+    open(ref + ".fai", "w").write("chr1\t4\t6\t4\t5\n")
+    open(os.path.join(data, "counts.jf"), "w").write("x\n")
+    open(os.path.join(data, "loci.bed"), "w").write("chr1\t0\t4\tX\n")
+
+    def sha(path):
+        return hashlib.sha256(open(path, "rb").read()).hexdigest()
+    good = f"fai={sha(ref + '.fai')} fasta={sha(ref)} image=x panel=y loci=z\n"
+    other = f"fai={sha(ref + '.fai')} fasta={hashlib.sha256(b'>chr1 TGCA').hexdigest()} image=x panel=y loci=z\n"
+    env = {"PATH": os.path.join(REPO, "scripts", "ci", "fake-docker") + os.pathsep + os.environ["PATH"], "HOME": w,
+           "GENOME_DIR": g, "FAKE_DOCKER_LOG": os.path.join(w, "docker.log")}
+
+    def step53(done):
+        if done is None:
+            if os.path.exists(os.path.join(data, "DONE")):
+                os.remove(os.path.join(data, "DONE"))
+        else:
+            open(os.path.join(data, "DONE"), "w").write(done)
+        p = subprocess.run(["bash", os.path.join(REPO, "scripts", "53-locityper.sh"), "s1"], env=env,
+                           stdin=subprocess.DEVNULL, capture_output=True, text=True)
+        return p.stdout + p.stderr
+
+    out = step53(other)
+    check("9. step 53: a DONE naming the same .fai and another FASTA stops it",
+          "built from another reference" in out and "Database:" not in out, out[-800:])
+    out = step53(good)
+    check("9. step 53: a DONE naming this .fai and FASTA passes the guard",
+          "another reference" not in out and "Database: " in out, out[-800:])
+    out = step53(None)
+    check("9. step 53: no DONE stops it", "database not found" in out and "Database:" not in out, out[-800:])
+
+    # The module's guard, as the task runs it: the lines of its script before Locityper, with
+    # Nextflow's ${...} filled in and \$ unescaped.
+    src = open(os.path.join(REPO, "modules", "local", "locityper", "main.nf")).read()
+    m = re.search(r'script:.*?"""\n(.*?)\n\s*mkdir -p work', src, re.S)
+    check("9. module: the guard is in the script block", m is not None)
+    if not m:
+        return
+    guard = (m.group(1).replace("${data}", data).replace("${fai}", ref + ".fai").replace("${fasta}", ref)
+             .replace("\\$", "$"))
+    if not shutil.which("sha256sum"):
+        print("[SKIP] 9. module guard: no sha256sum on this machine")
+        return
+
+    def module(done):
+        if done is None:
+            os.remove(os.path.join(data, "DONE"))
+        else:
+            open(os.path.join(data, "DONE"), "w").write(done)
+        p = subprocess.run(["bash", "-c", guard + "\necho GUARD-PASSED"], capture_output=True, text=True)
+        return p.stdout + p.stderr
+    check("9. module: another FASTA stops it", "built from another reference" in module(other))
+    check("9. module: this .fai and FASTA pass", "GUARD-PASSED" in module(good), module(good))
+    check("9. module: no DONE stops it with that reason", "not finished (no DONE)" in module(None))
+    check("9. control: a DONE naming only the .fai (the old guard's test) is not enough",
+          "GUARD-PASSED" not in module(f"fai={sha(ref + '.fai')} image=x\n"))
+
 
 def main():
     with tempfile.TemporaryDirectory() as t:
@@ -297,6 +411,9 @@ def main():
         # 7. run-all.sh plans step 53 only when TOOLS names it
         plan_checks(os.path.join(t, "plan"))
 
+        # 9. the DONE guards of step 53 and the module
+        done_guards(os.path.join(t, "guard"))
+
     print("\nRESULT:", "ALL PASS" if FAILS == 0 else f"{FAILS} FAILED")
     return 1 if FAILS else 0
 
diff --git a/tests/test_results_index.py b/tests/test_results_index.py
index 6b9b6247..8bcc6b3a 100755
--- a/tests/test_results_index.py
+++ b/tests/test_results_index.py
@@ -167,6 +167,39 @@ def main():
         check("a prose gene field: the lookup finds GC too",
               "none of the reports mentions GC" not in lookup_says(ppath, "GC"), lookup_says(ppath, "GC"))
 
+        # Sections that name their gene under "locus": the repeat expansions
+        # (flagged, key loci, the CFTR and UGT1A1 readings) and Locityper.
+        loc = copy.deepcopy(s_tsv)
+        loc["sections"]["expansions"].update(state="ok", data={
+            "flagged": [{"locus": "HTT", "status": "pathogenic"}],
+            "key_loci": [{"locus": "FMR1", "repeat_count": "."}],
+            "readings": [{"locus": "UGT1A1", "repeat": "promoter TA repeat"}]})
+        loc["sections"]["locityper"].update(state="ok", data={
+            "loci": [{"locus": "FCGR", "region": "chr1:161505430-161678654", "status": "called"},
+                     {"locus": "HLA-DRB1", "region": "chr6:32484348-32593538", "status": "called"}]})
+        lpath = os.path.join(work, "summary-locus.json")
+        with open(lpath, "w") as f:
+            json.dump(loc, f)
+        lt = results_index.build(loc, tsv, render_report.SECTION_OF)["terms"]
+        want = {"HTT": "expansions", "FMR1": "expansions", "UGT1A1": "expansions", "FCGR": "locityper",
+                "HLA-DRB1": "locityper"}
+        check("a 'locus' key: HTT, FMR1 and UGT1A1 are terms of the expansions section, FCGR and HLA-DRB1 of Locityper",
+              all(sec in [it["section"] for it in lt.get(t, [])] for t, sec in want.items()),
+              {t: lt.get(t) for t in want})
+        check("a 'locus' key: a region is not a term",
+              not [t for t in lt if t.startswith("CHR")], sorted(lt))
+        check("a 'locus' key: the lookup finds HTT and FCGR too",
+              all("none of the reports mentions" not in lookup_says(lpath, t) for t in ("HTT", "FCGR")),
+              [lookup_says(lpath, t) for t in ("HTT", "FCGR")])
+        keep = genome_lookup.GENE_KEYS
+        genome_lookup.GENE_KEYS = tuple(k for k in keep if k != "locus")
+        try:
+            gone = results_index.build(loc, tsv, render_report.SECTION_OF)["terms"]
+        finally:
+            genome_lookup.GENE_KEYS = keep
+        check("red control: without the 'locus' key none of them is a term",
+              not set(want) & set(gone), sorted(set(want) & set(gone)))
+
         # 2. bash and Nextflow forms
         nf = results_index.build(s_vcf, vcf, render_report.SECTION_OF)
         for key in results_index.TABLES: