Here is the tutorial (Analyzing WGBS data with bsseq)
I want to use this method to visualize differential dmr regions on the same gene.
Here is my code：
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("bsseq")
BiocManager::install("bsseqData")
library(bsseq)
library(bsseqData)
data(BS.cancer.ex.fit)
BS.cancer.ex.fit <- updateObject(BS.cancer.ex.fit)
BS.cancer.ex.fit
class(BS.cancer.ex.fit)
##
data(BS.cancer.ex.fit)
# BS.cancer.ex.fit <- updateObject(BS.cancer.ex.fit)
# BS.cancer.ex.fit
##
BS.cov <- getCoverage(BS.cancer.ex.fit)
keepLoci.ex <- which(rowSums(BS.cov[, BS.cancer.ex$Type == "cancer"] >= 2) >= 2 &
rowSums(BS.cov[, BS.cancer.ex$Type == "normal"] >= 2) >= 2)
length(keepLoci.ex)
# [1] 597371
BS.cancer.ex.fit <- BS.cancer.ex.fit[keepLoci.ex,]
length(BS.cancer.ex.fit)
## (the keepLoci.ex is also available for direct inspection in the bsseqData package.)
## We are now ready to compute t-statistics, by
BS.cancer.ex.tstat <- BSmooth.tstat(BS.cancer.ex.fit,
group1 = c("C1", "C2", "C3"),
group2 = c("N1", "N2", "N3"),
estimate.var = "group2",
local.correct = TRUE,
verbose = TRUE)
# BS.cancer.ex.tstat
# plot(BS.cancer.ex.tstat)
#####
dmrs0 <- dmrFinder(BS.cancer.ex.tstat, cutoff = c(-4.6, 4.6))
dmrs <- subset(dmrs0, n >= 3 & abs(meanDiff) >= 0.1)
nrow(dmrs)
head(dmrs, n = 3)
######
pData <- pData(BS.cancer.ex.fit)
pData$col <- rep(c("red", "blue"), each = 3)
pData(BS.cancer.ex.fit) <- pData
## Once this is setup, we can plot a single DMR like
plotRegion(BS.cancer.ex.fit, dmrs[1,], extend = 5000, addRegions = dmrs)
I successfully got the figure in the end, but my question is I think a little information is missing on my diagram.
Below is my figure produced by the code above:
The perfect picture below has the gene bar information. I don't know how to add this bar into my plot. Here is the article about this picture (figure 3b).
I don't think the gene bar below the picture is added manually. Because there are other precise information like promoter, exon and methylation regions.
Just like the description in the figure legend: below plots, gene diagrams, showing gene body (black bars), annotated promoter region (red bar), exons (blue bars), location of individual CpG sites (‘tick marks’) and gene direction (arrow).
I hope somebody could give me some advice or the right way to solve my problem. Very thankful.
