From fa610ddfab877c71d58873a85a6dc44a0c5d7f9a Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Wed, 4 Nov 2026 14:49:23 +0100 Subject: [PATCH 1/9] fix(cn): pass the interval table and CNVpytor calls to the couple screen and the Nextflow report Step 39 and COUPLE_CARRIER now hand couple_report.py each partner's step 35 interval table and step 18 CNVpytor VCF when they exist, so the carrier deletion alleles and AZF lines appear without a hand-run command. HTML_REPORT places _cn_intervals.tsv in paralogs/, where collect_summary.py reads it, and main.nf joins it into the report inputs. Bead: pgp-pfw6.36 --- docs/35-paralogs.md | 1 + docs/39-couple.md | 2 +- docs/cn-loci.md | 2 +- main.nf | 4 +- modules/local/couple/main.nf | 11 +- modules/local/html_report/main.nf | 1 + scripts/39-couple.sh | 8 +- tests/e2e/followups-cn-1-couple-report.sh | 225 ++++++++++++++++++++++ workflows/bam_analysis.nf | 3 + workflows/family.nf | 5 +- workflows/sv.nf | 1 + 11 files changed, 252 insertions(+), 11 deletions(-) create mode 100755 tests/e2e/followups-cn-1-couple-report.sh diff --git a/docs/35-paralogs.md b/docs/35-paralogs.md index 36f62cef..d9f54f9c 100644 --- a/docs/35-paralogs.md +++ b/docs/35-paralogs.md @@ -71,6 +71,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: `) and `.` in every other column | | `${SAMPLE}__copy_number.tsv` | The same, for every other catalog row (none in the catalog today) | | `${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}_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)) | ## Runtime diff --git a/docs/39-couple.md b/docs/39-couple.md index dbe6e0cf..d7952c3f 100644 --- a/docs/39-couple.md +++ b/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/_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/_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/_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/_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 passes both when a partner's folder holds them (`paralogs/_cn_intervals.tsv` of step 35, `cnvpytor/_cnvs.vcf.gz` of step 18), and the Nextflow couple screen passes them when `--tools` has `parascopy` and `cnvpytor`. ### The carrier list diff --git a/docs/cn-loci.md b/docs/cn-loci.md index 5e395313..300b41fa 100644 --- a/docs/cn-loci.md +++ b/docs/cn-loci.md @@ -83,7 +83,7 @@ The report does not show `_cn_dosage.tsv` yet; the table and the step lo A row on chrX or chrY is read against the sample's own count there: its controls sit on the same chromosome, so a male's one X or Y gives 1 copy, and `expected` says so. The sex is the declared one: step 35's second argument, and in Nextflow the samplesheet's sex column. Step 16's inferred sex is not used, so both give the same table. A chrY row on a female is `not applicable`; without a sex, chrX and chrY rows are `not assessed` (`sex_unknown`). Every row needs 8x of control depth (4x on chrY), flanks within 1.5x of each other, and an estimate within 0.3 of a whole number; a control with no read makes that row `not assessed` (`broken_control`) and the others still run. A failed bedcov pass writes every row as `not assessed` (`check_failed`). -The couple screen (step 39) reads the table with `couple_report.py screen --cn-intervals ID=_cn_intervals.tsv`: each carrier allele is present, absent or not assessable for each partner, and a male partner's AZF deletions are listed the same way (present means the deletion is there). With `--cnvpytor ID=_cnvs.vcf.gz` (step 18), a CNVpytor call that disagrees with the depth call is printed beside it. Step 39 and the Nextflow couple screen do not pass these two inputs yet. +The couple screen (step 39) reads the table with `couple_report.py screen --cn-intervals ID=_cn_intervals.tsv`: each carrier allele is present, absent or not assessable for each partner, and a male partner's AZF deletions are listed the same way (present means the deletion is there). With `--cnvpytor ID=_cnvs.vcf.gz` (step 18), a CNVpytor call that disagrees with the depth call is printed beside it. Step 39 and the Nextflow couple screen pass both when the partner has them. To change a row, change the script (its hand rows or a pinned download) and run `scripts/ci/build-cn-intervals.py`; `--check` fails when the catalog is not what the downloads give. diff --git a/main.nf b/main.nf index 40e447c4..db09c89b 100644 --- a/main.nf +++ b/main.nf @@ -853,7 +853,8 @@ workflow { "the couple carrier screen reads each partner's ClinVar hits, so it does not run." } def ch_couple_outputs = couple_ids && tools_list.contains('clinvar') ? Channel.empty().mix( - PGX.out.clinvar_dir, SV.out.delly_vcf, BAM_ANALYSIS.out.expansion_vcf, BAM_ANALYSIS.out.smn_copy_number + PGX.out.clinvar_dir, SV.out.delly_vcf, BAM_ANALYSIS.out.expansion_vcf, BAM_ANALYSIS.out.smn_copy_number, + BAM_ANALYSIS.out.cn_intervals, SV.out.cnvpytor_vcf ).filter { meta, f -> couple_ids.contains(meta.id) }.map { meta, f -> f } : Channel.empty() FAMILY( ch_gvcf.filter { meta, g, idx -> family_ids.contains(meta.id) }, @@ -1118,6 +1119,7 @@ workflow { .join(CLINICAL.out.ancestry_fine_png.map { meta, f -> [meta.id, f] }, remainder: true) .join(BAM_ANALYSIS.out.coverage.map { meta, f -> [meta.id, f] }, remainder: true) .join(BAM_ANALYSIS.out.callability.map { meta, f -> [meta.id, f] }, remainder: true) + .join(BAM_ANALYSIS.out.cn_intervals.map { meta, f -> [meta.id, f] }, remainder: true) .join(BAM_ANALYSIS.out.sample_qc.map { meta, f -> [meta.id, f] }, remainder: true) .join(BAM_ANALYSIS.out.telomere_results.map { meta, f -> [meta.id, f] }, remainder: true) .join(BAM_ANALYSIS.out.mito_vcf.map { meta, f -> [meta.id, f] }, remainder: true) diff --git a/modules/local/couple/main.nf b/modules/local/couple/main.nf index d85bb5b0..7c4e0edd 100644 --- a/modules/local/couple/main.nf +++ b/modules/local/couple/main.nf @@ -8,15 +8,17 @@ likely pathogenic hits (CLINVAR_SCREEN), Delly deletions over an exon of 50 bp or more (with --gencode_exons; one matching a common deletion of AnnotSV's benign loss set is set aside, the set read from - --annotsv_annotations when it is there), SMN1 copy number (PARASCOPY) and - the FMR1 and FXN repeats (EXPANSION_HUNTER). It writes couple__.tsv and .md, published + --annotsv_annotations when it is there), SMN1 copy number (PARASCOPY), + the FMR1 and FXN repeats (EXPANSION_HUNTER), and the carrier deletion + alleles and a male partner's AZF regions from PARASCOPY's interval table, + with CNVpytor's calls (CNVPYTOR_VCF) as their cross-check. It writes couple__.tsv and .md, published under the family's folder, // (--family_name, default family/), and one copy of the table per partner, _couple_carrier.tsv, published under //couple/ for that partner's report. The step outputs come as one list; each file is matched to its sample by its name (_clinvar_hits.vcf, _sv.vcf.gz, _eh.vcf, - _smn_copy_number.tsv). + _smn_copy_number.tsv, _cn_intervals.tsv, _cnvs.vcf.gz). Equivalent to: scripts/39-couple.sh ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -51,7 +53,8 @@ process COUPLE_CARRIER { script: def files = outputs instanceof List ? outputs : [outputs] - def kinds = [clinvar: '_clinvar_hits.vcf', delly: '_sv.vcf.gz', eh: '_eh.vcf', smn: '_smn_copy_number.tsv'] + def kinds = [clinvar: '_clinvar_hits.vcf', delly: '_sv.vcf.gz', eh: '_eh.vcf', smn: '_smn_copy_number.tsv', + 'cn-intervals': '_cn_intervals.tsv', cnvpytor: '_cnvs.vcf.gz'] def args = files.collect { f -> // A staged input's name is its stage path (in/): keep the file name only def n = f.name.tokenize('/')[-1] diff --git a/modules/local/html_report/main.nf b/modules/local/html_report/main.nf index 11acd25c..1b4f2e4e 100644 --- a/modules/local/html_report/main.nf +++ b/modules/local/html_report/main.nf @@ -72,6 +72,7 @@ process HTML_REPORT { ${id}_paraphase_copy_number.tsv|${id}_paraphase_variants.tsv|${id}_smn_copy_number.tsv|${id}_hba_copy_number.tsv) place paralogs "\$f" ;; ${id}.mosdepth.summary.txt|${id}_callability.bed.gz) place coverage "\$f" ;; ${id}_sample_qc.tsv) place qc "\$f" ;; + ${id}_cn_intervals.tsv) place paralogs "\$f" ;; ${id}_prs_summary.tsv) place prs "\$f" ;; ${id}_ancestry.tsv) place ancestry "\$f" ;; ${id}_ancestry_fine.tsv|${id}_painting.png) place ancestry_fine "\$f" ;; diff --git a/scripts/39-couple.sh b/scripts/39-couple.sh index 4b928d3e..8e2cfc50 100755 --- a/scripts/39-couple.sh +++ b/scripts/39-couple.sh @@ -10,7 +10,10 @@ # ($GENOME_DIR//clinvar/_clinvar_hits.vcf). Read when # they are there: Delly's calls (step 19, /delly/_sv.vcf.gz), # ExpansionHunter's (step 09, /expansion_hunter/_eh.vcf), -# Parascopy's SMN1 table (step 35, /paralogs/_smn_copy_number.tsv) +# Parascopy's SMN1 table (step 35, /paralogs/_smn_copy_number.tsv), +# step 35's interval table (/paralogs/_cn_intervals.tsv: +# the carrier deletion alleles and a male partner's AZF regions) and +# CNVpytor's calls, its cross-check (step 18, /cnvpytor/_cnvs.vcf.gz), # GENCODE's exons (setup.sh: reference/gencode.v.basic.exons.bed) # and AnnotSV's benign loss set (setup.sh --annotsv-data, with gnomAD-SV's # common deletions): a Delly deletion matching a common one is set aside. @@ -115,7 +118,8 @@ for s in ${NAMED[@]+"${NAMED[@]}"}; do fi [ -d "${GENOME_DIR}/${s}" ] || continue for kind in "clinvar:clinvar/${s}_clinvar_hits.vcf:06" "delly:delly/${s}_sv.vcf.gz:19" \ - "eh:expansion_hunter/${s}_eh.vcf:09" "smn:paralogs/${s}_smn_copy_number.tsv:35"; do + "eh:expansion_hunter/${s}_eh.vcf:09" "smn:paralogs/${s}_smn_copy_number.tsv:35" \ + "cn-intervals:paralogs/${s}_cn_intervals.tsv:35" "cnvpytor:cnvpytor/${s}_cnvs.vcf.gz:18"; do IFS=: read -r k rel step <<< "$kind" if [ -f "${GENOME_DIR}/${s}/${rel}" ]; then ARGS+=("--${k}" "${s}=/genome/${s}/${rel}") diff --git a/tests/e2e/followups-cn-1-couple-report.sh b/tests/e2e/followups-cn-1-couple-report.sh new file mode 100755 index 00000000..59570fd8 --- /dev/null +++ b/tests/e2e/followups-cn-1-couple-report.sh @@ -0,0 +1,225 @@ +#!/usr/bin/env bash +# Step 35's interval table (_cn_intervals.tsv) and step 18's CNVpytor +# calls reach the couple screen and the Nextflow report. +# 1. step 39 on a made-up couple, FCN_A (male) and FCN_B (female), each +# with an empty ClinVar screen, an interval table and a CNVpytor VCF in +# their folders: the couple report lists FCN_A's del(GJB6-D13S1830) and +# AZFc deletion as present, FCN_B's CFTRdele2,3 as absent with the +# CNVpytor call beside it, and FCN_B's AZF line as not applicable. +# Red control: with FCN_A's interval table gone the report has no line +# of FCN_A's alleles and names the missing table. +# 2. COUPLE_CARRIER, run for real in its image through a two-line pipeline +# on the same files: the same couple table and report as step 39. +# Red control: a copy of the module without the two input kinds writes +# no allele line. +# 3. HTML_REPORT, run for real through a two-line pipeline with the SMN +# table and the interval table: the task places the interval table and +# summary.json carries its rows with a change. Red control: a copy of +# the module without that place line stops on the input. +# 4. main.nf passes both files on: a stub run of the stub pedigree (two +# partners) with clinvar, parascopy, cnvpytor and html_report stages +# each partner's interval table and CNVpytor VCF in COUPLE_CARRIER and +# the interval table in each HTML_REPORT task. Red control: a copy of +# the pipeline without the two lines stages none. +. "$(dirname "$0")/lib.sh" + +command -v nextflow >/dev/null || { fail "nextflow is not on PATH"; finish; } + +VCF="${GENOME_DIR}/${SAMPLE}/vcf/${SAMPLE}.vcf.gz" +check "an earlier case left ${VCF#"$GENOME_DIR"/}" nonempty "${VCF#"$GENOME_DIR"/}" + +C="${E2E_WORK}/genome-followups-cn" +docker run --rm -v "${E2E_WORK}:/w" "$PYTHON_IMAGE" rm -rf /w/genome-followups-cn +A=FCN_A; B=FCN_B +for s in "$A" "$B"; do mkdir -p "${C}/${s}/clinvar" "${C}/${s}/paralogs" "${C}/${s}/cnvpytor"; done +G0="$GENOME_DIR" +export GENOME_DIR="$C" + +ivh() { printf 'locus\tname\tchrom\tstart\tend\tstatus\tcopies\texpected\tchange\testimate\treason\n'; } +iv() { # iv LOCUS CHROM START END STATUS COPIES EXPECTED CHANGE REASON + printf '%s\t%s name\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t.\t%s\n' "$1" "$1" "$2" "$3" "$4" "$5" "$6" "$7" "$8" "$9" +} +for s in "$A" "$B"; do + printf '##fileformat=VCFv4.2\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t%s\n' "$s" \ + > "${C}/${s}/clinvar/${s}_clinvar_hits.vcf" +done +{ + ivh + iv GJB6_D13S1830 chr13 20223037 20531806 called 1 2 deletion . + iv CFTR_dele2_3 chr7 117498312 117519392 called 2 2 "no change" . + iv HBB_exon3_del chr11 5225254 5225875 called 2 2 "no change" . + iv AZFa chrY 12537649 12920478 called 1 1 "no change" . + iv AZFb chrY 19702297 20593154 called 1 1 "no change" . + iv AZFc_DAZ chrY 23129354 24907040 called 2 4 deletion . +} > "${C}/${A}/paralogs/${A}_cn_intervals.tsv" +{ + ivh + iv GJB6_D13S1830 chr13 20223037 20531806 called 2 2 "no change" . + iv CFTR_dele2_3 chr7 117498312 117519392 called 2 2 "no change" . + iv HBB_exon3_del chr11 5225254 5225875 called 2 2 "no change" . + iv AZFa chrY 12537649 12920478 "not applicable" . . . no_y +} > "${C}/${B}/paralogs/${B}_cn_intervals.tsv" +# CNVpytor: FCN_A has no call, FCN_B a deletion over CFTRdele2,3 (absent by depth) +printf '##fileformat=VCFv4.2\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n' | gzip -n > "${C}/${A}/cnvpytor/${A}_cnvs.vcf.gz" +printf '##fileformat=VCFv4.2\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\nchr7\t117498001\tCNV1\tN\t\t.\tPASS\tSVTYPE=DEL;END=117520000\n' \ + | gzip -n > "${C}/${B}/cnvpytor/${B}_cnvs.vcf.gz" +printf 'sample\tsex\tfather\tmother\tpartner\n%s\tmale\t.\t.\t%s\n%s\tfemale\t.\t.\t%s\n' "$A" "$B" "$B" "$A" > "${C}/pedigree.tsv" + +GJB6_A="- ${A}, del(GJB6-D13S1830): \*\*present\*\* \(read depth, step 35: 1 of 2 copies\); CNVpytor \(step 18\) calls no deletion here" +AZFC_A="- ${A}, AZFc deletion \(DAZ copies; gr/gr and b2/b3 leave two of four\): \*\*present\*\* \(read depth, step 35: 2 of 4 copies\)" +CFTR_B="- ${B}, CFTRdele2,3: \*\*absent\*\* \(read depth, step 35: 2 of 2 copies\); CNVpytor \(step 18\) calls a deletion here" +AZF_B="- ${B}, AZF regions: not applicable \(no Y chromosome\)" + +# --- 1. step 39 ------------------------------------------------------------------- +run_step 39-couple.sh couple "${C}/pedigree.tsv" +check_step_exit 39-couple.sh +check "step 39 names both inputs of each partner" \ + sh -c "grep -q '${A}: paralogs/${A}_cn_intervals.tsv (step 35)' '$STEP_LOG' && grep -q '${B}: cnvpytor/${B}_cnvs.vcf.gz (step 18)' '$STEP_LOG'" +MD="${C}/couple/couple_${A}_${B}.md" +TXT=$(cat "$MD" 2>/dev/null) +grep -E "^- (${A}|${B}), " <<< "$TXT" +check "step 39: ${A}'s del(GJB6-D13S1830) present, CNVpytor disagreeing" has "$GJB6_A" "$TXT" +check "step 39: ${A}'s AZFc deletion present" has "$AZFC_A" "$TXT" +check "step 39: ${B}'s CFTRdele2,3 absent, the CNVpytor deletion beside it" has "$CFTR_B" "$TXT" +check "step 39: ${B}'s AZF line not applicable" has "$AZF_B" "$TXT" +cp "$MD" "${CASE_TMP}/step.md" +cp "${C}/couple/couple_${A}_${B}.tsv" "${CASE_TMP}/step.tsv" + +mv "${C}/${A}/paralogs/${A}_cn_intervals.tsv" "${CASE_TMP}/a_iv.tsv" +run_step 39-couple.sh couple "${C}/pedigree.tsv" +check_step_exit 39-couple.sh +TXT0=$(cat "$MD" 2>/dev/null) +check "red: without ${A}'s interval table no line of ${A}'s alleles" lacks "^- ${A}, (del\(GJB6|AZFc)" "$TXT0" +check "red: and the report names the missing table" has "Carrier deletion alleles of ${A}: no interval table \(step 35\)" "$TXT0" +mv "${CASE_TMP}/a_iv.tsv" "${C}/${A}/paralogs/${A}_cn_intervals.tsv" + +# --- 2. COUPLE_CARRIER -------------------------------------------------------------- +NP="${CASE_TMP}/project" +rm -rf "$NP"; mkdir -p "$NP" +cp -R "${REPO}/bin" "${NP}/bin" +mkdir -p "${NP}/couple-nokinds" +sed -e "/'cn-intervals': '_cn_intervals.tsv', cnvpytor: '_cnvs.vcf.gz'/d" \ + -e "s|smn: '_smn_copy_number.tsv',\$|smn: '_smn_copy_number.tsv']|" "${REPO}/modules/local/couple/main.nf" \ + > "${NP}/couple-nokinds/main.nf" +check_eq "control copy of the couple module: the two kinds are gone" \ + "$(grep -c "_cn_intervals.tsv'" "${NP}/couple-nokinds/main.nf")" 0 +nfconf() { # nfconf OUT + printf "params.outdir = '%s'\nparams.publish_dir_mode = 'copy'\nparams.family_name = 'couple'\ndocker.enabled = true\ndocker.runOptions = '--user %s:%s --network none'\n" \ + "$1" "$(id -u)" "$(id -g)" > "${NP}/nextflow.config" +} +nfrun() { # nfrun NAME: the pipeline in NP/main.nf, log CASE_TMP/NAME.log + echo "+ nextflow run ($1)" + ( cd "$CASE_TMP" && nextflow run "${NP}/main.nf" -ansi-log false -c "${REPO}/conf/containers.config" \ + -work-dir "${CASE_TMP}/work-$1" ) > "${CASE_TMP}/$1.log" 2>&1 + NF_RC=$? + [ "$NF_RC" = 0 ] || grep -E 'ERROR|Error executing process|Command error' -A 8 "${CASE_TMP}/$1.log" | head -40 +} +couple_run() { # couple_run NAME MODULE + rm -rf "${CASE_TMP:?}/$1" + nfconf "${CASE_TMP}/$1" + cat > "${NP}/main.nf" </dev/null)" +couple_run couple-nokinds "${NP}/couple-nokinds/main" +check "red: the module without the two kinds writes no allele line" \ + lacks "^- ${A}, del\(GJB6" "$(cat "${CASE_TMP}/couple-nokinds/couple/couple_${A}_${B}.md" 2>/dev/null)" + +# --- 3. HTML_REPORT ----------------------------------------------------------------- +T="${CASE_TMP}/tables" +rm -rf "$T"; mkdir -p "$T" +printf 'chrom\tstart\tend\tlocus\tagCN_filter\tagCN\tagCN_qual\tpsCN_filter\tpsCN\tpsCN_qual\thomologous_regions\n' \ + > "${T}/${SAMPLE}_smn_copy_number.tsv" +cp "${C}/${A}/paralogs/${A}_cn_intervals.tsv" "${T}/${SAMPLE}_cn_intervals.tsv" +mkdir -p "${NP}/report-noplace" +sed '/_cn_intervals.tsv) *place paralogs/d' "${REPO}/modules/local/html_report/main.nf" > "${NP}/report-noplace/main.nf" +check_eq "control copy of the report module: the place line is gone" \ + "$(grep -c '_cn_intervals.tsv)' "${NP}/report-noplace/main.nf")" 0 +report_run() { # report_run NAME MODULE + rm -rf "${CASE_TMP:?}/$1" + nfconf "${CASE_TMP}/$1" + cat > "${NP}/main.nf" </dev/null \ + | while read -r f; do grep -lE -- "$3" "$f"; done | wc -l | tr -d ' ' +} +for n in stub stub-nowire; do + rm -rf "${CASE_TMP:?}/${n}"; mkdir -p "${CASE_TMP}/${n}" + git -C "$REPO" ls-files -z | (cd "$REPO" && tar --null -T - -cf -) | tar -xf - -C "${CASE_TMP}/${n}" +done +W1='BAM_ANALYSIS.out.cn_intervals.map'; W2='BAM_ANALYSIS.out.cn_intervals, SV.out.cnvpytor_vcf' +check_eq "main.nf joins the interval table into the report inputs once" "$(grep -c "$W1" "${CASE_TMP}/stub/main.nf")" 1 +check_eq "main.nf hands both files to the couple screen once" "$(grep -c "$W2" "${CASE_TMP}/stub/main.nf")" 1 +sed -i -e "/${W1}/d" -e "s/BAM_ANALYSIS.out.smn_copy_number,\$/BAM_ANALYSIS.out.smn_copy_number/" -e "/${W2}/d" \ + "${CASE_TMP}/stub-nowire/main.nf" +check_eq "control copy: both lines are gone" "$(grep -c -e "$W1" -e "$W2" "${CASE_TMP}/stub-nowire/main.nf")" 0 +echo '{"kir": false}' > "${CASE_TMP}/no-kir.json" +for n in stub stub-nowire; do + echo "+ nextflow run main.nf -profile test_all -stub --pedigree assets/stub/pedigree.tsv (${n})" + ( cd "${CASE_TMP}/${n}" && nextflow run main.nf -profile test_all -stub -ansi-log false --outdir "${CASE_TMP}/${n}/out" \ + -work-dir "${CASE_TMP}/${n}/work" --pedigree assets/stub/pedigree.tsv --tools clinvar,parascopy,cnvpytor,html_report \ + -params-file "${CASE_TMP}/no-kir.json" ) > "${CASE_TMP}/${n}.log" 2>&1 + rc=$? + tail -n 6 "${CASE_TMP}/${n}.log" + check_eq "stub run (${n}) exits 0" "$rc" 0 +done +S="${CASE_TMP}/stub/work"; S0="${CASE_TMP}/stub-nowire/work" +check_eq "COUPLE_CARRIER stages the interval table" "$(staged_in "$S" COUPLE_CARRIER '_cn_intervals.tsv')" 1 +check_eq "COUPLE_CARRIER stages the CNVpytor VCF" "$(staged_in "$S" COUPLE_CARRIER '_cnvs.vcf.gz')" 1 +check_eq "both partners' tables reach it" \ + "$(find "$S" -name .command.run -exec grep -l '^### name: .*COUPLE_CARRIER' {} + 2>/dev/null | while read -r f; do grep -ho 'stub_[a-z]*_cn_intervals.tsv' "$f"; done | sort -u | wc -l | tr -d ' ')" 2 +check_eq "every HTML_REPORT task stages the interval table" \ + "$(staged_in "$S" HTML_REPORT '_cn_intervals.tsv')" "$(staged_in "$S" HTML_REPORT '_smn_copy_number.tsv')" +check_ge "HTML_REPORT tasks with the SMN table" "$(staged_in "$S" HTML_REPORT '_smn_copy_number.tsv')" 1 +check_ge "control: the stub run without the lines still ran COUPLE_CARRIER" "$(staged_in "$S0" COUPLE_CARRIER '_clinvar_hits.vcf')" 1 +check_eq "red: without the lines COUPLE_CARRIER stages neither file" \ + "$(staged_in "$S0" COUPLE_CARRIER '_cn_intervals[.]tsv|_cnvs[.]vcf[.]gz')" 0 +check_eq "red: and no HTML_REPORT task stages the interval table" "$(staged_in "$S0" HTML_REPORT '_cn_intervals.tsv')" 0 + +rm -rf "${CASE_TMP:?}/stub" "${CASE_TMP:?}/stub-nowire" "${CASE_TMP:?}"/work-* +docker run --rm -v "${E2E_WORK}:/w" "$PYTHON_IMAGE" rm -rf /w/genome-followups-cn +finish diff --git a/workflows/bam_analysis.nf b/workflows/bam_analysis.nf index 135f44da..2406ea7d 100644 --- a/workflows/bam_analysis.nf +++ b/workflows/bam_analysis.nf @@ -75,6 +75,7 @@ workflow BAM_ANALYSIS { ch_mito_annotated = Channel.empty() // [meta, _chrM_annotated.tsv] ch_mito_deletions = Channel.empty() // [meta, _chrM_deletions.tsv] ch_mito_copies = Channel.empty() // [meta, _mtdna_copies.tsv] + ch_cn_intervals = Channel.empty() // [meta, _cn_intervals.tsv] when Parascopy runs ch_kir_genotype = Channel.empty() ch_hla_optitype = Channel.empty() // [meta, _optitype_result.tsv] with --hla_qc ch_hla_concordance = Channel.empty() // [meta, _hla_concordance.tsv] when OptiType typed reads @@ -234,6 +235,7 @@ workflow BAM_ANALYSIS { // if (params.tools && params.tools.split(',').collect{it.trim()}.contains('parascopy')) { PARASCOPY(ch_bam, ch_reference, ch_reference_fai, ch_parascopy_data, ch_parascopy_bed) + ch_cn_intervals = PARASCOPY.out.intervals ch_smn_copy_number = PARASCOPY.out.copy_number ch_hba_copy_number = PARASCOPY.out.hba ch_versions = ch_versions.mix(PARASCOPY.out.versions) @@ -377,6 +379,7 @@ workflow BAM_ANALYSIS { telomere_results = ch_telomere_results coverage = ch_coverage callability = ch_callability // [meta, _callability.bed.gz] + cn_intervals = ch_cn_intervals // [meta, _cn_intervals.tsv]: bin/cn_loci.tsv's interval rows mito_vcf = ch_mito_vcf mito_annotated = ch_mito_annotated // [meta, _chrM_annotated.tsv] mito_deletions = ch_mito_deletions // [meta, _chrM_deletions.tsv] diff --git a/workflows/family.nf b/workflows/family.nf index 413c1c88..9e31a5fa 100644 --- a/workflows/family.nf +++ b/workflows/family.nf @@ -54,8 +54,9 @@ workflow FAMILY { ch_reference // channel: val(path) — reference FASTA ch_reference_fai // channel: val(path) — reference .fai ch_par_bed // channel: val(path) — assets/par_grch38.bed - ch_couple_outputs // channel: path — the partners' ClinVar hits, Delly VCF, ExpansionHunter VCF and SMN1 - // table, named _*; empty when the pedigree marks no partners + ch_couple_outputs // channel: path — the partners' ClinVar hits, Delly VCF, ExpansionHunter VCF, SMN1 + // table, step 35's interval table and CNVpytor VCF, named _*; empty when the + // pedigree marks no partners ch_carrier_genes // channel: val(path) — assets/carrier_genes.tsv ch_exons // channel: val(path) — GENCODE exons BED (--gencode_exons), or [] ch_ibd_panel // channel: val(path) — step 26b's panel folder (--ancestry_fine_panel); empty without it, diff --git a/workflows/sv.nf b/workflows/sv.nf index 4b8480ae..59b70e3b 100644 --- a/workflows/sv.nf +++ b/workflows/sv.nf @@ -152,6 +152,7 @@ workflow SV { manta_vcf = ch_manta_vcf delly_vcf = ch_delly_vcf cnvpytor_calls = ch_cnvpytor_calls + cnvpytor_vcf = ch_cnvpytor_vcf // [meta, _cnvs.vcf.gz] duphold_vcf = ch_duphold_vcf annotsv_tsv = ch_annotsv_tsv annotsv_class4_5 = ch_annotsv_checked // [meta, _sv_class4_5_checked.tsv] From 46a2bf89bca67c5648b4887fa63c835f231c78ad Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Wed, 4 Nov 2026 14:51:50 +0100 Subject: [PATCH 2/9] test(cn): assert HG002's AZFc_DAZ reading on real reads HG002's DAZ depth gives an estimate of 3.63 in CI run 41767270253, so the row is not assessed with reason between_counts. The case now checks that, with the estimate held to 3.3 to 3.9, and checks that the same assertion fails on the female table. Bead: pgp-pfw6.36 --- docs/cn-loci.md | 2 +- tests/e2e/cn-intervals-2-hg002.sh | 16 ++++++++++++++-- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/docs/cn-loci.md b/docs/cn-loci.md index 300b41fa..c54652e8 100644 --- a/docs/cn-loci.md +++ b/docs/cn-loci.md @@ -110,7 +110,7 @@ pypgx owns the genes it calls from the BAM, the `bam` rows of `assets/pgx_genes. ## How it is tested - `tests/test_copy_number.py` section 9 runs the converter on a small ClinGen BED whose chrX rows are spelled `chrx` (the rows load on chrX; the same rows left as `chrx` are refused; an unknown contig stops the build). Made-up depths call a 22q11.2 A-D deletion at 1 copy and a 17p12 duplication at 3, and 2 for both unaltered. AZFc at half its DAZ depth on a male is a deletion; every AZF row on a female is not applicable. No interval row's control overlaps another row's target. The interval block of step 35, cut out and run on made-up bedcov tables, writes one line per row, and a failed bedcov writes every row as not assessed. `tests/test_couple_report.py` section 11 checks the three states and the CNVpytor cross-check. -- The e2e case `cn-intervals-1-catalog` rebuilds the interval rows from the pinned downloads and fails on any difference with the committed catalog. `cn-intervals-2-hg002` runs the interval blocks of step 35 and of the PARASCOPY module on HG002's real reads at the AZF genes, their controls and 22q11.2: a male with one copy of AZFa and AZFb, no change at 22q11.2 A-B, no call on a row off chrY and chr22 (they have no read; the chr22 rows share A-B's control window, so they are left out of that check), and the AZF rows not applicable when the same reads are read as female. +- The e2e case `cn-intervals-1-catalog` rebuilds the interval rows from the pinned downloads and fails on any difference with the committed catalog. `cn-intervals-2-hg002` runs the interval blocks of step 35 and of the PARASCOPY module on HG002's real reads at the AZF genes, their controls and 22q11.2: a male with one copy of AZFa and AZFb, AZFc_DAZ not assessed (its estimate, between 3.3 and 3.9, is no whole count of the four DAZ copies, so the row says `between_counts` rather than a call), no change at 22q11.2 A-B, no call on a row off chrY and chr22 (they have no read; the chr22 rows share A-B's control window, so they are left out of that check), and the AZF rows not applicable when the same reads are read as female. - `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`. diff --git a/tests/e2e/cn-intervals-2-hg002.sh b/tests/e2e/cn-intervals-2-hg002.sh index 65b6bfb1..457cd324 100755 --- a/tests/e2e/cn-intervals-2-hg002.sh +++ b/tests/e2e/cn-intervals-2-hg002.sh @@ -4,14 +4,17 @@ # of four rows (AZFa, AZFb, AZFc_DAZ and 22q11.2 A-B, targets and controls). # - the interval block of step 35 (cut from scripts/35-paralogs.sh, run # with scripts/lib/common.sh's run_in) on the sample read as male: AZFa -# and AZFb 1 copy of 1, 22q11.2 A-B 2 of 2, no change; every row off +# and AZFb 1 copy of 1, 22q11.2 A-B 2 of 2, no change; AZFc_DAZ not +# assessed (between_counts: HG002's DAZ depth gives an estimate between +# 3.3 and 3.9, no whole count of the four copies); every row off # chrY and chr22 has no read and is not assessed, never a call (the # chr22 rows share A-B's fetched control window, so they are left out); # - the interval block of modules/local/parascopy/main.nf, cut from the # file, in PARASCOPY_IMAGE on the same BAM: the same table, byte for byte; # - the same reads read as female: the three AZF rows not applicable, # 22q11.2 A-B still 2 copies. -# Control: the male table and the female table differ. +# Control: the male table and the female table differ, and the AZFc check +# fails on the female table. . "$(dirname "$0")/lib.sh" SRC="https://storage.googleapis.com/deepvariant/case-study-testdata/HG002.novaseq.pcr-free.35x.dedup.grch38_no_alt.bam" @@ -62,6 +65,13 @@ for r in AZFa AZFb; do check_eq "male: ${r} 1 copy of 1, no change" "$(col "$T" "$r" copies)/$(col "$T" "$r" expected)/$(col "$T" "$r" change)" \ "1/1/no change" done +# azfc TABLE: "status/reason/in" when AZFc_DAZ's estimate is in 3.3 to 3.9, else ".../out" +azfc() { + echo "$(col "$1" AZFc_DAZ status)/$(col "$1" AZFc_DAZ reason)/$(col "$1" AZFc_DAZ estimate \ + | awk '{print ($1 ~ /^[0-9.]+$/ && $1 >= 3.3 && $1 <= 3.9) ? "in" : "out"}')" +} +check_eq "male: AZFc_DAZ not assessed, between two whole counts, estimate 3.3 to 3.9 ($(col "$T" AZFc_DAZ estimate))" \ + "$(azfc "$T")" "not assessed/between_counts/in" check_eq "male: 22q11.2 A-B 2 copies of 2, no change" \ "$(col "$T" 22q11.2_DGS_proximal_A-B copies)/$(col "$T" 22q11.2_DGS_proximal_A-B expected)/$(col "$T" 22q11.2_DGS_proximal_A-B change)" \ "2/2/no change" @@ -90,6 +100,8 @@ for r in AZFa AZFb AZFc_DAZ; do done check_eq "female: 22q11.2 A-B still 2 copies" "$(col "$T" 22q11.2_DGS_proximal_A-B copies)" 2 check "control: the male and the female tables differ" bash -c '! cmp -s "$1" "$2"' _ "$T" "${CASE_TMP}/male.tsv" +check "control: the AZFc check fails on the female table ($(azfc "$T"))" \ + test "$(azfc "$T")" != "not assessed/between_counts/in" rm -rf "$W" finish From 8b505f3d615975e7edb70683abf37a5ee53fd220 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Wed, 4 Nov 2026 14:58:15 +0100 Subject: [PATCH 3/9] feat(cn): read the HBB 619 bp and Sicilian delta-beta deletions from depth Two interval rows at their HbVar GRCh38 spans (979 and 1035), checked by placing the NG_000007.3 breakpoint flanks on NC_000011.10. The couple screen lists both as carrier deletion alleles. The step 35 interval table gains two lines; no existing row or column changes. Bead: pgp-pfw6.36 --- bin/cn_loci.tsv | 2 ++ bin/couple_report.py | 6 ++++-- docs/39-couple.md | 2 +- docs/cn-loci.md | 6 +++--- scripts/ci/build-cn-intervals.py | 14 ++++++++++++-- tests/test_copy_number.py | 12 +++++++++--- tests/test_couple_report.py | 18 ++++++++++++++---- 7 files changed, 45 insertions(+), 15 deletions(-) diff --git a/bin/cn_loci.tsv b/bin/cn_loci.tsv index 9e8898bf..fde09f0d 100644 --- a/bin/cn_loci.tsv +++ b/bin/cn_loci.tsv @@ -135,6 +135,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 diff --git a/bin/couple_report.py b/bin/couple_report.py index 7cd2d015..82dc85e2 100755 --- a/bin/couple_report.py +++ b/bin/couple_report.py @@ -19,7 +19,8 @@ and FXN repeats (ExpansionHunter, step 09). With --cn-intervals (step 35's _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 _cnvs.vcf.gz) a CNVpytor call that disagrees is printed beside it. It writes, per couple A and @@ -399,7 +400,8 @@ def smn_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)")) diff --git a/docs/39-couple.md b/docs/39-couple.md index d7952c3f..76b823ec 100644 --- a/docs/39-couple.md +++ b/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/_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/_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 passes both when a partner's folder holds them (`paralogs/_cn_intervals.tsv` of step 35, `cnvpytor/_cnvs.vcf.gz` of step 18), and the Nextflow couple screen passes them when `--tools` has `parascopy` and `cnvpytor`. +`couple_report.py screen` also takes step 35's interval table (`--cn-intervals ID=paralogs/_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. With `--cnvpytor ID=cnvpytor/_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 passes both when a partner's folder holds them (`paralogs/_cn_intervals.tsv` of step 35, `cnvpytor/_cnvs.vcf.gz` of step 18), and the Nextflow couple screen passes them when `--tools` has `parascopy` and `cnvpytor`. ### The carrier list diff --git a/docs/cn-loci.md b/docs/cn-loci.md index c54652e8..f13f8563 100644 --- a/docs/cn-loci.md +++ b/docs/cn-loci.md @@ -1,6 +1,6 @@ # Copy Number from Depth (the locus catalog) -Several steps count copies of a gene from short-read depth: HBA1 and HBA2 in step 35, the RHD deletion in step 41 and the LPA KIV-2 repeat in step 46. Step 35 also counts three gene-dosage loci that no pinned caller models (HP, LCE3B/LCE3C and CCL3L1) and the interval rows: the recurrent CNV regions ClinGen curates (22q11.2, 1q21.1, 15q11.2, 16p11.2, 17p12 and the rest), three carrier deletion alleles and the Y chromosome AZF regions. They all use one engine, `bin/copy_number.py`, and one table, `bin/cn_loci.tsv`, with one row per locus. A new locus is a new row, not a new script. +Several steps count copies of a gene from short-read depth: HBA1 and HBA2 in step 35, the RHD deletion in step 41 and the LPA KIV-2 repeat in step 46. Step 35 also counts three gene-dosage loci that no pinned caller models (HP, LCE3B/LCE3C and CCL3L1) and the interval rows: the recurrent CNV regions ClinGen curates (22q11.2, 1q21.1, 15q11.2, 16p11.2, 17p12 and the rest), five carrier deletion alleles and the Y chromosome AZF regions. They all use one engine, `bin/copy_number.py`, and one table, `bin/cn_loci.tsv`, with one row per locus. A new locus is a new row, not a new script. The table sits in `bin/` because `bin/` is the one repository folder that every container mounts: the bash steps mount it as `/pgp-bin`, and Nextflow puts it on every task's PATH. @@ -78,7 +78,7 @@ The report does not show `_cn_dosage.tsv` yet; the table and the step lo - **The ClinGen recurrent regions**: every region row of ClinGen Dosage Sensitivity's recurrent CNV BED, version 2.1, GRCh38 (CC0, pinned by `CLINGEN_RECURRENT_CNV_URL` and its sha256 in `versions.env`). The file's breakpoint and LCR rows are the segmental duplications that bound the regions. They are cut out of the targets, as are GIAB v3.6 segmental duplications, assembly gaps, satellites and the regions of every pinned Parascopy model (CHRNA7 in 15q13.3 stays with Parascopy). Targets are cut to pieces of at most 500 kb, and their depth is pooled. Controls are two 100 kb windows, one on each side of the cluster of ClinGen rows the region belongs to, clear of all of that and of the targets of the hand rows below (so no row's control sits on another row's target); a side with no clear window within 3 Mb (a chromosome end, as next to 22q11.2 and 15q11.2) is left out. The BED spells its chrX rows `chrx`; the script normalises them, and the engine refuses any contig name that is not a GRCh38 chromosome. The classification of a CNV found here follows the ACMG and ClinGen technical standards (Riggs 2020, PMID 31690835); a call from this step is a lead for a clinical test, not a classification. - **Two more ClinGen curations**: the 1p36 terminal region (ISCA-37434, haploinsufficiency score 3) and the Xq28 region that holds MECP2 (ISCA-46304, triplosensitivity score 3: the MECP2 duplication), from ClinGen's region curation list, read on 2026-11-04. -- **Carrier deletion alleles**, at their ClinVar GRCh38 spans: del(GJB6-D13S1830) (variation 5546; del Castillo 2002, PMID 11807148), CFTRdele2,3 (variation 66105) and a 621 bp HBB deletion that removes exon 3 (variation 3769040). The alpha-globin `--SEA` deletion is read by the HBA row: two alpha-globin genes, with the reading that names `--SEA`. The 619 bp HBB deletion and the Sicilian delta-beta deletion have no row: we could not confirm their GRCh38 spans in a primary source. +- **Carrier deletion alleles**, at their ClinVar GRCh38 spans: del(GJB6-D13S1830) (variation 5546; del Castillo 2002, PMID 11807148), CFTRdele2,3 (variation 66105) and a 621 bp HBB deletion that removes exon 3 (variation 3769040); and, at their HbVar GRCh38 spans, the 619 bp HBB deletion common in India and Pakistan (HbVar 979, `NG_000007.3:g.71609_72227del619`, chr11:5,225,389-5,226,007; Orkin 1979, PMID 287080) and the 13.4 kb Sicilian (delta-beta)0 deletion from HBD to 3' of HBB (HbVar 1035, `NG_000007.3:g.64336_77738del13403`, chr11:5,219,878-5,233,280; Henthorn 1990, PMID 2307466). We checked both spans on 2026-11-04 by placing 50 bp on each side of the NG_000007.3 breakpoints on NC_000011.10; they agree with HbVar's hg38 positions. The rows count copies over their span only: the Sicilian deletion covers both HBB rows, so a Sicilian carrier reads as a deletion on all three, and the 619 bp and 621 bp rows overlap over 487 bp, so either deletion can show on both. The alpha-globin `--SEA` deletion is read by the HBA row: two alpha-globin genes, with the reading that names `--SEA`. - **AZF** (Y chromosome microdeletions, Krausz 2024, PMID 37674303): AZFa as USP9Y and DDX3Y, AZFb as KDM5D and EIF1AY (single-copy genes the deletion removes), and AZFc as the four DAZ copies added up, against ZFY and NLGN4Y. A male has 1 copy of AZFa and AZFb and 4 DAZ copies; a complete AZFc deletion leaves 0, the partial gr/gr and b2/b3 deletions 2. Which copies are lost is not claimed. Gene spans from Ensembl, checked on 2026-11-04. A row on chrX or chrY is read against the sample's own count there: its controls sit on the same chromosome, so a male's one X or Y gives 1 copy, and `expected` says so. The sex is the declared one: step 35's second argument, and in Nextflow the samplesheet's sex column. Step 16's inferred sex is not used, so both give the same table. A chrY row on a female is `not applicable`; without a sex, chrX and chrY rows are `not assessed` (`sex_unknown`). Every row needs 8x of control depth (4x on chrY), flanks within 1.5x of each other, and an estimate within 0.3 of a whole number; a control with no read makes that row `not assessed` (`broken_control`) and the others still run. A failed bedcov pass writes every row as `not assessed` (`check_failed`). @@ -109,7 +109,7 @@ pypgx owns the genes it calls from the BAM, the `bam` rows of `assets/pgx_genes. ## How it is tested -- `tests/test_copy_number.py` section 9 runs the converter on a small ClinGen BED whose chrX rows are spelled `chrx` (the rows load on chrX; the same rows left as `chrx` are refused; an unknown contig stops the build). Made-up depths call a 22q11.2 A-D deletion at 1 copy and a 17p12 duplication at 3, and 2 for both unaltered. AZFc at half its DAZ depth on a male is a deletion; every AZF row on a female is not applicable. No interval row's control overlaps another row's target. The interval block of step 35, cut out and run on made-up bedcov tables, writes one line per row, and a failed bedcov writes every row as not assessed. `tests/test_couple_report.py` section 11 checks the three states and the CNVpytor cross-check. +- `tests/test_copy_number.py` section 9 runs the converter on a small ClinGen BED whose chrX rows are spelled `chrx` (the rows load on chrX; the same rows left as `chrx` are refused; an unknown contig stops the build). Made-up depths call a 22q11.2 A-D deletion at 1 copy and a 17p12 duplication at 3, and 2 for both unaltered. AZFc at half its DAZ depth on a male is a deletion; every AZF row on a female is not applicable. No interval row's control overlaps another row's target. The HBB 619 bp and Sicilian rows sit at HbVar's GRCh38 spans. The interval block of step 35, cut out and run on made-up bedcov tables, writes one line per row, and a failed bedcov writes every row as not assessed. `tests/test_couple_report.py` section 11 checks the three states and the CNVpytor cross-check. - The e2e case `cn-intervals-1-catalog` rebuilds the interval rows from the pinned downloads and fails on any difference with the committed catalog. `cn-intervals-2-hg002` runs the interval blocks of step 35 and of the PARASCOPY module on HG002's real reads at the AZF genes, their controls and 22q11.2: a male with one copy of AZFa and AZFb, AZFc_DAZ not assessed (its estimate, between 3.3 and 3.9, is no whole count of the four DAZ copies, so the row says `between_counts` rather than a call), no change at 22q11.2 A-B, no call on a row off chrY and chr22 (they have no read; the chr22 rows share A-B's control window, so they are left out of that check), and the AZF rows not applicable when the same reads are read as female. - `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. diff --git a/scripts/ci/build-cn-intervals.py b/scripts/ci/build-cn-intervals.py index eeb53499..c676deab 100755 --- a/scripts/ci/build-cn-intervals.py +++ b/scripts/ci/build-cn-intervals.py @@ -18,8 +18,14 @@ ClinGen_region_curation_list_GRCh38.tsv as read on 2026-11-04; - three carrier deletion alleles at their ClinVar GRCh38 spans (HAND_ROWS): del(GJB6-D13S1830) (VCV 5546), CFTRdele2,3 (VCV 66105) and a 621 bp HBB - deletion that removes exon 3 (VCV 3769040). The alpha-globin --SEA - deletion is read by the HBA row (step 35's HBA table), not here; + deletion that removes exon 3 (VCV 3769040); and two beta-globin + deletions at their HbVar GRCh38 spans, checked on 2026-11-04 by placing + 50 bp on each side of the NG_000007.3 breakpoints on NC_000011.10: the + 619 bp HBB deletion common in India and Pakistan (HbVar 979, + NG_000007.3:g.71609_72227del619, Orkin 1979) and the 13.4 kb Sicilian + (delta-beta)0 deletion from HBD to 3' of HBB (HbVar 1035, + NG_000007.3:g.64336_77738del13403, Henthorn 1990). The alpha-globin + --SEA deletion is read by the HBA row (step 35's HBA table), not here; - the Y chromosome AZF regions (HAND_ROWS), read only on a sample with a Y: AZFa as USP9Y and DDX3Y, AZFb as KDM5D and EIF1AY (single-copy genes each deletion removes), and AZFc as the four DAZ copies added up (a complete @@ -95,6 +101,10 @@ "."), ("HBB_exon3_del", "chr11", [("HBB deletion of exon 3 (621 bp), ClinVar 3769040", 5225254, 5225875)], None, "pooled", True, "."), + ("HBB_619bp_del", "chr11", [("HBB 619 bp deletion (South Asian), HbVar 979", 5225388, 5226007)], None, + "pooled", True, "287080"), + ("HBD_HBB_Sicilian_del", "chr11", [("Sicilian (delta-beta)0 deletion (13.4 kb, HBD to HBB 3'), HbVar 1035", + 5219877, 5233280)], None, "pooled", True, "2307466"), ("AZFa", "chrY", [("AZFa (USP9Y and DDX3Y)", 12537649, 12860852), ("DDX3Y", 12904107, 12920478)], Y_CONTROLS, "pooled", False, "37674303"), ("AZFb", "chrY", [("AZFb (KDM5D and EIF1AY)", 19702297, 19745148), ("EIF1AY", 20575693, 20593154)], diff --git a/tests/test_copy_number.py b/tests/test_copy_number.py index 31795645..0a583285 100644 --- a/tests/test_copy_number.py +++ b/tests/test_copy_number.py @@ -608,7 +608,8 @@ def dosage(cat, read): chrx\t6527561\t6537770\tXp22.31 distal breakpoint\t0\t.\t6527561\t6537770\t0,0,0 chrx\t6537771\t8156913\tXp22.31 recurrent region (includes STS)\t0\t.\t6537771\t8156913\t255,111,0 """ -HAND = ["1p36_terminal", "Xq28_MECP2", "GJB6_D13S1830", "CFTR_dele2_3", "HBB_exon3_del", "AZFa", "AZFb", "AZFc_DAZ"] +HAND = ["1p36_terminal", "Xq28_MECP2", "GJB6_D13S1830", "CFTR_dele2_3", "HBB_exon3_del", "HBB_619bp_del", + "HBD_HBB_Sicilian_del", "AZFa", "AZFb", "AZFc_DAZ"] STUB_IV = r""" run_in() { while [ $# -gt 0 ]; do case $1 in -v|--cpus|--memory) shift 2 ;; *) break ;; esac; done @@ -674,8 +675,13 @@ def intervals(cat, read): iv = [r for r in cat.values() if r["kind"] == "interval" and r["_passes"] == ["."]] names = [r["locus"] for r in iv] - check("interval rows: the ClinGen regions (over 70), then the hand rows", len(iv) > 78 and names[-8:] == HAND, - (len(iv), names[-8:])) + check("interval rows: the ClinGen regions (over 70), then the hand rows", len(iv) > 80 and names[-10:] == HAND, + (len(iv), names[-10:])) + # HbVar's GRCh38 spans (979, 1035), 0-based start: one target each, never cut + spans = {loc: [(s, e) for _, s, e, _ in cat[loc]["_targets"]] if loc in cat else None + for loc in ("HBB_619bp_del", "HBD_HBB_Sicilian_del")} + check("interval rows: the HBB 619 bp and Sicilian deletions at HbVar's GRCh38 spans", + spans == {"HBB_619bp_del": [(5225388, 5226007)], "HBD_HBB_Sicilian_del": [(5219877, 5233280)]}, spans) check("interval rows: each cross-checked by CNVpytor, step 18", all(r["cross_check"] == "cnvpytor:cnvs" for r in iv)) check("interval rows: the 15q13.3 rows keep CHRNA7 (a Parascopy model) out of their targets", diff --git a/tests/test_couple_report.py b/tests/test_couple_report.py index 2688e7f3..de28377a 100644 --- a/tests/test_couple_report.py +++ b/tests/test_couple_report.py @@ -42,7 +42,8 @@ 11. carrier deletion alleles and AZF (step 35's interval table): each allele of each partner is listed as present, absent or not assessable; the male partner's AZF deletions are listed by that name (absent on a - normal male), the female's are not applicable; + normal male), the female's are not applicable; the HBB 619 bp and the + Sicilian (delta-beta)0 deletions are two of the alleles; where CNVpytor (step 18) and the depth call disagree, both are printed; without --cn-intervals the section is not written; 12. a split half call (./1) is het with one copy, never hom. @@ -390,11 +391,15 @@ def ivrow(locus, chrom, s, e, status="called", copies="2", expected="2", change= + ivrow("CFTR_dele2_3", "chr7", 117498312, 117519392) + ivrow("HBB_exon3_del", "chr11", 5225254, 5225875, status="not assessed", copies=".", change=".", reason="between_counts") + + ivrow("HBB_619bp_del", "chr11", 5225388, 5226007) + + ivrow("HBD_HBB_Sicilian_del", "chr11", 5219877, 5233280) + ivrow("AZFa", "chrY", 12537649, 12920478, status="not applicable", copies=".", expected=".", change=".", reason="no_y")) b_iv = put(f"{t}/B_cn_intervals.tsv", ivh + ivrow("GJB6_D13S1830", "chr13", 20223037, 20531806) + ivrow("CFTR_dele2_3", "chr7", 117498312, 117519392) - + ivrow("HBB_exon3_del", "chr11", 5225254, 5225875) + + ivrow("HBB_exon3_del", "chr11", 5225254, 5225875, copies="1", change="deletion") + + ivrow("HBB_619bp_del", "chr11", 5225388, 5226007, copies="1", change="deletion") + + ivrow("HBD_HBB_Sicilian_del", "chr11", 5219877, 5233280, copies="1", change="deletion") + ivrow("AZFa", "chrY", 12537649, 12920478, copies="1", expected="1") + ivrow("AZFb", "chrY", 19702297, 20593154, copies="1", expected="1") + ivrow("AZFc_DAZ", "chrY", 23129354, 24907040, copies="2", expected="4", change="deletion")) @@ -407,14 +412,19 @@ def ivrow(locus, chrom, s, e, status="called", copies="2", expected="2", change= "--cnvpytor", f"A={a_cnv}", "--cnvpytor", f"B={b_cnv}", "--out-dir", f"{t}/cn") md = open(f"{t}/cn/couple_A_B.md").read() if os.path.exists(f"{t}/cn/couple_A_B.md") else "" lines = [x for x in md.splitlines() if x.startswith("- A,") or x.startswith("- B,")] - check("cn alleles: each of the three alleles of each partner has one of the three states", + check("cn alleles: each of the five alleles of each partner has one of the three states", rc == 0 and all(sum(f"- {p}, {n}: **{st}**" in md for st in ("present", "absent", "not assessable")) == 1 - for p in "AB" for n in ("del(GJB6-D13S1830)", "CFTRdele2,3", "HBB exon 3 deletion")), + for p in "AB" for n in ("del(GJB6-D13S1830)", "CFTRdele2,3", "HBB exon 3 deletion", + "HBB 619 bp deletion", "Sicilian (delta-beta)0 deletion")), (rc, log, lines)) check("cn alleles: A's GJB6 deletion present, A's HBB not assessable with its reason, B's CFTR absent", "- A, del(GJB6-D13S1830): **present** (read depth, step 35: 1 of 2 copies)" in md and "- A, HBB exon 3 deletion: **not assessable** (read depth, step 35: not assessed (between_counts))" in md and "- B, CFTRdele2,3: **absent**" in md, lines) + check("cn alleles: B's Sicilian and 619 bp deletions present, A's absent", + "- B, Sicilian (delta-beta)0 deletion: **present** (read depth, step 35: 1 of 2 copies)" in md + and "- B, HBB 619 bp deletion: **present** (read depth, step 35: 1 of 2 copies)" in md + and "- A, Sicilian (delta-beta)0 deletion: **absent** (read depth, step 35: 2 of 2 copies)" in md, lines) check("cn alleles: where CNVpytor and the depth call disagree, both are printed", "- A, del(GJB6-D13S1830): **present** (read depth, step 35: 1 of 2 copies); " "CNVpytor (step 18) calls no deletion here" in md From ba87cb7bbaf6efd77e5c68f88c9fa648a7384b52 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Wed, 4 Nov 2026 15:49:00 +0100 Subject: [PATCH 4/9] fix(cn): match the GJB6 line literally in the couple case and say the beta-globin rows overlap --- bin/couple_report.py | 5 ++++- docs/39-couple.md | 2 +- tests/e2e/followups-cn-1-couple-report.sh | 2 +- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/bin/couple_report.py b/bin/couple_report.py index 82dc85e2..e0e9b53c 100755 --- a/bin/couple_report.py +++ b/bin/couple_report.py @@ -615,7 +615,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 |", diff --git a/docs/39-couple.md b/docs/39-couple.md index 76b823ec..936849c0 100644 --- a/docs/39-couple.md +++ b/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/_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. With `--cnvpytor ID=cnvpytor/_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 passes both when a partner's folder holds them (`paralogs/_cn_intervals.tsv` of step 35, `cnvpytor/_cnvs.vcf.gz` of step 18), and the Nextflow couple screen passes them when `--tools` has `parascopy` and `cnvpytor`. +`couple_report.py screen` also takes step 35's interval table (`--cn-intervals ID=paralogs/_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/_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 passes both when a partner's folder holds them (`paralogs/_cn_intervals.tsv` of step 35, `cnvpytor/_cnvs.vcf.gz` of step 18), and the Nextflow couple screen passes them when `--tools` has `parascopy` and `cnvpytor`. ### The carrier list diff --git a/tests/e2e/followups-cn-1-couple-report.sh b/tests/e2e/followups-cn-1-couple-report.sh index 59570fd8..4b806910 100755 --- a/tests/e2e/followups-cn-1-couple-report.sh +++ b/tests/e2e/followups-cn-1-couple-report.sh @@ -65,7 +65,7 @@ printf '##fileformat=VCFv4.2\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\nchr | gzip -n > "${C}/${B}/cnvpytor/${B}_cnvs.vcf.gz" printf 'sample\tsex\tfather\tmother\tpartner\n%s\tmale\t.\t.\t%s\n%s\tfemale\t.\t.\t%s\n' "$A" "$B" "$B" "$A" > "${C}/pedigree.tsv" -GJB6_A="- ${A}, del(GJB6-D13S1830): \*\*present\*\* \(read depth, step 35: 1 of 2 copies\); CNVpytor \(step 18\) calls no deletion here" +GJB6_A="- ${A}, del\(GJB6-D13S1830\): \*\*present\*\* \(read depth, step 35: 1 of 2 copies\); CNVpytor \(step 18\) calls no deletion here" AZFC_A="- ${A}, AZFc deletion \(DAZ copies; gr/gr and b2/b3 leave two of four\): \*\*present\*\* \(read depth, step 35: 2 of 4 copies\)" CFTR_B="- ${B}, CFTRdele2,3: \*\*absent\*\* \(read depth, step 35: 2 of 2 copies\); CNVpytor \(step 18\) calls a deletion here" AZF_B="- ${B}, AZF regions: not applicable \(no Y chromosome\)" From 0c490e0e98aea5515707c3b55cae7d7a4378e7f3 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Wed, 4 Nov 2026 16:03:56 +0100 Subject: [PATCH 5/9] fix(cn): print no missing CNVpytor call for an allele under two of its 1 kb bins --- bin/couple_report.py | 5 ++++- docs/39-couple.md | 2 +- tests/test_couple_report.py | 5 +++++ 3 files changed, 10 insertions(+), 2 deletions(-) diff --git a/bin/couple_report.py b/bin/couple_report.py index e0e9b53c..22717f4b 100755 --- a/bin/couple_report.py +++ b/bin/couple_report.py @@ -407,6 +407,9 @@ def smn_finding(path, genes): ("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): @@ -455,7 +458,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" diff --git a/docs/39-couple.md b/docs/39-couple.md index 936849c0..466ebbb0 100644 --- a/docs/39-couple.md +++ b/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/_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/_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 passes both when a partner's folder holds them (`paralogs/_cn_intervals.tsv` of step 35, `cnvpytor/_cnvs.vcf.gz` of step 18), and the Nextflow couple screen passes them when `--tools` has `parascopy` and `cnvpytor`. +`couple_report.py screen` also takes step 35's interval table (`--cn-intervals ID=paralogs/_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/_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/_cn_intervals.tsv` of step 35, `cnvpytor/_cnvs.vcf.gz` of step 18), and the Nextflow couple screen passes them when `--tools` has `parascopy` and `cnvpytor`. ### The carrier list diff --git a/tests/test_couple_report.py b/tests/test_couple_report.py index de28377a..cb380e09 100644 --- a/tests/test_couple_report.py +++ b/tests/test_couple_report.py @@ -432,6 +432,11 @@ def ivrow(locus, chrom, s, e, status="called", copies="2", expected="2", change= "deletion here" in md and "- B, GJB6_D13S1830" not in md and "- B, del(GJB6-D13S1830): **absent** (read depth, step 35: 2 of 2 " "copies)\n" in md + "\n", lines) + check("cn alleles: no CNVpytor note under two of its 1 kb bins (B's 619 bp, 621 bp); the 13.4 kb Sicilian gets it", + "- B, HBB 619 bp deletion: **present** (read depth, step 35: 1 of 2 copies)\n" in md + and "- B, HBB exon 3 deletion: **present** (read depth, step 35: 1 of 2 copies)\n" in md + and "- B, Sicilian (delta-beta)0 deletion: **present** (read depth, step 35: 1 of 2 copies); " + "CNVpytor (step 18) calls no deletion here" in md, lines) check("cn alleles: the male partner's AZF regions, AZFc a deletion; the female's not applicable", "- B, AZFc deletion (DAZ copies; gr/gr and b2/b3 leave two of four): **present** (read depth, step 35: " "2 of 4 copies)" in md and "- B, AZFa deletion: **absent** (read depth, step 35: 1 of 1 copies)" in md From 75f5e762e1b23c287746a20f5c88a9ff3e8f2127 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Wed, 4 Nov 2026 18:53:39 +0100 Subject: [PATCH 6/9] test(cn): count the 90 interval rows in the catalog and HG002 cases --- tests/e2e/cn-intervals-1-catalog.sh | 2 +- tests/e2e/cn-intervals-2-hg002.sh | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/e2e/cn-intervals-1-catalog.sh b/tests/e2e/cn-intervals-1-catalog.sh index 640cb7c3..21b741f7 100755 --- a/tests/e2e/cn-intervals-1-catalog.sh +++ b/tests/e2e/cn-intervals-1-catalog.sh @@ -11,7 +11,7 @@ mkdir -p "$CACHE" echo "+ scripts/ci/build-cn-intervals.py --check" python3 "${REPO}/scripts/ci/build-cn-intervals.py" --cache "$CACHE" --check 2>&1 | tee "$STEP_LOG" check_eq "build-cn-intervals.py --check exits 0 on the committed catalog" "${PIPESTATUS[0]}" 0 -check "the log says all 88 interval rows were checked" grep -q "the 88 interval rows are what the pinned downloads give" "$STEP_LOG" +check "the log says all 90 interval rows were checked" grep -q "the 90 interval rows are what the pinned downloads give" "$STEP_LOG" # Control: one interval target moved by 1 kb EDITED="${CASE_TMP}/edited.tsv" diff --git a/tests/e2e/cn-intervals-2-hg002.sh b/tests/e2e/cn-intervals-2-hg002.sh index 457cd324..371861a0 100755 --- a/tests/e2e/cn-intervals-2-hg002.sh +++ b/tests/e2e/cn-intervals-2-hg002.sh @@ -59,7 +59,7 @@ step_block male T="${W}/paralogs/${P}_cn_intervals.tsv" cp "$T" "${CASE_TMP}/male.tsv" awk -F'\t' '$1 ~ /^(AZF|22q11.2_DGS_proximal_A-B)/' "$T" -check_ge "the table lists every interval row" "$(($(wc -l < "$T") - 1))" 88 +check_ge "the table lists every interval row" "$(($(wc -l < "$T") - 1))" 90 for r in AZFa AZFb; do check_eq "male: ${r} called" "$(col "$T" "$r" status)" called check_eq "male: ${r} 1 copy of 1, no change" "$(col "$T" "$r" copies)/$(col "$T" "$r" expected)/$(col "$T" "$r" change)" \ From 6b983577e683d4906a42d2900d8e96417dc33087 Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Wed, 4 Nov 2026 19:31:45 +0100 Subject: [PATCH 7/9] fix(cn): keep the interval rows in the summary when step 35 wrote no other table --- bin/collect_summary.py | 5 +++-- tests/test_copy_number.py | 5 +++++ 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/bin/collect_summary.py b/bin/collect_summary.py index f18c3d02..4565ab2b 100755 --- a/bin/collect_summary.py +++ b/bin/collect_summary.py @@ -842,7 +842,9 @@ 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") + 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)) if not p: return None, {} if "smn" in tables: @@ -865,7 +867,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") 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), diff --git a/tests/test_copy_number.py b/tests/test_copy_number.py index 0a583285..8a623e6d 100644 --- a/tests/test_copy_number.py +++ b/tests/test_copy_number.py @@ -801,6 +801,11 @@ def intervals(cat, read): ci = data.get("cn_intervals") or {} check("collect_summary: step 35's section lists the interval rows with a change (AZFc) and counts the rest", [c["locus"] for c in ci.get("changes", [])] == ["AZFc_DAZ"] and ci.get("rows") == len(iv), ci) + os.remove(os.path.join(sd, "paralogs", "S_smn_copy_number.tsv")) + _, data = cs.sec_parascopy(sd, "S") + ci = (data or {}).get("cn_intervals") or {} + check("collect_summary: the interval rows stay when step 35 wrote no other table", + [c["locus"] for c in ci.get("changes", [])] == ["AZFc_DAZ"] and ci.get("rows") == len(iv), ci) check("PARASCOPY module: the samplesheet sex goes to the interval rows", "def sex = meta.sex in ['male', 'female'] ? meta.sex : 'unknown'" in module and "--sex ${sex}" in module) From bebf10f10fa06c050ba76c1cb6c7efb2e2d911ed Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Thu, 5 Nov 2026 19:21:01 +0100 Subject: [PATCH 8/9] docs(cn): name the interval table among the report's paralog inputs and list the couple-report e2e case --- docs/cn-loci.md | 1 + modules/local/html_report/main.nf | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/cn-loci.md b/docs/cn-loci.md index f13f8563..b4192083 100644 --- a/docs/cn-loci.md +++ b/docs/cn-loci.md @@ -111,6 +111,7 @@ pypgx owns the genes it calls from the BAM, the `bam` rows of `assets/pgx_genes. - `tests/test_copy_number.py` section 9 runs the converter on a small ClinGen BED whose chrX rows are spelled `chrx` (the rows load on chrX; the same rows left as `chrx` are refused; an unknown contig stops the build). Made-up depths call a 22q11.2 A-D deletion at 1 copy and a 17p12 duplication at 3, and 2 for both unaltered. AZFc at half its DAZ depth on a male is a deletion; every AZF row on a female is not applicable. No interval row's control overlaps another row's target. The HBB 619 bp and Sicilian rows sit at HbVar's GRCh38 spans. The interval block of step 35, cut out and run on made-up bedcov tables, writes one line per row, and a failed bedcov writes every row as not assessed. `tests/test_couple_report.py` section 11 checks the three states and the CNVpytor cross-check. - The e2e case `cn-intervals-1-catalog` rebuilds the interval rows from the pinned downloads and fails on any difference with the committed catalog. `cn-intervals-2-hg002` runs the interval blocks of step 35 and of the PARASCOPY module on HG002's real reads at the AZF genes, their controls and 22q11.2: a male with one copy of AZFa and AZFb, AZFc_DAZ not assessed (its estimate, between 3.3 and 3.9, is no whole count of the four DAZ copies, so the row says `between_counts` rather than a call), no change at 22q11.2 A-B, no call on a row off chrY and chr22 (they have no read; the chr22 rows share A-B's control window, so they are left out of that check), and the AZF rows not applicable when the same reads are read as female. +- The e2e case `followups-cn-1-couple-report` runs step 39, COUPLE_CARRIER and HTML_REPORT on a made-up couple, plus a stub run of `main.nf`, and checks that the interval table and the CNVpytor calls reach the couple report and `summary.json`. Each wiring line has a red control: the case runs again without it and checks that the rows are gone. - `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`. diff --git a/modules/local/html_report/main.nf b/modules/local/html_report/main.nf index c13c7e02..ce2aa6ef 100644 --- a/modules/local/html_report/main.nf +++ b/modules/local/html_report/main.nf @@ -22,7 +22,8 @@ and go to assets/, which --assets names, as do assets/mito_confirmed.tsv (the confirmed mtDNA variants step 20's section names) and the paralog catalogs assets/paralog_loci.tsv and assets/paraphase_genes.tsv. Every - _
_copy_number.tsv of step 35 goes to paralogs/. Every gene and rsID the + _
_copy_number.tsv of step 35 goes to paralogs/, and so + does step 35's _cn_intervals.tsv. Every gene and rsID the results mention also goes to ${id}_search_index.json beside it (bin/results_index.py), for the live site's search, and the report's charts to ${id}_charts.json, for the live page (docs/24-html-report.md). From 2662e8f67986ec670983c28ea44a00b23456c47c Mon Sep 17 00:00:00 2001 From: GeiserX <9169332+GeiserX@users.noreply.github.com> Date: Thu, 5 Nov 2026 20:09:14 +0100 Subject: [PATCH 9/9] ci(e2e): list the files the couple-report case runs --- scripts/ci/e2e-case-paths.tsv | 1 + 1 file changed, 1 insertion(+) diff --git a/scripts/ci/e2e-case-paths.tsv b/scripts/ci/e2e-case-paths.tsv index 98e368d2..4c60e629 100644 --- a/scripts/ci/e2e-case-paths.tsv +++ b/scripts/ci/e2e-case-paths.tsv @@ -87,6 +87,7 @@ families-file-1 assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_g family-joint-1 assets/construct_variants.tsv assets/deep_intronic_sites.tsv assets/families.example.tsv assets/par_grch38.bed assets/pedigree.example.tsv assets/stub/pedigree.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 bin/collect_summary.py bin/family_report.py bin/family_segments.py bin/ibd_report.py bin/mosaic_sex_loss.py modules/local/align_minimap2/* modules/local/clinvar_screen/* modules/local/deepvariant/* modules/local/fastp/* modules/local/roh/* modules/local/somalier/* modules/local/vcf_precheck/* scripts/01b-fastp-qc.sh scripts/06-clinvar-screen.sh scripts/11-roh-analysis.sh scripts/33-sample-qc.sh scripts/38-family.sh scripts/ci/check-images.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh fixitfelix-chr21-1 assets/callability_genes.tsv scripts/16b-mosdepth.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh followups-clinvar-traits-1-rules-html assets/genotype_rules.tsv assets/traits.tsv bin/render_report.py bin/traits_report.py scripts/42-traits.sh scripts/ci/check-rule-alleles.sh +followups-cn-1-couple-report assets/carrier_genes.tsv assets/par_grch38.bed assets/stub/pedigree.tsv bin/collect_summary.py bin/couple_report.py bin/family_report.py bin/render_report.py bin/results_index.py modules/local/couple/* modules/local/html_report/* scripts/39-couple.sh scripts/ci/fake-docker/lib.sh genome-lookup-1 scripts/ci/fake-docker/lib.sh scripts/owngenome-lookup.sh genotype-rule-rows-1-str-extra assets/str_extra_loci.json assets/stub/db/catalog.json scripts/09-expansion-hunter.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh genotype-rules-1-phase assets/genotype_rules.tsv assets/traits.tsv bin/traits_report.py modules/local/traits/* modules/local/vcf_precheck/* scripts/42-traits.sh scripts/ci/check-rule-alleles.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh