### R code from vignette source 'ENmix.Rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: options ################################################### options(width=70) ################################################### ### code chunk number 2: example (eval = FALSE) ################################################### ## library(ENmix) ## #read in data ## require(minfiData) ## #data pre-processing ## beta=mpreprocess(RGsetEx,nCores=6) ## #or ## #read in IDAT files ## sheet <- read.metharray.sheet(file.path(find.package("minfiData"), ## "extdata"), pattern = "csv$") ## rgSet <- read.metharray.exp(targets = sheet,extended = TRUE) ## #quality control and data pre-processing ## beta=mpreprocess(rgSet,nCores=6,qc=TRUE,foutlier=TRUE, ## rmcr=TRUE,impute=TRUE) ################################################### ### code chunk number 3: example (eval = FALSE) ################################################### ## library(ENmix) ## #read in data ## sheet <- read.metharray.sheet(file.path(find.package("minfiData"), ## "extdata"), pattern = "csv$") ## rgSet <- read.metharray.exp(targets = sheet, extended = TRUE) ## #QC info ## qc<-QCinfo(rgSet) ## #background correction and dye bias correction ## #if provide qc info, the low quality samples and probes ## #will be excluded before background correction ## mdat<-preprocessENmix(rgSet, bgParaEst="oob", dyeCorr="RELIC", ## QCinfo=qc, nCores=6) ## #low quality samples and probes can also be excluded after ## #background correction using QCfilter. ## #mdat <- QCfilter(mdat,qcinfo=qc,outlier=TRUE) ## #inter-array normalization ## mdat<-norm.quantile(mdat, method="quantile1") ## #probe-type bias adjustment ## beta<-rcp(mdat,qcscore=qc) ## beta <- rm.outlier(beta,qcscore=qc,impute=TRUE,rmcr=TRUE) ################################################### ### code chunk number 4: example (eval = FALSE) ################################################### ## library(ENmix) ## #read in data ## sheet <- read.metharray.sheet(file.path(find.package("minfiData"), ## "extdata"), pattern = "csv$") ## rgSet <- read.metharray.exp(targets = sheet, extended = TRUE) ## #control plots ## plotCtrl(rgSet) ## #QC info ## qc<-QCinfo(rgSet) ## mraw <- preprocessRaw(rgSet) ## beta<-getBeta(mraw, "Illumina") ## #distribution plot ## multifreqpoly(beta,main="Methylation Beta value distribution") ## #Search for multimodal CpGs ## #sample size in this example data is too small for this purpose! ## #exclude low quality data first ## bb=beta; bb[qc$detP>0.05 | qc$nbead<3]=NA ## nmode<-nmode.mc(bb, minN = 3, modedist=0.2, nCores = 6) ## outCpG = names(nmode)[nmode>1] ## #background correction and dye bias correction ## mdat<-preprocessENmix(rgSet, bgParaEst="oob", dyeCorr="RELIC", ## QCinfo=qc, exCpG=outCpG, nCores=6) ## #inter-array normalization ## mdat<-norm.quantile(mdat, method="quantile1") ## #probe-type bias adjustment ## beta<-rcp(mdat,qcscore=qc) ## # Principal component regression analysis plot ## cov<-data.frame(group=pData(mdat)$Sample_Group, ## slide=factor(pData(mdat)$Slide)) ## pcrplot(beta, cov, npc=6) ## #filter out low quality and outlier data points for each probe; ## #rows and columns with too many missing value can be removed ## #if specify; Do imputation to fill missing data if specify. ## beta <- rm.outlier(beta,qcscore=qc,rmcr=TRUE,impute=TRUE) ## #Non-negative control surrogate variables ## sva<-ctrlsva(rgSet) ################################################### ### code chunk number 5: UnevaluatedCode (eval = FALSE) ################################################### ## library(ENmix) ## require(minfi) ## #see minfi user's guide for the format of sample_sheet.txt file ## targets <- read.table("./sample_sheet.txt", header=T) ## rgSet <- read.metharray.exp( targets = targets, extended = TRUE) ## # or read in all idat files under a directory ## rgSet <- read.metharray.exp(base = "path_to_directory_idat_files", ## targets = NULL, extended = TRUE, recursive=TRUE) ################################################### ### code chunk number 6: UnevaluatedCode (eval = FALSE) ################################################### ## M<-matrix_for_methylated_intensity ## U<-matrix_for_unmethylated_intensity ## pheno<-as.data.frame(cbind(colnames(M), colnames(M))) ## names(pheno)<-c("Basename","filenames") ## rownames(pheno)<-pheno$Basename ## pheno<-AnnotatedDataFrame(data=pheno) ## anno<-c("IlluminaHumanMethylation450k", "ilmn12.hg19") ## names(anno)<-c("array", "annotation") ## mdat<-MethylSet(Meth = M, Unmeth = U, annotation=anno, ## phenoData=pheno) ################################################### ### code chunk number 7: load (eval = FALSE) ################################################### ## library(ENmix) ## require(minfi) ## require(minfiData) ## sheet <- read.metharray.sheet(file.path(find.package("minfiData"), ## "extdata"), pattern = "csv$") ## rgSet <- read.metharray.exp(targets = sheet, extended = TRUE) ################################################### ### code chunk number 8: ctrlplot (eval = FALSE) ################################################### ## plotCtrl(rgSet) ################################################### ### code chunk number 9: ctrlplot (eval = FALSE) ################################################### ## pinfo=pData(rgSet) ## IDorder=rownames(pinfo)[order(pinfo$Slide,pinfo$Array)] ## plotCtrl(rgSet,IDorder) ################################################### ### code chunk number 10: ctrlplot (eval = FALSE) ################################################### ## mraw <- preprocessRaw(rgSet) ## #total intensity plot is userful for data quality inspection ## #and identification of outlier samples ## multifreqpoly(assayData(mraw)$Meth+assayData(mraw)$Unmeth, ## xlab="Total intensity") ## #Compare frequency polygon plot and density plot ## beta<-getBeta(mraw, "Illumina") ## anno=getAnnotation(rgSet) ## beta1=beta[anno$Type=="I",] ## beta2=beta[anno$Type=="II",] ## library(geneplotter) ## jpeg("dist.jpg",height=900,width=600) ## par(mfrow=c(3,2)) ## multidensity(beta,main="Multidensity") ## multifreqpoly(beta,main="Multifreqpoly",xlab="Beta value") ## multidensity(beta1,main="Multidensity: Infinium I") ## multifreqpoly(beta1,main="Multifreqpoly: Infinium I", ## xlab="Beta value") ## multidensity(beta2,main="Multidensity: Infinium II") ## multifreqpoly(beta2,main="Multifreqpoly: Infinium II", ## xlab="Beta value") ## dev.off() ################################################### ### code chunk number 11: filter (eval = FALSE) ################################################### ## qc<-QCinfo(rgSet) ## #exclude before backgroud correction ## mdat<-preprocessENmix(rgSet, bgParaEst="oob", dyeCorr="RELIC", ## QCinfo=qc, nCores=6) ## #Or exclude after background correction ## mdat <- QCfilter(mdat,qcinfo=qc,outlier=TRUE) ################################################### ### code chunk number 12: preprocessENmix (eval = FALSE) ################################################### ## #filter out outliers ## b1=rm.outlier(beta) ## #filter out low quality and outlier values ## b2=rm.outlier(beta,qcscore=qcscore) ## #filter out low quality and outlier values, remove rows and columns ## # with too many missing values ## b3=rm.outlier(beta,qcscore=qcscore,rmcr=TRUE) ## #filter out low quality and outlier values, remove rows and columns ## # with too many missing values, and then do imputation ## b3=rm.outlier(beta,qcscore=qcscore,rmcr=TRUE,impute=TRUE) ################################################### ### code chunk number 13: preprocessENmix (eval = FALSE) ################################################### ## qc=QCinfo(rgSet) ## mdat<-preprocessENmix(rgSet, bgParaEst="oob", dyeCorr="RELIC", ## QCinfo=qc, exCpG=NULL, nCores=6) ################################################### ### code chunk number 14: normalize.quantile.450k (eval = FALSE) ################################################### ## mdat<-norm.quantile(mdat, method="quantile1") ################################################### ### code chunk number 15: rcp (eval = FALSE) ################################################### ## beta<-rcp(mdat) ################################################### ### code chunk number 16: bmiq.mc (eval = FALSE) ################################################### ## beta<-bmiq.mc(mdat, nCores=6) ################################################### ### code chunk number 17: ctrlsva (eval = FALSE) ################################################### ## require(minfiData) ## sheet <- read.metharray.sheet(file.path(find.package("minfiData"), ## "extdata"), pattern = "csv$") ## rgSet <- read.metharray.exp(targets = sheet,extended = TRUE) ## sva<-ctrlsva(rgSet) ################################################### ### code chunk number 18: combat.mc (eval = FALSE) ################################################### ## batch<-factor(pData(mdat)$Slide) ## betaC<-ComBat.mc(beta, batch, nCores=6, mod=NULL) ################################################### ### code chunk number 19: pcrplot (eval = FALSE) ################################################### ## cov<-data.frame(group=pData(mdat)$Sample_Group, ## slide=factor(pData(mdat)$Slide)) ## pcrplot(beta, cov, npc=6) ################################################### ### code chunk number 20: nmode.mc (eval = FALSE) ################################################### ## nmode<- nmode.mc(beta, minN = 3, modedist=0.2, nCores = 5) ################################################### ### code chunk number 21: ENmixAndminfi (eval = FALSE) ################################################### ## library(ENmix) ## #minfi functions to read in data ## sheet <- read.metharray.sheet(file.path(find.package("minfiData"), ## "extdata"), pattern = "csv$") ## rgSet <- read.metharray.exp(targets = sheet, extended = TRUE) ## #ENmix function for control plot ## plotCtrl(rgSet) ## #minfi functions to extract methylation and annotation data ## mraw <- preprocessRaw(rgSet) ## beta<-getBeta(mraw, "Illumina") ## anno=getAnnotation(rgSet) ## #ENmix function for fast and accurate distribution plot ## multifreqpoly(beta,main="Data distribution") ## multifreqpoly(beta[anno$Type=="I",],main="Data distribution, type I") ## multifreqpoly(beta[anno$Type=="II",],main="Data distribution, type II") ## #ENmix background correction ## mset<-preprocessENmix(rgSet, bgParaEst="oob", dyeCorr="RELIC", nCores=6) ## #minfi functions for further preprocessing and analysis ## gmSet <- preprocessQuantile(mset) ## bumps <- bumphunter(gmSet, design = model.matrix(~ gmSet$status), B = 0, ## type = "Beta", cutoff = 0.25) ################################################### ### code chunk number 22: ENmixAndChAMP (eval = FALSE) ################################################### ## library(ENmix) ## library(ChAMP) ## testDir=system.file("extdata",package="ChAMPdata") ## myLoad=champ.load(directory=testDir) ## #ENmix background correction ## mset<-preprocessENmix(myLoad$rgSet,bgParaEst="oob", nCores=6) ## #remove probes filtered by champ.load() ## mset=mset[rownames(myLoad$beta),] ## #update myLoad object with background corrected intensity data ## myLoad$mset=mset ## myLoad$beta=getBeta(mset) ## myLoad$intensity=getMeth(mset)+getUnmeth(mset) ## #continue ChAMP pipeline ## myNorm=champ.norm() ################################################### ### code chunk number 23: sessionInfo ################################################### toLatex(sessionInfo())