Skip to content

fix(cn): stop calling a CYP21A2 deletion when the BAM dropped its multi-mapping reads #63

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
geiresx wants to merge 4 commits into from
Closed

fix(cn): stop calling a CYP21A2 deletion when the BAM dropped its multi-mapping reads #63

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
fea0cb4
fix(cn): stop calling a CYP21A2 deletion when the BAM dropped its mul…
GeiserX Nov 8, 2026
cd2e89c
test(e2e): run the CYP21A2 MAPQ 0 guard on HG002's reads with the MAP…
GeiserX Nov 8, 2026
3583bec
fix(cn): say a low MAPQ 0 share means few or no multi-mapping reads, …
GeiserX Nov 8, 2026
6c6684d
train: #334 on main and the PRs queued before it
GeiserX Nov 8, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 5 additions & 1 deletion bin/cn_loci.tsv
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 .
Expand All @@ -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 .
Expand Down
20 changes: 19 additions & 1 deletion bin/copy_number.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 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 '.'}"


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}")
Expand All @@ -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
Expand Down
Loading
Loading