Pipeline Steps¶

Each step is an independent module under scripts/. This page covers what each one does, what it writes, and what makes it skip itself.

Step 0: Verify Environment¶

scripts/0_verify_environment.nf

The gate for everything else. Nine stages run in parallel, each writing its own section, and the results are assembled into Output/Reports/0_verify_environment.txt. Nothing downstream starts until every one of them has passed.

The nine stages¶

Stage	Reports as	Fails when
`CheckReference`	`REFERENCE FILE CHECK`	The reference is missing. Gzipped or plain are both accepted
`CheckGFF`	`GFF FILE CHECK`	The annotation is missing while `annotate = true`
`CheckData`	`DATA SOURCE CHECK`, `DATA FOLDER CHECK`, `DATA FILES CHECK`	`dataSource` does not resolve, the directory is absent, or nothing matches `readPattern`
`CheckMetadataFile`	`METADATA FILE CHECK`, `METADATA CHECK`, `METADATA SAMPLE MATCH`, `METADATA CHANGE CHECK`	`metadata.csv` is missing or malformed, a FASTQ sample has no row, or the analysis-affecting columns differ from what produced your results
`CheckInstalledSoftware`	`SOFTWARE CHECK`	A tool named in `params.software` is not on `PATH`
`CheckTrimParameters`	`TRIM PARAMETER CHECK`	`autodetect = false` with `adapter1` or `adapter2` unset
`CheckRunParameters`	`PIPELINE VERSION`, `RUN PARAMETER CHECK`	The release differs from the one that produced your results, or an analysis-affecting parameter has changed
`CheckDirectories`	`DIRECTORY CHECK`	`mainDir` and `storageDir` are the same path, or either is the installation itself
`CheckMultiRun`	`MULTI-RUN CHECK`	The run table is missing, unparseable, or names a column that is not a parameter

CheckMetadataFile reports more than it refuses. It prints the pooling it worked out — which rows merged into which VCF column — and the pool sizes it will apply, before any compute is spent. A mistake there is visible in seconds rather than in a result months later.

The consistency guards¶

Three of those stages are not validating your input; they are checking that what you are about to run matches what produced the results already on disk. They exist because of how resume works: completed steps are skipped by looking for output files, not by checking what produced them, so a changed poolSize would leave one Frequencies/ folder holding tables computed under two different thresholds, invisibly.

Record	Contents
`.poolseqflow_version`	The release that produced these results. A mismatch is a hard stop — nothing else is compared, because the parameter set itself moves between releases
`.poolseqflow_params`	The analysis-affecting parameters, mirrored to a readable `run_parameters.txt`
`.parameters.config`	Your configuration file, copied verbatim
`.multirun.csv`	Your run table, copied verbatim, if you used one
`.poolseqflow_metadata`	The `RG_` and `param_` columns and the row order, kept beside the results they describe

Path and resource parameters are excluded, along with the software entries — they change where and how fast work happens rather than what the answer is.

When a guard trips, the run stops before any work happens and the report names what to delete. Deleting it is what clears the check. How much has to go depends on what changed: a reordered metadata.csv invalidates less than an edited tag value, and a changed pool size less again (table).

Step 1: Build Reference Dictionaries¶

scripts/1_build_dictionaries.nf

Puts an uncompressed reference into Reference/Dictionaries/ — decompressing it if it was gzipped, copying it if it was not — and builds three index sets there, in parallel:

Sub-step	Produces
`CreateBwaIndex`	`Ref.fasta.{amb,ann,bwt,pac,sa}`
`CreateSamtoolsFaiIndex`	`Ref.fasta.fai`
`BuildSnpEffDb`	`Reference/Dictionaries/snpEff/` and `snpEff.config` (only if `annotate = true`)

All of it lives under mainDir, beside the reference it was built from, and none of it is copied to storageDir. Dictionaries are working material: they are derived from your reference and can be rebuilt from it, so they are not a result to keep. Under a run table, runs sharing a reference build them once.

The SnpEff database name is derived from the GFF filename — reference.gff.gz becomes database reference.gff. The build copies the reference and GFF into a SnpEff data/ layout, generates a minimal config, and verifies that .bin files were produced before declaring success. Completion is marked by .build_complete, which is what the resume check looks for.

Build options are -gff3 -noCheckCds -noCheckProtein -v. The two -noCheck flags suppress SnpEff's protein and CDS consistency checks, which fail on many non-model GFFs for reasons that do not affect variant annotation.

Step 2: Trim & QC¶

scripts/2_trim_reads.nf

Two sub-steps.

TrimReads runs Trim Galore with --fastqc --paired --retain_unpaired -q 25, which removes adapters and low-quality 3′ ends and produces a FastQC report on the result.

ClipReads parses that FastQC report and derives per-sample clip points from the per-base composition table, then applies them with cutadapt and re-runs FastQC on the output.

The clipping algorithm, its failure modes and the one parameter worth tuning are covered in Trimming & Clipping.

Trimmed reads are deleted once clipping has consumed them. ClipReads is the only step configured to retry on failure.

Step 3: Align¶

scripts/3_align.nf

bwa mem -K 10000000 -T 30 -t <threads> reference R1_clipped R2_clipped \
  | samtools view -b -o <sample>.bam

Option	Parameter	Effect
`-T 30`	`bwa.minScoreOutput`	Minimum alignment score to report
`-K 10000000`	`bwa.batchSize`	Fixed bases per batch — makes output deterministic regardless of thread count

-K is worth knowing about. Without it, BWA processes a batch sized by thread count, so the same input aligned with different threads can produce slightly different output. Fixing the batch size makes a run reproducible across machines.

bwa.options is the whole flag string, and it ships commented out because the pipeline builds it for you from the two values above. There is nothing to set: change minScoreOutput or batchSize and the string follows.

// options      = "-K ${params.bwa.batchSize} -T ${params.bwa.minScoreOutput}"

Uncomment it only when you want flags those two cannot express — a different -A/-B scoring, say. From then on the string is used exactly as you write it and the two values above stop being read, so -K is gone unless you put it back yourself, and the run quietly loses the determinism the section above is about. -t is the exception: it is not in the string at either end, so the thread count still comes from the cores ladder whatever you pin here.

Like the values it replaces, it is analysis-affecting — step 0 refuses a run whose string differs from the one recorded beside your existing alignments.

Output goes to Output/Aligned/ as unsorted BAM.

Step 4: Clean BAM Files¶

scripts/4_clean.nf

A single streamed pipeline of seven operations: name-sort, fixmate, coordinate-sort, mark and remove duplicates, add read groups, filter, index. The @RG string is assembled per sample from that sample's row in metadata.csv, skipping empty fields.

Full detail, including why duplicates are removed rather than marked and what the MAPQ floor costs you: Alignment & Cleaning.

Output: Output/Ready/<sample>_ready.bam and .bai.

Step 5: Generate Reports¶

scripts/5_reports.nf

Three independent sub-steps per sample:

Sub-step	Command	Output
`AlignmentReport`	`bamtools stats`	`Output/Reports/Alignment/<sample>_alignment_report.txt`
`CoverageReport`	`samtools coverage`	`Output/Reports/Coverage/<sample>_coverage_report.txt`
`DepthProfile`	`samtools stats`	`Output/Reports/Depth/<sample>_depth_histogram.tsv` and `_depth_report.txt`

DepthProfile also decides the depth ceiling step 6 applies to that sample, which is why this step is no longer a leaf: it passes the BAM, its index and the chosen ceiling on to variant calling. It runs whatever capBAM.maxDepth is set to — the histogram is worth having even when nothing is capped, and the reports are what step 6 is judged by.

Reading the depth report on every run is the habit worth forming. The capped BAM step 6 builds is transient, so this file is the only record of which reads were actually called for that sample.

Step 6: Variant Calling¶

scripts/6_variant_call.nf

bcftools mpileup -B -C 50 -q 30 -Q 30 -d 0 -a AD,DP,SP,INFO/AD -Ou -f reference <all BAMs> \
  | bcftools call -m -A -v -Ov -o <name>.vcf

Before that, CapBAM truncates each BAM to the ceiling step 5 chose for it. It runs only for samples that have one: a sample the detector left alone, or a run with capBAM.maxDepth = 0, goes to bcftools with its ready BAM unchanged, and no CapBAM task appears in the trace. The capped BAM never reaches either storage root — it is built here, read by the pileup, and discarded, which is why step 5's reports are the durable record.

One joint task over every BAM, producing a multi-sample VCF with AD and DP FORMAT fields. Every flag and its consequence: Variant Calling.

The BAMs are sorted before being handed to bcftools, in metadata.csv row order. bcftools orders VCF columns by command-line order, so without this the column order would follow task-completion order — three consecutive runs on identical input gave three different orders.

A header fix is applied afterwards: INFO/MQ is declared Integer by bcftools but can carry a float, so the declaration is rewritten to Float. Without it, strict VCF parsers reject the file.

Step 7: VCF → Allele Frequency Tables¶

scripts/7_vcf2freq.nf

Five sub-steps in a serial chain, each deleting its input once its output is safe:

#	Sub-step	Does
1	`SortRefAltByFrequency`	Re-encodes so the most-read allele is `REF`; recomputes `DP` from `AD`; sets `GT` to `./.`
2	`FilterPotentialFalsePositives`	Splits multiallelics, applies the cross-sample support test, rejoins, re-normalizes
3	`DepthAndQualityFilter`	`bcftools view -e "FMT/DP<20"` then `vcftools --minQ 30`
4	`SplitSNPsAndINDELs`	Two vcftools passes into a SNP VCF and an INDEL VCF
5	`CalculateFrequencies`	Extracts `AD`, publishes it as the depth table, then divides to frequencies; runs once per split file

Sub-step 2 is the pool-aware core of the pipeline and is documented in full in The Filter Chain. The minimum credible frequency it enforces is

\[f_{\min} = \frac{1}{2 \times \text{ploidy} \times \text{poolSize}}\]

but the test is not a plain cutoff in either direction. An allele must reach that frequency in a fraction of samples, which is what lets the threshold sit so low — and each sample is judged against its own pool's threshold, taken from param_poolSize, rather than against one number for the whole run.

Output: Output/Frequencies/<name>_snp_freq.tsv and <name>_indel_freq.tsv, each beside a _depth.tsv holding the read counts it was computed from. Formats: Interpreting Results and The depth tables.

Step 8: Annotate Variants¶

scripts/8_annotate_variants.nf — optional, controlled by annotate.

bcftools norm -m - <name>.vcf | snpEff -v -stats snpeff_summary.html <db>

Splits multiallelic sites onto separate lines — SnpEff annotates one alternate allele per record — and annotates against the database built in step 1.

This runs on the raw call set

Step 8 takes step 6's output, not step 7's. It runs in parallel with the frequency branch, so <name>_annotated.vcf contains sites the step 7 filters removed, encoded against the original reference rather than the major allele. To attach annotations to your frequency tables, join on CHROM/POS and expect unmatched rows on the annotation side.

Output: Output/VCF/<name>_annotated.vcf, Output/Reports/snpeff_summary.html and Output/Reports/snpeff_summary.genes.txt. SnpEff names the gene table after the summary and writes the two together, so they arrive and are replaced as a pair.

Setting annotate = false skips this step and makes gffFile unnecessary — step 1 also stops building the SnpEff database.

The parameters that control it — annotate, gffFile and the snpEff block — are in Annotations.

Promotion¶

scripts/9_completion.nf

Not an analysis step, and it produces nothing of its own. It moves each finished artifact from the working volume to permanent storage once nothing needs it any more, so every byte crosses between the two exactly once.

An artifact enters mainDir/Utilized/ precisely when a later step will read it again. Things with no consumer at all — the unpaired reads, the trimming reports, the FastQC output, both step 5 reports — go straight to storageDir and never appear there.

What moves it is a completion signal from the last step that reads it, not the artifact itself. Which step that is can depend on your configuration: <name>.vcf is read by step 7 always and by step 8 only when annotate = true, with no ordering between them, so the gate has to wait for whichever set applies to this run.

If a run is interrupted, artifacts stay in Utilized/ and the skip checks find them there — a resumed run counts them as done and picks up where it stopped, rather than repeating the work because the file is not in its final place yet.