fix(cn): wire the copy-number intervals into the couple screen and the report, and close the review follow-ups

bin/cn_loci.tsv

@@ -148,6 +148,8 @@ Xq28_MECP2	interval	chrX	Xq28_MECP2_t1:154008528-154110279:Xq28 region (includes
 GJB6_D13S1830	interval	chr13	GJB6_D13S1830_t1:20223037-20531806:del(GJB6-D13S1830), ClinVar 5546	GJB6_D13S1830_L:20073037-20173037:flank left of the region;GJB6_D13S1830_R:20581806-20681806:flank right of the region	.	pooled	flanks	2	0	.	.	.	8	1.5	round:0.3	.	.	bedcov	cn-depth	cnvpytor:cnvs	11807148
 CFTR_dele2_3	interval	chr7	CFTR_dele2_3_t1:117498312-117519392:CFTRdele2,3 (21 kb), ClinVar 66105	CFTR_dele2_3_L:117348312-117448312:flank left of the region;CFTR_dele2_3_R:117569392-117669392:flank right of the region	.	pooled	flanks	2	0	.	.	.	8	1.5	round:0.3	.	.	bedcov	cn-depth	cnvpytor:cnvs	.
 HBB_exon3_del	interval	chr11	HBB_exon3_del_t1:5225254-5225875:HBB deletion of exon 3 (621 bp), ClinVar 3769040	HBB_exon3_del_L:5075254-5175254:flank left of the region;HBB_exon3_del_R:5275875-5375875:flank right of the region	.	pooled	flanks	2	0	.	.	.	8	1.5	round:0.3	.	.	bedcov	cn-depth	cnvpytor:cnvs	.
+HBB_619bp_del	interval	chr11	HBB_619bp_del_t1:5225388-5226007:HBB 619 bp deletion (South Asian), HbVar 979	HBB_619bp_del_L:5075388-5175388:flank left of the region;HBB_619bp_del_R:5276007-5376007:flank right of the region	.	pooled	flanks	2	0	.	.	.	8	1.5	round:0.3	.	.	bedcov	cn-depth	cnvpytor:cnvs	287080
+HBD_HBB_Sicilian_del	interval	chr11	HBD_HBB_Sicilian_del_t1:5219877-5233280:Sicilian (delta-beta)0 deletion (13.4 kb, HBD to HBB 3'), HbVar 1035	HBD_HBB_Sicilian_del_L:5069877-5169877:flank left of the region;HBD_HBB_Sicilian_del_R:5283280-5383280:flank right of the region	.	pooled	flanks	2	0	.	.	.	8	1.5	round:0.3	.	.	bedcov	cn-depth	cnvpytor:cnvs	2307466
 AZFa	interval	chrY	AZFa_t1:12537649-12860852:AZFa (USP9Y and DDX3Y);AZFa_t2:12904107-12920478:DDX3Y	AZFa_c1:2935200-2982579:ZFY, single-copy control on Yp;AZFa_c2:14522572-14845654:NLGN4Y, single-copy control on Yq	.	pooled	pooled	2	0	.	.	.	4	.	round:0.3	.	.	bedcov	cn-depth	cnvpytor:cnvs	37674303
 AZFb	interval	chrY	AZFb_t1:19702297-19745148:AZFb (KDM5D and EIF1AY);AZFb_t2:20575693-20593154:EIF1AY	AZFb_c1:2935200-2982579:ZFY, single-copy control on Yp;AZFb_c2:14522572-14845654:NLGN4Y, single-copy control on Yq	.	pooled	pooled	2	0	.	.	.	4	.	round:0.3	.	.	bedcov	cn-depth	cnvpytor:cnvs	37674303
 AZFc_DAZ	interval	chrY	AZFc_DAZ_t1:23129354-23199010:AZFc (DAZ1 to DAZ4 added up, gr/gr and b2/b3 leave two);AZFc_DAZ_t2:23219446-23291356:DAZ2;AZFc_DAZ_t3:24763068-24813492:DAZ3;AZFc_DAZ_t4:24833842-24907040:DAZ4	AZFc_DAZ_c1:2935200-2982579:ZFY, single-copy control on Yp;AZFc_DAZ_c2:14522572-14845654:NLGN4Y, single-copy control on Yq	.	sum	pooled	2	0	.	.	.	4	.	round:0.3	.	.	bedcov	cn-depth	cnvpytor:cnvs	37674303

bin/collect_summary.py

@@ -844,7 +844,10 @@ def sec_parascopy(d, s):
     the other loci."""
     tables = parascopy_tables(d, s)
     hp = os.path.join(d, "paralogs", f"{s}_hba_copy_number.tsv")
-    p = tables.get("smn") or next(iter(tables.values()), None) or (hp if os.path.isfile(hp) else None)
+    ip = os.path.join(d, "paralogs", f"{s}_cn_intervals.tsv")
+    cp = os.path.join(d, "paralogs", f"{s}_cyp21a2_psv.tsv")
+    p = (tables.get("smn") or next(iter(tables.values()), None) or (hp if os.path.isfile(hp) else None)
+         or (ip if os.path.isfile(ip) else None) or (cp if os.path.isfile(cp) else None))
     if not p:
         return None, {}
     if "smn" in tables:
@@ -867,8 +870,6 @@ def sec_parascopy(d, s):
         reading = f"{NO_SMN_TABLE}: no carrier reading"
     others = {n: t for n, t in tables.items() if n != "smn"}
     catalog = {r["table"]: r for r in asset_rows("paralog_loci.tsv")} if others else {}
-    ip = os.path.join(d, "paralogs", f"{s}_cn_intervals.tsv")
-    cp = os.path.join(d, "paralogs", f"{s}_cyp21a2_psv.tsv")
     return p, {"smn1": smn1, "smn2": {"cn": smn2["cn"], "text": smn2["text"]},
                "agcn": smn1["agcn"], "agcn_qual": smn1["agcn_qual"], "agcn_filter": smn1["agcn_filter"],
                "rows": len(rows),

bin/couple_report.py

@@ -20,7 +20,8 @@
       and FXN repeats (ExpansionHunter, step 09). With --cn-intervals (step
       35's <id>_cn_intervals.tsv) the report lists, for each partner, every
       carrier deletion allele of that table (del(GJB6-D13S1830), CFTRdele2,3,
-      the HBB exon 3 deletion) as present, absent or not assessable, and
+      the HBB exon 3 deletion, the HBB 619 bp deletion, the Sicilian
+      (delta-beta)0 deletion) as present, absent or not assessable, and
       a male partner's AZF deletions the same way; with --cnvpytor (step 18's
       <id>_cnvs.vcf.gz) a CNVpytor call that disagrees is printed beside it.
       With --cyp21a2 (step 35's <id>_cyp21a2_psv.tsv, short reads) a
@@ -448,12 +449,16 @@ def cyp21a2_finding(path, genes):
 
 # The rows of bin/cn_loci.tsv the couple report reads, with their names here.
 CN_ALLELES = (("GJB6_D13S1830", "del(GJB6-D13S1830)"), ("CFTR_dele2_3", "CFTRdele2,3"),
-              ("HBB_exon3_del", "HBB exon 3 deletion"))
+              ("HBB_exon3_del", "HBB exon 3 deletion"), ("HBB_619bp_del", "HBB 619 bp deletion"),
+              ("HBD_HBB_Sicilian_del", "Sicilian (delta-beta)0 deletion"))
 # Named as deletions, so present reads as "the deletion is there" (as for CN_ALLELES).
 CN_AZF = (("AZFa", "AZFa deletion"), ("AZFb", "AZFb deletion"),
           ("AZFc_DAZ", "AZFc deletion (DAZ copies; gr/gr and b2/b3 leave two of four)"))
 # A CNVpytor deletion counts for a row when it covers half of the row's span or more.
 CNV_COVER = 0.5
+# Step 18 calls on 1 kb bins, so a deletion shorter than two bins (the 619 bp
+# and 621 bp HBB rows) is one it cannot see: its silence there is no evidence.
+CNV_MIN_SPAN = 2000
 
 
 def read_cn_intervals(path):
@@ -502,7 +507,7 @@ def cnv_note(row, dels):
     s, e = int(row["start"]), int(row["end"])
     hit = any(c == row["chrom"] and min(e, de) - max(s, ds) >= CNV_COVER * (e - s) for c, ds, de in dels)
     state, _ = cn_state(row)
-    if state == "present" and not hit:
+    if state == "present" and not hit and e - s >= CNV_MIN_SPAN:
         return "CNVpytor (step 18) calls no deletion here"
     if state == "absent" and hit:
         return "CNVpytor (step 18) calls a deletion here"
@@ -672,7 +677,10 @@ def section(title, keep, none):
     if cn:
         out.extend(["## Carrier deletion alleles and AZF regions (read depth, step 35)", "",
                     "Each allele is present, absent or not assessable from the read depth over it. A present "
-                    "deletion needs a clinical test. The AZF regions are read on a male partner only.", ""])
+                    "deletion needs a clinical test. The AZF regions are read on a male partner only. The three "
+                    "beta-globin rows overlap: the Sicilian deletion covers both HBB rows, and the 619 bp and "
+                    "exon 3 rows share 487 bp, so one deletion can read as present on two or three of them; "
+                    "only a clinical test tells which allele it is.", ""])
         out.extend(cn)
         out.append("")
     out.extend(["## What was assessed", "", "| Partner | ClinVar hits (in listed genes / all) | Deletions | SMN1 | FMR1 and FXN repeats |",

docs/35-paralogs.md

@@ -76,6 +76,7 @@ All output is written to `${GENOME_DIR}/${SAMPLE}/paralogs/`.
 | `${SAMPLE}_smn_copy_number.tsv` | One row per region: `chrom`, `start`, `end`, `locus`, `agCN_filter`, `agCN`, `agCN_qual`, `psCN_filter`, `psCN`, `psCN_qual`, `homologous_regions`. A locus that was not assessed has one row over the catalog region (`chrom`, `start`, `end`), with `locus` (the model name), `agCN_filter` and `psCN_filter` (both `not assessed: <reason>`) and `.` in every other column |
 | `${SAMPLE}_<table>_copy_number.tsv` | The same, for every other catalog row (`gba1`, `pms2_pms2cl`, `c4` and the rest of the table above) |
 | `${SAMPLE}_hba_copy_number.tsv` | `metric`, `value`: `status` (`called` or `not assessed`), `alpha_globin_genes`, `estimate`, `hba2_depth`, `hba1_depth`, `flank_left_depth`, `flank_right_depth`, `mapq0_share` (the part of the depth over the two genes from reads of mapping quality 0), `message` |
+| `${SAMPLE}_cn_intervals.tsv` | One line per interval row of `bin/cn_loci.tsv` (ClinGen's recurrent CNV regions, the carrier deletion alleles, the Y chromosome AZF regions): `locus`, `name`, `chrom`, `start`, `end`, `status`, `copies`, `expected`, `change` (`deletion`, `duplication`, `no change`, or `.` when not called), `estimate`, `reason` ([the interval rows](cn-loci.md#the-interval-rows)). The couple screen (step 39) reads it for the carrier deletion alleles and AZF |
 | `${SAMPLE}_cyp21a2_psv.tsv` | `metric`, `value`: `status`, `cyp21a2_copies` (the least CYP21A2 copies over the low sites, else the median site's), `module_total` (CYP21A2 plus CYP21A1P copies), `expected_total` (4), `estimate`, `finding` (`none`, `deletion`, `conversion`, `pseudogene_deletion`), `conversion_sites`, one `psv_<site>` row per site (gene, pseudogene and other reads, copies, ratio to the other sites), `target_depth`, `control_depth`, `reason`, `reading` |
 | `${SAMPLE}_parascopy/<table>/` | Parascopy's own output for each catalog row (`smn/` for SMN1): `res.samples.bed.gz`, `res.paralog.bed.gz`, `psvs.vcf.gz` and the rest ([described upstream](https://github.com/tprodanov/parascopy/blob/main/docs/cn_output.md)) |
 

docs/39-couple.md

@@ -38,7 +38,7 @@ The table of what was assessed counts, for each partner, the deletions over a li
 
 Only the ClinVar screen is required. Each of the others is used when its file is there, and the report lists every one that was not.
 
-`couple_report.py screen` also takes step 35's interval table (`--cn-intervals ID=paralogs/<sample>_cn_intervals.tsv`, [the interval rows](cn-loci.md#the-interval-rows)). The report then lists, for each partner, del(GJB6-D13S1830), CFTRdele2,3 and the HBB exon 3 deletion as present, absent or not assessable from read depth, and a male partner's AZFa, AZFb and AZFc deletions (Y chromosome microdeletions, a cause of male infertility) the same way, where present means the deletion is there; a female partner's AZF line says not applicable. With `--cnvpytor ID=cnvpytor/<sample>_cnvs.vcf.gz` (step 18), a CNVpytor call that disagrees with the depth call (a deletion over half the allele's span or more where depth sees none, or none where depth sees one) is printed beside it. These lines are a list, not carrier findings: they do not enter the gene table. Step 39 and the Nextflow couple screen do not pass the two inputs yet; run the command by hand to get them.
+`couple_report.py screen` also takes step 35's interval table (`--cn-intervals ID=paralogs/<sample>_cn_intervals.tsv`, [the interval rows](cn-loci.md#the-interval-rows)). The report then lists, for each partner, del(GJB6-D13S1830), CFTRdele2,3, the HBB exon 3 deletion, the HBB 619 bp deletion and the Sicilian (delta-beta)0 deletion as present, absent or not assessable from read depth, and a male partner's AZFa, AZFb and AZFc deletions (Y chromosome microdeletions, a cause of male infertility) the same way, where present means the deletion is there; a female partner's AZF line says not applicable. The three beta-globin rows overlap (the Sicilian deletion covers both HBB rows, and the 619 bp and exon 3 rows share 487 bp), so one deletion can read as present on two or three of them; the report says so, and only a clinical test names the allele. With `--cnvpytor ID=cnvpytor/<sample>_cnvs.vcf.gz` (step 18), a CNVpytor call that disagrees with the depth call (a deletion over half the allele's span or more where depth sees none, or none where depth sees one) is printed beside it. A missing call is not printed for an allele shorter than 2 kb (the two HBB rows): step 18 counts in 1 kb bins and cannot see a deletion that short. These lines are a list, not carrier findings: they do not enter the gene table. Step 39 passes both when a partner's folder holds them (`paralogs/<sample>_cn_intervals.tsv` of step 35, `cnvpytor/<sample>_cnvs.vcf.gz` of step 18), and the Nextflow couple screen passes them when `--tools` has `parascopy` and `cnvpytor`.
 
 ### The carrier list