plot DMR using ggbio (ggplot type graphics, separate pane for every sample with a bar graph for methylation percent)

plotDMR_ggbio.R

library(GenomicRanges)
library(ggbio)
library(ChIPpeakAnno)

gr <- GRanges(seqnames = "NC_035082.1", IRanges(55355800, 55355870), strand = "*")

# Get methylation percents from metilene input and convert to GRanges object
input <- read.table("metilene_Control_Treated.input", sep = "\t", header = T, stringsAsFactors = F)
names(input)[1] <- "chr"
names(input)[2] <- "start"
input$end <- input$start
input <- toGRanges(input)

# Get DMRs and convert to GRanges object
list.files()
dmrs <- read.table("metilene_qval.0.05.out", sep = "\t", header = F, stringsAsFactors = F)
names(dmrs) <- c("chr", "start", "end", "qvalue", "meth_diff", "num_CpGs", "mean_Control", "mean_Selenium")
dmrs <- toGRanges(dmrs)
dmrs

# Create bar graphs showing methylation percent data for every samples
p1 <- ggplot(data.frame(cpgData@ranges, cpgData$Control), aes(x=start, y=cpgData.Control)) + 
                  geom_bar(stat="identity", fill = "dodgerblue") + ylim(0, 100) + ylab("methylation %")
p2 <- ggplot(data.frame(cpgData@ranges, cpgData$Control.1), aes(x=start, y=cpgData.Control.1)) + 
                  geom_bar(stat="identity", fill = "dodgerblue") + ylim(0, 100) + ylab("methylation %")
p3 <- ggplot(data.frame(cpgData@ranges, cpgData$Control.2), aes(x=start, y=cpgData.Control.2)) + 
                  geom_bar(stat="identity", fill = "dodgerblue") + ylim(0, 100) + ylab("methylation %")
p4 <- ggplot(data.frame(cpgData@ranges, cpgData$Selenium), aes(x=start, y=cpgData.Selenium)) + 
                  geom_bar(stat="identity", fill = "firebrick1") + ylim(0, 100) + ylab("methylation %")
p5 <- ggplot(data.frame(cpgData@ranges, cpgData$Selenium.1), aes(x=start, y=cpgData.Selenium.1)) + 
                  geom_bar(stat="identity", fill = "firebrick1") + ylim(0, 100) + ylab("methylation %")
p6 <- autoplot(cpgData) # Plot individual CpG positions

# Use tracks function to put it all together
tracks("Control" = p1, "Control" = p2, "Control" = p3, "Selenium" = p4, "Selenium" = p5, "CpGs" = p6)
ggsave("NC_035082.1_55355800_55355870_DMRplot.tiff", dpi = 200)