From ATAC-seq Libraries to a Testable Count Matrix
A differential accessibility analysis asks whether the amount of ATAC-seq signal at a genomic region differs systematically between biological groups. The analysis must distinguish a biological difference from a technical one: a sample with more usable fragments can have larger counts at many regions even when its underlying accessibility is unchanged. It must also distinguish a consistent group difference from variation among biological replicates.
The starting point is a common set of genomic regions, such as a consensus set of accessible regions obtained from peak calls across the samples or a set defined independently. For each sample and each region, count the eligible ATAC-seq fragments overlapping that region, using one consistent counting rule. The result is a matrix with regions as rows and samples as columns. Each entry is a nonnegative integer count. The comparison is made at the same regions in every sample; a peak being called in one group and not another is not, by itself, a statistical test of differential accessibility.
A biological replicate is an independently sampled biological unit, not an additional sequencing lane from the same library. Technical sequencing depth improves measurement precision but does not substitute for biological replication. A design with no independent replication cannot estimate between-unit biological variability reliably, so its apparently precise counts do not support the same population-level inference.
Negative-Binomial Model and Library-Size Offset
For region \(i\) and sample \(j\), a standard model treats the count as negative binomial with mean \(\mu_{ij}\) and dispersion \(\phi_i\). Under this parameterization, the variance is \(\mu_{ij}+\phi_i\mu_{ij}^{2}\). The extra variance term allows biological variability to exceed the variance of a Poisson model, whose variance equals its mean.
Here \(s_j\) is a positive library-size normalization factor for sample \(j\), and \(\log(s_j)\) is an offset: its coefficient is fixed at one rather than estimated. The vector \(\mathbf{x}_j\) records the design variables for that sample, and \(\boldsymbol{\beta}_i\) contains the region-specific regression coefficients. A factor of two in \(s_j\) represents twice the effective sequencing exposure, so equal underlying accessibility would be expected to produce twice the count. Factors are usually scaled to have a convenient common reference, since a shared multiplicative rescaling can be absorbed by the intercept.
Size factors are not necessarily raw total reads. They should reflect usable library depth and be estimated consistently across samples, often using a robust method based on the distribution of counts across regions. Such normalization relies on a reference assumption: most regions, or an appropriate set of reference regions, do not undergo a coordinated change that would make the typical sample-wide count shift biological rather than technical. If that assumption is implausible, ordinary library-size normalization can misrepresent the comparison. In particular, it does not by itself solve the problem of a global accessibility change.
The usual inferential assumptions are that biological samples are independent conditional on the design, that the count distribution is adequately represented by the negative-binomial mean-variance relationship, and that the design matrix can distinguish the condition effect from included covariates. The selected region set and counting rule should be consistent across samples. Dispersion estimates are generally stabilized using information across regions, because estimating one dispersion from only a few replicates at every region is noisy.
Worked Example: Correcting Counts for Different Library Sizes
Suppose two hypothetical samples have size factors \(s_1=0.80\) and \(s_2=1.20\). At one fixed region their counts are \(Y_{i1}=80\) and \(Y_{i2}=120\). Assume the factors are valid relative exposure measures and the samples are otherwise comparable.
The normalized counts, used here for a transparent descriptive comparison, are \(80/0.80=100\) and \(120/1.20=100\). Thus the raw counts differ by 50%, but after accounting for the 50% difference in size factors, the estimated relative count is equal. The corresponding log-rate contrast is \(\log(120/1.20)-\log(80/0.80)=\log(100)-\log(100)=0\).
This calculation illustrates the role of the offset; it is not a complete differential test. With replicated samples, the model also estimates uncertainty from biological variation. Nor should normalized counts be treated as independent observations in an ordinary test that ignores the count model.
Condition Effects, Covariates, and Tests
For a two-group comparison, let \(x_j=0\) for the reference group and \(x_j=1\) for the group of interest. With no other covariates, the model is \(\log(\mu_{ij})=\log(s_j)+\beta_{0i}+\beta_{1i}x_j\). The coefficient \(\beta_{1i}\) is the log ratio of expected counts between groups at equal library size. Its exponent, \(\exp(\beta_{1i})\), is the fold change. A positive coefficient indicates higher expected accessibility in the group coded one; a negative coefficient indicates lower expected accessibility.
If batch, sex, or another measured variable could explain systematic count differences, include it in the design when scientifically appropriate and when the study design permits its effect to be estimated. For example, adding a batch indicator gives \(\log(\mu_{ij})=\log(s_j)+\beta_{0i}+\beta_{1i}x_j+\beta_{2i}b_j\). Then \(\beta_{1i}\) compares conditions at the same modeled batch value. Covariate adjustment does not repair perfect confounding: if every control is in batch one and every treated sample is in batch two, condition and batch cannot be separated from these data.
A Wald test evaluates \(H_0:\beta_{1i}=0\) using the fitted coefficient divided by its standard error. A likelihood-ratio test instead compares the fitted model containing the condition term with a nested model that omits it. Both test whether the data support a condition effect under the specified model; the Wald test directly uses an estimated coefficient and standard error, while the likelihood-ratio test measures the loss of fit when the term is removed. In a genome-wide analysis, many regions are tested, so p-values should be adjusted for multiple testing. A false discovery rate procedure, such as Benjamini-Hochberg, controls the expected proportion of false discoveries among the rejected hypotheses under its stated conditions; an adjusted p-value is not the probability that a particular region is false.
Worked Example: Interpreting a Fitted Condition Effect
In a hypothetical replicated comparison, suppose the negative-binomial model gives a condition coefficient \(\hat{\beta}_{1i}=0.588\) with standard error \(0.200\). Assume the samples are independent biological replicates, the library-size factors are appropriate, and the design includes any necessary estimable covariates. The coefficient is on the natural-log scale.
The estimated fold change is \(\exp(0.588)=1.80\). Thus, at this region, the fitted expected fragment count in the group coded one is 1.80 times that in the reference group, or 80% higher, after accounting for the offset and modeled covariates.
The Wald statistic is \(z=0.588/0.200=2.94\). Using a two-sided standard normal reference distribution gives \(p\approx0.0033\). A 95% coefficient interval is \(0.588\pm1.96(0.200)\), or \((0.196,0.980)\). Exponentiating the endpoints gives an approximate 95% fold-change interval of \((1.22,2.66)\). Before multiple-testing adjustment, the interval excludes a fold change of one and the test provides evidence of a positive condition effect. It does not establish that the region is significant after genome-wide adjustment; that depends on the full set of tests.
Worked Example: Why an Unbalanced Covariate Can Change the Comparison
Consider four hypothetical samples. The two controls are in batch A and the two treated samples are in batch B. A model with both condition and batch terms would require separate coefficients for these factors. Here the condition indicator and batch-B indicator have exactly the same values in all four samples: both are zero for controls and one for treated samples.
Consequently, the model’s linear predictor contains \(\beta_{1i}x_j+\beta_{2i}b_j=(\beta_{1i}+\beta_{2i})x_j\). The data can estimate the sum of the condition and batch effects, but cannot determine the two coefficients separately. The design matrix is rank-deficient for these terms, so there is no unique adjusted condition estimate. Dropping batch would yield a group contrast, but that contrast would combine condition and batch rather than isolate condition. The appropriate conclusion is that this design cannot identify a condition-specific effect without additional assumptions or data that break the confounding.
Justification for the Offset and Diagnostics
The offset follows from a simple exposure argument. Let \(\lambda_{ij}\) be the expected count per unit of effective sequencing exposure at region \(i\) in sample \(j\). If counts scale proportionally with exposure, then \(\mu_{ij}=s_j\lambda_{ij}\). Taking logarithms gives \(\log(\mu_{ij})=\log(s_j)+\log(\lambda_{ij})\). A regression model for the log rate therefore includes \(\log(s_j)\) with its coefficient fixed at one. This adjusts the expected count for exposure without pretending that sequencing depth is a biological condition effect.
Diagnostics should address both the data and the model. Check that samples have adequate usable fragments, that the count matrix contains the intended regions and consistent counts, and that size factors are not driven by a few unusual samples. Examine whether condition groups separate because of batch or another covariate, and whether any group has too few independent replicates for reliable dispersion estimation. Very low-count regions carry little information about fold change; filtering them by a rule chosen without reference to the observed condition difference can improve power by removing tests that are unlikely to be informative.
After fitting, inspect dispersion estimates, fitted-versus-observed patterns, influential samples, and the distribution of residual or test statistics for signs of systematic misfit. Large changes in estimated effects after removing one sample warrant investigation, not automatic deletion. Compare raw and normalized sample-level summaries: a strong relationship between apparent group differences and library size suggests inadequate normalization or a design problem. Also verify effect direction against the stated reference group and check whether reported results are adjusted for multiple testing.
Alternatives, Limitations, and Common Pitfalls
A Poisson model is simpler, but its mean-equals-variance assumption is often too restrictive for biological count data and can produce standard errors that are too small when biological variation is present. A negative-binomial model accommodates extra-Poisson variation. Testing normalized counts with a simple t-test can be useful in settings where its assumptions are defensible, but it does not naturally model the discrete count mean-variance relationship and small-sample uncertainty as the negative-binomial framework does. Wald and likelihood-ratio tests are alternatives within a fitted regression model, not substitutes for a sound design or normalization.
The conclusions are conditional on the region definition, counting rule, normalization reference, and model. A consensus peak set is useful for comparing the same features, but regions chosen after looking for the strongest group separation can inflate evidence if that selection is not accounted for. Differential accessibility identifies a statistical change in measured fragment counts; it does not by itself show that a gene is regulated, establish a mechanism, or prove a causal effect. Effect size, uncertainty, replicate consistency, and multiplicity-adjusted evidence should all inform interpretation.
- Pitfall: treating raw count ratios as biological fold changes. A larger library can have larger counts across regions. Compare using the offset or a valid normalization method, not raw counts alone.
- Pitfall: treating a nonsignificant result as proof of no change. A large p-value can reflect low power or imprecise estimates. Inspect the effect estimate and interval as well.
- Pitfall: claiming an adjusted condition effect under perfect confounding. If condition and batch have identical patterns, their effects are not separately identifiable. State the limitation rather than assigning the combined difference to condition.
Check Your Understanding
Answer each question, then compare with the solutions below.
- A region has counts 60 and 90 in two samples with size factors 0.60 and 1.50. What are the normalized counts, and what does their comparison suggest descriptively?
- Why can sequencing two lanes from one biological sample not replace two independent biological replicates?
- A fitted condition coefficient is \(-0.693\). What fold change does it imply, and which group has lower expected counts if the coded group has \(x=1\)?
- What can be concluded about separate condition and batch effects when their indicators are identical for every sample?
- Why should a genome-wide set of region-level p-values be adjusted for multiple testing?
Solutions. (1) The normalized counts are \(60/0.60=100\) and \(90/1.50=60\); descriptively, the first sample has the larger count relative to its size factor, but this alone is not a replicated test. (2) Lanes from one library are technical measurements of the same biological unit and do not estimate between-unit biological variability. (3) \(\exp(-0.693)\approx0.500\), so the coded group has about half the expected count of the reference group, conditional on the model. (4) The effects cannot be separated because the design is rank-deficient; only their combined contrast is identifiable. (5) Testing many regions creates opportunities for small p-values by chance, so adjustment is needed to control a specified error rate across the collection of discoveries.