In perishky/meffil: Efficient algorithms for DNA methylation

library(knitr)
library(Cairo)
opts_chunk$set(warning=FALSE, fig.width=6, fig.height=6, dev="CairoPNG", stop=TRUE)

library(GEOquery)

Normalize a 450k/EPIC dataset

Download example data sets

path.450k <- download.450k.demo.dataset()

path.epic <- download.epic.demo.dataset()

Normalize dataset

Create samplesheet

library(meffil)
options(mc.cores=10)

samplesheet.epic <- meffil.read.samplesheet(path.epic, pattern="Demo_SampleSheet.csv")
samplesheet.450k <-  meffil.create.samplesheet(path.450k)[1:10,]

common.columns <- intersect(colnames(samplesheet.epic), colnames(samplesheet.450k))
samplesheet <- rbind(samplesheet.epic[,common.columns],
                     samplesheet.450k[,common.columns])

Parameters.

qc.file <- "450k-and-epic/qc-report.html"
author <- "Illumina, et al. and Prickett, et al."
study <- "Normalizing EPIC and 450K microarrays together"
number.pcs <- 5
norm.file <- "450k-and-epic/normalization-report.html"
cell.type.reference <- "blood gse35069"
featureset <- "common" ## select a featureset compatible with both 450k and epic

Generate quality control objects.

qc.objects <- meffil.qc(samplesheet, cell.type.reference=cell.type.reference,
                        featureset=featureset, verbose=T)

QC report.

qc.summary <- meffil.qc.summary(qc.objects, verbose=T)
meffil.qc.report(qc.summary,
                 output.file=qc.file,
                 author=author,
                 study=study)

Normalize dataset.

norm.objects <- meffil.normalize.quantiles(qc.objects, number.pcs=number.pcs, verbose=T)

norm.dataset <- meffil.normalize.samples(norm.objects,
                                        just.beta=F, 
                                        cpglist.remove=qc.summary$bad.cpgs$name,
                                        verbose=T)

Normalization report.

beta <- meffil.get.beta(norm.dataset$M, norm.dataset$U)
pcs <- meffil.methylation.pcs(beta, sites=meffil.get.autosomal.sites("common"), verbose=T)

parameters <- meffil.normalization.parameters(norm.objects)
parameters$batch.threshold <- 0.01
norm.summary <- meffil.normalization.summary(norm.objects=norm.objects,
                                             pcs=pcs,
                                             parameters=parameters, verbose=T)

meffil.normalization.report(norm.summary,
                            output.file=norm.file,
                            author=author,
                            study=study)