fix(cn): show the step 35 interval, CYP21A2 and gene-dosage results in both reports and the couple screen

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	-	-

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,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 rows and not {"locus", "name", "status", "copies", "expected", "change"} <= set(rows[0]):
+    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
@@ -904,6 +917,27 @@ 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")
+    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):
     """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

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,35 @@ def cnv_note(row, dels):
     return ""
 
 
+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; 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 {}
+            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 +682,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 +705,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 +790,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 +852,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, 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)" % (
             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)

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}
 

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 ".",
     }