Quality control: which samples can we trust?

The first step in methylation analysis is quality control. A failed array, a swapped sample, or sex-mismatch can cause serious trouble downstream if it slips through. This notebook runs the standard minfi (Aryee et al. 2014) QC process on the 96-sample subset and walks you through how to decide what (if anything) to drop.

We start from the raw two-channel object saved at the end of the Setup step — data/01_RGset.rds. Nothing in this chapter is pre-computed: every number and every figure below is derived from that object by the code you see, and the four QC measures are saved at the end so later work does not have to recompute them.

Code

library(minfi)     # detectionP(), getQC(), getSex(), getSnpBeta()
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.
rg <- readRDS("data/01_RGset.rds")

# Setup attached the sample sheet to the object, so the phenotype table comes
# back out of it here rather than being read from the CSV a second time. Its row
# order is guaranteed to match the column order of `rg`.
ss <- as.data.frame(pData(rg))

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

The four checks below are independent of one another, and each one reduces the million-probe object to a handful of numbers per sample. Running all four takes under a minute.

1. Detection p-values — did the probe beat background?

Every probe is compared to the array’s negative-control distribution. A large detection p-value means the signal is indistinguishable from background; the measurement is noise. We summarize per sample: the mean detection p, and the fraction of probes failing at p > 0.01.

Code

# detectionP() compares every probe on every array against that array's own
# distribution of negative-control probes, and returns a probe x sample matrix
# of p-values. This is the one expensive step in the chapter: roughly 40 seconds
# and 4 GB for 96 arrays, because the full 1,051,943 x 96 matrix is built.
detP <- detectionP(rg)

# Two per-sample summaries are all QC needs, so we collapse the big matrix
# immediately and keep only these:
det_mean <- colMeans(detP)          # average detection p-value across all probes
det_fail <- colMeans(detP > 0.01)   # fraction of probes that failed the 0.01 cut

# One row per array, ready to plot.
det <- data.frame(sample    = colnames(rg),
                  mean_detP = as.numeric(det_mean),
                  frac_fail = as.numeric(det_fail),
                  slide     = ss$slide)
head(det, 3)

Code

ggplot(det, aes(x = reorder(sample, mean_detP), y = mean_detP)) +
  geom_point(color = "#1A6B75", size = 1.6) +                                  # teal
  geom_hline(yintercept = 0.01, linetype = "dashed", color = "#8C3A4A") +      # plum
  labs(x = "sample (ranked)", y = "mean detection p-value",
       title = "Detection p-value per sample") +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

Per-sample detection p-value. All 96 samples are far below the 0.01 concern line — no failed arrays.

Code

cat("worst sample mean detP:", signif(max(det$mean_detP), 3), "\n")

worst sample mean detP: 0.00107

Code

cat("worst sample % probes p>0.01:", signif(max(det$frac_fail) * 100, 3), "%\n")

worst sample % probes p>0.01: 0.598 %

Code

cat("samples failing QC (mean detP > 0.01):", sum(det$mean_detP > 0.01), "\n")

samples failing QC (mean detP > 0.01): 0

What a failure would look like

On a real plate you will sometimes see one or two samples an order of magnitude above the rest, crossing the 0.01 line. Those get dropped here, before normalization.

2. Signal intensity — is the array bright enough?

getQC() reduces each sample to the median methylated (mMed) and unmethylated (uMed) log2 intensity. Weak arrays fall toward the bottom-left; minfi’s convention flags samples below a line at (mMed + uMed)/2 = 10.5.

Code

# getQC() needs the two channels separated into methylated and unmethylated
# signal, which is what preprocessRaw() does -- it is a re-arrangement of the
# raw intensities, not a normalization, so nothing is altered here.
mset  <- preprocessRaw(rg)
qc_df <- getQC(mset)           # one row per sample: mMed and uMed

qc <- as.data.frame(qc_df)
qc$sample <- rownames(qc)
# minfi's convention: flag an array whose two medians average below 10.5.
qc$flag <- ((qc$mMed + qc$uMed) / 2) < 10.5
head(qc, 3)

Code

ggplot(qc, aes(mMed, uMed, color = flag)) +
  geom_point(size = 2) +
  geom_abline(intercept = 2 * 10.5, slope = -1, linetype = "dashed", color = "gray40") +
  scale_color_manual(values = c(`FALSE` = "#1A6B75", `TRUE` = "#8C3A4A"),
                     labels = c("pass", "low signal"), name = NULL) +
  labs(x = "median methylated log2 int (mMed)",
       y = "median unmethylated log2 int (uMed)",
       title = "Signal-intensity QC") +
  coord_equal(xlim = c(9, 14), ylim = c(9, 14))

Median methylated vs unmethylated log2 intensity. All samples sit well above the 10.5 low-signal cutoff.

Code

cat("mMed range:", round(min(qc$mMed),2), "-", round(max(qc$mMed),2), "\n")

mMed range: 11.27 - 12.78

Code

cat("low-signal flagged samples:", sum(qc$flag), "\n")

low-signal flagged samples: 0

3. Predicted vs reported sex

Sex is the single most reliable phenotype to predict from methylation: the X and Y chromosome signal is unambiguous. getSex() compares median total intensity on X vs Y. A mismatch between predicted and reported sex is a strong signal of a sample swap or annotation error — and, because sex is a strong methylation covariate we must model, catching it here is important.

Code

# getSex() works on median intensity per chromosome, so the probes have to know
# where they sit in the genome first. mapToGenome() attaches those coordinates
# using the EPIC annotation package; it adds no new measurement.
gmset <- mapToGenome(mset)
psex  <- getSex(gmset)         # xMed, yMed and predictedSex per sample

sx <- data.frame(sample    = colnames(rg),
                 xMed      = psex$xMed,
                 yMed      = psex$yMed,
                 predicted = as.character(psex$predictedSex),
                 reported  = as.character(ss$sex))   # reported sex from the sample sheet
sx$mismatch <- sx$predicted != sx$reported
head(sx, 3)

Code

ggplot(sx, aes(xMed, yMed, color = reported, shape = mismatch)) +
  geom_point(size = 2.5) +
  scale_color_manual(values = c(F = "#B8873F", M = "#1A6B75"), name = "reported") +
  scale_shape_manual(values = c(`FALSE` = 16, `TRUE` = 4), guide = "none") +
  labs(x = "median X log2 int", y = "median Y log2 int",
       title = "Predicted vs reported sex")

Predicted sex from X/Y intensity, colored by reported sex. Clean separation, zero mismatches.

Code

tab <- table(reported = sx$reported, predicted = sx$predicted)
knitr::kable(tab, caption = "Reported vs predicted sex")

Reported vs predicted sex
	F	M
F	49	0
M	0	47

Code

cat("sex mismatches:", sum(sx$mismatch), "of", nrow(sx), "\n")

sex mismatches: 0 of 96

4. SNP-identity fingerprint

EPIC carries ~59 rs SNP probes that measure genotype, not methylation. Their beta values cluster near 0, 0.5, or 1 (the three genotypes). They are useful for two things: spotting sample duplicates / swaps (two arrays from the same person share a fingerprint) and, in longitudinal designs, confirming that repeat samples really are the same individual.

Code

# getSnpBeta() pulls just the rs probes out of the raw object and returns their
# beta values -- a 59 x 96 genotype matrix, small enough to handle directly.
snps <- getSnpBeta(rg)

# Cluster the samples by fingerprint so that any identical pair would sit
# side by side in the image below.
d   <- dist(t(snps))
ord <- hclust(d)$order
image(t(snps[, ord]), col = colorRampPalette(c("#1A6B75", "white", "#8C3A4A"))(64),
      axes = FALSE, main = "rs SNP-probe fingerprints (samples clustered)")
box()

The ~59 rs SNP probes across samples. Discrete banding near 0 / 0.5 / 1 is the expected genotype signal; no two columns are identical, so no duplicate samples.

Code

# nearest-neighbor distance: a near-zero value would flag a duplicate
dm <- as.matrix(dist(t(snps)))
diag(dm) <- NA
nn <- apply(dm, 1, min, na.rm = TRUE)
cat("SNP probes:", nrow(snps), "\n")

SNP probes: 59

Code

cat("smallest between-sample fingerprint distance:", round(min(nn), 3),
    "(a value near 0 would indicate a duplicate)\n")

smallest between-sample fingerprint distance: 2.224 (a value near 0 would indicate a duplicate)

Verdict

Code

# Each check contributes a list of sample IDs to drop.
drop_detp <- det$sample[det$mean_detP > 0.01]
drop_int  <- qc$sample[qc$flag]
drop_sex  <- sx$sample[sx$mismatch]     # sample IDs, to investigate rather than drop

# Only the first two are automatic. A sex mismatch is a signal to go and find
# out what happened, not grounds for silent removal, so it is reported and
# deliberately left out of `keep`.
keep <- setdiff(colnames(rg), unique(c(drop_detp, drop_int)))

cat("Samples in:", ncol(rg), "\n")

Samples in: 96

Code

cat("Dropped for detection p:", length(drop_detp), "\n")

Dropped for detection p: 0

Code

cat("Dropped for low intensity:", length(drop_int), "\n")

Dropped for low intensity: 0

Code

cat("Sex mismatches (investigate, don't auto-drop):", length(drop_sex), "\n")

Sex mismatches (investigate, don't auto-drop): 0

Code

cat("Samples passing QC:", length(keep), "\n")

Samples passing QC: 96

All 96 samples pass. We carry the full set forward to normalization.

Note that “96 samples pass QC” is not the same as “96 samples are analyzed”. Nine of them have no recorded PTSD status and will drop out when the exposure enters the model — see phenotype completeness below.

Always do QC. ALWAYS

Even though all samples passed QC thresholds in this subset, doing the QC was still important as now we have evidence supporting that these samples are ‘clean’. In a real dataset there are typically several samples that will not pass QC for one of the 4 checks.

Saving what we just computed

The four checks cost a minute of compute and reduce a 460 MB object to a few numbers per sample. Save them, so that revisiting a figure or a threshold later does not mean reading the IDATs and recomputing detection p-values again.

Two files, for two different readers:

Code

# 1. The full set of pieces, exactly as deposited in the Zenodo record. Keeping
#    the S4 objects (`qc` and `psex` are DFrames) means the saved file is the
#    same object the code above produced, with nothing coerced or rounded.
qc_pieces <- list(
  det_mean = det_mean,                 # named numeric, one per sample
  det_fail = det_fail,                 # named numeric, one per sample
  qc       = qc_df,                    # DFrame: mMed, uMed
  psex     = psex,                     # DFrame: xMed, yMed, predictedSex
  rep_sex  = as.character(ss$sex),     # reported sex, from the sample sheet
  snps     = snps                      # 59 rs probes x 96 samples
)
saveRDS(qc_pieces, "data/01_qc_pieces.rds")

# 2. A flat, per-sample table of the same measurements. This one is small enough
#    to live in the repository, so a reader who cannot hold the RGChannelSet in
#    memory can still reproduce every figure in this chapter (see the callout).
qc_metrics <- data.frame(
  sample         = colnames(rg),
  mean_detP      = as.numeric(det_mean),
  frac_detP_fail = as.numeric(det_fail),
  mMed           = as.numeric(qc_df$mMed),
  uMed           = as.numeric(qc_df$uMed),
  xMed           = as.numeric(psex$xMed),
  yMed           = as.numeric(psex$yMed),
  predicted_sex  = as.character(psex$predictedSex),
  reported_sex   = as.character(ss$sex),
  slide          = ss$slide,
  array_pos      = ss$array_pos,
  qc_pass        = colnames(rg) %in% keep,
  stringsAsFactors = FALSE)

write.table(qc_metrics, "data/01_qc_sample_metrics.tsv",
            sep = "\t", row.names = FALSE, quote = FALSE)

Picking the chapter up from the saved table instead

If reading the IDATs or holding the raw object in memory is not practical on your machine, you can skip straight to the figures. The per-sample table is in the repository, so it can be read directly over the network — no clone and no download step:

RStudio Console

library(data.table)
qc_metrics <- fread(paste0("https://raw.githubusercontent.com/krferrier/",
                           "Methylation-EWAS-tutorial/main/tutorial/data/",
                           "01_qc_sample_metrics.tsv"))

# The plotting frames used above, rebuilt from the table:
det <- data.frame(sample = qc_metrics$sample, mean_detP = qc_metrics$mean_detP,
                  frac_fail = qc_metrics$frac_detP_fail, slide = qc_metrics$slide)
qc  <- data.frame(sample = qc_metrics$sample, mMed = qc_metrics$mMed,
                  uMed = qc_metrics$uMed)
qc$flag <- ((qc$mMed + qc$uMed) / 2) < 10.5
sx  <- data.frame(sample = qc_metrics$sample, xMed = qc_metrics$xMed,
                  yMed = qc_metrics$yMed, predicted = qc_metrics$predicted_sex,
                  reported = qc_metrics$reported_sex)
sx$mismatch <- sx$predicted != sx$reported

Everything except the SNP-fingerprint heatmap will then run, because that one needs the 59 x 96 genotype matrix rather than a per-sample summary. What this route does not give you is the experience of watching the checks run — so use it because you have to, not to save a minute.

What the phenotype file is missing

QC checks whether the arrays worked. It is also the right moment to check whether the phenotype data you will eventually model is complete, because a covariate you discover is absent at modeling time is far more expensive than one you plan around now.

There are two distinct kinds of gap, and they have different consequences. A variable that is absent for every sample is unmeasured confounding — you cannot adjust for it. A variable that is present but incomplete costs you samples, because any model containing that term silently drops the rows where it is missing. Check for both.

Missing values in the variables we intend to model

Tabulate missingness for every term in the planned model, before you normalize:

Code

model_vars <- c("ptsd", "sex", "age", "childhood_abuse",
                "mergedcapsandpsswinthin30days")
miss <- data.frame(
  variable = model_vars,
  n_missing = sapply(model_vars, function(v)
    sum(is.na(ss[[v]]) | ss[[v]] == "")),
  row.names = NULL)
knitr::kable(miss, caption = "Missing values per candidate model variable, across all 96 QC-passing arrays.")

Missing values per candidate model variable, across all 96 QC-passing arrays.
variable	n_missing
ptsd	9
sex	0
age	0
childhood_abuse	2
mergedcapsandpsswinthin30days	9

PTSD status — the exposure — is missing for 9 of the 96 samples. The continuous CAPS/PSS severity score is missing for the same samples, so there is no alternative outcome to fall back on for them.

These nine samples are dropped later, and it is not a QC failure

Every one of the 96 arrays passed QC. But an EWAS model contains the exposure, so lm / limma will drop any sample with no PTSD status — the analysis in chapter 06 and the pipeline run in chapter 07 therefore fit n = 87, not 96.

This tutorial keeps all 96 samples through normalization (chapter 02), probe filtering (chapter 03), cell-composition estimation (chapter 04), and batch diagnosis (chapter 05) on purpose. None of those steps uses PTSD status: deconvolution and PCA describe the arrays, not the association model, so there is no reason to discard usable arrays before the step that actually needs the exposure. The count changes exactly once, at the point the exposure enters.

The distinction worth internalising: samples dropped for quality are samples you do not trust, and samples dropped for missing phenotype are samples you trust but cannot model. Report them separately. Collapsing them into one “n = 87 after QC” sentence — the most common way this gets written up — hides the fact that nine usable arrays were lost to an incomplete metadata table, which is a very different problem from nine failed arrays.

Because the nine are a loss from missingness rather than quality, it is worth a moment to ask whether they are missing at random or missing in a way that correlates with something. If the samples without PTSD status were, say, all male or all from one chip, dropping them would distort the design:

Code

no_ptsd <- is.na(ss$ptsd) | ss$ptsd == ""
cat("samples with no PTSD status:", sum(no_ptsd), "\n")

samples with no PTSD status: 9

Code

cat("sex:", paste(names(table(ss$sex[no_ptsd])),
                 table(ss$sex[no_ptsd]), collapse = " / "), "\n")

sex: F 4 / M 5

Code

cat("age range:", paste(range(ss$age[no_ptsd]), collapse = "-"), "\n")

age range: 22-65

Code

cat("distinct chips:", length(unique(ss$slide[no_ptsd])),
    "of", length(unique(ss$slide)), "\n")

distinct chips: 8 of 12

Four female and five male, ages spanning the full cohort range, spread across eight of the twelve chips — the missingness looks incidental rather than structured. That is reassuring, but note that it is a check you can only do for observed variables; it cannot rule out the samples differing on something unrecorded.

Variables absent for every sample

For blood methylation the most consequential gap is usually smoking. Smoking produces some of the largest and most reproducible methylation signals in blood (Joehanes et al. 2016), so an unmeasured smoking variable is unmeasured confounding in every model you fit. The Grady deposit gives PTSD, age, sex, race, and childhood trauma — but no smoking variable, which is common in public deposits rather than unusual.

That absence is worth noting here because it shapes decisions in the next two chapters. The full set of options, and the one this tutorial takes, is laid out in the EWAS chapter.

Inventory your covariates before you normalize

Write down the model you intend to fit, then check each term against the metadata you actually have. For each term ask two questions: is it present at all? and how many samples is it missing for? The first tells you what you cannot adjust for; the second tells you your real analysis n. Note both, and decide then how you will handle them — including whether a term is worth the samples it costs.

References

Aryee, Martin J, Andrew E Jaffe, Hector Corrada-Bravo, et al. 2014. “Minfi: A Flexible and Comprehensive Bioconductor Package for the Analysis of Infinium DNA Methylation Microarrays.” Bioinformatics 30 (10): 1363–69. https://doi.org/10.1093/bioinformatics/btu049.

Joehanes, Roby, Allan C Just, Riccardo E Marioni, et al. 2016. “Epigenetic Signatures of Cigarette Smoking.” Circulation: Cardiovascular Genetics 9 (5): 436–47. https://doi.org/10.1161/CIRCGENETICS.116.001506.

--- title: "Quality control: which samples can we trust?" --- ```{r} #| label: setup #| include: false source("_setup.R") ``` The first step in methylation analysis is quality control. A failed array, a swapped sample, or sex-mismatch can cause serious trouble downstream if it slips through. This notebook runs the standard `minfi` [@aryee2014minfi] QC process on the 96-sample subset and walks you through how to decide what (if anything) to drop. We start from the raw two-channel object saved at the end of the [Setup](00_setup.qmd) step — `data/01_RGset.rds`. Nothing in this chapter is pre-computed: every number and every figure below is derived from that object by the code you see, and the four QC measures are saved at the end so later work does not have to recompute them. ```{r} #| label: load library(minfi) # detectionP(), getQC(), getSex(), getSnpBeta() 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. rg <- readRDS("data/01_RGset.rds") # Setup attached the sample sheet to the object, so the phenotype table comes # back out of it here rather than being read from the CSV a second time. Its row # order is guaranteed to match the column order of `rg`. ss <- as.data.frame(pData(rg)) rg ``` The four checks below are independent of one another, and each one reduces the million-probe object to a handful of numbers per sample. Running all four takes under a minute. ## 1. Detection p-values — did the probe beat background? Every probe is compared to the array's negative-control distribution. A large **detection p-value** means the signal is indistinguishable from background; the measurement is noise. We summarize per sample: the mean detection p, and the fraction of probes failing at p > 0.01. ```{r} #| label: detp-compute # detectionP() compares every probe on every array against that array's own # distribution of negative-control probes, and returns a probe x sample matrix # of p-values. This is the one expensive step in the chapter: roughly 40 seconds # and 4 GB for 96 arrays, because the full 1,051,943 x 96 matrix is built. detP <- detectionP(rg) # Two per-sample summaries are all QC needs, so we collapse the big matrix # immediately and keep only these: det_mean <- colMeans(detP) # average detection p-value across all probes det_fail <- colMeans(detP > 0.01) # fraction of probes that failed the 0.01 cut # One row per array, ready to plot. det <- data.frame(sample = colnames(rg), mean_detP = as.numeric(det_mean), frac_fail = as.numeric(det_fail), slide = ss$slide) head(det, 3) ``` ```{r} #| label: detp-plot #| fig-cap: "Per-sample detection p-value. All 96 samples are far below the 0.01 concern line — no failed arrays." ggplot(det, aes(x = reorder(sample, mean_detP), y = mean_detP)) + geom_point(color = "#1A6B75", size = 1.6) + # teal geom_hline(yintercept = 0.01, linetype = "dashed", color = "#8C3A4A") + # plum labs(x = "sample (ranked)", y = "mean detection p-value", title = "Detection p-value per sample") + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) ``` ```{r} #| label: detp-table cat("worst sample mean detP:", signif(max(det$mean_detP), 3), "\n") cat("worst sample % probes p>0.01:", signif(max(det$frac_fail) * 100, 3), "%\n") cat("samples failing QC (mean detP > 0.01):", sum(det$mean_detP > 0.01), "\n") ``` ::: {.callout-note} ## What a failure would look like On a real plate you will sometimes see one or two samples an order of magnitude above the rest, crossing the 0.01 line. Those get dropped here, before normalization. ::: ## 2. Signal intensity — is the array bright enough? `getQC()` reduces each sample to the median methylated (`mMed`) and unmethylated (`uMed`) log2 intensity. Weak arrays fall toward the bottom-left; `minfi`'s convention flags samples below a line at `(mMed + uMed)/2 = 10.5`. ```{r} #| label: intensity-compute # getQC() needs the two channels separated into methylated and unmethylated # signal, which is what preprocessRaw() does -- it is a re-arrangement of the # raw intensities, not a normalization, so nothing is altered here. mset <- preprocessRaw(rg) qc_df <- getQC(mset) # one row per sample: mMed and uMed qc <- as.data.frame(qc_df) qc$sample <- rownames(qc) # minfi's convention: flag an array whose two medians average below 10.5. qc$flag <- ((qc$mMed + qc$uMed) / 2) < 10.5 head(qc, 3) ``` ```{r} #| label: intensity-plot #| fig-cap: "Median methylated vs unmethylated log2 intensity. All samples sit well above the 10.5 low-signal cutoff." ggplot(qc, aes(mMed, uMed, color = flag)) + geom_point(size = 2) + geom_abline(intercept = 2 * 10.5, slope = -1, linetype = "dashed", color = "gray40") + scale_color_manual(values = c(`FALSE` = "#1A6B75", `TRUE` = "#8C3A4A"), labels = c("pass", "low signal"), name = NULL) + labs(x = "median methylated log2 int (mMed)", y = "median unmethylated log2 int (uMed)", title = "Signal-intensity QC") + coord_equal(xlim = c(9, 14), ylim = c(9, 14)) ``` ```{r} #| label: intensity-summary cat("mMed range:", round(min(qc$mMed),2), "-", round(max(qc$mMed),2), "\n") cat("low-signal flagged samples:", sum(qc$flag), "\n") ``` ## 3. Predicted vs reported sex Sex is the single most reliable phenotype to predict from methylation: the X and Y chromosome signal is unambiguous. `getSex()` compares median total intensity on X vs Y. A mismatch between predicted and reported sex is a strong signal of a **sample swap or annotation error** — and, because sex is a strong methylation covariate we must model, catching it here is important. ```{r} #| label: sex-compute # getSex() works on median intensity per chromosome, so the probes have to know # where they sit in the genome first. mapToGenome() attaches those coordinates # using the EPIC annotation package; it adds no new measurement. gmset <- mapToGenome(mset) psex <- getSex(gmset) # xMed, yMed and predictedSex per sample sx <- data.frame(sample = colnames(rg), xMed = psex$xMed, yMed = psex$yMed, predicted = as.character(psex$predictedSex), reported = as.character(ss$sex)) # reported sex from the sample sheet sx$mismatch <- sx$predicted != sx$reported head(sx, 3) ``` ```{r} #| label: sex-plot #| fig-cap: "Predicted sex from X/Y intensity, colored by reported sex. Clean separation, zero mismatches." ggplot(sx, aes(xMed, yMed, color = reported, shape = mismatch)) + geom_point(size = 2.5) + scale_color_manual(values = c(F = "#B8873F", M = "#1A6B75"), name = "reported") + scale_shape_manual(values = c(`FALSE` = 16, `TRUE` = 4), guide = "none") + labs(x = "median X log2 int", y = "median Y log2 int", title = "Predicted vs reported sex") ``` ```{r} #| label: sex-table tab <- table(reported = sx$reported, predicted = sx$predicted) knitr::kable(tab, caption = "Reported vs predicted sex") cat("sex mismatches:", sum(sx$mismatch), "of", nrow(sx), "\n") ``` ## 4. SNP-identity fingerprint EPIC carries ~59 `rs` SNP probes that measure genotype, not methylation. Their beta values cluster near 0, 0.5, or 1 (the three genotypes). They are useful for two things: spotting **sample duplicates / swaps** (two arrays from the same person share a fingerprint) and, in longitudinal designs, confirming that repeat samples really are the same individual. ```{r} #| label: snp-heatmap #| fig-cap: "The ~59 rs SNP probes across samples. Discrete banding near 0 / 0.5 / 1 is the expected genotype signal; no two columns are identical, so no duplicate samples." # getSnpBeta() pulls just the rs probes out of the raw object and returns their # beta values -- a 59 x 96 genotype matrix, small enough to handle directly. snps <- getSnpBeta(rg) # Cluster the samples by fingerprint so that any identical pair would sit # side by side in the image below. d <- dist(t(snps)) ord <- hclust(d)$order image(t(snps[, ord]), col = colorRampPalette(c("#1A6B75", "white", "#8C3A4A"))(64), axes = FALSE, main = "rs SNP-probe fingerprints (samples clustered)") box() ``` ```{r} #| label: snp-dups # nearest-neighbor distance: a near-zero value would flag a duplicate dm <- as.matrix(dist(t(snps))) diag(dm) <- NA nn <- apply(dm, 1, min, na.rm = TRUE) cat("SNP probes:", nrow(snps), "\n") cat("smallest between-sample fingerprint distance:", round(min(nn), 3), "(a value near 0 would indicate a duplicate)\n") ``` ## Verdict ```{r} #| label: verdict # Each check contributes a list of sample IDs to drop. drop_detp <- det$sample[det$mean_detP > 0.01] drop_int <- qc$sample[qc$flag] drop_sex <- sx$sample[sx$mismatch] # sample IDs, to investigate rather than drop # Only the first two are automatic. A sex mismatch is a signal to go and find # out what happened, not grounds for silent removal, so it is reported and # deliberately left out of `keep`. keep <- setdiff(colnames(rg), unique(c(drop_detp, drop_int))) cat("Samples in:", ncol(rg), "\n") cat("Dropped for detection p:", length(drop_detp), "\n") cat("Dropped for low intensity:", length(drop_int), "\n") cat("Sex mismatches (investigate, don't auto-drop):", length(drop_sex), "\n") cat("Samples passing QC:", length(keep), "\n") ``` All 96 samples pass. We carry the full set forward to [normalization](02_normalization.qmd). Note that "96 samples pass QC" is not the same as "96 samples are analyzed". Nine of them have no recorded PTSD status and will drop out when the exposure enters the model — see [phenotype completeness](#what-the-phenotype-file-is-missing) below. ::: {.callout-tip} ## Always do QC. ALWAYS Even though all samples passed QC thresholds in this subset, doing the QC was still important as now we have evidence supporting that these samples are 'clean'. In a real dataset there are typically several samples that will not pass QC for one of the 4 checks. ::: ## Saving what we just computed The four checks cost a minute of compute and reduce a 460 MB object to a few numbers per sample. Save them, so that revisiting a figure or a threshold later does not mean reading the IDATs and recomputing detection p-values again. Two files, for two different readers: ```{r} #| label: save-qc # 1. The full set of pieces, exactly as deposited in the Zenodo record. Keeping # the S4 objects (`qc` and `psex` are DFrames) means the saved file is the # same object the code above produced, with nothing coerced or rounded. qc_pieces <- list( det_mean = det_mean, # named numeric, one per sample det_fail = det_fail, # named numeric, one per sample qc = qc_df, # DFrame: mMed, uMed psex = psex, # DFrame: xMed, yMed, predictedSex rep_sex = as.character(ss$sex), # reported sex, from the sample sheet snps = snps # 59 rs probes x 96 samples ) saveRDS(qc_pieces, "data/01_qc_pieces.rds") # 2. A flat, per-sample table of the same measurements. This one is small enough # to live in the repository, so a reader who cannot hold the RGChannelSet in # memory can still reproduce every figure in this chapter (see the callout). qc_metrics <- data.frame( sample = colnames(rg), mean_detP = as.numeric(det_mean), frac_detP_fail = as.numeric(det_fail), mMed = as.numeric(qc_df$mMed), uMed = as.numeric(qc_df$uMed), xMed = as.numeric(psex$xMed), yMed = as.numeric(psex$yMed), predicted_sex = as.character(psex$predictedSex), reported_sex = as.character(ss$sex), slide = ss$slide, array_pos = ss$array_pos, qc_pass = colnames(rg) %in% keep, stringsAsFactors = FALSE) write.table(qc_metrics, "data/01_qc_sample_metrics.tsv", sep = "\t", row.names = FALSE, quote = FALSE) ``` ::: {.callout-note collapse="true"} ## Picking the chapter up from the saved table instead If reading the IDATs or holding the raw object in memory is not practical on your machine, you can skip straight to the figures. The per-sample table is in the repository, so it can be read directly over the network — no clone and no download step: ```{.r filename="RStudio Console"} library(data.table) qc_metrics <- fread(paste0("https://raw.githubusercontent.com/krferrier/", "Methylation-EWAS-tutorial/main/tutorial/data/", "01_qc_sample_metrics.tsv")) # The plotting frames used above, rebuilt from the table: det <- data.frame(sample = qc_metrics$sample, mean_detP = qc_metrics$mean_detP, frac_fail = qc_metrics$frac_detP_fail, slide = qc_metrics$slide) qc <- data.frame(sample = qc_metrics$sample, mMed = qc_metrics$mMed, uMed = qc_metrics$uMed) qc$flag <- ((qc$mMed + qc$uMed) / 2) < 10.5 sx <- data.frame(sample = qc_metrics$sample, xMed = qc_metrics$xMed, yMed = qc_metrics$yMed, predicted = qc_metrics$predicted_sex, reported = qc_metrics$reported_sex) sx$mismatch <- sx$predicted != sx$reported ``` Everything except the SNP-fingerprint heatmap will then run, because that one needs the 59 x 96 genotype matrix rather than a per-sample summary. What this route does *not* give you is the experience of watching the checks run — so use it because you have to, not to save a minute. ::: ## What the phenotype file is missing QC checks whether the *arrays* worked. It is also the right moment to check whether the **phenotype data** you will eventually model is complete, because a covariate you discover is absent at modeling time is far more expensive than one you plan around now. There are two distinct kinds of gap, and they have different consequences. A variable that is **absent for every sample** is unmeasured confounding — you cannot adjust for it. A variable that is **present but incomplete** costs you samples, because any model containing that term silently drops the rows where it is missing. Check for both. ### Missing values in the variables we intend to model Tabulate missingness for every term in the planned model, before you normalize: ```{r} #| label: pheno-missing model_vars <- c("ptsd", "sex", "age", "childhood_abuse", "mergedcapsandpsswinthin30days") miss <- data.frame( variable = model_vars, n_missing = sapply(model_vars, function(v) sum(is.na(ss[[v]]) | ss[[v]] == "")), row.names = NULL) knitr::kable(miss, caption = "Missing values per candidate model variable, across all 96 QC-passing arrays.") ``` PTSD status — the exposure — is missing for `r sum(is.na(ss$ptsd) | ss$ptsd == "")` of the 96 samples. The continuous CAPS/PSS severity score is missing for the same samples, so there is no alternative outcome to fall back on for them. ::: {.callout-important} ## These nine samples are dropped later, and it is not a QC failure Every one of the 96 arrays passed QC. But an EWAS model contains the exposure, so `lm` / `limma` will drop any sample with no PTSD status — the analysis in [chapter 06](06_ewas.qmd) and the pipeline run in [chapter 07](07_pipeline.qmd) therefore fit **n = 87**, not 96. This tutorial keeps all 96 samples through normalization ([chapter 02](02_normalization.qmd)), probe filtering ([chapter 03](03_probe_filtering.qmd)), cell-composition estimation ([chapter 04](04_cell_composition.qmd)), and batch diagnosis ([chapter 05](05_batch_effects.qmd)) on purpose. None of those steps uses PTSD status: deconvolution and PCA describe the arrays, not the association model, so there is no reason to discard usable arrays before the step that actually needs the exposure. The count changes exactly once, at the point the exposure enters. The distinction worth internalising: **samples dropped for quality** are samples you do not trust, and **samples dropped for missing phenotype** are samples you trust but cannot model. Report them separately. Collapsing them into one "n = 87 after QC" sentence — the most common way this gets written up — hides the fact that nine usable arrays were lost to an incomplete metadata table, which is a very different problem from nine failed arrays. ::: Because the nine are a loss from missingness rather than quality, it is worth a moment to ask whether they are **missing at random** or missing in a way that correlates with something. If the samples without PTSD status were, say, all male or all from one chip, dropping them would distort the design: ```{r} #| label: missing-pattern no_ptsd <- is.na(ss$ptsd) | ss$ptsd == "" cat("samples with no PTSD status:", sum(no_ptsd), "\n") cat("sex:", paste(names(table(ss$sex[no_ptsd])), table(ss$sex[no_ptsd]), collapse = " / "), "\n") cat("age range:", paste(range(ss$age[no_ptsd]), collapse = "-"), "\n") cat("distinct chips:", length(unique(ss$slide[no_ptsd])), "of", length(unique(ss$slide)), "\n") ``` Four female and five male, ages spanning the full cohort range, spread across eight of the twelve chips — the missingness looks incidental rather than structured. That is reassuring, but note that it is a *check you can only do for observed variables*; it cannot rule out the samples differing on something unrecorded. ### Variables absent for every sample For blood methylation the most consequential gap is usually **smoking**. Smoking produces some of the largest and most reproducible methylation signals in blood [@joehanes2016smoking], so an unmeasured smoking variable is unmeasured confounding in every model you fit. The Grady deposit gives PTSD, age, sex, race, and childhood trauma — but no smoking variable, which is common in public deposits rather than unusual. That absence is worth noting here because it shapes decisions in the next two chapters. The full set of options, and the one this tutorial takes, is laid out in [the EWAS chapter](06_ewas.qmd). ::: {.callout-tip} ## Inventory your covariates before you normalize Write down the model you intend to fit, then check each term against the metadata you actually have. For each term ask two questions: *is it present at all?* and *how many samples is it missing for?* The first tells you what you cannot adjust for; the second tells you your real analysis n. Note both, and decide *then* how you will handle them — including whether a term is worth the samples it costs. ::: ## References {.unnumbered} ::: {#refs} :::