## ----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()