16  Case study: ATAC-seq

This chapter follows one ATAC-seq experiment through EpiNexus, from public files to biology. At each step it shows the result and how to read it: quality control, peaks, differential accessibility, motifs and Gene Ontology.

16.1 The question

SOX2 is a transcription factor that keeps mouse embryonic stem cells in their stem-cell state. It is a pioneer factor: it can bind DNA that is wrapped around nucleosomes, and help to open it (Maresca et al. 2023). If SOX2 keeps some sites open, those sites should close soon after SOX2 is removed. ATAC-seq measures where chromatin is open (Chapter 1), so it can test this.

We ask: which sites change their accessibility 2 hours after SOX2 is removed, and which transcription factors bind there?

16.2 The dataset

The data come from a study by Maresca and colleagues (GEO GSE209527) (Maresca et al. 2023). Their cells carry a SOX2 protein with a small tag (FKBPF36V). A molecule called dTAG-13 binds the tag and makes the cell destroy the tagged protein, within about an hour. We use four of their libraries (Table 16.1).

Table 16.1: The four ATAC-seq libraries used in this chapter. Read pairs from the European Nucleotide Archive.
Sample GEO Condition Read pairs
Timeseries_ATACseq_SOX2-FKBP_NT_Rep1 GSM6373853 untreated 27.0 million
Timeseries_ATACseq_SOX2-FKBP_NT_Rep2 GSM6373854 untreated 18.6 million
Timeseries_ATACseq_SOX2-FKBP_2H_Rep1 GSM6373859 2 h without SOX2 26.7 million
Timeseries_ATACseq_SOX2-FKBP_2H_Rep2 GSM6373860 2 h without SOX2 31.3 million

Each library was sequenced paired-end, 2 × 38 bp, in two runs. Together the files are about 4.9 GB. The study used the mouse genome mm10, so we use mm10 too.

This design is simple: two conditions, two replicates each. It also has a built-in check. We know which factor was removed. If the analysis works, the sites that close should carry the SOX2 motif.

NoteIn EpiNexus 1.2
  1. Create a project with the genome Mouse (mm10).
  2. On Upload Data, choose From a public repository and enter the four GSM numbers (Chapter 4). EpiNexus imports each sample’s two runs as one pair of FASTQ files.
  3. In the import table, rename the condition NT to untreated. EpiNexus treats labels such as untreated, control or WT as the reference, so fold changes read “2 h / untreated”.
  4. Start a run: ATAC-seq, Full Pipeline, genome mm10, and the default settings. These are MACS3 with q 0.05, IDR for the two replicates, and Library size (reads in peaks) normalisation. Regions count as differential at FDR < 0.1 and |log2 fold change| > 0.5.

The run used in this chapter is 1e7862c4, made with EpiNexus 1.2.1.

16.3 Quality control

Table 16.2: Quality control of the four libraries, from EpiNexus’s report (run 1e7862c4). Reads count both ends of each pair. TSS: TSS enrichment; Mito.: mitochondrial reads.
Sample Reads Mapped Duplicates FRiP TSS Mito. NRF PBC1
untreated, rep 1 53.8 M 93.6 % 27.0 % 45.1 % 9.5 1.5 % 0.58 0.65
untreated, rep 2 37.1 M 94.3 % 32.4 % 39.1 % 9.5 1.7 % 0.56 0.62
2 h, rep 1 53.1 M 94.2 % 33.6 % 43.3 % 9.7 2.0 % 0.52 0.59
2 h, rep 2 62.3 M 93.4 % 26.2 % 36.4 % 8.2 0.8 % 0.62 0.66

Some numbers look very good (Table 16.2):

  • FRiP is 36–45 %. ENCODE asks for more than 30 % (ENCODE Project Consortium 2026).
  • Mitochondrial reads are 0.8–2 %, far below the 15 % limit for the Omni-ATAC protocol (Chapter 6).
  • The fragment sizes show the expected pattern (Figure 16.1). The median fragment is 186–203 bp.
The ATAC-seq Analysis page: cards DARs 7,395 of 113,578 tested, more open 834, more closed 6,561; tab Fragment Sizes with a line chart of read count against fragment size from 10 to 990 bp for the four libraries. Each curve has a high peak near 60–70 bp, a second peak near 210 bp and a small bump near 390 bp.
Figure 16.1: The Fragment Sizes tab of the ATAC-seq page (run 1e7862c4). All four libraries show the nucleosome ladder: a peak below 100 bp from open DNA, a peak near 200 bp from DNA around one nucleosome, and a small one near 400 bp from two. The header counts 7,395 changed regions: 834 more open, 6,561 more closed. EpiNexus 1.2 interface; not covered by the book’s licence.

Other numbers look poor:

  • Duplicates are 26–34 %. More than 25 % “may suggest that the transposition was inefficient” (Grandi et al. 2022).
  • NRF and PBC1, the measures of library complexity (Chapter 6), are 0.52–0.66. ENCODE prefers values above 0.9 (ENCODE Project Consortium 2026).
  • TSS enrichment is 8.2–9.7. For mm10, ENCODE calls values below 10 “concerning” (ENCODE Project Consortium 2026). EpiNexus measures it with a different list of start sites than ENCODE, so the two scores are not directly comparable (Chapter 6).

EpiNexus 1.2.6 marks all four samples as failing QC, on duplicates, NRF and PBC1. It keeps them, and its report calls the downstream results “exploratory”.

Are the results unusable, then? Look at the other signs. The replicates agree: in the sample PCA (Figure 16.2), the first component holds 92.5 % of the variation and splits untreated from 2 h. The enrichment is strong. And, as the motif analysis below shows, the expected factor comes out on top. Low complexity means that the libraries hold fewer distinct DNA molecules than ideal. Deeper sequencing of these libraries would mostly read the same molecules again. It costs depth, not correctness.

Scatter plot of PC1 (92.5 %) against PC2 (4.7 %). Two blue points for 2 h without SOX2 sit on the left at PC1 about −58 and −39. Two orange points for untreated cells sit on the right at PC1 about 44 and 53.
Figure 16.2: Sample PCA of the four libraries, from the counts in the tested regions. The first component (92.5 % of the variation) separates the two conditions. Figure made by EpiNexus 1.2.1 (run 1e7862c4).
WarningPitfall: reading a failed metric as a failed experiment

A metric below its line is a question, not a verdict. Ask what the metric measures, and look at the other signs: enrichment (FRiP, TSS), replicate agreement (PCA) and a result you can predict. Here, low complexity sits next to strong enrichment, clear replicates and the expected motif.

16.4 Peaks

MACS3 found 83,954 to 112,309 peaks per sample. Two checks help to judge them.

Agreement with independent processing. ChIP-Atlas, a public database, processed the same four libraries with its own pipeline and a stricter threshold (q below 10−5) (Zou et al. 2024). It called 39,083 to 57,489 peaks per sample. EpiNexus found 89–91 % of them. The extra peaks in EpiNexus come from its looser threshold of 0.05.

Agreement between replicates. IDR keeps the peaks that rank high in both replicates (Chapter 7). ENCODE’s ATAC-seq standard says that “the number of peaks within an IDR peak file should be >70,000, though values >50,000 may be acceptable” (ENCODE Project Consortium 2026). EpiNexus’s IDR sets hold 32,934 peaks (untreated) and 34,098 (2 h). That is far below ENCODE’s line, although the replicates agree well:

  • 71,620 of the first untreated replicate’s peaks are also peaks in the second replicate. In the 2 h samples, the number is 78,413.
  • Take the strongest tenth of the first replicate’s peaks. Almost all of them (98.5 % untreated, 99.6 % at 2 h) are also found in the second replicate. Yet 17–18 % of these strong, shared peaks are missing from the IDR set.

Strong peaks found in both replicates should pass IDR. So, for this dataset, the IDR set in EpiNexus 1.2 misses many real sites.

NoteIn EpiNexus 1.2

The differential analysis does not use the IDR set by default. It builds its regions from each sample’s own peaks. For ATAC-seq data, keep that default: do not choose Replicate-reproducible peaks (IDR / overlap) under Consensus peaks for differential. Check the size of the IDR set against ENCODE’s numbers before you use it as your list of open sites.

16.5 Differential accessibility

EpiNexus made one set of regions from the four samples: windows of 501 bp around peak summits, found in at least two samples (Grandi et al. 2022). That gave 113,578 regions. It counted the reads in each region and tested the difference with DESeq2 (Love et al. 2014) (Chapter 9).

The result is clear (Figure 16.3). 6,561 regions lost accessibility and 834 gained it. Losses outnumber gains about eight to one. The lost regions also change more: 92 % of them lose more than half of their reads (log2 fold change below −1), and the strongest lose almost all of them. The gains are small: 98 % lie between log2 0.5 and 2.

This is what a pioneer factor predicts. Without SOX2, the sites it held open close. Few sites open, and the ones that do open change little.

Volcano plot of minus log10 FDR against log2 fold change. Thousands of blue points spread to the left, down to a fold change of about −8 and up to −log10 FDR of 28. A smaller group of orange points sits to the right, between fold changes 0.5 and 3, with −log10 FDR up to about 8. Grey points crowd the centre near zero.
Figure 16.3: EpiNexus’s differential result for 2 h without SOX2 against untreated cells: each point is a region of 501 bp. Blue: lost accessibility (6,561); orange: gained (834); grey: not significant (106,183; a sample of 5,000 is shown). Dashed lines: the thresholds (FDR 0.1, |log2 fold change| 0.5). Figure made by EpiNexus 1.2.1 (run 1e7862c4).

When most changes go one way: normalisation

Before comparing samples, a differential test divides each sample’s counts by a size factor, a number for how deeply that sample was sequenced. There are several ways to compute it:

  • Library size (reads in peaks), EpiNexus’s default for ATAC-seq: the number of reads in all tested regions.
  • Median-of-ratios (RLE), DESeq2’s default: for each region, the ratio of the sample’s count to a typical sample’s count, then the median of these ratios (Love et al. 2014).

The two agree when most regions do not change. They differ when many regions change in one direction, as here. Two hours after SOX2 removal, thousands of sites have closed. So the 2-hour samples have fewer reads in peaks. Scaling by reads in peaks lifts all their counts a little, to make the totals equal. That pushes every fold change slightly up. Weak gains near the threshold can then become significant. The median-of-ratios method looks at the typical region. A minority of changed regions does not move a median.

So, in an experiment like this one, treat the small set of gained regions with care. Some may come from the normalisation, not from the biology.

TipBest practice

Reske and colleagues compared eight ways to analyse differential ATAC-seq data. The normalisation method changed the number of differential regions and even the overall direction of change (Reske et al. 2020). Scaling by reads in peaks “should eliminate global differences in favor of reducing any technical biases”. They recommend that “researchers should systematically compare multiple normalization methods before continuing with differential accessibility analysis”. An MA plot (Chapter 9) shows quickly if one method shifts all regions up or down.

NoteIn EpiNexus 1.2

The normalisation is chosen on the Parameters step of the run wizard (Normalization). For ATAC-seq the default is Library size (reads in peaks). RLE (DESeq2 default) is the median-of-ratios method. If you expect many sites to change in one direction, use Run again with the other method, and compare the gained and lost regions of the two runs. The Differential page and the report state which normalisation a run used.

Figure 16.4 shows where the Differential page states it.

The Differential Analysis page for Timeseries SOX2-FKBP 2H versus untreated: a description naming DESeq2, normalisation by reads in peaks and 113,578 regions of 501 bp around peak summits, two orange notes about two replicates, four cards (113,578 total peaks, 7,395 significant, 834 gained, 6,561 lost) and an MA plot of log2 fold change against mean count, with a dense blue cloud below zero, grey points around zero and an orange band above.
Figure 16.4: The top of the Differential page for the ATAC-seq comparison (run 1e7862c4). The line Normalisation: reads in peaks names the method. In the MA plot, the lost regions (blue) form a large cloud below zero, and the gained regions (orange) a thin band above it. EpiNexus 1.2 interface; not covered by the book’s licence.

16.6 Motifs: which factors explain the changes?

EpiNexus looked for known motifs in the lost and gained regions. It compared them with regions of the same GC content that did not change (Chapter 12). The test is SEA (Bailey and Grant 2021), with the motifs of JASPAR 2024 (Rauluseviciute et al. 2024).

Table 16.3: Motifs in the lost and gained regions, from EpiNexus’s SEA results (q < 0.05 and enrichment above 1).
Regions Significant motifs Top motifs
Lost (6,561) 329 13 SOX and POU motifs share the top place, with q-values too small to compute. Among them are Pou5f1::Sox2 and POU2F1::SOX2. Pou5f1::Sox2 is in 42 % of lost regions and 11 % of unchanged ones. All of the top 20 motifs belong to the SOX or POU families.
Gained (834) 200 17 of the top 20 motifs are nuclear receptors (for example NR1H3, NR2F1 and ESRRB). No SOX motif is enriched.

The lost regions carry the SOX2 motif (Table 16.3). The strongest motifs include the joint motif of SOX2 and OCT4 (Pou5f1::Sox2). OCT4 is a factor that binds DNA together with SOX2. This is the built-in check we hoped for: the factor we removed explains the sites that closed. A second method in EpiNexus, monaLisa (Machlab et al. 2022), ranks all regions by their fold change. It finds the Pou5f1::Sox2 motif most enriched among the regions that lost the most.

The gained regions carry nuclear-receptor motifs. This is a lead, not a finding. The gained set is small, its changes are weak, and part of it may come from the normalisation (see above). A next step would be to check whether these regions still gain with the other normalisation, and whether nearby genes change their expression.

16.7 Gene Ontology: a pitfall with gene families

EpiNexus also tested the lost regions for Gene Ontology terms (Chapter 11). It found 89 significant terms. The top ones are “arachidonate epoxygenase activity” and “odorant binding”. These do not sound like stem-cell biology.

A look at the genes explains them. The first term comes from 25 Cyp genes, such as Cyp2c37 to Cyp2c40. The second comes from Mup and Vmn1r genes. These are gene families whose members lie side by side in the genome. GO enrichment of regions gives each region to the genes around it. So a handful of regions in one gene-rich area can light up a term shared by the whole family.

The GO table gives two tests for each term: one counts regions, and one counts genes. For lost regions, EpiNexus calls a term significant on the region test. Only 6 of the 89 terms also pass the gene test. “Odorant binding” does not: its gene-test FDR is 1.0.

WarningPitfall: GO terms from gene clusters

When a top GO term comes from one gene family, check the genes behind it, and check the gene-level test. A term carried by regions near one cluster of related genes says little about the biology of your experiment.

16.8 What this case study shows

  • The biology comes out clearly. Removing SOX2 closes thousands of sites within 2 hours, and the closed sites carry the SOX2 motif.
  • QC lines are guides. All four libraries failed EpiNexus’s complexity checks, yet the replicates agree and the expected factor comes out on top.
  • Normalisation matters when changes go one way. Treat a small set of changes in the other direction with care, and compare normalisation methods.
  • Check the IDR set before you use it. In EpiNexus 1.2 it was much smaller than ENCODE expects for this dataset.
  • Read GO terms through their genes. Gene families can produce enriched terms that say little about the experiment.

16.9 Further reading

  • Maresca et al. (2023): the study that produced these data, with a full time course and degron lines for OCT4 and NANOG.
  • Reske et al. (2020): how normalisation changes differential ATAC-seq results.
  • Grandi et al. (2022): the Omni-ATAC protocol and its quality checks.
  • ENCODE Project Consortium (2026): ENCODE’s ATAC-seq standards.