feat(report): print the LPA null allele beside the KIV-2 count and the SMN1 silent-carrier marker beside the SMN copy number

bin/collect_summary.py

@@ -763,6 +763,35 @@ def parascopy_cn_at(rows, chrom, pos):
     return {"cn": None, "text": "no row covers exon 7", "agcn": ".", "agcn_qual": ".", "agcn_filter": "."}
 
 
+MARKER_FIELDS = ("marker", "rsid", "status", "alt_reads", "depth", "fraction", "copies", "message")
+
+
+def read_marker_table(p, step, key):
+    """A copy-aware marker table (metric, value; bin/traits_report.py MARKERS):
+    step 46's KIV-2 4925G>A beside the count or step 35's SMN1 g.27134T>G
+    beside the psCN. A table that is not there or has no known status reads
+    as not assessed, with the reason in its message, so the report line
+    never goes missing beside the count."""
+    m = traits_report.MARKERS[key]
+    what = m["name"]
+    out = dict.fromkeys(MARKER_FIELDS, ".")
+    out.update(marker=what, rsid=m["dbsnp"][0], status="not assessed")
+    if not os.path.isfile(p):
+        out["message"] = (f"{what}: not assessed, step {step} wrote no {os.path.basename(p)} (a run from before "
+                          "this reader, or the check did not run).")
+        return out
+    try:
+        kv = dict(read_tsv_pairs(p))
+    except OSError as e:
+        kv = {"message": f"{os.path.basename(p)} could not be read: {e}"}
+    if kv.get("status") not in traits_report.MARKER_STATUSES:
+        out["message"] = (f"{what}: not assessed, {os.path.basename(p)} has no status row "
+                          f"({', '.join(traits_report.MARKER_STATUSES)}).")
+        return out
+    out.update({k: kv[k] for k in MARKER_FIELDS if kv.get(k)})
+    return out
+
+
 def smn_reading(smn1):
     """What an SMN1 copy number means for SMA, in one line."""
     if smn1 is None:
@@ -881,6 +910,9 @@ def sec_parascopy(d, s):
     catalog = {r["table"]: r for r in asset_rows("paralog_loci.tsv")} if others else {}
     return p, {"smn1": smn1, "smn2": {"cn": smn2["cn"], "text": smn2["text"]},
                "agcn": smn1["agcn"], "agcn_qual": smn1["agcn_qual"], "agcn_filter": smn1["agcn_filter"],
+               # SMN1 g.27134T>G beside the psCN (step 35, <id>_smn1_marker.tsv)
+               "smn1_marker": read_marker_table(os.path.join(d, "paralogs", f"{s}_smn1_marker.tsv"), "35",
+                                                "SMN1_27134"),
                "rows": len(rows),
                # HBA1/HBA2 from depth (bin/hba_copy_number.py); None when
                # step 35 wrote no HBA table (a run from before the check)
@@ -2786,6 +2818,9 @@ def sec_lpa(d, s):
         raise Unreadable(f"{os.path.basename(p)} has no status row (estimated, too little depth or not assessed)")
     data = {k: kv.get(k) or "." for k in ("status", "kiv2_copies", "ci95_low", "ci95_high", "class",
                                           "kiv2_depth", "flank_depth", "mapq0_share", "message", "advice")}
+    # KIV-2 4925G>A beside the count (step 46, lpa/<id>_lpa_kiv2_4925.tsv)
+    data["null_allele"] = read_marker_table(os.path.join(os.path.dirname(p), f"{s}_lpa_kiv2_4925.tsv"), "46",
+                                            "LPA_KIV2_4925")
     variants = []
     t = first_existing(d, [f"traits/{s}_traits.tsv"])
     for r in (read_tsv(t) if t else []):

bin/live_results.py

@@ -301,6 +301,8 @@ def a_lpa(S, s):
     d = S["lpa"]["data"]
     facts = [["KIV-2 repeats", f"{d.get('kiv2_copies')} (95% interval {d.get('ci95_low')} to {d.get('ci95_high')})"
               if d.get("status") == "estimated" else d.get("status")], ["Class", d.get("class")]]
+    if d.get("null_allele"):
+        facts.append(["KIV-2 4925G>A null allele", rr.marker_state(d["null_allele"])])
     return {"facts": facts, "lines": rr.lpa_lines(d), "attention": d.get("flag") == "high"}
 
 
@@ -357,15 +359,19 @@ def a_couple(S, s):
 
 
 def a_paralogs(S, s):
-    lines = []
+    lines, facts = [], []
     if (S.get("parascopy") or {}).get("state") in ("ok", "stale"):
         d = S["parascopy"]["data"]
         lines += rr.parascopy_lines(d) + rr.parascopy_locus_lines(d) + rr.cn_interval_lines(d.get("cn_intervals"))
+        if d.get("smn1_marker"):
+            # the silent-carrier marker beside the SMN1 psCN
+            facts = [["SMN1 copies", (d.get("smn1") or {}).get("text", ".")],
+                     ["SMN1 g.27134T>G silent-carrier marker", rr.marker_state(d["smn1_marker"])]]
     rows = []
     if (S.get("paralogs") or {}).get("state") in ("ok", "stale"):
         rows = [[g.get("gene"), g.get("copy_number"), g.get("status"), g.get("detail")]
                 for g in S["paralogs"]["data"].get("genes") or []]
-    return {"lines": lines, "table": table(["Gene", "Copies", "Status", "Detail"], rows)}
+    return {"facts": facts, "lines": lines, "table": table(["Gene", "Copies", "Status", "Detail"], rows)}
 
 
 def a_pharmcat(S, s):

bin/render_report.py

@@ -119,6 +119,23 @@ def ancient_lines(d):
                  "health on its own.")
 LPA_NOTE = ("The KIV-2 repeat count from read depth is a research estimate, not an Lp(a) level: the level is a "
             "blood test, and the count explains only part of the difference between people.")
+LPA_NULL_CAVEAT = ("A null allele lowers the Lp(a) made from the copies it sits on; the Lp(a) blood test is the "
+                   "measure.")
+SMN1_MARKER_CAVEAT = ("The marker raises the chance that two SMN1 copies sit on one chromosome (2+0, a silent carrier); "
+                      "it does not prove 2+0, and a carrier test in a clinical laboratory confirms it.")
+
+
+def marker_state(m):
+    """A copy-aware marker's state in a few words (collect_summary.read_marker_table):
+    present or uncertain with the reads and the copies that carry it, not seen
+    with the reads, too little depth, or not assessed."""
+    st = m.get("status") or "not assessed"
+    if st == "not assessed":
+        return st
+    out = f"{st}, {m.get('alt_reads', '.')} of {m.get('depth', '.')} reads carry it"
+    if m.get("copies") not in (None, "", "."):
+        out += f", about {m['copies']} copies"
+    return out
 
 
 def combined_cell(h):
@@ -1799,7 +1816,11 @@ def paralog_lines(d):
 def parascopy_lines(d):
     """Step 35's lines, the same in both reports."""
     smn1, smn2 = d.get("smn1") or {}, d.get("smn2") or {}
-    return [f"SMN1 copies: {smn1.get('text', '.')}",
+    mk = d.get("smn1_marker")
+    marker = ([f"SMN1 g.27134T>G silent-carrier marker ({mk.get('rsid') or '.'}): {marker_state(mk)}",
+               f"SMN1 marker reading: {mk.get('message', '.')}",
+               f"SMN1 marker caveat: {SMN1_MARKER_CAVEAT}"] if mk else [])
+    return [f"SMN1 copies: {smn1.get('text', '.')}"] + marker + [
             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', '.')})",
@@ -2146,6 +2167,10 @@ def lpa_lines(d):
     else:
         out = [f"KIV-2 repeats: {d.get('status', '.')}"]
     out.append(f"  {d.get('message', '.')}")
+    na = d.get("null_allele")
+    if na:
+        out += [f"KIV-2 4925G>A null allele ({(na.get('rsid') or '.').split(' ')[0]}): {marker_state(na)}",
+                f"  {na.get('message', '.')}", f"  {LPA_NULL_CAVEAT}"]
     for v in d.get("variants", []):
         n = v["risk_alleles"]
         out.append(f"{v['rsid']}: {v['genotype']} ({'not read' if n is None else str(n) + ' risk allele' + ('s' if n != 1 else '')})")
@@ -2897,7 +2922,8 @@ def a(key, lines):
     if S.get("parascopy", MISSING)["state"] != "missing":
         a("parascopy", card(S["parascopy"], "SMN1/SMN2 Copy Number (Parascopy)",
                [stat("Status", done_badge(S["parascopy"]))]
-               + [stat(*(E(x) for x in line.split(": ", 1))) for line in parascopy_lines(S["parascopy"]["data"])]
+               # stat escapes the label itself; the value is escaped here
+               + [stat(lab, E(val)) for lab, val 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)]))

docs/24-html-report.md

@@ -96,14 +96,15 @@ The cards, by chapter (cards not listed here are named in the chapter table abov
 - **Structural Variants** — Manta, Delly, CNVpytor and the SURVIVOR consensus (2+ callers) counts; then step 05's class 4 and 5 SVs (AnnotSV) with their caller support, exon-edge and frequency checks. A call step 05's rules moved to class 3 (a coding tandem repeat, an immunoglobulin or T-cell receptor locus, the MHC, a pseudogene named by a disease gene's alias, or read depth that does not show the change) is listed apart with AnnotSV's class, the rule and the reason, and the class 4-5 count leaves it out
 - **Repeat Expansions** — a table of every locus, the five key loci (HTT, FMR1, C9orf72, ATXN1, DMPK) first, with its repeat count, Stranger's status and notes; and counts of the loci outside the normal range (step 9b, Stranger), of those short reads could not size (not counted as normal), and the loci in an intermediate range (`bin/str_intermediate_ranges.tsv`). When Stranger flags RFC1, or another locus where the motif decides (`expansion_hunter/${SAMPLE}_str_motifs.tsv`), the card adds the verdict of the read-level motif count (step 9b). A locus the repeat catalog gives no thresholds (its `NoThresholds` field in `expansion_hunter/${SAMPLE}_eh_catalog.json`: MUC1, CEL, TYMS and PLIN4 with the STRchive catalog, among others) is marked not graded, with the catalog's reason, never normal. Any other locus outside Stranger's catalog has no status and is counted apart. Without step 9b the card says Stranger did not run; a Stranger file older than the ExpansionHunter result it came from is not read, and the card says to rerun step 9b. Repeat-expansion calls from short reads can be wrong at some loci
 - **Mitochondrial** — chrM PASS variants and heteroplasmic calls (allele fraction 0.05 to 0.95; below 5% NUMT reads and noise dominate)
+- **Lp(a): LPA KIV-2 Repeat Count** — when step 46 ran: the KIV-2 count with its interval and class, the KIV-2 4925G>A null allele under the count with its caveat (`not assessed` without step 46's 4925G>A table), step 42's two LPA risk variants, the flag and the blood-test advice
 - **Mosaic Loss of Chromosome Y or X** — step 33's call on chrY (men) or chrX (women), the reason for it, and a note that this is a population-level association, not a diagnosis
 
 **Carrier**
 
 - **Carrier List (this sample)** — the step 06 ClinVar hits in the recessive and X-linked genes of `assets/carrier_genes.tsv`, for this sample alone; hypomorph and pseudodeficiency alleles set aside
 - **Couple Carrier Screen** — when step 39 ran for this sample: the partner, the genes where both carry a change, the X-linked findings and the genes where only one partner does, each shared or X-linked gene with the chance for a child
 - **Paralog Genes (Paraphase)** — when step 35b ran: copy number and call per gene, SMN1 from Paraphase beside Parascopy, with whether they agree, when step 35 ran too, and the count of variants on phased haplotypes
-- **SMN1/SMN2 Copy Number (Parascopy)** — when step 35 ran: SMN1 and SMN2 copies at exon 7 with quality and filter, the two together (agCN), the SMA reading, and, when step 35 wrote its HBA table, the count of alpha-globin genes from depth or the reason it was not assessed
+- **SMN1/SMN2 Copy Number (Parascopy)** — when step 35 ran: SMN1 and SMN2 copies at exon 7 with quality and filter, the SMN1 g.27134T>G silent-carrier marker under the SMN1 copies with its caveat (`not assessed` without step 35's marker table), the two together (agCN), the SMA reading, and, when step 35 wrote its HBA table, the count of alpha-globin genes from depth or the reason it was not assessed
 - **Paralog Copy Number** — `not checked`, gray, when neither step 35 nor step 35b left a result, naming a step that failed
 
 **Medicines**

docs/35-paralogs.md

@@ -96,7 +96,7 @@ SMN2 holds the same 21 bases 875,420 bases upstream, so the reads cannot tell an
 | `too little depth` | under 10 reads over the two positions |
 | `not assessed` | the pileup could not be read, or the check failed |
 
-The message reads the marker beside the SMN1 psCN at exon 7 (filter `PASS` and quality 20 or more, the rule of the report's SMN1 line): with two SMN1 copies it says whether the marker raises or lowers the chance of 2+0; with another count, or a psCN that is not resolved, it says the marker refines a two-copy result only. Every message ends with "This marker does not prove 2+0, and its absence does not rule 2+0 out." The reports do not read this table yet; it sits beside the SMN table in the same folder.
+The message reads the marker beside the SMN1 psCN at exon 7 (filter `PASS` and quality 20 or more, the rule of the report's SMN1 line): with two SMN1 copies it says whether the marker raises or lowers the chance of 2+0; with another count, or a psCN that is not resolved, it says the marker refines a two-copy result only. Every message ends with "This marker does not prove 2+0, and its absence does not rule 2+0 out." Both reports and the live person page print the marker right under the SMN1 copies: its state (`present` or `uncertain` with the reads that carry G and the copies, `not seen` with the reads, `too little depth`, or `not assessed`), the reader's message, and the caveat that the marker raises the chance that two SMN1 copies sit on one chromosome, does not prove 2+0, and that a carrier test in a clinical laboratory confirms it. Without the table (a run from before this reader) or with a table that has no known status, the line says `not assessed` with the reason. The live card also lists the SMN1 copies and the marker state as facts. The Nextflow `HTML_REPORT` gets the table from the `PARASCOPY` task the same way.
 
 ## Runtime
 

docs/42-traits.md

@@ -115,8 +115,8 @@ The six rules after TAS2R38 need no phase. Their sites lie close together on one
 
 Two more markers of the planned list are not rules here, because each sits in sequence the genome holds more than once and a gVCF genotype at one position cannot read it. They are read from the BAM by the step that counts the copies, as a fraction of copies rather than a genotype (`MARKERS` in `bin/traits_report.py`, which pins each one's dbSNP b157 record and the GRCh38 sequence around every copy):
 
-- **LPA KIV-2 4925G>A** (Coassin et al. 2017, PMID 28444229), a splice variant that lowers the Lp(a) made from the chromosome that carries it. GRCh38 holds the site in each of its six KIV-2 units and in KIV-1. [Step 46](46-lpa.md) reads it beside the KIV-2 count.
-- **SMN1 g.27134T>G** (rs143838139; Luo et al. 2014, PMID 23788250), a marker of SMN1 duplication alleles that SMN2 holds too. [Step 35](35-paralogs.md) reads it beside the SMN1 copy number. It refines the chance of a silent 2+0 carrier for a two-copy result, and its line always says that it does not prove 2+0.
+- **LPA KIV-2 4925G>A** (Coassin et al. 2017, PMID 28444229), a splice variant that lowers the Lp(a) made from the chromosome that carries it. GRCh38 holds the site in each of its six KIV-2 units and in KIV-1. [Step 46](46-lpa.md) reads it beside the KIV-2 count, and the reports print it under the count.
+- **SMN1 g.27134T>G** (rs143838139; Luo et al. 2014, PMID 23788250), a marker of SMN1 duplication alleles that SMN2 holds too. [Step 35](35-paralogs.md) reads it beside the SMN1 copy number, and the reports print it under the SMN1 copies. It refines the chance of a silent 2+0 carrier for a two-copy result, and its line always says that it does not prove 2+0.
 
 Kim et al. (2003) found five TAS2R38 haplotypes worldwide. The rule reads the two common ones. A copy with AAV, AAI, PVI or a rarer one is `other_haplotype`, so the GRCh38 reference itself, which carries AAI (Ala49, Ala262, Ile296), is not read when every site matches it. The readings follow Bufe et al. (2005): PAV/PAV rated PTC and PROP the most bitter, PAV/AVI recognised them as bitter at nearly the same concentrations with a lower and more variable intensity, and AVI/AVI recognised PTC only at far higher concentrations.
 

docs/46-lpa.md

@@ -33,6 +33,7 @@ Pinned in `versions.env`; [Image versions](versions.md) lists the current tags.
 ## Input
 
 - `${SAMPLE}/aligned/${SAMPLE}_sorted.bam` with its `.bai` (`ALIGN_DIR` as for the other BAM steps), aligned to GRCh38 with `chr` contig names.
+- A failure of the 4925G>A reader never costs the count: step 46 and the `LPA` process then write a 4925G>A table that says `not assessed, the check failed`, print a warning, and go on.
 - The aligner must report a read that fits several places equally well, with mapping quality 0. minimap2 (step 2) and BWA-MEM do. novoalign by default leaves such reads unaligned: then the step finds almost no MAPQ 0 reads over KIV-2 and reports the count as not assessed rather than a count that comes out low.
 
 ## Command
@@ -105,7 +106,7 @@ What a short-read fraction cannot say: which chromosome or which repeat carries
 | `lpa/${SAMPLE}_lpa_kiv2_depth_q1.tsv` | the same, reads of MAPQ 1 or more |
 | `lpa/${SAMPLE}_lpa_kiv2_4925_pileup.tsv` | `samtools mpileup` at the seven 4925G>A positions |
 
-Both reports have an Lp(a) section: the count with its interval and class, the two risk variants from step 42 with their risk allele counts (rs10455872 G, rs3798220 C), a flag when there is at least one risk allele or the count is low, and the advice to measure Lp(a) in blood. Without step 42's table the two variants are shown as not read. The reports do not read the 4925G>A table yet; step 46 prints its message beside the count.
+Both reports have an Lp(a) section: the count with its interval and class, the two risk variants from step 42 with their risk allele counts (rs10455872 G, rs3798220 C), a flag when there is at least one risk allele or the count is low, and the advice to measure Lp(a) in blood. Without step 42's table the two variants are shown as not read. Right under the count, both reports and the live person page print the 4925G>A line: `present` or `uncertain` with the reads that carry it and the copies, `not seen` with the reads, `too little depth`, or `not assessed`; then the reader's message and the caveat that a null allele lowers the Lp(a) made from the copies it sits on and that the Lp(a) blood test is the measure. Without the 4925G>A table (a run from before this reader) or with a table that has no known status, the line says `not assessed` with the reason. The live card also lists the state as a fact. The Nextflow `HTML_REPORT` gets the table from the `LPA` process the same way.
 
 ## Runtime