From assay output to defensible interpretation
An epigenomic comparison can produce thousands of methylation probes and genomic regions with apparent differences between conditions. A small p-value is not, by itself, evidence that a biological difference is real: low-quality samples, technical structure, cell composition, and large numbers of simultaneous tests can all affect the result. A defensible analysis links sample-level diagnostics to a clearly specified model and reports effect size and uncertainty alongside statistical evidence.
The worked analysis below uses hypothetical samples from two conditions, measured with a bulk methylation array and a chromatin-sequencing assay. The assays measure different features: the array reports methylation proportions at probes, while sequencing produces read counts or normalized signal over genomic regions. As emphasized in the earlier material on batch effects in methylation arrays, chromatin-sequencing QC, and differential chromatin testing, QC, normalization, and model choice must reflect the assay and study design.
Start with sample and design checks
Before testing features, verify sample identities, condition labels, assay matching, and the design matrix. Check assay-specific QC summaries rather than relying on a single universal threshold. For methylation arrays, inspect detection failures, control-probe performance, intensity distributions, and sample-level clustering. For chromatin sequencing, inspect read depth, mapping and duplication summaries, library complexity, enrichment or signal metrics appropriate to the assay, and replicate concordance. Thresholds should be chosen for the platform and protocol before interpreting differential results; a metric that is unusual is a reason to investigate, not an automatic reason to discard a sample.
Also inspect whether condition is entangled with batch. If every control was processed in one batch and every treated sample in another, the condition effect and batch effect cannot be separated from these data. Including both labels in a regression does not repair perfect confounding: the design matrix is rank deficient, so the separate effects are not identifiable. Partial overlap can permit adjustment, although it may leave little information for a precise estimate.
Worked Example: Applying QC Without Overreacting to an Outlier
In a hypothetical study, nine matched tissue samples were attempted on both assays. Eight pass prespecified, assay-specific checks. One sample has methylation-array control metrics within the laboratory’s acceptable range, but its chromatin library has unusually low usable read depth and enrichment, and its replicate relationships are poor. These observations suggest a chromatin-library problem; they do not establish that the sample’s biology is unusual.
The appropriate response is to check sample identity, library records, and the protocol-defined QC criteria, then document the decision. If the chromatin library fails a prespecified technical criterion, omit that sample from the chromatin analysis and report the reason. The acceptable methylation-array data need not automatically be discarded from a methylation-only analysis. A cross-assay comparison requiring matched samples should use the matched set available for that comparison and state the resulting sample count.
The other eight samples are distributed across two processing batches: batch 1 has three controls and one treated sample; batch 2 has one control and three treated samples. Both conditions occur in both batches, so a batch-adjusted contrast is estimable. The imbalance means the estimate relies on an additive model and will be less secure than a design with stronger balance; it is not equivalent to perfect condition–batch confounding.
Adjusted methylation differences
At one methylation probe, let \(Y_i\) be the observed beta value for sample \(i\), a proportion on the scale from 0 to 1. For this worked comparison, use the linear model \(Y_i=\beta_0+\beta_1 T_i+\beta_2 B_i+\varepsilon_i\), where \(T_i=1\) for treated samples and \(0\) for controls, and \(B_i=1\) for batch 2 and \(0\) for batch 1. The target is \(\beta_1\), the adjusted treated-minus-control difference in beta value.
This model assumes independent biological samples, an additive batch effect, and residuals with a common variance and an approximately normal distribution for the usual small-sample t interval. It also requires that condition and batch are not perfectly confounded. The additive assumption says the condition contrast is the same in both batches. If that assumption is implausible, an interaction may be scientifically relevant, but this small, imbalanced design offers limited information to assess it.
Worked Example: A Batch-Adjusted Methylation Contrast
The eight hypothetical beta values are as follows: batch 1 controls are 0.18, 0.20, and 0.22; the batch 1 treated sample is 0.35; the batch 2 control is 0.30; and batch 2 treated samples are 0.43, 0.45, and 0.47. The values lie within the beta-value range. There are samples from both conditions in both batches, so the additive model is estimable.
The unadjusted means are \(0.225=(0.18+0.20+0.22+0.30)/4\) for controls and \(0.425=(0.35+0.43+0.45+0.47)/4\) for treated samples. Their raw difference is \(0.425-0.225=0.200\), or 20.0 percentage points. This contrast mixes condition and batch because the batch proportions differ between conditions.
The batch-specific means are 0.20 for batch 1 controls and 0.35 for the batch 1 treated sample, a difference of 0.15. In batch 2, the means are 0.30 for the control and 0.45 for treated samples, also a difference of 0.15. The fitted additive model therefore gives \(\hat{\beta}_1=0.150\), or an adjusted difference of 15.0 percentage points. The fitted batch 2 coefficient is \(0.10\).
Residuals around the fitted values are \(-0.02,0,0.02\) among batch 1 controls and \(-0.02,0,0.02\) among batch 2 treated samples; the other two observations have zero residual. Thus \(SSE=4(0.02^2)=0.0016\), and the residual variance estimate is \(s^2=0.0016/(8-3)=0.00032\), with \(df=5\). For this design, the estimated variance of the condition coefficient is \(s^2/1.5\), so \(SE(\hat{\beta}_1)=\sqrt{0.00032/1.5}=0.01461\).
Using the two-sided 95% t critical value \(t_{0.975,5}=2.571\), the interval is \(0.150\pm2.571(0.01461)=(0.1124,0.1876)\). At this probe, the estimated adjusted difference is 15.0 percentage points, with a model-based 95% confidence interval from 11.2 to 18.8 percentage points. This interval describes uncertainty under the stated model; it is not a 95% probability that the fixed parameter lies inside this particular interval.
A short justification explains why adjustment matters. In an additive model, the fitted contrast can be obtained from either batch-specific treated-minus-control difference when both are available: both differences here equal 0.15. The raw group means instead weight the batches differently: controls are mostly in batch 1, while treated samples are mostly in batch 2. The observed raw difference is therefore not the same estimand as the batch-adjusted contrast. With more complex imbalance, the model uses all observations according to the design rather than simply subtracting pooled means.
Beta values are bounded, so a linear model can be unsuitable when values cluster near 0 or 1, residual spread changes markedly with the mean, or fitted values fall outside the allowable range. Inspect residuals and influence, and compare the result with a method suited to the measurement distribution when diagnostics are poor. Do not add covariates mechanically: adjustment for a variable affected by condition can change the scientific question or introduce bias.
Adjusted chromatin-sequencing signals
For read counts over a region, a common model is a negative-binomial generalized linear model. Write \(Y_{ij}\sim NB(\mu_{ij},\alpha)\), with \(\log(\mu_{ij})=\log(L_j)+\gamma_0+\gamma_1T_j+\gamma_2B_j\). Here \(Y_{ij}\) is the count for region \(i\) in sample \(j\), \(L_j\) is a library-size or normalization factor used as an offset, \(T_j\) indicates condition, and \(B_j\) indicates batch. The negative-binomial variance is \(\mu_{ij}+\alpha\mu_{ij}^2\), allowing extra-Poisson variation. In practice, dispersion is estimated or stabilized across regions using the chosen differential-testing procedure.
The model assumes independent biological replicates conditional on the covariates and a suitable mean–variance relationship. Normalization must be appropriate to the assay and its global signal behavior; an offset does not correct an unsuitable normalization strategy. Inspect library diagnostics and replicate concordance, and check that the design is estimable. A region’s raw counts alone are not a differential test.
Worked Example: A Chromatin Region With an Adjusted Fold Change
At one hypothetical region, four controls have normalized counts of 60 each and four treated samples have normalized counts of 120 each. Assume equal size factors, condition and batch are modeled as above, batch has no fitted effect, and a dispersion estimate of \(\alpha=0.10\) is available from the analysis across regions. The mean counts are 60 and 120, so the fitted log mean ratio is \(\hat{\gamma}_1=\log(120/60)=\log(2)=0.6931\). The corresponding fold change is \(\exp(0.6931)=2.00\).
Using the approximate Wald standard error for the log ratio, \(SE(\hat{\gamma}_1)=\sqrt{(1/120+0.10)/4+(1/60+0.10)/4}=0.2372\). The Wald statistic is \(z=0.6931/0.2372=2.921\), giving an approximate two-sided \(p\)-value of 0.0035. A 95% interval on the log scale is \(0.6931\pm1.96(0.2372)=(0.2282,1.158)\). Exponentiating its endpoints gives a fold-change interval of approximately \((1.26,3.18)\).
Under this illustrative model, the treated group has an estimated twofold signal at the region, with a model-based 95% interval from about 1.26-fold to 3.18-fold. This is evidence for a condition-associated difference at this region before accounting for testing many regions. The calculation depends on the assumed dispersion and equal normalization factors; actual analysis should use the fitted model and its estimated uncertainty, not substitute this approximation for a full analysis.
Multiplicity, cross-assay evidence, and interpretation
A genome-wide analysis tests many probes or regions, so unadjusted p-values overstate the number of compelling findings. A standard option is the Benjamini–Hochberg procedure. If there are \(m\) tests, order the p-values as \(p_{(1)}\leq\cdots\leq p_{(m)}\), and compare each with \(i\alpha/m\) at rank \(i\). The procedure controls the false discovery rate under independence and certain positive-dependence conditions; its guarantees should not be assumed under every possible dependence structure. Report the adjustment method, the number of tests, and adjusted values or a clearly stated threshold. A p-value is not the probability that the null hypothesis is true.
Methylation and chromatin results can provide complementary context, but they are not interchangeable measurements. A methylation difference at a probe and a chromatin signal difference at a nearby region do not automatically establish a regulatory mechanism. Compare genomic locations using a prespecified mapping rule, consider whether the signals are consistent in direction with the biological question, and report discordant as well as concordant findings. If cell composition differs between bulk samples, apparent molecular differences may partly reflect different cell mixtures; adjustment is informative only when the relevant composition measures and assumptions are adequate.
| Analysis choice | Useful when | Important limitation |
|---|---|---|
| Batch-adjusted linear model for beta values | The beta-value scale and residual behavior support an additive mean comparison. | Bounded values and small samples can make normal, constant-variance assumptions poor. |
| Negative-binomial model for sequencing counts | Replicate counts show overdispersion relative to a Poisson model and normalization is appropriate. | Results depend on the mean–variance and normalization assumptions, and require biological replication. |
| Unadjusted group comparison | Design is balanced or technical covariates are irrelevant to the contrast. | Can mix condition differences with batch differences in an imbalanced design. |
| Technical correction or sample exclusion | Diagnostics identify a credible technical problem and the correction or exclusion is justified. | Cannot recover a condition effect that is perfectly confounded with batch; data-driven exclusion can bias results. |
The most important failure mode is to interpret technical separation as biology. If sample clustering follows batch, revisit labels, QC, normalization, and estimability before interpreting feature-level tests. A large condition coefficient with a wide interval may be compatible with several effect sizes; a small p-value with a negligible effect may have little biological importance. Report effect estimates, uncertainty, adjusted evidence, sample counts, exclusions, and diagnostic concerns together.
Check Your Understanding
Use the worked analysis to answer each question, then compare with the solutions.
- Why is the raw methylation difference of 20 percentage points not the same as the adjusted difference of 15 percentage points?
- What happens to separate condition and batch estimates if all controls are in batch 1 and all treated samples are in batch 2?
- In the chromatin example, what does the 95% fold-change interval of 1.26 to 3.18 describe?
- Does the chromatin p-value of 0.0035 alone establish a genome-wide discovery? Why or why not?
- If one sample fails a prespecified chromatin-library criterion but passes methylation-array QC, must it be removed from both assay analyses?
Solutions: (1) The groups contain different proportions of batches, and the raw contrast mixes the 0.15 condition difference with the 0.10 batch difference. The additive model estimates the condition contrast at a common batch composition. (2) Condition and batch are perfectly confounded, so their separate effects are not identifiable from these observations. (3) Under the stated negative-binomial model and dispersion assumptions, it is a 95% confidence interval for the treated-to-control signal fold change at that region, not a probability statement about the fixed parameter. (4) No. It is an unadjusted p-value for one region; genome-wide conclusions require multiplicity adjustment and consideration of effect size and diagnostics. (5) No. Exclusion should follow assay-specific evidence; the methylation data may remain usable in a methylation-only analysis, while a matched cross-assay analysis must state which samples it includes.