### R code from vignette source 'ENmix.Rnw' ################################################### ### 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 ## path <- file.path(find.package("minfiData"),"extdata") ## rgSet <- readidat(path = path,recursive = 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 ## path <- file.path(find.package("minfiData"),"extdata") ## rgSet <- readidat(path = path,recursive = 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 ## path <- file.path(find.package("minfiData"),"extdata") ## rgSet <- readidat(path = path,recursive = TRUE) ## sheet <- read.metharray.sheet(file.path(find.package("minfiData"),"extdata"), ## pattern = "csv$") ## rownames(sheet)=paste(sheet$Slide,sheet$Array,sep="_") ## colData(rgSet)=as(sheet[colnames(rgSet),],"DataFrame") ## #control plots ## plotCtrl(rgSet) ## #QC info ## qc<-QCinfo(rgSet) ## #calculate detection P values ## detp=calc_detP(rgSet,detPtype = "negative") ## detp2=calc_detP(rgSet,detPtype = "oob") ## #get methDataSet ## mraw <- getmeth(rgSet) ## #get raw methyaltion beta values ## beta<-getB(mraw) ## #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") ## colData(mdat)=as(sheet[colnames(mdat),],"DataFrame") ## #probe-type bias adjustment ## beta<-rcp(mdat,qcscore=qc) ## # Principal component regression analysis plot ## cov<-data.frame(group=colData(mdat)$Sample_Group, ## slide=factor(colData(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 ## csva<-ctrlsva(rgSet) ## #estimate cell type proportions ## #based on rgDataset ## celltype=estimateCellProp(userdata=rgSet,refdata="FlowSorted.Blood.450k", ## nonnegative = TRUE,normalize=TRUE) ## #using methDataSet ## celltype=estimateCellProp(userdata=mdat,refdata="FlowSorted.Blood.450k", ## nonnegative = TRUE,normalize=TRUE) ## #using beta value ## celltype=estimateCellProp(userdata=beta,refdata="FlowSorted.Blood.450k", ## nonnegative = TRUE,normalize=TRUE) ## #predic sex based on rgDataSet or methDataset ## sex=predSex(rgSet) ## sex=predSex(mdat) ## #Methylation age calculation ## mage=methyAge(beta) ## #ICC analysis, see manual ## #DMR analysis, see manual ################################################### ### code chunk number 5: UnevaluatedCode (eval = FALSE) ################################################### ## library(ENmix) ## rgSet <- readidat(path = "directory containing idat files",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: ctrlplot (eval = FALSE) ################################################### ## plotCtrl(rgSet) ################################################### ### code chunk number 8: ctrlplot (eval = FALSE) ################################################### ## ## path <- file.path(find.package("minfiData"),"extdata") ## rgSet <- readidat(path = path,recursive = TRUE) ## sheet <- read.metharray.sheet(file.path(find.package("minfiData"),"extdata"), ## pattern = "csv$") ## rownames(sheet)=paste(sheet$Slide,sheet$Array,sep="_") ## colData(rgSet)=as(sheet[colnames(rgSet),],"DataFrame") ## #control plots ## IDorder=rownames(colData(rgSet))[order(colData(rgSet)$Slide,colData(rgSet)$Array)] ## plotCtrl(rgSet,IDorder) ################################################### ### code chunk number 9: ctrlplot (eval = FALSE) ################################################### ## mraw <- getmeth(rgSet) ## #total intensity plot is userful for data quality inspection ## #and identification of outlier samples ## multifreqpoly(assays(mraw)$Meth+assays(mraw)$Unmeth, ## xlab="Total intensity") ## #Compare frequency polygon plot and density plot ## beta<-getB(mraw) ## anno=rowData(mraw) ## beta1=beta[anno$Infinium_Design_Type=="I",] ## beta2=beta[anno$Infinium_Design_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 10: 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 11: 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 12: preprocessENmix (eval = FALSE) ################################################### ## qc=QCinfo(rgSet) ## mdat<-preprocessENmix(rgSet, bgParaEst="oob", dyeCorr="RELIC", ## QCinfo=qc, exCpG=NULL, nCores=6) ################################################### ### code chunk number 13: normalize.quantile.450k (eval = FALSE) ################################################### ## mdat<-norm.quantile(mdat, method="quantile1") ################################################### ### code chunk number 14: rcp (eval = FALSE) ################################################### ## beta<-rcp(mdat) ################################################### ### code chunk number 15: ctrlsva (eval = FALSE) ################################################### ## sva<-ctrlsva(rgSet) ################################################### ### code chunk number 16: pcrplot (eval = FALSE) ################################################### ## cov<-data.frame(group=colData(mdat)$Sample_Group, ## slide=factor(colData(mdat)$Slide)) ## pcrplot(beta, cov, npc=6) ################################################### ### code chunk number 17: nmode.mc (eval = FALSE) ################################################### ## nmode<- nmode.mc(beta, minN = 3, modedist=0.2, nCores = 5) ################################################### ### code chunk number 18: 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<-getB(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 19: 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=getB(mset) ## myLoad$intensity=getMeth(mset)+getUnmeth(mset) ## #continue ChAMP pipeline ## myNorm=champ.norm() ################################################### ### code chunk number 20: sessionInfo ################################################### toLatex(sessionInfo())