Differential Histone Mark Enrichment

Why Differential Histone-Mark Enrichment Needs Its Own Design

A histone-mark assay measures where marked chromatin is enriched, but the genomic pattern can differ greatly by mark. Some marks produce relatively focal peaks; others extend across long domains or gene bodies. A differential analysis must therefore answer two questions together: what genomic features are being compared, and what is the unit of statistical replication? A test that treats adjacent genomic bins as independent biological samples can report impressive significance while providing no valid evidence about condition differences across samples.

The earlier material on chromatin sequencing established that sequencing counts and continuous signal tracks are related but not interchangeable summaries. Here the goal is to estimate and test condition-associated enrichment at a defined set of genomic regions. A significant difference is evidence of a change in measured enrichment under the model; it is not, by itself, proof that the mark caused a change in gene expression or phenotype.

Define the Features Before Testing

For a focal mark, such as promoter-associated H3K4me3, a region may be a peak interval or a fixed promoter window. For a broad mark, such as H3K27me3 or gene-body-associated H3K36me3, one peak boundary may not adequately represent the enrichment pattern. Possible choices include fixed-width bins, annotated gene bodies, or broad domains defined using a condition-blind or pooled procedure. The choice changes the biological question: a bin-level test asks about local enrichment, whereas a domain-level summary asks about enrichment over a larger interval.

The features should not be selected only because they look enriched in one condition and then tested on the same evidence as though the selection had been independent. That can bias effect estimates and p-values. A defensible analysis defines candidate regions from external annotations, a pooled set of samples, or another procedure that does not favor either condition; alternatively, selection and testing must be handled by a method that accounts for their dependence. The same region boundaries should then be used to count or summarize signal in every replicate.

The biological sample, not a read, bin, or technical library, is the replicate for inference about a condition. If a broad domain is divided into ten bins, those ten measurements from one sample do not become ten independent biological replicates. They can be combined into a region-level summary or modeled with a method that explicitly represents their dependence.

Count-Based Testing

A common approach is to count aligned fragments overlapping each fixed region in each biological sample. Counts are nonnegative integers and typically show more variability between biological replicates than a simple Poisson model allows. A negative-binomial model accommodates this extra variation. For region \(i\) and sample \(j\), one form is:

Here \(Y_{ij}\) is the count, \(\mu_{ij}\) its expected value, and \(\phi_i\) the region-specific dispersion. With a log link, a two-condition model can be written:

The size factor \(s_j\) is a sample-level exposure adjustment, and \(C_j\) indicates condition (for example, 0 for control and 1 for treatment). Covariates such as batch can be added when appropriate and supported by the design. Under this coding, \(\beta_{1i}\) is the estimated log rate ratio for treatment versus control at region \(i\); the estimated fold change is \(\exp(\beta_{1i})\). Testing no condition difference means testing \(H_0:\beta_{1i}=0\), against a specified alternative such as \(H_a:\beta_{1i}\ne 0\).

This model assumes independent biological samples conditional on the modeled covariates, comparable feature definitions, and a suitable mean–variance relationship. Its size-factor adjustment also relies on the normalization strategy being appropriate for the experiment. If a condition causes a genome-wide shift in chromatin signal, ordinary relative library scaling may not identify absolute changes; that issue requires special care rather than an automatic interpretation of the fitted fold changes.

Broad Marks and the Choice of Summary

For a broad mark, deciding how to represent the signal is part of the analysis rather than a cosmetic preprocessing choice. A region spanning a long domain may be divided into consistent bins and then summarized, or a pre-defined domain may be counted as a whole. Whole-domain counts can be easier to interpret as total enrichment over the interval, while bin-level estimates can reveal localized changes that an aggregate would hide. In either case, replicate-level measurements must be retained for inference.

Continuous Signal and Replicate-Level Comparisons

Some analyses use a continuous signal summary rather than raw region counts. For example, a fixed region may be summarized by mean signal, integrated signal, or a replicate-level score from a signal track. A comparison of such summaries can use a linear model if its assumptions are reasonable: independent biological replicates, approximately normal model residuals, and an adequate variance specification. A two-sample pooled t-test additionally assumes equal population variances. If that assumption is not reasonable, a Welch comparison or a suitable regression model with unequal-variance handling is preferable.

Signal tracks can be useful for visualization and for describing shape within a broad region, but values from tracks are often transformed, smoothed, or scaled. Their units and processing must be stated. They should not be treated as raw independent observations just because a track contains many genomic positions. For count-based inference, applying a test to plotted or rounded signal values discards the count model and can produce misleading standard errors.

Multiple Testing, Diagnostics, and Alternatives

A genome-wide analysis tests many regions. Even when every null hypothesis is true, some small p-values will occur by chance. The Benjamini–Hochberg procedure is commonly used to control the false discovery rate under independence or suitable positive dependence among tests: it adjusts the collection of p-values, and results should be reported with adjusted values (often called q-values) alongside effect estimates. Under other dependence structures, use a method with a guarantee appropriate to that structure. A threshold on unadjusted p-values alone does not control the expected proportion of false discoveries among reported regions.

Useful diagnostics include checking replicate-level counts or signal summaries, library quality, dispersion estimates, residual patterns, and whether results are dominated by one sample. Plotting the signal across a broad domain can reveal whether an aggregate difference reflects a coherent shift or a small subregion. Also inspect whether candidate regions were defined without condition-favoring selection and whether the design separates condition from batch. Earlier material on input controls is relevant: background enrichment and assay-specific artifacts can affect the measured signal, and subtracting input counts naively is not a general substitute for a model that handles background appropriately.

The main alternatives answer somewhat different questions. A negative-binomial count model is suited to replicate-level region counts and their extra-Poisson variation. A linear model for transformed counts or continuous signal can be useful when transformation, residual behavior, and variance modeling are justified. A permutation test may be attractive when the sample size and design permit valid label permutations, but it cannot repair confounding or a condition-dependent feature-selection step. Peak calling separately in each condition and comparing the resulting peak lists is another tempting approach, but it can make boundaries and selection depend on the same observed differences being tested. None of these approaches removes the need for biological replication or an appropriate feature definition.

Common Pitfalls and Interpretation

• Counting bins as replicates. Bins from the same sample are correlated measurements, not independent biological samples. Preserve sample-level replication in the model.
• Using condition-specific boundaries without accounting for selection. A region present only because one condition was used to call peaks may appear different by construction. Use comparable, defensible intervals.
• Confusing a large signal with a precise difference. Effect size and uncertainty are distinct. Report the fold change or signal difference, its uncertainty or test result, and any multiple-testing adjustment.
• Interpreting a p-value as biological importance. A small p-value does not show that the change is large, functional, or causal. Consider magnitude, reproducibility, genomic context, and independent evidence.
• Assuming relative normalization detects absolute global change. Standard size-factor methods rely on suitable reference behavior. If chromatin changes broadly across the genome, relative comparisons may conceal or distort that shift.