### R code from vignette source 'MMDiff.Rnw' ################################################### ### code chunk number 1: MMDiff.Rnw:142-147 ################################################### tmp <- tempfile(as.character(Sys.getpid())) pdf(tmp) savewarn <- options("warn") options(warn=-1) savewd <- getwd() ################################################### ### code chunk number 2: MMDiff.Rnw:180-188 ################################################### MMDiffBS_extdata <- system.file("extdata", package="MMDiffBamSubset") leftoverPlots <- dir(MMDiffBS_extdata, pattern=glob2rx("Rplots*.pdf")) if(length(leftoverPlots) > 0) { try({ file.remove(paste(MMDiffBS_extdata, leftoverPlots, sep=.Platform$file.sep)) }, silent=TRUE) } ################################################### ### code chunk number 3: MMDiff.Rnw:191-195 ################################################### library('MMDiff') library('MMDiffBamSubset') oldwd <- setwd(system.file("extdata", package="MMDiffBamSubset")) Cfp1 <- dba(sampleSheet="Cfp1.csv", minOverlap=3) ################################################### ### code chunk number 4: MMDiff.Rnw:209-210 ################################################### Cfp1 ################################################### ### code chunk number 5: MMDiff.Rnw:217-219 ################################################### Peaks <- dba.peakset(Cfp1,bRetrieve=TRUE) Peaks <- Peaks[1:1000] ################################################### ### code chunk number 6: MMDiff.Rnw:226-228 ################################################### Peaks.2 <- GRanges(seqnames = Rle('chr1'), ranges = IRanges(start=seq(3200000,3219900,100), width=100)) ################################################### ### code chunk number 7: MMDiff.Rnw:278-282 ################################################### Cfp1Profiles <- getPeakProfiles(Cfp1, Peaks, bin.length=50, save.files=FALSE,run.parallel=FALSE) setwd(oldwd) ################################################### ### code chunk number 8: MMDiff.Rnw:291-292 ################################################### names(Cfp1Profiles$MD) ################################################### ### code chunk number 9: MMDiff.Rnw:303-305 ################################################### PeakRawHists <- Cfp1Profiles$MD$PeakRawHists names(PeakRawHists)[1:10] ################################################### ### code chunk number 10: MMDiff.Rnw:313-315 ################################################### peak.id <- 5 dim(PeakRawHists[[peak.id]]) ################################################### ### code chunk number 11: MMDiff.Rnw:324-326 ################################################### colnames(PeakRawHists[[peak.id]]) rownames(PeakRawHists[[peak.id]]) ################################################### ### code chunk number 12: MMDiff.Rnw:331-334 ################################################### Peak.id <- 33 Sample.ids <- c("WT.AB2", "Null.AB2", "Resc.AB2") plotPeak(Cfp1Profiles, Peak.id, Sample.ids, NormMethod=NULL) ################################################### ### code chunk number 13: MMDiff.Rnw:345-350 ################################################### oldwd <- setwd(system.file("extdata", package="MMDiffBamSubset")) Cfp1Profiles2 <- getPeakProfiles(Cfp1, Peaks, bin.length=50, save.files=FALSE, keep.extra=TRUE) names(Cfp1Profiles2$MD$RawHists$WT_2) setwd(oldwd) ################################################### ### code chunk number 14: MMDiff.Rnw:380-382 ################################################### SampleIDs <- c("WT.AB2", "Null.AB2", "Resc.AB2") Cfp1Norm <- getNormFactors(Cfp1Profiles, method = "DESeq", SampleIDs=SampleIDs) ################################################### ### code chunk number 15: MMDiff.Rnw:387-388 ################################################### Cfp1Norm$MD$NormFactors$DESeq ################################################### ### code chunk number 16: MMDiff.Rnw:401-404 ################################################### Peak.id <- 33 Sample.ids <- c("WT.AB2", "Null.AB2", "Resc.AB2") plotPeak(Cfp1Norm, Peak.id, Sample.ids, NormMethod='DESeq') ################################################### ### code chunk number 17: MMDiff.Rnw:421-423 ################################################### Cfp1Dists <- compHistDists(Cfp1Norm, method='GMD', overWrite=FALSE, NormMethod='DESeq') ################################################### ### code chunk number 18: MMDiff.Rnw:429-434 ################################################### SampleIDs <- Cfp1Norm$samples$SampleID ID <- which(upper.tri(matrix(1, length(SampleIDs), length(SampleIDs)), diag = F), arr.ind=T) CompIDs <- rbind(SampleIDs[ID[,1]], SampleIDs[ID[,2]]) CompIDs ################################################### ### code chunk number 19: MMDiff.Rnw:439-441 ################################################### Cfp1Dists <- compHistDists(Cfp1Norm, method='GMD', CompIDs=CompIDs, overWrite=FALSE, NormMethod='DESeq') ################################################### ### code chunk number 20: MMDiff.Rnw:450-451 ################################################### Cfp1Dists$MD$DISTS$GMD[1:10,] ################################################### ### code chunk number 21: MMDiff.Rnw:457-459 ################################################### data(Cfp1Dists) names(Cfp1Dists$MD$DISTS) ################################################### ### code chunk number 22: MMDiff.Rnw:473-484 ################################################### data(Cfp1Dists) group1 <- c("WT.AB2","Resc.AB2") group2 <- c("Null.AB2") Cfp1Pvals <- detPeakPvals(Cfp1Dists, group1=group1, group2=group2, name1='Wt/Resc', name2='Null', method='MMD') Cfp1Pvals <- detPeakPvals(Cfp1Pvals, group1=group1, group2=group2, name1='Wt/Resc', name2='Null', method='GMD') Cfp1Pvals <- detPeakPvals(Cfp1Pvals, group1=group1, group2=group2, name1='Wt/Resc', name2='Null', method='Pearson') ################################################### ### code chunk number 23: MMDiff.Rnw:492-494 ################################################### idxMMD <- which(Cfp1Pvals$MD$Pvals$MMD<0.05) ################################################### ### code chunk number 24: MMDiff.Rnw:499-501 ################################################### idxGMD <- which(Cfp1Pvals$MD$Pvals$GMD<0.05) idxPearson <- which(Cfp1Pvals$MD$Pvals$Pearson<0.05) ################################################### ### code chunk number 25: MMDiff.Rnw:505-506 ################################################### rownames(Cfp1Pvals$MD$Pvals$MMD)[idxMMD[1:10]] ################################################### ### code chunk number 26: MMDiff.Rnw:513-519 ################################################### group1 <- c("WT.AB2", "Resc.AB2") group2 <- c("Null.AB2") plotHistDists(Cfp1Pvals, group1=group1, group2=group2, method='MMD') plotHistDists(Cfp1Pvals, group1=group1, group2=group2, method='GMD') plotHistDists(Cfp1Pvals, group1=group1, group2=group2, method='Pearson') ################################################### ### code chunk number 27: setup ################################################### sessionInfo() ################################################### ### code chunk number 28: MMDiff.Rnw:589-592 ################################################### dev.off() unlink(tmp) setwd(savewd)