DMRcate.R
In DMRcate: Methylation array and sequencing spatial analysis methods

## ----bioconductor, message=FALSE, warning=FALSE, eval=FALSE-------------------
#  if (!require("BiocManager"))
#  	install.packages("BiocManager")
#  BiocManager::install("DMRcate")

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

## ----tcells, message=FALSE----------------------------------------------------
library(ExperimentHub)
eh <- ExperimentHub()
FlowSorted.Blood.EPIC <- eh[["EH1136"]]  
tcell <- FlowSorted.Blood.EPIC[,colData(FlowSorted.Blood.EPIC)$CD4T==100 | 
                                 colData(FlowSorted.Blood.EPIC)$CD8T==100]

## ----detpnorm-----------------------------------------------------------------
detP <- detectionP(tcell)
remove <- apply(detP, 1, function (x) any(x > 0.01))
tcell <- preprocessFunnorm(tcell)
tcell <- tcell[!rownames(tcell) %in% names(which(remove)),]
tcellms <- getM(tcell)

## ----filter, message=FALSE----------------------------------------------------
nrow(tcellms)
tcellms.noSNPs <- rmSNPandCH(tcellms, dist=2, mafcut=0.05)
nrow(tcellms.noSNPs)

## ----avearrays----------------------------------------------------------------
tcell$Replicate
tcell$Replicate[tcell$Replicate==""] <- tcell$Sample_Name[tcell$Replicate==""]
tcellms.noSNPs <- limma::avearrays(tcellms.noSNPs, tcell$Replicate)
tcell <- tcell[,!duplicated(tcell$Replicate)]
tcell <- tcell[rownames(tcellms.noSNPs),]
colnames(tcellms.noSNPs) <- colnames(tcell)
assays(tcell)[["M"]] <- tcellms.noSNPs
assays(tcell)[["Beta"]] <- ilogit2(tcellms.noSNPs)

## ----annotate, message=FALSE--------------------------------------------------
type <- factor(tcell$CellType)
design <- model.matrix(~type) 
myannotation <- cpg.annotate("array", tcell, arraytype = "EPIC",
                             analysis.type="differential", design=design, coef=2)

## ----showmyannotation---------------------------------------------------------
myannotation

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

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

## ----plotting, message=FALSE--------------------------------------------------
groups <- c(CD8T="magenta", CD4T="forestgreen")
cols <- groups[as.character(type)]
cols

DMR.plot(ranges=results.ranges, dmr=1, CpGs=getBeta(tcell), what="Beta", 
         arraytype = "EPIC", phen.col=cols, genome="hg19")

## ----goregion-----------------------------------------------------------------
library(missMethyl)
enrichment_GO <- goregion(results.ranges[1:100], all.cpg = rownames(tcell), 
                          collection = "GO", array.type = "EPIC")
enrichment_GO <- enrichment_GO[order(enrichment_GO$P.DE),]
head(as.matrix(enrichment_GO), 10)

## ----loadeh, message=FALSE----------------------------------------------------
bis_1072 <- eh[["EH1072"]]
bis_1072
colnames(bis_1072)

## ----bisphen------------------------------------------------------------------
pData(bis_1072) <- data.frame(replicate=gsub(".*-", "", colnames(bis_1072)),
                              tissue=substr(colnames(bis_1072), 1, 
                                            nchar(colnames(bis_1072))-3), 
                              row.names=colnames(bis_1072))
colData(bis_1072)$tissue <- gsub("-", "_", colData(bis_1072)$tissue)
as.data.frame(colData(bis_1072))

## ----changeseqlevs------------------------------------------------------------
bis_1072 <- renameSeqlevels(bis_1072, mapSeqlevels(seqlevels(bis_1072), "UCSC"))

## ----chr2filter---------------------------------------------------------------
bis_1072 <- bis_1072[seqnames(bis_1072)=="chr19",]
bis_1072

## ----bsdesign, message=FALSE--------------------------------------------------
tissue <- factor(pData(bis_1072)$tissue)
tissue <- relevel(tissue, "Liver_Treg")

#Regular matrix design
design <- model.matrix(~tissue)
colnames(design) <- gsub("tissue", "", colnames(design))
colnames(design)[1] <- "Intercept"
rownames(design) <- colnames(bis_1072)
design

#Methylation matrix design
methdesign <- edgeR::modelMatrixMeth(design)
methdesign

## ----fitBSseq-----------------------------------------------------------------
cont.mat <- limma::makeContrasts(treg_vs_tcon=Lymph_N_Treg-Lymph_N_Tcon,
                                 fat_vs_ln=Fat_Treg-Lymph_N_Treg,
                                 skin_vs_ln=Skin_Treg-Lymph_N_Treg,
                                 fat_vs_skin=Fat_Treg-Skin_Treg,
                                 levels=methdesign)
cont.mat

## ----sequencingannotate-------------------------------------------------------
seq_annot <- sequencing.annotate(bis_1072, methdesign, all.cov = TRUE, 
                                   contrasts = TRUE, cont.matrix = cont.mat, 
                                   coef = "treg_vs_tcon", fdr=0.05)
seq_annot

## ----seqdmrcate---------------------------------------------------------------
dmrcate.res <- dmrcate(seq_annot, C=2, min.cpgs = 5)
dmrcate.res
treg_vs_tcon.ranges <- extractRanges(dmrcate.res, genome="mm10")
treg_vs_tcon.ranges

## ----seqDMRplot1, message=FALSE-----------------------------------------------
cols <- as.character(plyr::mapvalues(tissue, unique(tissue), 
                                     c("darkorange", "maroon", "blue", 
                                       "black", "magenta")))
names(cols) <- tissue

DMR.plot(treg_vs_tcon.ranges, dmr = 1, 
         CpGs=bis_1072[,tissue %in% c("Lymph_N_Tcon", "Lymph_N_Treg")], 
         phen.col = cols[tissue %in% c("Lymph_N_Tcon", "Lymph_N_Treg")], 
         genome="mm10")

## ----fatskin------------------------------------------------------------------
seq_annot <- sequencing.annotate(bis_1072, methdesign, all.cov = TRUE, 
                                   contrasts = TRUE, cont.matrix = cont.mat, 
                                   coef = "fat_vs_skin", fdr=0.05)

## ----redefinethresh-----------------------------------------------------------
seq_annot <- changeFDR(seq_annot, 0.25)

## ----dmrsfatskin--------------------------------------------------------------
dmrcate.res <- dmrcate(seq_annot, C=2, min.cpgs = 5)
fat_vs_skin.ranges <- extractRanges(dmrcate.res, genome="mm10")

## ----seqDMRplot2, message=FALSE-----------------------------------------------
cols
DMR.plot(fat_vs_skin.ranges, dmr = 1, CpGs=bis_1072, phen.col = cols, genome="mm10")

## ----DSS, message=FALSE-------------------------------------------------------
library(DSS)
DMLfit <- DMLfit.multiFactor(bis_1072, design=data.frame(tissue=tissue), formula=~tissue)
DSS_treg.vs.tcon <- DMLtest.multiFactor(DMLfit, Contrast=matrix(c(0, 0, -1, 1, 0)))
#Make sure to filter out all sites where the test statistic is NA
DSS_treg.vs.tcon <- DSS_treg.vs.tcon[!is.na(DSS_treg.vs.tcon$stat),]

seq_annot <- sequencing.annotate(obj=DSS_treg.vs.tcon, fdr=0.05)
seq_annot

dmrcate.res <- dmrcate(seq_annot, C=2, min.cpgs = 5)
DSS.treg_vs_tcon.ranges <- extractRanges(dmrcate.res, genome="mm10")

findOverlaps(treg_vs_tcon.ranges, DSS.treg_vs_tcon.ranges)

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