From Raw Profiles to a Diagnostic Check
Expression analysis is a sequence of linked decisions, not a single correction step. First inspect the raw profiles for sample-level anomalies; then choose a normalization compatible with the measurement platform; then represent known technical factors in a model that also preserves the biological comparison. Finally, examine the adjusted data to see whether the technical pattern has actually diminished without erasing the biological signal.
The sequence matters. A normalization method addresses particular systematic differences in measurement, such as array-wide intensity distributions or RNA-seq library composition. Batch adjustment addresses a different problem: expression differences associated with recorded processing groups. Neither operation guarantees that every technical pattern is gone, and neither should be judged only by whether a plot looks more uniform.
The assumptions behind the chosen procedures must be plausible. Common normalization methods rely, in different ways, on the premise that most genes are not changing in one direction across all samples, or that increases and decreases are sufficiently balanced. A simple additive batch model assumes that the batch effect is adequately represented by a common shift on the chosen expression scale. If biological condition and batch are completely confounded, their separate effects cannot be estimated from the data; no correction algorithm can recover information the design did not collect.
Inspect Before Normalizing
Begin by confirming sample identities, platform, units, and data status: raw array intensities are not interchangeable with background-corrected or log-transformed values, and RNA-seq counts should not be treated as continuous measurements on the original scale. Check for missing values, extreme library sizes or intensities, failed controls, and sample swaps. Plot raw distributions by sample, using boxplots or density curves, and inspect sample relationships with a suitable low-dimensional display. These diagnostics are described in the earlier material on diagnostic plots for expression normalization and principal-component batch detection.
A broad shift in one sample's raw distribution may reflect a technical scale difference, but it is not by itself proof that the sample should be removed. Compare the pattern with recorded processing information and other quality measures. Conversely, overlapping distributions do not prove that samples are free of batch effects: technical structure can be gene-specific or affect only a subset of genes.
For a microarray, inspect intensity distributions and spatial or control-probe diagnostics when available. For RNA-seq, check total reads, mapped reads, library composition, and the prevalence of low-count genes. A sample with an unusual profile may be a legitimate biological extreme, a technical failure, or both; removal requires a documented reason rather than a preference for a tidy plot.
Choose Normalization for the Platform
For many single-channel microarray workflows, quantile normalization makes the empirical distribution of probe intensities the same across arrays. It is appropriate when arrays are intended to be comparable and a common overall distribution is plausible. It can be problematic when a genuine, broad shift in expression distribution is scientifically expected, because forcing distributions to match may remove part of that signal.
RNA-seq begins with integer read counts. A common size-factor approach estimates relative sample scaling from the ratios of counts to gene-wise geometric means. The resulting factors can be used as offsets in a count model or to form normalized counts for descriptive plots. Normalized counts are useful for visualization, but the procedure does not make the counts normally distributed or replace a count model for inference. The earlier material on batch handling in RNA-seq count models treats that inferential distinction in more detail.
Worked Example: Quantile Normalization of Three Arrays
Suppose three hypothetical microarrays have raw log-scale intensities for four probes. The values are already ordered by probe within each array: Array A has \(2,4,6,8\), Array B has \(3,5,7,9\), and Array C has \(4,6,8,10\). Assume these arrays are intended to have comparable overall expression distributions, and that a shared distribution is scientifically reasonable.
At each rank, average the values across arrays. The rank-wise averages are \((2+3+4)/3=3\), \((4+5+6)/3=5\), \((6+7+8)/3=7\), and \((8+9+10)/3=9\). Quantile normalization replaces each value by the average for its rank within that array.
Thus Array A's values \(2,4,6,8\) become \(3,5,7,9\); Array B's values become \(3,5,7,9\); and Array C's values become \(3,5,7,9\). Each array retains its within-array rank ordering, while the three empirical distributions now match. This toy example has no ties; ties require a defined averaging convention for the tied ranks.
The numerical result is not evidence that the arrays are now biologically identical. Only their marginal distributions have been made equal. The method is defensible here because the example assumes comparable arrays and no meaningful global distribution shift; if a global biological shift were plausible, this normalization could distort the target signal.
Worked Example: Median-Ratio Size Factors for RNA-Seq
Consider hypothetical raw counts for three genes in three samples. Suppose most genes are stable in relative expression and the main sample difference is sequencing scale. The count table is:
| Gene | Sample 1 | Sample 2 | Sample 3 |
|---|---|---|---|
| Gene 1 | 100 | 200 | 400 |
| Gene 2 | 50 | 100 | 200 |
| Gene 3 | 25 | 50 | 100 |
For each gene, take the geometric mean across samples. The geometric means are 200 for Gene 1, 100 for Gene 2, and 50 for Gene 3. Divide each count by its gene's geometric mean; the ratios are \(0.5,0.5,0.5\) for Sample 1, \(1,1,1\) for Sample 2, and \(2,2,2\) for Sample 3. The median of each sample's ratios gives size factors \(0.5,1,2\), respectively.
Divide each count by the size factor for its sample. For example, Sample 1 has normalized counts \(100/0.5=200\), \(50/0.5=100\), and \(25/0.5=50\). Sample 2 has \(200/1=200\), \(100/1=100\), and \(50/1=50\). Sample 3 has \(400/2=200\), \(200/2=100\), and \(100/2=50\). Therefore, in every sample the normalized counts are 200 for Gene 1, 100 for Gene 2, and 50 for Gene 3—not 200 for every gene.
This result is consistent with the data: all three genes' raw counts double from one sample to the next, and the size factors remove that common scaling. In real data, the stability assumption can fail if a large fraction of genes changes strongly in one direction or if a small number of highly expressed genes dominate the library. Inspect library composition and normalized profiles; if the assumptions are doubtful, compare an appropriate alternative such as a trimmed-mean normalization rather than treating any one factor estimate as automatic.
Model Known Batch While Preserving Biology
After platform-appropriate normalization, represent known batch and biological condition together. For a gene \(g\) measured on a log expression scale, a basic additive model is:
Here \(x_i\) encodes the biological contrast of interest, \(\gamma_{g,b_i}\) is the batch contribution relative to a chosen reference, and \(\varepsilon_{gi}\) is residual variation. The model requires an identifiable design: the condition and batch columns cannot be perfectly redundant. For ordinary least squares inference, the usual assumptions include an adequate linear mean model and independent, mean-zero errors with appropriate variance assumptions; in this tutorial, adjustment is illustrated as estimation of the fitted batch contribution. If batch effects differ by condition, an additive model without a batch-by-condition interaction may be inadequate.
With a full-rank design matrix \(X\), least squares estimates are \(\hat{\beta}_g=(X^\mathsf{T}X)^{-1}X^\mathsf{T}Y_g\). To adjust expression while retaining the biological terms, subtract only the estimated batch contribution, not the condition contribution or the intercept. A chosen reference-batch parameterization leaves the reference batch unchanged and shifts other batches relative to it. Centered adjustments are also possible, but the convention must be stated because adjusted values can differ by a constant while fitted biological contrasts remain the same.
Worked Example: Adjusting a Balanced Two-Batch Design
Suppose a hypothetical study measures one gene on a normalized log2-expression scale. It has two control and two treated samples in each of two batches, so condition is represented in both batches. Values are: batch 1 controls \(8,9\), batch 1 treated \(10,11\), batch 2 controls \(10,11\), and batch 2 treated \(12,13\). The design is balanced, and an additive model with no batch-by-condition interaction is plausible for this illustration.
Use batch 1 and control as references, so \(x_i=1\) for treated samples and \(I(b_i=2)=1\) for batch 2. In each of the four condition-batch cells there are two observations. The cell means are 8.5, 10.5, 10.5, and 12.5, respectively. Thus the fitted intercept is the batch 1 control mean, \(\hat{\beta}_0=8.5\); the treatment contrast is \(10.5-8.5=2\) log2 units; and the batch 2 effect is \(10.5-8.5=2\) log2 units. The treated mean in batch 2, \(12.5\), is consistent with the additive fit \(8.5+2+2=12.5\).
The adjusted values are batch 1 controls \(8,9\), batch 1 treated \(10,11\), batch 2 controls \(8,9\), and batch 2 treated \(10,11\). The control and treated means remain 8.5 and 10.5 in each batch, so their difference remains 2 log2 units. Each observation's residual from its fitted cell mean is \(-0.5\) or \(0.5\); the residuals average to zero in every cell. In this deliberately small example, that is consistent with the fitted model, but it is not enough to establish that no other technical structure remains.
Because both conditions appear in both batches, the treatment contrast and batch shift are separately estimable. If all controls had been run in batch 1 and all treated samples in batch 2, the observed difference could not distinguish treatment from batch; subtracting a fitted batch effect would then require an untestable assumption about which source caused the difference.
Check the Adjusted Data, Not Just the Adjustment
A complete workflow compares relevant diagnostics before and after adjustment. Plot normalized expression distributions and sample-level ordination, then color or label samples by batch and biological condition. Examine whether samples still separate by batch after accounting for condition, and whether the biological grouping or contrast of interest remains. For a fitted gene-wise model, inspect residuals against fitted values and batch, and look for systematic trends, outliers, and unequal spread. Residual diagnostics assess what the model did not explain; they are not a guarantee that the model is correct.
A useful justification for checking residuals follows directly from the model. If the fitted conditional mean is correct, then the residual \(e_{gi}=Y_{gi}-\widehat{Y}_{gi}\) should not retain a systematic mean pattern by modeled batch. A plot of residual averages by batch can expose a remaining shift. However, centering residuals in a fitted model is partly built into least squares, so this check must be supplemented by design-aware plots and assessment of biology, influential samples, and possible interactions.
Confirm sample labels, platform, count or intensity scale, library or array quality, and any exclusions with documented reasons.
Use a method whose assumptions fit the platform and biological question; retain the raw data and record the normalization choices.
Inspect the condition-by-batch design and fit a model that includes the biological variables that must be preserved.
Compare batch and biology patterns before and after adjustment, inspect residuals, and investigate unexpected separation rather than applying another correction automatically.
The principal alternatives address different circumstances. Quantile normalization is a distribution-matching method commonly used for appropriate array data. RNA-seq size-factor and trimmed-mean approaches address sample-level count composition or scale. Linear covariate adjustment and empirical Bayes procedures model known batches on a suitable expression scale, while methods for latent factors target unrecorded structure under additional assumptions. These approaches are not interchangeable: the data type, design, scientific contrast, and assumptions determine which is reasonable.
| Approach | Typical role | Key caution |
|---|---|---|
| Quantile normalization | Aligns array-wide empirical distributions | Can remove a genuine global shift |
| RNA-seq size factors or trimmed-mean scaling | Accounts for relative library scale or composition | Depends on assumptions about the bulk of genes and composition |
| Known-batch model adjustment | Removes estimated effects of recorded technical groups | Cannot separate batch and biology under complete confounding |
| Latent-factor adjustment | Represents unrecorded technical variation | Factors may capture biology or unwanted variation depending on specification |
Common Failure Modes and Limitations
- Normalizing by appearance alone. Making all distributions look alike is not a universal goal. Choose a method based on platform and assumptions, and check whether the scientific signal could include a global shift.
- Adjusting raw counts as though they were log-scale measurements. RNA-seq count models handle counts and offsets directly. For descriptive transformed values, state the transformation and do not present the result as a substitute for count-based inference.
- Removing batch without preserving condition. A correction that ignores the biological design may erase the very contrast under study. Include estimable biological terms in the model and verify that the contrast remains.
- Interpreting a clean PCA plot as proof of success. A plot summarizes only some variation and may miss gene-specific residual batch effects. Use complementary diagnostics and report limitations.
- Overcorrecting or correcting confounded data. A method cannot identify the source of a perfectly confounded difference from expression values alone. In non-confounded data, aggressive adjustment can still suppress real biology if nuisance and biological signals are not distinguished.
No workflow can guarantee removal of every unmeasured technical effect. Normalization methods also have finite-sample behavior, can be sensitive to outliers or composition, and may behave poorly when their assumptions are violated. Treat adjusted values as outputs of an explicit model, retain an audit trail of raw data and choices, and carry appropriate uncertainty into downstream analysis rather than treating corrected expression as error-free.
Check Your Understanding
Answer each question, then compare with the brief solutions below.
- In the RNA-seq example, what is the normalized count for Gene 2 in Sample 3, and what does it represent?
- Why is quantile normalization not automatically suitable when a broad biological shift across arrays is expected?
- In the balanced batch example, what is the fitted treatment contrast, and why can it be estimated separately from batch?
- Suppose all controls were processed in batch 1 and all treated samples in batch 2. Can the treatment and batch effects be uniquely separated by the additive model?
Solutions. (1) The normalized count is \(200/2=100\). It is Gene 2's count after dividing Sample 3 by its size factor; it is not a raw count and does not itself provide a count-model test. (2) Quantile normalization forces the arrays to share an empirical distribution, which could remove a genuine global biological shift. (3) The fitted contrast is \(10.5-8.5=2\) log2-expression units. Both conditions occur in both batches, so the design contains information separating the two effects. (4) No. Condition and batch are perfectly confounded, so their effects are not separately identifiable from these data.