## ----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_KEGG <- goregion(results.ranges[results.ranges$meandiff < 0][1:100], all.cpg = rownames(tcell), collection = "KEGG", array.type = "EPIC") enrichment_KEGG <- enrichment_KEGG[order(enrichment_KEGG$P.DE),] head(as.matrix(enrichment_KEGG), 10) ## ----loadeh, message=FALSE---------------------------------------------------- bis_1072 <- eh[["EH1072"]] bis_1072 colnames(bis_1072) ## ----bisphen------------------------------------------------------------------ bsseq::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()