Hi,
Following code has been extracted from Daniel Hansen Lab Github course on 450K analysis (https://github.com/hansenlab/tutorial.450k/blob/master/vignettes/methylation450k.Rmd). I hope it helps.
Plotting DMRs
As with DMPs, it's always a good idea to plot DMRs. r Biocpkg("minfi")
does not currently include any functionality for doing this directly but will soon[^soon]. In the mean time, we can (with a little bit of a work) make a publication-quality figure using the r Biocpkg("Gviz")
.
[^soon]: Famous last words ...
NOTE: This figure takes a while to generate because r Biocpkg("Gviz")
needs to first download some data.
library(Gviz)
genome <- "hg19"
# NOTE: Using the non-permuted
dmr <- bh_dmrs$table[1, ]
chrom <- dmr$chr
start <- dmr$start
end <- dmr$end
minbase <- start - 0.25 * (end - start)
maxbase <- end + 0.25 * (end - start)
pal <- c("#E41A1C", "#377EB8")
# Start building the tracks
iTrack <- IdeogramTrack(genome = genome,
chromosome = dmr$chr,
name = "")
gTrack <- GenomeAxisTrack(col = "black",
cex = 1,
name = "",
fontcolor = "black")
# NOTE: This track takes a little while to create
rTrack <- UcscTrack(genome = genome,
chromosome = chrom,
track = "refGene",
from = minbase,
to = maxbase,
trackType = "GeneRegionTrack",
rstarts = "exonStarts",
rends = "exonEnds",
gene = "name",
symbol = "name2",
transcript = "name",
strand = "strand",
fill = "darkblue",
stacking = "squish",
name = "RefSeq",
showId = TRUE,
geneSymbol = TRUE)
# methylation data track
gr <- granges(GRset.funnorm)
gr$beta <- getBeta(GRset.funnorm)
methTrack <- DataTrack(range = gr,
groups = targets$status,
genome = genome,
chromosome = chrom,
ylim = c(-0.05, 1.05),
col = pal,
type = c("a","p"),
name = "DNA Meth.\n(beta value)",
background.panel = "white",
legend = TRUE,
cex.title = 0.8,
cex.axis = 0.8,
cex.legend = 0.8)
# DMR position data track
dmrTrack <- AnnotationTrack(start = start,
end = end,
genome = genome,
name = "DMR",
chromosom = chrom)
# Finally, plot the tracks
tracks <- list(iTrack, gTrack, methTrack, dmrTrack, rTrack)
sizes <- c(2, 2, 5, 2, 3) # set up the relative sizes of the tracks
plotTracks(tracks,
from = minbase,
to = maxbase,
showTitle = TRUE,
add53 = TRUE,
add35 = TRUE,
grid = TRUE,
lty.grid = 3,
sizes = sizes,
length(tracks))