ATAC-seq data analysis: from raw reads to accessible chromatin maps.

1. Overview: what is ATAC-seq data analysis?

ATAC-seq analysis converts sequencing reads from transposase-accessible chromatin into maps of regulatory regions. The analysis identifies open chromatin peaks, evaluates data quality, annotates regulatory elements, discovers enriched motifs and compares chromatin accessibility between conditions.

2. ATAC-seq assay principle

ATAC-seq uses Tn5 transposase to insert sequencing adapters preferentially into accessible DNA. Open chromatin regions generate more reads than compact or nucleosome-protected regions. The fragment-size distribution can reveal nucleosome-free and nucleosome-associated fragments.

3. ATAC-seq study types

4. Experimental design

ATAC-seq is sensitive to nuclei preparation, cell viability, mitochondrial DNA, sequencing depth and batch effects. Plan biological replicates and controls before sequencing.

Questions to answer early

• Is the project bulk ATAC-seq, low-input ATAC-seq, single-cell ATAC-seq or multiome?
• Which biological conditions, tissues, time points or treatments are compared?
• How many biological replicates are available per group?
• Are samples balanced across library-preparation batches and sequencing runs?
• What reference genome, blacklist and gene annotation will be used?
• Is the main goal peak discovery, differential accessibility, motif discovery, enhancer mapping or integration with RNA-seq?
• Are there expected cell-composition changes that may influence bulk ATAC-seq signal?

5. Input files and metadata

sample_id	group	replicate	batch	tissue	fastq_1	fastq_2
A_rep1	control	1	A	cells	A_rep1_R1.fastq.gz	A_rep1_R2.fastq.gz
A_rep2	control	2	A	cells	A_rep2_R1.fastq.gz	A_rep2_R2.fastq.gz
B_rep1	treated	1	B	cells	B_rep1_R1.fastq.gz	B_rep1_R2.fastq.gz
B_rep2	treated	2	B	cells	B_rep2_R1.fastq.gz	B_rep2_R2.fastq.gz

6. FASTQ quality control

Raw-read QC detects adapter contamination, low-quality bases, uneven read counts and sequencing problems before alignment.

mkdir -p results/qc/fastqc results/qc/multiqc

fastqc data/fastq/*.fastq.gz \
  --outdir results/qc/fastqc \
  --threads 8

multiqc results/qc \
  --outdir results/qc/multiqc

• Inspect per-base quality, adapter content and overrepresented sequences.
• Compare read counts across samples.
• Review duplication in context of library complexity and enrichment.
• Expect short fragments; adapter trimming is often relevant.

7. Adapter trimming

ATAC-seq libraries often contain short inserts, so adapter sequences may be present in reads. Trimming can improve alignment and downstream QC.

mkdir -p results/trimmed results/qc/fastp

sample="A_rep1"

fastp \
  --in1 "data/fastq/${sample}_R1.fastq.gz" \
  --in2 "data/fastq/${sample}_R2.fastq.gz" \
  --out1 "results/trimmed/${sample}_R1.trimmed.fastq.gz" \
  --out2 "results/trimmed/${sample}_R2.trimmed.fastq.gz" \
  --html "results/qc/fastp/${sample}.html" \
  --json "results/qc/fastp/${sample}.json" \
  --thread 8

8. Alignment

ATAC-seq reads are commonly aligned to the reference genome using Bowtie2, BWA or similar short-read aligners. Paired-end data are especially useful because fragment sizes are biologically informative.

mkdir -p results/alignments results/logs

sample="A_rep1"

bowtie2 \
  -x reference/bowtie2/genome \
  -1 "results/trimmed/${sample}_R1.trimmed.fastq.gz" \
  -2 "results/trimmed/${sample}_R2.trimmed.fastq.gz" \
  -p 16 \
  2> "results/logs/${sample}.bowtie2.log" | \
  samtools sort -@ 8 -o "results/alignments/${sample}.sorted.bam"

samtools index "results/alignments/${sample}.sorted.bam"

9. BAM filtering and blacklist removal

Filtering removes reads that are unmapped, low quality, duplicates, mitochondrial reads or located in known problematic regions. Filtering rules should be documented because they influence peaks and QC metrics.

mkdir -p results/filtered

sample="A_rep1"

# Keep high-quality, properly mapped nuclear reads.
samtools view -b -q 30 -F 1804 \
  "results/alignments/${sample}.sorted.bam" \
  > "results/filtered/${sample}.mapq30.bam"

samtools index "results/filtered/${sample}.mapq30.bam"

# Remove blacklist regions.
bedtools intersect \
  -v \
  -abam "results/filtered/${sample}.mapq30.bam" \
  -b reference/blacklist.bed \
  > "results/filtered/${sample}.filtered.bam"

samtools index "results/filtered/${sample}.filtered.bam"

10. Mitochondrial reads

ATAC-seq often contains mitochondrial DNA because mitochondria are accessible to Tn5. High mitochondrial fraction can reduce usable nuclear reads and indicate sample-preparation issues.

samtools idxstats results/alignments/A_rep1.sorted.bam \
  > results/qc/A_rep1.idxstats.txt

11. Tn5 shift correction

Tn5 transposase inserts adapters with a characteristic offset from the true cut site. Shift correction is often applied before generating cut-site tracks, peak summits or footprinting inputs.

12. Insert-size and nucleosome periodicity QC

Paired-end ATAC-seq fragment sizes should often show patterns corresponding to nucleosome-free, mono-nucleosome and multi-nucleosome fragments.

samtools stats results/filtered/A_rep1.filtered.bam \
  > results/qc/A_rep1.samtools.stats.txt

multiqc results/qc --outdir results/qc/multiqc_alignment

13. TSS enrichment and FRiP

TSS enrichment and FRiP are two of the most useful ATAC-seq quality metrics.

14. Peak calling

Peak calling identifies accessible regions enriched for ATAC-seq fragments. MACS2 or MACS3 are commonly used for bulk ATAC-seq peak calling.

mkdir -p results/peaks

sample="A_rep1"

macs2 callpeak \
  -t "results/filtered/${sample}.filtered.bam" \
  -f BAMPE \
  -g hs \
  -n "${sample}" \
  --outdir results/peaks \
  --nomodel \
  --shift -100 \
  --extsize 200 \
  -q 0.01

15. Replicates and reproducibility

Biological replicate consistency is essential for trustworthy ATAC-seq. Technical replicates can be useful, but biological replication is needed for condition-level conclusions.

16. Consensus peak sets

A consensus peak set is often used for read counting, annotation and differential accessibility. It can be created from merged peaks across replicates or conditions.

mkdir -p results/peaks/consensus

cat results/peaks/*_peaks.narrowPeak \
  | cut -f1-3 \
  | sort -k1,1 -k2,2n \
  | bedtools merge \
  > results/peaks/consensus/consensus_peaks.bed

17. Signal tracks and genome-browser visualization

bigWig signal tracks allow visualization of chromatin accessibility across the genome and are useful for reports, genome browsers and manual review of key loci.

mkdir -p results/tracks

bamCoverage \
  -b results/filtered/A_rep1.filtered.bam \
  -o results/tracks/A_rep1.RPGC.bw \
  --normalizeUsing RPGC \
  --effectiveGenomeSize 2913022398 \
  --binSize 10 \
  --numberOfProcessors 8

• Use consistent normalization when comparing samples visually.
• Inspect representative promoters, enhancers and positive-control regions.
• Always interpret browser snapshots together with genome-wide statistics.

18. Peak annotation

Peak annotation connects accessible regions to promoters, enhancers, genes, CpG islands, repeats or custom regulatory features. Nearby genes are useful hypotheses, not automatic regulatory targets.

bedtools intersect \
  -a results/peaks/consensus/consensus_peaks.bed \
  -b annotation/promoters.bed \
  -wa -wb \
  > results/annotation/atac_peaks_overlapping_promoters.tsv

19. Motif analysis

Motif enrichment analysis identifies transcription-factor binding motifs enriched in accessible regions. It is especially useful for interpreting regulatory programs and differential accessibility.

mkdir -p results/motifs

bedtools getfasta \
  -fi reference/genome.fa \
  -bed results/peaks/consensus/consensus_peaks.bed \
  -fo results/motifs/consensus_peak_sequences.fa

• Use appropriate background sequences matched for GC content and accessibility context when possible.
• Separate promoter and distal peaks if they represent different regulatory contexts.
• Motif enrichment suggests candidate regulators; it does not prove binding without additional evidence.

20. Transcription-factor footprinting

Footprinting aims to detect local depletion of Tn5 insertions at transcription-factor binding sites. It can provide regulatory hypotheses but is sensitive to bias, sequencing depth and normalization.

21. Differential accessibility analysis

Differential accessibility analysis tests whether chromatin accessibility differs between groups at peak regions. It usually uses read counts over a consensus peak set and count-based statistical models.

22. Integration with RNA-seq, ChIP-seq and methylation

ATAC-seq is most powerful when integrated with other regulatory and expression data.

23. Note on single-cell ATAC-seq

Single-cell ATAC-seq requires specialized processing because reads are assigned to cell barcodes and converted to sparse peak-by-cell or tile-by-cell matrices.

24. Example ATAC-seq analysis workflow

The following simplified workflow illustrates a common paired-end bulk ATAC-seq route. Real projects should adapt parameters to organism, sample type, replicates and validation requirements.

# 1. Raw-read QC
fastqc data/fastq/*.fastq.gz --outdir results/qc/fastqc --threads 8

# 2. Trim reads
fastp \
  --in1 data/fastq/A_rep1_R1.fastq.gz \
  --in2 data/fastq/A_rep1_R2.fastq.gz \
  --out1 results/trimmed/A_rep1_R1.fastq.gz \
  --out2 results/trimmed/A_rep1_R2.fastq.gz \
  --html results/qc/fastp/A_rep1.html \
  --json results/qc/fastp/A_rep1.json \
  --thread 8

# 3. Align
bowtie2 -x reference/bowtie2/genome \
  -1 results/trimmed/A_rep1_R1.fastq.gz \
  -2 results/trimmed/A_rep1_R2.fastq.gz \
  -p 16 2> results/logs/A_rep1.bowtie2.log | \
  samtools sort -@ 8 -o results/alignments/A_rep1.sorted.bam

samtools index results/alignments/A_rep1.sorted.bam

# 4. Filter
samtools view -b -q 30 -F 1804 results/alignments/A_rep1.sorted.bam \
  > results/filtered/A_rep1.mapq30.bam

bedtools intersect -v \
  -abam results/filtered/A_rep1.mapq30.bam \
  -b reference/blacklist.bed \
  > results/filtered/A_rep1.filtered.bam

samtools index results/filtered/A_rep1.filtered.bam

# 5. Peak calling
macs2 callpeak \
  -t results/filtered/A_rep1.filtered.bam \
  -f BAMPE -g hs -n A_rep1 \
  --outdir results/peaks \
  --nomodel --shift -100 --extsize 200 -q 0.01

# 6. Signal track
bamCoverage \
  -b results/filtered/A_rep1.filtered.bam \
  -o results/tracks/A_rep1.RPGC.bw \
  --normalizeUsing RPGC \
  --effectiveGenomeSize 2913022398 \
  --binSize 10 \
  --numberOfProcessors 8

# 7. Summarize
multiqc results --outdir results/qc/multiqc_final

25. Deliverables and reporting

• FASTQ QC, trimming and final MultiQC reports.
• Sorted, indexed and filtered BAM files.
• Mitochondrial read fraction, duplicate rate, mapping rate and blacklist fraction.
• Insert-size periodicity, TSS enrichment, FRiP and peak-count summaries.
• Peak files in BED/narrowPeak format.
• Consensus peak set for annotation and differential accessibility.
• Normalized bigWig signal tracks for genome-browser visualization.
• Peak annotation tables for promoters, enhancers, genes and custom features.
• Motif enrichment and candidate transcription-factor summaries.
• Differential accessibility tables and plots when comparing conditions.
• Methods section with software versions, genome assembly, blacklist version, parameters and limitations.

26. ATAC-seq analysis cheat sheet

Useful ATAC-seq tools and documentation

Frequently asked questions