Probe filtering

Not every one of the ~865,000 EPIC probes should be used in downstream analyses. Some probes need to be removed because they were designed for assessing quality control, as we did in QC. Other probes need to be removed because they have been found to be problematic in some way: cross-hybridization to more than one place in the genome; proximity to a common SNP (so what looks like a methylation difference is really a genotype difference); mapping to repetitive regions. These are properties of the probe, so these problematic probes are filtered out from everyone. This notebook explains the categories and applies Wanding Zhou’s curated masks (Zhou et al. 2017) for Illumina methylation arrays.

What this chapter needs

This chapter starts from data/02_funnorm_grs.rds, written at the end of normalization. If you did not run that chapter, fetch the published checkpoint instead — note it is a 1.25 GB download:

Terminal

./get_data.sh C_normalized

Code

library(minfi)       # getAnnotation() for probe coordinates, detectionP() for the p-value matrix
library(data.table)  # fread() for the mask tables
library(ggplot2)     # the ancestry-contrast figure

# The functional-normalized GenomicRatioSet written at the end of
# [normalization](02_normalization.qmd): 96 samples, beta and M-values for every
# probe that survived normalization. Nothing is filtered yet -- that is this
# chapter's job.
grs <- readRDS("data/02_funnorm_grs.rds")
grs

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

1. Why some probes are unreliable by design

There are three main types of problematic probes and each has been characterized in dedicated methods papers:

Cross-reactive / non-specific probes. A 50-mer probe sequence is not guaranteed to be unique in the genome. Some probes hybridize to multiple loci, so their signal is a blend from several genomic locations. Chen et al. (Chen et al. 2013) cataloged cross-reactive probes on the 450K; McCartney et al. (McCartney et al. 2016) extended this to the EPIC array. Zhou’s MASK_mapping flag captures probes that do not map uniquely.

SNP-affected probes. If a common SNP falls in the probe body — especially at the interrogated CpG or the single-base extension position — then the assay may reflect genotype, and not methylation. Two individuals with different alleles will show a large “methylation difference” that is really a difference in DNA sequence. This is a dangerous and often overlooked category for an EWAS, because it produces reproducible, highly significant, and completely spurious associations whenever the SNP allele frequency differs between your comparison groups. Zhou provides SNP masks recommendations built from 1000 Genomes allele frequencies. If you are studying a population that is not well represented by one of the 1000 Genomes Project, it is recommended to find or curate a list of population-specific SNPs and exclude any CpGs within the last 5bp of the CpG probe. If there is no appropriate 1000 Genomes population and you don’t have population/study-specific allele frequencies, you can use the general population SNP mask and note this as a limitation of the study. If you are performing a study on a multi-ancestry cohort or don’t know the ancestry of your cohort, then the general population masking is appropriate.

Repetitive / structurally masked probes. Probes overlapping RepeatMasker regions, copy-number-variable segments, or requiring an extension-base substitution have also been flagged historically — as MASK_rmsk15, MASK_sub25_copy and MASK_typeINextBaseSwitch in the legacy manifest files from Zhou’s annotation database. The current v8.1 mask retains only the extension-base category (as M_1baseSwitchSNP_*), and the repeat and copy-number terms are now handled by the MASK_mapping criteria (based on a mapping quality score).

The latest CpG masking recommendations from Zhou (EPIC.hg38.mask.cm, v8.1) include 24 separate mask terms. One of those terms is M_general, which combines the recommended masking for probes with problematic mapping (the probes that should be masked independent of population) with CpGs within 5bp of common SNPs (MAF ≥ 5%) for the general population. If you do have a specific population, the remaining terms give population-specific masking for the 5 super-populations from the 1000 Genomes. The population-specific masking is a bit more stringent, as it uses MAF ≥ 1% — which is what the _1pt suffix on those term names records.

The terms are compositional: rather than a single pre-baked per-population column, you select the mapping terms plus the SNP terms for each ancestry you need and take their union. The full set of 24 terms for EPIC is:

mapping / specificity — M_mapping, M_nonuniq
SNP at the interrogated CpG, per population — M_SNP_{AFR,EUR,EAS,SAS,AMR}_1pt, plus M_SNPcommon_1pt and M_SNPcommon_5pt for the general population
SNP at the extension base — M_2extBase_SNP_* (same population suffixes)
Type I extension-base switch — M_1baseSwitchSNP_*
the recommended general default — M_general

If you need more specific sub-population masking, there is an archived version of the recommended population masking (EPIC.hg38.manifest.pop.tsv.gz) that includes the 5 super-populations and 26 sub-populations. If using the archived population masking file, the MASK_general_<POP> column includes the recommended problematic probes as well as the population-specific probes. The MASK_snp5_<POP> columns only flag the population-specific masking.

Why use Zhou’s masks specifically

The Zhou lab re-mapped every Infinium probe against the genome from scratch and integrated SNP information across 1000 Genomes populations (Zhou et al. 2017). The masks are the annotation layer underneath the sesame package (Zhou et al. 2018) and are the current community standard.

Getting the mask

There are two routes to the same mask. Take the first one unless you have a reason not to.

Route 1 — the ready-made TSV (recommended)

The eight mask codes this chapter uses are shipped with the tutorial data as a plain tab-delimited table, so you can filter probes with fread() and nothing else. No new software to install:

Where the mask table comes from

data/EPIC.hg38.mask.v81.8code.tsv.gz is tracked in the repository, so a fresh clone already has it — 2.9 MB, no download needed. It is also inside the D_filtered data tier, which is where it lived before; if you have already fetched that tier you have two copies of the same file and can ignore this note.

It is a derived table, not a Zhou lab release: Route 2 below shows the yame commands that produce it from the published .cm mask, so you can rebuild it for another array or another mask version.

Code

# One row per probe on the array, one 0/1 column per mask code, plus the
# row-wise OR of the eight in `mask_v81_8code`.
mask <- fread("data/EPIC.hg38.mask.v81.8code.tsv.gz")

dim(mask)                            # 866,553 probes x 10 columns

[1] 866553     10

Code

sum(mask$mask_v81_8code)             # 90,832 probes flagged by the union

[1] 90832

Code

# the probes to drop
drop8 <- mask[mask_v81_8code == 1L, Probe_ID]

The table has one row per probe on the array and one 0/1 column per mask code, plus a mask_v81_8code column that is the row-wise OR of the eight — which is exactly the filter this chapter applies. Because the individual code columns are kept, you can also build a different combination: to mask on African-ancestry SNP flags only, OR the three *_AFR_1pt columns and ignore the rest.

This file was generated from the official .cm release with YAME (below) and verified to reproduce all eight per-code counts and the 90,832-probe union exactly. It is a convenience layer, not a separate source of truth.

Route 2 — from the official release with YAME

Use this if you want a mask term the TSV does not carry (it ships 8 of the 24), a different array, or a newer mask version than the one pinned here.

The v8.1 mask ships as EPIC.hg38.mask.cm — a bit-packed CX-format file rather than a table. It stores one bit per probe per mask term over the array’s canonical probe ordering, which is why the whole thing is only ~360 kB where the equivalent TSV is 2.8 MB. Because the file carries no probe names (row i means whatever row i of the ordering means), you fetch the ordering alongside it.

Read it with YAME, the Zhou lab’s encoder for these formats:

# install: conda install -c zhou-lab -c conda-forge yame
# (bioconda's linux build is pinned at 1.8, which predates `yame fetch` --
#  use the zhou-lab channel or build from source)

yame fetch EPIC/EPIC.hg38.mask.cm EPIC/EPIC.ordering.tsv.gz

yame info    EPIC.hg38.mask.cm   # 24 mask terms over 866,553 probes
yame summary EPIC.hg38.mask.cm   # how many probes each term flags
yame unpack -a EPIC.hg38.mask.cm > mask_bits.tsv   # all 24 terms, one column each

yame summary is the quickest way to see what a term includes you before you commit to it. Turning the bits into a probe list is then a row-wise OR across the columns you chose, indexed against the ordering:

RStudio Console

library(data.table)
# `yame unpack` writes the term names to a sidecar .idx, one per line, in the
# same order as the columns of mask_bits.tsv -- so read the column, not the header.
terms  <- fread("EPIC.hg38.mask.cm.idx", header = FALSE)$V1
probes <- fread("EPIC.ordering.tsv.gz")$Probe_ID          # canonical probe ordering
bits   <- fread("mask_bits.tsv", header = FALSE); setnames(bits, terms)

# Choose the terms to mask on. These eight are the ones this chapter uses, and
# the same eight the ready-made TSV carries.
mask_cols <- c("M_mapping", "M_nonuniq",
               "M_SNP_AFR_1pt", "M_1baseSwitchSNP_AFR_1pt", "M_2extBase_SNP_AFR_1pt",
               "M_SNP_EUR_1pt", "M_1baseSwitchSNP_EUR_1pt", "M_2extBase_SNP_EUR_1pt")

drop <- probes[Reduce(`|`, lapply(mask_cols, function(m) bits[[m]] == 1L))]

The per-population SNP counts tabulated in §2 need terms the ready-made TSV does not carry, so they come from this route. This is the one derivation in the chapter that cannot run at render time — it needs the yame binary and the .cm release fetched above, and neither ships with the repository — so it runs once and its five-row result is committed:

Code

# eval: false because this needs the yame CLI and the files `yame fetch` /
# `yame unpack -a` produce, which are not in the repository. Run it once after
# the shell block above, then commit data/03_snp_by_pop.tsv.
terms  <- fread("EPIC.hg38.mask.cm.idx", header = FALSE)$V1
bits   <- fread("mask_bits.tsv", header = FALSE); setnames(bits, terms)

# M_SNP_<POP>_1pt flags a common SNP (MAF >= 1%) at the interrogated CpG in that
# 1000 Genomes super-population. Counted over the whole array, not over this
# study's probes, because it is a property of the platform.
pops <- c("AFR", "AMR", "SAS", "EUR", "EAS")
snp_by_pop <- vapply(pops,
                     function(p) sum(bits[[paste0("M_SNP_", p, "_1pt")]] == 1L),
                     integer(1))

# Five rows of two columns: small and rectangular, so it is committed to the
# repository and every reader gets the §2 table whether or not they have yame.
write.table(data.frame(population = pops, n_probes = snp_by_pop),
            "data/03_snp_by_pop.tsv", sep = "\t", row.names = FALSE, quote = FALSE)

The archived manifest is still there when you need it

EPIC.hg38.manifest.pop.tsv.gz remains available and is the right choice when you need sub-population resolution: it carries MASK_general_<POP> and MASK_snp5_<POP> columns for the 5 super-populations and 26 sub-populations (GWD, YRI, ASW, CEU, FIN, …), where v8.1’s .cm stops at the 5 super-populations. Note also that the two generations do not mask the same categories — see the comparison at the end of this section.

The funnel below needs all three probe lists, so the two study-specific filters are computed here and explained in §4.

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

# The detection-p filter is per PROBE, so it needs the full probe x sample
# p-value matrix. The QC chapter kept only per-sample summaries, so the matrix is
# recomputed from the raw two-channel object here: about 40 seconds and ~4 GB for
# 96 arrays, which is the most expensive step in this chapter.
rg   <- readRDS("data/01_RGset.rds")
detP <- detectionP(rg)

# Fraction of samples in which each probe failed the p > 0.01 cut, over the
# probes that survived normalization -- those are the only ones the funnel can
# drop. A probe unreliable in more than 10% of samples goes.
frac_fail <- rowMeans(detP[rownames(grs), ] > 0.01)
drop_detp <- rownames(grs)[frac_fail > 0.10]

# Nothing below needs the raw object or the p-value matrix, and together they are
# the memory high-water mark of the chapter, so release them now.
rm(rg, detP); invisible(gc())
length(drop_detp)

[1] 1928

Code

# The eight-code union from Route 1, restricted to the probes actually present
# after normalization.
drop_mask <- intersect(rownames(grs), drop8)

# Sex-chromosome probes, from the manifest coordinates minfi attached to the
# object -- removed for an autosomal EWAS, see §4.
drop_sex <- rownames(grs)[getAnnotation(grs)$chr %in% c("chrX", "chrY")]

# A probe is dropped if ANY of the three lists flags it, so the three counts
# below overlap and do not sum to the total dropped.
keep <- setdiff(rownames(grs), Reduce(union, list(drop_mask, drop_detp, drop_sex)))

# The same filter with the EUR-only mask instead of the AFR + EUR union, for the
# ancestry contrast in §3. Detection-p and sex filters are held fixed.
eur_cols  <- c("M_mapping", "M_nonuniq",
               "M_SNP_EUR_1pt", "M_1baseSwitchSNP_EUR_1pt", "M_2extBase_SNP_EUR_1pt")
drop_eur8 <- mask$Probe_ID[Reduce(`|`, lapply(eur_cols, function(m) mask[[m]] == 1L))]
keep_eur  <- setdiff(rownames(grs),
                     Reduce(union, list(intersect(rownames(grs), drop_eur8),
                                        drop_detp, drop_sex)))

# Per-code probe sets, kept so §2 can ask how two codes overlap rather than just
# how many probes each one flags.
mask_cols <- c("M_mapping", "M_nonuniq",
               "M_SNP_AFR_1pt", "M_1baseSwitchSNP_AFR_1pt", "M_2extBase_SNP_AFR_1pt",
               "M_SNP_EUR_1pt", "M_1baseSwitchSNP_EUR_1pt", "M_2extBase_SNP_EUR_1pt")
per_code_sets <- lapply(mask_cols, function(m) mask$Probe_ID[mask[[m]] == 1L])
names(per_code_sets) <- mask_cols
mp <- list(per_code_sets = per_code_sets)

# Per-population SNP counts: derived by the `snp-by-pop` chunk above, read back
# here because that chunk cannot run at render time.
sp <- fread("data/03_snp_by_pop.tsv")

# Every count the funnel, the tables and the prose below quote, in one place.
funl <- list(
  n0               = nrow(grs),                        # probes entering the funnel
  drop_mask        = length(drop_mask),
  drop_detp        = length(drop_detp),
  drop_sex         = length(drop_sex),
  retained         = length(keep),                     # EWAS-ready probes
  retained_eur     = length(keep_eur),                 # under a EUR-only mask
  afr_not_eur_kept = length(setdiff(keep_eur, keep)),  # what EUR-only would leave in
  union8           = sum(mask$mask_v81_8code),         # array-wide, before any study filter
  snp_by_pop       = setNames(sp$n_probes, sp$population)
)

Broken down, the probe removal includes:

Code

funnel <- data.frame(
  step = c("Start (normalized probes)",
           "− Zhou v8.1 mask (8 codes, AFR + EUR)",
           "− Detection-p failures (>10% of samples)",
           "− Sex chromosomes (autosomal EWAS)",
           "= Retained (EWAS-ready)"),
  probes_dropped = c(NA, funl$drop_mask, funl$drop_detp, funl$drop_sex, NA),
  probes_remaining = c(funl$n0, NA, NA, NA, funl$retained)
)
knitr::kable(funnel, caption = "Probe-filtering funnel on the 96-sample subset (AFR + EUR aware).")

Probe-filtering funnel on the 96-sample subset (AFR + EUR aware).
step	probes_dropped	probes_remaining
Start (normalized probes)	NA	865859
− Zhou v8.1 mask (8 codes, AFR + EUR)	90832	NA
− Detection-p failures (>10% of samples)	1928	NA
− Sex chromosomes (autosomal EWAS)	19627	NA
= Retained (EWAS-ready)	NA	756273

Code

cat(sprintf("Retained %d of %d probes (%.1f%%).\n",
            funl$retained, funl$n0, 100 * funl$retained / funl$n0))

Retained 756273 of 865859 probes (87.3%).

2. The categories, quantified on EPIC

The SNP mask is not the same for every study population, because a “common SNP” depends on which population you ask about. The number of SNP-affected probes varies substantially by ancestry:

Code

snp <- data.frame(
  population  = c("AFR (African)", "AMR (Admixed American)", "SAS (South Asian)",
                  "EUR (European)", "EAS (East Asian)"),
  snp_masked  = funl$snp_by_pop[c("AFR", "AMR", "SAS", "EUR", "EAS")]
)
knitr::kable(snp, row.names = FALSE,
             caption = "Probes flagged by M_SNP_<POP>_1pt (common SNP at the interrogated CpG, MAF >= 1%) by 1000 Genomes super-population. African-ancestry samples have roughly twice as many SNP-affected probes as European.")

Probes flagged by M_SNP__1pt (common SNP at the interrogated CpG, MAF >= 1%) by 1000 Genomes super-population. African-ancestry samples have roughly twice as many SNP-affected probes as European.
population	snp_masked
AFR (African)	35420
AMR (Admixed American)	18629
SAS (South Asian)	17780
EUR (European)	17076
EAS (East Asian)	15281

African-ancestry populations carry the most common genetic variation, so more probes overlap a common SNP. Concretely, in these masks 35,420 probes are SNP-flagged for AFR versus 17,076 for EUR:

Code

# Both sets come from the per-code columns of the shipped mask table, so the
# comparison is between two mask terms over the same probe universe.
afr <- mp$per_code_sets$M_SNP_AFR_1pt
eur <- mp$per_code_sets$M_SNP_EUR_1pt
cat("masked in AFR but NOT EUR:", length(setdiff(afr, eur)), "\n")

masked in AFR but NOT EUR: 24435

Code

cat("masked in EUR but NOT AFR:", length(setdiff(eur, afr)), "\n")

masked in EUR but NOT AFR: 6091

Code

cat("masked in both:           ", length(intersect(afr, eur)), "\n")

masked in both:            10985

While there is some overlap between the AFR and EUR recommended SNP masks, there are distinct maskings for each population.

3. Ancestry matters

The Grady Trauma Project subset is ~95% African American. If we filtered with a general or only European-population mask (a common default), we would leave in probes that are SNP-affected in African-ancestry individuals. As African American individuals tend to have an admixture of both European and African Ancestry, it is important to include masking recommendations for both.

Code

contrast <- data.frame(
  mask = c("AFR + EUR union\n(used here)", "EUR only\n(wrong here)"),
  retained = c(funl$retained, funl$retained_eur)
)
ggplot(contrast, aes(mask, retained, fill = mask)) +
  geom_col(width = 0.6) +
  geom_text(aes(label = format(retained, big.mark = ",")), vjust = -0.4, size = 3.5) +
  scale_fill_manual(values = c("#1A6B75",      # teal: the mask this chapter applies
                               "#8C3A4A"),     # plum: the mask that would be wrong here
                    guide = "none") +
  labs(x = NULL, y = "probes retained",
       title = "Ancestry-aware probe filtering") +
  expand_limits(y = max(contrast$retained) * 1.06)

Probes retained under the AFR+EUR union mask used here versus a EUR-only mask. Filtering this African-American cohort on European SNP flags alone would keep tens of thousands of probes that are SNP-affected in African-ancestry individuals.

Code

cat("Probes kept by a EUR-only mask but dropped by the AFR + EUR union:",
    funl$afr_not_eur_kept, "\n")

Probes kept by a EUR-only mask but dropped by the AFR + EUR union: 36553

The union is what the code below applies, and it is the operational meaning of “include masking recommendations for both”.

So, if we are using the latest recommended masking file (v8.1), the columns we would want to include in the masking would be M_mapping, M_nonuniq, M_SNP_AFR_1pt,M_1baseSwitchSNP_AFR_1pt, M_2extBase_SNP_AFR_1pt, M_SNP_EUR_1pt, M_1baseSwitchSNP_EUR_1pt, and M_2extBase_SNP_EUR_1pt. If we were using the archived manifest population file, we would use the MASK_general_AFR and MASK_snp5_EUR.

What changing mask generation costs you

It is worth seeing the size of the version effect directly, because it runs the opposite way to intuition. The eight v8.1 codes cover two ancestries; the legacy MASK_general_AFR covers one. Yet the legacy mask is the more aggressive:

Code

gen <- data.frame(
  mask = c("Legacy MASK_general_AFR", "v8.1, 8 codes (AFR + EUR)"),
  # Legacy counts are quoted from the archived manifest; the v8.1 figures are the
  # ones derived above, so only the mask generation differs between the rows.
  probes_masked   = c(124132, funl$union8),
  probes_retained = c(724282, funl$retained)
)
knitr::kable(gen, row.names = FALSE,
             caption = "Same cohort, same detection-p and sex filters; only the mask generation differs.")

Same cohort, same detection-p and sex filters; only the mask generation differs.
mask	probes_masked	probes_retained
Legacy MASK_general_AFR	124132	724282
v8.1, 8 codes (AFR + EUR)	90832	756273

The difference is the repeat and copy-number terms described in §2: the legacy MASK_general_<POP> bundles them in, and v8.1 has no equivalent. So moving to v8.1 retains 31,991 more probes here — more sites tested, but also the repeat-overlapping probes no longer excluded. Neither choice is automatically right; what matters is that you state which mask generation and which terms you used, since the two are not interchangeable.

The teaching point

Probe filtering is not a fixed recipe — the SNP component depends on your study population. For a cohort of known, relatively homogeneous ancestry, use that population’s mask. For an admixed or mixed-ancestry cohort, the conservative choice is to take the union of the relevant population masks (drop a probe if it is SNP-affected in any represented ancestry). For a multi-ancestry cohort or cohort of unknown ancestry, the general population mask is the most appropriate.

Smoking-associated probes

Everything above masks probes that are unreliable by design. There is a separate, optional reason to drop probes: smoking is one of the largest environmental drivers of blood methylation, and if your dataset has no smoking variable — as the Grady deposit does not — one of the available strategies is to exclude CpGs known from the literature to be smoking-associated, at this filtering stage rather than after the fact. That is a decision about confounding, not about probe quality, so it does not belong in the funnel above. This tutorial does exclude smoking-associated CpGs, but at a later step as these probes are needed to create a proxy smoking variable. See the smoking section in the EWAS chapter for the full set of options.

4. Detection-p and sex-chromosome filtering

Two study-specific filters complete the funnel:

Detection-p failures. From QC we have the per-probe detection p-values. A probe that fails (p > 0.01) in more than ~10% of samples is unreliable in this dataset and is dropped. Here that removed 1928 probes.
Sex chromosomes. For an autosomal EWAS, X- and Y-chromosome probes are also removed (19627 probes): X-inactivation and sex differences would otherwise dominate, and they need sex-stratified handling if analyzed. They are removed from the autosomal analysis, not deleted from the study — a sex-specific analysis would use them deliberately.

The exact filtering code (run once, checkpointed):

Code

# `keep` is the set union of the three drop lists, taken out of the normalized
# probe set: the mask from Route 1 (`drop_mask`), the detection-p failures
# (`drop_detp`) and the sex chromosomes (`drop_sex`), all built above. Order does
# not matter -- a probe flagged by two filters is dropped once.
grs_filt <- grs[keep, ]

# Subsetting a GenomicRatioSet carries the beta/M matrices, the genome
# coordinates and the sample metadata along with it, so this object is
# analysis-ready as it stands.
grs_filt

class: GenomicRatioSet 
dim: 756273 96 
metadata(0):
assays(2): Beta CN
rownames(756273): cg26928153 cg16269199 ... cg19565306 cg22012975
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

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

Probe filtering is one of the few preprocessing steps that works identically whether your data came from IDATs or from a processed β-matrix — because the masks are keyed on probe IDs, not intensities. The workflow:

Load the matrix (see the normalization non-IDAT section) so you have a β-matrix with cg/ch probe IDs as rownames.

Apply the same masks by ID. Read Zhou’s mask/manifest for the correct array (450K vs EPIC vs EPICv2 — they have different probe sets), pick the mask for your study population, and subset:

RStudio Console

# `drop8` is the eight-code union read in Route 1 above. It works unchanged on
# a matrix, because the mask is keyed on probe IDs and nothing else.
beta_filt <- beta[!rownames(beta) %in% drop8, ]

# archived manifest, if you need sub-population resolution instead
library(data.table)
pop  <- fread("EPIC.hg38.manifest.pop.tsv.gz")
drop <- pop$probeID[pop$MASK_general_ASW]           # e.g. African-ancestry SW US
beta_filt <- beta[!rownames(beta) %in% drop, ]

You lose the detection-p filter unless the submitter deposited detection p-values alongside the matrix (some do, as a separate supplementary file — use it if present). Without it, rely on the mask plus removing probes with excessive NA values.
Sex chromosomes filter the same way, by intersecting with the manifest’s chromosome annotation. It is possible sex chromosome CpGs will have already been removed from a pre-processed dataset.

Match the mask to the array and the genome build

Zhou publishes masks per platform (450K / EPIC / EPICv2) and per genome build (hg19 / hg38). Use the manifest that matches the array the data were generated on; the probe IDs mostly overlap between 450K and EPIC but the mask contents differ. If you are combining a 450K series with an EPIC series, filter each on its own manifest, then intersect on the shared probe IDs.

We save the filtered, EWAS-ready object for cell-composition estimation next:

Code

# 1. The EWAS-ready GenomicRatioSet. An S4 object with four parallel matrices and
#    a GRanges, so RDS rather than a flat file. Read by cell-composition
#    estimation ([chapter 04](04_cell_composition.qmd)) and the batch-effect
#    diagnostics ([chapter 05](05_batch_effects.qmd)).
saveRDS(grs_filt, "data/03_grs_filtered.rds")

# 2. The funnel counts, so a later chapter or a methods section can quote them
#    without re-running the filter.
saveRDS(funl, "data/03_filter_funnel.rds")

# 3. The per-code probe sets behind §2. A list of eight character vectors of
#    unequal length, which a flat file would not hold.
saveRDS(mp, "data/03_mask_pieces.rds")

# 4. The funnel as a flat table as well: five rows, so it is committed to the
#    repository and readable straight from GitHub by anyone who cannot hold the
#    GenomicRatioSet in memory.
write.table(funnel, "data/03_filter_funnel.tsv",
            sep = "\t", row.names = FALSE, quote = FALSE, na = "")

If the detection-p step will not fit on your machine

Recomputing detectionP() holds the raw two-channel object and a probe × sample p-value matrix in memory at once — roughly 4 GB on top of the normalized object. If that is more than you have, fetch the filtered checkpoint this chapter writes and start from it instead:

Terminal

./get_data.sh D_filtered

RStudio Console

grs_filt <- readRDS("data/03_grs_filtered.rds")

What you give up is the funnel: every count above is computed from the three probe lists, not read from a file, so none of it appears if you skip the chapter’s code.

The result: 756,273 autosomal probes on 96 samples, filtered under the Zhou v8.1 mask for mapping quality, SNP masking appropriate to an admixed African-American cohort (AFR and EUR terms together), detection reliability, and sex chromosomes.

Nearly every probe that survives here is tested in the EWAS. Two small deductions happen later, both in chapter 05:

the 20 CpGs used to build the smoking proxy are excluded, because regressing those probes on a score derived from themselves is circular; and
two probes (cg17759086, cg01801182) have β pinned at 0 or 1 in the 87-sample analysis subset, so their M-values are ±∞ and they cannot be fitted.

So 756,273 probes pass filtering and 756,251 are actually tested.

Next: Cell composition — whole blood is a mixture of cell types with very different methylomes; we estimate the proportions and validate against the study’s deposited values.

References

Chen, Yi-an, Mathieu Lemire, Sanaa Choufani, et al. 2013. “Discovery of Cross-Reactive Probes and Polymorphic CpGs in the Illumina Infinium HumanMethylation450 Microarray.” Epigenetics 8 (2): 203–9. https://doi.org/10.4161/epi.23470.

McCartney, Daniel L, Rosie M Walker, Stewart W Morris, Andrew M McIntosh, David J Porteous, and Kathryn L Evans. 2016. “Identification of Polymorphic and Off-Target Probe Binding Sites on the Illumina Infinium MethylationEPIC BeadChip.” Genomics Data 9: 22–24. https://doi.org/10.1016/j.gdata.2016.05.012.

Zhou, Wanding, Peter W Laird, and Hui Shen. 2017. “Comprehensive Characterization, Annotation and Innovative Use of Infinium DNA Methylation BeadChip Probes.” Nucleic Acids Research 45 (4): e22. https://doi.org/10.1093/nar/gkw967.

Zhou, Wanding, Timothy J Triche Jr, Peter W Laird, and Hui Shen. 2018. “SeSAMe: Reducing Artifactual Detection of DNA Methylation by Infinium BeadChips in Genomic Deletions.” Nucleic Acids Research 46 (20): e123. https://doi.org/10.1093/nar/gky691.

--- title: "Probe filtering" --- ```{r} #| label: setup #| include: false source("_setup.R") ``` Not every one of the ~865,000 EPIC probes should be used in downstream analyses. Some probes need to be removed because they were designed for assessing quality control, as we did in [QC](01_qc.qmd). Other probes need to be removed because they have been found to be problematic in some way: cross-hybridization to more than one place in the genome; proximity to a common SNP (so what looks like a methylation difference is really a genotype difference); mapping to repetitive regions. These are properties of the **probe**, so these problematic probes are filtered out from everyone. This notebook explains the categories and applies Wanding Zhou's curated masks [@zhou2017comprehensive] for Illumina methylation arrays. ::: {.callout-note} ## What this chapter needs This chapter starts from `data/02_funnorm_grs.rds`, written at the end of [normalization](02_normalization.qmd). If you did not run that chapter, fetch the published checkpoint instead — note it is a 1.25 GB download: ```{.bash filename="Terminal"} ./get_data.sh C_normalized ``` ::: ```{r} #| label: load library(minfi) # getAnnotation() for probe coordinates, detectionP() for the p-value matrix library(data.table) # fread() for the mask tables library(ggplot2) # the ancestry-contrast figure # The functional-normalized GenomicRatioSet written at the end of # [normalization](02_normalization.qmd): 96 samples, beta and M-values for every # probe that survived normalization. Nothing is filtered yet -- that is this # chapter's job. grs <- readRDS("data/02_funnorm_grs.rds") grs ``` ## 1. Why some probes are unreliable by design {#sec-why} There are three main types of problematic probes and each has been characterized in dedicated methods papers: **Cross-reactive / non-specific probes.** A 50-mer probe sequence is not guaranteed to be unique in the genome. Some probes hybridize to multiple loci, so their signal is a blend from several genomic locations. Chen et al. [@chen2013crossreactive] cataloged cross-reactive probes on the 450K; McCartney et al. [@mccartney2016identification] extended this to the EPIC array. Zhou's `MASK_mapping` flag captures probes that do not map uniquely. **SNP-affected probes.** If a common SNP falls in the probe body — especially at the interrogated CpG or the single-base extension position — then the assay may reflect **genotype**, and not methylation. Two individuals with different alleles will show a large "methylation difference" that is really a difference in DNA sequence. This is a dangerous and often overlooked category for an EWAS, because it produces reproducible, highly significant, and completely spurious associations whenever the SNP allele frequency differs between your comparison groups. Zhou provides SNP masks recommendations built from 1000 Genomes allele frequencies. If you are studying a population that is not well represented by one of the 1000 Genomes Project, it is recommended to find or curate a list of population-specific SNPs and exclude any CpGs within the last 5bp of the CpG probe. If there is no appropriate 1000 Genomes population and you don't have population/study-specific allele frequencies, you can use the general population SNP mask and note this as a limitation of the study. If you are performing a study on a multi-ancestry cohort or don't know the ancestry of your cohort, then the general population masking is appropriate. **Repetitive / structurally masked probes.** Probes overlapping RepeatMasker regions, copy-number-variable segments, or requiring an extension-base substitution have also been flagged historically — as `MASK_rmsk15`, `MASK_sub25_copy` and `MASK_typeINextBaseSwitch` in the legacy manifest files from Zhou's annotation database. The current v8.1 mask retains only the extension-base category (as `M_1baseSwitchSNP_*`), and the repeat and copy-number terms are now handled by the `MASK_mapping` criteria (based on a mapping quality score). The latest CpG masking recommendations from Zhou (`EPIC.hg38.mask.cm`, v8.1) include 24 separate mask terms. One of those terms is `M_general`, which combines the recommended masking for probes with problematic mapping (the probes that should be masked independent of population) with CpGs within 5bp of common SNPs (MAF ≥ 5%) for the general population. If you do have a specific population, the remaining terms give population-specific masking for the 5 super-populations from the 1000 Genomes. The population-specific masking is a bit more stringent, as it uses MAF ≥ 1% — which is what the `_1pt` suffix on those term names records. The terms are compositional: rather than a single pre-baked per-population column, you select the mapping terms plus the SNP terms for each ancestry you need and take their union. The full set of 24 terms for EPIC is: - **mapping / specificity** — `M_mapping`, `M_nonuniq` - **SNP at the interrogated CpG**, per population — `M_SNP_{AFR,EUR,EAS,SAS,AMR}_1pt`, plus `M_SNPcommon_1pt` and `M_SNPcommon_5pt` for the general population - **SNP at the extension base** — `M_2extBase_SNP_*` (same population suffixes) - **Type I extension-base switch** — `M_1baseSwitchSNP_*` - **the recommended general default** — `M_general` If you need more specific sub-population masking, there is an archived version of the recommended population masking (`EPIC.hg38.manifest.pop.tsv.gz`) that includes the 5 super-populations and 26 sub-populations. If using the archived population masking file, the `MASK_general_<POP>` column includes the recommended problematic probes as well as the population-specific probes. The `MASK_snp5_<POP>` columns only flag the population-specific masking. ::: {.callout-note} ## Why use Zhou's masks specifically The Zhou lab re-mapped every Infinium probe against the genome from scratch and integrated SNP information across 1000 Genomes populations [@zhou2017comprehensive]. The masks are the annotation layer underneath the `sesame` package [@zhou2018sesame] and are the current community standard. ::: ### Getting the mask There are two routes to the same mask. Take the first one unless you have a reason not to. #### Route 1 — the ready-made TSV (recommended) The eight mask codes this chapter uses are shipped with the tutorial data as a plain tab-delimited table, so you can filter probes with `fread()` and nothing else. No new software to install: ::: {.callout-note} ## Where the mask table comes from `data/EPIC.hg38.mask.v81.8code.tsv.gz` is tracked in the repository, so a fresh clone already has it — 2.9 MB, no download needed. It is also inside the `D_filtered` data tier, which is where it lived before; if you have already fetched that tier you have two copies of the same file and can ignore this note. It is a derived table, not a Zhou lab release: [Route 2](#route-2) below shows the `yame` commands that produce it from the published `.cm` mask, so you can rebuild it for another array or another mask version. ::: ```{r} #| label: mask-read # One row per probe on the array, one 0/1 column per mask code, plus the # row-wise OR of the eight in `mask_v81_8code`. mask <- fread("data/EPIC.hg38.mask.v81.8code.tsv.gz") dim(mask) # 866,553 probes x 10 columns sum(mask$mask_v81_8code) # 90,832 probes flagged by the union # the probes to drop drop8 <- mask[mask_v81_8code == 1L, Probe_ID] ``` The table has one row per probe on the array and one 0/1 column per mask code, plus a `mask_v81_8code` column that is the row-wise OR of the eight — which is exactly the filter this chapter applies. Because the individual code columns are kept, you can also build a different combination: to mask on African-ancestry SNP flags only, OR the three `*_AFR_1pt` columns and ignore the rest. This file was generated from the official `.cm` release with YAME (below) and verified to reproduce all eight per-code counts and the 90,832-probe union exactly. It is a convenience layer, not a separate source of truth. #### Route 2 — from the official release with YAME {#route-2} Use this if you want a mask term the TSV does not carry (it ships 8 of the 24), a different array, or a newer mask version than the one pinned here. The v8.1 mask ships as `EPIC.hg38.mask.cm` — a bit-packed **CX-format** file rather than a table. It stores one bit per probe per mask term over the array's canonical probe ordering, which is why the whole thing is only ~360 kB where the equivalent TSV is 2.8 MB. Because the file carries no probe names (row *i* means whatever row *i* of the ordering means), you fetch the ordering alongside it. Read it with [YAME](https://github.com/zhou-lab/YAME), the Zhou lab's encoder for these formats: ```bash # install: conda install -c zhou-lab -c conda-forge yame # (bioconda's linux build is pinned at 1.8, which predates `yame fetch` -- # use the zhou-lab channel or build from source) yame fetch EPIC/EPIC.hg38.mask.cm EPIC/EPIC.ordering.tsv.gz yame info EPIC.hg38.mask.cm # 24 mask terms over 866,553 probes yame summary EPIC.hg38.mask.cm # how many probes each term flags yame unpack -a EPIC.hg38.mask.cm > mask_bits.tsv # all 24 terms, one column each ``` `yame summary` is the quickest way to see what a term includes you before you commit to it. Turning the bits into a probe list is then a row-wise OR across the columns you chose, indexed against the ordering: ```{.r filename="RStudio Console"} library(data.table) # `yame unpack` writes the term names to a sidecar .idx, one per line, in the # same order as the columns of mask_bits.tsv -- so read the column, not the header. terms <- fread("EPIC.hg38.mask.cm.idx", header = FALSE)$V1 probes <- fread("EPIC.ordering.tsv.gz")$Probe_ID # canonical probe ordering bits <- fread("mask_bits.tsv", header = FALSE); setnames(bits, terms) # Choose the terms to mask on. These eight are the ones this chapter uses, and # the same eight the ready-made TSV carries. mask_cols <- c("M_mapping", "M_nonuniq", "M_SNP_AFR_1pt", "M_1baseSwitchSNP_AFR_1pt", "M_2extBase_SNP_AFR_1pt", "M_SNP_EUR_1pt", "M_1baseSwitchSNP_EUR_1pt", "M_2extBase_SNP_EUR_1pt") drop <- probes[Reduce(`|`, lapply(mask_cols, function(m) bits[[m]] == 1L))] ``` The per-population SNP counts tabulated in §2 need terms the ready-made TSV does not carry, so they come from this route. This is the one derivation in the chapter that cannot run at render time — it needs the `yame` binary and the `.cm` release fetched above, and neither ships with the repository — so it runs once and its five-row result is committed: ```{r} #| label: snp-by-pop #| eval: false # eval: false because this needs the yame CLI and the files `yame fetch` / # `yame unpack -a` produce, which are not in the repository. Run it once after # the shell block above, then commit data/03_snp_by_pop.tsv. terms <- fread("EPIC.hg38.mask.cm.idx", header = FALSE)$V1 bits <- fread("mask_bits.tsv", header = FALSE); setnames(bits, terms) # M_SNP_<POP>_1pt flags a common SNP (MAF >= 1%) at the interrogated CpG in that # 1000 Genomes super-population. Counted over the whole array, not over this # study's probes, because it is a property of the platform. pops <- c("AFR", "AMR", "SAS", "EUR", "EAS") snp_by_pop <- vapply(pops, function(p) sum(bits[[paste0("M_SNP_", p, "_1pt")]] == 1L), integer(1)) # Five rows of two columns: small and rectangular, so it is committed to the # repository and every reader gets the §2 table whether or not they have yame. write.table(data.frame(population = pops, n_probes = snp_by_pop), "data/03_snp_by_pop.tsv", sep = "\t", row.names = FALSE, quote = FALSE) ``` ::: {.callout-note} ## The archived manifest is still there when you need it `EPIC.hg38.manifest.pop.tsv.gz` remains available and is the right choice when you need **sub-population** resolution: it carries `MASK_general_<POP>` and `MASK_snp5_<POP>` columns for the 5 super-populations *and* 26 sub-populations (GWD, YRI, ASW, CEU, FIN, …), where v8.1's `.cm` stops at the 5 super-populations. Note also that the two generations do not mask the same categories — see the comparison at the end of this section. ::: The funnel below needs all three probe lists, so the two study-specific filters are computed here and explained in [§4](#sec-detp-sex). ::: {.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: probe-detp # The detection-p filter is per PROBE, so it needs the full probe x sample # p-value matrix. The QC chapter kept only per-sample summaries, so the matrix is # recomputed from the raw two-channel object here: about 40 seconds and ~4 GB for # 96 arrays, which is the most expensive step in this chapter. rg <- readRDS("data/01_RGset.rds") detP <- detectionP(rg) # Fraction of samples in which each probe failed the p > 0.01 cut, over the # probes that survived normalization -- those are the only ones the funnel can # drop. A probe unreliable in more than 10% of samples goes. frac_fail <- rowMeans(detP[rownames(grs), ] > 0.01) drop_detp <- rownames(grs)[frac_fail > 0.10] # Nothing below needs the raw object or the p-value matrix, and together they are # the memory high-water mark of the chapter, so release them now. rm(rg, detP); invisible(gc()) length(drop_detp) ``` ```{r} #| label: probe-filters # The eight-code union from Route 1, restricted to the probes actually present # after normalization. drop_mask <- intersect(rownames(grs), drop8) # Sex-chromosome probes, from the manifest coordinates minfi attached to the # object -- removed for an autosomal EWAS, see §4. drop_sex <- rownames(grs)[getAnnotation(grs)$chr %in% c("chrX", "chrY")] # A probe is dropped if ANY of the three lists flags it, so the three counts # below overlap and do not sum to the total dropped. keep <- setdiff(rownames(grs), Reduce(union, list(drop_mask, drop_detp, drop_sex))) # The same filter with the EUR-only mask instead of the AFR + EUR union, for the # ancestry contrast in §3. Detection-p and sex filters are held fixed. eur_cols <- c("M_mapping", "M_nonuniq", "M_SNP_EUR_1pt", "M_1baseSwitchSNP_EUR_1pt", "M_2extBase_SNP_EUR_1pt") drop_eur8 <- mask$Probe_ID[Reduce(`|`, lapply(eur_cols, function(m) mask[[m]] == 1L))] keep_eur <- setdiff(rownames(grs), Reduce(union, list(intersect(rownames(grs), drop_eur8), drop_detp, drop_sex))) # Per-code probe sets, kept so §2 can ask how two codes overlap rather than just # how many probes each one flags. mask_cols <- c("M_mapping", "M_nonuniq", "M_SNP_AFR_1pt", "M_1baseSwitchSNP_AFR_1pt", "M_2extBase_SNP_AFR_1pt", "M_SNP_EUR_1pt", "M_1baseSwitchSNP_EUR_1pt", "M_2extBase_SNP_EUR_1pt") per_code_sets <- lapply(mask_cols, function(m) mask$Probe_ID[mask[[m]] == 1L]) names(per_code_sets) <- mask_cols mp <- list(per_code_sets = per_code_sets) # Per-population SNP counts: derived by the `snp-by-pop` chunk above, read back # here because that chunk cannot run at render time. sp <- fread("data/03_snp_by_pop.tsv") # Every count the funnel, the tables and the prose below quote, in one place. funl <- list( n0 = nrow(grs), # probes entering the funnel drop_mask = length(drop_mask), drop_detp = length(drop_detp), drop_sex = length(drop_sex), retained = length(keep), # EWAS-ready probes retained_eur = length(keep_eur), # under a EUR-only mask afr_not_eur_kept = length(setdiff(keep_eur, keep)), # what EUR-only would leave in union8 = sum(mask$mask_v81_8code), # array-wide, before any study filter snp_by_pop = setNames(sp$n_probes, sp$population) ) ``` Broken down, the probe removal includes: ```{r} #| label: funnel-table funnel <- data.frame( step = c("Start (normalized probes)", "− Zhou v8.1 mask (8 codes, AFR + EUR)", "− Detection-p failures (>10% of samples)", "− Sex chromosomes (autosomal EWAS)", "= Retained (EWAS-ready)"), probes_dropped = c(NA, funl$drop_mask, funl$drop_detp, funl$drop_sex, NA), probes_remaining = c(funl$n0, NA, NA, NA, funl$retained) ) knitr::kable(funnel, caption = "Probe-filtering funnel on the 96-sample subset (AFR + EUR aware).") cat(sprintf("Retained %d of %d probes (%.1f%%).\n", funl$retained, funl$n0, 100 * funl$retained / funl$n0)) ``` ## 2. The categories, quantified on EPIC {#sec-categories} The SNP mask is **not the same for every study population**, because a "common SNP" depends on which population you ask about. The number of SNP-affected probes varies substantially by ancestry: ```{r} #| label: pop-snp-table snp <- data.frame( population = c("AFR (African)", "AMR (Admixed American)", "SAS (South Asian)", "EUR (European)", "EAS (East Asian)"), snp_masked = funl$snp_by_pop[c("AFR", "AMR", "SAS", "EUR", "EAS")] ) knitr::kable(snp, row.names = FALSE, caption = "Probes flagged by M_SNP_<POP>_1pt (common SNP at the interrogated CpG, MAF >= 1%) by 1000 Genomes super-population. African-ancestry samples have roughly twice as many SNP-affected probes as European.") ``` African-ancestry populations carry the most common genetic variation, so more probes overlap a common SNP. Concretely, in these masks **`r format(funl$snp_by_pop[["AFR"]], big.mark = ",")`** probes are SNP-flagged for AFR versus **`r format(funl$snp_by_pop[["EUR"]], big.mark = ",")`** for EUR: ```{r} #| label: afr-eur-overlap # Both sets come from the per-code columns of the shipped mask table, so the # comparison is between two mask terms over the same probe universe. afr <- mp$per_code_sets$M_SNP_AFR_1pt eur <- mp$per_code_sets$M_SNP_EUR_1pt cat("masked in AFR but NOT EUR:", length(setdiff(afr, eur)), "\n") cat("masked in EUR but NOT AFR:", length(setdiff(eur, afr)), "\n") cat("masked in both: ", length(intersect(afr, eur)), "\n") ``` While there is some overlap between the AFR and EUR recommended SNP masks, there are distinct maskings for each population. ## 3. Ancestry matters {#sec-ancestry} The Grady Trauma Project subset is **~95% African American**. If we filtered with a general or only European-population mask (a common default), we would leave in probes that are SNP-affected in African-ancestry individuals. As African American individuals tend to have an admixture of both European and African Ancestry, it is important to include masking recommendations for both. ```{r} #| label: ancestry-contrast #| fig-cap: "Probes retained under the AFR+EUR union mask used here versus a EUR-only mask. Filtering this African-American cohort on European SNP flags alone would keep tens of thousands of probes that are SNP-affected in African-ancestry individuals." #| fig-width: 6.5 #| fig-height: 3 contrast <- data.frame( mask = c("AFR + EUR union\n(used here)", "EUR only\n(wrong here)"), retained = c(funl$retained, funl$retained_eur) ) ggplot(contrast, aes(mask, retained, fill = mask)) + geom_col(width = 0.6) + geom_text(aes(label = format(retained, big.mark = ",")), vjust = -0.4, size = 3.5) + scale_fill_manual(values = c("#1A6B75", # teal: the mask this chapter applies "#8C3A4A"), # plum: the mask that would be wrong here guide = "none") + labs(x = NULL, y = "probes retained", title = "Ancestry-aware probe filtering") + expand_limits(y = max(contrast$retained) * 1.06) ``` ```{r} #| label: ancestry-diff cat("Probes kept by a EUR-only mask but dropped by the AFR + EUR union:", funl$afr_not_eur_kept, "\n") ``` The union is what the code below applies, and it is the operational meaning of "include masking recommendations for both". So, if we are using the latest recommended masking file (v8.1), the columns we would want to include in the masking would be `M_mapping`, `M_nonuniq`, `M_SNP_AFR_1pt`,`M_1baseSwitchSNP_AFR_1pt`, `M_2extBase_SNP_AFR_1pt`, `M_SNP_EUR_1pt`, `M_1baseSwitchSNP_EUR_1pt`, and `M_2extBase_SNP_EUR_1pt`. If we were using the archived manifest population file, we would use the `MASK_general_AFR` and `MASK_snp5_EUR`. ### What changing mask generation costs you It is worth seeing the size of the version effect directly, because it runs the opposite way to intuition. The eight v8.1 codes cover **two** ancestries; the legacy `MASK_general_AFR` covers one. Yet the legacy mask is the more aggressive: ```{r} #| label: mask-generation-contrast gen <- data.frame( mask = c("Legacy MASK_general_AFR", "v8.1, 8 codes (AFR + EUR)"), # Legacy counts are quoted from the archived manifest; the v8.1 figures are the # ones derived above, so only the mask generation differs between the rows. probes_masked = c(124132, funl$union8), probes_retained = c(724282, funl$retained) ) knitr::kable(gen, row.names = FALSE, caption = "Same cohort, same detection-p and sex filters; only the mask generation differs.") ``` The difference is the repeat and copy-number terms described in §2: the legacy `MASK_general_<POP>` bundles them in, and v8.1 has no equivalent. So moving to v8.1 retains **`r format(funl$retained - 724282, big.mark = ",")`** more probes here — more sites tested, but also the repeat-overlapping probes no longer excluded. Neither choice is automatically right; what matters is that you state which mask generation and which terms you used, since the two are not interchangeable. ::: {.callout-important} ## The teaching point Probe filtering is not a fixed recipe — the SNP component depends on your study population. For a cohort of known, relatively homogeneous ancestry, use that population's mask. For an admixed or mixed-ancestry cohort, the conservative choice is to take the **union** of the relevant population masks (drop a probe if it is SNP-affected in *any* represented ancestry). For a multi-ancestry cohort or cohort of unknown ancestry, the general population mask is the most appropriate. ::: ::: {.callout-note} ## Smoking-associated probes Everything above masks probes that are unreliable *by design*. There is a separate, optional reason to drop probes: smoking is one of the largest environmental drivers of blood methylation, and if your dataset has no smoking variable — as the Grady deposit does not — one of the available strategies is to exclude CpGs known from the literature to be smoking-associated, at this filtering stage rather than after the fact. That is a decision about confounding, not about probe quality, so it does not belong in the funnel above. This tutorial does exclude smoking-associated CpGs, but at a later step as these probes are needed to create a proxy smoking variable. See [the smoking section in the EWAS chapter](06_ewas.qmd) for the full set of options. ::: ## 4. Detection-p and sex-chromosome filtering {#sec-detp-sex} Two study-specific filters complete the funnel: - **Detection-p failures.** From [QC](01_qc.qmd) we have the per-probe detection p-values. A probe that fails (p > 0.01) in more than ~10% of samples is unreliable *in this dataset* and is dropped. Here that removed `r funl$drop_detp` probes. - **Sex chromosomes.** For an autosomal EWAS, X- and Y-chromosome probes are also removed (`r funl$drop_sex` probes): X-inactivation and sex differences would otherwise dominate, and they need sex-stratified handling if analyzed. They are *removed from the autosomal analysis*, not deleted from the study — a sex-specific analysis would use them deliberately. The exact filtering code (run once, checkpointed): ```{r} #| label: apply-filter # `keep` is the set union of the three drop lists, taken out of the normalized # probe set: the mask from Route 1 (`drop_mask`), the detection-p failures # (`drop_detp`) and the sex chromosomes (`drop_sex`), all built above. Order does # not matter -- a probe flagged by two filters is dropped once. grs_filt <- grs[keep, ] # Subsetting a GenomicRatioSet carries the beta/M matrices, the genome # coordinates and the sample metadata along with it, so this object is # analysis-ready as it stands. grs_filt ``` ## 5. When you don't have IDATs: filtering a processed matrix {#sec-no-idats} Probe filtering is one of the few preprocessing steps that works **identically** whether your data came from IDATs or from a processed β-matrix — because the masks are keyed on **probe IDs**, not intensities. The workflow: 1. **Load the matrix** (see the [normalization](02_normalization.qmd) non-IDAT section) so you have a β-matrix with `cg`/`ch` probe IDs as rownames. 2. **Apply the same masks by ID.** Read Zhou's mask/manifest for the correct array (450K vs EPIC vs EPICv2 — they have different probe sets), pick the mask for your study population, and subset: ```{.r filename="RStudio Console"} # `drop8` is the eight-code union read in Route 1 above. It works unchanged on # a matrix, because the mask is keyed on probe IDs and nothing else. beta_filt <- beta[!rownames(beta) %in% drop8, ] # archived manifest, if you need sub-population resolution instead library(data.table) pop <- fread("EPIC.hg38.manifest.pop.tsv.gz") drop <- pop$probeID[pop$MASK_general_ASW] # e.g. African-ancestry SW US beta_filt <- beta[!rownames(beta) %in% drop, ] ``` 3. **You lose the detection-p filter** unless the submitter deposited detection p-values alongside the matrix (some do, as a separate supplementary file — use it if present). Without it, rely on the mask plus removing probes with excessive `NA` values. 4. **Sex chromosomes** filter the same way, by intersecting with the manifest's chromosome annotation. It is possible sex chromosome CpGs will have already been removed from a pre-processed dataset. ::: {.callout-tip} ## Match the mask to the array *and* the genome build Zhou publishes masks per platform (450K / EPIC / EPICv2) and per genome build (hg19 / hg38). Use the manifest that matches the array the data were generated on; the probe IDs mostly overlap between 450K and EPIC but the mask *contents* differ. If you are combining a 450K series with an EPIC series, filter each on its own manifest, then intersect on the shared probe IDs. ::: We save the filtered, EWAS-ready object for cell-composition estimation next: ```{r} #| label: save-filtered # 1. The EWAS-ready GenomicRatioSet. An S4 object with four parallel matrices and # a GRanges, so RDS rather than a flat file. Read by cell-composition # estimation ([chapter 04](04_cell_composition.qmd)) and the batch-effect # diagnostics ([chapter 05](05_batch_effects.qmd)). saveRDS(grs_filt, "data/03_grs_filtered.rds") # 2. The funnel counts, so a later chapter or a methods section can quote them # without re-running the filter. saveRDS(funl, "data/03_filter_funnel.rds") # 3. The per-code probe sets behind §2. A list of eight character vectors of # unequal length, which a flat file would not hold. saveRDS(mp, "data/03_mask_pieces.rds") # 4. The funnel as a flat table as well: five rows, so it is committed to the # repository and readable straight from GitHub by anyone who cannot hold the # GenomicRatioSet in memory. write.table(funnel, "data/03_filter_funnel.tsv", sep = "\t", row.names = FALSE, quote = FALSE, na = "") ``` ::: {.callout-note collapse="true"} ## If the detection-p step will not fit on your machine Recomputing `detectionP()` holds the raw two-channel object and a probe × sample p-value matrix in memory at once — roughly 4 GB on top of the normalized object. If that is more than you have, fetch the filtered checkpoint this chapter writes and start from it instead: ```{.bash filename="Terminal"} ./get_data.sh D_filtered ``` ```{.r filename="RStudio Console"} grs_filt <- readRDS("data/03_grs_filtered.rds") ``` What you give up is the funnel: every count above is computed from the three probe lists, not read from a file, so none of it appears if you skip the chapter's code. ::: The result: **`r format(funl$retained, big.mark=",")`** autosomal probes on 96 samples, filtered under the Zhou v8.1 mask for mapping quality, SNP masking appropriate to an admixed African-American cohort (AFR and EUR terms together), detection reliability, and sex chromosomes. Nearly every probe that survives here is tested in the EWAS. Two small deductions happen later, both in [chapter 05](05_batch_effects.qmd): - the **20 CpGs used to build the smoking proxy** are excluded, because regressing those probes on a score derived from themselves is circular; and - **two probes** (`cg17759086`, `cg01801182`) have β pinned at 0 or 1 in the 87-sample analysis subset, so their M-values are ±∞ and they cannot be fitted. So `r format(funl$retained, big.mark=",")` probes pass filtering and **756,251** are actually tested. --- **Next:** [Cell composition](04_cell_composition.qmd) — whole blood is a mixture of cell types with very different methylomes; we estimate the proportions and validate against the study's deposited values. ## References {.unnumbered} ::: {#refs} :::