Rohit GuptaA DIY guide to creating a brain tumor classifier
Brain tumor classification from DNA methylation is one of the clearer success stories of machine learning in pathology. The Capper et al. (2018) classifier showed that a random forest trained on Illumina methylation arrays can sort central nervous system tumors into a large set of molecularly defined classes, and that this often refines or corrects the histological diagnosis. The reference data are public, the arrays are standardized, and the core model is a random forest. That combination makes it a great project for learning how to build a classifier properly.
The word “properly” is doing a lot of work there. With tens of thousands of features, a few thousand samples, and dozens of classes, it is very easy to build a model that looks excellent in cross-validation and falls apart on the next hospital’s samples. Most of this post is about avoiding that.
What follows is a tutorial, not a clinical tool. Anything used for patient care needs validation, quality management and regulatory work far beyond a blog post.
What the input looks like
The Illumina Infinium HumanMethylation450 (450k) and MethylationEPIC (EPIC) arrays measure methylation at several hundred thousand CpG sites. Each probe gives a methylated intensity (M) and an unmethylated intensity (U), and the usual summary is the beta value:
beta = M / (M + U + offset)
Beta is roughly the fraction of copies methylated at that CpG, between 0 and 1. The raw data come as pairs of IDAT files per sample (one for the red channel, one for green), named after the slide barcode and array position, for example 200123450001_R01C01_Grn.idat.
For a classifier, you want a samples by probes matrix of beta values after normalization and probe filtering. I do that part in R with either minfi or sesame.
Getting beta values with minfi
# Schematic: needs your own IDATs or a large GEO download, so it isn't run here
library(minfi)
library(IlluminaHumanMethylation450kanno.ilmn12.hg19) # EPIC needs the EPIC manifest/annotation packages
# Placeholder path: a directory containing *_Grn.idat / *_Red.idat pairs
rgSet <- read.metharray.exp(base = "path/to/idats", recursive = TRUE)
# Detection p-values: probe signal versus background
detP <- detectionP(rgSet)
# Drop samples where many probes fail
bad_sample <- colMeans(detP > 0.01) > 0.05
rgSet <- rgSet[, !bad_sample]
detP <- detP[, !bad_sample]
# Normalize (Noob background and dye-bias correction), map to the genome
mSet <- preprocessNoob(rgSet)
grSet <- ratioConvert(mapToGenome(mSet), what = "beta")
# Probes that failed detection in more than 1% of samples
detP <- detP[rownames(grSet), ]
grSet <- grSet[rowMeans(detP > 0.01) < 0.01, ]
# Probes with a SNP at the CpG or single-base extension site
grSet <- dropLociWithSnps(grSet, snps = c("SBE", "CpG"), maf = 0)
# Sex chromosomes and non-CpG ("ch.") probes
ann <- getAnnotation(grSet)
grSet <- grSet[!(ann$chr %in% c("chrX", "chrY")), ]
grSet <- grSet[grepl("^cg", rownames(grSet)), ]
# Cross-reactive probes: placeholder file, one probe ID per line,
# taken from a published list (Chen et al. 2013 for 450k, Pidsley et al. 2016 for EPIC)
cross_reactive <- readLines("path/to/cross_reactive_probes.txt")
grSet <- grSet[!(rownames(grSet) %in% cross_reactive), ]
beta <- getBeta(grSet)
data.table::fwrite(
data.table::data.table(probe = rownames(beta), beta),
"betas.tsv.gz", sep = "\t"
)
Each filter has a reason. Sex chromosome probes mostly encode the patient’s sex (and X inactivation), which you do not want the model to learn as a proxy for tumor type, especially if sex is unevenly distributed across classes. SNP-affected probes measure genotype as much as methylation. Cross-reactive probes hybridize to more than one genomic location, so their signal is ambiguous. Probes that fail detection are noise.
Or with sesame
# Schematic: needs your own IDATs or a large GEO download, so it isn't run here
library(sesame)
sesameDataCache() # one-time download of annotation data
betas <- openSesame("path/to/idats") # probes x samples matrix
openSesame applies its own preprocessing, including pOOBAH detection masking, and returns NA for masked probes. Sesame also ships its own probe masks based on Zhou et al. (2017), which cover much of the SNP and cross-reactivity filtering above. Either pipeline is fine. Pick one and use it for every sample, including future ones. Mixing preprocessing pipelines between training data and new cases is a quiet way to introduce a batch effect.
The reference cohort
The Capper et al. reference cohort is deposited in GEO as GSE90496, about 2,800 samples profiled on the 450k array, spanning 82 CNS tumor methylation classes plus several control tissue classes. The companion validation cohort is GSE109379. Look at the supplementary files on the GEO series page: there is a raw IDAT archive, which I recommend over any processed matrix because it lets you run the same preprocessing you will apply to new samples.
The class labels live in the sample characteristics. In R, GEOquery gets them:
# Schematic: needs your own IDATs or a large GEO download, so it isn't run here
library(GEOquery)
gse <- getGEO("GSE90496", GSEMatrix = TRUE)
meta <- pData(gse[[1]])
# Find the column holding the methylation class, then check it
grep("class", colnames(meta), value = TRUE, ignore.case = TRUE)
Inspect the column names rather than trusting my memory of them, and check that the sample identifiers in the metadata match the IDAT file names or matrix columns you end up with. Save a two-column table (sample ID, class) as labels.tsv.
Building the training set
Load the beta matrix and labels in Python:
import numpy as np
import pandas as pd
# Placeholder paths. betas.tsv.gz is probes x samples, written by the R step above.
betas = pd.read_csv("betas.tsv.gz", sep="\t", index_col=0).T.astype(np.float32)
labels = pd.read_csv("labels.tsv", sep="\t", index_col=0).iloc[:, 0]
labels = labels.loc[betas.index] # same sample order; fails loudly if IDs don't match
A full 450k matrix for a few thousand samples is several gigabytes in float32, so keep an eye on memory.
Class imbalance and minimum class size
Methylation classes are very unequal in size. Common entities like glioblastoma subclasses have many reference samples; rare entities have a handful. Three decisions follow:
- Set a minimum class size. With 5-fold outer and 3-fold inner cross-validation, a class needs enough samples to appear in every training and validation split. I use at least 10 as a floor for this kind of tutorial. Classes below the floor are either merged into a sensible parent (if one exists biologically) or excluded and documented.
- Stratify every split by class.
- Weight or resample. In scikit-learn,
class_weight="balanced_subsample"reweights classes inside each bootstrap sample. The alternative is to downsample to equal class sizes per tree, which is whatBalancedRandomForestClassifierin imbalanced-learn does.
from sklearn.model_selection import train_test_split
counts = labels.value_counts()
keep_classes = counts[counts >= 10].index
mask = labels.isin(keep_classes)
X = betas.loc[mask]
y = labels.loc[mask]
print(f"{X.shape[0]} samples, {y.nunique()} classes, "
f"dropped {counts.size - keep_classes.size} small classes")
A held-out test set
Before anything else, split off a test set and do not look at it until the very end:
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.2, stratify=y, random_state=2023
)
Nested cross-validation (below) already gives an honest estimate of generalization, so why also hold out a test set? Because you are going to make choices beyond hyperparameters: which probes to filter, which classes to merge, which calibration method to use, where to put the threshold. Every one of those decisions is made by looking at cross-validation results. The test set is the one number that no decision touched. If you look at it, change something and look again, it becomes a second validation set and stops being a test set.
Even better than a random split is a truly external set: samples from a different center or a later time period, such as the validation cohort, or your own institution’s cases if you have appropriate approval.
Feature selection, and why it lives inside the loop
You cannot reasonably feed 400,000 probes into a random forest with a few thousand samples, and most probes are near-constant across CNS tumors anyway. The standard first step is to keep the most variable CpGs, for example the top 5,000 to 20,000 by variance or standard deviation.
Here is the trap. If you compute probe variances on the full dataset, select the top probes, and then run cross-validation on that reduced matrix, information from the validation folds has already shaped the features. With an unsupervised filter like variance, the leak is modest. With any supervised selection (t-tests between classes, random forest importance, anything that uses labels) the leak can be severe, and on high-dimensional data it can produce impressive accuracy from pure noise. Varma and Simon (2006) is the classic demonstration of how far optimistic these estimates can get.
The rule: anything that learns from data is part of the model, and is refit inside each training fold. That includes imputation, feature selection, scaling, and batch correction if you do it.
scikit-learn makes this straightforward with a Pipeline. There is no built-in “top k by variance” selector, but SelectKBest accepts any scoring function, and a scoring function that ignores y gives you exactly that:
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.feature_selection import SelectKBest
from sklearn.ensemble import RandomForestClassifier
def variance_score(X, y=None):
"""Score each probe by its variance. Ignores y on purpose."""
return np.nanvar(X, axis=0)
pipe = Pipeline([
("impute", SimpleImputer(strategy="median")),
("select", SelectKBest(score_func=variance_score, k=10_000)),
("rf", RandomForestClassifier(
n_estimators=1000,
class_weight="balanced_subsample",
n_jobs=-1,
random_state=0,
)),
])
When this pipeline is fit on a training fold, the imputer learns medians from that fold, the selector picks the top variable probes from that fold, and the forest trains on those probes. The validation fold is then transformed with the fold’s medians and the fold’s probe list. Nothing about the validation samples influences the model.
A practical note: running the variance filter on all 400,000 probes inside every fold is slow and memory hungry. A reasonable compromise is a coarse, label-free pre-filter computed on X_train only (say, the top 50,000 probes by variance), followed by the in-pipeline selection of k from those. That pre-filter never sees the test set, so the test estimate stays clean. It does see all the cross-validation folds, so the cross-validation estimate gets a small optimistic bias, which is a trade I’ll accept for an unsupervised filter. I would not make that trade for a supervised one.
train_var = X_train.var(axis=0)
prefilter = train_var.sort_values(ascending=False).index[:50_000]
X_train = X_train[prefilter]
X_test = X_test[prefilter]
Beta versus M-values does not matter much here. Random forest splits are invariant to monotone transformations of individual features, so the trees would be identical. The variance ranking is different between beta and M-values, though, so the selected probes would change. I use betas for this.
Nested cross-validation
Hyperparameter tuning is also a form of learning from data. If you tune on the same folds you report, the reported score is the best of many tries and is biased upward. Nested cross-validation separates the two jobs:
X_train, y_train (test set already locked away)
│
├── OUTER loop: StratifiedKFold(5) → estimates generalization
│ │
│ ├── outer fold 1: train on folds 2-5, validate on fold 1
│ │ │
│ │ └── INNER loop: StratifiedKFold(3) on folds 2-5 only
│ │ ├── for each candidate (k, n_estimators, max_features, min_samples_leaf):
│ │ │ fit pipeline on 2/3, score on 1/3, repeat 3x, average
│ │ └── pick best candidate, refit it on all of folds 2-5
│ │
│ │ score the refit model on outer fold 1 → balanced accuracy #1
│ │
│ ├── outer fold 2: same procedure, fresh inner search → #2
│ ├── ...
│ └── outer fold 5 → #5
│
├── report mean and spread of #1-#5
│
└── final model: run the inner search once on all of X_train, refit best,
then evaluate ONCE on X_test
The outer score estimates how well the whole procedure (including tuning and feature selection) generalizes. It does not belong to any single set of hyperparameters, and different outer folds may well pick different ones. That is fine and expected.
The inner search
from sklearn.model_selection import StratifiedKFold, RandomizedSearchCV
param_distributions = {
"select__k": [2_000, 5_000, 10_000, 20_000],
"rf__n_estimators": [500, 1000],
"rf__max_features": ["sqrt", "log2", 0.01],
"rf__min_samples_leaf": [1, 2, 5],
}
inner_cv = StratifiedKFold(n_splits=3, shuffle=True, random_state=1)
outer_cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=2)
search = RandomizedSearchCV(
pipe,
param_distributions=param_distributions,
n_iter=20,
scoring="balanced_accuracy",
cv=inner_cv,
refit=True,
random_state=3,
n_jobs=1, # the forest already uses all cores
)
A few comments on the grid. n_estimators is not really a hyperparameter you need to tune: more trees reduce variance and never make a random forest meaningfully worse, they only cost time. I include it to show the syntax; in practice I pick a large number and leave it. max_features controls how many probes each split considers and is usually the most influential setting. min_samples_leaf smooths the trees and the probability estimates. select__k tunes the feature filter itself, which is only possible because the filter is inside the pipeline.
The scoring metric is balanced accuracy, the average of per-class recall. With imbalanced classes, plain accuracy rewards a model that is good at the big classes and ignores the small ones.
The outer loop
The short version uses cross_validate, which treats the whole search object as the estimator:
from sklearn.model_selection import cross_validate
nested = cross_validate(
search, X_train, y_train,
cv=outer_cv,
scoring=["balanced_accuracy", "accuracy"],
return_estimator=True,
)
print(nested["test_balanced_accuracy"])
print([est.best_params_ for est in nested["estimator"]])
I usually write the outer loop explicitly so I can keep the out-of-fold predictions for a confusion matrix:
from sklearn.base import clone
from sklearn.metrics import balanced_accuracy_score
outer_scores = []
oof_pred = pd.Series(index=y_train.index, dtype=object)
for fold, (tr, va) in enumerate(outer_cv.split(X_train, y_train)):
X_tr, X_va = X_train.iloc[tr], X_train.iloc[va]
y_tr, y_va = y_train.iloc[tr], y_train.iloc[va]
fold_search = clone(search)
fold_search.fit(X_tr, y_tr)
pred = fold_search.predict(X_va)
oof_pred.iloc[va] = pred
score = balanced_accuracy_score(y_va, pred)
outer_scores.append(score)
print(f"fold {fold}: balanced accuracy {score:.3f} {fold_search.best_params_}")
print(f"nested CV balanced accuracy: "
f"{np.mean(outer_scores):.3f} +/- {np.std(outer_scores):.3f}")
This is computationally heavy: 5 outer folds times 20 candidates times 3 inner folds is 300 forest fits, plus refits. Start with a small n_iter and a smaller k range to check that everything runs, then scale up on a machine with plenty of cores and memory.
The final model
Once you are happy with the procedure, run it once more on all of the training data:
search.fit(X_train, y_train)
final_model = search.best_estimator_
print(search.best_params_)
Evaluation
The headline is balanced accuracy, but the information is in the per-class numbers. Evaluate on the test set once, at the end:
from sklearn.metrics import classification_report, confusion_matrix, recall_score
classes = final_model.classes_
test_pred = final_model.predict(X_test)
print("balanced accuracy:", balanced_accuracy_score(y_test, test_pred))
print(classification_report(y_test, test_pred, zero_division=0))
per_class_recall = pd.Series(
recall_score(y_test, test_pred, labels=classes, average=None, zero_division=0),
index=classes,
).sort_values()
print(per_class_recall.head(10)) # the classes the model struggles with
With dozens of classes, a full confusion matrix plot is hard to read. Listing the largest off-diagonal cells is more useful:
cm = pd.DataFrame(
confusion_matrix(y_test, test_pred, labels=classes),
index=pd.Index(classes, name="true"),
columns=pd.Index(classes, name="predicted"),
)
off_diag = cm.where(~np.eye(len(classes), dtype=bool)).stack()
print(off_diag[off_diag > 0].sort_values(ascending=False).head(15))
You can run the same code on oof_pred versus y_train for the cross-validation view. Expect most confusions to be between biologically related classes, such as subclasses within a tumor family. Those are much less worrying than confusions between unrelated entities. Capper et al. grouped related classes into families for exactly this reason, and summing calibrated probabilities over a family gives a family-level score. A confusion between unrelated classes is worth investigating sample by sample, because it often points to a labeling problem, a low-purity sample, or a technical artifact.
Calibration and the “no match” option
A random forest’s predict_proba returns the fraction of trees voting for each class (averaged over leaf class frequencies, to be exact). Those numbers are scores, not probabilities. With many related classes, votes get split between neighbors, so a correctly classified sample might get a top score of 0.4. A score of 0.4 for one sample and 0.4 for another do not necessarily mean the same level of confidence across classes.
Capper et al. addressed this by adding a calibration step: a multinomial logistic regression that takes the vector of random forest scores and maps it to calibrated class probabilities. The idea is that after calibration, among all samples given a probability of 0.9 for some class, about 90% truly belong to it. That is what lets you set a meaningful threshold.
The key detail is that the calibrator must be trained on out-of-fold forest scores. Scores for samples the forest was trained on are inflated, because each tree has seen most of them.
from sklearn.model_selection import GridSearchCV, cross_val_predict
from sklearn.linear_model import LogisticRegression
calib_cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=4)
# Out-of-fold random forest scores for every training sample
rf_scores_oof = cross_val_predict(
final_model, X_train, y_train, cv=calib_cv, method="predict_proba"
)
# Multinomial logistic regression on those scores, with the
# regularization strength chosen by cross-validated log loss
calibrator = GridSearchCV(
LogisticRegression(max_iter=5000),
param_grid={"C": np.logspace(-2, 3, 11)},
scoring="neg_log_loss",
cv=calib_cv,
)
calibrator.fit(rf_scores_oof, y_train)
print("chosen C:", calibrator.best_params_["C"])
assert list(calibrator.classes_) == list(final_model.classes_)
def predict_calibrated(model, calibrator, X_new):
rf_scores = model.predict_proba(X_new)
return pd.DataFrame(
calibrator.predict_proba(rf_scores),
index=X_new.index,
columns=calibrator.classes_,
)
cross_val_predict refits clones of final_model with its chosen hyperparameters, so the variance filter is refit per fold here too. Tune the regularization rather than leaving C at its default. The inputs are vote fractions between 0 and 1, so the coefficients need to be fairly large to produce confident probabilities, and with the default C=1.0 the penalty can flatten everything toward uniform. On simulated data I used to test this code, the default gave far less confident (and worse calibrated, by log loss) probabilities than the tuned value.
If you want something off the shelf, CalibratedClassifierCV does a similar job:
from sklearn.calibration import CalibratedClassifierCV
calibrated = CalibratedClassifierCV(
estimator=clone(final_model), method="sigmoid", cv=5
)
calibrated.fit(X_train, y_train)
For multiclass problems it calibrates each class one-vs-rest and then renormalizes, which is a little different from a single multinomial model over the whole score vector. method="isotonic" is more flexible but needs more data per class than most rare tumor classes have.
Choosing a threshold
Calibrated probabilities let you say “no match” when the top class is not convincing. Capper et al. used a cutoff of 0.9 on their calibrated scores, but that number belongs to their model and data. Choose yours by looking at the trade-off between how many samples get a call and how accurate those calls are, using out-of-fold data, never the test set:
# Out-of-fold calibrated probabilities: refit the calibrator (with its
# own C search) inside each fold of the out-of-fold forest scores
calib_oof = cross_val_predict(
clone(calibrator), rf_scores_oof, y_train,
cv=calib_cv, method="predict_proba",
)
top_prob = calib_oof.max(axis=1)
top_class = calibrator.classes_[calib_oof.argmax(axis=1)]
correct = top_class == y_train.to_numpy()
rows = []
for t in [0.3, 0.5, 0.7, 0.8, 0.9, 0.95]:
called = top_prob >= t
rows.append({
"threshold": t,
"fraction_called": called.mean(),
"accuracy_when_called": correct[called].mean() if called.any() else np.nan,
})
print(pd.DataFrame(rows))
Pick the threshold based on what an error costs in your setting, write it down, and then apply it once to the test set:
THRESHOLD = 0.9 # whatever you chose from the table above
test_probs = predict_calibrated(final_model, calibrator, X_test)
test_call = test_probs.idxmax(axis=1).where(
test_probs.max(axis=1) >= THRESHOLD, "no match"
)
called = test_call != "no match"
print(f"called {called.mean():.1%} of test samples")
print("accuracy among called:", (test_call[called] == y_test[called]).mean())
Report both numbers together. A classifier that is right 99% of the time on the 40% of cases it is willing to call is a very different tool from one that calls 95% of cases at slightly lower accuracy.
One more caveat: in this layout the calibrator and the threshold were chosen outside the nested loop, so the nested cross-validation score does not include them. The held-out test set is what covers the full pipeline.
Pitfalls that cost more than any hyperparameter
Batch effects. Arrays processed on different days, in different labs, or from different sample types (formalin-fixed paraffin-embedded versus frozen tissue) differ systematically. If a class was mostly profiled in one batch, the model can learn the batch. Capper et al. included an adjustment for material type in their pipeline. At minimum, check whether your top principal components track slide, center or material type, and when you validate, try to hold out whole batches or centers rather than random samples (StratifiedGroupKFold does this if you have a group label).
Tumor purity. A tumor sample is a mixture of tumor cells and normal brain, blood and immune cells. Low-purity samples drift toward the control tissue classes or get low scores everywhere. This is one reason the reference includes control tissue classes, and one reason the “no match” outcome is useful. If you have purity estimates, look at whether your errors concentrate in low-purity samples.
Array type differences. The reference is 450k; most new cases now run on EPIC or EPIC v2. These arrays share most but not all probes, and some probes behave differently across versions. Restrict training features to probes present (and passing QC) on every array you will use. In minfi, combineArrays can merge 450k and EPIC data onto their common probes. EPIC v2 probe IDs carry suffixes, so map them back to the original IDs carefully. Then check that EPIC samples of a known class score the way 450k samples do.
Label noise. Reference labels come from a mix of histology, molecular markers and methylation-based clustering. Some are wrong and some classes are defined partly by the method you are reproducing. Persistent cross-validation errors are worth reviewing case by case before you blame the model. The random forest is fairly tolerant of a little label noise; your evaluation is not.
Rare classes. A class with 10 samples gets maybe 2 in each test fold. Its recall estimate is coarse, and a single sample swings it a lot. Report per-class sample counts next to per-class metrics, and be careful about claims for the small ones.
Reusing the test set. Covered above, but it is the one I see most often. If you evaluated on the test set and then changed something, you need a new test set, or at least an honest note that the number is optimistic.
Wrong unit of independence. If the same patient contributes several samples (a primary and a recurrence, or two regions), keep them in the same fold. Otherwise the model can recognize the patient instead of the tumor class. Again, group-aware splitting handles this.
An R route
If you prefer to stay in R, the same structure works with ranger (fast, multithreaded) or randomForest. ranger(x = x, y = y, num.trees = 1000, mtry = 100, probability = TRUE) fits a probability forest, and the nested loop can be written with caret, tidymodels, or mlr3, or by hand with explicit fold indices. The same rules apply: feature selection inside the folds, calibration on out-of-fold scores, and the test set touched once.
Checklist
Before trusting any number from your classifier:
- Same preprocessing pipeline (minfi or sesame, same version and settings) for reference and new samples.
- Probe QC: detection p-values, sex chromosomes, SNP and cross-reactive probes removed; features restricted to probes shared by every array type you use.
- Classes below a minimum size merged or excluded, and documented.
- Test set split off first and used once.
- Imputation and feature selection inside a
Pipeline, refit per fold. - Hyperparameters tuned in an inner loop; generalization estimated in an outer loop.
- Balanced accuracy plus per-class recall and the largest confusions, with per-class sample counts.
- Calibrator trained on out-of-fold scores; threshold chosen without the test set; coverage reported alongside accuracy.
- Top principal components checked against batch, center, array type and material.
- Patients, and ideally centers, kept within a single fold.
The Capper et al. paper and its supplementary methods are worth reading in full once you have built your own version. Many of their design choices make a lot more sense after you have tripped over the problems they solve.
References
- Capper D, et al. DNA methylation-based classification of central nervous system tumours. Nature 555, 469–474 (2018).
- Breiman L. Random Forests. Machine Learning 45, 5–32 (2001).
- Varma S, Simon R. Bias in error estimation when using cross-validation for model selection. BMC Bioinformatics 7, 91 (2006).
- Aryee MJ, et al. Minfi: a flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays. Bioinformatics 30, 1363–1369 (2014).
- Zhou W, Triche TJ, Laird PW, Shen H. SeSAMe: reducing artifactual detection of DNA methylation by Infinium BeadChips in genomic deletions. Nucleic Acids Research 46, e123 (2018).
- Zhou W, Laird PW, Shen H. Comprehensive characterization, annotation and innovative use of Infinium DNA methylation BeadChip probes. Nucleic Acids Research 45, e22 (2017).
- Chen YA, et al. Discovery of cross-reactive probes and polymorphic CpGs in the Illumina Infinium HumanMethylation450 microarray. Epigenetics 8, 203–209 (2013).
- Pidsley R, et al. Critical evaluation of the Illumina MethylationEPIC BeadChip microarray for whole-genome DNA methylation profiling. Genome Biology 17, 208 (2016).