fix(report): read every repeat, Y and phase result truthfully, and bring the docs and pins in line with 1.0.0

bin/collect_summary.py

@@ -2100,6 +2100,10 @@ def str_not_graded(d, s):
 
 # The repeats STRchive leaves out (assets/str_extra_loci.json), read from
 # their counts when step 09 sized them (STR_EXTRA_LOCI=yes; docs/09-str-expansions.md).
+STR_EXTRA_SETUP = "STR_EXTRA_LOCI=yes ./scripts/setup.sh --strchive-data <genome_dir>"
+# What a section state other than ok says about a table another step wrote.
+OTHER_STEP_STATE = {"stale": "from an earlier run", "failed": "its step failed in the last run",
+                    "unreadable": "unreadable"}
 # CFTR intron 9 (legacy IVS8): ExpansionHunter writes the TG and the poly-T
 # tract as two records, so which TG count sits on the copy of which T count
 # is known only when one of the two is the same on both copies. R117H
@@ -2140,13 +2144,25 @@ def trait_genotype(d, s, rsid):
     return (row.get("status") or "", row.get("genotype") or "") if row else None
 
 
-def cftr_reading(by, d, s):
+def r117h_genotype(d, s, state):
+    """(genotype or None, why not read, note) of R117H in step 42's traits
+    table, given the traits section's state: a table of a failed or
+    unreadable step is not read, one from an earlier run is read and said so."""
+    if state in ("failed", "unreadable"):
+        return None, f"step 42's traits table is {OTHER_STEP_STATE[state]}", ""
+    r = trait_genotype(d, s, CFTR_R117H)
+    if not r or r[0] != "called":
+        return None, "no step 42 traits row for rs78655421 was called", ""
+    return r[1], "", (" (R117H from step 42's traits table of an earlier run)" if state == "stale" else "")
+
+
+def cftr_reading(by, d, s, traits_state="ok"):
     tg, t = by.get("CFTR_TG"), by.get("CFTR_T")
     out = {"locus": "CFTR", "repeat": "intron 9 (legacy IVS8) TG and poly-T tract", "source": CFTR_SOURCE}
     if not tg or not t:
         return dict(out, status="not_sized", genotype=".", phase=".",
-                    reading="not in step 09's output: STR_EXTRA_LOCI=yes setup.sh --strchive-data adds the tract "
-                            "to the catalog (docs/09-str-expansions.md).")
+                    reading=f"not in step 09's output: {STR_EXTRA_SETUP} adds the tract to the catalog "
+                            "(docs/09-str-expansions.md).")
     a, b = str_counts(tg), str_counts(t)
     if not a or not b:
         return dict(out, status="not_sized", genotype=".", phase=".",
@@ -2162,10 +2178,24 @@ def cftr_reading(by, d, s):
     other = sorted({n for n in b if n not in (5, 7, 9)})
     odd = (f" A poly-T length of {', '.join(f'{n}T' for n in other)} is none of the usual 5T, 7T and 9T."
            if other else "")
+    r117, why, note = r117h_genotype(d, s, traits_state)
     if not five:
-        return dict(out, five_t=False, phase="not_needed",
-                    reading=f"No 5T copy ({b[0]}T/{b[1]}T), so the tract does not change how a CFTR finding is "
-                            f"read.{odd}")
+        # No 5T risk of its own, but the poly-T still sets how R117H reads:
+        # with no 5T copy, an R117H copy sits on a 7T or 9T (or other) copy.
+        line = (f"No 5T copy ({b[0]}T/{b[1]}T): the tract adds no 5T allele of its own. The poly-T length still "
+                "matters beside R117H, which is milder and of variable penetrance on a 7T or 9T copy and CF-causing "
+                "with 5T in cis.")
+        if r117 is None:
+            line += f" R117H was not read ({why})."
+        elif r117.split("/").count("A") == 2:
+            copies = f"both {b[0]}T copies" if b[0] == b[1] else f"both copies ({b[0]}T and {b[1]}T)"
+            line += f" R117H ({r117}) is read: it sits on {copies}.{note}"
+        elif r117.split("/").count("A"):
+            copies = f"{b[0]}T" if b[0] == b[1] else f"{b[0]}T or {b[1]}T"
+            line += f" R117H ({r117}) is read: with no 5T copy it sits on a {copies} copy.{note}"
+        else:
+            line += f" No R117H copy.{note}"
+        return dict(out, five_t=False, phase="not_needed", r117h=r117 or ".", reading=line + odd)
     if len(pairs) == 1:
         # One TG count on both copies, or 5T on both copies with their two TG counts.
         tgs = [a[0]] if a[0] == a[1] else a
@@ -2175,36 +2205,34 @@ def cftr_reading(by, d, s):
     line = (f"5T: {five} copy{'s' if five > 1 else ''} with the 5T poly-T tract, {tg_text}. 5T lowers the splicing "
             "of exon 10 and is a variable-penetrance allele: with TG12 or TG13 on the same copy it more often "
             "causes CBAVD or another CFTR-related disorder when the other copy carries a CF-causing variant.")
-    r = trait_genotype(d, s, CFTR_R117H)
-    if not r or r[0] != "called":
+    if r117 is None:
         phase = "not_read"
-        r117 = (" R117H was not read (no step 42 traits row for rs78655421 was called), so whether it sits on the "
-                "5T copy is not known.")
+        text = f" R117H was not read ({why}), so whether it sits on the 5T copy is not known."
     else:
-        n = r[1].split("/").count("A")
+        n = r117.split("/").count("A")
         if n == 0:
             phase = "not_needed"
-            r117 = " No R117H copy: the 5T is read on its own."
+            text = " No R117H copy: the 5T is read on its own."
         elif n == 2 or five == 2:
             phase = "in_cis"
-            r117 = (f" R117H ({r[1]}) and 5T share a copy: R117H with 5T in cis is a CF-causing allele of "
+            text = (f" R117H ({r117}) and 5T share a copy: R117H with 5T in cis is a CF-causing allele of "
                     "variable penetrance.")
         else:
             phase = "unknown"
-            r117 = (" One R117H copy and one 5T copy: phase unknown. R117H lies 17.6 kb from the tract, farther "
+            text = (" One R117H copy and one 5T copy: phase unknown. R117H lies 17.6 kb from the tract, farther "
                     "than a read pair spans, so whether they share a copy (in cis, a CF-causing allele of variable "
                     "penetrance) or not (in trans) is not read; testing a parent tells them apart.")
-    return dict(out, five_t=True, phase=phase, r117h=r[1] if r and r[0] == "called" else ".",
-                reading=line + r117 + odd)
+        text += note
+    return dict(out, five_t=True, phase=phase, r117h=r117 or ".", reading=line + text + odd)
 
 
-def ugt1a1_reading(by, d, s):
+def ugt1a1_reading(by, d, s, cpic_state="ok"):
     item = by.get("UGT1A1")
     out = {"locus": "UGT1A1", "repeat": "promoter TA repeat", "source": UGT1A1_SOURCE}
     if not item:
         return dict(out, status="not_sized", genotype=".", pharmcat=".", agreement="not_compared",
-                    reading="not in step 09's output: STR_EXTRA_LOCI=yes setup.sh --strchive-data adds the repeat "
-                            "to the catalog (docs/09-str-expansions.md). PharmCAT (step 7) still reads UGT1A1 *28.")
+                    reading=f"not in step 09's output: {STR_EXTRA_SETUP} adds the repeat to the catalog "
+                            "(docs/09-str-expansions.md). PharmCAT (step 7) still reads UGT1A1 *28.")
     c = str_counts(item)
     if not c:
         return dict(out, status="not_sized", genotype=".", pharmcat=".", agreement="not_compared",
@@ -2214,15 +2242,20 @@ def ugt1a1_reading(by, d, s):
     genotype = "/".join(names)
     dec = sum(n in UGT1A1_DECREASED for n in names)
     ta = "/".join(f"TA{n - 1}" for n in c)
-    reading = (f"{genotype} ({ta}, {c[0]}/{c[1]} TA units as ExpansionHunter counts them): "
-               + ("no decreased-function TA allele." if not dec else
-                  f"{dec} decreased-function allele{'s' if dec > 1 else ''} (*28 and *37 are decreased function).")
-               + ("" if all(n in UGT1A1_ALLELES for n in c) else " A TA length outside TA5 to TA8 is not a named allele here."))
     out.update(status="read", genotype=genotype, ta_units=c)
+    reading = f"{genotype} ({ta}, {c[0]}/{c[1]} TA units as ExpansionHunter counts them): "
+    if not all(n in UGT1A1_ALLELES for n in c):
+        # A length with no named allele has no function to read, and so no
+        # count of decreased-function copies to set beside PharmCAT's.
+        return dict(out, pharmcat=".", agreement="not_compared",
+                    reading=reading + "a TA length outside TA5 to TA8 is not a named allele here, so its function "
+                                      "is not read and it is not compared with PharmCAT (step 7).")
+    reading += ("no decreased-function TA allele." if not dec else
+                f"{dec} decreased-function allele{'s' if dec > 1 else ''} (*28 and *37 are decreased function).")
     # PharmCAT's UGT1A1 call (step 27's table, with the reference-block fill
     # of bin/pgx_refblock_fill.py): its count of *28 and *37 copies.
     pc, agreement = ".", "not_compared"
-    p = first_existing(d, [f"cpic/{s}_phenotypes.tsv"])
+    p = first_existing(d, [f"cpic/{s}_phenotypes.tsv"]) if cpic_state not in ("failed", "unreadable") else None
     row = None
     if p:
         try:
@@ -2236,20 +2269,26 @@ def ugt1a1_reading(by, d, s):
             # A combination allele reads "*80+*28" or "[*80 + *28]".
             n = sum(any(x in re.split(r"[+ \[\]]", a) for x in UGT1A1_DECREASED) for a in alleles)
             agreement = "agrees" if n == dec else "differs"
+    old = " (step 27's table is from an earlier run)" if cpic_state == "stale" and agreement != "not_compared" else ""
     if agreement == "agrees":
-        reading += f" PharmCAT (step 7) called {pc}: the same number of *28 and *37 copies."
+        reading += f" PharmCAT (step 7) called {pc}: the same number of *28 and *37 copies{old}."
     elif agreement == "differs":
-        reading += (f" PharmCAT (step 7) called {pc}, a different number of *28 and *37 copies: read the reads at "
-                    "chr2:233760234 before using either.")
+        reading += (f" PharmCAT (step 7) called {pc}, a different number of *28 and *37 copies{old}: read the reads "
+                    "at chr2:233760234 before using either.")
+    elif cpic_state in ("failed", "unreadable"):
+        reading += f" Not compared with PharmCAT (step 27's table is {OTHER_STEP_STATE[cpic_state]})."
     else:
         reading += " Not compared with PharmCAT (no unambiguous UGT1A1 call from step 27)."
     return dict(out, pharmcat=pc, agreement=agreement, reading=reading)
 
 
-def str_readings(by, d, s):
+def str_readings(by, d, s, states=None):
     """What the counts of the extra repeats mean, one entry per locus, and the
-    5-HTTLPR line that no step can test."""
-    return [cftr_reading(by, d, s), ugt1a1_reading(by, d, s),
+    5-HTTLPR line that no step can test. states: the traits and cpic section
+    states, once collect() knows them (a stale, failed or unreadable table is
+    said so, not read as current)."""
+    states = states or {}
+    return [cftr_reading(by, d, s, states.get("traits", "ok")), ugt1a1_reading(by, d, s, states.get("cpic", "ok")),
             {"locus": "SLC6A4", "repeat": "5-HTTLPR promoter length polymorphism", "status": "not_tested",
              "genotype": ".", "reading": HTTLPR_READING, "source": "Bousman et al. 2023, PMID 37032427"}]
 
@@ -2639,9 +2678,14 @@ def sec_y_tree(d, s):
            "explore": haplogroup_explore("y", hg)}
     lvl = first_existing(d, [f"y_haplogroup/{s}_y_tree_path.tsv"])
     if lvl:
-        prows = read_tsv(lvl)
-        if not prows or "Branch" not in prows[0]:
-            raise Unreadable(f"{os.path.basename(lvl)} has no Branch column or no rows")
+        # bin/y_tree_place.py writes the header alone when no branch is derived:
+        # an empty path, not an unreadable file.
+        with open_text(lvl) as f:
+            rd = csv.DictReader(f, delimiter="\t")
+            prows = list(rd)
+            header = rd.fieldnames or []
+        if "Branch" not in header:
+            raise Unreadable(f"{os.path.basename(lvl)} has no Branch column")
         out["path"] = [{k: (row.get(k) or "") for k in
                         ("Level", "Branch", "SNPs", "On_tables", "Derived", "Ancestral", "No_call",
                          "Ancestral_SNPs")} for row in prows]
@@ -3966,10 +4010,15 @@ def over(chrom, a, b):
 
 def phase_haps(gt):
     """The copies (haplotype indexes) a phased GT puts the alt allele on, or
-    None when the GT is not phased: '0|1' -> {1}, '1|0' -> {0}, '0/1' -> None."""
+    None when the GT is not phased or a copy is not called: '0|1' -> {1},
+    '1|0' -> {0}, '0/1' -> None, '.|1' -> None (the missing copy could carry
+    the other allele, so no cis or trans can be read from it)."""
     if not gt or "|" not in gt:
         return None
-    return frozenset(i for i, a in enumerate(gt.split("|")) if a not in ("0", "."))
+    alleles = gt.split("|")
+    if "." in alleles:
+        return None
+    return frozenset(i for i, a in enumerate(alleles) if a != "0")
 
 
 def pair_phase(dl, v):
@@ -4855,6 +4904,13 @@ def collect(sample, sample_dir, declared_sex=None):
         sections[key] = sec
 
     cyp2d6 = cyp2d6_block(sample, sample_dir, sections)
+    # The repeat readings lean on the traits (R117H) and CPIC (UGT1A1)
+    # tables, whose sections come after the expansions one: read them again
+    # with those sections' states.
+    ex = sections["expansions"]
+    if ex["state"] in ("ok", "stale"):
+        ex["data"]["readings"] = str_readings({l["locus"]: l for l in ex["data"]["loci"]}, sample_dir, sample,
+                                              {k: sections[k]["state"] for k in ("traits", "cpic")})
     if sections["locityper"]["state"] in ("ok", "stale"):
         primary = locityper_primary(sections)
         for x in sections["locityper"]["data"]["loci"]:

docs/09-str-expansions.md

@@ -49,9 +49,9 @@ Each entry has a `NoThresholds` reason and no Stranger thresholds: step 9b gives
 
 The reports read the counts (`bin/collect_summary.py`, under the repeat expansions section, `readings` in `summary.json`):
 
-- **CFTR**: the TG and poly-T counts of both copies, such as `TG11-7T/TG11-7T`. ExpansionHunter writes the two tracts as two records, so when both the TG and the poly-T counts differ between the copies, which TG sits with which poly-T is not known: the line says `TG11/TG12 5T/7T` and names both pairings. A 5T copy gets the 5T line: 5T is a variable-penetrance allele, more often disease-causing with TG12 or TG13 on the same copy, and it matters when the other copy carries a CF-causing variant. The line then reads R117H (rs78655421), a row of step 42's traits table: no R117H copy; R117H and 5T in cis when one of them is on both copies; otherwise phase unknown, because R117H lies 17.6 kb from the tract, farther than a read pair spans. No 5T copy gives no 5T line.
-- **UGT1A1**: the TA count of each copy as its allele (`*36`, `*1`, `*28`, `*37`), and the number of decreased-function copies (`*28`, `*37`) compared with PharmCAT's UGT1A1 call of step 27, which reads `*28` from the small-variant calls with the reference-block fill of `bin/pgx_refblock_fill.py`. The line says whether the two agree; when they differ, read the reads before using either. PharmCAT's call is still the one step 27's drug guidance uses.
-- Without the opt-in, both lines say the repeat was not sized and how to add it.
+- **CFTR**: the TG and poly-T counts of both copies, such as `TG11-7T/TG11-7T`. ExpansionHunter writes the two tracts as two records, so when both the TG and the poly-T counts differ between the copies, which TG sits with which poly-T is not known: the line says `TG11/TG12 5T/7T` and names both pairings. A 5T copy gets the 5T line: 5T is a variable-penetrance allele, more often disease-causing with TG12 or TG13 on the same copy, and it matters when the other copy carries a CF-causing variant. The line then reads R117H (rs78655421), a row of step 42's traits table: no R117H copy; R117H and 5T in cis when one of them is on both copies; otherwise phase unknown, because R117H lies 17.6 kb from the tract, farther than a read pair spans. No 5T copy gives no 5T line, but the poly-T still matters beside R117H (milder on a 7T or 9T copy, CF-causing with 5T in cis), so the line still reads R117H and says which poly-T copy it sits on. A traits table from an earlier run is read and said to be so; one of a failed or unreadable step 42 is not read.
+- **UGT1A1**: the TA count of each copy as its allele (`*36`, `*1`, `*28`, `*37`), and the number of decreased-function copies (`*28`, `*37`) compared with PharmCAT's UGT1A1 call of step 27, which reads `*28` from the small-variant calls with the reference-block fill of `bin/pgx_refblock_fill.py`. The line says whether the two agree; when they differ, read the reads before using either. A TA length outside TA5 to TA8 has no named allele, so the line makes no function claim and does not compare it. A step 27 table from an earlier run is compared and said to be so; one of a failed or unreadable step is not compared. PharmCAT's call is still the one step 27's drug guidance uses.
+- Without the opt-in, both lines say the repeat was not sized and give the command that adds it: `STR_EXTRA_LOCI=yes ./scripts/setup.sh --strchive-data <genome_dir>`.
 
 `tests/test_str_catalog.py` checks each region against the GRCh38 sequence, `tests/test_collect_summary.py` checks the readings on made-up counts (a TG11-7T/TG11-7T sample is the red control of the 5T line), the e2e case `genotype-rule-rows-1-str-extra` sizes both tracts on made-up reads with ExpansionHunter 5.0.0, and the e2e case `str-extra-repeats-1-readings` runs step 09 and both reports on made-up reads and checks the lines.