oxo-flow-atacseq — ATAC-seq: peak calling and QC

★ Verified · ⇄ Official port of nf-core/atacseq @ 2.1.2 — same tools, same versions, same commands. Part of the oxo-flow-community catalog.

ATAC-seq (Assay for Transposase-Accessible Chromatin) analysis for a cohort of samples: raw-read QC (FastQC), adapter trimming (Trim Galore), BWA-MEM alignment, Picard mark-duplicates, BAMTools filtering, MACS2 broad-peak calling, HOMER peak annotation, FRiP scoring, normalised bigWig tracks, deepTools QC plots (profile, heatmap and fingerprint) and a single MultiQC report. Input is a directory of <sample>.fastq.gz reads plus pre-built reference files; output lands in results/ with per-step subdirectories and a combined MultiQC HTML report.

The default plan is the upstream single-end aligner=bwa main path (15 rules, 29 rule instances for the fixture cohort). Upstream branches that are off the main path — paired-end, alternative aligners (Bowtie2/Chromap/STAR), reference preparation, mitochondrial filtering, consensus peaks / DESeq2, preseq, Picard metrics, ataqv, IGV, R QC plots — are ported as when-gated branches (see Gated branches); each activates only when its config key is set, so the default plan is unchanged.

Requires oxo-flow >= 0.12.0. Recommended — prebuilt release binary:

curl -fL -o oxo-flow.tar.gz \
  https://github.com/Traitome/oxo-flow/releases/latest/download/oxo-flow-latest-x86_64-unknown-linux-gnu.tar.gz
tar xzf oxo-flow.tar.gz
sudo mv oxo-flow /usr/local/bin/

Alternative via conda:

conda install -c bioconda oxo-flow-cli

Note the bioconda package may lag the latest release; other platform binaries are on the releases page.

git clone https://github.com/oxo-flow-community/oxo-flow-atacseq.git
cd oxo-flow-atacseq

Reference data (pre-built files, declared under [config] in main.oxoflow):

The port does not auto-do what upstream prepare_genome does on the default config: reference/annotation files must be pre-built uncompressed GTF/BED/FASTA and unpacked index files. The convenience branches are available when you need them — set gtf to a .gz source? Leave it canonical and point gtf_src at the archive. The khmer genome-size estimate is on independently: set macs_gsize = "" to auto-estimate from the FASTA at read_length (default 2.7e9 = GRCh37/38 @ 50 bp — see Fidelity).

To do upstream prepare_genome:	set these config keys	happens
decompress a .gz GTF	gtf_src (canonical gtf stays fixed)	ref::gunzip_gtf (upstream GUNZIP_GTF)
GFF3 → GTF	gff with gtf left empty	ref::gffread (upstream GFFREAD); GTF goes to results/genome/index/ref.gtf
decompress a .gz gene/TSS BED	gene_bed_src / tss_bed_src	ref::gunzip_gene_bed / ref::gunzip_tss_bed
unpack a BWA index tarball	bwa_index_src (a .tar.gz)	ref::untar_bwa_index; the extracted index is found at runtime (find -name "*.amb"), exactly like upstream BWA_MEM

The when-gated convenience rules live in modules/reference.oxoflow; the consumer rules (HOMER, deepTools, BWA-MEM) pick them up via depends_on (a closed gate counts as satisfied, so the default pre-built path is unchanged) and resolve the file at shell level, mirroring the khmer genome-size wiring.

Input data: single-end raw/<sample>.fastq.gz reads, named in [[sample_groups]] (see test/fixtures/ for a tiny working set); the paired-end branch (config.paired=true) reads raw/<sample>_1.fastq.gz + raw/<sample>_2.fastq.gz.

Compute: up to 12 CPUs / 72 GB per rule (trimming, alignment and deepTools rules are the heaviest); a few rules need as little as 1 CPU / 6 GB.

Tools: a mixed delivery — 46 of 49 rules run in pinned Docker images (biocontainers/* tags, e.g. biocontainers/macs2:2.2.7.1--py38h4a8c8d9_3), executed by oxo-flow via Docker or Singularity; the remaining three rules (picard_mergesamfiles, picard_markduplicates, merge_replicates) share one pinned conda env at envs/picard-samtools.yaml (picard 3.0.0, samtools 1.17), which requires conda/mamba at runtime.

# 1. install oxo-flow (see Installation)
# 2. prepare data: <raw_dir>/<sample>.fastq.gz (single-end), see fixtures/
# 3. preview the plan
oxo-flow dry-run main.oxoflow
# 4. run
oxo-flow run main.oxoflow -j 8
# 5. run a subset
oxo-flow run main.oxoflow -t multiqc --samples first:1

Sample names are declared in [[sample_groups]] (S1, S2 in the fixture set) — add your own names there or point raw_dir at your data. Pipeline behaviour is tuned through [config]: macs_gsize (default 2.7e9; "" = khmer auto-estimate), narrow_peak, broad_cutoff, fragment_size, min_trimmed_reads, out_dir, and the skip_* toggles (skip_fastqc, skip_qc, skip_trimming, skip_plot_profile, skip_plot_fingerprint, skip_multiqc, skip_peak_annotation). All can be overridden on the CLI.

Each non-default branch below is [[include]]d from modules/ and gated by a when condition on a [config] key (defaults in parentheses). Nothing from a gated branch runs unless its key is set, and a dry-run/run shows exactly which branch a toggle activates.

Where an upstream branch runs by default (--save_reference, ataqv, Picard metrics, IGV, consensus/DESeq2, R QC plots), this port ships it off so the default plan stays the minimal main path; enabling it is one config flag. See Known divergences.

Upstream: nf-core/atacseq @ 2.1.2 (commit 1a1dbe52ffbd82256c941a032b0e22abbd925b8a), MIT license. Created 2026-08-15; this workflow may lag behind upstream releases. Upstream attribution in NOTICE.md.

Ported with upstream defaults: aligner=bwa, single-end, narrow_peak=false (broad peaks), no control. One row per upstream process; steps not ported are listed with reasons. when-gated rules carry the gate in the Notes column.

• Branches that upstream runs by default ship OFF here: ataqv, Picard metrics, IGV, consensus peaks/DESeq2 and the R QC plots all run on the upstream default path; this port gates them behind skip_ataqv/skip_picard_metrics/skip_igv/skip_consensus_peaks/ skip_peak_qc (default true = off) so the default plan stays the minimal main path. skip_preseq = true matches the upstream default (preseq is off there too). Enabling a branch is a single config flag.
• Reference inputs: upstream builds/derives the reference from --genome (iGenomes) at runtime; this port consumes pre-built files (reference, bwa_index prefix, gtf, gene_bed, tss_bed, chrom_sizes, optional blacklist). prepare_reference = true instead generates the BWA index and chrom sizes from the FASTA (the fixture already ships them, so the rules report up-to-date there).
• ataqv and Picard metrics are SE-only: both rules consume the single-end filtered mLb.clN BAM; the paired-end variants (over the mLb.flT orphan-removed BAM) are not ported — with paired=true the rules are skipped even when their skip_* gates are off.
• Alternative aligners are single-end only (when adds !config.paired): the paired branch is bwa-only, matching upstream's paired alignment options. Misconfigurations (e.g. paired=true with aligner="star") surface as validation warnings via missing inputs.
• PE requires bwa_index and produces bwa/library/ outputs: the paired branch's bwa_mem_pe emits the full read group (-R '@RG\tID:{sample}\tSM:$SM\tPL:ILLUMINA\tLB:{sample}\tPU:1', where SM strips a _T[0-9]* suffix from the sample) rather than upstream's minimal '@RG\tID:{sample}\tSM:{sample}'; the default SE path keeps the original @RG handling.
• Broad peaks are hardcoded: the port's macs2_callpeak and all consumers use --broad; narrow_peak = true is honoured by macs2_callpeak but the downstream rules in this port read *_peaks.broadPeak only.
• MultiQC config: path_filters/module_order were trimmed to the ported steps; preseq/featureCounts/ataqv/DESeq2 sections appear when their branches are enabled (the PE MultiQC adds the PE fastqc/trimgalore logs). Report comment points at the upstream pipeline.
• macs_gsize: upstream keys it off the read-length genome block and auto-estimates via khmer when the value is empty; the port exposes it as config.macs_gsize (default 2.7e9, the upstream GRCh37/38 @ 50 bp value). Setting macs_gsize = "" ports the upstream empty-value semantics: ref::khmer_uniquekmers estimates the size from the FASTA at config.read_length (upstream KHMER_UNIQUEKMERS) and macs2_callpeak reads it into --gsize. Divergence: upstream additionally consults the read-length keyed genome table (e.g. a non-empty macs_gsize there overrides khmer); the port uses the explicit config.macs_gsize value or the khmer estimate, nothing else.
• IGV session lists one merged mLb_* track set: upstream separates per-replicate and merged-replicate bigWigs/peaks into two IGV track sets (mLb_* vs mLb_clN_*); this port's igv rule scans merged_library/{bigwig,macs2/broad_peak} once, so merged and per-replicate files land in the same set (the merged files are included when config.merged_samples names them).
• featureCounts is unstranded: -s 0 is passed explicitly (the upstream default); a [SCI-FEATURECOUNTS-STRAND] preflight hint may appear for the gated rule — informational.
• size_factors/ stays in the workdir: upstream DESeq2_QC publishes the size_factors directory; the port leaves it in the rule workdir (not declared as an output) — documented, not lost.
• Legacy image pin: qce::multiqc_custom_peaks (modules/qc_extra.oxoflow) still pins the retired quay.io/nf-core/ubuntu:20.04 image (upstream's shell-only rule used a plain ubuntu container). It is kept as-is for byte-identical default behavior; the rule is gated off by default (multiqc_custom_peaks = false).
• Dry-run prints benign input ✗ lines on the default plan: {sample}.mLb.clN.sorted.bam is declared as an output of both bamtools_filter (per-sample chain, runs) and merge_replicates (merged-replicate chain, off with the default empty merged_samples), so a dry-run with no files on disk reports that path unresolved for 20 instances (per-sample consumers such as macs2_callpeak, bedtools_genomecov, frip_score, plotfingerprint, their transitive bigWig/peak consumers, and gated-off pe::/qce:: copies). oxo-flow dry-run still exits 0; a real run resolves the files.
• Consensus branch needs ≥ 2 samples: upstream filters the peak channel to size() > 1; with one sample the consensus rules would produce degenerate output.

Upstream's MERGED_REPLICATE_* branch (default ON) folds per-replicate libraries into one merged BAM per sample and drives the whole downstream chain (markdup-free: the merged BAM is NOT re-marked; bigWig, MACS2/HOMER, consensus/DESeq2) from it. This port implements the same analysis with oxo-flow's input_groups primitive — a declared pattern plus group_by that folds matching files the same way groupTuple() does.

How to use it — declare replicate samples with the upstream _REP\d+ suffix convention (S1_REP1, S1_REP2, …) in your samplesheet, then pass the base names of the replicate groups through merged_samples so the MultiQC/IGV aggregation rules wait for the merged files:

# samples_replicates.csv: name,samples rows — one group per sample is fine
#   cohort,S1_REP1
#   cohort,S1_REP2
#   cohort,S2_REP1
#   cohort,S2_REP2
oxo-flow run main.oxoflow -j 4 \
    --samples @test/fixtures/samples_replicates.csv \
    --arg merged_samples=S1,S2

The merge is gated by skip_merge_replicates (default false = ON, matching upstream) AND a non-empty merged_samples — both are needed for the merged chain to run, so a plain samplesheet with the default config never instantiates it. With no _REP\d+ samples the merge_replicates rule instantiates nothing and the default plan is untouched.

Engine requirements: the merged-replicate mode uses the input_groups primitive (Traitome/oxo-flow#231, unreleased at time of writing). On released engines the gate is inert — the default path runs unchanged (the downstream rules fall back to their declared per-sample inputs), and a replicate run fails fast with a clear message instead of misbehaving.

What happens: merge_replicates folds the per-replicate filtered BAMs ({sample}_REP{rep}.mLb.clN.sorted.bam) by base id and writes the merged BAM + index + samtools stats to the canonical results/bwa/merged_library/{sample}.mLb.clN.sorted.bam path — the same location the per-sample chain writes — so the regular chain's rules (macs2_callpeak, homer_annotatepeaks, frip_score, bedtools_genomecov, ucsc_bedgraphtobigwig, deeptools_plots, plotfingerprint, multiqc_custom_peaks, and the PE equivalents) run on the merged BAM with their identical commands, exactly like upstream's MERGED_REPLICATE_CALL_ANNOTATE_PEAKS / ..._BAM_BEDGRAPH_BIGWIG / ..._PLOT_DEEPTOOLS chain. multiqc/igv aggregate both the per-replicate and the merged files (upstream does the same; see the MultiQC/IGV divergences).

Deviations from upstream (all documented, none change the default non-merge path):

• No re-markdup on the merged BAM: upstream's MERGED_REPLICATE_MARKDUPLICATES_PICARD re-runs MarkDuplicates with --REMOVE_DUPLICATES true on the merged BAM. Here the merged BAM keeps the per-replicate markdup state (mLb.clN, --REMOVE_DUPLICATES false), matching the port's single-library handling.
• Single-replicate groups are merged too: upstream keeps only groups with ≥ 2 replicates (groupTuple() size() > 1); input_groups instantiates every group it finds, so a lone S1_REP1 still produces a merged S1 BAM (MergeSamFiles handles a single input fine). Declare _REP\d+ names only when you intend the merge.
• Merged naming uses the library suffix (mLb), not mRp: upstream names merged-replicate files {id}.mRp.clN.sorted.bam; here the merged BAM reuses the canonical {sample}.mLb.clN.sorted.bam name, which is what lets the regular chain consume it unchanged.
• One consensus analysis: upstream runs the consensus/DESeq2 chain twice (per-replicate and merged-replicate peak sets). Here the glob-based cons::* rules aggregate per-replicate AND merged peaks in a single run — two sample sets in one consensus matrix instead of two separate analyses.
• Do not mix a plain sample with its replicate base: declaring both S1 and S1_REP1 makes bamtools_filter and merge_replicates write the same S1.mLb.clN.sorted.bam — a declared-sample conflict (the duplicate would be caught at validation). Use _REP names OR the plain name for one sample, never both.

bash test/run.sh

Runs oxo-flow validate + lint + dry-run (sample selection comes from the workflow's [[sample_groups]] declaration) against the fixture data; a debug pass additionally asserts that no literal {wildcards} leak into expanded commands. See test/run.sh for details.

Apache-2.0. Copyright (c) 2026 oxo-flow-community. Upstream (nf-core/atacseq) is MIT — see NOTICE.md and LICENSE.upstream.