## ----bioconductor, message=FALSE, warning=FALSE, eval=FALSE--------------
#  source("http://bioconductor.org/biocLite.R")
#  biocLite("DMRcate")

## ----libr, message=FALSE, warning=FALSE----------------------------------
library(DMRcate)

## ----loaddata------------------------------------------------------------
data(dmrcatedata)
myMs <- logit2(myBetas)

## ----filter--------------------------------------------------------------
nrow(snpsall)
nrow(myMs)
myMs.noSNPs <- rmSNPandCH(myMs, dist=2, mafcut=0.05)
nrow(myMs.noSNPs)

## ----annotate------------------------------------------------------------
patient <- factor(sub("-.*", "", colnames(myMs)))
type <- factor(sub(".*-", "", colnames(myMs)))
design <- model.matrix(~patient + type) 
myannotation <- cpg.annotate("array", myMs.noSNPs, what="M", arraytype = "450K",
                             analysis.type="differential", design=design, coef=39)
#Or, alternatively
grset <- makeGenomicRatioSetFromMatrix(myMs.noSNPs, array = "IlluminaHumanMethylation450k", 
                                       annotation = "ilmn12.hg19", mergeManifest = TRUE, 
                                       what = "M")
myannotation <- cpg.annotate("array", grset, analysis.type="differential", design=design, 
                             coef=39)


## ----dmrcate, warning=FALSE----------------------------------------------
dmrcoutput <- dmrcate(myannotation, lambda=1000, C=2)

## ----ranges--------------------------------------------------------------
results.ranges <- extractRanges(dmrcoutput, genome = "hg19")
results.ranges

## ----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, what="Beta", arraytype = "450K", 
         phen.col=cols, genome="hg19", samps=samps)

## ----wgbssuitecpgs-------------------------------------------------------
CpGs

## ----prepareDSS, warning=FALSE-------------------------------------------
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) 


## ----wgbsDMRcate, warning=FALSE------------------------------------------
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")

## ----sessionInfo---------------------------------------------------------
sessionInfo()