ATAC-seq Analysis Workflow: From Raw Reads to Regulatory Hypotheses

ATAC seq (Assay for Transposase Accessible Chromatin with high throughput sequencing) maps open chromatin genome wide by probing DNA accessibility with a hyperactive Tn5 transposase. This guide walks you through the core analysis workflow from raw FASTQ files to testable regulatory hypotheses. It is intended for experimental biologists planning their first ATAC seq project, bioinformatics beginners seeking a structured pipeline, and experienced analysts who want a concise reference. The focus is on bulk ATAC seq from a model organism or human sample, though many principles transfer to single cell approaches. For a thorough introduction to the biology and method, see the NCBI Bookshelf [1]. The EMBL EBI Training [2] provides additional context on sequencing data processing and quality assurance.

At a Glance

Quality Control and Preprocessing

Raw sequencing data from the NCBI Sequence Read Archive [5] or your own sequencer must pass basic quality screens. Run FastQC on the original FASTQ files to check per base quality scores, GC content, adapter contamination, and duplicate levels. If adapter sequences are present (common in short fragment ATAC seq libraries) trim them with Cutadapt or Trimmomatic. Perform a second FastQC round on trimmed reads to confirm improvement.

A key decision is whether to remove reads mapping to the mitochondrial genome. Mitochondrial chromatin is highly accessible and can consume 50% or more of reads. The Galaxy Training Network [3] recommends filtering out mitochondrial reads after alignment to reduce noise. Some pipelines remove them at the FASTQ stage by aligning to a combined reference without the mitochondrial chromosome, but post alignment filtering with SAMtools is safer. At this stage also check for excessive PCR duplicates not attributable to biological copy number. High duplication rates suggest low library complexity and can bias downstream signals. Record the proportion of reads retained after each filter.

Read Alignment

Align quality trimmed reads to the appropriate reference genome using Bowtie2, which is optimized for short reads and accepts the paired end format typical of ATAC seq. Use the --very sensitive preset and allow for discordant alignments if fragment sizes vary widely. The Bioconductor project [4] offers extensive documentation on handling aligned reads in R, but alignment itself is best done in a command line environment.

Remove duplicate reads with Picard MarkDuplicates. Deduplication is essential because Tn5 preferentially inserts into open chromatin and the same insertion site can be sequenced multiple times. Retain only properly paired reads with mapping quality (MAPQ) at least 30. Use SAMtools to filter and sort the resulting BAM file. For single end data (less common) treat each read independently but still enforce a minimum MAPQ.

Decision criteria: If your experiment uses paired end 50 bp reads you can safely use Bowtie2 default settings. For longer reads (100 bp or more) increase the seed length. If you plan to call peaks with MACS2, keep the alignment in BAM format with read pairs properly ordered.

Accessibility Signal and Quality Metrics

Before calling peaks, evaluate the quality of your ATAC seq library using fragment size distribution and transcription start site (TSS) enrichment. Use deepTools plotFingerprint or bamPEFragmentSize to generate a fragment size histogram. A healthy ATAC seq library shows a strong peak at around 50 to 100 bp (nucleosome free fragments) and smaller peaks at ~200 bp and ~400 bp representing mono and dinucleosomal fragments. Paired end reads allow precise fragment length calculation from the inferred insert size.

Compute a TSS enrichment score by collecting reads within a window around annotated TSSs and normalizing to local background. As described in bioRxiv [6], this metric quantifies signal at regulatory regions and should exceed 5 for most cell types. Low TSS enrichment indicates poor chromatin accessibility signal or high background. If your enrichment is below 2, consider whether the sample is degraded or the transposition reaction was suboptimal. The Galaxy Training Network [3] provides ready to use workflows for generating these metrics.

Peak Calling

Peak calling identifies regions of significantly enriched Tn5 insertions relative to a background model. The most widely used tool is MACS2, which handles paired end ATAC seq data well when you set the --shift parameter to compensate for the 9 bp duplication introduced by Tn5. For paired end data MACS2 automatically computes the fragment size from the alignment, do not use --shift and --extsize as you would for ChIP seq. Instead use MACS2 with -f BAMPE to treat each read pair as a single fragment. The EMBL EBI Training [2] has a dedicated module on peak calling with MACS2.

Other callers include Genrich (which integrates ATAC seq specific settings) and HMMRATAC (which models nucleosome positioning). For most bulk projects MACS2 with a q value cutoff of 0.05 produces reproducible narrow peaks. If you are studying broad accessible domains (e.g. in heterochromatin) consider using MACS2 --broad mode or Genrich. Always compare peak numbers across replicates and calculate overlaps using the Irreproducible Discovery Rate (IDR) method, especially if you plan to submit data to ENCODE or similar consortia.

Decision criteria: Use MACS2 narrow peaks for transcription factor motif analysis and broad peaks for enhancer landscapes. For samples with low read depth (< 20 million unique fragments) consider using Genrich or relaxed q values, but be aware of increased false positives.

Motif and Footprinting Analysis

Once you have a confident peak set, ask which transcription factors (TFs) may bind these accessible regions. Run motif enrichment with HOMER or MEME on the peak sequences. HOMER performs de novo motif discovery and matches known motifs against databases such as JASPAR. For network based transcription factor enrichment the Nucleic Acids Research [10] paper describes ChEA KG, which integrates knowledge graphs to prioritize TFs. The RAGER platform [9] provides an integrated workbench for combining ATAC seq and RNA seq data, making motif to target connections more robust.

Footprinting goes a step further by detecting short regions (6 to 20 bp) within peaks that are protected from Tn5 cleavage by bound proteins. Tools such as TOBIAS or HINT ATAC require high coverage (100 million reads or more) and precise alignment. As shown in bioRxiv [6], careful footprinting can uncover novel cis regulatory elements in Drosophila. For human or mouse, compare footprints with existing chromatin immunoprecipitation data to validate.

Biological Interpretation

Annotate peaks to nearby genes using ChIPseeker or GREAT. Define regulatory hypotheses based on proximity to TSSs, overlap with enhancer marks (H3K27ac, H3K4me1), and correlation with gene expression. The Cell [8] paper demonstrates single cell mapping of regulatory DNA protein interactions, which parallels bulk analysis at higher resolution. However, for bulk data you can integrate with public ChIP seq tracks in the UCSC Genome Browser.

Cautious interpretation is critical. Open chromatin indicates potential regulatory activity, not proof. A peak near a gene does not imply regulation of that gene. Consider co accessibility across a locus, evolutionary conservation, and functional validation. For instance, a motif for a repressor TF does not guarantee repression. Always express results as hypotheses: "The promoter of gene X becomes accessible upon stimulation, suggesting that TF Y may bind there and activate transcription." The RAGER platform [9] helps formalize such integration with statistical support.

Common Mistakes

• Ignoring mitochondrial reads. They can constitute the majority of aligned reads. Filter them out early.
• Using ChIP seq peak calling parameters. MACS2 default (-f AUTO) may misestimate fragment size for ATAC seq. Use -f BAMPE for paired end data.
• Overinterpreting motif presence. A motif may be bound in vitro but not in vivo. Combine footprinting and expression data.
• Skipping replicate analysis. Single sample peak lists often have high false discovery. Use IDR or at least a consensus approach.
• Neglecting batch effects. If samples are processed on different days or by different people, include batch variables in differential accessibility testing.

Limits and Uncertainty

ATAC seq signal depends on cell type, transposition efficiency, sequencing depth, and library preparation. Low coverage samples may miss accessible regions, especially in heterochromatin. The method does not distinguish between active enhancers and poised or repressed accessible sites. Co accessibility and chromatin state integration are necessary to infer functional activity.

Biological interpretation is also limited by the reference genome and annotation quality. For non model organisms, motif databases may be incomplete. Finally, ATAC seq does not provide single base resolution of protein binding, that requires footprinting or complementary methods like ChIP seq. Always treat peak calls as regions of interest, not definitive regulatory elements.

Frequently Asked Questions

How many reads do I need for a typical bulk ATAC seq experiment?

Aim for at least 50 million total reads per sample. After quality filtering and deduplication you should retain 20 to 30 million unique fragments. Low complexity samples (e.g. from tissues with many mitochondrial reads) may require more sequencing.

What is the best peak caller for ATAC seq?

MACS2 with -f BAMPE and default q value is a reliable starting point for most species. For broad domains or non model organisms, consider Genrich which has ATAC seq specific settings. Compare outputs from two callers if you have concerns.

Can I use this workflow for single cell ATAC seq?

No. Single cell ATAC seq requires specialized pipelines (e.g. ArchR, Signac, Cell Ranger ATAC) to handle sparse matrices, fragment files, and clustering. However, the interpretation of peaks, motifs, and integration with expression data is conceptually similar.

How do I handle low sample number (fewer than three replicates)?

With only one or two replicates you lose statistical power for differential accessibility. You can still call peaks and explore motifs, but avoid making strong comparisons. Consider using public data as biological replicates if available.

References and Further Reading

• NCBI Bookshelf ATAC seq overview and chromatin biology
• EMBL EBI Training Sequence analysis fundamentals
• Galaxy Training Network ATAC seq tutorial and workflows
• Bioconductor Genomic analysis packages for R
• NCBI Sequence Read Archive Public ATAC seq data sets
• bioRxiv Novel Drosophila cis regulatory elements from ATAC seq footprinting
• Nucleic Acids Research ChEA KG transcription factor enrichment tool
• PLoS One RAGER platform for integrated RNA and ATAC seq
• Cell Single cell mapping of regulatory DNA protein interactions
• Genome Med circVDJ seq for T cell clonotype detection

Related Articles

• RNA Sequencing Analysis: From FASTQ Files to Biological Questions
• RNA-seq Quality Control: What to Check Before Differential Expression
• How to Plan a Bulk RNA-seq Differential Expression Study
• Single-Cell RNA-seq Workflow: A Practical Analysis Roadmap
• Single-Cell RNA-seq Quality Control: Cells, Genes, and Mitochondrial Reads