wateRmelon.R
In wateRmelon: Illumina 450 methylation array normalization and metrics

### R code from vignette source 'wateRmelon.Rnw'

###################################################
### code chunk number 1: UnevaluatedCode (eval = FALSE)
###################################################
## install.packages('ROCR', 'matrixStats')
## if (!requireNamespace("BiocManager", quietly=TRUE))
##     install.packages("BiocManager")
## BiocManager::install( 'limma', 'minfi',
##    'IlluminaHumanMethylation450kmanifest',
##    'methylumi', 'lumi')


###################################################
### code chunk number 2: UnevaluatedCode (eval = FALSE)
###################################################
## install.packages('wateRmelon_0.9.9.tar.gz', repos=NULL, type='source')


###################################################
### code chunk number 3: code-block
###################################################
library('wateRmelon')
# load in melon dataset
data (melon)
# display dimensions of data matrix
dim(melon)
# quality filter using default thresholds
melon.pf<-pfilter(melon)
# preprocess using our best method
melon.dasen.pf <- dasen(melon.pf)


###################################################
### code chunk number 4: dmrse
###################################################
# calculate iDMR metrics on QC'd betas
dmrse_row(melon.pf)
# calculate iDMR metrics on QC'd and preprocessed betas
dmrse_row(melon.dasen.pf)
# slightly lower (better) standard errors


###################################################
### code chunk number 5: genki
###################################################
# calculate SNP metrics on QC'd betas
genki(melon.pf)
# calculate SNP metrics on QC'd and preprocessed betas
genki(melon.dasen.pf)
# slightly lower (better) standard errors


###################################################
### code chunk number 6: seabi
###################################################
# calculate X-chromosome metrics on QC'd betas
seabi(melon.pf, sex=pData(melon.pf)$sex, X=fData(melon.pf)$CHR=='X')
# calculate X-chromosome metrics on QC'd and preprocessed betas
seabi(melon.dasen.pf,
   sex=pData(melon.dasen.pf)$sex,
   X=fData(melon.dasen.pf)$CHR=='X'
)


###################################################
### code chunk number 7: UnevaluatedCode (eval = FALSE)
###################################################
## library(methylumi)
## melon <- methyLumiR('finalreport.txt')


###################################################
### code chunk number 8: UnevaluatedCodeecc (eval = FALSE)
###################################################
## mlumi <- readEPIC('path/to/directory')


###################################################
### code chunk number 9: IncludeGraphic
###################################################
boxplot(log(methylated(melon)), las=2, cex.axis=0.8 )


###################################################
### code chunk number 10: IncludeGraphic
###################################################
boxplot(log(unmethylated(melon)), las=2, cex.axis=0.8 )


###################################################
### code chunk number 11: outlyx
###################################################
outlyx(melon) # Can take some time on large data-sets


###################################################
### code chunk number 12: bscon
###################################################
bsc <- bscon(melon)
hist(bsc)


###################################################
### code chunk number 13: pfilter
###################################################
melon.pf <- pfilter(melon)


###################################################
### code chunk number 14: dasen
###################################################
melon.dasen.pf <- dasen(melon.pf)


###################################################
### code chunk number 15: agep
###################################################
data(coef)
agep(melon.dasen.pf, coeff= coef, method='horvath')


###################################################
### code chunk number 16: UnevaluatedCodeecc (eval = FALSE)
###################################################
## # Code will not work with melon as melon only has a subset of probes
## estimateCellCounts.wmln(melon.dasen.pf)


###################################################
### code chunk number 17: qual
###################################################
normv <- qual(betas(melon.dasen.pf), betas(melon.pf))
plot(normv[,1:2])


###################################################
### code chunk number 18: pwod
###################################################
melon.pwod.dasen.pf <- pwod(betas(melon.dasen.pf))


###################################################
### code chunk number 19: workflow
###################################################
data(melon)
# load in melon dataset
# Optional read-in data using:
# melon <- readEPIC('path/to/idats')

outliers <- outlyx(melon, plot=F)
sum(outliers$out)
sum(bscon(melon)<85)
# Check for outliers

melon.pf<-pfilter(melon)
# perform QC on raw data matrix using default thresholds
melon.dasen.pf<-dasen(melon.pf)
# preprocess using our best method

qual <- qual(betas(melon.dasen.pf), betas(melon.pf))
# Check for bad samples (again)

sex  <- pData(melon.dasen.pf)$sex
# extract phenotypic information for test
bet<-betas(melon.dasen.pf)
# extract processed beta values
melon.sextest<-sextest(bet,sex)
# run t-test to idenitify sex difference

agep(melon.dasen.pf)
# Check ages (see if they match up e.t.c)

melon.pwod <- pwod(bet)
# Clean up data for statistical testing