Cell composition: Whole Blood

Whole blood is not one thing. It is a mixture of neutrophils, monocytes, natural killer (NK) cells, B cells, and CD4+/CD8+ T cells, and each cell type has a distinctly different methylome. When you draw blood from two people, differences in their cell-type proportions — driven by age, infection, stress, time of day — can produce large methylation differences that are unrelated to the biology you’re studying. For most blood EWAS, cell composition is a large source of confounding (Houseman et al. 2012; Jaffe and Irizarry 2014). This notebook estimates the proportions, explains the reference-based method, and validates our estimates against the proportions the study deposited.

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)                  # estimateCellCounts2(), getBeta()
library(FlowSorted.Blood.EPIC)  # the flow-sorted EPIC reference and the IDOL probe library
library(GEOquery)               # getGEO(), for the study's own deposited proportions
library(data.table)             # fread() for the deposited table, rbindlist() for the figure
library(ggplot2)                # every figure in this chapter

# The raw two-channel object written by the Setup chapter: 1,051,943 addresses
# across 96 arrays, with no normalization applied. Deconvolution is one of the
# few steps that needs the raw intensities rather than the normalized object
# from chapter 02 -- the reason is in the callout below.
rg <- readRDS("data/01_RGset.rds")

1. Why cell composition matters

A neutrophil and a CD4+ T cell carry the same genome but very different DNA methylation. A whole-blood methylation measurement at any CpG is therefore a weighted average across the cell types present, with weights equal to the cell proportions:

\[\beta_{\text{observed}} \approx \sum_{k} w_k \, \beta_k\]

where $w_k$ is the fraction of cell type $k$ and $\beta_k$ is that cell type’s methylation at the CpG. If cases have more neutrophils than controls — even for reasons unrelated to the exposure — every cell-type-discriminating CpG will look “differentially methylated.” Neutrophil fraction alone routinely explains the largest chunk of genome-wide methylation variance in blood.

2. Reference-based deconvolution (Houseman)

The standard approach is the Houseman method (Houseman et al. 2012):

Start from a reference of purified, flow-sorted blood cell types, each profiled on the same array platform. For each cell type we know its methylation at a panel of cell-type-discriminating CpGs.
Choose the discriminating CpGs. The original method picked them by F-test; the IDOL optimization (Salas et al. 2018) selects a library that minimizes deconvolution error and is the current best practice for EPIC. We use the IDOL-optimized 450-CpG library from FlowSorted.Blood.EPIC.
For each sample, solve a constrained regression (proportions ≥ 0, summing to ~1) of the observed methylation at those CpGs against the reference profiles. The fitted coefficients are the estimated cell proportions.

Crucially, the sample and reference data must be normalized together so that technical differences don’t confound composition. minfi’s estimateCellCounts2 handles this by combining RGChannelSet with the reference RGChannelSet, normalizing jointly (here with noob), then deconvolving:

Code

# The IDOL-optimized probe library ships with FlowSorted.Blood.EPIC as a dataset
# rather than a function, so it has to be attached before it can be passed in.
data("IDOLOptimizedCpGs")

# The published proportions were produced with this seed set; keeping it means a
# re-render reproduces the same numbers rather than something merely similar.
set.seed(1)

# This is the most memory-hungry step in the whole tutorial: roughly 1 min 54 s
# and a 15.0 GB peak for 96 arrays. The cost is structural, not wasteful -- the
# 96 arrays and the flow-sorted reference are stacked into one object and noob-
# normalized together (see the callout below), so the peak holds both datasets
# at full probe resolution at once. An 8 GB laptop cannot run this chunk; the
# collapsed callout further down is the route around it.
cc <- estimateCellCounts2(
  rg,                                   # raw RGChannelSet (not yet normalized)
  compositeCellType   = "Blood",
  processMethod       = "preprocessNoob",
  probeSelect         = "IDOL",
  cellTypes           = c("CD8T","CD4T","NK","Bcell","Mono","Neu"),
  referencePlatform   = "IlluminaHumanMethylationEPIC",
  IDOLOptimizedCpGs   = IDOLOptimizedCpGs,
  returnAll           = FALSE)

# `prop` is the estimated proportion per cell type; the sibling `counts` element
# is only filled in when you supply absolute cell counts to compare against, so
# on this call it comes back empty. Read `prop`.
props <- cc$prop                        # 96 samples x 6 cell types
dim(props)

[1] 96  6

Reference-based needs the raw data

estimateCellCounts2 starts from the raw RGChannelSet, not the normalized object from normalization, because it must normalize the samples and reference. This is one of the few steps that reaches back to the raw IDATs. The reference dataset itself is downloaded once from Bioconductor’s ExperimentHub and cached.

Running out of memory on the chunk above

A 15.0 GB peak is more than a 16 GB laptop has spare, and far more than an 8 GB one has at all. The proportions this chapter writes at the end are committed to the repository, so you can skip the deconvolve chunk and read them instead — every section below it then runs:

RStudio Console

props <- readRDS(url(paste0("https://raw.githubusercontent.com/krferrier/",
                            "Methylation-EWAS-tutorial/main/tutorial/data/",
                            "04_cellcounts.rds")))

What this route costs you is the one thing this chapter is about: you do not see the reference being fetched, you do not see the joint normalization, and if your own study ever needs different cell types or a different reference, you will not have run the call that changes. It also leaves cc undefined, so the save of the full deconvolution return at the end of the chapter cannot run. Use it because your machine forces you to.

3. The composition of this subset

Code

# Collapse the 96 x 6 matrix to one row per cell type: the mean proportion and
# the observed range across samples. Percentages read more easily than
# proportions in a table, hence the 100x.
comp <- data.frame(cell_type = colnames(props),
                   mean_pct   = round(100 * colMeans(props), 1),
                   min_pct    = round(100 * apply(props, 2, min), 1),
                   max_pct    = round(100 * apply(props, 2, max), 1))
comp <- comp[order(-comp$mean_pct), ]
knitr::kable(comp, row.names = FALSE,
             caption = "Estimated cell-type composition across the 96 samples. Neutrophils dominate, as expected for whole blood, but the range is wide — this between-sample variation is exactly what must be modeled.")

Estimated cell-type composition across the 96 samples. Neutrophils dominate, as expected for whole blood, but the range is wide — this between-sample variation is exactly what must be modeled.
cell_type	mean_pct	min_pct	max_pct
Neu	53.6	17.4	83.7
CD4T	16.4	2.9	32.3
CD8T	11.7	2.6	50.2
Mono	8.5	2.8	27.0
Bcell	7.0	1.5	13.8
NK	5.2	1.7	13.7

Code

# reshape() wants a data.frame with an id column, so the matrix goes wide-to-long
# by hand here: one row per sample x cell type, which is what geom_col stacks.
pl <- as.data.frame(props); pl$sample_id <- rownames(pl)
plm <- reshape(pl, varying = colnames(props), v.names = "prop",
               timevar = "cell", times = colnames(props),
               direction = "long")
# Stacking order is the factor order, so set it explicitly: largest fraction at
# the bottom of the bar, rarest at the top.
plm$cell <- factor(plm$cell, levels = c("Neu","CD4T","CD8T","Mono","Bcell","NK"))
# Order the bars by neutrophil fraction so the gradient across samples is
# visible instead of being scattered by array order.
ord <- order(props[,"Neu"]); plm$sample_id <- factor(plm$sample_id,
                                                     levels = rownames(props)[ord])
ggplot(plm, aes(sample_id, prop, fill = cell)) +
  geom_col(width = 1) +
  scale_fill_brewer(palette = "Set2") +
  labs(x = "sample (ordered by neutrophil fraction)", y = "estimated proportion",
       fill = "cell type", title = "Blood cell composition varies widely between samples") +
  theme_bw(base_size = 11) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

Estimated cell-type proportions per sample (stacked). The neutrophil fraction (largest band) varies from 17% to 84% across individuals — a swing that would swamp most disease effects if left unmodeled.

4. Validation against the deposited answer key

The Grady Trauma Project deposited its own cell-type proportion estimates with the GEO record. That gives us the opportunity to check our deconvolution.

Code

# The deposited proportions are part of the GEO record's *phenotype*, not of the
# IDATs, so they come out of the series matrix rather than out of `rg`. GEOquery
# parses each "label: value" characteristic into its own column named
# "<label>:ch1", which is where the six proportions arrive.
dep_file <- "data/04_deposited_cellcounts.tsv"

if (file.exists(dep_file)) {
  # The series matrix is ~30 MB and the deposited values never change, so once
  # the small table below exists it is what later renders read. Delete it to
  # force the fetch again.
  dep <- as.data.frame(fread(dep_file))
} else {
  # 60 seconds is not enough to pull 30 MB from GEO on a slow link.
  options(timeout = max(1200, getOption("timeout")))

  # getGPL = FALSE skips the platform manifest, which we do not need here.
  gse <- getGEO("GSE132203", GSEMatrix = TRUE, getGPL = FALSE)[[1]]
  pd  <- Biobase::pData(gse)

  # Our cell-type names on the left, GEO's column names on the right.
  dep_cols <- c(CD8T = "cd8t:ch1", CD4T = "cd4t:ch1", NK = "nk:ch1",
                Bcell = "bcell:ch1", Mono = "mono:ch1", Neu = "neu:ch1")

  # Everything out of a series matrix is text, including the numbers.
  dep <- as.data.frame(lapply(dep_cols,
                              function(k) as.numeric(as.character(pd[[k]]))))
  dep$sample_id <- as.character(pd$geo_accession)

  # The series carries 795 samples; match() keeps the 96 on our plate subset,
  # in the row order of `props` so the merge below cannot silently misalign.
  dep <- dep[match(rownames(props), dep$sample_id), c("sample_id", names(dep_cols))]
}

Code

cells <- c("CD8T","CD4T","NK","Bcell","Mono","Neu")

# One row per sample with both sets of numbers side by side. Both frames carry
# the same six column names, so merge() appends the suffixes that tell them
# apart -- .est for ours, .dep for the study's.
est <- as.data.frame(props); est$sample_id <- rownames(props)
m   <- merge(est, dep, by = "sample_id", suffixes = c(".est", ".dep"))

# Three complementary summaries per cell type: correlation says whether the
# between-sample variation agrees (the part that matters for confounding
# adjustment), mean absolute error says how far apart the absolute values are,
# and the two means expose a systematic offset the correlation would hide.
res <- data.frame(cell = cells, r = NA_real_, mae = NA_real_,
                  mean_est = NA_real_, mean_dep = NA_real_)
for (i in seq_along(cells)) {
  e <- m[[paste0(cells[i], ".est")]]
  d <- m[[paste0(cells[i], ".dep")]]
  res$r[i]        <- cor(e, d, use = "complete.obs")
  res$mae[i]      <- mean(abs(e - d), na.rm = TRUE)
  res$mean_est[i] <- mean(e, na.rm = TRUE)
  res$mean_dep[i] <- mean(d, na.rm = TRUE)
}

# The per-sample comparison and the per-cell-type summary travel together: this
# is the object saved at the end of the chapter, and the one the paragraph below
# the figure reads its numbers from.
val <- list(merged = m, validation = res)

Code

disp <- data.frame(
  `cell type`        = res$cell,
  `Pearson r`        = sprintf("%.3f", res$r),
  `mean abs. error`  = sprintf("%.4f", res$mae),
  `our mean`         = sprintf("%.1f%%", 100*res$mean_est),
  `deposited mean`   = sprintf("%.1f%%", 100*res$mean_dep),
  check.names = FALSE)
knitr::kable(disp, caption = "Our estimates vs the study's deposited proportions, per cell type (n = 96).")

Our estimates vs the study’s deposited proportions, per cell type (n = 96).
cell type	Pearson r	mean abs. error	our mean	deposited mean
CD8T	0.993	0.0066	11.7%	11.8%
CD4T	0.993	0.0100	16.4%	15.5%
NK	0.953	0.0095	5.2%	4.2%
Bcell	0.963	0.0140	7.0%	5.6%
Mono	0.989	0.0044	8.5%	8.3%
Neu	0.997	0.0112	53.6%	54.5%

Code

# Six panels, one per cell type, so stack the pairs into a long table.
long <- rbindlist(lapply(cells, function(k)
  data.table(cell = k, est = m[[paste0(k, ".est")]], dep = m[[paste0(k, ".dep")]])))
long$cell <- factor(long$cell, levels = cells)

# The correlation from the table above, printed in each panel's corner: with
# free axes per panel a reader cannot judge agreement by eye alone.
labs_r <- data.frame(cell = factor(res$cell, levels = cells),
                     lab  = sprintf("r = %.3f", res$r))

p_val <- ggplot(long, aes(dep, est)) +
  # y = x first, so the points sit on top of it.
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "#6E675B") +  # gray
  geom_point(alpha = .6, size = 1.3, color = "#1A6B75") +                          # teal
  facet_wrap(~cell, scales = "free") +
  geom_text(data = labs_r, aes(x = -Inf, y = Inf, label = lab),
            hjust = -0.15, vjust = 1.5, size = 3, inherit.aes = FALSE) +
  labs(x = "deposited proportion (study answer key)", y = "our Houseman/IDOL estimate",
       title = "Deconvolution validation: estimated vs deposited cell proportions")
p_val

Estimated (y) vs deposited (x) proportions for each cell type; dashed line is perfect agreement. Every cell type correlates > 0.95 with the answer key, with mean absolute error under 1.5 percentage points — this pipeline reproduces the study’s own estimates.

Correlations range from 0.953 (NK, the rarest and hardest-to-estimate type) to 0.997 (neutrophils), and the mean absolute error never exceeds 1.4 percentage points. Small systematic offsets (e.g. our B-cell estimates run slightly high) are expected when the reference library or normalization differs from the original authors’, and they matter far less than the near-perfect ranking — because in the EWAS model it is the between-sample variation, not the absolute level, that does the confounding adjustment.

Practical note

Many studies do not provide cell proportions. In a real analysis you trust the method and sanity-check the output: neutrophils should dominate whole blood (~50–70%), proportions should sum to ~1, and the distribution should look biologically plausible.

5. When you don’t have IDATs: deconvolution from a processed matrix

Reference-based deconvolution normally wants the raw RGChannelSet, so a β-matrix-only deposit is the harder case. Two routes:

projectCellType on a β-matrix. The Houseman projection itself only needs methylation values at the reference CpGs, not intensities. Pull the reference β-values and its cell-type means, subset your matrix to the shared reference CpGs, and run the constrained projection directly:

library(FlowSorted.Blood.EPIC)
library(genefilter)
# `beta` is a CpG x sample matrix -- in the deposit-only case it is whatever
# the submitter posted. Locally, the filtered object from chapter 03 gives one:
beta <- minfi::getBeta(readRDS("data/03_grs_filtered.rds"))
# reference cell-type means at the IDOL CpGs (from the package)
ref <- ...   # matrix: IDOL CpGs x 6 cell types
shared <- intersect(rownames(beta), rownames(ref))
props  <- minfi:::projectCellType(beta[shared, ], ref[shared, ])

The estimates are slightly noisier than the joint-normalization route because the sample and reference weren’t normalized together, but for a well-processed matrix they are close.

EpiDISH (reference-based, matrix-native). The EpiDISH package (Teschendorff et al. 2017) is designed to take a β-matrix plus a reference centroid matrix and return proportions via robust partial correlation (RPC), CIBERSORT, or constrained projection. It ships blood reference centroids and is the most convenient option when you only have processed data:
```
library(EpiDISH)
library(minfi)
beta <- minfi::getBeta(readRDS("data/03_grs_filtered.rds"))   # CpG x sample, as above
data(centDHSbloodDMC.m)                     # blood reference centroids
out <- epidish(beta, ref.m = centDHSbloodDMC.m, method = "RPC")
props <- out$estF
```

Reference-free alternatives

When no suitable reference exists (unusual tissues, or you worry the reference doesn’t match your population), reference-free methods estimate latent composition-like components directly from the data: RefFreeEWAS (Houseman et al. 2016), SVA (Leek et al. 2012), and RUV. These are covered in batch effects; the trade-off is that the recovered components are not labeled cell types, only surrogate variables to adjust for.

We carry the estimated proportions forward as EWAS covariates:

Code

# 1. The proportions themselves: a 96 x 6 matrix with GSM ids as rownames. This
#    is the file the later chapters read when they need cell-type covariates.
saveRDS(props, "data/04_cellcounts.rds")

# 2. The whole estimateCellCounts2() return, as deposited in the Zenodo record.
#    Keeping the full object means the file on record is the object the code
#    above produced, with nothing dropped or coerced.
saveRDS(cc, "data/04_cc_full.rds")

# 3. The validation: the per-sample comparison (`merged`) and the per-cell-type
#    summary (`validation`) that the table, the figure and the paragraph above
#    are all drawn from.
saveRDS(val, "data/04_validation.rds")

# 4. The deposited proportions, flat and small enough to live in the repository
#    rather than in a Zenodo tier -- which is what lets the fetch above run once
#    instead of on every render.
write.table(dep, "data/04_deposited_cellcounts.tsv",
            sep = "\t", row.names = FALSE, quote = FALSE)

# 5. The validation figure at the size the rendered page uses. The site serves
#    the committed copy, so it is written by the same code that drew it above.
ggsave("data/04_validation_scatter.png", p_val, width = 8, height = 5.2, dpi = 130)

Proportions are estimated for all 96 QC-passing arrays. Deconvolution does not use the exposure, so there is no reason to restrict it to the modeled subset here — the nine samples with no recorded PTSD status (chapter 01) drop out later, when the exposure enters the model, leaving n = 87 for the EWAS.

Next: Batch effects — technical variation from chips, array positions, and plates, and how to detect it before deciding whether to correct.

References

Houseman, Eugene Andres, William P Accomando, Devin C Koestler, et al. 2012. “DNA Methylation Arrays as Surrogate Measures of Cell Mixture Distribution.” BMC Bioinformatics 13: 86. https://doi.org/10.1186/1471-2105-13-86.

Houseman, Eugene Andres, Molly L Kile, David C Christiani, Tan A Ince, Karl T Kelsey, and Carmen J Marsit. 2016. “Reference-Free Deconvolution of DNA Methylation Data and Mediation by Cell Composition Effects.” BMC Bioinformatics 17: 259. https://doi.org/10.1186/s12859-016-1140-4.

Jaffe, Andrew E, and Rafael A Irizarry. 2014. “Accounting for Cellular Heterogeneity Is Critical in Epigenome-Wide Association Studies.” Genome Biology 15 (2): R31. https://doi.org/10.1186/gb-2014-15-2-r31.

Leek, Jeffrey T, W Evan Johnson, Hilary S Parker, Andrew E Jaffe, and John D Storey. 2012. “The Sva Package for Removing Batch Effects and Other Unwanted Variation in High-Throughput Experiments.” Bioinformatics 28 (6): 882–83. https://doi.org/10.1093/bioinformatics/bts034.

Salas, Lucas A, Devin C Koestler, Rondi A Butler, et al. 2018. “An Optimized Library for Reference-Based Deconvolution of Whole-Blood Biospecimens Assayed Using the Illumina HumanMethylationEPIC BeadArray.” Genome Biology 19 (1): 64. https://doi.org/10.1186/s13059-018-1448-7.

Teschendorff, Andrew E, Charles E Breeze, Shijie C Zheng, and Stephan Beck. 2017. “A Comparison of Reference-Based Algorithms for Correcting Cell-Type Heterogeneity in Epigenome-Wide Association Studies.” BMC Bioinformatics 18 (1): 105. https://doi.org/10.1186/s12859-017-1511-5.

--- title: "Cell composition: Whole Blood" --- ```{r} #| label: setup #| include: false source("_setup.R") ``` Whole blood is not one thing. It is a mixture of neutrophils, monocytes, natural killer (NK) cells, B cells, and CD4+/CD8+ T cells, and each cell type has a distinctly different methylome. When you draw blood from two people, differences in their cell-type *proportions* — driven by age, infection, stress, time of day — can produce large methylation differences that are unrelated to the biology you're studying. For most blood EWAS, cell composition is a large source of confounding [@houseman2012cellcomp; @jaffe2014celltype]. This notebook estimates the proportions, explains the reference-based method, and validates our estimates against the proportions the study deposited. ::: {.callout-note} ## What this chapter needs This chapter starts from `data/01_RGset.rds`, written by [chapter 00](00_setup.qmd). If you did not run that chapter, fetch the published checkpoint instead: ```{.bash filename="Terminal"} ./get_data.sh B_qc ``` ::: ```{r} #| label: load library(minfi) # estimateCellCounts2(), getBeta() library(FlowSorted.Blood.EPIC) # the flow-sorted EPIC reference and the IDOL probe library library(GEOquery) # getGEO(), for the study's own deposited proportions library(data.table) # fread() for the deposited table, rbindlist() for the figure library(ggplot2) # every figure in this chapter # The raw two-channel object written by the Setup chapter: 1,051,943 addresses # across 96 arrays, with no normalization applied. Deconvolution is one of the # few steps that needs the raw intensities rather than the normalized object # from chapter 02 -- the reason is in the callout below. rg <- readRDS("data/01_RGset.rds") ``` ## 1. Why cell composition matters A neutrophil and a CD4+ T cell carry the same genome but very different DNA methylation. A whole-blood methylation measurement at any CpG is therefore a **weighted average** across the cell types present, with weights equal to the cell proportions: $$\beta_{\text{observed}} \approx \sum_{k} w_k \, \beta_k$$ where $w_k$ is the fraction of cell type $k$ and $\beta_k$ is that cell type's methylation at the CpG. If cases have more neutrophils than controls — even for reasons unrelated to the exposure — every cell-type-discriminating CpG will look "differentially methylated." Neutrophil fraction alone routinely explains the largest chunk of genome-wide methylation variance in blood. ## 2. Reference-based deconvolution (Houseman) The standard approach is the **Houseman method** [@houseman2012cellcomp]: 1. Start from a **reference** of purified, flow-sorted blood cell types, each profiled on the same array platform. For each cell type we know its methylation at a panel of cell-type-discriminating CpGs. 2. Choose the discriminating CpGs. The original method picked them by F-test; the **IDOL** optimization [@salas2018idol] selects a library that minimizes deconvolution error and is the current best practice for EPIC. We use the IDOL-optimized 450-CpG library from `FlowSorted.Blood.EPIC`. 3. For each sample, solve a **constrained regression** (proportions ≥ 0, summing to ~1) of the observed methylation at those CpGs against the reference profiles. The fitted coefficients are the estimated cell proportions. Crucially, the sample and reference data must be **normalized together** so that technical differences don't confound composition. `minfi`'s `estimateCellCounts2` handles this by combining `RGChannelSet` with the reference `RGChannelSet`, normalizing jointly (here with noob), then deconvolving: ```{r} #| label: deconvolve # The IDOL-optimized probe library ships with FlowSorted.Blood.EPIC as a dataset # rather than a function, so it has to be attached before it can be passed in. data("IDOLOptimizedCpGs") # The published proportions were produced with this seed set; keeping it means a # re-render reproduces the same numbers rather than something merely similar. set.seed(1) # This is the most memory-hungry step in the whole tutorial: roughly 1 min 54 s # and a 15.0 GB peak for 96 arrays. The cost is structural, not wasteful -- the # 96 arrays and the flow-sorted reference are stacked into one object and noob- # normalized together (see the callout below), so the peak holds both datasets # at full probe resolution at once. An 8 GB laptop cannot run this chunk; the # collapsed callout further down is the route around it. cc <- estimateCellCounts2( rg, # raw RGChannelSet (not yet normalized) compositeCellType = "Blood", processMethod = "preprocessNoob", probeSelect = "IDOL", cellTypes = c("CD8T","CD4T","NK","Bcell","Mono","Neu"), referencePlatform = "IlluminaHumanMethylationEPIC", IDOLOptimizedCpGs = IDOLOptimizedCpGs, returnAll = FALSE) # `prop` is the estimated proportion per cell type; the sibling `counts` element # is only filled in when you supply absolute cell counts to compare against, so # on this call it comes back empty. Read `prop`. props <- cc$prop # 96 samples x 6 cell types dim(props) ``` ::: {.callout-note} ## Reference-based needs the *raw* data `estimateCellCounts2` starts from the raw `RGChannelSet`, not the normalized object from [normalization](02_normalization.qmd), because it must normalize the samples *and* reference. This is one of the few steps that reaches back to the raw IDATs. The reference dataset itself is downloaded once from Bioconductor's ExperimentHub and cached. ::: ::: {.callout-note collapse="true"} ## Running out of memory on the chunk above A 15.0 GB peak is more than a 16 GB laptop has spare, and far more than an 8 GB one has at all. The proportions this chapter writes at the end are committed to the repository, so you can skip the `deconvolve` chunk and read them instead — every section below it then runs: ```{.r filename="RStudio Console"} props <- readRDS(url(paste0("https://raw.githubusercontent.com/krferrier/", "Methylation-EWAS-tutorial/main/tutorial/data/", "04_cellcounts.rds"))) ``` What this route costs you is the one thing this chapter is about: you do not see the reference being fetched, you do not see the joint normalization, and if your own study ever needs different cell types or a different reference, you will not have run the call that changes. It also leaves `cc` undefined, so the save of the full deconvolution return at the end of the chapter cannot run. Use it because your machine forces you to. ::: ## 3. The composition of this subset ```{r} #| label: comp-summary # Collapse the 96 x 6 matrix to one row per cell type: the mean proportion and # the observed range across samples. Percentages read more easily than # proportions in a table, hence the 100x. comp <- data.frame(cell_type = colnames(props), mean_pct = round(100 * colMeans(props), 1), min_pct = round(100 * apply(props, 2, min), 1), max_pct = round(100 * apply(props, 2, max), 1)) comp <- comp[order(-comp$mean_pct), ] knitr::kable(comp, row.names = FALSE, caption = "Estimated cell-type composition across the 96 samples. Neutrophils dominate, as expected for whole blood, but the range is wide — this between-sample variation is exactly what must be modeled.") ``` ```{r} #| label: comp-stack #| fig-cap: "Estimated cell-type proportions per sample (stacked). The neutrophil fraction (largest band) varies from 17% to 84% across individuals — a swing that would swamp most disease effects if left unmodeled." #| fig-width: 8 #| fig-height: 3.4 # reshape() wants a data.frame with an id column, so the matrix goes wide-to-long # by hand here: one row per sample x cell type, which is what geom_col stacks. pl <- as.data.frame(props); pl$sample_id <- rownames(pl) plm <- reshape(pl, varying = colnames(props), v.names = "prop", timevar = "cell", times = colnames(props), direction = "long") # Stacking order is the factor order, so set it explicitly: largest fraction at # the bottom of the bar, rarest at the top. plm$cell <- factor(plm$cell, levels = c("Neu","CD4T","CD8T","Mono","Bcell","NK")) # Order the bars by neutrophil fraction so the gradient across samples is # visible instead of being scattered by array order. ord <- order(props[,"Neu"]); plm$sample_id <- factor(plm$sample_id, levels = rownames(props)[ord]) ggplot(plm, aes(sample_id, prop, fill = cell)) + geom_col(width = 1) + scale_fill_brewer(palette = "Set2") + labs(x = "sample (ordered by neutrophil fraction)", y = "estimated proportion", fill = "cell type", title = "Blood cell composition varies widely between samples") + theme_bw(base_size = 11) + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) ``` ## 4. Validation against the deposited answer key The Grady Trauma Project deposited its own cell-type proportion estimates with the GEO record. That gives us the opportunity to check our deconvolution. ```{r} #| label: deposited-fetch # The deposited proportions are part of the GEO record's *phenotype*, not of the # IDATs, so they come out of the series matrix rather than out of `rg`. GEOquery # parses each "label: value" characteristic into its own column named # "<label>:ch1", which is where the six proportions arrive. dep_file <- "data/04_deposited_cellcounts.tsv" if (file.exists(dep_file)) { # The series matrix is ~30 MB and the deposited values never change, so once # the small table below exists it is what later renders read. Delete it to # force the fetch again. dep <- as.data.frame(fread(dep_file)) } else { # 60 seconds is not enough to pull 30 MB from GEO on a slow link. options(timeout = max(1200, getOption("timeout"))) # getGPL = FALSE skips the platform manifest, which we do not need here. gse <- getGEO("GSE132203", GSEMatrix = TRUE, getGPL = FALSE)[[1]] pd <- Biobase::pData(gse) # Our cell-type names on the left, GEO's column names on the right. dep_cols <- c(CD8T = "cd8t:ch1", CD4T = "cd4t:ch1", NK = "nk:ch1", Bcell = "bcell:ch1", Mono = "mono:ch1", Neu = "neu:ch1") # Everything out of a series matrix is text, including the numbers. dep <- as.data.frame(lapply(dep_cols, function(k) as.numeric(as.character(pd[[k]])))) dep$sample_id <- as.character(pd$geo_accession) # The series carries 795 samples; match() keeps the 96 on our plate subset, # in the row order of `props` so the merge below cannot silently misalign. dep <- dep[match(rownames(props), dep$sample_id), c("sample_id", names(dep_cols))] } ``` ```{r} #| label: validation-compute cells <- c("CD8T","CD4T","NK","Bcell","Mono","Neu") # One row per sample with both sets of numbers side by side. Both frames carry # the same six column names, so merge() appends the suffixes that tell them # apart -- .est for ours, .dep for the study's. est <- as.data.frame(props); est$sample_id <- rownames(props) m <- merge(est, dep, by = "sample_id", suffixes = c(".est", ".dep")) # Three complementary summaries per cell type: correlation says whether the # between-sample variation agrees (the part that matters for confounding # adjustment), mean absolute error says how far apart the absolute values are, # and the two means expose a systematic offset the correlation would hide. res <- data.frame(cell = cells, r = NA_real_, mae = NA_real_, mean_est = NA_real_, mean_dep = NA_real_) for (i in seq_along(cells)) { e <- m[[paste0(cells[i], ".est")]] d <- m[[paste0(cells[i], ".dep")]] res$r[i] <- cor(e, d, use = "complete.obs") res$mae[i] <- mean(abs(e - d), na.rm = TRUE) res$mean_est[i] <- mean(e, na.rm = TRUE) res$mean_dep[i] <- mean(d, na.rm = TRUE) } # The per-sample comparison and the per-cell-type summary travel together: this # is the object saved at the end of the chapter, and the one the paragraph below # the figure reads its numbers from. val <- list(merged = m, validation = res) ``` ```{r} #| label: validation-table disp <- data.frame( `cell type` = res$cell, `Pearson r` = sprintf("%.3f", res$r), `mean abs. error` = sprintf("%.4f", res$mae), `our mean` = sprintf("%.1f%%", 100*res$mean_est), `deposited mean` = sprintf("%.1f%%", 100*res$mean_dep), check.names = FALSE) knitr::kable(disp, caption = "Our estimates vs the study's deposited proportions, per cell type (n = 96).") ``` ```{r} #| label: validation-scatter #| fig-cap: "Estimated (y) vs deposited (x) proportions for each cell type; dashed line is perfect agreement. Every cell type correlates > 0.95 with the answer key, with mean absolute error under 1.5 percentage points — this pipeline reproduces the study's own estimates." #| fig-width: 8 #| fig-height: 5.2 # Six panels, one per cell type, so stack the pairs into a long table. long <- rbindlist(lapply(cells, function(k) data.table(cell = k, est = m[[paste0(k, ".est")]], dep = m[[paste0(k, ".dep")]]))) long$cell <- factor(long$cell, levels = cells) # The correlation from the table above, printed in each panel's corner: with # free axes per panel a reader cannot judge agreement by eye alone. labs_r <- data.frame(cell = factor(res$cell, levels = cells), lab = sprintf("r = %.3f", res$r)) p_val <- ggplot(long, aes(dep, est)) + # y = x first, so the points sit on top of it. geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "#6E675B") + # gray geom_point(alpha = .6, size = 1.3, color = "#1A6B75") + # teal facet_wrap(~cell, scales = "free") + geom_text(data = labs_r, aes(x = -Inf, y = Inf, label = lab), hjust = -0.15, vjust = 1.5, size = 3, inherit.aes = FALSE) + labs(x = "deposited proportion (study answer key)", y = "our Houseman/IDOL estimate", title = "Deconvolution validation: estimated vs deposited cell proportions") p_val ``` Correlations range from **`r sprintf("%.3f", min(val$validation$r))`** (NK, the rarest and hardest-to-estimate type) to **`r sprintf("%.3f", max(val$validation$r))`** (neutrophils), and the mean absolute error never exceeds **`r sprintf("%.1f", 100*max(val$validation$mae))` percentage points**. Small systematic offsets (e.g. our B-cell estimates run slightly high) are expected when the reference library or normalization differs from the original authors', and they matter far less than the near-perfect *ranking* — because in the EWAS model it is the between-sample variation, not the absolute level, that does the confounding adjustment. ::: {.callout-tip} ## Practical note Many studies do not provide cell proportions. In a real analysis you trust the method and **sanity-check the output**: neutrophils should dominate whole blood (~50–70%), proportions should sum to ~1, and the distribution should look biologically plausible. ::: ## 5. When you don't have IDATs: deconvolution from a processed matrix Reference-based deconvolution normally wants the raw `RGChannelSet`, so a β-matrix-only deposit is the harder case. Two routes: 1. **`projectCellType` on a β-matrix.** The Houseman projection itself only needs methylation *values* at the reference CpGs, not intensities. Pull the reference β-values and its cell-type means, subset your matrix to the shared reference CpGs, and run the constrained projection directly: ```r library(FlowSorted.Blood.EPIC) library(genefilter) # `beta` is a CpG x sample matrix -- in the deposit-only case it is whatever # the submitter posted. Locally, the filtered object from chapter 03 gives one: beta <- minfi::getBeta(readRDS("data/03_grs_filtered.rds")) # reference cell-type means at the IDOL CpGs (from the package) ref <- ... # matrix: IDOL CpGs x 6 cell types shared <- intersect(rownames(beta), rownames(ref)) props <- minfi:::projectCellType(beta[shared, ], ref[shared, ]) ``` The estimates are slightly noisier than the joint-normalization route because the sample and reference weren't normalized together, but for a well-processed matrix they are close. 2. **`EpiDISH` (reference-based, matrix-native).** The `EpiDISH` package [@teschendorff2017reference] is designed to take a β-matrix plus a reference centroid matrix and return proportions via robust partial correlation (RPC), CIBERSORT, or constrained projection. It ships blood reference centroids and is the most convenient option when you only have processed data: ```r library(EpiDISH) library(minfi) beta <- minfi::getBeta(readRDS("data/03_grs_filtered.rds")) # CpG x sample, as above data(centDHSbloodDMC.m) # blood reference centroids out <- epidish(beta, ref.m = centDHSbloodDMC.m, method = "RPC") props <- out$estF ``` ::: {.callout-note} ## Reference-free alternatives When no suitable reference exists (unusual tissues, or you worry the reference doesn't match your population), **reference-free** methods estimate latent composition-like components directly from the data: RefFreeEWAS [@houseman2014reffree], SVA [@leek2012sva], and RUV. These are covered in [batch effects](05_batch_effects.qmd); the trade-off is that the recovered components are not labeled cell types, only surrogate variables to adjust for. ::: We carry the estimated proportions forward as EWAS covariates: ```{r} #| label: save-cellcounts # 1. The proportions themselves: a 96 x 6 matrix with GSM ids as rownames. This # is the file the later chapters read when they need cell-type covariates. saveRDS(props, "data/04_cellcounts.rds") # 2. The whole estimateCellCounts2() return, as deposited in the Zenodo record. # Keeping the full object means the file on record is the object the code # above produced, with nothing dropped or coerced. saveRDS(cc, "data/04_cc_full.rds") # 3. The validation: the per-sample comparison (`merged`) and the per-cell-type # summary (`validation`) that the table, the figure and the paragraph above # are all drawn from. saveRDS(val, "data/04_validation.rds") # 4. The deposited proportions, flat and small enough to live in the repository # rather than in a Zenodo tier -- which is what lets the fetch above run once # instead of on every render. write.table(dep, "data/04_deposited_cellcounts.tsv", sep = "\t", row.names = FALSE, quote = FALSE) # 5. The validation figure at the size the rendered page uses. The site serves # the committed copy, so it is written by the same code that drew it above. ggsave("data/04_validation_scatter.png", p_val, width = 8, height = 5.2, dpi = 130) ``` Proportions are estimated for all **96** QC-passing arrays. Deconvolution does not use the exposure, so there is no reason to restrict it to the modeled subset here — the nine samples with no recorded PTSD status ([chapter 01](01_qc.qmd)) drop out later, when the exposure enters the model, leaving n = 87 for the EWAS. --- **Next:** [Batch effects](05_batch_effects.qmd) — technical variation from chips, array positions, and plates, and how to detect it before deciding whether to correct. ## References {.unnumbered} ::: {#refs} :::