2026-08-14
This also happens to be the first step in the TCGA processing pipeline. The complete pipeline is available on GitHub as the EGC.tools project. Ten samples from the TCGA breast cancer (BRCA) project are included in the TCGAMethylation450k package, which should be installed for this step.
suppressPackageStartupMessages({
library(methylumi)
library(TCGAMethylation450k)
library(FDb.InfiniumMethylation.hg19)
})## read in 10 BRCA IDATs
idatPath <- system.file("extdata/idat", package = "TCGAMethylation450k")
mset450k <- methylumIDAT(getBarcodes(path = idatPath), idatPath = idatPath)
#> 0 HumanMethylation27 samples found
#> 10 HumanMethylation450 samples found
#> Attempting to extract protocolData() from list...
#> Determining chip type from IDAT protocolData...
sampleNames(mset450k) <- paste0("TCGA", seq_along(sampleNames(mset450k)))
show(mset450k)
#>
#> Object Information:
#> MethyLumiSet (storageMode: lockedEnvironment)
#> assayData: 485577 features, 10 samples
#> element names: betas, methylated, methylated.OOB, pvals, unmethylated, unmethylated.OOB
#> protocolData: none
#> phenoData
#> sampleNames: TCGA1 TCGA2 ... TCGA10 (10 total)
#> varLabels: barcode
#> varMetadata: labelDescription
#> featureData
#> featureNames: cg00000029 cg00000108 ... rs9839873 (485577 total)
#> fvarLabels: Probe_ID DESIGN COLOR_CHANNEL
#> fvarMetadata: labelDescription
#> experimentData: use 'experimentData(object)'
#> Annotation: IlluminaHumanMethylation450k
#> Major Operation History:
#> submitted finished
#> 1 2026-08-14 18:52:58.867331 2026-08-14 18:53:20.8319
#> 2 2026-08-14 18:53:21.827672 2026-08-14 18:53:23.290412
#> command
#> 1 methylumIDAT(barcodes = getBarcodes(path = idatPath), idatPath = idatPath)
#> 2 Subset of 485577 features.Note that the default is to collect opposite-channel fluorescence from Type I methylation probes (which are paired and designed to fluoresce in one channel) in the matrices methylated.OOB and unmethylated.OOB (OOB, as in out-of-band) for use in background correction and perhaps additional steps. This also allows a user to coerce the resulting object into minfi’s RGChannelSet if desired, since all of the signal information in the IDATs is thus retained.
library(ggplot2)
## for larger datasets, the by.type argument can be set to FALSE;
## positional effects will manifest as a wave-like pattern
p <- qc.probe.plot(mset450k, by.type = TRUE)
print(p)
After importing the data from IDATs, the next step is to background correct and dye bias equalize the data. The default for background correction is a normal-exponential model which uses the out-of-band intensities as control probes. Dye bias correction is performed by picking the least-biased sample and using it as a reference for red:green intensity ratio adjustments based on the normalization controls. Other approaches to preprocessing (as implemented in minfi and lumi) include various flavors of quantile normalization and smoothing spline fits.
After preprocessing, we can reduce the size of the resulting MethyLumiSet substantially by dropping the out-of-band intensities with stripOOB(). This frees up some memory, but precludes later coercion to an RGChannelSet.
mset450k.proc <- stripOOB(normalizeMethyLumiSet(methylumi.bgcorr(mset450k)))
#> Background mean & SD estimated from 178406 probes
#> Background mean & SD estimated from 92596 probes
#> Normalizing via Illumina controls...
#> Using sample number 5 as reference level...Now we compare the post-processing controls with those from the figure above.
p2 <- qc.probe.plot(mset450k.proc, by.type = TRUE)
print(p2)
Coercions are provided to and from various data structures in the lumi and minfi packages. Each provides various functionality and exhibits different design decisions. One may be more appropriate than the other for some needs. Preprocessing in methylumi retains SNP probes, which can identify label swaps, but is less efficient than preprocessing in minfi, and cannot use shinyMethyl.
Coercing to lumi (e.g. for lumiMethyN or similar):
suppressPackageStartupMessages(library(lumi))
#> No methods found in package 'RSQLite' for request: 'dbListFields' when loading 'lumi'
mset450k.lumi <- as(mset450k.proc, "MethyLumiM")
show(mset450k.lumi)
#> MethyLumiM (storageMode: lockedEnvironment)
#> assayData: 485577 features, 10 samples
#> element names: detection, exprs, methylated, unmethylated
#> protocolData: none
#> phenoData
#> sampleNames: TCGA1 TCGA2 ... TCGA10 (10 total)
#> varLabels: barcode mu.Cy3 ... offset.Cy5 (9 total)
#> varMetadata: labelDescription
#> featureData
#> featureNames: cg00000029 cg00000108 ... rs9839873 (485577 total)
#> fvarLabels: Probe_ID DESIGN COLOR_CHANNEL
#> fvarMetadata: labelDescription
#> experimentData: use 'experimentData(object)'
#> Annotation: IlluminaHumanMethylation450kCoercing back to a MethyLumiSet:
mset450k.andBack <- as(mset450k.lumi, "MethyLumiSet")
show(mset450k.andBack)
#>
#> Object Information:
#> MethyLumiSet (storageMode: lockedEnvironment)
#> assayData: 485577 features, 10 samples
#> element names: betas, methylated, pvals, unmethylated
#> protocolData: none
#> phenoData
#> sampleNames: TCGA1 TCGA2 ... TCGA10 (10 total)
#> varLabels: barcode mu.Cy3 ... offset.Cy5 (9 total)
#> varMetadata: labelDescription
#> featureData
#> featureNames: cg00000029 cg00000108 ... rs9839873 (485577 total)
#> fvarLabels: Probe_ID DESIGN COLOR_CHANNEL
#> fvarMetadata: labelDescription
#> experimentData: use 'experimentData(object)'
#> Annotation: IlluminaHumanMethylation450k
#> Major Operation History:
#> submitted finished
#> 1 2026-08-14 18:52:58.867331 2026-08-14 18:53:20.8319
#> 2 2026-08-14 18:53:21.827672 2026-08-14 18:53:23.290412
#> 3 2026-08-14 18:53:34.24768 18:53:47
#> 4 2026-08-14 18:53:47.38769 18:53:51
#> 5 2026-08-14 18:53:51.578525 2026-08-14 18:53:51.579526
#> 6 2026-08-14 18:53:59.349332 2026-08-14 18:54:00.019334
#> command
#> 1 methylumIDAT(barcodes = getBarcodes(path = idatPath), idatPath = idatPath)
#> 2 Subset of 485577 features.
#> 3 methylumi.bgcorr(x = mset450k)
#> 4 normalizeViaControls(x = x)
#> 5 stripOOB(object = normalizeMethyLumiSet(methylumi.bgcorr(mset450k)))
#> 6 asMethod(from = object)
#> lumiVersion
#> 1 <NA>
#> 2 <NA>
#> 3 <NA>
#> 4 <NA>
#> 5 <NA>
#> 6 2.59.3MethyLumiSet objects with OOB matrices can be coerced to RGChannelSet objects for further processing using functions found in the minfi or ChAMP packages.
suppressPackageStartupMessages(library(FDb.InfiniumMethylation.hg19))
rgSet450k <- as(mset450k, "RGChannelSet")
#> Fetching coordinates for hg19...
show(rgSet450k)
#> class: RGChannelSet
#> dim: 621928 10
#> metadata(0):
#> assays(2): Green Red
#> rownames(621928): 14782418 12709357 ... 28673402 13742412
#> rowData names(0):
#> colnames(10): TCGA1 TCGA2 ... TCGA9 TCGA10
#> colData names(1): barcode
#> Annotation
#> array: IlluminaHumanMethylation450kThe above will not work for the processed data, but only because we called stripOOB() on the resulting object to reduce its size. If you plan on using a preprocessed MethyLumiSet in minfi for further processing, don’t strip it.
The GenomicMethylSet and GenomicRatioSet classes in minfi inherit from the RangedSummarizedExperiment class, which has some particularly useful features:
suppressPackageStartupMessages({
library(minfi)
library(IlluminaHumanMethylation450kanno.ilmn12.hg19)
})
grSet450k <- mapToGenome(mset450k.andBack)
sexChroms <- GRanges(seqnames = c("chrX", "chrY"),
IRanges(start = c(1, 1),
end = c(155270560, 59373566)),
strand = c("*", "*"))
summary(subsetByOverlaps(grSet450k, sexChroms))
#> [1] "GenomicMethylSet object of length 11648 with 0 metadata columns"
dim(subsetByOverlaps(grSet450k, sexChroms))
#> [1] 11648 10These SummarizedExperiment-derived objects can be subsetted by nearly anything that has an interval-based representation. Here we extract some promoters, but one could just as easily use AnnotationHub resources to find CTCF peaks or, say, H3K4me1 peaks in ChIP-seq data (often associated with transcriptional enhancers; the presence or absence of DNA methylation may help determine their activity).
## perhaps more topical:
suppressPackageStartupMessages({
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(Homo.sapiens)
})
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
KDM6AEntrezID <- org.Hs.egSYMBOL2EG[["KDM6A"]]
txs.KDM6A <- transcriptsBy(txdb, "gene")[[KDM6AEntrezID]]
tss.KDM6A <- unique(resize(txs.KDM6A, 1, fix = "start")) ## two start sites
promoters.KDM6A <- flank(tss.KDM6A, 100) ## an arbitrary distance upstream
show(subsetByOverlaps(grSet450k, promoters.KDM6A)) ## probes in this window
#> class: GenomicMethylSet
#> dim: 6 10
#> metadata(0):
#> assays(2): Meth Unmeth
#> rownames(6): cg14384228 cg07167981 ... cg17824914 cg06877198
#> rowData names(0):
#> colnames(10): TCGA1 TCGA2 ... TCGA9 TCGA10
#> colData names(9): barcode mu.Cy3 ... alpha.Cy5 offset.Cy5
#> Annotation
#> array: IlluminaHumanMethylation450k
#> annotation: ilmn12.hg19
#> Preprocessing
#> Method: methylumi, background corrected, dye bias equalized
#> minfi version: 1.59.1
#> Manifest version: 0.4Consult the AnnotationHub package vignette for some other possibilities. If you are unfamiliar with the powerful GenomicRanges and GenomicFeatures packages, you may want to familiarize yourself with them as well.
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] parallel stats4 stats graphics grDevices utils datasets
#> [8] methods base
#>
#> other attached packages:
#> [1] Homo.sapiens_1.3.1
#> [2] GO.db_3.23.1
#> [3] OrganismDbi_1.55.1
#> [4] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1
#> [5] lumi_2.65.0
#> [6] TCGAMethylation450k_1.49.0
#> [7] methylumi_2.59.3
#> [8] minfi_1.59.1
#> [9] bumphunter_1.55.1
#> [10] locfit_1.5-9.12
#> [11] iterators_1.0.14
#> [12] foreach_1.5.2
#> [13] Biostrings_2.81.6
#> [14] XVector_0.53.0
#> [15] SummarizedExperiment_1.43.0
#> [16] MatrixGenerics_1.25.0
#> [17] FDb.InfiniumMethylation.hg19_2.2.0
#> [18] org.Hs.eg.db_3.23.1
#> [19] TxDb.Hsapiens.UCSC.hg19.knownGene_3.22.1
#> [20] GenomicFeatures_1.65.0
#> [21] AnnotationDbi_1.75.2
#> [22] GenomicRanges_1.65.1
#> [23] Seqinfo_1.3.0
#> [24] IRanges_2.47.2
#> [25] S4Vectors_0.51.6
#> [26] matrixStats_1.5.0
#> [27] ggplot2_4.0.3
#> [28] reshape2_1.4.5
#> [29] scales_1.4.0
#> [30] Biobase_2.73.2
#> [31] BiocGenerics_0.59.12
#> [32] generics_0.1.4
#>
#> loaded via a namespace (and not attached):
#> [1] RColorBrewer_1.1-3 jsonlite_2.0.0
#> [3] magrittr_2.0.5 farver_2.1.2
#> [5] rmarkdown_2.31 BiocIO_1.23.3
#> [7] vctrs_0.7.3 multtest_2.69.0
#> [9] memoise_2.0.1 Rsamtools_2.29.0
#> [11] DelayedMatrixStats_1.35.0 RCurl_1.98-1.19
#> [13] askpass_1.2.1 htmltools_0.5.9
#> [15] S4Arrays_1.13.0 BiocBaseUtils_1.15.1
#> [17] curl_7.1.0 Rhdf5lib_2.1.0
#> [19] SparseArray_1.13.2 rhdf5_2.57.10
#> [21] KernSmooth_2.23-27 nor1mix_1.3-3
#> [23] plyr_1.8.9 cachem_1.1.0
#> [25] GenomicAlignments_1.49.1 lifecycle_1.0.5
#> [27] pkgconfig_2.0.3 Matrix_1.7-6
#> [29] R6_2.6.1 fastmap_1.2.0
#> [31] digest_0.6.39 siggenes_1.87.0
#> [33] reshape_0.8.10 RSQLite_3.53.3
#> [35] base64_2.0.2 mgcv_1.9-4
#> [37] httr_1.4.8 abind_1.4-8
#> [39] compiler_4.6.1 beanplot_1.3.1
#> [41] rngtools_1.5.2 bit64_4.8.2
#> [43] withr_3.0.3 S7_0.2.2
#> [45] BiocParallel_1.47.0 DBI_1.3.0
#> [47] HDF5Array_1.41.2 MASS_7.3-66
#> [49] openssl_2.4.2 DelayedArray_0.39.5
#> [51] rjson_0.2.23 tools_4.6.1
#> [53] otel_0.2.0 rentrez_1.2.4
#> [55] glue_1.8.1 quadprog_1.5-8
#> [57] h5mread_1.5.0 restfulr_0.0.17
#> [59] nlme_3.1-170 rhdf5filters_1.25.4
#> [61] grid_4.6.1 gtable_0.3.6
#> [63] tzdb_0.5.0 preprocessCore_1.75.0
#> [65] tidyr_1.3.2 data.table_1.18.4
#> [67] hms_1.1.4 xml2_1.6.0
#> [69] pillar_1.11.1 stringr_1.6.0
#> [71] limma_3.69.4 genefilter_1.95.0
#> [73] splines_4.6.1 dplyr_1.2.1
#> [75] lattice_0.23-1 survival_3.8-9
#> [77] rtracklayer_1.73.0 bit_4.6.0
#> [79] GEOquery_2.81.28 annotate_1.91.0
#> [81] RBGL_1.89.0 tidyselect_1.2.1
#> [83] knitr_1.51 xfun_0.60
#> [85] scrime_1.3.7 statmod_1.5.2
#> [87] UCSC.utils_1.9.0 stringi_1.8.9
#> [89] yaml_2.3.12 evaluate_1.0.5
#> [91] codetools_0.2-20 cigarillo_1.3.1
#> [93] tibble_3.3.1 graph_1.91.0
#> [95] BiocManager_1.30.27 affyio_1.83.0
#> [97] cli_3.6.6 xtable_1.8-8
#> [99] dichromat_2.0-1 Rcpp_1.1.2
#> [101] GenomeInfoDb_1.49.1 png_0.1-9
#> [103] XML_3.99-0.23 readr_2.2.0
#> [105] blob_1.3.0 mclust_6.1.3
#> [107] doRNG_1.8.6.3 sparseMatrixStats_1.25.0
#> [109] bitops_1.1-0 affy_1.91.0
#> [111] illuminaio_0.55.0 nleqslv_3.3.7
#> [113] purrr_1.2.2 crayon_1.5.3
#> [115] rlang_1.3.0 KEGGREST_1.53.6