Nothing
## ----style-Sweave, eval=TRUE, echo=FALSE, results='asis'-------------------
BiocStyle::latex()
## ----eval=FALSE------------------------------------------------------------
# rrbs <- readBedFiles(files, colData, bed_type = 'Encode', eData=NaN)
## ----eval=FALSE------------------------------------------------------------
# files <- c('wgEncodeHaibMethylRrbsH1hescHaibSitesRep1.bed.gz',
# 'wgEncodeHaibMethylRrbsH1hescHaibSitesRep2.bed.gz',
# 'wgEncodeHaibMethylRrbsK562HaibSitesRep1.bed.gz',
# 'wgEncodeHaibMethylRrbsK562HaibSitesRep2.bed.gz')
# group <- factor(c('H1-hESC','H1-hESC','K562','K562'))
# samples <- c('H1-hESC1','H1-hESC2','K562-1','K562-2')
# colData <- DataFrame(group,row.names= samples)
## ----setup,message=FALSE---------------------------------------------------
library(BiSeq)
library(M3D)
## --------------------------------------------------------------------------
data(rrbsDemo)
rrbsDemo
## --------------------------------------------------------------------------
data(CpGsDemo)
CpGsDemo
## --------------------------------------------------------------------------
data(CpGsDemo)
data(rrbsDemo)
overlaps <- findOverlaps(CpGsDemo,rowRanges(rrbsDemo))
## ----eval=FALSE------------------------------------------------------------
# MMDlistDemo <- M3D_Wrapper(rrbsDemo, overlaps)
## --------------------------------------------------------------------------
data(MMDlistDemo)
## --------------------------------------------------------------------------
# the full MMD
head(MMDlistDemo$Full)
# the coverage only MMD
head(MMDlistDemo$Coverage)
## --------------------------------------------------------------------------
M3Dstat <- MMDlistDemo$Full-MMDlistDemo$Coverage
## ----M3DstatPlot-----------------------------------------------------------
repCols <- c(1,6)
plot(as.vector(MMDlistDemo$Full[,repCols]),
as.vector(MMDlistDemo$Coverage[,repCols]),
xlab='Full MMD',ylab='Coverage MMD')
title('Test Statistics: Replicate Comparison')
abline(a=0,b=1,col='red',lwd=2)
## ----M3DstatBetweenPlot----------------------------------------------------
groupCols <- 2:5
plot(as.vector(MMDlistDemo$Full[,groupCols]),
as.vector(MMDlistDemo$Coverage[,groupCols]),
xlab='Full MMD',ylab='Coverage MMD')
title('Test Statistics: Inter-Group Comparison')
abline(a=0,b=1,col='red',lwd=2)
## ----M3DhistPlot-----------------------------------------------------------
repCols <- c(1,6)
groupCols <- 2:5
M3Dstat <- MMDlistDemo$Full - MMDlistDemo$Coverage
breaks <- seq(-0.2,1.2,0.1)
# WT reps
hReps <- hist(M3Dstat[,repCols], breaks=breaks,plot=FALSE)
# mean between groups
hGroups <- hist(rowMeans(M3Dstat[,groupCols]),breaks=breaks,plot=FALSE)
col1 <- "#0000FF40"
col2 <- "#FF000040"
plot(hReps,col=col1, freq=FALSE, ylab='Density',
xlab='MMD statistic', main= 'M3D Stat Histogram')
plot(hGroups, col=col2, freq=FALSE, add=TRUE)
legend(x=0.5, y =3, legend=c('Replicates', 'Groups'), fill=c(col1, col2))
## --------------------------------------------------------------------------
data(PDemo)
## --------------------------------------------------------------------------
group1 <- unique(colData(rrbsDemo)$group)[1]
group2 <-unique(colData(rrbsDemo)$group)[2]
PDemo <- pvals(rrbsDemo, CpGsDemo, M3Dstat,
group1, group2, smaller=FALSE, comparison='allReps', method='empirical', closePara=0.005)
head(PDemo$FDRmean)
## --------------------------------------------------------------------------
called <- which(PDemo$FDRmean<=0.01)
head(called)
head(CpGsDemo[called])
## --------------------------------------------------------------------------
par(mfrow=c(2,2))
for (i in 1:4){
CpGindex <- called[i]
plotMethProfile(rrbsDemo, CpGsDemo, 'H1-hESC', 'K562', CpGindex)
}
## ----eval=FALSE------------------------------------------------------------
# MMDlistDemo <- M3D_Wrapper(rrbsDemo, overlaps)
## ----eval=FALSE------------------------------------------------------------
# MMDlistDemo <- M3D_Wrapper_lite(rrbsDemo, overlaps)
## ----eval=FALSE------------------------------------------------------------
# MMDlistDemo <- M3D_Wrapper(rrbsDemo, overlaps)
# M3Dstat <- MMDlistDemo$Full-MMDlistDemo$Coverage
# PDemo <- pvals(rrbsDemo, CpGsDemo, M3Dstat,
# group1, group2, smaller=FALSE, comparison='allReps', method='empirical', closePara=0.005)
## ----eval=FALSE------------------------------------------------------------
# MMDlistDemo <- M3D_Wrapper_lite(rrbsDemo, overlaps)
# PDemo <- pvals_lite(rrbsDemo, CpGsDemo, M3Dstat,
# group1, group2, smaller=FALSE, comparison='allReps', method='empirical', closePara=0.005)
## ----eval=FALSE------------------------------------------------------------
# # change working directory to the location of the files
# group <- factor(c('H1-hESC','H1-hESC','K562','K562'))
# samples <- c('H1-hESC1','H1-hESC2','K562-1','K562-2')
# colData <- DataFrame(group,row.names= samples)
# files <- c('wgEncodeHaibMethylRrbsH1hescHaibSitesRep1.bed.gz',
# 'wgEncodeHaibMethylRrbsH1hescHaibSitesRep2.bed.gz',
# 'wgEncodeHaibMethylRrbsK562HaibSitesRep1.bed.gz',
# 'wgEncodeHaibMethylRrbsK562HaibSitesRep2.bed.gz')
# rrbs <- readBedFiles(files,colData, bed_type = ENCODE)
## ----eval=FALSE------------------------------------------------------------
# # Create the CpGs
# rrbs.CpGs <- clusterSites(object=rrbs,groups=unique(colData(rrbs)$group),
# perc.samples = 3/4, min.sites = 20, max.dist=100)
# CpGs <- clusterSitesToGR(rrbs.CpGs)
## ----eval=FALSE------------------------------------------------------------
# inds <- vector()
# overlaps <- findOverlaps(CpGs,rowRanges(rrbs.CpGs))
# for (i in 1:length(CpGs)){
# covs <- colSums(totalReads(rrbs.CpGs)[subjectHits(
# overlaps[queryHits(overlaps)==i]),])
# if (!any(covs<=100)){
# inds <- c(inds,i)
# }
# }
# CpGs <- CpGs[inds]
# rm(inds)
## ----eval=FALSE------------------------------------------------------------
# # reduce the rrbs to only the cytosine sites within regions
# rrbs <- rrbs.CpGs
# CpGs <- CpGs[1:1000] # first 1000 regions
# overlaps <- findOverlaps(CpGs,rowRanges(rrbs))
# rrbs <- rrbs[subjectHits(overlaps)]
# overlaps <- findOverlaps(CpGs,rowRanges(rrbs))
## --------------------------------------------------------------------------
sessionInfo()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.