diff --git a/assets/acmg_sf.tsv b/assets/acmg_sf.tsv index b4856e36..9d7d27e6 100644 --- a/assets/acmg_sf.tsv +++ b/assets/acmg_sf.tsv @@ -9,7 +9,10 @@ # report_if when a P/LP variant (ClinVar, or a rare HIGH-impact record) is a secondary finding: # any one variant is enough # biallelic two variants: a homozygote, or two different variants even when -# their phase is unknown (Table 1 footnote a) +# their phase is unknown (Table 1 footnote a); a called deletion +# of the gene (AnnotSV, step 35's copy number) counts as one of the +# two unless read-backed phase puts it on the variant's copy +# (bin/collect_summary.py two_hits); so it does for hemizygous # hom_variant only a variant of the variants column, homozygous # hemizygous a hemizygous or homozygous variant, or two heterozygous ones; a # single heterozygous variant is not reported diff --git a/assets/carrier_genes.tsv b/assets/carrier_genes.tsv index 46668dca..7926c39b 100644 --- a/assets/carrier_genes.tsv +++ b/assets/carrier_genes.tsv @@ -14,7 +14,9 @@ # severity: alleles of the gene that are not a full loss of function, separated by ";": each is # hypomorph or pseudodeficiency, then chrom:pos:ref:alt (GRCh38), then a name; "." for none. Such # an allele never makes a carrier finding on its own (the couple screen, step 39, and the -# carrier list of the sample report). +# carrier list of the sample report). With a called deletion of the gene it is a second hit (the +# two-hit lines of the sample report, bin/collect_summary.py two_hits). In the couple screen a +# hypomorph beside a full loss of function in the other partner's same gene is a shared gene. # Rebuilt with: couple_report.py genes (docs/39-couple.md says how). gene chrom start end inheritance lists condition severity ISG15 chr1 1001138 1014540 AR Orphanet Mendelian susceptibility to mycobacterial diseases due to complete ISG15 deficiency . @@ -186,7 +188,7 @@ HSD3B2 chr1 119413725 119453802 AR Orphanet Congenital adrenal hyperplasia due t PHGDH chr1 119648411 119744227 AR Orphanet 3-phosphoglycerate dehydrogenase deficiency, infantile/juvenile form; and 1 other . HMGCS2 chr1 119747979 119769092 AR Orphanet 3-hydroxy-3-methylglutaryl-CoA synthase deficiency . PEX11B chr1 145911348 145918747 AR Orphanet Infantile Refsum disease; and 2 others . -RBM8A chr1 145920935 145928067 AR Orphanet Thrombocytopenia-absent radius syndrome . +RBM8A chr1 145920935 145928067 AR Orphanet Thrombocytopenia-absent radius syndrome hypomorph chr1:145927447:C:T c.-21G>A (rs139428292), a low-expression allele: TAR syndrome with a null allele such as the 1q21.1 deletion on the other copy; hypomorph chr1:145927328:C:G c.67+32G>C (rs201779890), the other low-expression allele of TAR syndrome HJV chr1 146017468 146036746 AR Orphanet HJV or HAMP-related hemochromatosis . GJA8 chr1 147902795 147909269 AR Orphanet Cataract-microcornea syndrome; and 3 others . NOTCH2NLC chr1 149390534 149471833 AR Orphanet Oculopharyngodistal myopathy . diff --git a/bin/collect_summary.py b/bin/collect_summary.py index 1431e5cf..7605d801 100755 --- a/bin/collect_summary.py +++ b/bin/collect_summary.py @@ -1777,7 +1777,10 @@ def sec_annotsv(d, s): except ValueError: af = None genes = [g for g in (r.get("Gene_name") or "").split(";") if g] + fmt = annotsv_postfilter.sample_format(r) out.append({"sv": f"{r['SV_chrom']}:{r['SV_start']}-{r['SV_end']}", "type": r["SV_type"], + # the call's genotype and phase set (PS), for the two-hit rule (two_hits) + "gt": fmt.get("GT") or None, "ps": fmt.get("PS") or None, "length": r.get("SV_length") or "", "acmg_class": r["ACMG_class"], "genes": genes, "caller": r["Caller"], "support": sp["Support"], "support_callers": sp["Support_callers"].split(",") if sp["Support_callers"] else [], @@ -2557,7 +2560,7 @@ def sec_carrier(d, s): p = first_existing(d, [f"clinvar/{s}_clinvar_hits.vcf"]) if not p: return None, {} - return p, carrier_list_or_note(p) + return p, carrier_list_or_note(p, d, s) def sec_exomiser(d, s): @@ -3236,7 +3239,7 @@ def csq_high_impact_genes(path): gt = r[9].split(":")[0] if len(r) > 9 else "" weak, why, gq, dp, ad = call_check(r[8], r[9], floor) if len(r) > 9 else (False, [], None, None, None) out[(r[0], int(r[1]), r[3], r[4])] = { - "genotype": gt, "genes": genes, "transcript": cls, "frequency": freq, + "genotype": gt, "ps": gt_field(r, "PS"), "genes": genes, "transcript": cls, "frequency": freq, "weak": weak, "why": why, "gq": gq, "dp": dp, "ad": ad} except (OSError, EOFError, gzip.BadGzipFile, zlib.error) as e: raise Unreadable(f"{os.path.basename(path)}: {e}") from e @@ -3327,7 +3330,8 @@ def carrier_genes(): for n, r in enumerate(asset_rows("carrier_genes.tsv"), 1): try: out[r["gene"]] = {"inheritance": r["inheritance"], "lists": r["lists"], "condition": r["condition"], - "severity": parse_alleles(r.get("severity"), kinds=True)} + "severity": parse_alleles(r.get("severity"), kinds=True), + "span": (r["chrom"], int(r["start"]), int(r["end"]))} except (KeyError, ValueError, TypeError) as e: raise Unreadable(f"carrier_genes.tsv row {n}: {e}") from e return out @@ -3347,6 +3351,169 @@ def zygosity_word(copies, ploidy): return "hom" if copies >= 2 else "het" +# --- two hits in one gene across callers: a deletion plus a small variant -------- +# +# Small variants alone miss a gene hit twice when one hit is a deletion: an +# exon deletion of one copy of MUTYH plus a pathogenic variant on the other, an +# HBA deletion plus an HBA point variant, one SMN1 copy plus an SMN1 variant, a +# 1q21.1 deletion plus the RBM8A low-expression allele (TAR syndrome: Albers et +# al., Nat Genet 2012;44(4):435-439, PMID 22366785). The deletions come from the +# engines that call them: AnnotSV's class 4 and 5 deletions (step 05), step 35's +# interval table (bin/cn_loci.tsv), its HBA count and Parascopy's SMN1 count. +# Two hits are in trans (biallelic) only when the reads show it: the variant sits +# inside the deletion and is called homozygous (only one copy is left, and it +# carries the variant), or both calls carry a read-backed phase (GT with "|" and +# the same PS) on different copies. Otherwise they are possibly in trans; on the +# same copy (cis) they are no second hit. +TWO_HIT_ORDER = ("biallelic", "possibly biallelic") + + +def gene_deletions(d, s, spans): + """The genes a called deletion reaches: ({gene: [deletion]}, engines read, + engines not read). spans: {gene: (chrom, start, end)}, 1-based inclusive + (assets/carrier_genes.tsv). A deletion: chrom, start, end (1-based, + inclusive), text, inside (a variant inside its span is on the one copy left; + False for the HBA count, which cannot say which gene of the pair is gone), + gt and ps (its phase, AnnotSV only).""" + out, read, unread = {}, [], [] + + def add(g, x): + out.setdefault(g, []).append(x) + + def over(chrom, a, b): + return [g for g, (c, gs, ge) in spans.items() if c == chrom and gs <= b and ge >= a] + try: + p, data = sec_annotsv(d, s) + except (Unreadable, OSError, ValueError, KeyError) as e: + p, data = None, {} + unread.append(f"AnnotSV (step 05): {e}") + if p: + read.append("AnnotSV class 4 and 5 deletions (step 05)") + for x in data["rows"]: + chrom, _, span = x["sv"].partition(":") + a, _, b = span.partition("-") + if x["type"] != "DEL" or x["common"] or not (a.isdigit() and b.isdigit()): + continue + chrom = chrom if chrom.startswith("chr") else "chr" + chrom + # SV_start is the base before the deleted bases (annotsv_postfilter.sv_span): + # they are SV_start+1 to SV_end, 1-based. The text keeps AnnotSV's coordinates. + first = int(a) + 1 + hit = set(over(chrom, first, int(b))) + # A listed gene AnnotSV names counts only when its span reaches the deleted bases: one + # that ends at the padding base loses none. A gene with no span here is taken as named. + named = {g for g in x["genes"] if g not in spans} + for g in sorted(hit | named): + add(g, {"chrom": chrom, "start": first, "end": int(b), "inside": True, + "text": f"DEL {chrom}:{a}-{b} (AnnotSV class {x['acmg_class']}, {x['caller']})", + "gt": x.get("gt"), "ps": x.get("ps")}) + elif not unread: + unread.append("AnnotSV (step 05): no table") + try: + p, data = sec_parascopy(d, s) + why = "none" + except (Unreadable, OSError, ValueError, KeyError) as e: + p, data, why = None, {}, str(e) + if not p: + unread.append(f"step 35's copy-number tables: {why}") + return out, read, unread + iv = data.get("cn_intervals") + if iv is not None: + read.append("step 35's interval table (read depth)") + for r in iv["changes"]: + if r["change"] != "deletion" or not (r["start"].isdigit() and r["end"].isdigit()): + continue + a, b = int(r["start"]) + 1, int(r["end"]) # the table is BED: 0-based start + for g in over(r["chrom"], a, b): + add(g, {"chrom": r["chrom"], "start": a, "end": b, "inside": True, "gt": None, "ps": None, + "text": f"{r['locus']} deletion, {r['copies']} of {r['expected']} copies " + "(read depth, step 35)"}) + else: + unread.append("step 35's interval table: none") + hba = data.get("hba") or {} + n = hba.get("alpha_globin_genes", "") + if hba.get("status") == "called" and n.isdigit(): + read.append("step 35's HBA1/HBA2 count") + # 0 genes: no alpha-globin gene is left to carry a variant, so a hit there pairs with nothing. + if 1 <= int(n) < 4: + for g in ("HBA1", "HBA2"): + if g in spans: + c, a, b = spans[g] + add(g, {"chrom": c, "start": a, "end": b, "inside": False, "gt": None, "ps": None, + "text": f"{n} of 4 alpha-globin genes (read depth, step 35)"}) + else: + unread.append(f"step 35's HBA1/HBA2 count: {hba.get('message') or 'not called'}") + smn1 = (data.get("smn1") or {}).get("cn") + if smn1 is not None: + read.append("Parascopy's SMN1 count (step 35)") + if smn1 == 1 and "SMN1" in spans: + c, a, b = spans["SMN1"] + add("SMN1", {"chrom": c, "start": a, "end": b, "inside": True, "gt": None, "ps": None, + "text": "one SMN1 copy (Parascopy, step 35)"}) + else: + unread.append("Parascopy's SMN1 count: not resolved") + return out, read, unread + + +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.""" + if not gt or "|" not in gt: + return None + return frozenset(i for i, a in enumerate(gt.split("|")) if a not in ("0", ".")) + + +def pair_phase(dl, v): + """(verdict, why) of a deletion and a small variant in one gene: verdict + biallelic, possibly biallelic, or cis (no second hit).""" + inside = dl["inside"] and dl["chrom"] == v["chrom"] and dl["start"] <= v["pos"] <= dl["end"] + # Read-backed phase first: on the same copy (cis) is no second hit, wherever the variant sits. + hd, hv = phase_haps(dl.get("gt")), phase_haps(v.get("gt")) + phased = None + if hd and hv and dl.get("ps") and dl.get("ps") == v.get("ps"): + if hd == hv: + return "cis", f"in cis: read-backed phase puts both on the same copy (phase set {v['ps']})" + phased = (f"in trans: read-backed phase puts them on different copies (phase set {v['ps']})" if not hd & hv + else f"read-backed phase puts the variant on both copies (phase set {v['ps']}); check the reads") + if inside and v["copies"] >= 2: + return "biallelic", ("the variant is called homozygous inside the deletion: it is hemizygous, on the one " + "copy the deletion leaves") + if inside and v["ploidy"] == 2: + return "possibly biallelic", ("a heterozygous call inside a called deletion: one of the two calls may be " + "wrong; check the reads") + if phased: + return ("biallelic" if phased.startswith("in trans") else "possibly biallelic"), phased + return "possibly biallelic", "phase unknown, possibly in trans (no read-backed phase links the two calls)" + + +def two_hits(variants, dels): + """The genes hit twice: ({gene: entry}, {(gene, key): deletion text} of + the variants called homozygous inside a deletion, shown as hemizygous). + variants: [{gene, key, chrom, pos, text, genotype, copies, ploidy, gt, ps, + allele}] (allele: the severity name of a hypomorph, or None). An entry is + the best pair of the gene: verdict, why, deletion, variant, genotype, and + the allele name; a gene whose every pair is in cis has none. pairs: + {(gene, key): " with : "}, the best pair of each + variant that has one.""" + entries, hemi, pairs = {}, {}, {} + for v in variants: + for dl in dels.get(v["gene"], []): + verdict, why = pair_phase(dl, v) + if verdict == "biallelic" and why.startswith("the variant is called homozygous"): + hemi[(v["gene"], v["key"])] = dl["text"] + if verdict == "cis": + continue + e = {"gene": v["gene"], "verdict": verdict, "why": why, "deletion": dl["text"], "variant": v["text"], + "genotype": "hemizygous" if (v["gene"], v["key"]) in hemi else v["genotype"], + "allele": v.get("allele")} + old = entries.get(v["gene"]) + if old is None or TWO_HIT_ORDER.index(verdict) < TWO_HIT_ORDER.index(old["verdict"]): + entries[v["gene"]] = e + k = (v["gene"], v["key"]) + if k not in pairs or TWO_HIT_ORDER.index(verdict) < TWO_HIT_ORDER.index(pairs[k][0]): + pairs[k] = (verdict, f"{verdict} with {dl['text']}: {why}") + return entries, hemi, {k: t for k, (_, t) in pairs.items()} + + def sf_verdicts(gene, rule, alleles): """Which of a gene's variants the ACMG rule reports. alleles: {key: (copies, ploidy)} of the gene's distinct small variants, structural ones @@ -3443,7 +3610,8 @@ def acmg_sf_tier(d, s, clinical_vcf): "genotype": row["genotype"], "significance": row["significance"], "stars": row["stars"], "low_confidence": row["low_confidence"], "key": (row["chrom"], row["pos"], row["ref"], row["alt"]), - "copies": copies, "ploidy": ploidy}) + "copies": copies, "ploidy": ploidy, + "gt": r[9].split(":")[0] if len(r) > 9 else "", "ps": gt_field(r, "PS")}) if clinical_vcf: # A common record is in the clinical VCF through the ClinVar tier, # not as a rare HIGH impact; one of unknown frequency is listed and @@ -3461,7 +3629,8 @@ def acmg_sf_tier(d, s, clinical_vcf): "frequency_unknown": x["frequency"] == "unknown", "weak_call": x["weak"], "call_note": "; ".join(x["why"]) or None, "gq": x["gq"], "dp": x["dp"], "ad": x["ad"], - "key": (c, pos, ref, alt), "copies": copies, "ploidy": ploidy}) + "key": (c, pos, ref, alt), "copies": copies, "ploidy": ploidy, + "gt": gt, "ps": x.get("ps")}) svs, sv_checked, sv_note = sf_structural(d, s, rules) # One verdict per gene, over its distinct variants: a ClinVar hit that is # also a HIGH-impact record is one variant. @@ -3470,6 +3639,30 @@ def acmg_sf_tier(d, s, clinical_vcf): alleles.setdefault(x["gene"], {}).setdefault(x["key"], (max(x["copies"], 1), x["ploidy"] or 2)) for x in svs: alleles.setdefault(x["gene"], {}).setdefault(x["key"], (1, 2)) + # A called deletion of a gene whose rule needs two variants counts as one + # of them (two_hits), unless read-backed phase puts it on the variant's copy. + two_hit_note = None + try: + spans = {g: x["span"] for g, x in carrier_genes().items()} + except Unreadable as e: + spans = {} + two_hit_note = (f"A deletion as the second hit: only AnnotSV's gene names were read; step 35's deletions " + f"need the gene spans of {e}") + dels, _, _ = gene_deletions(d, s, spans) + seen, variants = set(), [] + for x in clinvar + high: + pseudo = rules[x["gene"]]["severity"].get(x["key"], {}).get("kind") == "pseudodeficiency" + if (x["gene"], x["key"]) not in seen and not pseudo: + seen.add((x["gene"], x["key"])) + variants.append({"gene": x["gene"], "key": x["key"], "chrom": x["key"][0], "pos": x["key"][1], + "text": x["variant"], "genotype": x["genotype"], "copies": x["copies"], + "ploidy": x["ploidy"] or 2, "gt": x["gt"], "ps": x["ps"], "allele": None}) + entries, hemi, pairs = two_hits(variants, dels) + two = {} + for g in entries: + if rules[g]["report_if"] in ("biallelic", "hemizygous"): + alleles[g][("deletion", g)] = (1, 2) + two.update({k: t for k, t in pairs.items() if k[0] == g}) verdict = {g: sf_verdicts(g, rules[g], a) for g, a in alleles.items()} reported = {"clinvar": [], "high": [], "sv": []} not_reported = [] @@ -3477,10 +3670,15 @@ def acmg_sf_tier(d, s, clinical_vcf): for x in rows: why, partial = verdict[x["gene"]][0][x["key"]], verdict[x["gene"]][1] rule = rules[x["gene"]] - for k in ("key", "copies", "ploidy"): + key = x["key"] + if (x["gene"], key) in hemi: + x["genotype"] = "hemizygous" + for k in ("key", "copies", "ploidy", "gt", "ps"): x.pop(k, None) if why is None: x.update(rule=rule["table_text"], rule_note=rule["note"], severity=partial) + if kind != "sv" and (x["gene"], key) in two: + x["two_hit"] = two[(x["gene"], key)] reported[kind].append(x) else: not_reported.append({"gene": x["gene"], "variant": x["variant"], "genotype": x.get("genotype"), @@ -3494,21 +3692,33 @@ def acmg_sf_tier(d, s, clinical_vcf): "high_impact_weak": sum(1 for x in high if x["weak_call"]), "high_impact_frequency_unknown": sum(1 for x in high if x["frequency_unknown"]), "structural": reported["sv"], "structural_checked": sv_checked, "structural_note": sv_note, - "not_reported": not_reported} + "not_reported": not_reported, "two_hit_note": two_hit_note} + + +def gt_field(r, key): + """One FORMAT field of a VCF record's first sample, or None.""" + if len(r) < 10: + return None + v = dict(zip(r[8].split(":"), r[9].split(":"))).get(key) + return None if v in (None, "", ".") else v -def carrier_list(p): +def carrier_list(p, d=None, s=None): """The sample's own carrier list: the ClinVar hits of step 6 (p) in the genes of assets/carrier_genes.tsv, one person, no partner needed. An allele the list marks hypomorph or pseudodeficiency is set aside, never a - carrier finding on its own.""" + carrier finding on its own. With the sample folder (d, s), the genes of + autosomal recessive conditions that a called deletion and a hit both reach + (two_hits), and a hit called homozygous inside a deletion is shown as + hemizygous.""" genes = carrier_genes() - carriers, aside = [], [] + carriers, aside, variants = [], [], [] for r in vcf_records(p): if len(r) < 8: continue row = clinvar_row(r) - copies, ploidy = allele_copies(r[9].split(":")[0] if len(r) > 9 else "") + gt = r[9].split(":")[0] if len(r) > 9 else "" + copies, ploidy = allele_copies(gt) key = (row["chrom"], row["pos"], row["ref"], row["alt"]) for g in row["gene"].split(","): if g not in genes: @@ -3523,8 +3733,27 @@ def carrier_list(p): aside.append(item) else: carriers.append(item) - return {"genes_on_list": len(genes), "list": "assets/carrier_genes.tsv", - "carriers": carriers, "set_aside": aside} + if not sev or sev["kind"] == "hypomorph": + variants.append({"gene": g, "key": key, "chrom": row["chrom"], "pos": row["pos"], + "text": item["variant"], "genotype": item["genotype"], "copies": copies, + "ploidy": ploidy, "gt": gt, "ps": gt_field(r, "PS"), + "allele": sev["name"] if sev else None, "item": item}) + out = {"genes_on_list": len(genes), "list": "assets/carrier_genes.tsv", + "carriers": carriers, "set_aside": aside} + if d is not None: + dels, read, unread = gene_deletions(d, s, {g: x["span"] for g, x in genes.items()}) + entries, hemi, _ = two_hits(variants, dels) + for v in variants: + if (v["gene"], v["key"]) in hemi: + v["item"].update(genotype="hemizygous", + hemizygous=f"called homozygous inside {hemi[(v['gene'], v['key'])]}") + rows = [] + for g, e in sorted(entries.items()): + if genes[g]["inheritance"] == "AR": + e.update(condition=genes[g]["condition"], inheritance="AR") + rows.append(e) + out["two_hits"] = {"genes": rows, "deletions_read": read, "deletions_not_read": unread} + return out def acmg_sf_or_note(d, s, clinical_vcf): @@ -3537,11 +3766,11 @@ def acmg_sf_or_note(d, s, clinical_vcf): "structural": [], "not_reported": []} -def carrier_list_or_note(p): +def carrier_list_or_note(p, d=None, s=None): """carrier_list, or the reason it could not run, which both reports print: the ClinVar hits stay.""" try: - return carrier_list(p) + return carrier_list(p, d, s) except Unreadable as e: return {"unavailable": str(e)} diff --git a/bin/couple_report.py b/bin/couple_report.py index 93fcb820..9ae97440 100755 --- a/bin/couple_report.py +++ b/bin/couple_report.py @@ -234,9 +234,12 @@ def common_match(benign, chrom, s0, e0): # --- findings ---------------------------------------------------------------- -def finding(gene, text, copies): - """copies: how many of the person's copies the finding is on (1 or 2).""" - return {"gene": gene, "text": text, "copies": copies} +def finding(gene, text, copies, hypomorph=None): + """copies: how many of the person's copies the finding is on (1 or 2). + hypomorph: the name of a low-function allele of the list's severity + column, which counts only beside a full loss of function in the partner + (screen_couple).""" + return {"gene": gene, "text": text, "copies": copies, "hypomorph": hypomorph} def zygosity(gt, chrom, sex, par): @@ -259,7 +262,11 @@ def stars_text(n): def clinvar_findings(path, genes, sex, par, aside=None): """Step 06's hits in genes of the list, and the count of all hits. A hit the list marks as a hypomorph or pseudodeficiency allele (its severity - column) is no carrier finding: it goes to aside as a text.""" + column) is no carrier finding on its own: it goes to aside as a text. A + hypomorph is also returned as a finding marked hypomorph: beside a full + loss of function in the other partner's same gene it is the second hit of + the child (TAR syndrome is a null allele plus an RBM8A low-expression + allele); screen_couple drops it otherwise. total counts every hit.""" out, total = [], 0 for r in collect_summary.vcf_records(path): if len(r) < 8: @@ -276,6 +283,10 @@ def clinvar_findings(path, genes, sex, par, aside=None): if aside is not None: aside.append("%s:%d %s>%s %s in %s: %s, a %s allele" % ( h["chrom"], h["pos"], h["ref"], h["alt"], z, g, sev["name"], sev["kind"])) + if sev["kind"] == "hypomorph": + text = "%s:%d %s>%s %s, %s, a hypomorph allele" % (h["chrom"], h["pos"], h["ref"], h["alt"], + z, sev["name"]) + out.append(finding(g, text, 1 if z == "het" else 2, sev["name"])) continue if g in genes: ident = "" if r[2] in (".", "") else ", ClinVar %s" % r[2] @@ -485,6 +496,11 @@ def ar_row(g, a, b, fa, fb): more = [n for n, f in ((a, fa), (b, fb)) if len(f) > 1 and all(x["copies"] == 1 for x in f)] why = ("Both partners carry a change in %s, linked to %s (autosomal recessive). A child who inherits " "the change of each parent can have the condition: %s." % (g["gene"], g["condition"], risk)) + low = [n for n, f in ((a, fa), (b, fb)) if f and all(x.get("hypomorph") for x in f)] + if low: + why += (" The change of %s is a low-function (hypomorph) allele, not a full loss of function: with the " + "other parent's loss-of-function change a child can have the condition, for some genes in a " + "milder form; two such alleles alone are set aside." % " and ".join(low)) if more: why += (" %s %s two changes in the gene; when they sit on different copies, every child gets one, " "and the chance is higher." % (" and ".join(more), "has" if len(more) == 1 else "have")) @@ -520,10 +536,13 @@ def xl_lines(g, who, sex, f): def screen_couple(a, b, sexes, found, genes): """Rows of the couple table, shared genes first, then X-linked, then the rest.""" rows = [] - for name in sorted(set(x["gene"] for x in found[a] + found[b])): + # A hypomorph (finding marked hypomorph) counts only when the other + # partner has a full loss of function in the same gene. + full = {s: set(x["gene"] for x in found[s] if not x.get("hypomorph")) for s in (a, b)} + for name in sorted(full[a] | full[b]): g = genes[name] - fa = [x for x in found[a] if x["gene"] == name] - fb = [x for x in found[b] if x["gene"] == name] + fa = [x for x in found[a] if x["gene"] == name and (not x.get("hypomorph") or name in full[b])] + fb = [x for x in found[b] if x["gene"] == name and (not x.get("hypomorph") or name in full[a])] if g["inheritance"] == "XL": parts = [xl_lines(g, w, sexes.get(w), f) for w, f in ((a, fa), (b, fb)) if f] result, risk, why = "x_linked", "; ".join(p[0] for p in parts), " ".join(p[1] for p in parts) @@ -686,7 +705,7 @@ def cmd_screen(a): gaps.append("Carrier deletion alleles of %s: no interval table (step 35)." % s) aside = [] f, total = clinvar_findings(clinvar[s], genes, sexes[s], par, aside) - info = {"clinvar": "%d / %d" % (len(f), total)} + info = {"clinvar": "%d / %d" % (sum(1 for x in f if not x["hypomorph"]), total)} for note in aside: gaps.append("Set aside for %s, not a carrier finding on its own: %s." % (s, note)) if s not in delly: @@ -860,7 +879,9 @@ def place(g): "# severity: alleles of the gene that are not a full loss of function, separated by \";\": each is", "# hypomorph or pseudodeficiency, then chrom:pos:ref:alt (GRCh38), then a name; \".\" for none. Such", "# an allele never makes a carrier finding on its own (the couple screen, step 39, and the", - "# carrier list of the sample report).", + "# carrier list of the sample report). With a called deletion of the gene it is a second hit (the", + "# two-hit lines of the sample report, bin/collect_summary.py two_hits). In the couple screen a", + "# hypomorph beside a full loss of function in the other partner's same gene is a shared gene.", "# Rebuilt with: couple_report.py genes (docs/39-couple.md says how).", ] with open(a.out, "w", encoding="utf-8") as f: diff --git a/bin/render_report.py b/bin/render_report.py index 9286436d..c8898ee1 100755 --- a/bin/render_report.py +++ b/bin/render_report.py @@ -1635,7 +1635,7 @@ def sf_rows(sf): the same in both reports: ClinVar hits, rare HIGH-impact records, then class 4 and 5 SVs.""" def extra(x): - return "".join(f"; {v}" for v in (x.get("severity"), x.get("rule_note")) if v) + return "".join(f"; {v}" for v in (x.get("severity"), x.get("rule_note"), x.get("two_hit")) if v) rows = [(x["gene"], x["variant"], x["genotype"], "ClinVar " + x["significance"] + ("; low-confidence call: check the reads" if x.get("low_confidence") else "") + extra(x)) for x in sf["clinvar_hits"]] @@ -1655,8 +1655,9 @@ def sf_apart_lines(sf): no rule could be applied.""" if sf.get("unavailable"): return [f"Secondary findings: not assessed, {sf['unavailable']}"] - return [f"{x['gene']} {x['variant']} {x.get('genotype') or '.'} ({x['source']}): not reported, {x['reason']}" - for x in sf.get("not_reported", [])] + return ([f"{x['gene']} {x['variant']} {x.get('genotype') or '.'} ({x['source']}): not reported, {x['reason']}" + for x in sf.get("not_reported", [])] + + ([sf["two_hit_note"]] if sf.get("two_hit_note") else [])) def carrier_lines(c): @@ -1665,12 +1666,31 @@ def carrier_lines(c): return [f"Carrier list: not made, {c['unavailable']}"] lines = [f"Carrier list ({c['genes_on_list']} genes): {len(c['carriers'])} finding(s)" + (f", {len(c['set_aside'])} set aside" if c["set_aside"] else "")] + def hemi(x): + return f" ({x['hemizygous']})" if x.get("hemizygous") else "" for x in c["carriers"]: - lines.append(f"{x['gene']} ({x['condition']}; {x['inheritance']}): {x['variant']} {x['genotype']}, " + lines.append(f"{x['gene']} ({x['condition']}; {x['inheritance']}): {x['variant']} {x['genotype']}{hemi(x)}, " f"ClinVar {x['significance']}") for x in c["set_aside"]: - lines.append(f"Set aside: {x['gene']} {x['variant']} {x['genotype']}: {x['name']} is a {x['kind']} allele, " - "not a carrier finding on its own") + lines.append(f"Set aside: {x['gene']} {x['variant']} {x['genotype']}{hemi(x)}: {x['name']} is a {x['kind']} " + "allele, not a carrier finding on its own") + return lines + two_hit_lines(c.get("two_hits")) + + +def two_hit_lines(t): + """The genes a called deletion and a ClinVar hit both reach (collect_summary.two_hits), and which deletion + callers were read; none for a summary written before the rule.""" + if not t: + return [] + lines = [] + for x in t["genes"]: + allele = f" ({x['allele']})" if x.get("allele") else "" + lines.append(f"{x['verdict'].capitalize()}: {x['gene']} ({x['condition']}; AR), both copies may be hit: " + f"{x['deletion']} and {x['variant']} {x['genotype']}{allele}; {x['why']}") + lines.append("Two hits across callers (a deletion plus a variant on the other copy): " + + (f"{len(t['genes'])} gene(s); " if t["genes"] else "none found; ") + + "deletions read from " + ("; ".join(t["deletions_read"]) or "no caller") + + (f"; not read: {'; '.join(t['deletions_not_read'])}" if t["deletions_not_read"] else "")) return lines diff --git a/docs/interpreting-results.md b/docs/interpreting-results.md index 435bc36e..af3b9631 100644 --- a/docs/interpreting-results.md +++ b/docs/interpreting-results.md @@ -76,7 +76,15 @@ The most common "pathogenic" finding in any genome is **heterozygous carrier sta - BTD with p.Asp444His: reported as partial, not profound, biotinidase deficiency; - AnnotSV class 4 and 5 structural variants (step 5) in BRCA1, BRCA2, MLH1, MSH2 and LDLR reach the same card, and so does a deletion of the 3' end of EPCAM (exons 8 and 9), which silences MSH2 (Kovacs et al., Hum Mutat 2009;30(2):197-203). -A deletion in trans with a small variant in a recessive gene is not counted as two variants yet. The TTN rule (truncating variants only) and the LMNA rule (missense only with cardiac evidence) are printed beside the finding, not applied. It is still a list to review with a clinician. +**Two hits across callers: a deletion plus a variant on the other copy.** Small-variant callers see one hit when the other copy of the gene is deleted. So the carrier list also joins the deletions the other engines call with the ClinVar hits, gene by gene (`bin/collect_summary.py`, `two_hits`). The deletions come from AnnotSV's class 4 and 5 deletions (step 5), step 35's interval table (the recurrent regions of `bin/cn_loci.tsv`, such as 1q21.1), its HBA1/HBA2 count and Parascopy's SMN1 count. A recessive gene that a deletion and a hit both reach gets one line: + +- **biallelic** when the reads show the two hits on different copies. Either the hit is called homozygous inside the deletion, so it sits on the one copy left (the line and the carrier list call it hemizygous), or both calls carry read-backed phase (a `|` genotype with the same phase set, PS) on different copies; +- **possibly biallelic** otherwise: the phase is unknown. Most short-read calls carry no phase, so this is the usual line. Step 05's checked table keeps no genotype column yet, and step 35's counts have no phase, so in a pipeline run a deletion pairs as possibly biallelic unless the hit is called homozygous inside it. A parent's genome (step 38, [family analysis](multi-sample.md#family-analysis-step-38-opt-in)) or long reads can settle it; +- no line when read-backed phase puts both on the same copy (cis): the other copy is intact. + +Examples: an exon deletion of MUTYH plus a pathogenic MUTYH variant; three alpha-globin genes plus an HBA2 point variant; one SMN1 copy plus an SMN1 variant (short reads place few SMN1 variants reliably, so check with a targeted test); and the 1q21.1 deletion plus one of the two RBM8A low-expression alleles (c.-21G>A, rs139428292, or c.67+32G>C, rs201779890), the cause of TAR syndrome (thrombocytopenia with absent radius; Albers et al., Nat Genet 2012;44(4):435-439, PMID 22366785). On its own an RBM8A low-expression allele stays set aside. The [couple screen](39-couple.md) (step 39) applies the same pairing between partners: a hypomorph allele of the list in one partner and a full loss of function in the same gene in the other (a ClinVar pathogenic variant, a Delly deletion) is a shared gene with its risk for a child; the allele in both partners, or in one alone, stays set aside. In the secondary-findings card a called deletion counts as one of the two variants a recessive gene needs, unless read-backed phase puts it on the variant's copy. The line under the carrier list names the deletion callers it read and the ones it could not read. A deletion AnnotSV ranks below class 4 is not counted. + +The TTN rule (truncating variants only) and the LMNA rule (missense only with cardiac evidence) are printed beside the finding, not applied. It is still a list to review with a clinician. --- diff --git a/tests/test_cross_engine_biallelic.py b/tests/test_cross_engine_biallelic.py new file mode 100644 index 00000000..777232b4 --- /dev/null +++ b/tests/test_cross_engine_biallelic.py @@ -0,0 +1,330 @@ +#!/usr/bin/env python3 +"""Two hits in one gene across callers (bin/collect_summary.py two_hits): a +called deletion plus a small variant, on synthetic sample folders. Nothing +here describes a person: the calls are made up at real gene positions. + +Checks: + 1. MUTYH: an AnnotSV exon deletion plus a heterozygous ClinVar P/LP hit + elsewhere in the gene gives a "possibly biallelic" line, and the ACMG SF + tier now reports the hit (MUTYH needs two variants: the deletion is one); + 2. the same two with read-backed phase on different copies (GT with "|", + same PS) give "biallelic"; + 3. the same two in cis (same copy, same PS) give no two-hit line, and the SF + tier keeps the one hit apart; + 4. TAR syndrome: step 35's 1q21.1 deletion plus the RBM8A low-expression + allele (rs139428292, assets/carrier_genes.tsv) gives the TAR line; + 5. a hit called 1/1 inside a called deletion is shown as hemizygous, in the + carrier list and in the SF tier; a 1/1 hit at AnnotSV's SV_start (the + base before the deleted bases) stays homozygous; + 6. SMN1: one SMN1 copy (Parascopy) plus an SMN1 hit, and HBA: three alpha + genes plus an HBA2 hit, each give a "possibly biallelic" line; + 7. without any deletion caller's output the hit gives no line, and the + summary names the callers it could not read; + 8. both reports print the two-hit line and the hemizygous genotype; + 9. the couple screen (step 39): the RBM8A low-expression allele in one + partner and a deletion of RBM8A in the other is a shared gene (TAR risk + for a child); the allele in both partners, or in one alone, is set aside; + the assessed table counts the allele as no finding (0 of 1 hit); + 10. edges: a gene AnnotSV names that ends at the padding base gets no line; + 0 alpha-globin genes pair with nothing; a het hit inside a deletion + phased in cis is no second hit; a HIGH-impact record keeps its phase set + (trans gives biallelic); with a cis and a trans hit in one gene only the + trans one carries the SF two-hit note; an unreadable carrier list is a + note on the SF card. +""" +import gzip +import os +import shutil +import subprocess +import sys +import tempfile + +REPO = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) +sys.path.insert(0, os.path.join(REPO, "bin")) +import collect_summary # noqa: E402 +import couple_report # noqa: E402 +import render_report # noqa: E402 + +FAILS = [] + + +def check(name, ok, detail=""): + print(("[PASS] " if ok else "[FAIL] ") + name + ("" if ok else " -- " + str(detail)[:1500])) + if not ok: + FAILS.append(name) + + +def put(path, text, gz=False): + os.makedirs(os.path.dirname(path), exist_ok=True) + with (gzip.open(path, "wt") if gz else open(path, "w")) as f: + f.write(text) + return path + + +VCF_HEAD = "##fileformat=VCFv4.2\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tS\n" +ANNOTSV_COLS = ("SV_chrom\tSV_start\tSV_end\tSV_length\tSV_type\tGene_name\tACMG_class\tCaller\tExon_edge\t" + "Exon_edge_detail\tPopulation_AF\tPopulation_source\tCommon\tFORMAT\tS\n") +# MUTYH spans chr1:45329163-45340893 (assets/carrier_genes.tsv). +SNV = ("chr1", 45332000, "C", "T") +DEL = ("chr1", 45335000, 45337000) +INSIDE = ("chr1", 45336000, "G", "A") +# RBM8A c.-21G>A (rs139428292), on the plus strand C>T. +TAR_ALLELE = ("chr1", 145927447, "C", "T") + + +def hit(v, gene, gt, ps=None): + chrom, pos, ref, alt = v + fmt, val = ("GT:GQ:DP:PS", f"{gt}:50:30:{ps}") if ps else ("GT:GQ:DP", f"{gt}:50:30") + return (f"{chrom}\t{pos}\t.\t{ref}\t{alt}\t50\tPASS\tGENEINFO={gene}:1;CLNSIG=Pathogenic;" + f"CLNREVSTAT=criteria_provided,_multiple_submitters,_no_conflicts\t{fmt}\t{val}\n") + + +def sv(span, genes, gt="0/1", ps=None): + chrom, a, b = span + fmt, val = ("GT:PS", f"{gt}:{ps}") if ps else ("GT", gt) + return f"{chrom}\t{a}\t{b}\t-{b - a}\tDEL\t{genes}\t4\tmanta\tyes\t.\t\t\tno\t{fmt}\t{val}\n" + + +SMN_COLS = "chrom\tstart\tend\tagCN\tpsCN\tpsCN_qual\tpsCN_filter\tagCN_qual\tagCN_filter\thomologous_regions\n" +INTERVAL_COLS = "locus\tname\tchrom\tstart\tend\tstatus\tcopies\texpected\tchange\testimate\treason\n" + + +def step35(d, smn1="2", intervals="", hba_genes=None): + """Step 35's tables: Parascopy's SMN table (SMN1 exon 7 at chr5:70951946, + SMN2's at chr5:70076526), the interval table and the HBA count.""" + put(f"{d}/paralogs/S_smn_copy_number.tsv", SMN_COLS + + f"chr5\t70951000\t70952000\t4\t{smn1},2\t40,40\tPASS\t40\tPASS\tchr5:70076000-70077000\n") + if intervals: + put(f"{d}/paralogs/S_cn_intervals.tsv", INTERVAL_COLS + intervals) + if hba_genes is not None: + put(f"{d}/paralogs/S_hba_copy_number.tsv", + f"metric\tvalue\nstatus\tcalled\nalpha_globin_genes\t{hba_genes}\nestimate\t{hba_genes}.02\n" + "message\tone gene lost on one chromosome\n") + + +def sample(work, name, hits, annotsv=None): + d = os.path.join(work, name, "S") + put(f"{d}/clinvar/S_clinvar_hits.vcf", VCF_HEAD + "".join(hits)) + if annotsv is not None: + put(f"{d}/annotsv/S_sv_class4_5_checked.tsv", ANNOTSV_COLS + annotsv) + return d + + +def carrier(d): + return collect_summary.sec_carrier(d, "S")[1] + + +def two(d, gene): + return [x for x in (carrier(d).get("two_hits") or {}).get("genes", []) if x["gene"] == gene] + + +def sf(d): + return collect_summary.acmg_sf_tier(d, "S", None) + + +def main(): + work = tempfile.mkdtemp(prefix="two_hits_") + try: + # 1. a MUTYH exon deletion plus a heterozygous hit: possibly biallelic + d = sample(work, "mutyh", [hit(SNV, "MUTYH", "0/1")], sv(DEL, "MUTYH")) + t = two(d, "MUTYH") + check("MUTYH deletion + het hit: one possibly biallelic line", + len(t) == 1 and t[0]["verdict"] == "possibly biallelic", t) + check("MUTYH deletion + het hit: the line names the deletion and the hit, phase unknown", + t and "DEL chr1:45335000-45337000" in t[0]["deletion"] and t[0]["variant"] == "chr1:45332000 C>T" + and "possibly in trans" in t[0]["why"], t) + s = sf(d) + rep = [x for x in s["clinvar_hits"] if x["gene"] == "MUTYH"] + check("MUTYH deletion + het hit: the SF tier reports the hit, counted with the deletion", + len(rep) == 1 and "possibly biallelic with DEL" in (rep[0].get("two_hit") or ""), rep) + check("MUTYH deletion + het hit: nothing left apart", + not [x for x in s["not_reported"] if x["gene"] == "MUTYH"], s["not_reported"]) + + # 2. read-backed phase on different copies: biallelic + d = sample(work, "trans", [hit(SNV, "MUTYH", "0|1", ps=45330000)], sv(DEL, "MUTYH", "1|0", ps=45330000)) + t = two(d, "MUTYH") + check("phased in trans: biallelic", len(t) == 1 and t[0]["verdict"] == "biallelic" + and "in trans" in t[0]["why"], t) + # a different phase set links nothing: possibly + d = sample(work, "otherps", [hit(SNV, "MUTYH", "0|1", ps=45330000)], sv(DEL, "MUTYH", "1|0", ps=45000000)) + t = two(d, "MUTYH") + check("phased in two different blocks: possibly biallelic", t and t[0]["verdict"] == "possibly biallelic", t) + + # 3. in cis: no line, and the SF tier keeps the one hit apart + d = sample(work, "cis", [hit(SNV, "MUTYH", "1|0", ps=45330000)], sv(DEL, "MUTYH", "1|0", ps=45330000)) + check("phased in cis: no two-hit line", two(d, "MUTYH") == [], two(d, "MUTYH")) + s = sf(d) + check("phased in cis: the SF tier lists the one hit apart, as before", + not [x for x in s["clinvar_hits"] if x["gene"] == "MUTYH"] + and [x for x in s["not_reported"] if x["gene"] == "MUTYH"], s) + + # 4. TAR: the 1q21.1 deletion of step 35 plus the RBM8A low-expression allele + tar_row = ("1q21.1_TAR_syndrome_proximal_BP2-BP3\t1q21.1_TAR_syndrome_proximal_BP2-BP3_t1\tchr1\t145686999\t" + "146048495\tcalled\t1\t2\tdeletion\t1.04\t.\n") + d = sample(work, "tar", [hit(TAR_ALLELE, "RBM8A", "1/1")]) + step35(d, intervals=tar_row) + c = carrier(d) + t = [x for x in (c.get("two_hits") or {}).get("genes", []) if x["gene"] == "RBM8A"] + check("TAR: the 1q21.1 deletion + the RBM8A allele give one RBM8A line", + len(t) == 1 and t[0]["verdict"] == "biallelic", c.get("two_hits")) + line = [x for x in render_report.carrier_lines(c) if "RBM8A" in x and "both copies" in x] + check("TAR: the line names the condition, the deletion and the allele", + len(line) == 1 and "Thrombocytopenia-absent radius syndrome" in line[0] + and "1q21.1_TAR_syndrome" in line[0] and "rs139428292" in line[0], line) + check("TAR: the allele alone stays set aside, never a carrier finding", + [x for x in c["set_aside"] if x["gene"] == "RBM8A"] and not + [x for x in c["carriers"] if x["gene"] == "RBM8A"], c) + # the allele without the deletion: no line + d = sample(work, "tar_alone", [hit(TAR_ALLELE, "RBM8A", "0/1")]) + step35(d, intervals=tar_row.replace("\t1\t2\tdeletion\t1.04", "\t2\t2\tno change\t2.01")) + check("TAR allele without a deletion: no line", two(d, "RBM8A") == [], two(d, "RBM8A")) + + # 5. a hit called 1/1 inside the deletion: hemizygous + d = sample(work, "hemi", [hit(INSIDE, "MUTYH", "1/1")], sv(DEL, "MUTYH")) + c = carrier(d) + item = [x for x in c["carriers"] if x["gene"] == "MUTYH"] + check("1/1 inside a deletion: the carrier list says hemizygous", + len(item) == 1 and item[0]["genotype"] == "hemizygous" and "inside DEL" in item[0].get("hemizygous", ""), + item) + t = two(d, "MUTYH") + check("1/1 inside a deletion: biallelic, the variant on the one copy left", + t and t[0]["verdict"] == "biallelic" and t[0]["genotype"] == "hemizygous", t) + rep = [x for x in sf(d)["clinvar_hits"] if x["gene"] == "MUTYH"] + check("1/1 inside a deletion: the SF tier says hemizygous", rep and rep[0]["genotype"] == "hemizygous", rep) + d = sample(work, "hom_out", [hit(SNV, "MUTYH", "1/1")], sv(DEL, "MUTYH")) + item = [x for x in carrier(d)["carriers"] if x["gene"] == "MUTYH"] + check("1/1 outside the deletion stays hom", item and item[0]["genotype"] == "hom", item) + # AnnotSV's SV_start is the base before the deleted bases: a 1/1 hit there is not inside + d = sample(work, "hom_edge", [hit(("chr1", DEL[1], "G", "A"), "MUTYH", "1/1")], sv(DEL, "MUTYH")) + item = [x for x in carrier(d)["carriers"] if x["gene"] == "MUTYH"] + check("1/1 at AnnotSV's SV_start (the padding base) stays hom", item and item[0]["genotype"] == "hom", item) + d = sample(work, "hom_first", [hit(("chr1", DEL[1] + 1, "G", "A"), "MUTYH", "1/1")], sv(DEL, "MUTYH")) + item = [x for x in carrier(d)["carriers"] if x["gene"] == "MUTYH"] + check("1/1 at the first deleted base is hemizygous", item and item[0]["genotype"] == "hemizygous", item) + + # 6. SMN1 and HBA + d = sample(work, "smn", [hit(("chr5", 70951950, "G", "A"), "SMN1", "0/1")]) + step35(d, smn1="1") + t = two(d, "SMN1") + check("one SMN1 copy + an SMN1 hit: a line", t and "one SMN1 copy" in t[0]["deletion"], t) + d = sample(work, "smn2", [hit(("chr5", 70951950, "G", "A"), "SMN1", "0/1")]) + step35(d, smn1="2") + check("two SMN1 copies + an SMN1 hit: no line", two(d, "SMN1") == [], two(d, "SMN1")) + d = sample(work, "hba", [hit(("chr16", 173600, "C", "T"), "HBA2", "0/1")]) + step35(d, hba_genes="3") + t = two(d, "HBA2") + check("three alpha genes + an HBA2 hit: possibly biallelic", + t and t[0]["verdict"] == "possibly biallelic" and "3 of 4 alpha-globin genes" in t[0]["deletion"], t) + d = sample(work, "hba4", [hit(("chr16", 173600, "C", "T"), "HBA2", "0/1")]) + step35(d, hba_genes="4") + check("four alpha genes + an HBA2 hit: no line", two(d, "HBA2") == [], two(d, "HBA2")) + + # 7. no deletion caller: no line, and the gap is named + d = sample(work, "none", [hit(SNV, "MUTYH", "0/1")]) + t = carrier(d).get("two_hits") or {"genes": None, "deletions_not_read": []} + check("no deletion caller: no line", t["genes"] == [], t) + check("no deletion caller: the callers not read are named", + any("AnnotSV" in x for x in t["deletions_not_read"]) + and any("step 35" in x for x in t["deletions_not_read"]), t) + check("no deletion caller: the report line says so", + any("not read: AnnotSV" in x for x in render_report.carrier_lines(carrier(d))), + render_report.carrier_lines(carrier(d))) + + # 8. both reports + d = sample(work, "report", [hit(INSIDE, "MUTYH", "1/1"), hit(("chr1", 45330500, "A", "G"), "MUTYH", "0/1")], + sv(DEL, "MUTYH")) + # step 23's clinical VCF, with no record of its own: the SF card is drawn from it + put(f"{d}/clinical/S_clinical.vcf.gz", VCF_HEAD, gz=True) + summ = collect_summary.collect("S", d) + txt, html = render_report.text_report(summ), render_report.html_report(summ) + check("text report: the two-hit line", "Biallelic: MUTYH" in txt, txt[-3000:]) + check("HTML report: the two-hit line", "Biallelic: MUTYH" in html) + check("text report: the hemizygous genotype", "45336000 G>A" not in txt and "hemizygous (called homozygous " + "inside DEL" in txt, txt[-3000:]) + check("text report: the SF row carries the two-hit detail", "biallelic with DEL chr1:45335000-45337000" in txt, + txt[-3000:]) + check("HTML report: the SF row carries the two-hit detail", "biallelic with DEL chr1:45335000-45337000" in html) + + # 9. the couple screen: hypomorph x null is a shared gene, hypomorph x hypomorph is not + genes = couple_report.read_genes(os.path.join(REPO, "assets", "carrier_genes.tsv")) + tar = sample(work, "couple_tar", [hit(TAR_ALLELE, "RBM8A", "0/1")]) + fa, _ = couple_report.clinvar_findings(f"{tar}/clinvar/S_clinvar_hits.vcf", genes, "female", {}, []) + null = couple_report.finding("RBM8A", "deletion chr1:145700000-146000000 (300.0 kb) het, 6 exons of the " + "gene (Delly)", 1) + sexes = {"A": "female", "B": "male"} + rows = couple_report.screen_couple("A", "B", sexes, {"A": fa, "B": [null]}, genes) + r = [x for x in rows if x["gene"] == "RBM8A"] + check("couple: RBM8A allele in A + RBM8A deletion in B is a shared gene, 1 in 4", + len(r) == 1 and r[0]["result"] == "shared" and r[0]["offspring_risk"] == "1 in 4 per pregnancy" + and "rs139428292" in r[0]["findings_a"] and "hypomorph" in r[0]["explanation"], rows) + rows = couple_report.screen_couple("A", "B", sexes, {"A": fa, "B": list(fa)}, genes) + check("couple: the RBM8A allele in both partners is set aside, no row", + not [x for x in rows if x["gene"] == "RBM8A"], rows) + rows = couple_report.screen_couple("A", "B", sexes, {"A": fa, "B": []}, genes) + check("couple: the RBM8A allele in one partner alone is set aside, no row", + not [x for x in rows if x["gene"] == "RBM8A"], rows) + ped = put(f"{work}/couple_run/pedigree.tsv", "sample\tsex\tfather\tmother\tpartner\nA\tfemale\t.\t.\tB\n" + "B\tmale\t.\t.\tA\n") + empty = sample(work, "couple_empty", []) + r = subprocess.run([sys.executable, os.path.join(REPO, "bin", "couple_report.py"), "screen", "--pedigree", ped, + "--genes", os.path.join(REPO, "assets", "carrier_genes.tsv"), + "--clinvar", f"A={tar}/clinvar/S_clinvar_hits.vcf", + "--clinvar", f"B={empty}/clinvar/S_clinvar_hits.vcf", + "--out-dir", f"{work}/couple_run/out"], capture_output=True, text=True) + md = f"{work}/couple_run/out/couple_A_B.md" + md_text = open(md).read() if os.path.isfile(md) else "" + check("couple run: the RBM8A allele counts as no finding of the listed genes (0 / 1) and is a gap line", + r.returncode == 0 and "| A | 0 / 1 |" in md_text and "Set aside for A" in md_text, + (r.returncode, r.stderr[-500:], md_text[-1500:])) + + # 10. edges + # MUTYH ends at chr1:45340893: a deletion whose padding base is that one removes none of it + d = sample(work, "pad_gene", [hit(SNV, "MUTYH", "0/1")], sv(("chr1", 45340893, 45342000), "MUTYH")) + check("a gene AnnotSV names that ends at the padding base: no line", two(d, "MUTYH") == [], two(d, "MUTYH")) + d = sample(work, "hba0", [hit(("chr16", 173600, "C", "T"), "HBA2", "0/1")]) + step35(d, hba_genes="0") + check("no alpha-globin gene left + an HBA2 hit: no line", two(d, "HBA2") == [], two(d, "HBA2")) + d = sample(work, "cis_inside", [hit(INSIDE, "MUTYH", "1|0", ps=45330000)], + sv(DEL, "MUTYH", "1|0", ps=45330000)) + check("a het hit inside a deletion, phased in cis: no line", two(d, "MUTYH") == [], two(d, "MUTYH")) + # a HIGH-impact record (no ClinVar hit) phased in trans with the deletion + csq = ("##fileformat=VCFv4.2\n##INFO=\n" + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tS\n" + f"chr1\t{SNV[1]}\t.\tC\tT\t50\tPASS\tCSQ=T|stop_gained|HIGH|MUTYH|0.0001\tGT:GQ:DP:AD:PS\t" + "0|1:50:30:15,15:45330000\n") + d = sample(work, "high_trans", [], sv(DEL, "MUTYH", "1|0", ps=45330000)) + clin = put(f"{d}/clinical/S_clinical.vcf", csq) + rep = [x for x in collect_summary.acmg_sf_tier(d, "S", clin)["high_impact"] if x["gene"] == "MUTYH"] + check("a HIGH-impact hit phased in trans with the deletion: biallelic in the SF tier", + len(rep) == 1 and (rep[0].get("two_hit") or "").startswith("biallelic with DEL"), rep) + # one hit in cis, one in trans: only the trans one carries the note + cis_hit = ("chr1", 45330500, "A", "G") + d = sample(work, "cis_trans", [hit(cis_hit, "MUTYH", "1|0", ps=45330000), + hit(SNV, "MUTYH", "0|1", ps=45330000)], sv(DEL, "MUTYH", "1|0", ps=45330000)) + rows = {x["variant"]: x.get("two_hit") for x in sf(d)["clinvar_hits"] if x["gene"] == "MUTYH"} + check("a cis and a trans hit in one gene: only the trans hit carries the two-hit note", + (rows.get("chr1:45332000 C>T") or "").startswith("biallelic with DEL") + and "chr1:45330500 A>G" in rows and rows["chr1:45330500 A>G"] is None, rows) + # an unreadable carrier list: the SF card says what was not read + d = sample(work, "no_list", [hit(SNV, "MUTYH", "0/1")], sv(DEL, "MUTYH")) + real = collect_summary.carrier_genes + + def broken(): + raise collect_summary.Unreadable("carrier_genes.tsv row 1: broken") + collect_summary.carrier_genes = broken + try: + s = sf(d) + finally: + collect_summary.carrier_genes = real + check("an unreadable carrier list: a note on the SF card", + "carrier_genes.tsv" in (s.get("two_hit_note") or "") + and any("carrier_genes.tsv" in x for x in render_report.sf_apart_lines(s)), s.get("two_hit_note")) + finally: + shutil.rmtree(work, ignore_errors=True) + print(f"\n{len(FAILS)} failed" if FAILS else "\nall passed") + return 1 if FAILS else 0 + + +if __name__ == "__main__": + sys.exit(main())