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