From fea0cb43e80fc731317f23b739f943e135e9b263 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sun, 8 Nov 2026 15:47:08 +0100 Subject: [PATCH 1/3] fix(cn): stop calling a CYP21A2 deletion when the BAM dropped its multi-mapping reads The CYP21A2 row added up the depth over CYP21A2 and CYP21A1P with no mapping quality check. An aligner that leaves out the reads fitting both copies (MAPQ 0) lowers the module total, which can read as the 30 kb chimeric deletion. The row now runs a second bedcov pass with -Q 1 and mapq0_min 0.02, as HBA does, in step 35 and the PARASCOPY module. Under that share it is not assessed (low_mapq0_share). The table gains two rows, mapq0_share and message; the reports print the message. HG002's 35x reads give a share of 0.099 over the two targets from BWA-MEM and from minimap2 -ax sr alike. Bead: pgp-5bja.21 --- bin/cn_loci.tsv | 6 ++- bin/copy_number.py | 20 ++++++++- docs/35-paralogs.md | 10 ++--- docs/cn-loci.md | 4 +- docs/lessons-learned.md | 7 +++ modules/local/parascopy/main.nf | 9 ++-- scripts/35-paralogs.sh | 12 +++-- tests/test_copy_number.py | 78 ++++++++++++++++++++++++++++++--- 8 files changed, 124 insertions(+), 22 deletions(-) diff --git a/bin/cn_loci.tsv b/bin/cn_loci.tsv index 49e60f23..aa1b5d4c 100644 --- a/bin/cn_loci.tsv +++ b/bin/cn_loci.tsv @@ -53,6 +53,10 @@ # In2G rs6467, I172N rs6475, V237E rs12530380 (the exon 6 cluster), V281L rs6471, Q318X rs7755898, # R356W rs7769409 (dbSNP positions, the CYP21A1P bases from the reference, checked 2026-11-04). # Coordinates from UCSC's ncbiRefSeq and genomicSuperDups tracks, checked 2026-11-04. +# MAPQ 0 guard as HBA's (0.02): HG002's 35x reads give a MAPQ 0 share of 0.099 over the two +# targets from BWA-MEM and from minimap2 -ax sr alike (measured 2026-11-08); an aligner that +# leaves those reads out gives about 0 and a lower module total (4.22 to 3.81 on HG002), a deletion +# on a genome whose depth runs lower. # psv:0.4, wider than HBA's 0.3: HG002's 35x reads give a module total of 4.22. locus kind chrom targets controls details target_depth control_depth scale offset bedcov mapq0_min mapq0_max min_control_depth max_control_ratio call interval readings input owner cross_check pmids HBA gene chr16 HBA2:172875-173710;HBA1:176679-177522 flank_left:100000-140000;flank_right:190000-230000 . sum flanks 2 0 .,-Q 1 0.02 . 8 1.5 round:0.3 . 0=no alpha-globin gene: the pattern of Hb Bart's hydrops fetalis (--/--); unconfirmed: this short-read count has been checked on one real genome only (HG002, four genes), so a count under four needs a clinical test before it means anything|1=one alpha-globin gene: the pattern of HbH disease (--/-a); unconfirmed: this short-read count has been checked on one real genome only (HG002, four genes), so a count under four needs a clinical test before it means anything|2=two alpha-globin genes: alpha-thalassaemia trait, either both genes lost on one chromosome (--/aa, as --SEA; a child with a partner who carries the same can have Hb Bart's) or one on each (-a/-a); depth cannot tell these apart; unconfirmed: this short-read count has been checked on one real genome only (HG002, four genes), so a count under four needs a clinical test before it means anything|3=three alpha-globin genes: one gene lost on one chromosome (-a/aa, as -a3.7 or -a4.2), a silent carrier; unconfirmed: this short-read count has been checked on one real genome only (HG002, four genes), so a count under four needs a clinical test before it means anything|4=four alpha-globin genes (aa/aa): no deletion seen|more={n} alpha-globin genes: more than four, a triplication (as aaa/aa) bedcov cn-depth paraphase:hba . @@ -61,7 +65,7 @@ RHD gene chr1 RHD_GENE:25272392-25330445:RHD, first to last base (ENST0000032866 HP gene chr16 HP_UNIT1:72056529-72057466:HP exons 3 and 4, the first copy of the Hp2 duplication (ENST00000355906);HP_UNIT2:72058253-72059190:HP exons 5 and 6, its second copy HP_FLANK_L:72005000-72045000:flank upstream of HP;HP_FLANK_R:72085000-72125000:flank downstream of HPR . sum flanks 2 0 .,-Q 1 . . 8 1.5 round:0.3 . 0=no copy of the duplicated HP exons: no HP gene found on either chromosome|1=one copy of the duplicated HP exons: Hp1 on one chromosome and no HP gene found on the other|2=two copies of the duplicated HP exons: Hp1-1, or rarely Hp2 with no HP gene on the other chromosome|3=three copies of the duplicated HP exons: Hp2-1|4=four copies of the duplicated HP exons: Hp2-2|more={n} copies of the duplicated HP exons: more than Hp2-2 carries, an uncommon allele with a further copy bedcov cn-depth none 26901066 LCE3BC interval chr1 LCE3BC_DEL:152600233-152614098:LCE3C to LCE3B, inside the LCE3C_LCE3B deletion LCE3_FLANK_L:152520000-152560000:flank on the LCE3D side, outside the deletion;LCE3_FLANK_R:152640000-152680000:flank on the LCE3A side, outside the deletion . sum flanks 2 0 .,-Q 1 . . 8 1.5 round:0.3 . 0=LCE3B and LCE3C deleted on both chromosomes (LCE3C_LCE3B-del on both), a common genotype that goes with a higher psoriasis risk|1=LCE3B and LCE3C deleted on one chromosome (LCE3C_LCE3B-del on one)|2=no LCE3B/LCE3C deletion seen|more={n} copies of the LCE3B/LCE3C stretch: more than two, a duplication bedcov cn-depth none 19169253 CCL3L1 gene chr17 CCL3L3:36194868-36196758:CCL3L3, the one copy of CCL3L1/CCL3L3 on the reference (ENST00000619989) CCL3L_FLANK_L:36112000-36148000:single-copy sequence between the CCL3 and CCL3L3 duplications;CCL3L_FLANK_R:36460000-36500000:single-copy sequence past the duplications . sum flanks 2 0 .,-Q 1 . . 8 1.5 round:0.3 . 0=no copy of CCL3L1 or CCL3L3|1=one copy of CCL3L1/CCL3L3 in all|more={n} copies of CCL3L1/CCL3L3 in all: the count differs widely between people, and no reading of it is given bedcov cn-depth none 15637236 -CYP21A2 gene chr6 CYP21A2_BODY:32038950-32041644:CYP21A2 from the end of the C4A model to its last base (NM_000500.9);CYP21A1P_BODY:32006220-32008909:the same stretch of CYP21A1P (NR_040090.1) RCCX_FLANK_L:31945000-31975000:single-copy sequence before the RCCX duplication;RCCX_FLANK_R:32070000-32100000:single-copy sequence in TNXB past the duplication PSV_In2G:32039080-32039081:32006352 AC>G c.293-13A/C>G (In2G, intron 2);PSV_I172N:32039425-32039426:32006689 T>A c.518T>A (p.I172N, exon 4);PSV_V237E:32039809-32039810:32007073 T>A c.713T>A (p.V237E, exon 6 cluster);PSV_V281L:32040109-32040110:32007373 G>T c.844G>T (p.V281L, exon 7);PSV_Q318X:32040420-32040421:32007685 C>T c.955C>T (p.Q318X, exon 8);PSV_R356W:32040534-32040535:32007799 C>T c.1069C>T (p.R356W, exon 8) sum flanks 2 0 . . . 8 1.5 psv:0.4 . none=no sign of the 30 kb chimeric deletion or of a pseudogene conversion at the six sites read; short reads see neither a change elsewhere in CYP21A2 nor two copies on one chromosome and none on the other|deletion=one RCCX module fewer and one CYP21A2 copy fewer at the differing sites: the pattern of the 30 kb CYP21A1P/CYP21A2 chimeric deletion, a common cause of carrier status for 21-hydroxylase deficiency (congenital adrenal hyperplasia); unconfirmed: a clinical test (MLPA or long-range PCR) confirms it|conversion=fewer CYP21A2 copies carry the gene's own base at {sites} than at the other sites: a possible conversion of one CYP21A2 copy to the pseudogene sequence there, a carrier pattern for 21-hydroxylase deficiency (congenital adrenal hyperplasia) when it sits on one chromosome; which chromosome carries it is not known; unconfirmed: a clinical test confirms it|pseudogene_deletion=one RCCX module fewer with CYP21A2 kept at every site read: a CYP21A1P (pseudogene) deletion, common and not a carrier finding bedcov cn-depth paraphase:rccx 34285390,27228465 +CYP21A2 gene chr6 CYP21A2_BODY:32038950-32041644:CYP21A2 from the end of the C4A model to its last base (NM_000500.9);CYP21A1P_BODY:32006220-32008909:the same stretch of CYP21A1P (NR_040090.1) RCCX_FLANK_L:31945000-31975000:single-copy sequence before the RCCX duplication;RCCX_FLANK_R:32070000-32100000:single-copy sequence in TNXB past the duplication PSV_In2G:32039080-32039081:32006352 AC>G c.293-13A/C>G (In2G, intron 2);PSV_I172N:32039425-32039426:32006689 T>A c.518T>A (p.I172N, exon 4);PSV_V237E:32039809-32039810:32007073 T>A c.713T>A (p.V237E, exon 6 cluster);PSV_V281L:32040109-32040110:32007373 G>T c.844G>T (p.V281L, exon 7);PSV_Q318X:32040420-32040421:32007685 C>T c.955C>T (p.Q318X, exon 8);PSV_R356W:32040534-32040535:32007799 C>T c.1069C>T (p.R356W, exon 8) sum flanks 2 0 .,-Q 1 0.02 . 8 1.5 psv:0.4 . none=no sign of the 30 kb chimeric deletion or of a pseudogene conversion at the six sites read; short reads see neither a change elsewhere in CYP21A2 nor two copies on one chromosome and none on the other|deletion=one RCCX module fewer and one CYP21A2 copy fewer at the differing sites: the pattern of the 30 kb CYP21A1P/CYP21A2 chimeric deletion, a common cause of carrier status for 21-hydroxylase deficiency (congenital adrenal hyperplasia); unconfirmed: a clinical test (MLPA or long-range PCR) confirms it|conversion=fewer CYP21A2 copies carry the gene's own base at {sites} than at the other sites: a possible conversion of one CYP21A2 copy to the pseudogene sequence there, a carrier pattern for 21-hydroxylase deficiency (congenital adrenal hyperplasia) when it sits on one chromosome; which chromosome carries it is not known; unconfirmed: a clinical test confirms it|pseudogene_deletion=one RCCX module fewer with CYP21A2 kept at every site read: a CYP21A1P (pseudogene) deletion, common and not a carrier finding bedcov cn-depth paraphase:rccx 34285390,27228465 # --- interval rows: written by scripts/ci/build-cn-intervals.py; do not edit by hand --- 1q21.1_TAR_syndrome_proximal_BP2-BP3 interval chr1 1q21.1_TAR_syndrome_proximal_BP2-BP3_t1:145686999-145809382:1q21.1 recurrent (TAR syndrome) region (proximal, BP2-BP3) (includes RBM8A);1q21.1_TAR_syndrome_proximal_BP2-BP3_t2:145823828-145923649;1q21.1_TAR_syndrome_proximal_BP2-BP3_t3:145927476-146048495 1q21.1_TAR_syndrome_proximal_BP2-BP3_L:143034587-143134587:flank left of the region;1q21.1_TAR_syndrome_proximal_BP2-BP3_R:149926380-150026380:flank right of the region . pooled flanks 2 0 . . . 8 1.5 round:0.3 . . bedcov cn-depth cnvpytor:cnvs . 1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4 interval chr1 1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4_t1:145687000-145809382:1q21.1 recurrent (TAR syndrome) region (proximal and distal, BP1-BP4) (includes RBM8A and GJA5);1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4_t2:145823828-145923649;1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4_t3:145927476-146048496;1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4_t4:147105903-147473859;1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4_t5:147473859-147841815;1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4_t6:147842996-147917509 1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4_L:143034587-143134587:flank left of the region;1q21.1_TAR_syndrome_proximal_and_distal_BP1-BP4_R:149926380-150026380:flank right of the region . pooled flanks 2 0 . . . 8 1.5 round:0.3 . . bedcov cn-depth cnvpytor:cnvs . diff --git a/bin/copy_number.py b/bin/copy_number.py index 0eeb2e6c..92bd31f1 100755 --- a/bin/copy_number.py +++ b/bin/copy_number.py @@ -580,6 +580,22 @@ def psv_assess(row, b, counts): return b +def psv_message(row, b): + """'.' for a called psv row, else 'not assessed: ' and why. A low MAPQ 0 + share is spelled out: the reports print this message, and the reason + code alone would not say that the BAM, not the genome, stopped the call.""" + if b["status"] == "called": + return "." + if b["reason"] == "low_mapq0_share": + targets = " and ".join(n.removesuffix("_BODY") for n, *_ in row["_targets"]) + return (f"not assessed: the BAM has no multi-mapping reads here (reads of mapping quality 0 make up " + f"{100 * b['mapq0_share']:.1f}% of the depth over {targets}, " + f"under {100 * row['_mapq0_min']:g}%), so the aligner left out the reads that fit both copies (as " + "novoalign does by default) and the module total would come out low and read as a deletion; " + "realign with minimap2 (step 02) or BWA-MEM") + return f"not assessed: {b['reason'] or '.'}" + + def psv_table(row, b): """([(metric, value)], one line) of a psv row's step table.""" f = (lambda x: "." if x is None else f"{x:.2f}") @@ -591,7 +607,9 @@ def psv_table(row, b): rows.append((f"psv_{name}", f"{s['gene_reads']} gene, {s['pseudo_reads']} pseudogene, {s['other_reads']} other" f" reads; copies {d(s['copies'])} ({f(s['estimate'])}); ratio {f(s['ratio'])}")) rows += [("target_depth", f(b["target_depth"])), ("control_depth", f(b["control_depth"])), - ("reason", b["reason"] or "."), ("reading", b["reading"] or ".")] + ("reason", b["reason"] or "."), ("reading", b["reading"] or "."), + ("mapq0_share", "." if b.get("mapq0_share") is None else f"{b['mapq0_share']:.3f}"), + ("message", psv_message(row, b))] line = (f"{row['locus']}: {b['status']} ({b['reason'] or b.get('finding')}; {d(b['copies'])} gene copies, " f"module total {d(b.get('module_total'))} of {d(b['expected'])})") return rows, line diff --git a/docs/35-paralogs.md b/docs/35-paralogs.md index 51a450e5..79ff7035 100644 --- a/docs/35-paralogs.md +++ b/docs/35-paralogs.md @@ -64,7 +64,7 @@ With Nextflow: `--tools ...,parascopy --parascopy_data ${GENOME_DIR}/reference/p 3. Writes one row per region of each locus profile to `${SAMPLE}__copy_number.tsv` (`${SAMPLE}_smn_copy_number.tsv` for the SMN1 row). 4. `samtools bedcov` (in the Parascopy image) sums the depth over HBA2, HBA1 and two flanks, once over every read and once over reads of mapping quality 1 or more (`-Q 1`), and `bin/copy_number.py call --locus HBA` (the copy-number engine, on the HBA row of `bin/cn_loci.tsv`; see [Copy number from depth](cn-loci.md)) writes `${SAMPLE}_hba_copy_number.tsv`, laid out by `bin/hba_copy_number.py`. A failure here prints a warning and writes a table that says not assessed; the SMN result stands. 5. The gene-dosage rows (HP, LCE3BC, CCL3L1) the same way, into `${SAMPLE}_cn_dosage.tsv` ([the gene-dosage rows](cn-loci.md#the-gene-dosage-rows)). -6. `samtools bedcov` over the CYP21A2 row's segments (one pass, every read), `samtools mpileup` without a reference at its PSV sites, then `bin/copy_number.py call --locus CYP21A2 --pileup ...` writes `${SAMPLE}_cyp21a2_psv.tsv`. A failure prints a warning and writes a table that says not assessed (`check_failed`). +6. `samtools bedcov` over the CYP21A2 row's segments (two passes, all reads and mapping quality 1 or more, as HBA), `samtools mpileup` without a reference at its PSV sites, then `bin/copy_number.py call --locus CYP21A2 --bedcov-mapq1 ... --pileup ...` writes `${SAMPLE}_cyp21a2_psv.tsv`. A failure prints a warning and writes a table that says not assessed (`check_failed`). 7. The interval rows into `${SAMPLE}_cn_intervals.tsv` ([the interval rows](cn-loci.md#the-interval-rows)). 8. SMN1 g.27134T>G beside the SMN1 psCN (before the HBA count): `samtools mpileup` without a reference at its SMN1 and SMN2 positions, then `bin/traits_report.py --marker SMN1_27134` with the SMN table writes `${SAMPLE}_smn1_marker.tsv` ([the SMN1 silent-carrier marker](#the-smn1-silent-carrier-marker)). A failure prints a warning and writes a table that says not assessed. @@ -79,7 +79,7 @@ All output is written to `${GENOME_DIR}/${SAMPLE}/paralogs/`. | `${SAMPLE}_smn1_marker.tsv` | `metric`, `value`: `marker`, `rsid`, `status` (`present`, `not seen`, `uncertain`, `too little depth`, `not assessed`), `alt_reads`, `depth`, `fraction`, `copies` (the fraction times agCN), `smn1_pscn` (the SMN1 psCN at exon 7 with its quality and filter), `agcn`, `message`. The step's log prints the message after the SMN table | | `${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)). Both reports show it, and 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_` row per site (gene, pseudogene and other reads, copies, ratio to the other sites), `target_depth`, `control_depth`, `reason`, `reading` | +| `${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_` row per site (gene, pseudogene and other reads, copies, ratio to the other sites), `target_depth`, `control_depth`, `reason`, `reading`, `mapq0_share` (the part of the depth over the two targets from reads of mapping quality 0), `message` (`.` when called, else `not assessed:` and why; the reports print it) | | `${SAMPLE}_parascopy/
/` | 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)) | ## The SMN1 silent-carrier marker @@ -228,7 +228,7 @@ The `finding`: | `conversion` | the module total is 4 or more and a site is low: its estimate is more than 0.6 copies under the larger of the median site's count and 2 | a possible conversion at that site; `conversion_sites` names it | | `pseudogene_deletion` | the module total is under 4 and no site is low | a CYP21A1P deletion, common and not a carrier finding | -The ratio of a site is its estimate over the median of the sites, so one copy converted out of two reads 0.5. A site needs 30 reads (CYP21A2 and CYP21A1P bases together); with no site there the row is not assessed (`psv_low_depth`). Which chromosome carries a change is never said, and a change away from the six sites is not seen. The row has no mapping quality 0 check yet, as HBA has: a BAM from an aligner that leaves out the reads fitting both genes (novoalign by default) would give a low module total, which can read as a deletion. Use minimap2 or BWA-MEM. A PacBio HiFi BAM gets CYP21A2 from Paraphase (step 35b) instead, which phases the copies. +The ratio of a site is its estimate over the median of the sites, so one copy converted out of two reads 0.5. A site needs 30 reads (CYP21A2 and CYP21A1P bases together); with no site there the row is not assessed (`psv_low_depth`). Which chromosome carries a change is never said, and a change away from the six sites is not seen. The row has a mapping quality 0 check, as HBA has. A read that fits CYP21A2 and CYP21A1P equally gets mapping quality 0, and HG002's 35x reads give about 10% of the depth over the two targets from such reads, with BWA-MEM and with minimap2 alike. An aligner that leaves those reads out (novoalign by default) gives a lower module total (on HG002's reads the estimate falls from 4.22 to 3.81), which can read as a deletion. So when they make up under 2% of that depth, the row is not assessed (`low_mapq0_share`), and the `message` row and both reports say the BAM has no multi-mapping reads there. Realign with minimap2 or BWA-MEM. A PacBio HiFi BAM gets CYP21A2 from Paraphase (step 35b) instead, which phases the copies. **A research estimate, not a diagnostic test.** A deletion or a conversion needs a clinical test (MLPA, long-range PCR or targeted sequencing). @@ -240,8 +240,8 @@ The ratio of a site is its estimate over the median of the sites, so one copy co - `tests/test_copy_number.py` checks the HBA row and that these depth tables give, through the engine, the same tables as the code before it. `tests/test_hba_copy_number.py` gives `bin/hba_copy_number.py` made-up depth tables: four genes at even depth, three with HBA2 at half depth, two, a triplication, an estimate between two whole numbers, low flanks, flanks that disagree, and no read of mapping quality 0 over the genes (the depths of the novoalign BAM, and a clean three), which must be not assessed. - The e2e case `paralog-checks-and-doc-fixes-1-hba` writes synthetic 30x reads over chr16:95,001-235,000 for two samples, one with every second read over HBA2 left out (a planted one-gene loss), runs step 35 on each and checks four and three genes. A synthetic read whose bases occur more than once in that stretch gets mapping quality 0, as an aligner gives it. The HBA block of `modules/local/parascopy/main.nf`, cut from the file, then runs on the same BAMs in the same image and must write the same table. Last, the block runs on HG002's real Illumina reads at the cluster, three ways: GIAB's 60x novoalign BAM sampled to about 30x must be not assessed, naming the missing mapping quality 0 reads; Google's 35x NovaSeq BWA-MEM BAM must count four; and the same reads realigned with `minimap2 -a -x sr` to chr16:1-300,000 must count four. The truth is four, from Paraphase on HG002's HiFi reads. The run's notes record each estimate, the four depths and the share of mapping quality 0. - The e2e case `paralog-rows-1` runs the whole catalog on the fixture's slice plus GIAB's HG002 reads over every catalog locus on the fixture's chromosomes, sampled to about 30x, timed (it fails past 30 minutes): every row's table must hold an integer agCN or say not assessed with the reason, and a header-only table fails that check. It then runs the step with a catalog that has no SMN1 row and checks that both reports still show the HBA count and the other loci. `tests/test_collect_summary.py` holds the same case on made-up tables, with the loci read in the catalog's order. -- `tests/test_copy_number.py` section 10 runs the CYP21A2 row on made-up bedcov and mpileup tables (`tests/fixtures/cn/CYP21A2_*`): unmodified, 2 copies and no finding; a planted chimeric deletion, 1 copy with the module total 3; a planted In2G conversion on one copy, a ratio of 0.5 at In2G and a possible conversion there. The same depths with every site even give no finding. The CYP21A2 block of step 35, cut out, writes those tables and a not-assessed one when mpileup fails, and the couple screen reads them. -- The e2e case `cyp21a2-short-reads-1-hg002` runs the CYP21A2 block of step 35 on HG002's Illumina reads over the row, unmodified and with a planted In2G conversion and a planted chimeric deletion made from the same reads, and the block of the PARASCOPY module on the unmodified BAM, which must write the same table. +- `tests/test_copy_number.py` section 10 runs the CYP21A2 row on made-up bedcov and mpileup tables (`tests/fixtures/cn/CYP21A2_*`): unmodified, 2 copies and no finding; a planted chimeric deletion, 1 copy with the module total 3; a planted In2G conversion on one copy, a ratio of 0.5 at In2G and a possible conversion there. The same depths with every site even give no finding. The deletion's tables with no MAPQ 0 read (a share under 0.02) read not assessed (`low_mapq0_share`), and with a share of 0.021 or HG002's 0.1 still give the deletion. The CYP21A2 block of step 35, cut out, writes those tables, a not-assessed one when mpileup fails and a `low_mapq0_share` one when its MAPQ 1 pass equals the first, and the couple screen and the report read them. +- The e2e case `cyp21a2-short-reads-1-hg002` runs the CYP21A2 block of step 35 on HG002's Illumina reads over the row, unmodified and with a planted In2G conversion and a planted chimeric deletion made from the same reads, and the block of the PARASCOPY module on the unmodified BAM, which must write the same table. The e2e case `cn-guard-1-mapq0` runs both blocks on the same reads with every MAPQ 0 read removed (`samtools view -q 1`): not assessed, `low_mapq0_share`, never a deletion, while the unmodified reads keep a share over 0.05. - The e2e case `hba-short-read-accuracy-1-nf-report` runs the Nextflow `HTML_REPORT` with and without an HBA table and checks the HBA line and the not-assessed list, then stub-runs every tool and checks that `main.nf` passes the table to each `HTML_REPORT` task that got the SMN table, against a copy of the pipeline with that join removed. - `tests/test_collect_summary.py` (`cn_report_section`) reads made-up CYP21A2, dosage and interval tables into the summary and both reports: a CYP21A2 deletion, HP 3, LCE3BC 1, CCL3L1 not assessed, a 22q11.2 deletion with its counts; without the tables each line says the table is not in the output; a header without the copy-number columns is unreadable for the dosage, interval, catalog and SMN tables, also with no row. - The e2e case `followups-nf-1-nf-wiring` stub-runs the pipeline on the stub pedigree and checks that `main.nf` passes the CYP21A2 and dosage tables to each `HTML_REPORT` task that got the interval table, and each partner's CYP21A2 table to `COUPLE_CARRIER`, against a copy of the pipeline without those lines. diff --git a/docs/cn-loci.md b/docs/cn-loci.md index 6c31e174..6c352885 100644 --- a/docs/cn-loci.md +++ b/docs/cn-loci.md @@ -85,7 +85,7 @@ The call has two parts: The sites, each CYP21A2's base in the reference with CYP21A1P's base at the matching position (dbSNP positions, the CYP21A1P bases read from the reference, checked on 2026-11-04): In2G c.293-13A/C>G (rs6467, chr6:32,039,081; CYP21A2 carries A or C there, CYP21A1P G), I172N c.518T>A (rs6475), V237E c.713T>A (rs12530380, the exon 6 cluster), V281L c.844G>T (rs6471), Q318X c.955C>T (rs7755898) and R356W c.1069C>T (rs7769409). -A site is low when its estimate is more than 0.6 copies under the larger of the median site's count and 2. `cyp21a2_copies` is the least count over the low sites, or the median site's count when no site is low. The `finding` is `deletion` when a site is low and the module total is under 4 (the 30 kb chimeric deletion leaves one module fewer and a chimera with CYP21A1P's sequence at its 5' sites), `conversion` at the low sites when the module total is 4 or more, `pseudogene_deletion` when the module total is under 4 with no low site, and `none` otherwise. At 30x a site has about 70 reads, and read sampling alone gives an estimate of 2 copies there a standard deviation of about 0.25. On HG002's 35x reads the six sites read 1.62 to 2.82 copies, a wider spread, so the 0.6 margin is a screen setting, not a measured false-call rate. Which chromosome carries a change is not claimed, and a two-copy result does not rule out two copies on one chromosome and none on the other. The row has no MAPQ 0 guard yet: an aligner that leaves out the reads that fit both genes (novoalign by default) would lower the module total and could read as a deletion, and no such BAM has been measured here to set a threshold. Use a BAM from minimap2 or BWA-MEM. Step 35 writes the table to `_cyp21a2_psv.tsv` ([step 35](35-paralogs.md#cyp21a2-from-depth-and-psv-sites)); the couple screen takes it with `--cyp21a2`. +A site is low when its estimate is more than 0.6 copies under the larger of the median site's count and 2. `cyp21a2_copies` is the least count over the low sites, or the median site's count when no site is low. The `finding` is `deletion` when a site is low and the module total is under 4 (the 30 kb chimeric deletion leaves one module fewer and a chimera with CYP21A1P's sequence at its 5' sites), `conversion` at the low sites when the module total is 4 or more, `pseudogene_deletion` when the module total is under 4 with no low site, and `none` otherwise. At 30x a site has about 70 reads, and read sampling alone gives an estimate of 2 copies there a standard deviation of about 0.25. On HG002's 35x reads the six sites read 1.62 to 2.82 copies, a wider spread, so the 0.6 margin is a screen setting, not a measured false-call rate. Which chromosome carries a change is not claimed, and a two-copy result does not rule out two copies on one chromosome and none on the other. The row has the MAPQ 0 guard of the HBA row: a second bedcov pass with `-Q 1`, and `mapq0_min` 0.02. HG002's 35x reads give a MAPQ 0 share of 0.099 over the two targets, from BWA-MEM and from minimap2 `-ax sr` alike (measured on 2026-11-08). An aligner that leaves out the reads that fit both genes (novoalign by default) gives about 0 and a lower module total estimate (4.22 to 3.81 on HG002's reads with every MAPQ 0 read removed), which on a genome whose depth runs lower reads as a deletion; under 0.02 the row is not assessed (`low_mapq0_share`), and its `message` row says the BAM has no multi-mapping reads there. The PSV sites still count every read whatever its mapping quality. Step 35 writes the table to `_cyp21a2_psv.tsv` ([step 35](35-paralogs.md#cyp21a2-from-depth-and-psv-sites)); the couple screen takes it with `--cyp21a2`. ### The interval rows @@ -130,5 +130,5 @@ pypgx owns the genes it calls from the BAM, the `bam` rows of `assets/pgx_genes. - `tests/test_copy_number.py` checks every row and refuses broken ones. For the three gene-dosage rows it calls made-up depths at 0 to 4 copies, cuts the dosage block out of step 35 and runs it on the made-up bedcov tables in `tests/fixtures/cn/` (HP 3, LCE3BC 1, CCL3L1 5 copies), checks that a row shifted by 1 kb no longer matches its table, and that the PARASCOPY module runs the same list. It runs the owner check with planted AMY1C, C4B and DEFB4 rows. It feeds made-up bedcov tables at 0, 1, 2, 3 and 4 copies and a MAPQ 0 share over the ceiling. It checks that a control with no read stops the run with `BROKEN_CONTROL`. It also checks that the depth tables of `tests/test_hba_copy_number.py` and `tests/test_lpa_kiv2.py` give, through the engine, the same tables the code before the engine wrote. - The fake-docker case `cn-engine-calls` runs steps 35, 41 and 46 with made-up depths. Each step calls `copy_number.py` once for its regions, step 35 and 46 once more for the call, and none calls the old per-gene scripts. A control exon with no read fails step 46 with `BROKEN_CONTROL`. -- `tests/test_copy_number.py` section 10 checks the CYP21A2 row: its six sites, the pileup parser, and the fixtures `tests/fixtures/cn/CYP21A2_{normal,deletion,conversion}` (2 copies and no finding; 1 copy with the module total 3; an In2G ratio of 0.5 and a possible conversion), with the conversion's depths read with every site even as the red control (no finding). The e2e case `cyp21a2-short-reads-1-hg002` runs the step's and the module's blocks on HG002's reads, unmodified and with a planted conversion and a planted deletion. +- `tests/test_copy_number.py` section 10 checks the CYP21A2 row: its six sites, the pileup parser, and the fixtures `tests/fixtures/cn/CYP21A2_{normal,deletion,conversion}` (2 copies and no finding; 1 copy with the module total 3; an In2G ratio of 0.5 and a possible conversion), with the conversion's depths read with every site even as the red control (no finding). The e2e case `cyp21a2-short-reads-1-hg002` runs the step's and the module's blocks on HG002's reads, unmodified and with a planted conversion and a planted deletion. The same section feeds the deletion's tables with MAPQ 0 shares of 0, 0.019, 0.021 and 0.1 (not assessed twice, then the deletion), and `cn-guard-1-mapq0` runs both blocks on HG002's reads with every MAPQ 0 read removed: not assessed (`low_mapq0_share`), never a deletion. - The e2e case `cn-engine-1-parity` runs steps 41 and 46 and the HBA block of step 35 on HG002's real reads at the three loci, at full depth and at an eighth of it. The code before the engine is fetched by its commit, 6056452, so the comparison still means something after this change is on main. The LPA and HBA tables and the RhD rows of step 41's table and evidence file must be byte for byte what it wrote. diff --git a/docs/lessons-learned.md b/docs/lessons-learned.md index 7d68e08c..8576ac66 100644 --- a/docs/lessons-learned.md +++ b/docs/lessons-learned.md @@ -620,6 +620,13 @@ Most bioinformatics containers run as non-root users. If writing to bind-mounted - **Why:** `ubuntu-latest` gives some legs a runner whose root disk is 72 GB. `/mnt` is a plain directory on that disk, so the work area pick found no second disk. The runner does have one: a blank NVMe disk (220.7 GB, no partition, no file system, nothing mounted), seen with `lsblk` in probe run 42037810203. - **Fix:** the leg's first step formats a blank whole disk and mounts it at `/mnt`; the work-area pick then takes `/mnt`, and the "Docker image store" step sets Docker's `data-root` to `/mnt/docker`. +## Copy number from depth (2026-11) + +### A BAM without MAPQ 0 reads could read as a CYP21A2 deletion +- **Failed:** the CYP21A2 row of `bin/cn_loci.tsv` added up the depth over CYP21A2 and CYP21A1P with no mapping quality check. The two are about 98% alike, so a read that fits both gets MAPQ 0. An aligner that leaves those reads out (novoalign does by default) lowers the module total. On HG002's reads, removing them takes the estimate from 4.22 to 3.81; on a genome whose depth runs lower it rounds to 3, which with a low PSV site reads as the 30 kb chimeric deletion: a false carrier finding for congenital adrenal hyperplasia. +- **Why:** HBA had the guard (a second `samtools bedcov -Q 1` pass and `mapq0_min`); the CYP21A2 row was added later without it. On HG002's 35x reads, 9.9% of the depth over the two targets comes from MAPQ 0 reads, with BWA-MEM and with minimap2 `-ax sr` alike, so a share near 0 means the aligner dropped them. +- **Fix:** the row runs both bedcov passes with `mapq0_min` 0.02. Under that share it is not assessed (`low_mapq0_share`), and the table's `message` row, which both reports print, says the BAM has no multi-mapping reads there. `tests/test_copy_number.py` and the e2e case `cn-guard-1-mapq0` check it. + ## Downloads (2026-11) ### One bad transfer of a large file failed the whole setup diff --git a/modules/local/parascopy/main.nf b/modules/local/parascopy/main.nf index 6263776b..16da880f 100644 --- a/modules/local/parascopy/main.nf +++ b/modules/local/parascopy/main.nf @@ -212,19 +212,22 @@ process PARASCOPY { if copy_number.py bed --locus CYP21A2 > cyp21a2_regions.bed \\ && copy_number.py bed --psv --locus CYP21A2 > cyp21a2_psv.bed \\ && samtools bedcov cyp21a2_regions.bed ${bam} > cyp21a2_bedcov.tsv \\ + && samtools bedcov -Q 1 cyp21a2_regions.bed ${bam} > cyp21a2_bedcov_q1.tsv \\ && samtools mpileup -A -B -q 0 -Q 13 -d 100000 \\ -r "\$(awk 'NR == 1 {c = \$1; s = \$2 + 1} {e = \$3} END {print c ":" s "-" e}' cyp21a2_psv.bed)" \\ -l cyp21a2_psv.bed ${bam} > cyp21a2_pileup.tsv \\ - && copy_number.py call --locus CYP21A2 --bedcov cyp21a2_bedcov.tsv --pileup cyp21a2_pileup.tsv \\ + && copy_number.py call --locus CYP21A2 --bedcov cyp21a2_bedcov.tsv --bedcov-mapq1 cyp21a2_bedcov_q1.tsv \\ + --pileup cyp21a2_pileup.tsv \\ --table ${meta.id}_cyp21a2_psv.tsv.tmp; then mv ${meta.id}_cyp21a2_psv.tsv.tmp ${meta.id}_cyp21a2_psv.tsv else rm -f ${meta.id}_cyp21a2_psv.tsv.tmp echo "WARNING: the CYP21A2 depth and PSV check failed; ${meta.id}_cyp21a2_psv.tsv says not assessed." >&2 - printf 'metric\\tvalue\\nstatus\\tnot assessed\\ncyp21a2_copies\\t.\\nmodule_total\\t.\\nexpected_total\\t.\\nestimate\\t.\\nfinding\\t.\\nconversion_sites\\t.\\ntarget_depth\\t.\\ncontrol_depth\\t.\\nreason\\tcheck_failed\\nreading\\t%s\\n' \\ + printf 'metric\\tvalue\\nstatus\\tnot assessed\\ncyp21a2_copies\\t.\\nmodule_total\\t.\\nexpected_total\\t.\\nestimate\\t.\\nfinding\\t.\\nconversion_sites\\t.\\ntarget_depth\\t.\\ncontrol_depth\\t.\\nreason\\tcheck_failed\\nreading\\t%s\\nmapq0_share\\t.\\nmessage\\tnot assessed: %s\\n' \\ + 'the CYP21A2 depth and PSV check failed; see the log of step 35' \\ 'the CYP21A2 depth and PSV check failed; see the log of step 35' > ${meta.id}_cyp21a2_psv.tsv fi - rm -f cyp21a2_regions.bed cyp21a2_psv.bed cyp21a2_bedcov.tsv cyp21a2_pileup.tsv + rm -f cyp21a2_regions.bed cyp21a2_psv.bed cyp21a2_bedcov.tsv cyp21a2_bedcov_q1.tsv cyp21a2_pileup.tsv # Interval rows of bin/cn_loci.tsv (as scripts/35-paralogs.sh): one bedcov over every row's segments, # then the verdict; a failed pass writes every row as not assessed (check_failed) diff --git a/scripts/35-paralogs.sh b/scripts/35-paralogs.sh index e984661b..64d609f9 100755 --- a/scripts/35-paralogs.sh +++ b/scripts/35-paralogs.sh @@ -336,7 +336,9 @@ rm -rf "$DOSAGE_WORK" # single-copy flanks: the RCCX module total), samtools mpileup without a # reference at its paralogous sequence variant sites (each on CYP21A2 and on # CYP21A1P, every read whatever its MAPQ), then copy_number.py call: the -# CYP21A2 copies at each site, and a deletion or a possible conversion. A +# CYP21A2 copies at each site, and a deletion or a possible conversion. The +# bedcov runs twice, all MAPQ and MAPQ 1 or more (as HBA): a BAM without the +# MAPQ 0 reads that fit both copies reads not assessed, never a deletion. A # failure writes a table that says not assessed, and the step goes on. CYP21_OUT="${OUTDIR}/${SAMPLE}_cyp21a2_psv.tsv" CYP21_WORK="${OUTDIR}/CYP21A2.part" @@ -348,18 +350,22 @@ if run_in -v "${PGP_ROOT}/bin:/pgp-bin:ro" "${PARASCOPY_IMAGE}" \ python3 /pgp-bin/copy_number.py bed --psv --locus CYP21A2 > "${CYP21_WORK}/psv.bed" \ && run_in --cpus 2 --memory 2g "${PARASCOPY_IMAGE}" \ samtools bedcov "$(cpath "${CYP21_WORK}/regions.bed")" "$(cpath "$BAM")" > "${CYP21_WORK}/bedcov.tsv" \ + && run_in --cpus 2 --memory 2g "${PARASCOPY_IMAGE}" \ + samtools bedcov -Q 1 "$(cpath "${CYP21_WORK}/regions.bed")" "$(cpath "$BAM")" > "${CYP21_WORK}/bedcov_q1.tsv" \ && run_in --cpus 2 --memory 2g "${PARASCOPY_IMAGE}" \ samtools mpileup -A -B -q 0 -Q 13 -d 100000 \ -r "$(awk 'NR == 1 {c = $1; s = $2 + 1} {e = $3} END {print c ":" s "-" e}' "${CYP21_WORK}/psv.bed")" \ -l "$(cpath "${CYP21_WORK}/psv.bed")" "$(cpath "$BAM")" > "${CYP21_WORK}/pileup.tsv" \ && run_in -v "${PGP_ROOT}/bin:/pgp-bin:ro" "${PARASCOPY_IMAGE}" \ python3 /pgp-bin/copy_number.py call --locus CYP21A2 --bedcov "$(cpath "${CYP21_WORK}/bedcov.tsv")" \ - --pileup "$(cpath "${CYP21_WORK}/pileup.tsv")" --table "$(cpath "${CYP21_OUT}.tmp")"; then + --bedcov-mapq1 "$(cpath "${CYP21_WORK}/bedcov_q1.tsv")" --pileup "$(cpath "${CYP21_WORK}/pileup.tsv")" \ + --table "$(cpath "${CYP21_OUT}.tmp")"; then mv "${CYP21_OUT}.tmp" "$CYP21_OUT" else rm -f "${CYP21_OUT}.tmp" echo "WARNING: the CYP21A2 depth and PSV check failed; ${CYP21_OUT##*/} says not assessed." >&2 - printf 'metric\tvalue\nstatus\tnot assessed\ncyp21a2_copies\t.\nmodule_total\t.\nexpected_total\t.\nestimate\t.\nfinding\t.\nconversion_sites\t.\ntarget_depth\t.\ncontrol_depth\t.\nreason\tcheck_failed\nreading\t%s\n' \ + printf 'metric\tvalue\nstatus\tnot assessed\ncyp21a2_copies\t.\nmodule_total\t.\nexpected_total\t.\nestimate\t.\nfinding\t.\nconversion_sites\t.\ntarget_depth\t.\ncontrol_depth\t.\nreason\tcheck_failed\nreading\t%s\nmapq0_share\t.\nmessage\tnot assessed: %s\n' \ + 'the CYP21A2 depth and PSV check failed; see the log of step 35' \ 'the CYP21A2 depth and PSV check failed; see the log of step 35' > "$CYP21_OUT" fi rm -rf "$CYP21_WORK" diff --git a/tests/test_copy_number.py b/tests/test_copy_number.py index bb812ac2..b887c92a 100644 --- a/tests/test_copy_number.py +++ b/tests/test_copy_number.py @@ -495,7 +495,8 @@ def read(p): check(f"{path}: copy_number.py call --locus {loc}", f"copy_number.py call --locus {loc}" in text) passes = ["" if p == "." else p for p in cat[loc]["_passes"]] if loc == "HBA": # and the dosage loop's two passes (section 8), the interval pass (section 9) - passes += ["", "-Q 1", "", ""] # and the CYP21A2 pass (section 10) + passes += ["", "-Q 1", "", ""] # and the CYP21A2 passes (section 10) + passes += ["-Q 1"] # the commands only: options, then a path, a quoted argument or a variable calls = [o for o, arg in re.findall(r"samtools bedcov((?: -[A-Za-z](?: \d+)?)*) (\S+)", text) if re.match(r'["$/]|\S+\.bed\b', arg)] @@ -864,13 +865,17 @@ def cyp_fixture(r, total, a2): return bed, "".join(x for _, x in sorted(lines)) -def cyp_call(r, bedcov_text, pileup_text): +def cyp_call(r, bedcov_text, pileup_text, mapq0=0.1): + """The block of a CYP21A2 genome; mapq0 is the share of every segment's + depth from MAPQ 0 reads (0.1 as HG002's BWA-MEM and minimap2 reads give + over the two targets; 0 for an aligner that leaves those reads out).""" with tempfile.TemporaryDirectory() as t: b, p = os.path.join(t, "b.tsv"), os.path.join(t, "p.tsv") open(b, "w").write(bedcov_text) open(p, "w").write(pileup_text) m = cn.read_bedcov(b, r) - return cn.psv_assess(r, cn.assess(r, m), cn.read_pileup(p, r)) + q1 = {n: v * (1 - mapq0) for n, v in m.items()} + return cn.psv_assess(r, cn.assess(r, m, q1), cn.read_pileup(p, r)) # run_in and cpath for the CYP21A2 block of step 35: copy_number.py runs for real, @@ -881,7 +886,14 @@ def cyp_call(r, bedcov_text, pileup_text): shift case "$1 $2" in "python3 /pgp-bin/copy_number.py") shift 2; python3 "${REPO}/bin/copy_number.py" "$@" ;; - "samtools bedcov") cat "${FIX}/CYP21A2_${CYP_CASE}.bedcov.tsv" ;; + "samtools bedcov") + # the MAPQ >= 1 pass: 90% of each depth (HG002's share), or all of it for + # an aligner that left the MAPQ 0 reads out (CYP_NO_MAPQ0) + if [ "${3:-}" = "-Q" ] && [ -z "${CYP_NO_MAPQ0:-}" ]; then + awk -F'\t' -v OFS='\t' '{$5 = int($5 * 0.9); print}' "${FIX}/CYP21A2_${CYP_CASE}.bedcov.tsv" + else + cat "${FIX}/CYP21A2_${CYP_CASE}.bedcov.tsv" + fi ;; "samtools mpileup") [ -z "${MPILEUP_FAILS:-}" ] && cat "${FIX}/CYP21A2_${CYP_CASE}.pileup.tsv" ;; *) return 1 ;; esac @@ -993,17 +1005,46 @@ def psv(cat, read): f"\t{round(30.0 * 2689)}\n", f"\t{round(37.5 * 2689)}\n"), pile) check("CYP21A2: a module total between counts (4.5) is not assessed, no finding", (b["status"], b["reason"], b["finding"]) == ("not assessed", "between_counts", None), b) + # The MAPQ 0 guard: an aligner that left out the reads fitting both copies lowers + # the module total, so the deletion's tables from such a BAM (no MAPQ 0 read) read + # not assessed, not a deletion; with HG002's 10% share they still give the deletion. + check("CYP21A2 row: a MAPQ 0 guard (two bedcov passes, all reads and MAPQ >= 1; minimum share 0.02)", + r["_passes"] == [".", "-Q 1"] and r["_mapq0_min"] == 0.02 and r["_mapq0_max"] is None, + (r["_passes"], r["_mapq0_min"], r["_mapq0_max"])) + total, a2 = CYP_CASES["deletion"] + bed, pile = cyp_fixture(r, total, a2) + for share, want in ((0.0, ("not assessed", "low_mapq0_share", None, None)), + (0.019, ("not assessed", "low_mapq0_share", None, None)), + (0.021, ("called", None, "deletion", 1)), (0.1, ("called", None, "deletion", 1))): + b = cyp_call(r, bed, pile, share) + check(f"CYP21A2: the deletion's tables with a MAPQ 0 share of {share}: {want[0]} ({want[1] or want[2]})", + (b["status"], b["reason"], b["finding"], b["copies"]) == want, b) + b = cyp_call(r, bed, pile, 0.0) + rows, line = cn.psv_table(r, b) + tab = dict(rows) + check("CYP21A2 without MAPQ 0 reads: the table says so (message, share 0.000), no copies, module total or finding", + tab["message"].startswith("not assessed: the BAM has no multi-mapping reads here") + and "CYP21A2 and CYP21A1P" in tab["message"] and tab["mapq0_share"] == "0.000" + and (tab["cyp21a2_copies"], tab["module_total"], tab["finding"]) == (".", ".", ".") + and "low_mapq0_share" in line, (tab, line)) + bed, pile = cyp_fixture(r, 4, {}) + tab = dict(cn.psv_table(r, cyp_call(r, bed, pile))[0]) + check("CYP21A2 unmodified with HG002's MAPQ 0 share: called, share 0.100, message '.'", + (tab["status"], tab["finding"], tab["mapq0_share"], tab["message"]) == ("called", "none", "0.100", "."), tab) step, module = read("scripts/35-paralogs.sh"), read("modules/local/parascopy/main.nf") blk = re.search(r'\n(CYP21_OUT=.*?\nrm -rf "\$CYP21_WORK"\n)', step, re.S) check("step 35 holds the CYP21A2 block", blk is not None) tabs = {} with tempfile.TemporaryDirectory() as t: - for case in list(CYP_CASES) + ["fail"]: + for case in list(CYP_CASES) + ["fail", "no_mapq0"]: env = dict(os.environ, REPO=REPO, FIX=CYP_FIX, OUTDIR=t, SAMPLE="S", PGP_ROOT=REPO, PARASCOPY_IMAGE="img", - BAM=os.path.join(t, "S.bam"), CYP_CASE="normal" if case == "fail" else case) + BAM=os.path.join(t, "S.bam"), + CYP_CASE={"fail": "normal", "no_mapq0": "deletion"}.get(case, case)) if case == "fail": env["MPILEUP_FAILS"] = "1" + if case == "no_mapq0": # the deletion's reads from a BAM with no MAPQ 0 read + env["CYP_NO_MAPQ0"] = "1" out = subprocess.run(["bash", "-euo", "pipefail", "-c", STUB_CYP + (blk.group(1) if blk else "false\n")], capture_output=True, text=True, env=env) tab = os.path.join(t, "S_cyp21a2_psv.tsv") @@ -1014,13 +1055,25 @@ def psv(cat, read): if os.path.exists(tab): os.replace(tab, tabs[case]) want = {"normal": ("called", "2", "4", "none"), "deletion": ("called", "1", "3", "deletion"), - "conversion": ("called", "1", "4", "conversion"), "fail": ("not assessed", ".", ".", ".")}[case] + "conversion": ("called", "1", "4", "conversion"), "fail": ("not assessed", ".", ".", "."), + "no_mapq0": ("not assessed", ".", ".", ".")}[case] check(f"step 35 CYP21A2 block ({case}): status, copies, module total, finding {want}", tuple(rows.get(k) for k in ("status", "cyp21a2_copies", "module_total", "finding")) == want, rows) check("step 35 CYP21A2 block with a failed mpileup: reason check_failed", open(tabs["fail"]).read().count("reason\tcheck_failed") == 1 if os.path.exists(tabs["fail"]) else False) + rd = (lambda c: dict(x.split("\t", 1) for x in open(tabs[c]).read().splitlines()[1:]) + if os.path.exists(tabs[c]) else {}) + nm, nl = rd("no_mapq0"), rd("normal") + check("step 35 CYP21A2 block on the deletion's reads without MAPQ 0 reads: low_mapq0_share, the message says why", + (nm.get("reason"), nm.get("mapq0_share")) == ("low_mapq0_share", "0.000") + and nm.get("message", "").startswith("not assessed: the BAM has no multi-mapping reads here"), nm) + check("step 35 CYP21A2 block, unmodified: the MAPQ >= 1 pass ran (share 0.100), message '.'", + (nl.get("mapq0_share"), nl.get("message")) == ("0.100", "."), nl) + check("step 35 CYP21A2 block with a failed mpileup: the rows of a called table but the sites'", + list(rd("fail")) == [k for k in nl if not k.startswith("psv_")], (list(rd("fail")), list(nl))) for path, text, q in (("scripts/35-paralogs.sh", step, True), ("modules/local/parascopy/main.nf", module, False)): for cmd in ("copy_number.py bed --locus CYP21A2", "copy_number.py bed --psv --locus CYP21A2", + "--bedcov-mapq1", "cyp21a2_bedcov_q1.tsv" if not q else 'bedcov_q1.tsv")"', "samtools mpileup -A -B -q 0 -Q 13 -d 100000", "copy_number.py call --locus CYP21A2 --bedcov", "--pileup", "_cyp21a2_psv.tsv"): check(f"{path}: {cmd}", cmd in text) @@ -1039,6 +1092,9 @@ def psv(cat, read): check("couple screen: the unmodified table is no finding, 2 copies read", got.get("normal", (["x"], ""))[0] == [] and "CYP21A2 2 copies" in got.get("normal", ([], ""))[1], got.get("normal")) + check("couple screen: the table without MAPQ 0 reads is no finding and says not assessed (low_mapq0_share)", + got.get("no_mapq0", (["x"], ""))[0] == [] and got.get("no_mapq0", ([], ""))[1] + == "CYP21A2 not assessed (low_mapq0_share)", got.get("no_mapq0")) check("couple screen: the failed table is no finding and says not assessed", got.get("fail", (["x"], ""))[0] == [] and got.get("fail", ([], ""))[1].startswith("CYP21A2 not assessed"), got.get("fail")) @@ -1083,6 +1139,14 @@ def psv(cat, read): c = data.get("cyp21a2") or {} check("collect_summary: step 35's section carries the CYP21A2 table (deletion, 1 copy)", (c.get("finding"), c.get("cyp21a2_copies")) == ("deletion", "1"), c) + # and the reports print why a table without MAPQ 0 reads is not assessed + import render_report + if os.path.exists(tabs["no_mapq0"]): + os.replace(tabs["no_mapq0"], os.path.join(sd, "paralogs", "S_cyp21a2_psv.tsv")) + _, data = cs.sec_parascopy(sd, "S") + got = render_report.cyp21a2_line(data.get("cyp21a2")) + check("report: the CYP21A2 line without MAPQ 0 reads is not assessed and says the BAM has no multi-mapping reads", + ": not assessed (the BAM has no multi-mapping reads here" in got and "deletion," not in got, got) if __name__ == "__main__": From cd2e89c478ce73b4d909f6e1ebac0b4307f3de92 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sun, 8 Nov 2026 15:47:15 +0100 Subject: [PATCH 2/3] test(e2e): run the CYP21A2 MAPQ 0 guard on HG002's reads with the MAPQ 0 reads removed Bead: pgp-5bja.21 --- scripts/ci/e2e-case-paths.tsv | 1 + scripts/ci/e2e-case-times.tsv | 1 + tests/e2e/cn-guard-1-mapq0.sh | 96 +++++++++++++++++++++++++++++++++++ 3 files changed, 98 insertions(+) create mode 100755 tests/e2e/cn-guard-1-mapq0.sh diff --git a/scripts/ci/e2e-case-paths.tsv b/scripts/ci/e2e-case-paths.tsv index 32f1f190..7008d1f1 100644 --- a/scripts/ci/e2e-case-paths.tsv +++ b/scripts/ci/e2e-case-paths.tsv @@ -73,6 +73,7 @@ clinvar-clip-check-1 assets/construct_variants.tsv assets/deep_intronic_sites.ts clinvar-combined-allele-1 assets/construct_variants.tsv assets/deep_intronic_sites.tsv assets/stub/sample.bam assets/stub/sample.bam.bai bin/cdna_scan.py bin/clinvar_hits.awk bin/clinvar_join_adjacent.awk bin/clip_check.py modules/local/clinvar_screen/* modules/local/vcf_precheck/* scripts/06-clinvar-screen.sh scripts/24-html-report.sh scripts/ci/fake-docker/lib.sh scripts/generate-report.sh scripts/validate-setup.sh tests/fixtures/vcf/combined_allele_clinvar.vcf tests/fixtures/vcf/combined_allele_sample.vcf clinvar-cpsr-quality-1 assets/construct_variants.tsv assets/deep_intronic_sites.tsv bin/cdna_scan.py bin/clinvar_hits.awk bin/clinvar_join_adjacent.awk bin/clip_check.py bin/render_report.py scripts/06-clinvar-screen.sh scripts/ci/fake-docker/lib.sh cn-engine-1-parity assets/bloodtype_markers.tsv bin/bloodtype_report.py bin/cn_loci.tsv bin/copy_number.py bin/hba_copy_number.py bin/lpa_kiv2.py scripts/41-bloodtype.sh scripts/46-lpa.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh +cn-guard-1-mapq0 bin/cn_loci.tsv bin/copy_number.py modules/local/parascopy/* scripts/35-paralogs.sh cn-intervals-1-catalog bin/cn_loci.tsv scripts/ci/build-cn-intervals.py scripts/ci/fake-docker/lib.sh cn-intervals-2-hg002 assets/par_grch38.bed assets/paralog_loci.tsv bin/cn_loci.tsv bin/copy_number.py bin/hba_copy_number.py bin/lpa_kiv2.py modules/local/deepvariant/* modules/local/parascopy/* modules/local/vcf_precheck/* scripts/35-paralogs.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh cn-report-1-tables assets/carrier_genes.tsv assets/par_grch38.bed bin/collect_summary.py bin/couple_report.py bin/family_report.py bin/paraphase_report.py bin/render_report.py bin/results_index.py scripts/24-html-report.sh scripts/39-couple.sh scripts/ci/fake-docker/lib.sh scripts/generate-report.sh diff --git a/scripts/ci/e2e-case-times.tsv b/scripts/ci/e2e-case-times.tsv index 9354daf2..3ff5ba01 100644 --- a/scripts/ci/e2e-case-times.tsv +++ b/scripts/ci/e2e-case-times.tsv @@ -72,6 +72,7 @@ clinvar-clip-check-1 shard 8 - clinvar-combined-allele-1 shard 25 - clinvar-cpsr-quality-1 shard 5 clinvar-clip-check-1 cn-engine-1-parity shard 33 - +cn-guard-1-mapq0 shard 15 - cn-intervals-1-catalog shard 0 - cn-intervals-2-hg002 shard 19 - cn-report-1-tables shard 60 - diff --git a/tests/e2e/cn-guard-1-mapq0.sh b/tests/e2e/cn-guard-1-mapq0.sh new file mode 100755 index 00000000..957bf473 --- /dev/null +++ b/tests/e2e/cn-guard-1-mapq0.sh @@ -0,0 +1,96 @@ +#!/usr/bin/env bash +# The MAPQ 0 guard of the CYP21A2 row of bin/cn_loci.tsv on real reads: HG002 +# from the 35x BWA-MEM GRCh38 no-ALT BAM of DeepVariant's case study, over the +# row's segments (CYP21A2, CYP21A1P, the two flanks). Two BAMs: +# - unmodified: the MAPQ 0 share over the two targets is over 0.05 (HG002's +# reads give about 0.1), and the row is never not assessed for it; +# - every MAPQ 0 read removed (samtools view -q 1), as an aligner that leaves +# out the reads fitting both copies writes it: not assessed, +# low_mapq0_share, no copies, no module total, no finding, and a message +# that says the BAM has no multi-mapping reads. +# Red control: the BAM without MAPQ 0 reads gives a lower estimate (the +# table keeps it when not assessed), so the guard stops a real drift. +# The block of step 35 runs on each (cut from scripts/35-paralogs.sh, with +# scripts/lib/common.sh's run_in), and the PARASCOPY module's block, cut from +# its main.nf, writes the same table on the BAM without MAPQ 0 reads, byte for +# byte. +. "$(dirname "$0")/lib.sh" + +SRC="https://storage.googleapis.com/deepvariant/case-study-testdata/HG002.novaseq.pcr-free.35x.dedup.grch38_no_alt.bam" +P=HG002cnguard +W="${GENOME_DIR}/${P}" +rm -rf "$W"; mkdir -p "${W}/aligned" "${W}/paralogs" +python3 "${REPO}/bin/copy_number.py" bed --locus CYP21A2 | sort -k1,1 -k2,2n > "${W}/fetch.bed" +check_ge "segments of the CYP21A2 row to fetch" "$(wc -l < "${W}/fetch.bed" | tr -d ' ')" 10 + +retry() { local i; for i in 1 2 3 4; do "$@" && return 0; [ "$i" -lt 4 ] && { echo " try ${i} of 4 failed; again in $((i * 15)) s"; sleep $((i * 15)); }; done; return 1; } +fetch_reads() { + rm -f "${W}/unmodified.bam" "${W}/$(basename "$SRC").bai" + docker run --rm -u "$(id -u):$(id -g)" -e HOME=/tmp -v /etc/ssl/certs:/etc/ssl/certs:ro \ + -e CURL_CA_BUNDLE=/etc/ssl/certs/ca-certificates.crt -v "${GENOME_DIR}:/genome" -w "/genome/${P}" "$SAMTOOLS_IMAGE" \ + samtools view -b -M -L "/genome/${P}/fetch.bed" -o "/genome/${P}/unmodified.bam" "$SRC" +} +check "HG002's reads over the CYP21A2 row download" retry fetch_reads +rm -f "${W}/$(basename "$SRC").bai" +sam index "${P}/unmodified.bam" +sam view -b -q 1 -o "${P}/no_mapq0.bam" "${P}/unmodified.bam" && sam index "${P}/no_mapq0.bam" +all=$(sam view -c "${P}/unmodified.bam" 2>/dev/null || echo 0) +kept=$(sam view -c "${P}/no_mapq0.bam" 2>/dev/null || echo 0) +check_ge "HG002 reads over the CYP21A2 row" "$all" 5000 +check_eq "no read of MAPQ 0 is left in the BAM without them" \ + "$(sam view -c -e 'mapq == 0' "${P}/no_mapq0.bam" 2>/dev/null | tr -d ' ')" 0 +check "the BAM without MAPQ 0 reads lost some (${kept} of ${all} kept)" awk -v a="$all" -v k="$kept" 'BEGIN {exit !(k > 0 && k < a)}' + +col() { awk -F'\t' -v k="$2" '$1 == k {print $2}' "$1" 2>/dev/null; } + +# --- step 35's CYP21A2 block, with common.sh's run_in --------------------------------- +sed -n '/^CYP21_OUT=/,/^rm -rf "\$CYP21_WORK"$/p' "${REPO}/scripts/35-paralogs.sh" > "${CASE_TMP}/step_body.sh" +check_ge "step 35's CYP21A2 block was cut from the script" "$(wc -l < "${CASE_TMP}/step_body.sh" | tr -d ' ')" 20 +check "step 35's CYP21A2 block runs the MAPQ 1 pass" grep -q -- '--bedcov-mapq1' "${CASE_TMP}/step_body.sh" +step_block() { # step_block NAME: the block on ${P}/.bam, its table copied to CASE_TMP/.tsv + { echo 'set -euo pipefail' + printf 'SAMPLE=%q GENOME_DIR=%q\n' "$P" "$GENOME_DIR" + printf '. %q\n' "${REPO}/scripts/lib/common.sh" + printf 'OUTDIR=%q BAM=%q\n' "${W}/paralogs" "${W}/$1.bam" + cat "${CASE_TMP}/step_body.sh"; } > "${CASE_TMP}/block.sh" + echo "+ step 35's CYP21A2 block ($1)" + bash "${CASE_TMP}/block.sh" 2>&1 | tail -n 5 + cp "${W}/paralogs/${P}_cyp21a2_psv.tsv" "${CASE_TMP}/$1.tsv" 2>/dev/null || : > "${CASE_TMP}/$1.tsv" + cat "${CASE_TMP}/$1.tsv" +} +for v in unmodified no_mapq0; do step_block "$v"; done +U="${CASE_TMP}/unmodified.tsv" N="${CASE_TMP}/no_mapq0.tsv" + +check "unmodified: the MAPQ 0 share ($(col "$U" mapq0_share)) is over 0.05" \ + awk -v s="$(col "$U" mapq0_share)" 'BEGIN {exit !(s != "" && s != "." && s + 0 > 0.05)}' +check "unmodified: never not assessed for the MAPQ 0 share" lacks 'low_mapq0_share' "$(col "$U" reason)" +check_eq "without MAPQ 0 reads: not assessed (low_mapq0_share), share 0.000, no copies, module total or finding" \ + "$(col "$N" status)/$(col "$N" reason)/$(col "$N" mapq0_share)/$(col "$N" cyp21a2_copies)/$(col "$N" module_total)/$(col "$N" finding)" \ + "not assessed/low_mapq0_share/0.000/././." +check "without MAPQ 0 reads: the message says the BAM has no multi-mapping reads" \ + has '^not assessed: the BAM has no multi-mapping reads here' "$(col "$N" message)" + +# Red control: the estimate the guard withheld. Removing the MAPQ 0 reads lowers +# the module total estimate (4.22 to 3.81 on these reads, measured 2026-11-08), +# the drift that, on a genome whose depth runs lower, reads as one module fewer. +check "red control: without MAPQ 0 reads the estimate ($(col "$N" estimate)) is at least 0.2 under the unmodified one ($(col "$U" estimate))" \ + awk -v n="$(col "$N" estimate)" -v u="$(col "$U" estimate)" 'BEGIN {exit !(n != "" && n != "." && u != "" && u != "." && n + 0.2 <= u + 0)}' +echo "cn-guard: HG002 MAPQ 0 share $(col "$U" mapq0_share) ($(col "$U" status), module total $(col "$U" module_total)," \ + "estimate $(col "$U" estimate)); without MAPQ 0 reads ${kept} of ${all} reads, $(col "$N" status) ($(col "$N" reason));" \ + "estimate $(col "$N" estimate)" >> "$E2E_NOTES" + +# --- the PARASCOPY module's CYP21A2 block, in its image ------------------------------- +MOD="${REPO}/modules/local/parascopy/main.nf" +sed -n '/# CYP21A2 from depth and PSVs (as scripts\/35-paralogs.sh)/,/rm -f cyp21a2_regions.bed/p' "$MOD" \ + | sed -e 's/\\\\/\\/g' -e 's/\\\$/$/g' > "${CASE_TMP}/module_body.sh" +check_ge "the module's CYP21A2 block was cut from main.nf" "$(wc -l < "${CASE_TMP}/module_body.sh" | tr -d ' ')" 15 +M="${W}/module" +rm -rf "$M"; mkdir -p "$M" +sed -e "s|\${bam}|/genome/${P}/no_mapq0.bam|g" -e "s|\${meta.id}|${P}|g" "${CASE_TMP}/module_body.sh" > "${M}/block.sh" +docker run --rm --user "$(id -u):$(id -g)" -v "${GENOME_DIR}:/genome" -v "${REPO}/bin:/pgp-bin:ro" -w "/genome/${P}/module" \ + "$PARASCOPY_IMAGE" sh -c 'PATH="/pgp-bin:$PATH"; sh -e block.sh' 2>&1 | tail -n 3 +check "the module's block writes step 35's table on the BAM without MAPQ 0 reads, byte for byte" \ + cmp -s "$N" "${M}/${P}_cyp21a2_psv.tsv" + +rm -rf "$W" +finish From 3583bec2fd9baffa2c2ccceaf8a45191d570ecfb Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Sun, 8 Nov 2026 18:20:04 +0100 Subject: [PATCH 3/3] fix(cn): say a low MAPQ 0 share means few or no multi-mapping reads, not none The guard trips under a 2% share, so the message names the measured share and the likely cause instead of claiming the BAM holds no such read. The e2e row takes the 5 s the case took on its first run. Bead: pgp-5bja.21 --- bin/copy_number.py | 10 +++++----- docs/35-paralogs.md | 2 +- docs/cn-loci.md | 2 +- docs/lessons-learned.md | 2 +- scripts/ci/e2e-case-times.tsv | 2 +- tests/e2e/cn-guard-1-mapq0.sh | 6 +++--- tests/test_copy_number.py | 8 ++++---- 7 files changed, 16 insertions(+), 16 deletions(-) diff --git a/bin/copy_number.py b/bin/copy_number.py index 92bd31f1..177ecd02 100755 --- a/bin/copy_number.py +++ b/bin/copy_number.py @@ -588,11 +588,11 @@ def psv_message(row, b): return "." if b["reason"] == "low_mapq0_share": targets = " and ".join(n.removesuffix("_BODY") for n, *_ in row["_targets"]) - return (f"not assessed: the BAM has no multi-mapping reads here (reads of mapping quality 0 make up " - f"{100 * b['mapq0_share']:.1f}% of the depth over {targets}, " - f"under {100 * row['_mapq0_min']:g}%), so the aligner left out the reads that fit both copies (as " - "novoalign does by default) and the module total would come out low and read as a deletion; " - "realign with minimap2 (step 02) or BWA-MEM") + return (f"not assessed: the BAM has few or no multi-mapping reads here (reads of mapping quality 0 make " + f"up {100 * b['mapq0_share']:.1f}% of the depth over {targets}, under " + f"{100 * row['_mapq0_min']:g}%), the sign of an aligner that leaves out the reads that fit both " + "copies (as novoalign does by default); the module total may then come out low and read as a " + "deletion; realign with minimap2 (step 02) or BWA-MEM") return f"not assessed: {b['reason'] or '.'}" diff --git a/docs/35-paralogs.md b/docs/35-paralogs.md index 79ff7035..f7751387 100644 --- a/docs/35-paralogs.md +++ b/docs/35-paralogs.md @@ -228,7 +228,7 @@ The `finding`: | `conversion` | the module total is 4 or more and a site is low: its estimate is more than 0.6 copies under the larger of the median site's count and 2 | a possible conversion at that site; `conversion_sites` names it | | `pseudogene_deletion` | the module total is under 4 and no site is low | a CYP21A1P deletion, common and not a carrier finding | -The ratio of a site is its estimate over the median of the sites, so one copy converted out of two reads 0.5. A site needs 30 reads (CYP21A2 and CYP21A1P bases together); with no site there the row is not assessed (`psv_low_depth`). Which chromosome carries a change is never said, and a change away from the six sites is not seen. The row has a mapping quality 0 check, as HBA has. A read that fits CYP21A2 and CYP21A1P equally gets mapping quality 0, and HG002's 35x reads give about 10% of the depth over the two targets from such reads, with BWA-MEM and with minimap2 alike. An aligner that leaves those reads out (novoalign by default) gives a lower module total (on HG002's reads the estimate falls from 4.22 to 3.81), which can read as a deletion. So when they make up under 2% of that depth, the row is not assessed (`low_mapq0_share`), and the `message` row and both reports say the BAM has no multi-mapping reads there. Realign with minimap2 or BWA-MEM. A PacBio HiFi BAM gets CYP21A2 from Paraphase (step 35b) instead, which phases the copies. +The ratio of a site is its estimate over the median of the sites, so one copy converted out of two reads 0.5. A site needs 30 reads (CYP21A2 and CYP21A1P bases together); with no site there the row is not assessed (`psv_low_depth`). Which chromosome carries a change is never said, and a change away from the six sites is not seen. The row has a mapping quality 0 check, as HBA has. A read that fits CYP21A2 and CYP21A1P equally gets mapping quality 0, and HG002's 35x reads give about 10% of the depth over the two targets from such reads, with BWA-MEM and with minimap2 alike. An aligner that leaves those reads out (novoalign by default) gives a lower module total (on HG002's reads the estimate falls from 4.22 to 3.81), which can read as a deletion. So when they make up under 2% of that depth, the row is not assessed (`low_mapq0_share`), and the `message` row and both reports say the BAM has few or no multi-mapping reads there. Realign with minimap2 or BWA-MEM. A PacBio HiFi BAM gets CYP21A2 from Paraphase (step 35b) instead, which phases the copies. **A research estimate, not a diagnostic test.** A deletion or a conversion needs a clinical test (MLPA, long-range PCR or targeted sequencing). diff --git a/docs/cn-loci.md b/docs/cn-loci.md index 6c352885..6c1d898a 100644 --- a/docs/cn-loci.md +++ b/docs/cn-loci.md @@ -85,7 +85,7 @@ The call has two parts: The sites, each CYP21A2's base in the reference with CYP21A1P's base at the matching position (dbSNP positions, the CYP21A1P bases read from the reference, checked on 2026-11-04): In2G c.293-13A/C>G (rs6467, chr6:32,039,081; CYP21A2 carries A or C there, CYP21A1P G), I172N c.518T>A (rs6475), V237E c.713T>A (rs12530380, the exon 6 cluster), V281L c.844G>T (rs6471), Q318X c.955C>T (rs7755898) and R356W c.1069C>T (rs7769409). -A site is low when its estimate is more than 0.6 copies under the larger of the median site's count and 2. `cyp21a2_copies` is the least count over the low sites, or the median site's count when no site is low. The `finding` is `deletion` when a site is low and the module total is under 4 (the 30 kb chimeric deletion leaves one module fewer and a chimera with CYP21A1P's sequence at its 5' sites), `conversion` at the low sites when the module total is 4 or more, `pseudogene_deletion` when the module total is under 4 with no low site, and `none` otherwise. At 30x a site has about 70 reads, and read sampling alone gives an estimate of 2 copies there a standard deviation of about 0.25. On HG002's 35x reads the six sites read 1.62 to 2.82 copies, a wider spread, so the 0.6 margin is a screen setting, not a measured false-call rate. Which chromosome carries a change is not claimed, and a two-copy result does not rule out two copies on one chromosome and none on the other. The row has the MAPQ 0 guard of the HBA row: a second bedcov pass with `-Q 1`, and `mapq0_min` 0.02. HG002's 35x reads give a MAPQ 0 share of 0.099 over the two targets, from BWA-MEM and from minimap2 `-ax sr` alike (measured on 2026-11-08). An aligner that leaves out the reads that fit both genes (novoalign by default) gives about 0 and a lower module total estimate (4.22 to 3.81 on HG002's reads with every MAPQ 0 read removed), which on a genome whose depth runs lower reads as a deletion; under 0.02 the row is not assessed (`low_mapq0_share`), and its `message` row says the BAM has no multi-mapping reads there. The PSV sites still count every read whatever its mapping quality. Step 35 writes the table to `_cyp21a2_psv.tsv` ([step 35](35-paralogs.md#cyp21a2-from-depth-and-psv-sites)); the couple screen takes it with `--cyp21a2`. +A site is low when its estimate is more than 0.6 copies under the larger of the median site's count and 2. `cyp21a2_copies` is the least count over the low sites, or the median site's count when no site is low. The `finding` is `deletion` when a site is low and the module total is under 4 (the 30 kb chimeric deletion leaves one module fewer and a chimera with CYP21A1P's sequence at its 5' sites), `conversion` at the low sites when the module total is 4 or more, `pseudogene_deletion` when the module total is under 4 with no low site, and `none` otherwise. At 30x a site has about 70 reads, and read sampling alone gives an estimate of 2 copies there a standard deviation of about 0.25. On HG002's 35x reads the six sites read 1.62 to 2.82 copies, a wider spread, so the 0.6 margin is a screen setting, not a measured false-call rate. Which chromosome carries a change is not claimed, and a two-copy result does not rule out two copies on one chromosome and none on the other. The row has the MAPQ 0 guard of the HBA row: a second bedcov pass with `-Q 1`, and `mapq0_min` 0.02. HG002's 35x reads give a MAPQ 0 share of 0.099 over the two targets, from BWA-MEM and from minimap2 `-ax sr` alike (measured on 2026-11-08). An aligner that leaves out the reads that fit both genes (novoalign by default) gives about 0 and a lower module total estimate (4.22 to 3.81 on HG002's reads with every MAPQ 0 read removed), which on a genome whose depth runs lower reads as a deletion; under 0.02 the row is not assessed (`low_mapq0_share`), and its `message` row says the BAM has few or no multi-mapping reads there. The PSV sites still count every read whatever its mapping quality. Step 35 writes the table to `_cyp21a2_psv.tsv` ([step 35](35-paralogs.md#cyp21a2-from-depth-and-psv-sites)); the couple screen takes it with `--cyp21a2`. ### The interval rows diff --git a/docs/lessons-learned.md b/docs/lessons-learned.md index 8576ac66..6f286775 100644 --- a/docs/lessons-learned.md +++ b/docs/lessons-learned.md @@ -625,7 +625,7 @@ Most bioinformatics containers run as non-root users. If writing to bind-mounted ### A BAM without MAPQ 0 reads could read as a CYP21A2 deletion - **Failed:** the CYP21A2 row of `bin/cn_loci.tsv` added up the depth over CYP21A2 and CYP21A1P with no mapping quality check. The two are about 98% alike, so a read that fits both gets MAPQ 0. An aligner that leaves those reads out (novoalign does by default) lowers the module total. On HG002's reads, removing them takes the estimate from 4.22 to 3.81; on a genome whose depth runs lower it rounds to 3, which with a low PSV site reads as the 30 kb chimeric deletion: a false carrier finding for congenital adrenal hyperplasia. - **Why:** HBA had the guard (a second `samtools bedcov -Q 1` pass and `mapq0_min`); the CYP21A2 row was added later without it. On HG002's 35x reads, 9.9% of the depth over the two targets comes from MAPQ 0 reads, with BWA-MEM and with minimap2 `-ax sr` alike, so a share near 0 means the aligner dropped them. -- **Fix:** the row runs both bedcov passes with `mapq0_min` 0.02. Under that share it is not assessed (`low_mapq0_share`), and the table's `message` row, which both reports print, says the BAM has no multi-mapping reads there. `tests/test_copy_number.py` and the e2e case `cn-guard-1-mapq0` check it. +- **Fix:** the row runs both bedcov passes with `mapq0_min` 0.02. Under that share it is not assessed (`low_mapq0_share`), and the table's `message` row, which both reports print, says the BAM has few or no multi-mapping reads there. `tests/test_copy_number.py` and the e2e case `cn-guard-1-mapq0` check it. ## Downloads (2026-11) diff --git a/scripts/ci/e2e-case-times.tsv b/scripts/ci/e2e-case-times.tsv index 3ff5ba01..6595b720 100644 --- a/scripts/ci/e2e-case-times.tsv +++ b/scripts/ci/e2e-case-times.tsv @@ -72,7 +72,7 @@ clinvar-clip-check-1 shard 8 - clinvar-combined-allele-1 shard 25 - clinvar-cpsr-quality-1 shard 5 clinvar-clip-check-1 cn-engine-1-parity shard 33 - -cn-guard-1-mapq0 shard 15 - +cn-guard-1-mapq0 shard 5 - cn-intervals-1-catalog shard 0 - cn-intervals-2-hg002 shard 19 - cn-report-1-tables shard 60 - diff --git a/tests/e2e/cn-guard-1-mapq0.sh b/tests/e2e/cn-guard-1-mapq0.sh index 957bf473..e1c5ca2a 100755 --- a/tests/e2e/cn-guard-1-mapq0.sh +++ b/tests/e2e/cn-guard-1-mapq0.sh @@ -7,7 +7,7 @@ # - every MAPQ 0 read removed (samtools view -q 1), as an aligner that leaves # out the reads fitting both copies writes it: not assessed, # low_mapq0_share, no copies, no module total, no finding, and a message -# that says the BAM has no multi-mapping reads. +# that says the BAM has few or no multi-mapping reads. # Red control: the BAM without MAPQ 0 reads gives a lower estimate (the # table keeps it when not assessed), so the guard stops a real drift. # The block of step 35 runs on each (cut from scripts/35-paralogs.sh, with @@ -67,8 +67,8 @@ check "unmodified: never not assessed for the MAPQ 0 share" lacks 'low_mapq0_sha check_eq "without MAPQ 0 reads: not assessed (low_mapq0_share), share 0.000, no copies, module total or finding" \ "$(col "$N" status)/$(col "$N" reason)/$(col "$N" mapq0_share)/$(col "$N" cyp21a2_copies)/$(col "$N" module_total)/$(col "$N" finding)" \ "not assessed/low_mapq0_share/0.000/././." -check "without MAPQ 0 reads: the message says the BAM has no multi-mapping reads" \ - has '^not assessed: the BAM has no multi-mapping reads here' "$(col "$N" message)" +check "without MAPQ 0 reads: the message says the BAM has few or no multi-mapping reads" \ + has '^not assessed: the BAM has few or no multi-mapping reads here' "$(col "$N" message)" # Red control: the estimate the guard withheld. Removing the MAPQ 0 reads lowers # the module total estimate (4.22 to 3.81 on these reads, measured 2026-11-08), diff --git a/tests/test_copy_number.py b/tests/test_copy_number.py index b887c92a..9665db09 100644 --- a/tests/test_copy_number.py +++ b/tests/test_copy_number.py @@ -1023,7 +1023,7 @@ def psv(cat, read): rows, line = cn.psv_table(r, b) tab = dict(rows) check("CYP21A2 without MAPQ 0 reads: the table says so (message, share 0.000), no copies, module total or finding", - tab["message"].startswith("not assessed: the BAM has no multi-mapping reads here") + tab["message"].startswith("not assessed: the BAM has few or no multi-mapping reads here") and "CYP21A2 and CYP21A1P" in tab["message"] and tab["mapq0_share"] == "0.000" and (tab["cyp21a2_copies"], tab["module_total"], tab["finding"]) == (".", ".", ".") and "low_mapq0_share" in line, (tab, line)) @@ -1066,7 +1066,7 @@ def psv(cat, read): nm, nl = rd("no_mapq0"), rd("normal") check("step 35 CYP21A2 block on the deletion's reads without MAPQ 0 reads: low_mapq0_share, the message says why", (nm.get("reason"), nm.get("mapq0_share")) == ("low_mapq0_share", "0.000") - and nm.get("message", "").startswith("not assessed: the BAM has no multi-mapping reads here"), nm) + and nm.get("message", "").startswith("not assessed: the BAM has few or no multi-mapping reads here"), nm) check("step 35 CYP21A2 block, unmodified: the MAPQ >= 1 pass ran (share 0.100), message '.'", (nl.get("mapq0_share"), nl.get("message")) == ("0.100", "."), nl) check("step 35 CYP21A2 block with a failed mpileup: the rows of a called table but the sites'", @@ -1145,8 +1145,8 @@ def psv(cat, read): os.replace(tabs["no_mapq0"], os.path.join(sd, "paralogs", "S_cyp21a2_psv.tsv")) _, data = cs.sec_parascopy(sd, "S") got = render_report.cyp21a2_line(data.get("cyp21a2")) - check("report: the CYP21A2 line without MAPQ 0 reads is not assessed and says the BAM has no multi-mapping reads", - ": not assessed (the BAM has no multi-mapping reads here" in got and "deletion," not in got, got) + check("report: the CYP21A2 line without MAPQ 0 reads is not assessed and says the BAM has few or no multi-mapping reads", + ": not assessed (the BAM has few or no multi-mapping reads here" in got and "deletion," not in got, got) if __name__ == "__main__":