diff --git a/assets/output_contract.tsv b/assets/output_contract.tsv index 4eb1563f..5477c732 100644 --- a/assets/output_contract.tsv +++ b/assets/output_contract.tsv @@ -107,6 +107,7 @@ file {sample}/cpic/{sample}_cpic_recommendations.txt 27 CPIC_LOOKUP both text st file {sample}/cpic/{sample}_phenotypes.tsv 27 CPIC_LOOKUP both tsv stable both - file {sample}/cpsr/{sample}.cpsr.grch38.classification.tsv.gz 17 CPSR nextflow tsv stable - - file {sample}/cpsr/{sample}.cpsr.grch38.html 17 CPSR nextflow html stable - - +file {sample}/cpsr/{sample}_cpsr_cdna_check.tsv 17 CLINVAR_CDNA_CHECK both tsv stable - - file {sample}/cyrius/{sample}_cyp2d6.tsv 21 CYRIUS nextflow tsv stable - - file {sample}/cyrius/{sample}_cyp2d6_depth_check.tsv 21 CYRIUS nextflow tsv stable - - file {sample}/delly/{sample}_sv.vcf.gz 19 DELLY_BCF2VCF nextflow vcf stable nextflow INFO/SVTYPE diff --git a/bin/owngenome_live.py b/bin/owngenome_live.py index a28d9f6c..b8576a73 100755 --- a/bin/owngenome_live.py +++ b/bin/owngenome_live.py @@ -678,7 +678,29 @@ def run_files(self, run): return out def run_file_allowed(self, run, rel): - return rel in self.run_files(run) or any(fnmatch.fnmatch(rel, g) for g in self.app.run_globs) + files = self.run_files(run) + if rel in files or any(fnmatch.fnmatch(rel, g) for g in self.app.run_globs): + return True + # The one HTML file directly in a folder the report links (cpsr/, which + # holds CPSR's own report): the folder request redirects to it (run_path). + folder = os.path.dirname(rel) + return bool(folder) and self.sole_html(run, folder, files) == os.path.basename(rel) + + def sole_html(self, run, folder, files=None): + """The name of the only HTML file directly in a run folder the report links, inside the + run and the root; None when the folder is not linked or holds none or several.""" + d = run["dir"] + if not d or folder not in (self.run_files(run) if files is None else files): + return None + p = os.path.join(d, folder) + try: + names = os.listdir(p) if os.path.isdir(p) else [] + except OSError: + return None + top = os.path.realpath(d) + pages = [f for f in names if f.lower().endswith(".html") and os.path.isfile(os.path.join(p, f)) + and Live.inside(os.path.join(p, f), top) and Live.inside(os.path.join(p, f), self.root)] + return pages[0] if len(pages) == 1 else None def index_page(self): assets = self.cfg.assets @@ -952,6 +974,14 @@ def run_path(self, run, rest): self.log("file", name, rel, 404) return self.deny(404) p = os.path.join(d, rel) + if os.path.isdir(p) and Live.inside(p, os.path.realpath(d)) and Live.inside(p, self.live.root): + # The site lists no folder: a linked folder with one HTML file in it + # (cpsr/ with CPSR's report) opens that file. + page = self.live.sole_html(run, rel) + if page: + self.log("file", name, rel, 302) + return self.send_body(302, "text/plain", "", [("Location", "/runs/" + urllib.parse.quote(name) + "/" + + urllib.parse.quote(rel + "/" + page))]) if not os.path.isfile(p) or not Live.inside(p, os.path.realpath(d)) or not Live.inside(p, self.live.root): self.log("file", name, rel, 404) return self.deny(404) diff --git a/docs/17-cpsr.md b/docs/17-cpsr.md index f4f97027..a5d59a2b 100644 --- a/docs/17-cpsr.md +++ b/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. The Nextflow `CPSR` process does not run it: it gets no BAM, as `CLINVAR_SCREEN` does not. +- `${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. ### What the reports show diff --git a/docs/live-site.md b/docs/live-site.md index b9602775..26018b78 100644 --- a/docs/live-site.md +++ b/docs/live-site.md @@ -63,11 +63,11 @@ The reads cross the network to the browser that logged in: only the byte ranges ## Privacy and what is stored -- The site serves only its own pages, the run's state, the report, the files the report links to by a relative link, and the genome viewer's files of the run (see above). Any other file is 404, and so is a link that leaves the run folder or a symlink out of `GENOME_DIR`. It never lists a folder. +- The site serves only its own pages, the run's state, the report, the files the report links to by a relative link, and the genome viewer's files of the run (see above). Any other file is 404, and so is a link that leaves the run folder or a symlink out of `GENOME_DIR`. It never lists a folder: a link to a folder that holds one HTML file (the report's `cpsr/`, with CPSR's own report) opens that file, and any other folder link is 404. - Each message is cut to what a progress view needs before it is stored: the task, its process and status, times, CPU and memory figures. Nextflow's parameters, command line and task scripts are dropped. Three paths stay on the server for later views (the output folder, the CPU and memory caps, the work directory); none of them reaches a page. - Messages go to `/runs//events.jsonl` (state dir: `GENOME_DIR/.owngenome-live`, or `OWNGENOME_LIVE_STATE_DIR`). The site writes nowhere else, so `GENOME_DIR` can be mounted read-only. - The access log, `/access.log`, has one tab-separated row per page or file read and per login (the genome viewer's files excepted: one row per page open, see above): time, client, the first 8 characters of a hash of the session, action, run, file, status. It never holds a query string, a request body, a result or the password. It moves to `access.log.1` past `--access-log-max-mb` (10). -- Every page carries a Content-Security-Policy that lets it talk to this server only, and the page loads no script from anywhere else. The run's own HTML (the report, CPSR, PharmCAT, MultiQC) is sent sandboxed: no script runs and nothing outside is fetched. +- Every page carries a Content-Security-Policy that lets it talk to this server only, and the page loads no script from anywhere else. The run's own HTML (the report, CPSR, PharmCAT, MultiQC) is sent sandboxed: no script runs and nothing outside is fetched. CPSR's report draws its variant tables with scripts, so on the live site those tables stay empty: open the file from `/cpsr/` on disk to see them. ## The password and the limits diff --git a/docs/nextflow.md b/docs/nextflow.md index de59d6cb..b15ae58d 100644 --- a/docs/nextflow.md +++ b/docs/nextflow.md @@ -6,7 +6,7 @@ The pipeline is a [Nextflow](https://www.nextflow.io/) DSL2 pipeline, `main.nf`. - **a BAM or a CRAM** without a VCF: it is called; - **a VCF** from any caller (nf-core/sarek, DRAGEN, a provider), with an optional BAM or CRAM. -Every BAM then goes through a sex check (indexcov), and the pipeline runs pharmacogenomics, variant annotation, clinical screening, BAM analyses, structural variant calling and reporting: 8 workflows, 85 processes in 47 module files under `modules/local/`. A VCF given in the samplesheet needs FILTER=PASS records and GRCh38 contig names with chr; see [FILTER=PASS required](#filterpass-required) and [Contig names and gVCF input](#contig-names-and-gvcf-input). Starting from a provider's VCF: [Starting from a Vendor VCF](vcf-first.md). +Every BAM then goes through a sex check (indexcov), and the pipeline runs pharmacogenomics, variant annotation, clinical screening, BAM analyses, structural variant calling and reporting: 8 workflows under `workflows/`, built from the processes under `modules/local/` (one folder per tool). A VCF given in the samplesheet needs FILTER=PASS records and GRCh38 contig names with chr; see [FILTER=PASS required](#filterpass-required) and [Contig names and gVCF input](#contig-names-and-gvcf-input). Starting from a provider's VCF: [Starting from a Vendor VCF](vcf-first.md). > **Nextflow is the pipeline; the scripts are single steps.** Each numbered script in `scripts/` runs one step on its own and takes its image tags from the same `versions.env` and its helpers from `scripts/lib/common.sh`; `run-all.sh` is a launcher: it writes a one-row samplesheet and starts this pipeline with `-resume`, then runs the few script-only steps you ask for and the reports. CI runs the scripts and the pipeline on the same reads and fails when their results differ (see [Bash vs Nextflow parity](#bash-vs-nextflow-parity)). The Singularity profile is untested (see [Profiles](#profiles)). @@ -207,7 +207,8 @@ results/ │ ├── vep/ # VEP VCF, and the vcfanno-enriched VCF when a score file is set │ ├── slivar/ # Prioritized variants + compound hets │ ├── clinical/ # Clinically relevant variant subset -│ ├── cpsr/ # Cancer predisposition report +│ ├── cpsr/ # Cancer predisposition report: CPSR's HTML and classification table, and the cDNA +│ │ # read check of its P/LP rows (_cpsr_cdna_check.tsv) │ ├── roh/ # Runs of homozygosity │ ├── prs/ # Polygenic risk scores │ ├── ancestry/ # Ancestry: projection onto the panel (with --ancestry_ref) @@ -215,13 +216,22 @@ results/ │ ├── mito/ # Mitochondrial variant calls with their ClinVar and gnomAD annotation, the haplogroup │ │ # from them and haplocheck's contamination check │ ├── y_haplogroup/ # Y-chromosome haplogroup of male samples (Yleaf, opt-in: y_haplogroup) +│ ├── y_str/ # _y_str.tsv, .csv and the calls VCF: Y-STR repeat counts of male samples (opt-in: y_str) │ ├── bloodtype/ # Blood group antigens predicted from the genome (opt-in: bloodtype) │ ├── lpa/ # LPA KIV-2 repeat count from read depth (opt-in: lpa) │ ├── chip/ # _chip.tsv, _chip_depth.tsv and Mutect2's filtered calls: clonal │ │ # haematopoiesis in 13 genes (opt-in: chip) +│ ├── vntyper/ # _vntyper.tsv and VNtyper's run folder: the MUC1 VNTR screen (opt-in: vntyper) +│ ├── viral/ # _viral.tsv, the coverage of each virus and of the human genome, and the viral reads BAM (opt-in: viral) +│ ├── exomiser/ # _exomiser_top.tsv, the text summary, the HPO terms and Exomiser's own results: +│ │ # variants ranked against the symptoms (--hpo, exomiser) +│ ├── mei/ # _mei.tsv, _mei_sites.tsv, the depth table and MEGANE's VCF: mobile element insertions (opt-in: mei) +│ ├── mosaic/ # _mocha.tsv and _mocha_summary.tsv: mosaic chromosomal changes (MoChA, opt-in: mocha) +│ ├── locityper/ # _locityper.tsv and Locityper's work folder: genotypes of hard genes (opt-in: locityper) │ ├── archaic/ # Neanderthal and Denisovan share and segments (IBDmix, opt-in: archaic) │ ├── ancient/ # The nearest ancient groups (table and text), their info and the PCA plot (AADR, opt-in: ancient) │ ├── hla/ # HLA typing +│ ├── hla_optitype/ # OptiType's HLA-A and HLA-B, and _hla_concordance.tsv against T1K (--hla_qc) │ ├── expansion_hunter/ # Repeat expansion calls, Stranger's annotation and, when Stranger flags RFC1, its motif count │ ├── telomere/ # Telomere length estimation │ ├── coverage/ # Coverage statistics (mosdepth) @@ -347,7 +357,7 @@ Where a module and its script differ on purpose: | Mito variants (20) | extracts the chrM reads with `samtools view`; annotates the calls with `clinvar/clinvar_chr.vcf.gz` and the gnomAD chrM sites under `gnomad_mito/` when they are installed | with GATK `PrintReads`. Mutect2 applies its own read filters to either, so the calls are expected to match; CI does not compare them. Both mark possible NuMTs at the median autosomal depth from mosdepth (step 16b, or `MOSDEPTH` when `mosdepth` is in `--tools`). `MITO_VARIANTS_ANNOTATE` runs the same annotation (`bin/collect_summary.py mito-annotate`) with the files of `--mito_clinvar` and `--gnomad_mito`; without them the table names each file as not given and leaves its columns `.`. `run-all.sh` passes `--mito_clinvar` when `clinvar/clinvar_chr.vcf.gz` and its `.tbi` are installed, and `--gnomad_mito` when the gnomAD chrM sites are | | AnnotSV (05) | annotates step 15's duphold-filtered calls; without them, or when Manta's calls are newer, Manta's calls, and says so | `ANNOTSV` always reads `DUPHOLD_FILTER`'s output: `annotsv` in `--tools` needs `duphold` and `manta` | | Y haplogroup (37) | reads the sex step 16 (indexcov) infers | runs on the rows whose samplesheet sex is male, which `INDEXCOV` has checked against the reads | -| HTML report (24) | renders every section from `bin/collect_summary.py`'s summary | `HTML_REPORT` runs the same code on the outputs of this run's steps, among them the cDNA read check of the ClinVar hits, the CYP2D6 depth check, the phased compound-het pairs, the RFC1 motif count, AnnotSV's class 4 and 5 table, the chrM annotation and step 35's HBA1/HBA2 count. The clinical filter and slivar cards show counts only, beside slivar's pairs by phase (the module gets their VCFs, not their summary tables). For every section, run `GENOME_DIR= scripts/24-html-report.sh ` on the Nextflow output | +| HTML report (24) | renders every section from `bin/collect_summary.py`'s summary | `HTML_REPORT` runs the same code on the outputs of this run's steps, among them the cDNA read check of the ClinVar hits, the CYP2D6 depth check, the phased compound-het pairs, the RFC1 motif count, AnnotSV's class 4 and 5 table, the chrM annotation and step 35's HBA1/HBA2 count. It also gets the VEP VCF (each ClinVar hit's gnomAD AF and rsIDs), CPSR's classification table and the cDNA read check of its P/LP rows, `ALIGN_MINIMAP2`'s read count, and the clinical filter's and slivar's summary tables. The slivar card has no compound-het variant and gene counts: the Nextflow slivar module writes no `_compound_hets.tsv`. For every section, run `GENOME_DIR= scripts/24-html-report.sh ` on the Nextflow output | | CNVpytor (18) | mounts each resource file over the image's data folder | copies the files into the image's `site-packages`, so it needs a writable container | | Cyrius (21) | runs Cyrius from the install `setup.sh --cyrius` made | the same install, given as `--cyrius_install` (see [Cyrius, opt-in](#cyrius-opt-in)) | @@ -433,7 +443,7 @@ Cyrius (CYP2D6 star alleles) is under a non-commercial licence and is not in the ### CI validation scope -`nextflow.yml` lints the pipeline and stub-runs every process, on the pinned Nextflow release and the newest 26.10.x, with one samplesheet row per starting point; it also checks that `T1K_BUILD` ran once for all samples, and that a VCF+BAM-only samplesheet against a reference with no `.sr.mmi` beside it finishes without building an index or aligning anything. A job that fails prints which tasks were running or waiting and uploads the Nextflow logs as the artifact `nextflow-logs-`. The E2E workflow runs the pipeline with real containers on a slice of the public HG002 genome: from the FASTQ pair through alignment, the sex check and DeepVariant to the default tools but cyrius (the fixture's BAM lacks the autosomal bins Cyrius normalises on), plus clinvar and hla_typing; through the VCF+BAM entry, over several cases, with clinvar, mosdepth, delly, manta, vcfanno, roh, pharmcat, cpic, pypgx, telomere_hunter, mito_haplogroup, survivor_merge, y_haplogroup, sample_qc, prs, cram_archive and html_report (the prs case's row names a gVCF too); and the sex-check stop (see [Testing](testing.md#the-e2e-job)). It then runs the parity check above. It cannot run the tools whose databases do not fit a CI runner (vep's offline cache, cpsr) or the ones it has no data for. Before trusting results from a new installation, run the pipeline on a known sample and compare key outputs (PharmCAT star alleles, ClinVar hit counts) against expected values. +`nextflow.yml` lints the pipeline and stub-runs every process, on the pinned Nextflow release and the newest 26.10.x, with one samplesheet row per starting point; it also checks that `T1K_BUILD` ran once for all samples, and that a VCF+BAM-only samplesheet against a reference with no `.sr.mmi` beside it finishes without building an index or aligning anything. A job that fails prints which tasks were running or waiting and uploads the Nextflow logs as the artifact `nextflow-logs-`. The E2E workflow runs the pipeline with real containers on a slice of the public HG002 genome: from the FASTQ pair through alignment, the sex check and DeepVariant to the default tools but cyrius (the fixture's BAM lacks the autosomal bins Cyrius normalises on), plus clinvar and hla_typing; through the VCF+BAM entry, over several cases, with clinvar, mosdepth, delly, manta, vcfanno, roh, pharmcat, cpic, pypgx, telomere_hunter, mito_haplogroup, survivor_merge, y_haplogroup, sample_qc, prs, cram_archive and html_report (the prs case's row names a gVCF too), and duphold, expansion_hunter and cnvpytor (case `nf-wiring-2-sv-str-tools`); and the sex-check stop (see [Testing](testing.md#the-e2e-job)). It then runs the parity check above. It cannot run the tools whose databases do not fit a CI runner (vep's offline cache, cpsr) or the ones it has no data for. Before trusting results from a new installation, run the pipeline on a known sample and compare key outputs (PharmCAT star alleles, ClinVar hit counts) against expected values. --- diff --git a/main.nf b/main.nf index 6345bf79..7160fb9d 100644 --- a/main.nf +++ b/main.nf @@ -59,6 +59,22 @@ include { MOCHA_PREP; MOCHA_PHASE; MOCHA_CALL; MOCHA_REPORT } from './modules/lo workflow { // ─── Validate inputs ──────────────────────────────────────────────── + // Every parameter must be one nextflow_schema.json declares: Nextflow + // takes any --name as a new parameter, so a misspelt option (--outdirr, + // --max-cpus) would be ignored with no word. --pedigree is read with + // containsKey and has no default; config_profile_* are the test + // profiles' labels. + def schema_params = new groovy.json.JsonSlurper() + .parseText(file("${projectDir}/nextflow_schema.json").text)['definitions'] + .values().collectMany { group -> (group['properties'] ?: [:]).keySet() as List } + def known_params = schema_params + ['pedigree', 'config_profile_name', 'config_profile_description'] + def unknown_params = params.keySet().findAll { k -> !known_params.contains(k) }.sort() + if (unknown_params) { + error "unknown parameter${unknown_params.size() > 1 ? 's' : ''} " + + "${unknown_params.collect { k -> '--' + k }.join(', ')}: not an option of this pipeline " + + "(nextflow_schema.json lists them; options take underscores, as --max_cpus). Check the spelling." + } + if (!params.input) { error "Please provide a samplesheet with --input " } @@ -635,8 +651,8 @@ workflow { ch_parascopy_bed = Channel.value(params.parascopy_depth_bed ? file(params.parascopy_depth_bed, checkIfExists: true) : []) ch_yleaf_data = Channel.value(params.yleaf_data ? file(params.yleaf_data, checkIfExists: true) : []) // GENCODE's exons (setup.sh: reference/gencode.v.basic.exons.bed): - // the cDNA read check of CLINVAR_SCREEN's hits, and the couple screen's - // deletions over an exon + // the cDNA read check of CLINVAR_SCREEN's hits and of CPSR's P/LP rows, + // and the couple screen's deletions over an exon ch_gencode_exons = Channel.value(params.gencode_exons ? file(params.gencode_exons, checkIfExists: true) : []) // Step 20's annotation of the chrM calls (MITO_VARIANTS_ANNOTATE): the // whole ClinVar VCF with chr names, its .tbi beside it (setup.sh installs @@ -747,7 +763,9 @@ workflow { ch_reference, ch_reference_fai, BAM_ANALYSIS.out.mito_vcf, - ch_ancestry_fine_panel + ch_ancestry_fine_panel, + ch_bam, + ch_gencode_exons ) // ─── Traits panel (traits, step 42) ───────────────────────────────── @@ -1118,10 +1136,17 @@ workflow { .join(PGX.out.pharmcat_html.map { meta, f -> [meta.id, f] }, remainder: true) .join(PGX.out.cpic_phenotypes.map { meta, f -> [meta.id, f] }, remainder: true) .join(PGX.out.cpic_recommendations.map { meta, f -> [meta.id, f] }, remainder: true) + .join(ANNOTATION.out.vep_vcf.map { meta, f -> [meta.id, f] }, remainder: true) + .join(ANNOTATION.out.vep_vcf_index.map { meta, f -> [meta.id, f] }, remainder: true) .join(ANNOTATION.out.clinical_vcf.map { meta, f -> [meta.id, f] }, remainder: true) + .join(ANNOTATION.out.clinical_summary.map { meta, f -> [meta.id, f] }, remainder: true) .join(ANNOTATION.out.slivar_vcf.map { meta, f -> [meta.id, f] }, remainder: true) + .join(ANNOTATION.out.slivar_summary.map { meta, f -> [meta.id, f] }, remainder: true) .join(ANNOTATION.out.comphet_pairs.map { meta, f -> [meta.id, f] }, remainder: true) .join(CLINICAL.out.cpsr_html.map { meta, f -> [meta.id, f] }, remainder: true) + .join(CLINICAL.out.cpsr_tsv.map { meta, f -> [meta.id, f] }, remainder: true) + .join(CLINICAL.out.cpsr_cdna.map { meta, f -> [meta.id, f] }, remainder: true) + .join(UPSTREAM.out.read_count.map { meta, f -> [meta.id, f] }, remainder: true) .join(CLINICAL.out.roh_regions.map { meta, f -> [meta.id, f] }, remainder: true) .join(CLINICAL.out.haplogroup.map { meta, f -> [meta.id, f] }, remainder: true) .join(CLINICAL.out.haplocheck.map { meta, f -> [meta.id, f] }, remainder: true) diff --git a/modules/local/clinvar_screen/main.nf b/modules/local/clinvar_screen/main.nf index bcc0db14..6d6afc82 100644 --- a/modules/local/clinvar_screen/main.nf +++ b/modules/local/clinvar_screen/main.nf @@ -263,37 +263,55 @@ process CLINVAR_SCREEN { _clinvar_cdna_scan.tsv, the loci with reads across two or more exon junctions; ext.cdna_scan = false skips it (the table then says so). - Equivalent to: scripts/06-clinvar-screen.sh, step 4 + source 'cpsr' runs the same check on CPSR's P/LP rows instead (the + CLINICAL workflow): hits_tsv is CPSR's classification table, read with + clip_check.py --cpsr, and the task writes cpsr/_cpsr_cdna_check.tsv + and no scan, as scripts/17 does. construct is [] there, as scripts/17 + passes no construct catalog. + + Equivalent to: scripts/06-clinvar-screen.sh, step 4, and + scripts/17-cpsr.sh, the cDNA read check */ process CLINVAR_CDNA_CHECK { tag "$meta.id" label 'process_single' - // The table only: CLINVAR_SCREEN's versions.yml is the one in clinvar/. + // The tables only: CLINVAR_SCREEN's (CPSR's) versions.yml is the one in clinvar/ (cpsr/). publishDir { "${params.outdir}/${meta.id}/clinvar" }, mode: params.publish_dir_mode, pattern: '*_clinvar_cdna_*.tsv' + publishDir { "${params.outdir}/${meta.id}/cpsr" }, mode: params.publish_dir_mode, pattern: '*_cpsr_cdna_check.tsv' input: // Fixed names: a samplesheet path may hold any character, and these go on a command line. tuple val(meta), path(hits_tsv), path(bam, stageAs: 'sample.bam'), path(bai, stageAs: 'sample.bam.bai') path(exons, stageAs: 'exons.bed') path(construct, stageAs: 'construct_variants.tsv') + val(source) // 'clinvar' (step 06's hits) or 'cpsr' (step 17's classification table) output: - tuple val(meta), path("${meta.id}_clinvar_cdna_check.tsv"), emit: cdna_check - tuple val(meta), path("${meta.id}_clinvar_cdna_scan.tsv"), emit: cdna_scan + // One source per task: the clinvar tables, or the cpsr one. + tuple val(meta), path("${meta.id}_clinvar_cdna_check.tsv"), emit: cdna_check, optional: true + tuple val(meta), path("${meta.id}_clinvar_cdna_scan.tsv"), emit: cdna_scan, optional: true + tuple val(meta), path("${meta.id}_cpsr_cdna_check.tsv"), emit: cpsr_cdna_check, optional: true path "versions.yml", emit: versions when: task.ext.when == null || task.ext.when script: + if (!(source in ['clinvar', 'cpsr'])) { + error "CLINVAR_CDNA_CHECK: source must be 'clinvar' or 'cpsr', not '${source}'" + } def reads = !bam ? "--skip 'no BAM for this sample in the samplesheet'" : !exons ? "--skip 'no GENCODE exons (--gencode_exons)'" : "--bam ${bam} --exons ${exons}" def scan = task.ext.cdna_scan == false ? "--skip 'ext.cdna_scan = false'" : reads + def rows = source == 'cpsr' ? "--cpsr ${hits_tsv}" : "--hits ${hits_tsv}" + def cons = construct ? "--construct ${construct}" : '' + // The scan reads every exon, whatever the rows: once per sample, with step 06. + def scan_cmd = source == 'clinvar' ? "cdna_scan.py --out ${meta.id}_clinvar_cdna_scan.tsv ${scan}" : '' """ - clip_check.py --hits ${hits_tsv} --out ${meta.id}_clinvar_cdna_check.tsv --construct ${construct} ${reads} - cdna_scan.py --out ${meta.id}_clinvar_cdna_scan.tsv ${scan} + clip_check.py ${rows} --out ${meta.id}_${source}_cdna_check.tsv ${cons} ${reads} + ${scan_cmd} cat <<-END_VERSIONS > versions.yml "${task.process}": @@ -303,7 +321,8 @@ process CLINVAR_CDNA_CHECK { stub: """ - touch ${meta.id}_clinvar_cdna_check.tsv ${meta.id}_clinvar_cdna_scan.tsv + touch ${meta.id}_${source}_cdna_check.tsv + ${source == 'clinvar' ? "touch ${meta.id}_clinvar_cdna_scan.tsv" : ''} cat <<-END_VERSIONS > versions.yml "${task.process}": diff --git a/modules/local/cpsr/main.nf b/modules/local/cpsr/main.nf index 2f8a5895..cbf0e538 100644 --- a/modules/local/cpsr/main.nf +++ b/modules/local/cpsr/main.nf @@ -12,9 +12,9 @@ and its files are renamed back to .cpsr.*. CPSR's own HTML title then shows that id. - Unlike scripts/17, this process does not run the cDNA read check on - CPSR's P/LP rows (bin/clip_check.py --cpsr): it gets no BAM, as - CLINVAR_SCREEN does not for step 06's check. + The cDNA read check on CPSR's P/LP rows (bin/clip_check.py --cpsr) runs + after it, in CLINVAR_CDNA_CHECK with source 'cpsr' (workflows/clinical.nf), + which gets the BAM and the exons and writes cpsr/_cpsr_cdna_check.tsv. Equivalent to: scripts/17-cpsr.sh ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/modules/local/html_report/main.nf b/modules/local/html_report/main.nf index baf425a3..8473e5f6 100644 --- a/modules/local/html_report/main.nf +++ b/modules/local/html_report/main.nf @@ -10,7 +10,9 @@ The inputs are the outputs of the steps that ran for this sample, as one list (main.nf joins them, so the report waits for each). They are linked into the folder layout collect_summary.py reads (clinvar/, pharmcat/, - cpic/, roh/, mito/, telomere/, hla_t1k/, hla_optitype/, kir/, pypgx/, annotsv/, traits/, chip/, couple/, ...); a step that did + cpic/, vep/, clinical/, slivar/, cpsr/, qc/, roh/, mito/, telomere/, + hla_t1k/, hla_optitype/, kir/, pypgx/, annotsv/, traits/, chip/, couple/, + ...); a step that did not run, or whose task failed and was ignored, has no file, and its card says "Not run". The called or given VCF comes separately, because its name is the user's. indexcov's folder becomes indexcov/ itself (the @@ -65,16 +67,17 @@ process HTML_REPORT { ${id}_clinvar_hits.vcf|${id}_clinvar_cdna_check.tsv|${id}_clinvar_cdna_scan.tsv|${id}_clinvar_deep_intronic.tsv) place clinvar "\$f" ;; ${id}.report.json|${id}.report.html) place pharmcat "\$f" ;; ${id}_phenotypes.tsv|${id}_cpic_recommendations.txt) place cpic "\$f" ;; - ${id}_clinical.vcf.gz) place clinical "\$f" ;; - ${id}.cpsr.grch38.html) place cpsr "\$f" ;; - ${id}_prioritized.vcf.gz|${id}_compound_het_pairs.tsv) place slivar "\$f" ;; + ${id}_vep.vcf.gz|${id}_vep.vcf.gz.tbi) place vep "\$f" ;; + ${id}_clinical.vcf.gz|${id}_clinical_summary.tsv) place clinical "\$f" ;; + ${id}.cpsr.grch38.html|${id}.cpsr.grch38.classification.tsv.gz|${id}_cpsr_cdna_check.tsv) place cpsr "\$f" ;; + ${id}_prioritized.vcf.gz|${id}_compound_het_pairs.tsv|${id}_slivar_summary.tsv) place slivar "\$f" ;; ${id}_roh.txt) place roh "\$f" ;; ${id}_haplogroup.txt|${id}_haplocheck.txt) place mito "\$f" ;; ${id}_y_haplogroup.txt|${id}_y_tree.txt|${id}_y_tree_path.tsv) place y_haplogroup "\$f" ;; ${id}_y_str.tsv) place y_str "\$f" ;; ${id}_paraphase_copy_number.tsv|${id}_paraphase_variants.tsv|${id}_*_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}_sample_qc.tsv|${id}_read_count.tsv) place qc "\$f" ;; ${id}_cn_intervals.tsv) place paralogs "\$f" ;; ${id}_prs_summary.tsv) place prs "\$f" ;; ${id}_ancestry.tsv) place ancestry "\$f" ;; diff --git a/scripts/ci/e2e-case-paths.tsv b/scripts/ci/e2e-case-paths.tsv index 0c1fc6ce..5dc19dc9 100644 --- a/scripts/ci/e2e-case-paths.tsv +++ b/scripts/ci/e2e-case-paths.tsv @@ -110,6 +110,8 @@ nextflow-from-fastq-3-sexcheck assets/callability_genes.tsv assets/live/index.ht nextflow-hardening-1-manta assets/par_grch38.bed assets/sv_rule_regions.tsv bin/annotsv_postfilter.py bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/manta_inv_filter.py bin/sv_svlen_fix.py scripts/04-manta.sh scripts/05-annotsv.sh scripts/15-duphold.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/ci/fake-docker/lib.sh nextflow-hardening-2-nf assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/construct_variants.tsv assets/deep_intronic_sites.tsv assets/live/index.html assets/mito_confirmed.tsv assets/par_grch38.bed assets/pgx_genes.tsv assets/stub/sample.bam assets/stub/sample.bam.bai assets/stub/sample.vcf.gz assets/sv_rule_regions.tsv assets/viewer/index.html bin/annotsv_postfilter.py bin/cdna_scan.py bin/clinvar_hits.awk bin/clinvar_join_adjacent.awk bin/clip_check.py bin/cn_loci.tsv bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/copy_number.py bin/cyp2d6_depth_check.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/manta_inv_filter.py bin/pgx_mito_merge.py bin/pgx_outside_calls.py bin/pgx_refblock_fill.py bin/render_report.py bin/results_index.py bin/sv_svlen_fix.py bin/y_tree_place.py bin/yleaf_run.py modules/local/clinvar_screen/* modules/local/delly/* modules/local/html_report/* modules/local/indexcov/* modules/local/manta/* modules/local/mito_haplogroup/* modules/local/mosdepth/* modules/local/pgx_consensus/* modules/local/pharmcat/* modules/local/pypgx/* modules/local/survivor_merge/* modules/local/telomere_hunter/* modules/local/vcf_precheck/* modules/local/vcfanno/* modules/local/y_haplogroup/* scripts/04-manta.sh scripts/05-annotsv.sh scripts/06-clinvar-screen.sh scripts/07-pharmacogenomics.sh scripts/10-telomere-hunter.sh scripts/12-mito-haplogroup.sh scripts/15-duphold.sh scripts/16-indexcov.sh scripts/16b-mosdepth.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/20-mtoolbox.sh scripts/21-cyrius.sh scripts/22-survivor-merge.sh scripts/24-html-report.sh scripts/30-vcfanno.sh scripts/32-pypgx.sh scripts/36-pgx-consensus.sh scripts/37-y-haplogroup.sh scripts/ci/fake-docker/lib.sh scripts/ci/parity-diff.sh scripts/cyrius-constraints.txt scripts/setup.sh scripts/validate-setup.sh nextflow-report-parity-1-sections assets/acmg_sf.tsv assets/callability_genes.tsv assets/carrier_genes.tsv assets/construct_variants.tsv assets/deep_intronic_sites.tsv assets/mito_confirmed.tsv assets/par_grch38.bed assets/pgx_genes.tsv assets/str_extra_loci.json assets/stub/db/catalog.json assets/stub/sample.bam assets/stub/sample.bam.bai assets/stub/sample.vcf.gz assets/sv_rule_regions.tsv bin/annotsv_postfilter.py bin/cdna_scan.py bin/clinvar_hits.awk bin/clinvar_join_adjacent.awk bin/clip_check.py bin/cn_loci.tsv bin/cnvpytor_to_vcf.py bin/collect_summary.py bin/comphet_phase.py bin/constraint_join.awk bin/copy_number.py bin/csq_rank.awk bin/cyp2d6_depth_check.py bin/export_consumer.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/pgx_outside_calls.py bin/render_report.py bin/results_index.py bin/rfc1_motif.py bin/strchive_convert.py bin/sv_svlen_fix.py modules/local/annotsv/* modules/local/clinvar_screen/* modules/local/expansion_hunter/* modules/local/export_consumer/* modules/local/html_report/* modules/local/mito_variants/* modules/local/pypgx/* modules/local/slivar/* modules/local/stranger/* modules/local/vcf_precheck/* modules/local/vep/* scripts/05-annotsv.sh scripts/06-clinvar-screen.sh scripts/09-expansion-hunter.sh scripts/09b-stranger.sh scripts/13-vep-annotation.sh scripts/15-duphold.sh scripts/16b-mosdepth.sh scripts/18-cnvpytor.sh scripts/19-delly.sh scripts/20-mtoolbox.sh scripts/24-html-report.sh scripts/31-slivar.sh scripts/32-pypgx.sh scripts/36-pgx-consensus.sh scripts/40-export-consumer.sh scripts/ci/fake-docker/lib.sh scripts/validate-setup.sh tests/fixtures/consumer_export/sites.tsv.gz +nf-wiring-1-report-tables assets/acmg_sf.tsv assets/test_samplesheet.csv assets/stub/* bin/clip_check.py bin/collect_summary.py bin/cpsr_sample_id bin/render_report.py bin/results_index.py conf/test_all.config conf/test.config modules/local/align_minimap2/* modules/local/clinical_filter/* modules/local/clinvar_screen/* modules/local/cpsr/* modules/local/html_report/* modules/local/slivar/* modules/local/vep/* scripts/24-html-report.sh tests/fixtures/consumer_export/sites.tsv.gz +nf-wiring-2-sv-str-tools bin/cnvpytor_to_vcf.py bin/manta_inv_filter.py bin/strchive_convert.py bin/sv_svlen_fix.py modules/local/cnvpytor/* modules/local/duphold/* modules/local/expansion_hunter/* modules/local/indexcov/* modules/local/manta/* modules/local/vcf_precheck/* output-contract-1-hg002 assets/* bin/* modules/* scripts/* paralog-catalog-1-gba assets/paralog_loci.tsv assets/parascopy_models_v1.7.txt bin/cn_loci.tsv bin/copy_number.py bin/hba_copy_number.py bin/lpa_kiv2.py modules/local/parascopy/* modules/local/vcf_precheck/* scripts/35-paralogs.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh scripts/validate-setup.sh paralog-checks-and-doc-fixes-1-hba assets/par_grch38.bed assets/paralog_loci.tsv bin/cn_loci.tsv bin/collect_summary.py bin/copy_number.py bin/hba_copy_number.py bin/lpa_kiv2.py bin/paraphase_report.py modules/local/align_minimap2/* modules/local/deepvariant/* modules/local/fastp/* modules/local/paraphase/* modules/local/parascopy/* modules/local/vcf_precheck/* scripts/01b-fastp-qc.sh scripts/35-paralogs.sh scripts/35b-paraphase.sh scripts/ci/check-images.sh scripts/ci/fake-docker/lib.sh scripts/setup.sh scripts/validate-setup.sh diff --git a/scripts/ci/e2e-case-times.tsv b/scripts/ci/e2e-case-times.tsv index ef818c27..6ad857e5 100644 --- a/scripts/ci/e2e-case-times.tsv +++ b/scripts/ci/e2e-case-times.tsv @@ -109,6 +109,8 @@ nextflow-from-fastq-3-sexcheck shard 97 - nextflow-hardening-1-manta shard 21 - nextflow-hardening-2-nf shard 643 bash-step-19-delly,nextflow-hardening-1-manta nextflow-report-parity-1-sections shard 206 clinvar-clip-check-1,comphet-phasing-1-reads +nf-wiring-1-report-tables shard 300 clinvar-clip-check-1 +nf-wiring-2-sv-str-tools shard 900 - output-contract-1-hg002 shard 30 - paralog-catalog-1-gba shard 114 - paralog-checks-and-doc-fixes-1-hba shard 57 - diff --git a/tests/e2e/nf-wiring-1-report-tables.sh b/tests/e2e/nf-wiring-1-report-tables.sh new file mode 100755 index 00000000..af473ba2 --- /dev/null +++ b/tests/e2e/nf-wiring-1-report-tables.sh @@ -0,0 +1,302 @@ +#!/usr/bin/env bash +# The Nextflow report gets every table the bash report reads, and CPSR's cDNA +# read check runs in Nextflow too. +# 1. CLINVAR_CDNA_CHECK with source 'cpsr', run for real (PYPGX_IMAGE), on a +# CPSR table holding the planted hit as a P/LP row and a benign row, with +# the spliced reads clinvar-clip-check-1 left: it publishes +# cpsr/_cpsr_cdna_check.tsv, the P/LP row is flagged with the counts +# bin/clip_check.py --cpsr gives step 17, and the benign row is not +# checked. Control: no exon file, and the row is not checked. +# 2. Sample NFW: step 24 and HTML_REPORT on the same tables: the ClinVar +# hits, a VEP file with a gnomAD AF on each hit (3% on the planted one), +# CPSR's table and its cDNA read check, step 02's read count and the +# clinical filter's and slivar's summary tables. The read_count, clinvar, +# cpsr, clinical and slivar sections of the two summary.json files must +# be equal, with the AF, the CPSR classes and the flagged row in them. +# Control: the module as it was (no place for these tables) stops. +# 3. main.nf passes them on: a stub run of every tool (-profile test_all) +# stages each table in each HTML_REPORT task whose sample has the step's +# main output. A misspelt option (--outdirr) stops the run and names it. +# Control: a copy without the new joins and without the option check +# stages none and takes --outdirr in silence. +. "$(dirname "$0")/lib.sh" + +command -v nextflow >/dev/null || { fail "nextflow is not on PATH"; finish; } + +S=NFW +SD="${GENOME_DIR}/${S}" +C="${SAMPLE}cdna" +CHR=$(planted chrom) +POS=$(planted pos) +REF=$(planted ref) +ALT=$(planted alt) +HG_VCF="${GENOME_DIR}/${SAMPLE}/vcf/${SAMPLE}.vcf.gz" +HITS_VCF="${GENOME_DIR}/${SAMPLE}/clinvar/${SAMPLE}_clinvar_hits.vcf" +HITS_TSV="${GENOME_DIR}/${SAMPLE}/clinvar/${SAMPLE}_clinvar_hits.tsv" +SPLICED="${GENOME_DIR}/${C}/aligned/${C}_sorted.bam" +for f in "$HG_VCF" "$HITS_VCF" "$HITS_TSV" "$SPLICED" "${SPLICED}.bai" \ + "${GENOME_DIR}/${SAMPLE}/qc/${SAMPLE}_read_count.tsv" \ + "${GENOME_DIR}/${SAMPLE}/clinical/${SAMPLE}_clinical_summary.tsv" \ + "${GENOME_DIR}/${SAMPLE}/clinical/${SAMPLE}_clinical.vcf.gz" \ + "${GENOME_DIR}/${SAMPLE}/slivar/${SAMPLE}_slivar_summary.tsv" \ + "${GENOME_DIR}/${SAMPLE}/slivar/${SAMPLE}_prioritized.vcf.gz"; do + check "an earlier case left ${f#"$GENOME_DIR"/}" test -s "$f" +done + +# A small pipeline in its own project folder, which gets the repository's bin/ +NP="${CASE_TMP}/project" +rm -rf "$NP" +mkdir -p "$NP" +cp -R "${REPO}/bin" "${NP}/bin" +# nf_run NAME: run NP/main.nf with the real containers, output under CASE_TMP/NAME +nf_run() { + local name=$1 out="${CASE_TMP}/$1" + rm -rf "$out" + { + echo "params.outdir = '${out}/results'" + echo "params.publish_dir_mode = 'copy'" + echo "docker.enabled = true" + echo "docker.runOptions = '--user $(id -u):$(id -g) --network none'" + } > "${NP}/nextflow.config" + echo "+ nextflow run (${name})" + ( + cd "$CASE_TMP" && nextflow run "${NP}/main.nf" -ansi-log false \ + -c "${REPO}/conf/containers.config" -work-dir "${CASE_TMP}/work-${name}" + ) > "${CASE_TMP}/${name}.log" 2>&1 + NF_RC=$? +} + +# --- 1. CLINVAR_CDNA_CHECK on CPSR's rows ----------------------------------------- +mkdir -p "${CASE_TMP}/in" +CLS="${CASE_TMP}/in/${S}.cpsr.grch38.classification.tsv.gz" +{ + printf 'VAR_ID\tGENOMIC_CHANGE\tSYMBOL\tCLINVAR_CLASSIFICATION\tCPSR_CLASSIFICATION\tCLASSIFICATION\n' + printf '%s_%s_%s_%s\t%s:g.%s%s>%s\t%s\tPathogenic\tPathogenic\tPathogenic\n' "${CHR#chr}" "$POS" "$REF" "$ALT" \ + "${CHR#chr}" "$POS" "$REF" "$ALT" "$(planted gene)" + printf '%s_%s_A_G\t%s:g.%sA>G\tOTHER\tBenign\tBenign\tBenign\n' "${CHR#chr}" "$((POS + 500))" "${CHR#chr}" "$((POS + 500))" +} | gzip -c > "$CLS" +# The exon file of clinvar-cpsr-quality-1: the planted hit's exon and the next one. +EXONS="${CASE_TMP}/in/exons.bed" +P0=$((POS - 1)) +{ + printf '%s\t%s\t%s\n' "$CHR" "$((P0 - 30))" "$((P0 + 30))" + printf '%s\t%s\t%s\n' "$CHR" "$((P0 + 3000))" "$((P0 + 3150))" +} > "$EXONS" +cdna_nf() { # cdna_nf NAME EXONS-EXPRESSION + cat > "${NP}/main.nf" < /dev/null +ROW=$(awk -F'\t' -v p="$POS" 'NR > 1 && $2 == p' "$NF_CDNA" 2>/dev/null) +echo "Nextflow CPSR row: ${ROW}" +check_eq "only the P/LP row is in the table" "$(awk 'NR > 1' "$NF_CDNA" 2>/dev/null | wc -l | tr -d ' ')" 1 +check_eq "the planted P/LP row is flagged" "$(cut -f9 <<< "$ROW")" flagged +check "the table is the one bin/clip_check.py --cpsr writes for step 17, byte for byte" cmp -s "${CASE_TMP}/in/step17.tsv" "$NF_CDNA" +cdna_nf cdna-no-exons "[]" +OFF=$(awk -F'\t' -v p="$POS" 'NR > 1 && $2 == p' "${CASE_TMP}/cdna-no-exons/results/${S}/cpsr/${S}_cpsr_cdna_check.tsv" 2>/dev/null) +check_eq "control: without exons the row is not checked, so the flagged check can fail" "$(cut -f9 <<< "$OFF")" not_checked + +# --- 2. step 24 and HTML_REPORT on the same tables --------------------------------- +rm -rf "$SD" +mkdir -p "${SD}/vcf" "${SD}/clinvar" "${SD}/vep" "${SD}/cpsr" "${SD}/qc" "${SD}/clinical" "${SD}/slivar" +ln -f "$HG_VCF" "${SD}/vcf/${S}.vcf.gz" +ln -f "${HG_VCF}.tbi" "${SD}/vcf/${S}.vcf.gz.tbi" +cp "$HITS_VCF" "${SD}/clinvar/${S}_clinvar_hits.vcf" +cp "$CLS" "${SD}/cpsr/" +cp "$NF_CDNA" "${SD}/cpsr/" +cp "${GENOME_DIR}/${SAMPLE}/qc/${SAMPLE}_read_count.tsv" "${SD}/qc/${S}_read_count.tsv" +cp "${GENOME_DIR}/${SAMPLE}/clinical/${SAMPLE}_clinical_summary.tsv" "${SD}/clinical/${S}_clinical_summary.tsv" +cp "${GENOME_DIR}/${SAMPLE}/clinical/${SAMPLE}_clinical.vcf.gz" "${SD}/clinical/${S}_clinical.vcf.gz" +cp "${GENOME_DIR}/${SAMPLE}/slivar/${SAMPLE}_slivar_summary.tsv" "${SD}/slivar/${S}_slivar_summary.tsv" +cp "${GENOME_DIR}/${SAMPLE}/slivar/${SAMPLE}_prioritized.vcf.gz" "${SD}/slivar/${S}_prioritized.vcf.gz" +# A VEP file with one record per hit: a gnomAD AF on each, 3% on the planted one. +python3 - "$HITS_TSV" "${SD}/vep/plain.vcf" "$CHR" "$POS" <<'PY' +import sys +hits, out, pchrom, ppos = sys.argv[1], sys.argv[2], sys.argv[3], int(sys.argv[4]) +rows = [l.rstrip("\n").split("\t") for l in open(hits)][1:] +recs = {} +for i, r in enumerate(rows): + chrom, pos, ref, alt = r[0], int(r[1]), r[2], r[3] + af = 0.03 if (chrom, pos) == (pchrom, ppos) else round(0.0001 * (i + 1), 6) + snv = len(ref) == 1 and len(alt) == 1 + allele = alt if snv or alt[:1] != ref[:1] else (alt[1:] or "-") + recs[(chrom, pos, ref, alt)] = f"{chrom}\t{pos}\trs{1000 + i}\t{ref}\t{alt}\t50\tPASS\tCSQ={allele}|x|LOW|G|{af}|{af / 2}|rs{1000 + i}\tGT\t0/1" +order = {c: i for i, c in enumerate(sorted({k[0] for k in recs}, key=lambda c: (len(c), c)))} +with open(out, "w") as f: + f.write("##fileformat=VCFv4.2\n") + for c in order: + f.write(f"##contig=\n") + f.write('##INFO=\n') + f.write('##FORMAT=\n') + f.write("#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tS\n") + for k in sorted(recs, key=lambda k: (order[k[0]], k[1])): + f.write(recs[k] + "\n") +print(f"{len(recs)} VEP records") +PY +in_user() { docker run --rm -u "$(id -u):$(id -g)" -v "${GENOME_DIR}:/genome" -w /genome "$BCFTOOLS_IMAGE" "$@"; } +in_user bcftools view --no-version -Oz -o "${S}/vep/${S}_vep.vcf.gz" "${S}/vep/plain.vcf" +in_user bcftools index -t "${S}/vep/${S}_vep.vcf.gz" +rm -f "${SD}/vep/plain.vcf" +check "the VEP file is bgzipped and indexed" test -s "${SD}/vep/${S}_vep.vcf.gz.tbi" +run_step 24-html-report.sh "$S" +check_step_exit 24-html-report.sh + +# The report assets main.nf passes HTML_REPORT (report_assets) +ASSET_FILES=$(for n in acmg_sf.tsv carrier_genes.tsv callability_genes.tsv mito_confirmed.tsv paralog_loci.tsv paraphase_genes.tsv; do + printf "file('%s/assets/%s'), " "$REPO" "$n"; done | sed 's/, $//') +report_nf() { # report_nf NAME MODULE + cat > "${NP}/main.nf" </dev/null; } +check_eq "the ClinVar hits' frequencies come from the VEP file" \ + "$(nf_py "d['clinvar']['data']['frequency']['source']")" "vep/${S}_vep.vcf.gz" +check_eq "the planted hit has its 3% AF and is marked common" \ + "$(nf_py "[(h['af'], h['common']) for h in d['clinvar']['data']['hits'] if h['chrom'] == '${CHR}' and h['pos'] == ${POS}]")" \ + "[(0.03, True)]" +check_eq "the hits carry VEP's rsIDs" \ + "$(nf_py "sum(1 for h in d['clinvar']['data']['hits'] if h['rsids'])")" "$(awk 'NR > 1' "$HITS_TSV" | wc -l | tr -d ' ')" +check_eq "the CPSR card counts the classes" "$(nf_py "d['cpsr']['data']['classification']")" "{'Pathogenic': 1, 'Benign': 1}" +check_eq "and flags the P/LP row the read check flagged" "$(nf_py "d['cpsr']['data']['cdna_flagged']")" 1 +check_eq "the read count section is there" "$(nf_py "d['read_count']['data']['read_count_check']")" \ + "$(awk -F'\t' '$1 == "read_count_check" {print $2}' "${SD}/qc/${S}_read_count.tsv")" +check_eq "the clinical card counts the summary table's rows" "$(nf_py "d['clinical']['data']['variants']")" \ + "$(awk 'NR > 1' "${SD}/clinical/${S}_clinical_summary.tsv" | wc -l | tr -d ' ')" +check "and has its impact breakdown" test -n "$(nf_py "d['clinical']['data']['by_impact']")" +check_eq "the slivar card counts the summary table's rows" "$(nf_py "d['slivar']['data']['prioritized']")" \ + "$(awk 'NR > 1' "${SD}/slivar/${S}_slivar_summary.tsv" | wc -l | tr -d ' ')" +# Control: the module as it was, without a place for the new tables. +OLD="${CASE_TMP}/old-module/main.nf" +mkdir -p "$(dirname "$OLD")" +sed -e '/_vep\.vcf\.gz|/d' -e 's/|\${id}_clinical_summary\.tsv//' -e 's/|\${id}_slivar_summary\.tsv//' \ + -e 's/|\${id}\.cpsr\.grch38\.classification\.tsv\.gz|\${id}_cpsr_cdna_check\.tsv//' \ + -e 's/|\${id}_read_count\.tsv//' "${REPO}/modules/local/html_report/main.nf" > "$OLD" +check_eq "control copy: the module names none of the new tables" \ + "$(grep -cE 'vep\.vcf|clinical_summary|slivar_summary|classification\.tsv|cpsr_cdna_check|read_count' "$OLD")" 0 +report_nf report-old "${OLD%.nf}" +check "control: the module without the places stops on the first new table" test "$NF_RC" -ne 0 +check "control: and says it has no place for it" grep -q 'HTML_REPORT got an input it has no place for' -r "${CASE_TMP}/work-report-old" --include=.command.err + +# --- 3. main.nf passes every table on (stub run of every tool) --------------------- +# PAIRS: each new table and the step output (or published file, with "out:") +# that says the step ran for the sample +PAIRS="_vep.vcf.gz:_clinical.vcf.gz _vep.vcf.gz.tbi:_clinical.vcf.gz _clinical_summary.tsv:_clinical.vcf.gz \ +_slivar_summary.tsv:_prioritized.vcf.gz .cpsr.grch38.classification.tsv.gz:.cpsr.grch38.html \ +_cpsr_cdna_check.tsv:.cpsr.grch38.html _read_count.tsv:out:qc" +# stub_staged DIR: one line per pair, " " +stub_staged() { + local d id p new old hits got + for p in $PAIRS; do + new=${p%%:*} old=${p#*:} hits=0 got=0 + while IFS= read -r d; do + d=$(dirname "$d") + id=$(sed -n "s/^### name: [\"'].*HTML_REPORT (\(.*\))[\"']\$/\1/p" "${d}/.command.run") + case "$old" in + out:*) [ -e "$1/out/${id}/${old#out:}/${id}${new}" ] || continue ;; + *) [ -e "${d}/input/${id}${old}" ] || continue ;; + esac + hits=$((hits + 1)) + [ ! -e "${d}/input/${id}${new}" ] || got=$((got + 1)) + done < <(grep -l "^### name: [\"'].*HTML_REPORT (" -r "$1/work" --include=.command.run 2>/dev/null) + echo "$new $hits $got" + done +} +stub_run() { # stub_run NAME [OPTION...] + local d="${CASE_TMP}/$1"; shift + echo "+ nextflow run main.nf -profile test_all -stub ($(basename "$d") $*)" + (cd "$d" && nextflow run main.nf -profile test_all -stub -ansi-log false --outdir "${d}/out" -work-dir "${d}/work" "$@") \ + > "${d}.log" 2>&1 +} +for n in stub stub-control; do + rm -rf "${CASE_TMP:?}/${n}" + mkdir -p "${CASE_TMP}/${n}" + cp -R "${REPO}/main.nf" "${REPO}/nextflow.config" "${REPO}/nextflow_schema.json" "${REPO}/versions.env" \ + "${REPO}/conf" "${REPO}/modules" "${REPO}/workflows" "${REPO}/bin" "${REPO}/assets" "${REPO}/scripts" "${CASE_TMP}/${n}/" + mkdir -p "${CASE_TMP}/${n}/tests/fixtures" # test_all's consumer_sites + cp -R "${REPO}/tests/fixtures/consumer_export" "${CASE_TMP}/${n}/tests/fixtures/" +done +JOINS='ANNOTATION.out.vep_vcf.map|ANNOTATION.out.vep_vcf_index.map|ANNOTATION.out.clinical_summary.map|ANNOTATION.out.slivar_summary.map|CLINICAL.out.cpsr_tsv.map|CLINICAL.out.cpsr_cdna.map|UPSTREAM.out.read_count.map' +check_eq "main.nf joins the seven tables into the report inputs" "$(grep -cE "$JOINS" "${CASE_TMP}/stub/main.nf")" 7 +sed -i -E "/${JOINS}/d" "${CASE_TMP}/stub-control/main.nf" +sed -i 's/^ if (unknown_params) {$/ if (false) {/' "${CASE_TMP}/stub-control/main.nf" +check_eq "control copy: the joins are gone" "$(grep -cE "$JOINS" "${CASE_TMP}/stub-control/main.nf")" 0 +check_eq "control copy: the option check is off" "$(grep -c '^ if (false) {$' "${CASE_TMP}/stub-control/main.nf")" 1 +stub_run stub +check_eq "stub run exits 0" "$?" 0 +tail -n 5 "${CASE_TMP}/stub.log" +stub_run stub-control --outdirr x +check_eq "control: the stub run without the option check takes --outdirr and exits 0" "$?" 0 +while read -r new hits got; do + check_ge "stub run: tasks with the step before ${new}" "$hits" 1 + check_eq "stub run: each of them gets ${new}" "$got" "$hits" +done < <(stub_staged "${CASE_TMP}/stub") +while read -r new hits got; do + check_ge "control: tasks with the step before ${new}" "$hits" 1 + check_eq "control: none gets ${new}" "$got" 0 +done < <(stub_staged "${CASE_TMP}/stub-control") +check_ge "stub run: CPSR's read check ran for a sample" \ + "$(find "${CASE_TMP}/stub/out" -path '*/cpsr/*_cpsr_cdna_check.tsv' | wc -l | tr -d ' ')" 1 +rm -rf "${CASE_TMP:?}/stub/work" "${CASE_TMP:?}/stub/out" +stub_run stub --outdirr x --max-cpus 2 +check "a misspelt option stops the run" test "$?" -ne 0 +check "and the message names --outdirr" grep -q 'unknown parameters .*--outdirr' "${CASE_TMP}/stub.log" +check "and the kebab-case --max-cpus" grep -Eq 'unknown parameters .*--(max-cpus|maxCpus)' "${CASE_TMP}/stub.log" +grep -m1 'unknown parameter' "${CASE_TMP}/stub.log" + +rm -rf "$SD" "${CASE_TMP:?}/stub" "${CASE_TMP:?}/stub-control" "${CASE_TMP:?}"/work-* +finish diff --git a/tests/e2e/nf-wiring-2-sv-str-tools.sh b/tests/e2e/nf-wiring-2-sv-str-tools.sh new file mode 100644 index 00000000..e79af992 --- /dev/null +++ b/tests/e2e/nf-wiring-2-sv-str-tools.sh @@ -0,0 +1,132 @@ +#!/usr/bin/env bash +# The Nextflow modules a default run-all.sh run plans for the BAM, with their +# real images: MANTA, DUPHOLD (and DUPHOLD_FILTER), EXPANSION_HUNTER and +# CNVPYTOR (and CNVPYTOR_VCF), through main.nf on the fixture's VCF and BAM. +# Each must publish a readable output with content: Manta's diploid SVs, +# duphold's depth tags on them, ExpansionHunter's repeat calls, and CNVpytor's +# VCF with its index. ExpansionHunter gets the bundled catalog cut to the +# fixture's contigs (as bash-step-09-str does); CNVpytor the pinned v1.3.2 +# GC and mask files docs/00-reference-setup.md downloads. +# Control: a copy of the pipeline whose EXPANSION_HUNTER passes a wrong flag +# fails that task, so a broken module command line cannot ship green. +# nextflow.config ignores a failed EXPANSION_HUNTER (report-only), so the run +# still exits 0: failed_tasks.tsv names it and no repeat VCF is published. +. "$(dirname "$0")/lib.sh" + +command -v nextflow >/dev/null || { fail "nextflow is not on PATH"; finish; } + +G="$GENOME_DIR" +REF="${G}/reference/GRCh38_no_alt_analysis_set.fasta" +VCF="${G}/${SAMPLE}/vcf/${SAMPLE}.vcf.gz" +BAM="${G}/${SAMPLE}/aligned/${SAMPLE}_sorted.bam" +for f in "$REF" "${REF}.fai" "$VCF" "${VCF}.tbi" "$BAM" "${BAM}.bai"; do + check "the base cases left ${f#"$G"/}" test -s "$f" +done + +# ExpansionHunter's bundled GRCh38 catalog, cut to the loci on the fixture's +# contigs: it stops on a locus whose contig the reference lacks. +CAT="${CASE_TMP}/eh_catalog.json" +docker run --rm "$EXPANSIONHUNTER_IMAGE" cat /usr/local/share/ExpansionHunter/variant_catalog/grch38/variant_catalog.json \ + > "${CASE_TMP}/catalog_full.json" +python3 - "${CASE_TMP}/catalog_full.json" "${REF}.fai" "$CAT" <<'PY' +import json, sys +loci = json.load(open(sys.argv[1])) +contigs = {line.split("\t")[0] for line in open(sys.argv[2])} +def regions(locus): + out = [] + for key in ("ReferenceRegion", "OfftargetRegions"): + r = locus.get(key, []) + out += r if isinstance(r, list) else [r] + return out +def known(c): + return c in contigs or "chr" + c in contigs +kept = [l for l in loci if all(known(r.split(":")[0]) for r in regions(l))] +print(f"catalog loci: {len(loci)}, on the fixture's contigs: {len(kept)}") +json.dump(kept, open(sys.argv[3], "w"), indent=1) +PY +check_ge "catalog loci on the fixture's contigs" \ + "$(python3 -c 'import json, sys; print(len(json.load(open(sys.argv[1]))))' "$CAT" 2>/dev/null || echo 0)" 5 + +# CNVpytor's pinned resources (docs/00-reference-setup.md): hg38's GC and mask +# files are read; the others only have to exist. +RES="${CASE_TMP}/cnvpytor" +mkdir -p "$RES" +BASE=https://github.com/abyzovlab/CNVpytor/raw/v1.3.2/cnvpytor/data +for f in gc_hg38.pytor mask_hg38.pytor gc_hg19.pytor mask_hg19.pytor gc_chm13v2.0.pytor gc_chm13v1.1.pytor gc_kn99.pytor; do + curl -fsSL --retry 3 -o "${RES}/${f}" "${BASE}/${f}" || fail "download of ${f}" +done +for f in gc_hg38.pytor mask_hg38.pytor; do + check "${f} is an HDF5 file" sh -c "head -c 8 '${RES}/${f}' | od -An -c | tr -d ' ' | grep -q '^211HDF'" +done + +SHEET="${CASE_TMP}/sheet.csv" +printf 'sample,vcf,vcf_index,bam,bam_index,sex\n%s,%s,%s,%s,%s,male\n' "$SAMPLE" "$VCF" "${VCF}.tbi" "$BAM" "${BAM}.bai" > "$SHEET" +OUT_REL="nf-wiring-sv" +OUT="${G}/${OUT_REL}" +# nf_run DIR NAME TOOLS: main.nf of DIR with the real images on the fixture. The +# fixture BAM is a slice, so indexcov's sex call is noise: --sex_check warn keeps +# the declared male (ExpansionHunter reads it). +nf_run() { + rm -rf "$OUT" + echo "+ nextflow run main.nf -profile docker --tools $3 ($2)" + ( + cd "$CASE_TMP" && nextflow run "$1/main.nf" -profile docker -ansi-log false \ + -work-dir "${E2E_WORK}/nf-work-$2" --input "$SHEET" --reference "$REF" --tools "$3" \ + --expansion_catalog "$CAT" --cnvpytor_resources "$RES" --sex_check warn \ + --outdir "$OUT" --max_cpus 4 --max_memory 14.GB + ) > "${CASE_TMP}/$2.log" 2>&1 +} +nf_run "$REPO" tools manta,duphold,expansion_hunter,cnvpytor +NF_RC=$? +grep -vE 'Pulling|Waiting|Verifying|Download complete|Pull complete|Already exists' "${CASE_TMP}/tools.log" | tail -n 30 +check_eq "nextflow run exits 0" "$NF_RC" 0 +[ "$NF_RC" = 0 ] || grep -E 'Error executing process|Command exit status|Command error' -A 10 "${CASE_TMP}/tools.log" | head -60 +cp "${CASE_TMP}/.nextflow.log" "${E2E_WORK}/logs/${CASE_NAME}.nextflow.log" 2>/dev/null +for p in MANTA DUPHOLD DUPHOLD_FILTER EXPANSION_HUNTER CNVPYTOR CNVPYTOR_VCF; do + check "${p} ran" grep -Eq "Submitted process > [A-Z_:]*:${p} \(${SAMPLE}\)" "${CASE_TMP}/tools.log" +done + +R="${OUT_REL}/${SAMPLE}" +MANTA_VCF="${R}/manta/results/variants/diploidSV.vcf.gz" +check "MANTA writes a readable diploid SV VCF" vcf_ok "$MANTA_VCF" +N_SV=$(vcf_count "$MANTA_VCF") +echo "Manta SVs: ${N_SV}" +check_ge "Manta calls SVs on the fixture" "$N_SV" 1 +DH="${G}/${R}/sv_duphold/${SAMPLE}_sv_duphold.vcf" +check "DUPHOLD declares its depth tag DHFFC" grep -q '^##FORMAT=/dev/null || echo 0)" "$N_SV" +check_ge "and tags them with DHFFC" "$(grep -v '^#' "$DH" 2>/dev/null | grep -c 'DHFFC' || true)" 1 +check "DUPHOLD_FILTER writes a readable VCF" vcf_ok "${R}/sv_filtered/${SAMPLE}_sv_filtered.vcf.gz" +check "and its index" nonempty "${R}/sv_filtered/${SAMPLE}_sv_filtered.vcf.gz.tbi" +EH="${G}/${R}/expansion_hunter/${SAMPLE}_eh.vcf" +check_ge "EXPANSION_HUNTER calls repeats" "$(grep -vc '^#' "$EH" 2>/dev/null || true)" 1 +check_ge "with a repeat count (REPCN) on them" "$(grep -v '^#' "$EH" 2>/dev/null | grep -c 'REPCN' || true)" 1 +check "EXPANSION_HUNTER publishes the catalog it ran" cmp -s "$CAT" "${G}/${R}/expansion_hunter/${SAMPLE}_eh_catalog.json" +CNV="${R}/cnvpytor/${SAMPLE}_cnvs.vcf.gz" +check "CNVPYTOR_VCF writes a readable VCF" vcf_ok "$CNV" +check "and its index" nonempty "${CNV}.tbi" +check "CNVPYTOR writes its call table" test -e "${G}/${R}/cnvpytor/${SAMPLE}_cnvs.txt" +check "the CNV VCF names the fixture's contigs" grep -q 'contig=/dev/null)" +echo "CNVpytor calls: $(vcf_count "$CNV"), table rows: $(wc -l < "${G}/${R}/cnvpytor/${SAMPLE}_cnvs.txt" 2>/dev/null | tr -d ' ')" >> "$E2E_NOTES" +echo "Manta SVs: ${N_SV}, EH loci: $(grep -vc '^#' "$EH" 2>/dev/null || true)" >> "$E2E_NOTES" + +# --- control: a broken EXPANSION_HUNTER command line fails the run --------------- +BROKEN="${CASE_TMP}/broken" +rm -rf "$BROKEN" +mkdir -p "$BROKEN" +cp -R "${REPO}/main.nf" "${REPO}/nextflow.config" "${REPO}/nextflow_schema.json" "${REPO}/versions.env" \ + "${REPO}/conf" "${REPO}/modules" "${REPO}/workflows" "${REPO}/bin" "${REPO}/assets" "$BROKEN/" +sed -i 's/--reads \${bam}/--readz ${bam}/' "${BROKEN}/modules/local/expansion_hunter/main.nf" +check_eq "control copy: EXPANSION_HUNTER passes --readz" \ + "$(grep -c -- '--readz \${bam}' "${BROKEN}/modules/local/expansion_hunter/main.nf")" 1 +rm -f "${CASE_TMP}/failed_tasks.tsv" +nf_run "$BROKEN" broken expansion_hunter +echo "control run exit: $?" +check "control: failed_tasks.tsv names EXPANSION_HUNTER as failed and ignored" \ + grep -qE $'(^|:)EXPANSION_HUNTER\tignored$' "${CASE_TMP}/failed_tasks.tsv" +check "control: no repeat VCF is published" test ! -e "$EH" + +# The tasks' files are root's: remove them from a container. +docker run --rm -v "${E2E_WORK}:/w" "$BCFTOOLS_IMAGE" rm -rf /w/nf-work-tools /w/nf-work-broken "/w/genome/${OUT_REL}" +rm -rf "$BROKEN" "$RES" +finish diff --git a/tests/e2e/short-sample-id-stub.sh b/tests/e2e/short-sample-id-stub.sh index b04a6b53..89f0fdf8 100755 --- a/tests/e2e/short-sample-id-stub.sh +++ b/tests/e2e/short-sample-id-stub.sh @@ -47,7 +47,9 @@ check_eq "S1 with TELOMERE_HUNTER failing: the run exits 0" "$RC" 0 O="${DIR}/out/S1" check "CPSR's HTML is named after the sample" test -e "${O}/cpsr/S1.cpsr.grch38.html" check "CPSR's classification is named after the sample" test -e "${O}/cpsr/S1.cpsr.grch38.classification.tsv.gz" -check "no CPSR file is named after the padded id" test -z "$(find "${O}/cpsr" -name 'S1_cpsr*' 2>/dev/null)" +# CPSR's own files under the padded id would read S1_cpsr.cpsr.grch38.*; the read check's +# S1_cpsr_cdna_check.tsv is named after the sample on purpose and must not match. +check "no CPSR file is named after the padded id" test -z "$(find "${O}/cpsr" -name 'S1_cpsr.*' 2>/dev/null)" CPSR_SH=$(find "${DIR}/work" -name .command.sh -exec grep -l 'cpsr_sample_id' {} + 2>/dev/null | head -n 1) check "the CPSR task ran bin/cpsr_sample_id" test -n "$CPSR_SH" check "the planted failure ran (TelomereHunter's task failed)" has 'TELOMERE_HUNTER.*(Error is ignored|terminated with an error)' "$LOG" diff --git a/tests/test_owngenome_live.py b/tests/test_owngenome_live.py index aebbc859..f7ce603d 100644 --- a/tests/test_owngenome_live.py +++ b/tests/test_owngenome_live.py @@ -109,9 +109,20 @@ def make_root(work): with open(outside, "w") as f: f.write("outside the root\n") os.symlink(outside, os.path.join(d, "clinvar", "escape.tsv")) + # The report links CPSR's folder, which holds its HTML report and its table + os.makedirs(os.path.join(d, "cpsr")) + with open(os.path.join(d, "cpsr", f"{RUN}.cpsr.grch38.html"), "w") as f: + f.write("CPSR\n") + with open(os.path.join(d, "cpsr", f"{RUN}.cpsr.grch38.classification.tsv.gz"), "w") as f: + f.write("not linked\n") + os.makedirs(os.path.join(d, "twopages")) + for n in ("a.html", "b.html"): + with open(os.path.join(d, "twopages", n), "w") as f: + f.write("\n") with open(os.path.join(d, f"{RUN}_report.html"), "w") as f: f.write('' 'hits escape ' + 'the cpsr folder two pages ' 'outside\n') return root @@ -843,6 +854,20 @@ def stored_leaks(text): c4, _, _ = request(port, "GET", f"/runs/{RUN}/..%2F..%2Foutside.txt", headers=auth) check("an unlinked file, a symlink out of the root and a path out of the run are 404", (c1, c2, c3, c4) == (404, 404, 404, 404), (c1, c2, c3, c4)) + code, hd, _ = request(port, "GET", f"/runs/{RUN}/cpsr/", headers=auth) + page = f"/runs/{RUN}/cpsr/{RUN}.cpsr.grch38.html" + check("the report's cpsr/ folder link opens the CPSR HTML report in it", + code == 302 and hd.get("location") == page, (code, hd.get("location"))) + code, hd, body = request(port, "GET", page, headers=auth) + check("the CPSR HTML report in the linked folder is served sandboxed", + code == 200 and b"CPSR" in body + and hd.get("content-security-policy") == owngenome_live.RUN_HTML_CSP, (code, hd)) + c1, _, _ = request(port, "GET", f"/runs/{RUN}/cpsr/{RUN}.cpsr.grch38.classification.tsv.gz", headers=auth) + c2, _, _ = request(port, "GET", f"/runs/{RUN}/twopages/", headers=auth) + c3, _, _ = request(port, "GET", f"/runs/{RUN}/clinvar/", headers=auth) + c4, _, _ = request(port, "GET", f"/runs/{RUN}/twopages/a.html", headers=auth) + check("a linked folder serves only its one HTML; a folder with two HTML files (itself or either file) " + "and an unlinked folder are 404", (c1, c2, c3, c4) == (404, 404, 404, 404), (c1, c2, c3, c4)) code, _, body = request(port, "GET", "/", headers=auth) check("the index lists the run", code == 200 and f'/runs/{RUN}/'.encode() in body, code) log = open(os.path.join(state_dir, "access.log")).read() diff --git a/workflows/annotation.nf b/workflows/annotation.nf index 6acefa59..11c1d100 100644 --- a/workflows/annotation.nf +++ b/workflows/annotation.nf @@ -59,7 +59,8 @@ workflow ANNOTATION { // VEP_NORM first splits and left-aligns the PASS records, as step 13 // does, so an indel finds its gnomAD and ClinVar record. // - ch_vep_vcf = Channel.empty() + ch_vep_vcf = Channel.empty() + ch_vep_vcf_index = Channel.empty() if (params.tools && params.tools.split(',').collect{it.trim()}.contains('vep')) { VEP_NORM(ch_current_vcf, ch_reference, ch_reference_fai) ch_versions = ch_versions.mix(VEP_NORM.out.versions) @@ -73,6 +74,7 @@ workflow ANNOTATION { ) ch_versions = ch_versions.mix(VEP.out.versions) ch_vep_vcf = VEP.out.vcf + ch_vep_vcf_index = VEP.out.vcf_index // Combine VCF + index for downstream ch_current_vcf = VEP.out.vcf @@ -123,6 +125,7 @@ workflow ANNOTATION { // (slivar image). Sequential from VCFANNO (or VEP, or raw VCF) // ch_slivar_vcf = Channel.empty() + ch_slivar_summary = Channel.empty() ch_comphet_pairs = Channel.empty() // [meta, <id>_compound_het_pairs.tsv] if (params.tools && params.tools.split(',').collect{it.trim()}.contains('slivar')) { if (!params.tools.split(',').collect{it.trim()}.contains('vep')) { @@ -144,6 +147,7 @@ workflow ANNOTATION { ) ch_versions = ch_versions.mix(SLIVAR_PRIORITIZE.out.versions, SLIVAR.out.versions, COMPHET_PHASE.out.versions) ch_slivar_vcf = SLIVAR_PRIORITIZE.out.vcf + ch_slivar_summary = SLIVAR_PRIORITIZE.out.summary_tsv ch_comphet_pairs = COMPHET_PHASE.out.pairs } @@ -151,7 +155,8 @@ workflow ANNOTATION { // STEP 3b: CLINICAL_FILTER — Extract clinically relevant variants // Branches from VCFANNO output (parallel with SLIVAR) // - ch_clinical_vcf = Channel.empty() + ch_clinical_vcf = Channel.empty() + ch_clinical_summary = Channel.empty() if (params.tools && params.tools.split(',').collect{it.trim()}.contains('clinical_filter')) { if (!params.tools.split(',').collect{it.trim()}.contains('vep')) { error "clinical_filter requires VEP-annotated input (IMPACT field). Add 'vep' to --tools or remove 'clinical_filter'." @@ -161,13 +166,17 @@ workflow ANNOTATION { ) ch_versions = ch_versions.mix(CLINICAL_FILTER.out.versions) ch_clinical_vcf = CLINICAL_FILTER.out.vcf + ch_clinical_summary = CLINICAL_FILTER.out.summary_tsv } emit: vep_vcf = ch_vep_vcf // channel: [meta, vcf] + vep_vcf_index = ch_vep_vcf_index // channel: [meta, tbi] enriched_vcf = ch_enriched_vcf // channel: [meta, vcf] slivar_vcf = ch_slivar_vcf // channel: [meta, vcf] + slivar_summary = ch_slivar_summary // channel: [meta, <id>_slivar_summary.tsv] comphet_pairs = ch_comphet_pairs // channel: [meta, <id>_compound_het_pairs.tsv] clinical_vcf = ch_clinical_vcf // channel: [meta, vcf] + clinical_summary = ch_clinical_summary // channel: [meta, <id>_clinical_summary.tsv] versions = ch_versions // channel: versions.yml } diff --git a/workflows/clinical.nf b/workflows/clinical.nf index a6f463b1..220ca95c 100644 --- a/workflows/clinical.nf +++ b/workflows/clinical.nf @@ -2,7 +2,8 @@ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ CLINICAL — Clinical Screening & Population Genetics Workflow ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - Runs cancer predisposition (CPSR), runs of homozygosity, polygenic scores + Runs cancer predisposition (CPSR, with step 17's cDNA read check of its + P/LP rows), runs of homozygosity, polygenic scores (pgsc_calc; with an ancestry panel also the sample's projection onto it, step 26), fine-grained ancestry against the panel of setup.sh --ancestry-fine-panel (step 26b), and the mitochondrial haplogroup (from MITO_VARIANTS' Mutect2 @@ -15,6 +16,7 @@ */ include { CPSR } from '../modules/local/cpsr/main' +include { CLINVAR_CDNA_CHECK } from '../modules/local/clinvar_screen/main' include { ROH } from '../modules/local/roh/main' include { PRS_PREPARE } from '../modules/local/prs/main' include { PRS_SCORE_SITES } from '../modules/local/prs/main' @@ -46,12 +48,16 @@ workflow CLINICAL { ch_reference_fai // channel: val(path) — reference .fai ch_mito_vcf // channel: [meta, chrM VCF] — MITO_VARIANTS' Mutect2 calls, for the samples it ran for ch_ancestry_fine_panel // channel: path — the folder setup.sh --ancestry-fine-panel writes, or [] + ch_bam // channel: [meta, bam, bai] — every sample with a BAM (CPSR's cDNA read check) + ch_gencode_exons // channel: path — GENCODE's exons (--gencode_exons), or [] main: ch_versions = Channel.empty() // Initialise output channels with empty defaults ch_cpsr_html = Channel.empty() + ch_cpsr_tsv = Channel.empty() + ch_cpsr_cdna = Channel.empty() ch_roh_regions = Channel.empty() ch_prs_scores = Channel.empty() ch_ancestry_results = Channel.empty() @@ -67,7 +73,23 @@ workflow CLINICAL { if (params.tools && params.tools.split(',').collect{it.trim()}.contains('cpsr')) { CPSR(ch_vcf, ch_pcgr_data, ch_vep_cache_cpsr) ch_cpsr_html = CPSR.out.html_report + ch_cpsr_tsv = CPSR.out.tsv_report ch_versions = ch_versions.mix(CPSR.out.versions) + // Step 17's cDNA read check of CPSR's P/LP rows, with the BAM and + // exons of step 06's check, or [] for a sample without a BAM, so its + // table says "not checked" and why. No construct catalog, as step 17. + CLINVAR_CDNA_CHECK( + CPSR.out.tsv_report + .map { meta, tsv -> [meta.id, meta, tsv] } + .join(ch_bam.map { meta, bam, bai -> [meta.id, [bam, bai]] }, remainder: true) + .filter { row -> row[1] != null } + .map { row -> tuple(row[1], row[2], row[3] ? row[3][0] : [], row[3] ? row[3][1] : []) }, + ch_gencode_exons, + [], + 'cpsr' + ) + ch_cpsr_cdna = CLINVAR_CDNA_CHECK.out.cpsr_cdna_check + ch_versions = ch_versions.mix(CLINVAR_CDNA_CHECK.out.versions) } // @@ -250,6 +272,8 @@ workflow CLINICAL { emit: cpsr_html = ch_cpsr_html + cpsr_tsv = ch_cpsr_tsv // [meta, <id>.cpsr.grch38.classification.tsv.gz] + cpsr_cdna = ch_cpsr_cdna // [meta, <id>_cpsr_cdna_check.tsv] roh_regions = ch_roh_regions prs_scores = ch_prs_scores ancestry_results = ch_ancestry_results diff --git a/workflows/pgx.nf b/workflows/pgx.nf index 70ef5abb..d0abde2d 100644 --- a/workflows/pgx.nf +++ b/workflows/pgx.nf @@ -87,7 +87,8 @@ workflow PGX { .filter { row -> row[1] != null } .map { row -> tuple(row[1], row[2], row[3] ? row[3][0] : [], row[3] ? row[3][1] : []) }, ch_gencode_exons, - Channel.value(file("${projectDir}/assets/construct_variants.tsv", checkIfExists: true)) + Channel.value(file("${projectDir}/assets/construct_variants.tsv", checkIfExists: true)), + 'clinvar' ) ch_clinvar_cdna = CLINVAR_CDNA_CHECK.out.cdna_check ch_clinvar_scan = CLINVAR_CDNA_CHECK.out.cdna_scan diff --git a/workflows/upstream.nf b/workflows/upstream.nf index abfbfb51..578cf2a9 100644 --- a/workflows/upstream.nf +++ b/workflows/upstream.nf @@ -113,6 +113,7 @@ workflow UPSTREAM { emit: bam = ch_bam // [meta, bam, bai]: every BAM, after the sex check indexcov = INDEXCOV.out.indexcov // [meta, indexcov/]: goleft's folder, for the report's sex check + read_count = ALIGN_MINIMAP2.out.read_count // [meta, <id>_read_count.tsv]: the FASTQ rows only, for the report vcf = DEEPVARIANT.out.vcf // [meta, vcf, tbi] gvcf = DEEPVARIANT.out.gvcf // [meta, g.vcf.gz, tbi] versions = ch_versions