Normalization: making arrays comparable

The next step in methylation data pre-processing is normalization. There is variation in intensities from the different probe chemistries, as well as from technical variation in the assay, which needs to be accounted for before the probe intensities can be compared across samples in downstream analyses. Normalization is the method used to do this. This notebook walks through the standard normalization strategies used for methylation data.

Since all of the samples from our small example subset passed QC, we will start normaliztion from the raw two-channel object produced in QC. If you did identify samples that failed to meet QC thresholds, those samples should be removed prior to normalization.

What this chapter needs

This chapter starts from data/01_RGset.rds, written by chapter 00. If you did not run that chapter, fetch the published checkpoint instead:

Terminal

./get_data.sh B_qc

Code

library(minfi)                                        # preprocessRaw/Noob/Funnorm(), getBeta()
library(IlluminaHumanMethylationEPICanno.ilm10b4.hg19) # EPIC probe annotation: design type, coordinates
library(ggplot2)                                       # every figure in this chapter

# The RGChannelSet written by the Setup chapter: raw red and green intensities
# for 1,051,943 addresses across 96 arrays, with no normalization applied. This
# is the same object the QC chapter worked from -- all 96 samples passed, so
# nothing has been dropped from it.
rg <- readRDS("data/01_RGset.rds")

rg

class: RGChannelSet 
dim: 1051943 96 
metadata(0):
assays(2): Green Red
rownames(1051943): 1600101 1600111 ... 99810990 99810992
rowData names(0):
colnames(96): GSM3853168 GSM3853169 ... GSM3853880 GSM3853881
colData names(15): sample_id series ... childhood_abuse sex
Annotation
  array: IlluminaHumanMethylationEPIC
  annotation: ilm10b4.hg19

1. Two sources of unwanted variation

Within-array, technical. Every Infinium measurement is a fluorescence intensity, and fluorescence has a background floor (autofluorescence, non-specific hybridization, optical noise) and a dye bias: the Cy3 (green) and Cy5 (red) channels are not excited or detected with equal efficiency, so the same amount of methylated and unmethylated DNA does not produce the same raw signal. Left uncorrected, this shifts β-values systematically.

Within-array, by probe design. The EPIC array (like the 450K before it) mixes two assay chemistries (Pidsley et al. 2016):

Type I probes use two bead types per CpG (one for methylated, one for unmethylated), both read in the same color channel. They tolerate a wider dynamic range but use up two bead addresses per site.
Type II probes use a single bead per CpG, with methylated and unmethylated states read in different channels. They are twice as space-efficient — which is why 84% of EPIC probes are Type II (723,722 of 865,859; the comparable 450K figure is 72%) — but their β-value distribution is compressed toward the center: they under-call fully methylated and fully unmethylated sites relative to Type I.

We can see this compression directly in the subset. Below are the median β and the between-probe spread (SD of per-probe mean β) for each design, on raw data:

Code

# Which design each CpG uses comes from the EPIC annotation package, not from the
# data: the `Type` column is "I" or "II" for each of the 865,859 annotated probes.
# These are the probes that carry hg19 coordinates, which is also what the genome
# mapping inside preprocessFunnorm needs -- so sampling from this set guarantees
# the same probes exist in all three normalized matrices built below.
ann    <- getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
typeI  <- rownames(ann)[ann$Type == "I"]
typeII <- rownames(ann)[ann$Type == "II"]

# preprocessRaw() only re-arranges the two raw channels into methylated and
# unmethylated signal -- no correction of any kind is applied, which is exactly
# what we want as the "before" picture. It takes a few seconds.
mset_raw     <- preprocessRaw(rg)
beta_raw_all <- getBeta(mset_raw)

# The figures and tables in this chapter describe the shape of the β
# distribution, and 20,000 probes describe that shape as well as 865,859 do at a
# thirtieth of the memory (a full matrix is ~665 MB; three of them will not fit
# comfortably alongside normalization). The seed is fixed so the sample -- and
# therefore every number in the tables below -- is the same on every render.
set.seed(1)
probe_sample <- sample(intersect(rownames(ann), rownames(beta_raw_all)), 20000)

beta_raw <- beta_raw_all[probe_sample, ]

# Drop the full-size objects as soon as the subset exists; the peak memory of
# this chapter is set by the normalization step below, not by this one.
rm(beta_raw_all, mset_raw)
invisible(gc())

dim(beta_raw)

[1] 20000    96

Code

# One summary row per probe design. `spread_sd` is the SD of the per-probe mean β
# across samples, i.e. how much the probes of a design differ from one another --
# not within-probe noise, which is why rowMeans() comes first.
type_summ <- function(beta, typeI, typeII) {
  bi  <- beta[rownames(beta) %in% typeI, , drop = FALSE]
  bii <- beta[rownames(beta) %in% typeII, , drop = FALSE]
  data.frame(
    design      = c("Type I", "Type II"),
    n_probes    = c(nrow(bi), nrow(bii)),
    median_beta = round(c(median(bi, na.rm = TRUE), median(bii, na.rm = TRUE)), 3),
    spread_sd   = round(c(sd(rowMeans(bi, na.rm = TRUE)),
                          sd(rowMeans(bii, na.rm = TRUE))), 3)
  )
}
knitr::kable(type_summ(beta_raw, typeI, typeII),
             caption = "Raw β by probe design (20k-probe sample). Type II sits closer to 0.5 with a narrower spread — the design-bias normalization must correct.")

Raw β by probe design (20k-probe sample). Type II sits closer to 0.5 with a narrower spread — the design-bias normalization must correct.
design	n_probes	median_beta	spread_sd
Type I	3184	0.092	0.373
Type II	16816	0.745	0.286

Why this matters for an EWAS

If Type I and Type II probes are on different scales, a differential-methylation test will partly reflect which chemistry a probe uses rather than biology, and effect sizes are not comparable across probe types.

2. Normalization methods

Methods fall into two groups by what they try to make comparable.

Within-array methods (fix probe-design bias per sample)

noob (normal-exponential out-of-band background correction) models the background as a normal-exponential mixture using the out-of-band Type I signals, then applies a dye-bias correction (Triche et al. 2013). It is a background/dye step, not a design-bias step, but it is the standard first move and is what most pipelines build on.
BMIQ (beta-mixture quantile) fits a three-state beta mixture (unmethylated / hemimethylated / methylated) separately to Type I and Type II probes, then transforms the Type II distribution onto the Type I scale (Teschendorff et al. 2013). It directly targets the design bias.
SWAN (subset-quantile within-array normalization) matches Type I and Type II intensity distributions using a subset of probes with matched underlying biology (Maksimovic et al. 2012).

Between-array methods (make samples comparable to each other)

Quantile normalization forces every sample’s intensity distribution to a common reference. Powerful, but it assumes the global methylation distribution is the same across samples — a poor assumption when groups differ biologically (e.g. tumor vs normal), where it can erase real signal.
Functional normalization (funnorm) is a between-array method designed for exactly that case (Fortin et al. 2014). Instead of forcing distributions to match, it regresses out variation explained by the array’s internal control probes — bisulfite conversion, staining, extension, hybridization controls — which capture technical variation without being tied to biological signal. It is the recommended default when you expect large biological differences between groups, and it wraps noob as its first step.

Choosing a method

There is no single winner. A defensible default for blood EWAS is noob + functional normalization (what preprocessFunnorm does in one call), because it corrects background/dye and removes technical variation without assuming samples share a methylation distribution. If your study is a within-tissue comparison with small expected differences, quantile-based methods (or preprocessQuantile) are also reasonable.

3. Running it on the subset

We compute three versions and compare: raw (preprocessRaw, no correction), noob (preprocessNoob), and functional (preprocessFunnorm). The functional result is a GenomicRatioSet — CpGs mapped to the genome with β and M-values ready for analysis.

The code to produce these subsets is shown below. It may take several minutes to complete.

Code

# noob on its own: background correction from the out-of-band Type I signals plus
# a dye-bias correction, one array at a time. Roughly a minute and several GB for
# 96 arrays. Subset to the same 20,000 probes immediately and drop the full
# MethylSet, so the funnorm step below starts from a clean slate.
mset_noob <- preprocessNoob(rg)
beta_noob <- getBeta(mset_noob)[probe_sample, ]

rm(mset_noob)
invisible(gc())

Code

# This is the expensive step in the chapter: about 2 minutes and a peak of
# ~12.5 GB at 96 arrays, because funnorm runs noob first and then fits the
# control-probe PCA across the full 865,859-probe matrix. If your machine cannot
# hold that, take the resume path in the callout below.
grs_fun <- preprocessFunnorm(rg)

# `grs_fun` is the handoff object for the rest of the tutorial -- keep it whole.
# Only the β matrix gets subset, for the figures in this chapter.
beta_fun <- getBeta(grs_fun)[probe_sample, ]

# A design comparison is only meaningful if all three matrices describe the same
# probes; fail loudly here rather than silently comparing different probe sets.
stopifnot(identical(rownames(beta_raw), rownames(beta_noob)),
          identical(rownames(beta_raw), rownames(beta_fun)))

grs_fun

class: GenomicRatioSet 
dim: 865859 96 
metadata(0):
assays(2): Beta CN
rownames(865859): cg14817997 cg26928153 ... cg07587934 cg16855331
rowData names(0):
colnames(96): GSM3853168 GSM3853169 ... GSM3853880 GSM3853881
colData names(18): sample_id series ... yMed predictedSex
Annotation
  array: IlluminaHumanMethylationEPIC
  annotation: ilm10b4.hg19
Preprocessing
  Method: NA
  minfi version: NA
  Manifest version: NA

If normalization does not fit on your machine

The two chunks above need the 460 MB RGChannelSet in memory and peak around 12.5 GB. If that is not practical, the published checkpoint holds exactly the pieces the figures below use — the three 20,000-probe β matrices and the two design lists — so you can fetch it and carry on:

Terminal

./get_data.sh C_normalized

RStudio Console

np <- readRDS("data/02_norm_pieces.rds")
beta_raw  <- np$beta_raw
beta_noob <- np$beta_noob
beta_fun  <- np$beta_fun
typeI     <- np$typeI
typeII    <- np$typeII

Every figure and table in this chapter will then run. What you give up is grs_fun itself, which probe filtering reads — for that chapter, fetch ./get_data.sh D_filtered instead. You also give up watching the normalization run, so take this route because you have to, not to save two minutes.

Then we can compare the distribution of signals across the three versions of the data and the probe types.

Code

# Reshape one β matrix into the long form ggplot needs. 8,000 probes x 96 samples
# is already ~768,000 density points per method, which is plenty -- the extra
# thinning keeps the three panels quick to draw.
mk_long <- function(beta, label, typeII) {
  b <- beta[sample(nrow(beta), min(8000, nrow(beta))), ]
  design <- ifelse(rownames(b) %in% typeII, "Type II", "Type I")
  # as.vector() unstacks a matrix column by column, so the probe labels repeat
  # once per sample -- rep(..., times = ncol(b)) keeps them aligned to `beta`.
  data.frame(beta = as.vector(b),
             design = rep(design, times = ncol(b)),
             method = label)
}
long <- rbind(
  mk_long(beta_raw,  "raw",     typeII),
  mk_long(beta_noob, "noob",    typeII),
  mk_long(beta_fun,  "funnorm", typeII)
)
# Panels in processing order, not alphabetical: raw is the "before" picture.
long$method <- factor(long$method, levels = c("raw", "noob", "funnorm"))
ggplot(long, aes(x = beta, color = design, linetype = design)) +
  geom_density(na.rm = TRUE, linewidth = 0.7) +
  facet_wrap(~method) +
  scale_color_manual(values = c("Type I"  = "#1A6B75",     # teal
                                "Type II" = "#8C3A4A")) +  # plum
  labs(x = "β-value", y = "density",
       title = "Probe-design bias, before and after normalization") +
  theme(legend.position = "bottom")

β-value density by probe design, before and after normalization (20k-probe sample). Raw Type II (dashed) is visibly pulled toward the center; noob and funnorm restore a more bimodal shape and bring the two designs closer together.

Code

# The same two-row summary as above, stacked across the three methods, so the
# Type I / Type II gap can be read down the column.
tab <- rbind(
  cbind(method = "raw",     type_summ(beta_raw,  typeI, typeII)),
  cbind(method = "noob",    type_summ(beta_noob, typeI, typeII)),
  cbind(method = "funnorm", type_summ(beta_fun,  typeI, typeII))
)
knitr::kable(tab, row.names = FALSE,
             caption = "Type I vs Type II β summaries across methods. The Type II spread widens toward the Type I spread after normalization. The medians do not converge, and are not expected to — see the text.")

Type I vs Type II β summaries across methods. The Type II spread widens toward the Type I spread after normalization. The medians do not converge, and are not expected to — see the text.
method	design	n_probes	median_beta	spread_sd
raw	Type I	3184	0.092	0.373
raw	Type II	16816	0.745	0.286
noob	Type I	3184	0.059	0.391
noob	Type II	16816	0.807	0.318
funnorm	Type I	3184	0.058	0.386
funnorm	Type II	16816	0.790	0.313

Read the spread_sd column, not the medians. Type II’s spread rises from 0.286 (raw) to 0.313 after functional normalization, closing about a sixth of the gap to the Type I spread of 0.385. That expansion of the compressed Type II range is what design-bias correction does, and it is the same thing the density curves above show.

The medians move the other way: Type II sits at 0.745 raw and 0.790 after correction, so the median gap between the two designs widens slightly rather than closing. That is not a failure of the correction, and it is worth understanding why. Type I probes are enriched in CpG islands, which are predominantly unmethylated, while Type II probes cover proportionally more open sea. The two rows are therefore summarising different populations of CpGs, and the distance between their medians is a fact about where each design sits in the genome rather than a measure of design bias. A statistic that mixes the two cannot tell you whether normalization worked.

Functional normalization additionally removes between-array technical variance that these density plots do not show; its payoff appears later, in the batch-effects notebook, as reduced technical structure in the principal components.

4. β-values vs M-values

Normalization outputs β-values (proportion methylated, bounded 0–1), which are interpretable but heteroscedastic — their variance is small near 0 and 1 and large near 0.5, which violates the constant-variance assumption of linear models. For statistical testing we use the M-value, the log₂ ratio of methylated to unmethylated signal (Du et al. 2010):

\[ M = \log_2\!\left(\frac{\beta}{1-\beta}\right) \]

M-values are roughly homoscedastic and better behaved for the linear models used in an EWAS. The convention is to test on M-values and report/interpret effect sizes on β-values (a β difference of 0.05 is a “5% methylation change,” which is biologically legible; an M-value difference is not). It is becoming more common to report/interpret from M-values instead of beta values, so here we only show analysis and reporting from M-values, but to report effect sizes from Beta values you would follow the same steps, just with the Beta value matrix instead of the M-value matrix.

Code

# The transform itself, drawn over the open interval: β = 0 and β = 1 map to
# -Inf and +Inf, which is the practical reason M-values need a bounded β.
b <- seq(0.001, 0.999, length.out = 500)
m <- log2(b / (1 - b))
dfm <- data.frame(beta = b, M = m)
ggplot(dfm, aes(beta, M)) +
  geom_line(color = "#1A6B75", linewidth = 0.9) +        # teal
  geom_hline(yintercept = 0, color = "gray70") +
  geom_vline(xintercept = 0.5, color = "gray70") +
  labs(x = "β-value", y = "M-value", title = "β → M-value transform")

The β→M transform. β is bounded and heteroscedastic; M is unbounded and closer to constant-variance.

5. Normalizing from a processed matrix

Many GEO series do not deposit raw IDATs, but do provide a processed matrix: a table of β-values (or sometimes M-values) with CpGs in rows and samples in columns, often as a GSExxxxx_series_matrix.txt.gz or a supplementary file. You cannot run minfi’s IDAT-based normalization on these, because the raw two-channel intensities are gone. Here is how to work with them.

Step 1 — get the matrix into R. For a series matrix:

RStudio Console

library(GEOquery)
gse  <- getGEO("GSExxxxxx", GSEMatrix = TRUE)      # returns a list of ExpressionSets
eset <- gse[[1]]
beta <- Biobase::exprs(eset)                        # CpG x sample β matrix
pheno <- Biobase::pData(eset)                       # sample metadata

For a supplementary flat file (getGEOSuppFiles downloads it), read it with data.table::fread() and coerce to a matrix with CpG IDs as rownames.

Step 2 — find out what was already done. Read the series’ processing description (experimentData() / pData(eset)$data_processing) to learn whether the matrix is already normalized, and how. Common cases:

Already noob + BMIQ (or funnorm) normalized — the most common. Do not re-normalize; treat the matrix as analysis-ready β and proceed to probe filtering.
Raw/unnormalized β — you can apply a β-space normalization such as BMIQ (wateRmelon::BMIQ()), which only needs β and probe design, not intensities.
Already batch-corrected (e.g. ComBat) or cell-composition-adjusted — be cautious: further correction can over-adjust, and you have lost the ability to model those effects yourself. Note it as a limitation.

Step 3 — QC on what you have. Without intensities you lose detection p-values and the intensity/sex-signal checks from the QC notebook. You can still:

check for missing values and impute or drop probes with excessive NA;
run PCA/MDS to spot outlier samples and technical structure;
predict sex from the X/Y-chromosome β distributions (or wateRmelon’s estimateSex) and compare to the reported sex, catching label swaps;
verify the probe count and platform match what you expect (EPIC ≈ 865k).

Step 4 — harmonize probe IDs. Processed matrices vary: some are already masked/filtered, some use different ID conventions, some are 450K and some EPIC. Before combining or comparing series, intersect on common cg/ch probe IDs and apply the same probe-filtering (next notebook) you would apply to IDAT-derived data.

Practical rule

From IDATs you control normalization; from a processed matrix you inherit it. The job shifts from “normalize” to “audit what was done and decide whether it is adequate,” then continue from probe filtering. Always record, in your methods, the processing state of every series you use — reviewers will ask.

We save the functional-normalized GenomicRatioSet for the downstream steps:

Code

# 1. The handoff object: 865,859 CpGs x 96 samples, genome-mapped, noob +
#    functional normalized. This is what probe filtering (chapter 03) reads, and
#    it is the file deposited as Zenodo tier C_normalized. ~1 GB on disk.
saveRDS(grs_fun, "data/02_funnorm_grs.rds")

# 2. The pieces the figures above are drawn from: the three 20,000-probe β
#    matrices and the two design lists. Saving them means the resume callout in
#    section 3 restores exactly the objects this chapter computed, and that a
#    reader who cannot run normalization still gets every figure. Small next to
#    the GenomicRatioSet (~46 MB), and an RDS rather than a flat file because it
#    is three matrices and two character vectors, not one rectangle.
norm_pieces <- list(
  beta_raw  = beta_raw,    # 20,000 x 96, no correction
  beta_noob = beta_noob,   # 20,000 x 96, noob only
  beta_fun  = beta_fun,    # 20,000 x 96, noob + functional normalization
  typeI     = typeI,       # Type I probe IDs, whole array
  typeII    = typeII       # Type II probe IDs, whole array
)
saveRDS(norm_pieces, "data/02_norm_pieces.rds")

Next: Probe filtering — removing probes that are unreliable by design (cross-reactive, SNP-affected, sex-chromosome), using Zhou’s curated masks with ancestry awareness.

References

Du, Pan, Xiao Zhang, Chiang-Ching Huang, et al. 2010. “Comparison of Beta-Value and m-Value Methods for Quantifying Methylation Levels by Microarray Analysis.” BMC Bioinformatics 11: 587. https://doi.org/10.1186/1471-2105-11-587.

Fortin, Jean-Philippe, Aurélie Labbe, Mathieu Lemire, et al. 2014. “Functional Normalization of 450k Methylation Array Data Improves Replication in Large Cancer Studies.” Genome Biology 15 (12): 503. https://doi.org/10.1186/s13059-014-0503-2.

Maksimovic, Jovana, Lavinia Gordon, and Alicia Oshlack. 2012. “SWAN: Subset-Quantile Within Array Normalization for Illumina Infinium HumanMethylation450 BeadChips.” Genome Biology 13 (6): R44. https://doi.org/10.1186/gb-2012-13-6-r44.

Pidsley, Ruth, Elena Zotenko, Timothy J Peters, et al. 2016. “Critical Evaluation of the Illumina MethylationEPIC BeadChip Microarray for Whole-Genome DNA Methylation Profiling.” Genome Biology 17 (1): 208. https://doi.org/10.1186/s13059-016-1066-1.

Teschendorff, Andrew E, Francesco Marabita, Matthias Lechner, et al. 2013. “A Beta-Mixture Quantile Normalization Method for Correcting Probe Design Bias in Illumina Infinium 450k DNA Methylation Data.” Bioinformatics 29 (2): 189–96. https://doi.org/10.1093/bioinformatics/bts680.

Triche, Timothy J, Daniel J Weisenberger, David Van Den Berg, Peter W Laird, and Kimberly D Siegmund. 2013. “Low-Level Processing of Illumina Infinium DNA Methylation BeadArrays.” Nucleic Acids Research 41 (7): e90. https://doi.org/10.1093/nar/gkt090.