fix(report): show the LPA null-allele and SMN1 marker lines everywhere, and say why CYP21A2 was not assessed on the couple screen

bin/couple_report.py

@@ -413,17 +413,34 @@ def smn_finding(path, genes):
     return [finding("SMN1", text, 1 if p["cn"] == 1 else 2)], "SMN1 %d cop%s" % (p["cn"], "y" if p["cn"] == 1 else "ies")
 
 
-def cyp21a2_finding(path, genes):
-    """CYP21A2 from depth and PSV sites (step 35's <id>_cyp21a2_psv.tsv):
-    ([findings], text of what was read). A deletion or a possible
-    conversion is a finding on one copy, on two when a site has no CYP21A2
-    copy left."""
+def cyp21a2_table(path):
+    """Step 35's CYP21A2 table (metric, value) as a dict."""
     t = {}
     with open(path, encoding="utf-8", errors="replace") as f:
         for line in f:
             c = line.rstrip("\n").split("\t", 1)
             if len(c) == 2 and c[0] != "metric":
                 t[c[0]] = c[1]
+    return t
+
+
+def cyp21a2_gap(path, txt):
+    """Why step 35 did not assess CYP21A2, for the not-assessed list: the
+    table's message (it says, for a BAM without multi-mapping reads, that the
+    aligner and not the genome stopped the call), or the short text with the
+    reason code for a table from before the message row."""
+    msg = cyp21a2_table(path).get("message", ".")
+    if msg.startswith("not assessed: "):
+        return "CYP21A2 " + msg
+    return txt
+
+
+def cyp21a2_finding(path, genes):
+    """CYP21A2 from depth and PSV sites (step 35's <id>_cyp21a2_psv.tsv):
+    ([findings], text of what was read). A deletion or a possible
+    conversion is a finding on one copy, on two when a site has no CYP21A2
+    copy left."""
+    t = cyp21a2_table(path)
     if t.get("status") not in ("called", "not assessed"):
         die(["%s: not step 35's CYP21A2 table (bin/copy_number.py call --locus CYP21A2)" % path])
     if t["status"] != "called":
@@ -835,7 +852,7 @@ def cmd_screen(a):
                 m, txt = cyp21a2_finding(cyp[s], genes)
                 f += m
                 if txt.startswith("CYP21A2 not assessed"):
-                    gaps.append("CYP21A2 deletion and conversions of %s: %s." % (s, txt))
+                    gaps.append("CYP21A2 deletion and conversions of %s: %s." % (s, cyp21a2_gap(cyp[s], txt)))
             else:
                 gaps.append("CYP21A2 deletion and conversions of %s: step 35's CYP21A2 table was not passed "
                             "to the couple screen (--cyp21a2)." % s)

docs/39-couple.md

@@ -25,7 +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 and the Nextflow couple screen (`COUPLE_CARRIER`) both pass it |
+| 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. A table whose `status` is `not assessed` gives no finding, and the list of what was not assessed prints its `message` (why, such as a BAM without multi-mapping reads), or its `reason` for a table without that row. Step 39 and the Nextflow couple screen (`COUPLE_CARRIER`) both pass it |
 | 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.