fix(nextflow): give the Nextflow report and couple screen the CYP21A2, dosage, gVCF and release data the bash path already uses

README.md

@@ -30,7 +30,7 @@ This pipeline takes raw sequencing data (FASTQ/BAM/VCF) from any vendor and runs
 - Calls SNPs and indels with DeepVariant, and structural and copy number variants with Manta, Delly and CNVpytor merged into a consensus.
 - Screens ClinVar and runs CPSR cancer predisposition panels, VEP annotation with CADD, SpliceAI, REVEL and AlphaMissense, and slivar prioritization.
 - Pharmacogenomics with PharmCAT, pypgx (29 genes, 13 of them with copy number from the BAM, so gene deletions, duplications and hybrid genes) and CPIC drug recommendations. PharmCAT gets T1K's HLA types, and a CYP2D6 call only when pypgx and Cyrius (opt-in, non-commercial licence) agree.
-- Repeat expansions, HLA typing (KIR opt-in), SMN1/SMN2 copy number (opt-in), and from PacBio HiFi reads HBA1/HBA2, CYP21A2, GBA1, PMS2, STRC, NCF1 and OTOA with Paraphase (opt-in), telomere length, mitochondrial haplogroup, heteroplasmy and an mtDNA contamination check, Y-chromosome haplogroup (opt-in), mosaic loss of chrY or chrX with its uncertainty, ROH, ancestry, admixture fractions with a chromosome painting (opt-in) and polygenic risk scores, and the blood group antigens (ABO, Rh, Kell, Duffy, Kidd and the other systems, opt-in).
+- Repeat expansions, HLA typing (KIR opt-in), SMN1/SMN2 copy number and, from short reads, CYP21A2 with its 30 kb deletion and pseudogene conversions (opt-in, step 35), and from PacBio HiFi reads HBA1/HBA2, CYP21A2, GBA1, PMS2, STRC, NCF1 and OTOA with Paraphase (opt-in), telomere length, mitochondrial haplogroup, heteroplasmy and an mtDNA contamination check, Y-chromosome haplogroup (opt-in), mosaic loss of chrY or chrX with its uncertainty, ROH, ancestry, admixture fractions with a chromosome painting (opt-in) and polygenic risk scores, and the blood group antigens (ABO, Rh, Kell, Duffy, Kidd and the other systems, opt-in).
 - A raw-data file in the 23andMe layout from your genome (opt-in step 40), for genealogy and health services that read chip files. The pipeline writes it; uploading it is up to you.
 - A traits panel: lactase persistence, caffeine, alcohol flush, bitter taste, earwax, APOE and other replicated variants, read from the gVCF, each with the primary paper it comes from.
 - Every tool runs in a Docker container with a hard memory limit and no network; CPU is a hard cap in the single-step scripts and a Docker share in the pipeline ([a shared host](https://geiserx.github.io/OwnGenome/hardware-requirements/#a-shared-host)). Each image is pinned by tag or digest in `versions.env`. The opt-in Cyrius is installed once by `setup.sh --cyrius` from PyPI, every file checked against its hash, and then runs offline too. No script uploads your data; [a few steps download public files](https://geiserx.github.io/OwnGenome/why-local/#network-calls-during-a-run) during a run.

assets/steps.tsv

@@ -37,7 +37,7 @@ step	label	tool	recorded	processes	sections	upstream	hours	hours_from	hours_low
 15	duphold	duphold	yes	DUPHOLD,DUPHOLD_FILTER		04	0.07	measured	0.05	
 16	indexcov (sex check)		yes	INDEXCOV	sex_check	02	0.01	measured	0.0002	
 16b	mosdepth	mosdepth	yes	MOSDEPTH	coverage,callability	02	0.05	measured	0.0333	
-17	CPSR	cpsr	yes	CPSR	cpsr	03	0.6	measured	0.2833	PCGR_DATA_BUNDLE,cpsr_clinvar_release
+17	CPSR	cpsr	yes	CPSR,CPSR_CDNA_CHECK	cpsr	03	0.6	measured	0.2833	PCGR_DATA_BUNDLE,cpsr_clinvar_release
 18	CNVpytor	cnvpytor	yes	CNVPYTOR,CNVPYTOR_VCF	cnvpytor	02	0.33	measured	0.2333	
 19	Delly	delly	yes	DELLY,DELLY_BCF2VCF	delly	02	3.32	measured	0.35	
 20	Mito variants (Mutect2)	mito_variants	yes	MITO_VARIANTS,MITO_VARIANTS_ANNOTATE,MITO_DELETIONS	mito	02	0.03	measured	0.0166	clinvar_file_date@self

bin/write_manifest.sh

@@ -5,7 +5,8 @@
 #
 # Writes ${GENOME_DIR}/<sample>/run_manifest.tsv, three tab-separated columns
 # (section, key, value):
-#   run      written_utc, written_by, git_commit, declared_sex
+#   run      written_utc, written_by, git_commit, pipeline_version (the OwnGenome
+#            release: the version in nextflow.config's manifest block), declared_sex
 #   version  every NAME=VALUE line of versions.env
 #   image    each *_IMAGE reference and the digest Docker resolved it to
 #            ("not pulled" when the image is not on this machine)
@@ -53,6 +54,19 @@ if git -C "$PGP_ROOT" rev-parse --verify -q HEAD >/dev/null 2>&1; then
   fi
 fi
 row run git_commit "$commit"
+# The OwnGenome release: the version = '...' line inside the manifest { } block
+# of nextflow.config (the release tag is v<version>). No row when there is none.
+release=$(awk '
+  /^[[:space:]]*manifest[[:space:]]*\{/ {inside = 1; next}
+  inside && /^[[:space:]]*\}/ {inside = 0}
+  inside && /^[[:space:]]*version[[:space:]]*=/ {
+    s = $0; sub(/^[^=]*=[[:space:]]*/, "", s); sub(/[[:space:]]*(\/\/.*)?$/, "", s)
+    gsub(/^["\047]|["\047]$/, "", s); print s; exit
+  }' "${PGP_ROOT}/nextflow.config" 2>/dev/null || true)
+if [ ! -r "${PGP_ROOT}/nextflow.config" ]; then
+  echo "WARNING: cannot read ${PGP_ROOT}/nextflow.config: run_manifest.tsv gets no pipeline_version row" >&2
+fi
+[ -z "$release" ] || row run pipeline_version "$release"
 [ -z "$DECLARED_SEX" ] || row run declared_sex "$DECLARED_SEX"
 
 # versions.env as written (comments and blank lines dropped)

conf/containers.config

@@ -33,6 +33,7 @@ process {
     withName: 'COUPLE_CARRIER' { container = 'python:3.15.2@sha256:d7a732b20de2f5d3cf977cbcf69aacf30dc928a37908ff9fdb7bb3503024c91a' }
     withName: 'CPIC_LOOKUP' { container = 'python:3.15.2@sha256:d7a732b20de2f5d3cf977cbcf69aacf30dc928a37908ff9fdb7bb3503024c91a' }
     withName: 'CPSR' { container = 'ghcr.io/sigven/pcgr:2.3.2' }
+    withName: 'CPSR_CDNA_CHECK' { container = 'quay.io/biocontainers/pypgx:0.26.0--pyh7e72e81_0' }
     withName: 'CRAM_ARCHIVE' { container = 'staphb/samtools:1.24' }
     withName: 'CRAM_TO_BAM' { container = 'staphb/samtools:1.24' }
     withName: 'CYP2D6_DEPTH' { container = 'quay.io/biocontainers/mosdepth:0.3.14--h87be163_2' }

docs/17-cpsr.md

@@ -121,7 +121,7 @@ Run by hand like this, `--sample_id` must be 3 to 40 characters long. The script
 - `${SAMPLE}.cpsr.grch38.classification.tsv.gz` — Tab-separated variant classifications (gzipped; read it with `zcat`)
 - Every variant in the panel genes gets one of the five ACMG/AMP classes: Pathogenic, Likely pathogenic, VUS, Likely benign, Benign
 - In the TSV, `CLASSIFICATION` is the final class and `ASSERTION_AUTHORITY` says whose it is: ClinVar's, unless ClinVar has no record or a conflicted one, then CPSR's own (`CPSR_CLASSIFICATION`, which every variant has). PCGR 2.2.5 named the final column `FINAL_CLASSIFICATION`; the reports (`bin/collect_summary.py`) read either. The file names are the same in both versions.
-- `${SAMPLE}_cpsr_cdna_check.tsv` — step 6's cDNA read check on CPSR's P/LP rows (final class or CPSR's own), in the columns of step 6's `clinvar/${SAMPLE}_clinvar_cdna_check.tsv`. The rows are read from `VAR_ID` (or `GENOMIC_CHANGE`) and keep CPSR's contig names (`13`, not `chr13`); the BAM is read under either name. It needs the sample's BAM, `aligned/${SAMPLE}_sorted.bam` with its `.bai` (`ALIGN_DIR` picks another folder), and the GENCODE exons `setup.sh` installs; without either, every row is "not checked" with the reason. In Nextflow the `CLINVAR_CDNA_CHECK` process runs it after `CPSR` (source `cpsr`), with the samplesheet's BAM and `--gencode_exons`, and publishes the table to `cpsr/`; the HTML report puts it beside each P/LP row there too.
+- `${SAMPLE}_cpsr_cdna_check.tsv` — step 6's cDNA read check on CPSR's P/LP rows (final class or CPSR's own), in the columns of step 6's `clinvar/${SAMPLE}_clinvar_cdna_check.tsv`. The rows are read from `VAR_ID` (or `GENOMIC_CHANGE`) and keep CPSR's contig names (`13`, not `chr13`); the BAM is read under either name. It needs the sample's BAM, `aligned/${SAMPLE}_sorted.bam` with its `.bai` (`ALIGN_DIR` picks another folder), and the GENCODE exons `setup.sh` installs; without either, every row is "not checked" with the reason. In Nextflow step 6's `CLINVAR_CDNA_CHECK` process runs it after `CPSR` (source `cpsr`) under the name `CPSR_CDNA_CHECK`, so the run's process table and the live view count it under step 17, with the samplesheet's BAM and `--gencode_exons`, and publishes the table to `cpsr/`; the HTML report puts it beside each P/LP row there too. Its image comes from a `withName: 'CPSR_CDNA_CHECK'` selector (a setting for step 6's `CLINVAR_CDNA_CHECK` does not reach it).
 
 ### What the reports show
 

docs/24-html-report.md

@@ -24,7 +24,7 @@ All output directories from previous pipeline steps. The script detects which st
 
 Two more files, when present:
 
-- `${SAMPLE}/run_manifest.tsv` (`bin/write_manifest.sh`): the pipeline commit, `versions.env`, the digest of every image, the ClinVar file date, the VEP cache, the PCGR and pypgx bundles, whether pgsc_calc's ancestry panel is installed, the IPD-IMGT/HLA release of the `hla.dat` the run used (`run-all.sh` passes it as `HLA_DAT`; without it the one `setup.sh` installs, then the release in the T1K index folder name, and `unknown` when HLA calls exist but no database is found) and the header of each PGS scoring file. `run-all.sh` writes it when a run starts and again just before this step, so an image pulled during the run gets its digest; this step writes one when the sample has none.
+- `${SAMPLE}/run_manifest.tsv` (`bin/write_manifest.sh`): the pipeline commit, the OwnGenome release (the `manifest` version in `nextflow.config`, row `run pipeline_version`), `versions.env`, the digest of every image, the ClinVar file date, the VEP cache, the PCGR and pypgx bundles, whether pgsc_calc's ancestry panel is installed, the IPD-IMGT/HLA release of the `hla.dat` the run used (`run-all.sh` passes it as `HLA_DAT`; without it the one `setup.sh` installs, then the release in the T1K index folder name, and `unknown` when HLA calls exist but no database is found) and the header of each PGS scoring file. `run-all.sh` writes it when a run starts and again just before this step, so an image pulled during the run gets its digest; this step writes one when the sample has none.
 - `${SAMPLE}/logs/run_status.tsv`: when the latest `run-all.sh` run started, then one line per step. A step it did not start is `skipped (<reason>)`, for example `skipped (no BAM)` or `skipped (not in TOOLS)`. Each step it starts gets a `planned` row just before the pipeline starts, which stands until the step's result row replaces it; a `meta finished_utc` row closes the file when `run-all.sh` ends, and a step still `planned` then becomes `not finished` ([Where a run stands](run-state.md)). Once the pipeline has ended, each step it started is one of:
     - `ok`: the pipeline finished. After a pipeline that another step stopped, a step is still `ok` when each of its processes has a task that completed in this run or came from Nextflow's cache (`-resume`), and none was left unfinished.
     - `failed`: one of the step's tasks failed. A report-only tool's failure is ignored by the pipeline, which goes on without it. Ancestry (step 26) is computed inside the PRS run, so it is `failed` when PRS is.

docs/35-paralogs.md

@@ -244,6 +244,7 @@ The ratio of a site is its estimate over the median of the sites, so one copy co
 - The e2e case `cyp21a2-short-reads-1-hg002` runs the CYP21A2 block of step 35 on HG002's Illumina reads over the row, unmodified and with a planted In2G conversion and a planted chimeric deletion made from the same reads, and the block of the PARASCOPY module on the unmodified BAM, which must write the same table.
 - The e2e case `hba-short-read-accuracy-1-nf-report` runs the Nextflow `HTML_REPORT` with and without an HBA table and checks the HBA line and the not-assessed list, then stub-runs every tool and checks that `main.nf` passes the table to each `HTML_REPORT` task that got the SMN table, against a copy of the pipeline with that join removed.
 - `tests/test_collect_summary.py` (`cn_report_section`) reads made-up CYP21A2, dosage and interval tables into the summary and both reports: a CYP21A2 deletion, HP 3, LCE3BC 1, CCL3L1 not assessed, a 22q11.2 deletion with its counts; without the tables each line says the table is not in the output; a header without the copy-number columns is unreadable for the dosage, interval, catalog and SMN tables, also with no row.
+- The e2e case `followups-nf-1-nf-wiring` stub-runs the pipeline on the stub pedigree and checks that `main.nf` passes the CYP21A2 and dosage tables to each `HTML_REPORT` task that got the interval table, and each partner's CYP21A2 table to `COUPLE_CARRIER`, against a copy of the pipeline without those lines.
 - The e2e case `cn-report-1-tables` runs steps 24, `generate-report.sh` and 39 in their images on made-up tables: both reports show the CYP21A2 deletion, the dosage loci and the 22q11.2 deletion; step 39 passes each partner's CYP21A2 table, and CFTRdele2,3 present in both partners is named under the shared-risk heading. Red controls: the CYP21A2 table removed, for the report and for one partner.
 
 ## In the Reports
@@ -258,7 +259,7 @@ The same section shows the tables no Parascopy model makes, each read on its own
 
 A table that is not in the run's output reads `not assessed (step 35's <id>_... is not in this run's output)`, so a missing table never reads as a normal result. An SMN, interval or dosage table with a header but without its columns makes the section unreadable, and a CYP21A2 table without a status row reads as not assessed; a catalog table without the copy-number columns makes that locus unreadable, also with no row under the header.
 
-The couple screen reads the CYP21A2 table with `couple_report.py screen --cyp21a2 ID=<id>_cyp21a2_psv.tsv`, where a deletion or a possible conversion is a CYP21A2 carrier finding. Step 39 passes it when a partner's folder holds it. The Nextflow `PARASCOPY` task writes the CYP21A2 and dosage tables, but `main.nf` does not pass them to `HTML_REPORT` or the couple screen yet, so a Nextflow run's reports say they are not in the output.
+The couple screen reads the CYP21A2 table with `couple_report.py screen --cyp21a2 ID=<id>_cyp21a2_psv.tsv`, where a deletion or a possible conversion is a CYP21A2 carrier finding. Step 39 passes it when a partner's folder holds it. In Nextflow, `main.nf` passes the `PARASCOPY` task's CYP21A2 and dosage tables to each sample's `HTML_REPORT`, and the CYP21A2 table of each partner to `COUPLE_CARRIER`, so both reports and the couple screen read them as the bash steps do.
 
 When neither step 35 nor step 35b left a result, both reports print `Paralog copy number: not checked (neither step 35 nor step 35b ran)`, naming a step that failed instead, so a missing check never reads as a normal result. Unless step 35's depth count made a call, or step 35b called both HBA1 and HBA2 (a run limited to other regions with `--paraphase_regions` has neither, and a `no call` row does not count), the "Not Assessed" list also names HBA1/HBA2 copy number.
 

docs/39-couple.md

@@ -25,7 +25,7 @@ Before this step a couple had to compare two ClinVar tables by hand ([Multi-Samp
 | ClinVar pathogenic or likely pathogenic hit | step 06, `clinvar/<sample>_clinvar_hits.vcf` (required) | the genotype carries the alternate allele and ClinVar's `GENEINFO` names a gene of the carrier list |
 | Deletion | step 19 (Delly), `delly/<sample>_sv.vcf.gz` | `SVTYPE=DEL`, `FILTER` PASS and a per-sample `FT` that is not a failed filter, the genotype carries it, it is 50 bp or longer, it removes at least one exon of a listed gene (GENCODE exons), when it is longer than 1,000 bp read depth shows the loss, and it does not match a common deletion (see below) |
 | SMN1 copy number | step 35 (Parascopy), `paralogs/<sample>_smn_copy_number.tsv` | SMN1 exon 7 at 1 copy (a carrier) or 0, with filter PASS and quality 20 or more, by step 35's rule |
-| CYP21A2 deletion or conversion | step 35 (depth and six PSV sites), `paralogs/<sample>_cyp21a2_psv.tsv` | the table's `finding` is `deletion` or `conversion` (on both copies when a site has no CYP21A2 copy left); unconfirmed, short reads. Step 39 passes it; the Nextflow couple screen does not yet |
+| CYP21A2 deletion or conversion | step 35 (depth and six PSV sites), `paralogs/<sample>_cyp21a2_psv.tsv` | the table's `finding` is `deletion` or `conversion` (on both copies when a site has no CYP21A2 copy left); unconfirmed, short reads. Step 39 and the Nextflow couple screen (`COUPLE_CARRIER`) both pass it |
 | FMR1 and FXN repeats | step 09 (ExpansionHunter), `expansion_hunter/<sample>_eh.vcf` | an FMR1 allele of 55 CGG or more (a premutation; more than 200 is a full mutation), an FXN allele of 66 GAA or more |
 
 Short-read Delly can mark PASS a long deletion, even one of many megabases, that read depth contradicts, and such a call removes an exon of hundreds of listed genes at once. So a deletion longer than 1,000 bp counts only when Delly's own read-depth copy number (`FORMAT/RDCN`, the depth inside the deletion against its flanks) is below a normal region's copies: 1 on a male's chrX or chrY outside the pseudoautosomal regions, 2 everywhere else. When the pedigree gives no sex, every region is judged against 2. Delly infers the sex itself and scales RDCN to it, so a long artefact on an unsexed male's chrX or chrY reads 1 and still counts: give the sex in the pedigree and in the samplesheet. A heterozygous loss reads 1, a homozygous or hemizygous loss 0. A long deletion with RDCN at the normal count or more, or with no estimate (`.` or a negative value, no reads in the flanks), is set aside, and the table of what was assessed gives the count. 1,000 bp is Delly's own line between small and large SVs in its read counts (`indelsize` in Delly's source); a deletion that size or smaller rests on its split reads and read pairs alone. A deletion VCF without the RDCN field (not written by Delly 2.x) is read as before, without the depth check.

docs/42-traits.md

@@ -108,7 +108,7 @@ The rules today:
 
 The six rules after TAS2R38 need no phase. Their sites lie close together on one copy (the two APOL1 G1 changes were in perfect linkage disequilibrium in Genovese et al. 2010; the 9p21 and 17q21.31 tags each mark one haplotype), or, for HFE, a double heterozygote reads as one change on each copy. A copy the rule does not describe, such as a 9p21 copy with only one of the two risk alleles, is `other_haplotype`, not read. What each rule leaves out, said in its note:
 
-- APOL1 G2 (rs71785313) is a 6-base deletion. The rule table reads single-base sites only, so the rule itself reads G1, and a G2 copy reads as G0 there. The HTML and text reports read G2 beside the rule (`bin/genotype_rules_indels.py`, through `bin/collect_summary.py`) from the sample's gVCF, or its VCF when there is no gVCF, and show the G1 plus G2 result under "With G2": G1/G2, G2/G2 and G1/G1 are the high-risk genotype. Two G2 copies read G2/G2 even when the G1 sites were not read. A VCF without reference blocks cannot tell a reference base from a position with no reads, so there a result with no G2 copy is marked `(no G2 call)`. The Nextflow report stages the called VCF only, so a Nextflow run gets that mark where a run-all.sh run, which reads the gVCF, does not.
+- APOL1 G2 (rs71785313) is a 6-base deletion. The rule table reads single-base sites only, so the rule itself reads G1, and a G2 copy reads as G0 there. The HTML and text reports read G2 beside the rule (`bin/genotype_rules_indels.py`, through `bin/collect_summary.py`) from the sample's gVCF, or its VCF when there is no gVCF, and show the G1 plus G2 result under "With G2": G1/G2, G2/G2 and G1/G1 are the high-risk genotype. Two G2 copies read G2/G2 even when the G1 sites were not read. A VCF without reference blocks cannot tell a reference base from a position with no reads, so there a result with no G2 copy is marked `(no G2 call)`. The Nextflow report gets the sample's gVCF too (DeepVariant's, or the samplesheet's `gvcf` column), so both runs read the gVCF when there is one.
 - SERPINA1: M is any copy with neither S nor Z; rarer deficiency and null alleles are not read. A serum level confirms a deficiency.
 - HFE: the C282Y homozygote is also a secondary finding (`assets/acmg_sf.tsv`, [step 23](23-clinical-filter.md)); the rule gives the reading of every genotype.
 - GC: the protein levels are those of Powe et al. (2013), whose monoclonal assay was questioned in later letters. A protein type, not a vitamin D deficiency.