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

###################################################
### code chunk number 1: bioconductor (eval = FALSE)
###################################################
## source("http://bioconductor.org/biocLite.R")
## biocLite("DMRcate")


###################################################
### code chunk number 2: libr
###################################################
library(DMRcate)


###################################################
### code chunk number 3: loaddata
###################################################
data(dmrcatedata)
myMs <- logit2(myBetas)


###################################################
### code chunk number 4: filter
###################################################
nrow(illuminaSNPs)
nrow(myMs)
myMs.noSNPs <- rmSNPandCH(myMs, dist=2, mafcut=0.05)
nrow(myMs.noSNPs)


###################################################
### code chunk number 5: annotate
###################################################
patient <- factor(sub("-.*", "", colnames(myMs)))
type <- factor(sub(".*-", "", colnames(myMs)))
design <- model.matrix(~patient + type) 
myannotation <- cpg.annotate("array", myMs.noSNPs, analysis.type="differential",
    design=design, coef=39)


###################################################
### code chunk number 6: dmrcate
###################################################
dmrcoutput <- dmrcate(myannotation, lambda=1000, C=2)


###################################################
### code chunk number 7: ranges
###################################################
results.ranges <- extractRanges(dmrcoutput, genome = "hg19")
results.ranges


###################################################
### code chunk number 8: plotting
###################################################
groups <- c(Tumour="magenta", Normal="forestgreen")
cols <- groups[as.character(type)]
samps <- c(1:6, 38+(1:6))
DMR.plot(ranges=results.ranges, dmr=1, CpGs=myBetas, phen.col=cols, genome="hg19", 
         samps=samps)


###################################################
### code chunk number 9: wgbssuitecpgs
###################################################
CpGs


###################################################
### code chunk number 10: prepareDSS
###################################################
meth <- as.data.frame(CpGs)[,c(1:2, grep(".C$", colnames(as.data.frame(CpGs))))]
coverage <- as.data.frame(CpGs)[,c(1:2, grep(".cov$", colnames(as.data.frame(CpGs))))]

treat1 <- data.frame(chr=coverage$seqnames, pos=coverage$start, 
                     N=coverage$Treatment1.cov, X=meth$Treatment1.C)

treat2 <- data.frame(chr=coverage$seqnames, pos=coverage$start, 
                     N=coverage$Treatment2.cov, X=meth$Treatment2.C)

treat3 <- data.frame(chr=coverage$seqnames, pos=coverage$start, 
                     N=coverage$Treatment3.cov, X=meth$Treatment3.C)

ctrl1 <- data.frame(chr=coverage$seqnames, pos=coverage$start, 
                     N=coverage$Control1.cov, X=meth$Control1.C)

ctrl2 <- data.frame(chr=coverage$seqnames, pos=coverage$start, 
                     N=coverage$Control2.cov, X=meth$Control2.C)

ctrl3 <- data.frame(chr=coverage$seqnames, pos=coverage$start, 
                     N=coverage$Control3.cov, X=meth$Control3.C)

samples <- list(treat1, treat2, treat3, ctrl1, ctrl2, ctrl3)
sampnames <- sub("\\..*", "", colnames(meth))[-c(1:2)]

obj_bsseq <- makeBSseqData(samples, sampnames)
DSSres <- DMLtest(obj_bsseq, group1=sampnames[1:3], group2=sampnames[4:6], smoothing=FALSE) 



###################################################
### code chunk number 11: wgbsDMRcate
###################################################
wgbsannot <- cpg.annotate("sequencing", DSSres)
wgbs.DMRs <- dmrcate(wgbsannot, lambda = 1000, C = 50, pcutoff = 0.05, mc.cores = 1)
wgbs.ranges <- extractRanges(wgbs.DMRs, genome = "hg19")
groups <- c(Treatment="darkorange", Control="blue")
cols <- groups[sub("[0-9]", "", sampnames)]
DMR.plot(ranges=wgbs.ranges, dmr=1, CpGs=CpGs, phen.col=cols, genome="hg19")


###################################################
### code chunk number 12: sessionInfo
###################################################
sessionInfo()