Nextflow Execution¶

The pipeline is a Nextflow DSL2 pipeline, main.nf. Each samplesheet row starts from one of three places:

FASTQ: the reads are trimmed (fastp), aligned (minimap2, read group, duplicates marked) and called (DeepVariant, a VCF and a gVCF, chrX and chrY haploid for a male sample);
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 and Contig names and gVCF input. Starting from a provider's VCF: Starting from a Vendor VCF.

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). The Singularity profile is untested (see Profiles).

Quick Start¶

Prerequisites¶

Docker (already required for the bash pipeline)
Java 17 or later (Nextflow 26.10 runtime requirement; CI runs Java 17)
Nextflow, the release CI validates (NEXTFLOW_VERSION in versions.env). Pin it when installing, because the plain installer fetches the newest release:
```
curl -s https://get.nextflow.io | NXF_VER=26.10.5 bash
sudo mv nextflow /usr/local/bin/
```

Run the Pipeline¶

# 1. Create a samplesheet CSV: one row per sample, from FASTQ, a BAM, or a VCF
cat > samplesheet.csv << 'EOF'
sample,fastq_1,fastq_2,bam,bam_index,vcf,vcf_index,sex
sample1,/path/to/sample1_R1.fastq.gz,/path/to/sample1_R2.fastq.gz,,,,,female
sample2,,,/path/to/sample2_sorted.bam,/path/to/sample2_sorted.bam.bai,/path/to/sample2.vcf.gz,/path/to/sample2.vcf.gz.tbi,male
EOF

# 2. Run (default tools need no external databases; prs and vcfanno are
#    skipped with a warning until --pgs_scoring or a score file is set)
nextflow run main.nf \
    --input samplesheet.csv \
    --reference /path/to/GRCh38_no_alt_analysis_set.fasta \
    --outdir ./results \
    -profile docker

# 3. To enable database-requiring tools, add them to --tools with their flags:
#    --tools '...,vep,slivar,clinical_filter'  + --vep_cache /path/to/vep_cache
#    --tools '...,cpsr'                        + --pcgr_data + --vep_cache_cpsr
#    --tools '...,clinvar'                     + --clinvar + --clinvar_index (optional --gencode_exons
#                                                <reference/gencode.v50.basic.exons.bed>: the cDNA read check
#                                                of each hit with the sample's BAM, as step 06 does)
#    --tools '...,expansion_hunter'            + --expansion_catalog (and a sex column)
#    --tools '...,hla_typing'                  + --hla_dat + --hla_genes (hla.dat and the GENCODE gene lines setup.sh installs)
#    --tools '...,annotsv'                     + --annotsv_annotations
#    --tools '...,cnvpytor'                    + --cnvpytor_resources
#    --tools '...,delly'                       (optional --delly_exclude <excl.tsv>, passed as delly sr -x)
#    --tools '...,manta'                       (optional --manta_call_regions <regions.bed.gz>, its .tbi beside it)
#    --tools '...,sample_qc'                   + --somalier_sites + --verifybamid2_panel (setup.sh --sample-qc-data
#                                                installs both; see docs/33-sample-qc.md)
#    --tools '...,cram_archive'                writes a checked CRAM beside each BAM (docs/34-cram-archive.md)
#    --tools '...,cyrius'                      + --cyrius_install (setup.sh --cyrius; non-commercial licence)
#    --tools '...,parascopy'                   + --parascopy_data (setup.sh --parascopy-data; docs/35-paralogs.md)
#    --tools '...,bloodtype'                   no data to install; RhD needs the BAM (docs/41-bloodtype.md)
#    --tools '...,lpa'                         no data to install; needs the BAM (docs/46-lpa.md)
#    --tools '...,chip'                        + --chip_gff (setup.sh --chip-data), optional --chip_germline_resource
#                                                and --chip_pon (each with its .tbi beside it); needs the BAM
#                                                (docs/49-chip.md)
#    --tools '...,archaic'                     + --archaic_data (setup.sh --archaic), --archaic_pop, --hla_genes for the
#                                                genes of each segment (docs/43-archaic.md)
#    --tools '...,ancient'                     + --ancient_panel (setup.sh --ancient-panel); reads the gVCF, so a sample
#                                                without one is skipped with a warning (docs/44-ancient.md)
#    --hpo HP:0001250,HP:0001263               + --exomiser_data (setup.sh --exomiser-data): Exomiser ranks the
#                                                variants against these symptoms; --exomiser_proband names
#                                                the sample when there are several (docs/48-exomiser.md)
#    mito_variants (a default tool) takes --mito_clinvar <clinvar/clinvar_chr.vcf.gz, its .tbi beside it>
#                                                and --gnomad_mito <gnomad_mito/gnomad.genomes.v3.1.sites.chrM.vcf.bgz>
#                                                (setup.sh installs both): ClinVar and gnomAD for each chrM call,
#                                                as step 20 does; without them those columns are '.'
#    --kir true (with hla_typing)              + --kir_dat (setup.sh --kir-data; KIR genes, a second T1K pass)
#    --hla_qc true (with hla_typing)           OptiType as a second HLA caller; HLA-A and HLA-B reach PharmCAT
#                                                only when it agrees with T1K (docs/08-hla-typing.md)
#    With pharmcat, PGX_CONSENSUS gives PharmCAT the HLA types of hla_typing and a CYP2D6 call only when
#    pypgx and cyrius agree (docs/36-pgx-consensus.md); PharmCAT then waits for those BAM steps.
#    telomere_hunter (a default tool) takes --cytoband <cytoBand.hg38.txt> for GRCh38 bands; without it,
#    TelomereHunter uses its own hg19 bands and the run logs a warning
#    An unknown name in --tools stops the run.
#    --pedigree pedigree.tsv                   + --somalier_sites: the family workflow over the samplesheet rows
#                                                it names (two or more, each with a gVCF; docs/multi-sample.md)

Family workflow (opt-in). --pedigree (tab separated: sample, sex, father, mother, partner; assets/pedigree.example.tsv) runs workflows/family.nf over the samplesheet rows it names. --pedigree may be a families file (a leading family column, assets/families.example.tsv). --family_name names the family to run, which a file with more than one family needs, and the folder its results go to (<outdir>/<name>/ instead of <outdir>/family/). The workflow does the same work as step 38: FAMILY_PREPARE checks the pedigree and that each member's gVCF names its own sample, GLNEXUS and GLNEXUS_VCF joint-call the gVCFs (GLNEXUS asks for 16 GB and gives GLnexus half of its memory as --mem-gbytes), FAMILY_RELATE and FAMILY_CHECK stop the run when somalier's relatedness of a pair disagrees with the pedigree (a sample swap), WHATSHAP_PED phases the children of complete trios by inheritance, and FAMILY_REPORT writes each child's allele origin, de novo candidates when both parents are sequenced, and compound heterozygotes among the clinical filter's candidates. With --ancestry_fine_panel (the panel folder of step 26b, setup.sh --ancestry-fine-panel), four more processes find the DNA each pair shares by descent, once FAMILY_CHECK has passed: FAMILY_IBD_SITES keeps the joint call's genotypes at the panel's SNVs on each autosome (a chromosome with fewer than 100 is left out), FAMILY_IBD_PHASE phases them with Beagle against the panel, FAMILY_IBD_SEGMENTS runs hap-ibd (its default 2 cM minimum) with the panel's genetic maps, and FAMILY_IBD_REPORT writes ibd_pairs.tsv, ibd_segments.tsv and ibd_report.md. They feed only the report, so a failure of one is logged and the run goes on. Without the panel none of them runs. Each member needs a gVCF: DEEPVARIANT's for a FASTQ or BAM row, the gvcf column for a VCF row (the run stops on a VCF row without one). Fewer than two named samples, or no --somalier_sites, stops the run at the start. When the pedigree marks two of those samples as each other's partner, COUPLE_CARRIER runs the couple carrier screen of step 39 on them (Couple Carrier Screen): it reads each partner's ClinVar hits, so it needs clinvar in --tools (without it the run logs a warning and skips the screen), uses deletions over an exon only with --gencode_exons (the GENCODE exons BED setup.sh installs as reference/gencode.v50.basic.exons.bed), reads Delly's calls, ExpansionHunter's and Parascopy's SMN1 table when those steps ran, and, when --annotsv_annotations holds AnnotSV's benign loss set (benign_Loss_SV_GRCh38.sorted.bed, with gnomAD-SV's common deletions merged in by setup.sh --annotsv-data), sets aside a Delly deletion that matches a common one there. Without --pedigree none of these processes exists and a run is unchanged. What it checks and what the files hold: Multi-sample comparison.

--reference is the FASTA the reads are aligned to and the BAMs were aligned to, the GRCh38 no-ALT analysis set that setup.sh installs as reference/GRCh38_no_alt_analysis_set.fasta. The Nextflow pipeline does not compare a given BAM's header with it; ./scripts/validate-setup.sh <sample> does, and Realigning after a reference change covers a BAM aligned to another reference.

From reads¶

Parameter	Default	Description
`--skip_trim`	false	Align the raw reads, without fastp (`SKIP_TRIM=true` of step 01b)
`--minimap2_index`	`<reference base>.sr.mmi` beside the FASTA	The minimap2 index built with `-x sr`, the file step 02 builds. When it is not there, the run builds one (1-3 minutes for GRCh38 in the measured runs, cached by `-resume`)
`--intervals`	whole genome	Space-separated regions DeepVariant calls (`INTERVALS` of step 03)
`--sex_check`	`fail`	When the sex indexcov infers from a BAM differs from the samplesheet's: `fail` stops the run before any BAM step starts, `warn` logs both and goes on with the declared sex. With `sample_qc`, the same applies to the sex somalier infers from the reads

The trimmed reads stay in the work directory; the fastp reports are published. ALIGN_MINIMAP2 pipes minimap2 into samtools fixmate, sort and markdup in one task, as step 02 does, so no SAM file is written. It asks for the process_high label's 8 CPUs and 32 GB, and minimap2 and samtools share them: samtools sort takes a quarter of the task's memory (1 GB per thread at 32 GB and 8 CPUs, step 02's -m 1G), and the whole task peaked at 18.3-18.4 GB on GRCh38 in two measured 30-40x runs (runtime). A kill for memory anywhere in the pipe ends the task with that exit status, so it is retried once with twice the memory. DeepVariant also asks for 32 GB by default; --max_memory caps both.

Resume After Failure¶

-resume reuses every task that finished. If a step fails, fix the issue and resume:

nextflow run main.nf -resume [same params as before]

The failed tasks and the ones downstream of them run again. A task that was running when the run stopped starts over: DeepVariant is a single task, so a run stopped during it loses that task's progress. Run the pipeline inside tmux or screen, or with nohup, so closing the terminal does not stop it. Nextflow keys an input file by its path, size and modification time, not by its content: a file rewritten in place with the same size and time is not seen as changed.

One change reruns the alignment once: since alignment and duplicate marking became one task (ALIGN_MINIMAP2, which no longer writes a .sam.gz), a direct nextflow run -resume on the FASTQ samplesheet of a run started with an earlier version reuses fastp (and the minimap2 index, when the run built one) but aligns again, and every task after the alignment runs again too; run-all.sh starts the call after such a run, once it ended, from the BAM and VCF it published and does not align again.

One exception, once: Nextflow 26.04 hashes a map input by its keys and values where 25.10 hashed it differently, so the first -resume of a work directory written by 25.10 reruns every process that takes the meta map, which is every process here. Start a fresh work directory for the first run after the upgrade instead of waiting on a resume that caches nothing.

Samplesheet Format¶

Column	Required	Description
`sample`	Yes	Sample identifier (used as output directory name)
`fastq_1`, `fastq_2`	One of three*	Paired gzipped FASTQ: trimmed, aligned and called
`bam`, `bam_index`	One of three*	Aligned BAM and its index (`.bam.bai`): called when the row has no VCF
`cram`, `crai`	Instead of a BAM	Aligned CRAM and its index, read with `--reference` (the FASTA it was written with). `CRAM_TO_BAM` writes it out as a BAM in the work directory, checked against it, and the row goes on as a BAM row
`vcf`, `vcf_index`	One of three*	Bgzipped VCF (`.vcf.gz`) and its tabix index from any caller; a BAM on the same row is optional
`gvcf`, `gvcf_index`	No	The gVCF of the row's VCF (step 03 writes one) and its index. PharmCAT and PRS read it as they read the gVCF DEEPVARIANT writes; the row is not called again. Needs `vcf` on the row; `run-all.sh` fills it from `vcf/<sample>.g.vcf.gz`
`sex`	On called rows**	`male` or `female`

* A row starts from FASTQ, or from a BAM or CRAM, or from a VCF (with or without a BAM or CRAM); a row with FASTQ and a BAM, CRAM or VCF stops the run. A VCF-only row is valid for annotation and PGx, but most default tools (mosdepth, telomere_hunter, mito_variants) and opt-in tools (expansion_hunter, hla_typing, pypgx, cyrius, parascopy) need a BAM. Provide reads or a BAM for full analysis. A row the pipeline calls gets a VCF and a gVCF, and a VCF row can give its gVCF in the gvcf column; PharmCAT and PRS then read the sites where the sample matches the reference from the gVCF.

The VCF must name its contigs the GRCh38 way with chr (chr1 to chr22, chrX, chrY, chrM); a VCF named 1, MT stops the run with the rename command. A gVCF given in the vcf column stops the run with pharmcat selected: the pipeline expands the reference blocks of the gVCF DeepVariant writes for a called row, or of the one in the gvcf column, not of a gVCF given as the VCF, and a variants-only VCF leaves about half of PharmCAT's genes Unknown. Starting from a Vendor VCF has the commands for both.

** sex is required on every row the pipeline calls (FASTQ, or a BAM or CRAM without a VCF): for a male sample DeepVariant calls chrX and chrY haploid outside the pseudoautosomal regions. It is required on every row with a BAM or CRAM when expansion_hunter is in --tools, where it sets the chrX ploidy (ExpansionHunter's default is female). A row that needs it and lacks it stops the run at parse time.

Sex check. With sample_qc in --tools, somalier also infers the sex from the reads (chrX heterozygosity) and the same rule and --sex_check apply; see Step 33. INDEXCOV (goleft indexcov, seconds per sample: it reads only the .bai) infers each BAM's sex from the chrX and chrY copy numbers and writes it to <sample>_sex_check.tsv. When the row declares a sex and indexcov infers another (or cannot tell), the run stops before any step reads the BAM, and the message gives both values and the copy numbers: the sample is not the one you think, the declared sex is wrong, or the sample has a sex-chromosome aneuploidy. --sex_check warn logs it and goes on with the declared sex. On a small region slice, like the test fixture, indexcov's call is not reliable.

Each sample value must appear once; a repeated id stops the run, because the id names the output directory and keys every per-sample join.

Using Sarek Output¶

If you ran nf-core/sarek for alignment and variant calling, point the samplesheet at sarek's output files. Sarek 3.x writes CRAM by default: give it in the cram and crai columns, with sarek's reference as --reference, or run sarek with --save_output_as_bam to get the .recal.bam files below:

sample,vcf,vcf_index,bam,bam_index
sample1,results/variant_calling/deepvariant/sample1/sample1.deepvariant.vcf.gz,results/variant_calling/deepvariant/sample1/sample1.deepvariant.vcf.gz.tbi,results/preprocessing/recalibrated/sample1/sample1.recal.bam,results/preprocessing/recalibrated/sample1/sample1.recal.bam.bai

Profiles¶

Profile	Description
`docker`	Run with Docker containers (default for local). Every container runs with `--network none`.
`singularity`	Singularity/Apptainer. Untested. One module writes inside its image and needs a writable container: `cnvpytor` (copies resources into its `site-packages`). This profile does not cut the network.
`test`	Minimal test with reduced resources
`test_full`	Full-size test with real WGS data

Combine profiles: -profile docker,test

Resource Configuration¶

Each process asks for a fixed CPU count, a memory and a time limit, from its label or from its own block in conf/base.config. After an exit that usually means a kill for memory (137 and the others conf/base.config lists), the task is retried once with its memory doubled. Time doubles too, except DeepVariant's, which gets --max_time from its first attempt. On Nextflow before 26.10.0, a task that runs past its time limit can stop the whole run instead of being retried (nextflow-io/nextflow#7569, fixed in 26.10.0); we read this in the source and have not reproduced it. These limits cap each request:

Parameter	Default	Description
`--max_cpus`	16	Maximum CPUs per process
`--max_memory`	64.GB	Maximum memory per process
`--max_time`	48.h	Maximum wall time per process

--max_cpus and --max_memory cap each task. They do not limit how much runs at once: Nextflow fills the machine's CPUs and RAM with as many tasks as fit. On Nextflow before 26.10.0 that count of free CPUs and RAM can go wrong and leave a task waiting for room it will never be given; we have seen it in the CI stub runs only (Troubleshooting). 26.10.0, the pinned release, fixes it (nextflow-io/nextflow#7816).

Override for smaller machines, with --max_memory below the RAM the machine reports. On a 32 GB Linux machine the kernel reports a little less than 32 GiB, so a 32 GB cap leaves every 32 GB task (DeepVariant, the minimap2 index and ALIGN_MINIMAP2) larger than the machine, and Nextflow refuses to start it:

nextflow run main.nf --max_cpus 8 --max_memory 30.GB [other params]

Output Structure¶

results/
├── sample1/
│   ├── fastq_trimmed/      # fastp reports, JSON + HTML (FASTQ rows)
│   ├── aligned/            # <sample>_sorted.bam + .bai: minimap2, duplicates marked (FASTQ rows);
│   │                       #   <sample>_sorted.cram + .crai, checked against the BAM (cram_archive)
│   ├── vcf/                # <sample>.vcf.gz and <sample>.g.vcf.gz + .tbi: DeepVariant (called rows)
│   ├── indexcov/           # goleft indexcov coverage plots and .ped (rows with a BAM)
│   ├── <sample>_sex_check.tsv  # the sex indexcov infers, with CNchrX and CNchrY
│   ├── qc/                 # <sample>_sample_qc.tsv: somalier's sex, FREEMIX, same person as (sample_qc);
│   │                       #   somalier/ and verifybamid2/ hold each tool's own files;
│   │                       #   <sample>_read_count.tsv: reads in and primary records in the BAM (ALIGN_MINIMAP2)
│   ├── pharmcat/           # PharmCAT PGx reports (HTML + JSON)
│   ├── clinvar/            # ClinVar pathogenic variant screen: hits as VCF and TSV, and the cDNA read check of each hit
│   ├── pypgx/              # pypgx star allele calling (optional)
│   ├── cpic/               # CPIC drug-gene recommendations (optional)
│   ├── 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
│   ├── roh/                # Runs of homozygosity
│   ├── prs/                # Polygenic risk scores
│   ├── ancestry/           # Ancestry: projection onto the panel (with --ancestry_ref)
│   ├── ancestry_fine/      # Fine-grained ancestry: fractions, painting BED and PNG (opt-in: ancestry_fine)
│   ├── 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)
│   ├── bloodtype/          # Blood group antigens predicted from the genome (opt-in: bloodtype)
│   ├── lpa/                # LPA KIV-2 repeat count from read depth (opt-in: lpa)
│   ├── chip/               # <sample>_chip.tsv, <sample>_chip_depth.tsv and Mutect2's filtered calls: clonal
│   │                       #   haematopoiesis in 13 genes (opt-in: chip)
│   ├── 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
│   ├── expansion_hunter/   # Repeat expansion calls, Stranger's annotation and, when Stranger flags RFC1, its motif count
│   ├── telomere/           # Telomere length estimation
│   ├── coverage/           # Coverage statistics (mosdepth)
│   ├── cyrius/             # CYP2D6 star allele (Cyrius, opt-in) and the CYP2D6 depth check
│   ├── pgx_consensus/      # PharmCAT's outside calls (HLA-A/B, an agreed CYP2D6) and who said what
│   ├── kir/                # KIR genotypes and the IPD-KIR release (--kir)
│   ├── paralogs/           # SMN1/SMN2 copy number (Parascopy, opt-in) and HBA1/HBA2 from depth; Paraphase's copy number and variant
│   │                       #   tables and its own folder (opt-in: paraphase, PacBio HiFi BAMs only)
│   ├── traits/             # <sample>_traits.tsv, <sample>_genotype_rules.tsv and <sample>_traits_genotypes.tsv: the traits panel (traits)
│   ├── consumer_export/    # 23andMe-layout and AncestryDNA-layout raw-data zips and the summary table
│   │                       #   (opt-in: consumer_export)
│   ├── couple/             # <sample>_couple_carrier.tsv: the couple carrier screen for this partner (--pedigree)
│   ├── manta/              # SV calling (optional): diploidSV.vcf.gz with inversions as SVTYPE=INV,
│   │                       #   except a lone inversion junction of 1 Mb or more, kept as its breakend pair;
│   │                       #   diploidSV.raw.vcf.gz as Manta wrote it (an inversion is two BND records)
│   ├── sv_duphold/         # Manta SVs with duphold depth tags (optional)
│   ├── sv_filtered/        # Manta SVs after the duphold depth filter (optional)
│   ├── annotsv/            # AnnotSV ACMG classification of the filtered SVs, and its class 4-5 SVs with their caller, exon-edge and frequency checks (optional)
│   ├── delly/              # SV calling (optional)
│   ├── cnvpytor/           # CNV calling (optional)
│   ├── sv_merged/          # SV consensus of two or more callers, SURVIVOR merge (optional)
│   ├── summary.json        # The numbers the report is rendered from (html_report)
│   ├── <sample>_delta.*    # .txt, .html and .json: what changed since the previous run (run-all.sh only:
│   │                       #   generate-report.sh writes them after the pipeline; delta-report.md)
│   ├── history/            # each run's summary.json and run_manifest.tsv, the next delta's baseline (run-all.sh only)
│   └── *_report.html       # Summary HTML report (published to sample root): QC, ClinVar, PharmCAT, CPIC,
│                           #   CPSR, clinical filter (with its ACMG SF tier), slivar, ROH, mito haplogroup and haplocheck,
│                           #   Y haplogroup and every other step that ran; "Not run" for a tool not selected
├── somalier/               # somalier relate over every sample of the run: samples, pairs, HTML (sample_qc)
├── family/                 # --pedigree (or <name>/ with --family_name): family_plan.tsv, pedigree.tsv,
│                           #   family.vcf.gz (joint call), relationships.tsv, somalier/, family.phased.vcf.gz,
│                           #   <child>_inheritance.tsv.gz, de_novo.tsv, compound_het.tsv, family_report.md;
│                           #   ibd_pairs.tsv, ibd_segments.tsv and ibd_report.md with --ancestry_fine_panel (FAMILY_IBD_REPORT);
│                           #   couple_<a>_<b>.tsv and .md for two partners (COUPLE_CARRIER)
├── multiqc/                # MultiQC report across samples (reads mosdepth: needs a BAM; a VCF-only run logs the skip)
└── pipeline_info/
    ├── timeline_*.html
    ├── report_*.html
    ├── trace_*.txt       # one row per task: the 14 default columns, then process, tag, attempt, cpus,
    │                     #   memory, start, complete (trace.fields in nextflow.config)
    ├── step_history/     # <box>.tsv: bin/run_eta.py's per-task history of this box, for the step times
    │                     #   and the hours left (Where a Run Stands); no sample id, no path
    └── dag_*.svg

Watching a run live¶

Nextflow's weblog, which posts each task's progress to a URL, is the nf-weblog plugin, and Nextflow 26.10 loads it only when a config names it. So nextflow.config has no weblog block: with OWNGENOME_LIVE_URL set, run-all.sh writes <sample>/nextflow/live.config, which loads the pinned plugin and points the weblog at OWNGENOME_LIVE_URL/ingest/<sample>, and passes it with -c. A run without OWNGENOME_LIVE_URL loads no plugin. Nextflow downloads the plugin on the first live run and keeps it in ~/.nextflow/plugins. For a bare nextflow run, write the same file and pass -c; Watch a Run Live has it and the live site that reads the messages.

The pipeline does not anonymise anything: the outputs carry whatever identified you in the input. Before you share a file, know what it holds:

The VCF's sample name (the last column of its #CHROM line) is repeated in the ClinVar VCFs (clinvar/<sample>_clinvar_hits.vcf, <sample>_pass.vcf.gz), on every line of roh/<sample>_roh.txt, in mito/<sample>_haplogroup.txt and as sampleId in the PharmCAT JSON files.
The input's header lines pass through into the ClinVar VCFs, including the provider's and bcftools' command lines, which often name the sample or a file.
A BAM the pipeline aligned names the sample label in its @RG line, and its @PG lines hold the command lines with your file paths.
The input file name is in the ##bcftools_viewCommand and ##bcftools_normCommand lines of the ClinVar VCFs, and can be in the command lines other tools print into their outputs.
pipeline_info/ (report, timeline, trace) holds absolute paths of your machine, and the trace names each sample in its name and tag columns. step_history/ holds neither, but it does say what ran on your box and when.
The samplesheet's sample label names every output folder and file.

A neutral input removes the metadata that names you. Use a label that does not name you, name the file after it, and keep only the header lines the tools read. The recipe below lists the lines to keep and renames the sample column. A list of keys to delete would miss lines, because callers name them differently.

It does not anonymise the data. Genotype files (VCF, gVCF, BAM, CRAM) can identify you and your relatives through genealogy databases: Erlich et al. 2018 (Science 362:690) projected a third-cousin-or-closer match for about 60% of searches for people of European descent, and Gymrek et al. 2013 (Science 339:321) recovered surnames from Y-chromosome STRs.

bcftools view --no-version -h in.vcf.gz | grep -E '^(##(fileformat|FILTER|INFO|FORMAT|ALT|contig)=|#CHROM)' > h.txt
echo SAMPLE > names.txt
bcftools reheader -h h.txt -s names.txt -o SAMPLE.vcf.gz in.vcf.gz
bcftools index -t SAMPLE.vcf.gz

bcftools annotate -x cannot remove these lines (it exits with "No matching tag"). The same recipe, with bcftools from the pinned image, is step 4 of Starting from a Vendor VCF. Do not share pipeline_info/.

The pipeline and the single-step scripts¶

Run the pipeline for a sample: it runs every step a default run selects, in parallel where the inputs allow, and -resume reruns only what changed. Run a numbered script to run or rerun one step by hand, on the layout run-all.sh and setup.sh use; the scripts and the modules run the same commands at the same tool releases (see the parity section below). The pipeline needs Java 17 or later and Nextflow beside Docker; the scripts need Docker only. The HPC/Singularity profile exists but is untested.

Known Limitations & Design Decisions¶

Scope: from reads to report¶

The pipeline starts from paired short-read FASTQ, from a BAM or from a VCF. The reads go through fastp, minimap2 with the sr preset (samtools fixmate, sort, markdup) and DeepVariant's WGS model, with the flags of steps 01b, 02 and 03. Other starting points are scripts (below): ORA files, long reads, other aligners and callers, chip data.

Bash vs Nextflow parity¶

The pipeline and the scripts run the same commands at the same tool releases, so on the same input they must give the same answers. The E2E workflow checks this on the test fixture: the scripts 01b, 02, 03, 06, 07, 11 and 25, and the pipeline from the same FASTQ pair, then scripts/ci/parity-diff.sh compares, item by item:

Item	Compared
alignment	the two BAMs record by record: the md5 of every record (`samtools view`), the md5 of the header without its `@PG` lines (each side's command lines), `samtools flagstat` and `samtools idxstats` (read through the index)
variants	`bcftools isec` of the two VCFs: records in one only, and shared records with another FILTER or GT
gvcf	every gVCF record
clinvar	the rows of `<sample>_clinvar_hits.tsv`
roh	the `bcftools roh` segments and the 5 Mb summary
pharmcat	the diplotype of each gene
prs	each score's sum, matched variant count and input (gVCF or VCF)
sv_merged	each SV consensus record (position, end, type, length, support and each caller's GT); on the Nextflow hardening run, against steps 22 and 37 on the same calls and BAM
y_haplogroup	Yleaf's prediction table; on the same hardening run
traits	`<sample>_traits.tsv` and `<sample>_traits_genotypes.tsv`, byte for byte; in case `traits-panel-1`, which runs step 42 and `--tools traits,html_report` on the same VCF and gVCF. `<sample>_genotype_rules.tsv` with WhatsHap's phase from a BAM, byte for byte, in case `genotype-rules-1-phase`
archaic	`<sample>_archaic.tsv`, `<sample>_archaic_segments.tsv` and `<sample>_archaic.txt`, byte for byte; in case `archaic-ancestry-1`, which runs step 43 and `--tools archaic` on the same VCF and table
ancient	the ranked groups of `<sample>_ancient.tsv` (`rank`, `group`, `n`) and `sites_used`; in case `ancient-similarity-1`, which runs step 44 and `--tools ancient,html_report` on the same gVCF and panel

A difference fails the check unless the script lists it with its reason; none is listed today. Output file names and folders differ (the scripts write under $GENOME_DIR/<sample>/, the pipeline under <outdir>/<sample>/, mapped in the script).

Where a module and its script differ on purpose:

Step	Script	Module
Alignment (02)	pipes minimap2 into samtools, each in its own image (`MINIMAP2_IMAGE`, `SAMTOOLS_IMAGE`)	`ALIGN_MINIMAP2` runs the same pipe in one task, so in one image: `MINIMAP2_SAMTOOLS_IMAGE`, a BioContainers image with both tools at the releases of those two. `scripts/ci/check-images.sh` fails when the releases part, and the parity check above compares the BAMs record by record. The task writes no SAM file
Sex check (16)	runs beside the other steps and stops itself on a mismatch	`INDEXCOV` runs before every BAM step, and a mismatch stops the run before DeepVariant starts
DeepVariant (03)	`MODEL_TYPE` picks WGS, WES, PACBIO or ONT_R104	the WGS model only: the pipeline takes paired short reads
HLA typing (08)	keeps the T1K index under `t1k_idx/`, named after the T1K version, the IPD-IMGT/HLA release and the GENCODE release	`T1K_BUILD` builds it once per run for every sample; the task hash covers the same three, and `-resume` reuses it
PRS (25)	scores the list in `assets/pgs_scores.tsv`, downloading a missing file	scores every file `--pgs_scoring` holds, labelled from `assets/pgs_scores.tsv` (an id not in it is labelled with its file's `trait_reported`); `PRS` runs pgsc_calc on the host, see No network inside the containers
ExpansionHunter (09)	uses STRchive's catalog when `setup.sh` built it, else the GRCh38 catalog inside the image, or `EH_CATALOG`	needs `--expansion_catalog` (`run-all.sh` passes the same one)
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=<outdir> scripts/24-html-report.sh <sample>` 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)

Scripts with no module, and why:

Script	Why it stays bash-only
`01-ora-to-fastq.sh`	`orad` is a native binary that runs on the host, on the raw reads
`02a-alignment-bwamem2.sh`, `03a-gatk-haplotypecaller.sh`, `03b-freebayes.sh`, `03c-strelka2-germline.sh`, `03d-octopus.sh`, `03e-clair3.sh`	alternative and legacy aligners and callers, kept for benchmarking: the pipeline aligns with minimap2 and calls with DeepVariant, and takes the VCF of any other caller in the `vcf` column
`02b-alignment-longread.sh`, `04c-sniffles2.sh`	the long-read path: the samplesheet holds one short-read BAM per sample
`04a-tiddit.sh`	an alternative to Manta; its local assembly needs the classic BWA index
`04b-gridss.sh`	needs a 31 GB Java heap and the BWA index, more than the SV modules are sized for
`14-imputation-prep.sh`	writes per-chromosome files for an upload to an imputation server, a manual step
`29-mutect2-somatic.sh`	experimental tumor-only somatic calling
`chip-to-vcf.sh`	turns consumer array data into a VCF: an input for the pipeline, not a step on one
`benchmark-variants.sh`	compares callers with each other or with a truth set: a check of the calling, not an analysis of the sample
`generate-report.sh`	the text report. Run `GENOME_DIR=<outdir> scripts/generate-report.sh <sample>` on the Nextflow output
`run-all.sh`, `setup.sh`, `validate-setup.sh`	running, installing and checking the bash path

Reference databases not auto-downloaded¶

Several tools require large reference databases that are not automatically downloaded by the pipeline. You must obtain and provide paths for these yourself:

Parameter	Required by	Size
`--vep_cache`	VEP annotation	~28 GB download (release 116)
`--pcgr_data`	CPSR cancer predisposition	~7 GB bundle (20260620)
`--vep_cache_cpsr`	CPSR	~24 GB download (release 115)
`--pypgx_bundle`	PyPGx star allele calling	~370 MB
`--annotsv_annotations`	AnnotSV SV classification; with `--pedigree`, also the couple screen's common-deletion filter (its benign loss set)	~5.3 GB download
`--cadd_snv`, `--spliceai_snv`, etc.	vcfanno score annotation	~100 GB total
`--gnomad_constraint`	Slivar gene constraint	~95 MB
`--pgs_scoring`	Polygenic risk scores	~400 MB (`setup.sh` fills `<genome_dir>/prs_scores`)
`--ancestry_ref`	PRS percentiles and ancestry (step 26)	~7 GB (`setup.sh --ancestry-panel`)
`--ancestry_fine_panel`	Fine-grained ancestry (step 26b, `ancestry_fine`); with `--pedigree`, also the family's shared segments (`FAMILY_IBD_*`)	the folder `setup.sh --ancestry-fine-panel` writes

Tools that require external databases (VEP, slivar, clinvar, CPSR, ExpansionHunter, HLA typing, pypgx, AnnotSV, CNVpytor) will fail at startup if enabled in --tools without their required parameters. vcfanno and prs are in the default tools but are skipped, with a warning, until a score file or --pgs_scoring is set. The gnomAD constraint table is optional for slivar; when it is set and no gene matches it, the task fails.

Ancestry reference panel¶

--ancestry_ref is pgsc_calc's panel, <genome_dir>/reference/pgsc_calc/pgsc_1000G_v1.tar.zst, with the pgsc_1000G_v1_GRCh38_sites.tsv that setup.sh --ancestry-panel writes beside it (the run stops when the list is missing). With it, prs genotypes the panel's SNVs from the gVCF too, pgsc_calc projects the sample onto the panel, each score gets a percentile among the most similar population, and PRS_SUMMARY publishes step 26's table in ancestry/. ancestry in --tools needs prs and --ancestry_ref; without them it is skipped with a warning. See step 26.

SV consensus merge¶

survivor_merge runs SURVIVOR merge (step 22's parameters: breakpoints within 1 kb, same type and strands, at least 50 bp, two or more callers) over the callers this run selected among manta, delly and cnvpytor, each cut to its PASS records first (SURVIVOR_PREP, SURVIVOR_MERGE, SURVIVOR_SORT). CNVpytor's depth-only calls count as one caller like the others; their coarse breakpoints often lie more than 1 kb from the paired-end callers'. TIDDIT and Sniffles2 have no module, so the consensus of the script (step 22) can hold those two callers as well. GRIDSS is in neither: its breakend (BND) records never match the DEL, DUP and INV records of the other callers. See step 22.

FILTER=PASS required¶

ClinVar screen, clinical filter and slivar keep only records with FILTER=PASS. Before any analysis, VCF_PRECHECK counts the FILTER values of each sample. A VCF with no PASS record at all (for example unfiltered GATK HaplotypeCaller or FreeBayes output, where FILTER is .) stops the run with a message naming the sample, because every PASS-only step would report zero hits. Filter it with your caller's recommended filters, or add --allow_unfiltered to treat FILTER . as PASS for that file. A VCF with any PASS record is used as given.

One sample per VCF¶

One samplesheet row is one sample, so a VCF given on a row must hold one sample column. VCF_PRECHECK stops the run, before any analysis, on a VCF with more than one (a joint-called family file, for example): the message names the sample count and the first five sample names, and prints the bcftools view -s <name> -a -c 1 command that keeps one sample's column. That command writes a variant-only VCF, which is what the vcf column takes. A gvcf for the same row is split with -s <name> alone, so its reference blocks stay. validate-setup.sh <sample> applies the same rule to vcf/<sample>.vcf.gz, and fails when it cannot read that file's header.

Contig names and gVCF input¶

VCF_PRECHECK also stops the run, before any analysis, in four cases, and the message names the sample and the fix:

No contig that holds records is named the chr way (1, MT instead of chr1, chrM). Without this stop the mito haplogroup file comes out empty and chrX segments leak into the autosomal ROH summary, with exit 0. The message prints the bcftools annotate --rename-chrs command with its 25-line map.
pharmcat is selected and the VCF is a gVCF (a ##GVCFBlock header line, or reference-block records: ALT <*>, <NON_REF> or . with INFO/END), or only its name says so (.g.vcf, .genomic.vcf). PharmCAT refuses both. Without pharmcat, a gVCF runs: ROH, the mito haplogroup and the ClinVar screen give the same results as on the matching variants-only file.
A ##contig line for chr1-22, chrX or chrY (or 1-22, X, Y) gives a length GRCh38 does not have, for example chr1 length 249250621 (GRCh38: 248956422), which is GRCh37 (hg19). The message lists up to five such contigs. Renaming the contigs does not fix this: get a GRCh38 VCF, or start from the FASTQ.
The header gives no length for chr1-22, chrX or chrY (no ##contig lines, or lines without lengths), and the REF bases of the file's first 1000 records (less any on a contig the reference lacks) differ from --reference in more than 5% of them. bcftools norm -c x does the comparison, as in validate-setup.sh; the message gives the count and the rate, for example they differ in 45 of the first 60 records (75.0%). On GRCh38 none or almost none differ, on another build most do; 5% is a line drawn, not measured. When no record is left to compare, or bcftools norm fails, the task log says so and the run goes on.

Starting from a Vendor VCF has the commands that fix the first two.

Security model¶

This pipeline is designed for personal, single-user use on trusted data. Sample labels are restricted to [A-Za-z0-9._-] (they name folders and go into shell commands); this is not anonymisation, see Before you share outputs. HTML report fields from VCF INFO are escaped to prevent XSS. However, it is not hardened for multi-tenant or untrusted-input scenarios. Do not expose the pipeline or its outputs as a web service without additional security review.

Failed commands fail the task¶

Every task script runs under bash -euo pipefail (process.shell in conf/base.config). In a pipe such as bcftools view -f PASS in.vcf.gz | bcftools +split-vep ..., a failure of the first command, a truncated input for example, fails the task. Under Nextflow's default bash -ue only the last command counted, and the task wrote a short, valid file and passed.

No network inside the containers¶

With -profile docker every container runs with --network none (process.containerOptions in nextflow.config). The steps read only their inputs, so a tool that tries to download something at run time fails instead of fetching an unpinned file. The singularity profile does not cut the network.

One process runs on the host instead of in a container: PRS starts pgsc_calc, a Nextflow pipeline of its own that starts its own containers. It gets the images of versions.env through conf/containers.config (ext.pipeline_images, written by scripts/ci/gen-containers-config.sh), and gives every one of its containers -u with your ids and --network none, followed by the run's own docker.runOptions and the PRS process's containerOptions less their user and network flags (step 25). With --pgsc_calc <genome_dir>/tools/pgsc_calc-<release> (setup.sh installs it with its nf-schema plugin) it runs offline; without, the task first fetches GitHub's archive of the pinned release, checked against PGSC_CALC_SHA256, and Nextflow the plugin.

Cyrius, opt-in¶

Cyrius (CYP2D6 star alleles) is under a non-commercial licence and is not in the default tools. scripts/setup.sh --cyrius <genome_dir> installs it once into tools/cyrius-1.1.1/, with pip --require-hashes --no-deps --only-binary :all: from scripts/cyrius-constraints.txt, which pins every wheel by its sha256. The CYRIUS module runs it from there (--tools ...,cyrius --cyrius_install <genome_dir>/tools/cyrius-1.1.1) with no network, like every other task. See step 21.

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-<line>. 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). 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.

Relationship to nf-core¶

This pipeline uses nf-core template patterns and tooling for code quality, but is not an official nf-core pipeline (it uses a GPL-3.0 license; nf-core requires MIT).

Individual modules (PharmCAT, pypgx, slivar) will be contributed to nf-core/modules under MIT license for use by the broader community.

Acknowledgement¶

This pipeline was created using tools and best practices from the nf-core community (Ewels et al., 2020, Nat Biotechnol). nf-core components used here are released under the MIT license.