inst/doc/ENmix.R

### 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())

Try the ENmix package in your browser

Any scripts or data that you put into this service are public.

ENmix documentation built on May 2, 2018, 4:24 a.m.