Differential methylation calling

DSS

A very important step of a WGBS data analysis is the differential methylation calling. During this step, we would like to identify changes between groups or conditions first on the site level. Then, we will need to identify larger regions (DMR-s, differentially methylated regions) affected by methylation changes that are more likely to represent functional alterations. Methrix doesn’t have a DMR caller, therefore we will use DSS. At this stage, DSS is not yet able to work directly with methrix, therefore we transform our methrix object to bsseq for the sake of DMR calling.


if(!requireNamespace("DSS")) {
BiocManager::install("DSS")
}
library(DSS)

For more detailed description, options and recommendations on smoothing, please refer to the ´DSS´ vignette

bs_obj <- methrix::methrix2bsseq(meth)


dmlTest = DSS::DMLtest(bs_obj, group1=rownames(sample_anno)[sample_anno$Day=="D6"], group2=rownames(sample_anno)[sample_anno$Day=="D14"], 
                       smoothing = TRUE)

Call DMRs and DMLs


dmls = DSS::callDML(dmlTest, p.threshold = 0.05)
dmrs = DSS::callDMR(dmlTest, p.threshold = 0.05)

DMRs:


kable(dmrs) %>%
  kable_styling(bootstrap_options = "striped", full_width = F) %>%
  scroll_box(width = "100%", height = "200px")

DMLs:


kable(dmls) %>%
  kable_styling(bootstrap_options = "striped", full_width = F) %>%
  scroll_box(width = "100%", height = "200px")

DMR calling with `dmrseq`

A popular R package for DMR calling is dmrseq. For more detailed description of the possibilites, see the vignette.



if(!requireNamespace("dmrseq")) {
BiocManager::install("dmrseq")
}
library(dmrseq)

register(MulticoreParam(2))

testCovariate <- "Day"
regions <- dmrseq(bs=bs_obj[240001:260000,],
                  cutoff = 0.05,
                  testCovariate=testCovariate, verbose = T)
saveRDS(regions, "dmrseq_regions.RDS")

Visualizing DMRs and candidate regions

Visualizing a region

We can visualize certain regions and/or DMRs using the Gviz package. The plotting function is included in the best_practices_WGBS.Rmd file.

dmrs <- makeGRangesFromDataFrame(dmrs, keep.extra.columns = T)

genome <- "hg19"
mart <- useMart(biomart="ENSEMBL_MART_ENSEMBL", host="feb2014.archive.ensembl.org", dataset="hsapiens_gene_ensembl")



i <- 5

#candidate_region <- GRanges(paste0(seqlevels(dmrs)[i], ":", start(dmrs)[i], "-", end(dmrs)[i]))
candidate_region <- dmrs[i]
region_plot(m = meth, region=candidate_region, mart=mart, genome=genome, groups=meth$Day, dmrs=dmrs)
#the region_plot function is available in the best_practices.Rmd file.

Heatmap of DMRs


heatmap_data <- as.data.frame(methrix::get_region_summary(meth, makeGRangesFromDataFrame(dmrs)))
heatmap_data <- heatmap_data[complete.cases(heatmap_data),]
rownames(heatmap_data) <- paste0(heatmap_data$chr,":", heatmap_data$start, "-", heatmap_data$end)
pheatmap::pheatmap(heatmap_data[,-(1:5)], show_rownames = FALSE, annotation_col = as.data.frame(meth@colData), 
                   fontsize = 8)

Number and size of DMRs

if(!requireNamespace("plotly")) {
install.packages("plotly")
}
if(!requireNamespace("ggplot2")) {
install.packages("ggplot2")
}
if(!requireNamespace("scales")) {
install.packages("scales")
}
library(plotly)
library(ggplot2)
library(scales)


count <- c("Higher methylation in day 14" = length(dmrs[dmrs$diff.Methy<0,]),
           "Higher methylation in day 6" = length(dmrs[dmrs$diff.Methy>0,]))
length <- c("Higher methylation in day 14" = sum(width(dmrs[dmrs$diff.Methy<0,])), 
            "Higher methylation in day 6" = sum(width(dmrs[dmrs$diff.Methy<0,]))) 


data <- data.frame(Direction=c("Gain", "Loss"), Count=count, Length=length)

   g <- ggplot(data=data)+geom_col(aes(x=factor(1), y=Count, fill=Direction), position = "dodge")+theme_bw()+theme(axis.title.x = element_blank(), axis.text.x = element_blank())+scale_fill_brewer(palette = "Dark2")+ggtitle("Number of differentially methylated regions")+scale_y_continuous(labels = comma)
ggplotly(g)
   p <- ggplot(data=data)+geom_col(aes(x=factor(1), y=Length, fill=Direction), position = "dodge")+theme_bw()+theme(axis.title.x = element_blank(), axis.text.x = element_blank())+scale_fill_brewer(palette = "Dark2")+ggtitle("Length of differentially methylated regions")+scale_y_continuous(labels = comma)
ggplotly(p)

Region annotation

To describe and interpret the differentially methylated regions, we can use the ChIPseeker package. The location of the DMRs can already give us a hint the involved processes, but once we assign the regions to genes, a pathway enrichment analysis is also possible. For additional details, see the vignette of ChIPseeker:

## loading packages
if(!requireNamespace("ChIPseeker")) {
BiocManager::install("ChIPseeker")
}
if(!requireNamespace("TxDb.Hsapiens.UCSC.hg19.knownGene")) {
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")
}
if(!requireNamespace("ReactomePA")) {
BiocManager::install("ReactomePA")
}
if(!requireNamespace("clusterProfiler")) {
BiocManager::install("clusterProfiler")
}
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
library(ReactomePA)
library(clusterProfiler)

 

  peakAnnoList <- lapply(list("Higher in day 14"=dmrs[dmrs$diff.Methy<0,], "Higher in day 6" = dmrs[dmrs$diff.Methy>0,]), annotatePeak, TxDb=txdb, tssRegion=c(-3000, 3000), verbose=FALSE)

plotAnnoBar(peakAnnoList)
plotDistToTSS(peakAnnoList)



genes <- lapply(list("Higher in day 14"=dmrs[dmrs$diff.Methy<0,], "Higher in day 6" = dmrs[1:3,]), seq2gene, tssRegion = c(-1000, 1000), flankDistance = 3000, TxDb=txdb)

names(genes) <- gsub(" ", "\n", names(genes))
comppthw <- compareCluster(geneCluster   = genes,
                         fun           = "enrichPathway",
                         pvalueCutoff  = 0.05,
                         pAdjustMethod = "BH")
dotplot(comppthw, font.size = 8)

gc()

03: Indentifying DMRs