gynecoloji/snakemake_ATACseq_spikein

Reproducible Snakemake workflow for spike-in–normalized ATAC-seq — concatenated-genome alignment, MACS2 peaks, a consensus fragment-count matrix, an interactive QC report, and DESeq2 differential binding. Docker/Apptainer-ready.

Overview

Latest release: v1.4.6, Last update: 2026-08-03

Share link: https://snakemake.github.io/snakemake-workflow-catalog?wf=gynecoloji/snakemake_ATACseq_spikein

Quality control: linting: failed formatting: failed

Topics: apptainer atac-seq atacseq bioinformatics deeptools docker epigenomics reproducible-research snakemake spike-in-normalization deseq2 differential-binding peak-calling

Deployment

Step 1: Install Snakemake and Snakedeploy

Snakemake and Snakedeploy are best installed via the Conda package manager. It is recommended to install conda via Miniforge. Run

conda create -c conda-forge -c bioconda -c nodefaults --name snakemake snakemake snakedeploy

to install both Snakemake and Snakedeploy in an isolated environment. For all following commands ensure that this environment is activated via

conda activate snakemake

For other installation methods, refer to the Snakemake and Snakedeploy documentation.

Step 2: Deploy workflow

With Snakemake and Snakedeploy installed, the workflow can be deployed as follows. First, create an appropriate project working directory on your system and enter it:

mkdir -p path/to/project-workdir
cd path/to/project-workdir

In all following steps, we will assume that you are inside of that directory. Then run

snakedeploy deploy-workflow https://github.com/gynecoloji/snakemake_ATACseq_spikein . --tag v1.4.6

Snakedeploy will create two folders, workflow and config. The former contains the deployment of the chosen workflow as a Snakemake module, the latter contains configuration files which will be modified in the next step in order to configure the workflow to your needs.

Step 3: Configure workflow

To configure the workflow, adapt config/config.yml to your needs following the instructions below.

Step 4: Run workflow

The deployment method is controlled using the --software-deployment-method (short --sdm) argument.

To run the workflow using apptainer/singularity, use

snakemake --cores all --sdm apptainer

To run the workflow using a combination of conda and apptainer/singularity for software deployment, use

snakemake --cores all --sdm conda apptainer

To run the workflow with automatic deployment of all required software via conda/mamba, use

snakemake --cores all --sdm conda

Snakemake will automatically detect the main Snakefile in the workflow subfolder and execute the workflow module that has been defined by the deployment in step 2.

For further options such as cluster and cloud execution, see the docs.

Step 5: Generate report

After finalizing your data analysis, you can automatically generate an interactive visual HTML report for inspection of results together with parameters and code inside of the browser using

snakemake --report report.zip

Configuration

The following section is imported from the workflow’s config/README.md.

Configuration

This workflow is configured through two files in this directory:

  • config.yaml — all workflow parameters (see below)

  • samples.csv — the sample sheet

plus reference data you download into ref/ (not tracked in git; see Reference data).

Sample sheet (config/samples.csv)

CSV with one row per sample and these columns:

column

description

sample_id

Sample name. Raw reads must be data/<sample_id>_R1_001.fastq.gz / _R2_001.fastq.gz.

type

Free-text condition label (e.g. Control, NICD3).

group

Replicate group. Reproducibility handling is chosen from group size (below).

Example:

sample_id,type,group
GSF4007-Control_1_S11,Control,group1
GSF4007-Control_2_S13,Control,group1
GSF4007-Control_3_S15,Control,group1
GSF4007-NICD3-V5_1_S12,NICD3,group2
GSF4007-NICD3-V5_2_S14,NICD3,group2
GSF4007-NICD3-V5_3_S16,NICD3,group2

Per-group reproducibility is derived automatically from the number of replicates in each group:

  • ≥ 3 replicates → majority vote (a peak is kept if it recurs in ≥ consensus_min_replicates replicates).

  • exactly 2 replicates → IDR (idr_threshold).

  • 1 replicate → the sample’s own peaks are used as-is.

Parameters (config/config.yaml)

Every parameter — with its type, default, and description — is defined once in the config schema, workflow/schemas/config.schema.yaml. That schema is the single source of truth: the workflow validates config.yaml against it on every run (and fills in defaults for anything you omit), and the Snakemake Workflow Catalog renders it as a parameter table on the workflow page.

To configure a run, edit config.yaml directly — it ships with working defaults and an inline comment on every parameter. At minimum, point the reference-file paths (human_fasta, spikein_fasta, blacklist, gtf, promoter_bed, enhancer_bed) at the files you provide (see Reference data).

Reference data

Genomes, indexes and large annotations are not shipped in the repo (they are .gitignored). Download / place them under ref/ before running, matching the paths in config.yaml:

  • ref/hg38.fa — chr-prefixed UCSC human genome

  • ref/dm6.fa (or another spike-in genome)

  • ref/hg38_blacklist_regions.bed — ENCODE hg38 blacklist (shipped)

  • ref/gencode.v36.annotation.gtf — GENCODE annotation (for TSS QC)

  • ref/hg38.2bit — for computeGCBias

  • ref/picard.jar — Picard (used by MarkDuplicates)

The combined Bowtie2 index (ref/COMBINED/) is built automatically by the build_combined_genome rule from human_fasta + spikein_fasta.

See the top-level README.md for full setup and run instructions.

Workflow parameters

The following table is automatically parsed from the workflow’s config.schema.y(a)ml file.

Parameter

Type

Description

Required

Default

samples_table

string

Path to the sample sheet CSV. Columns: sample_id, type, group. Per-group reproducibility (majority vote / IDR / single) is derived from group size.

yes

config/samples.csv

adapter_r1

string

Optional. Explicit R1 adapter sequence that OVERRIDES fastp auto-detection. Leave unset to auto-detect adapters for paired-end reads (–detect_adapter_for_pe).

adapter_r2

string

Optional. Explicit R2 adapter sequence (used together with adapter_r1).

human_fasta

string

Human genome FASTA. Must be chr-prefixed UCSC (chr1..chrX) to match the blacklist.

yes

ref/hg38.fa

spikein_fasta

string

Spike-in genome FASTA (any species, e.g. Drosophila dm6, yeast, E. coli).

yes

ref/dm6.fa

spikein_prefix

string

Prefix prepended to spike-in chromosome names before concatenation, so spike-in reads can be split out after alignment.

yes

spikein_

combined_index

string

Bowtie2 index prefix for the built human + spike-in reference (created automatically by the build_combined_genome rule).

yes

ref/COMBINED/genome

align_chroms

array

Human chromosomes kept when building the index (reads then align only to these plus the spike-in). Use an empty list to keep all human chromosomes.

yes

keep_chroms

array

Analysis keep-set for the final human BAM (mito-% QC is recorded first). Must be a subset of align_chroms and consistent with keep_chroms_regex.

yes

blacklist

string

ENCODE-style blacklist BED (chr-prefixed).

yes

ref/hg38_blacklist_regions.bed

peak_types

array

Peak types to analyze (narrowPeak for standard ATAC-seq; add broadPeak if needed).

[‘narrowPeak’]

effective_genome_size

integer

Effective genome size for deepTools RPGC normalization (hg38 default).

yes

2913022398

bin_size

integer

bigWig bin size in bp.

yes

25

consensus_window

integer

Fixed consensus peak width around each summit, in bp.

yes

500

consensus_min_replicates

integer

Majority-vote threshold for conditions with >=3 replicates.

yes

2

idr_threshold

number

IDR threshold for conditions with exactly 2 replicates.

yes

0.05

idr_relaxed_pvalue

number

MACS2 -p value for the relaxed peak calls used as IDR input.

yes

0.1

idr_top_n_peaks

integer

Number of top relaxed peaks retained per replicate for IDR.

yes

150000

keep_chroms_regex

string

Regex used by the consensus step to filter chromosomes; keep consistent with keep_chroms.

yes

^chr([1-9]

macs2_genome

string

MACS2 -g effective genome preset (e.g. hs, mm, ce, dm).

yes

hs

gtf

string

GENCODE GTF (chr-prefixed) used for TSS-enrichment QC.

yes

ref/gencode.v36.annotation.gtf

promoter_bed

string

Promoter BED used for the reads-in-annotation QC.

yes

ref/promoter_chr1-22X.bed

enhancer_bed

string

Enhancer BED used for the reads-in-annotation QC.

yes

ref/enhancer_chr1-22X.bed

spikein_pct_min

number

Lower bound (%) of the expected spike-in read fraction; drives the QC status flag.

yes

2

spikein_pct_max

number

Upper bound (%) of the expected spike-in read fraction; drives the QC status flag.

yes

10

spikein_min_reads

integer

Hard floor on spike-in reads per sample. compute_spikein_factors raises an error rather than emitting normalization factors if any sample falls below this. Factors derived from counts this low are dominated by Poisson noise (SD = sqrt(n), so ~22% per-sample error at 20 reads) and produce large numbers of false differential peaks. Unlike spikein_pct_min/spikein_pct_max, which only set a QC status flag, this gates the primary and differential stages. Set 0 to disable.

100000

spikein_max_within_spread

number

Maximum tolerated fold-spread of spike-in reads between replicates of one condition. Replicates of a condition received the same spike-in, so spread among them is measurement error rather than biology. The spikein and anchor_shape modes refuse to run above this, because size factors that noisy destroy statistical power: on a benchmark dataset with a 3.5x within-condition spread the spikein mode found 0 differential peaks where the spike-in-free baseline found 1,488. A read-depth floor does not catch this case – that dataset had 5.9-21.0 M spike-in reads (11.8-36.4% of the library).

2

ctcf_bed

string

Constitutive-CTCF anchor BED used by the ctcf differential-openness normalization and the anchor+shape hybrid. Built by workflow/scripts/build_constitutive_ctcf.py: the union of ENCODE CTCF ChIP-seq peaks across 59 cell types, restricted to regions bound in >=90% of them. Genome coordinates only, so it is reusable for any human sample.

ref/constitutive_ctcf_hg38.bed

diffopen_modes

array

Normalization modes run by the opt-in diffopen_all target, each written to its own results/diffopen// directory: none = DESeq2 median-of-ratios over all peaks, spikein = size factors from Drosophila spike-in depth, ctcf = median-of-ratios restricted to CTCF anchors (spike-in free), rnastable = median-of-ratios restricted to promoter-class peaks over the TSS window of RNA-seq-stable genes (needs diffopen_rna_table).

[‘none’, ‘spikein’, ‘ctcf’]

diffopen_ref_label

string

Reference level of the sample sheet’s type column (the denominator of the contrast), e.g. “Control” for a NICD3-vs-Control comparison.

Control

ctcf_trim_k

number

ctcf mode: MAD multiplier for the invariance trim — anchors whose between-condition shift is more than this many MADs from the median shift are dropped before re-estimating size factors. Large values disable it.

2.5

ctcf_trim_iter

integer

ctcf mode: number of trim / re-estimate iterations (0 = no trimming).

2

diffopen_rna_table

string

Path to the RNA-seq DESeq2/edgeR results table (TSV, or CSV by extension) used by the rnastable normalization. Required when rnastable is in diffopen_modes. Must contain the gene, log2FoldChange, padj and baseMean columns named by the diffopen_rna_*_col parameters.

diffopen_rna_gene_col

string

Column in diffopen_rna_table holding the gene symbol (matched to GTF gene_name).

gene

diffopen_rna_lfc_col

string

Column in diffopen_rna_table holding the log2 fold change.

log2FoldChange

diffopen_rna_padj_col

string

Column in diffopen_rna_table holding the adjusted p-value.

padj

diffopen_rna_basemean_col

string

Column in diffopen_rna_table holding mean expression.

baseMean

diffopen_rna_basemean_min

number

rnastable: minimum baseMean for a gene to be an eligible stable anchor (excludes underpowered low-count genes).

10

diffopen_rna_padj_min

number

rnastable: a gene counts as non-significant (stable) when padj >= this (or padj is NA). Higher = stricter “clearly unchanged” bar.

0.5

diffopen_rna_lfc_max

number

rnastable: maximum

log2FoldChange

for a gene to count as stable.

diffopen_rna_tss_window

integer

rnastable: half-width (bp) of the TSS window used to link a stable gene to peaks.

2000

diffopen_rna_min_anchors

integer

rnastable: refuse to normalize on fewer than this many anchor peaks.

100

diffopen_rna_promoter_class_required

boolean

rnastable: require anchor peaks to be promoter-class (Ensembl Regulatory Build) in addition to overlapping a stable-gene TSS window. Set false to relax to TSS-window overlap only when the strict set is too small.

true

rnastable_trim_k

number

rnastable: MAD multiplier for the invariance trim (as ctcf_trim_k).

2.5

rnastable_trim_iter

integer

rnastable: number of trim / re-estimate iterations (0 = no trimming).

2

anchor_shape_span

number

loess span for the intensity-dependent shape curve (anchor+shape hybrid).

0.6

anchor_shape_trim_k

number

MAD multiplier for trimming CTCF anchors that move between conditions, so only invariant sites shape the correction (anchor+shape hybrid).

2.5

anchor_shape_iter

integer

Shape-fit / anchor-trim refinement iterations (anchor+shape hybrid).

2

Linting and formatting

Linting results
1WorkflowError in file "/tmp/tmpt2z49o6a/gynecoloji-snakemake_ATACseq_spikein-114cd6b/workflow/Snakefile", line 48:
2Snakemake version mismatch: the driver running this workflow is 9.21.0, but workflow/envs/snakemake.yaml pins snakemake=9.3.2.
3Rules that use `script:` will fail with an empty log file.
4
5Fix by installing a matching driver, e.g.
6    mamba create -n atacseq -c conda-forge -c bioconda snakemake-minimal=9.3.2 pandas
7or by driving the workflow with the environment's own Snakemake.
8  File "/tmp/tmpt2z49o6a/gynecoloji-snakemake_ATACseq_spikein-114cd6b/workflow/Snakefile", line 59, in <module>
9  File "/tmp/tmpt2z49o6a/gynecoloji-snakemake_ATACseq_spikein-114cd6b/workflow/Snakefile", line 48, in _check_driver_version
Formatting results
 1[DEBUG] 
 2[DEBUG] 
 3[DEBUG] 
 4[DEBUG] 
 5[DEBUG] In file "/tmp/tmpt2z49o6a/gynecoloji-snakemake_ATACseq_spikein-114cd6b/workflow/rules/diffopen.smk":  Formatted content is different from original
 6[DEBUG] 
 7[INFO] 1 file(s) would be changed 😬
 8[INFO] 4 file(s) would be left unchanged 🎉
 9
10snakefmt version: 0.11.5