### R code from vignette source 'vignettes/missMethyl/inst/doc/missMethyl.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: style-Sweave
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: reading_data
###################################################
library(missMethyl)
library(limma)
library(minfi)
library(minfiData)
baseDir <- system.file("extdata", package = "minfiData")
targets <- read.450k.sheet(baseDir)
targets[,1:9]
targets[,10:12]
rgSet <- read.450k.exp(base = baseDir,targets = targets)


###################################################
### code chunk number 3: swan
###################################################
mSet <- preprocessRaw(rgSet)
mSetSw <- SWAN(mSet,verbose=TRUE)


###################################################
### code chunk number 4: betasByType
###################################################
par(mfrow=c(1,2), cex=1.25)
densityByProbeType(mSet[,1], main = "Raw")
densityByProbeType(mSetSw[,1], main = "SWAN")


###################################################
### code chunk number 5: filtering
###################################################
detP <- detectionP(rgSet)
keep <- rowSums(detP < 0.01) == ncol(rgSet)
mSet <- mSet[keep,]


###################################################
### code chunk number 6: extraction
###################################################
mset_reduced <- mSet[sample(1:nrow(mSet), 20000),]
meth <- getMeth(mset_reduced)
unmeth <- getUnmeth(mset_reduced)
Mval <- log2((meth + 100)/(unmeth + 100))
beta <- getBeta(mset_reduced)
dim(Mval)


###################################################
### code chunk number 7: mdsplot
###################################################
par(mfrow=c(1,1))
plotMDS(Mval, labels=targets$Sample_Name, col=as.integer(factor(targets$status)))
legend("topleft",legend=c("Cancer","Normal"),pch=16,cex=1.2,col=1:2)


###################################################
### code chunk number 8: design
###################################################
library(limma)
group <- factor(targets$status,levels=c("normal","cancer"))
id <- factor(targets$person)
design <- model.matrix(~id + group)
design


###################################################
### code chunk number 9: diffmeth
###################################################
fit <- lmFit(Mval,design)
fit <- eBayes(fit)


###################################################
### code chunk number 10: results
###################################################
summary(decideTests(fit))
top<-topTable(fit,coef=4)
top


###################################################
### code chunk number 11: top4
###################################################
cpgs <- rownames(top)
par(mfrow=c(2,2))
for(i in 1:4){
stripchart(beta[rownames(beta)==cpgs[i],]~design[,4],method="jitter",
group.names=c("Normal","Cancer"),pch=16,cex=1.5,col=c(4,2),ylab="Beta values",
vertical=TRUE,cex.axis=1.5,cex.lab=1.5)
title(cpgs[i],cex.main=1.5)
}


###################################################
### code chunk number 12: diffvar
###################################################
fitvar <- varFit(Mval, design = design, coef = c(1,4))


###################################################
### code chunk number 13: results
###################################################
summary(decideTests(fitvar))
topDV <- topVar(fitvar, coef=4)
topDV


###################################################
### code chunk number 14: alternative
###################################################
design2 <- model.matrix(~0+group+id)
fitvar.contr <- varFit(Mval, design=design2, coef=c(1,2))
contr <- makeContrasts(groupcancer-groupnormal,levels=colnames(design2))
fitvar.contr <- contrasts.varFit(fitvar.contr,contrasts=contr)


###################################################
### code chunk number 15: altresults
###################################################
summary(decideTests(fitvar.contr))
topVar(fitvar.contr,coef=1)


###################################################
### code chunk number 16: top4DV
###################################################
cpgsDV <- rownames(topDV)
par(mfrow=c(2,2))
for(i in 1:4){
stripchart(beta[rownames(beta)==cpgsDV[i],]~design[,4],method="jitter",
group.names=c("Normal","Cancer"),pch=16,cex=1.5,col=c(4,2),ylab="Beta values",
vertical=TRUE,cex.axis=1.5,cex.lab=1.5)
title(cpgsDV[i],cex.main=1.5)
}


###################################################
### code chunk number 17: loadingdata
###################################################
library(tweeDEseqCountData)
data(pickrell1)
counts<-exprs(pickrell1.eset)
dim(counts)
gender <- pickrell1.eset$gender
table(gender)
rm(pickrell1.eset)
data(genderGenes)
data(annotEnsembl63)
annot <- annotEnsembl63[,c("Symbol","Chr")]
rm(annotEnsembl63)


###################################################
### code chunk number 18: dgelist
###################################################
library(edgeR)
y <- DGEList(counts=counts, genes=annot[rownames(counts),])


###################################################
### code chunk number 19: filtering
###################################################
isexpr <- rowSums(cpm(y)>1) >= 20
hasannot <- rowSums(is.na(y$genes))==0
y <- y[isexpr & hasannot,,keep.lib.sizes=FALSE]
dim(y)
y <- calcNormFactors(y)


###################################################
### code chunk number 20: testhapmap
###################################################
design.hapmap <- model.matrix(~gender)
fitvar.hapmap <- varFit(y, design = design.hapmap)
fitvar.hapmap$genes <- y$genes


###################################################
### code chunk number 21: resultshapmap
###################################################
summary(decideTests(fitvar.hapmap))
topDV.hapmap <- topVar(fitvar.hapmap,coef=ncol(design.hapmap))
topDV.hapmap


###################################################
### code chunk number 22: top4DVhapmap
###################################################
genesDV <- rownames(topDV.hapmap)
par(mfrow=c(2,2))
for(i in 1:4){
stripchart(cpm(y,log=TRUE)[rownames(y)==genesDV[i],]~design.hapmap[,ncol(design.hapmap)],method="jitter",
group.names=c("Female","Male"),pch=16,cex=1.5,col=c(4,2),ylab="Log counts per million",
vertical=TRUE,cex.axis=1.5,cex.lab=1.5)
title(genesDV[i],cex.main=1.5)
}


###################################################
### code chunk number 23: sessionInfo
###################################################
toLatex(sessionInfo())