### R code from vignette source 'vignettes/MMDiff/inst/doc/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:181-185 ################################################### library('MMDiff') library('MMDiffBamSubset') oldwd <- setwd(system.file("extdata", package="MMDiffBamSubset")) Cfp1 <- dba(sampleSheet="Cfp1.csv", minOverlap=3) ################################################### ### code chunk number 3: MMDiff.Rnw:199-200 ################################################### Cfp1 ################################################### ### code chunk number 4: MMDiff.Rnw:207-209 ################################################### Peaks <- dba.peakset(Cfp1,bRetrieve=TRUE) Peaks <- Peaks[1:1000] ################################################### ### code chunk number 5: MMDiff.Rnw:216-218 ################################################### Peaks.2 <- GRanges(seqnames = Rle('chr1'), ranges = IRanges(start=seq(3200000,3219900,100), width=100)) ################################################### ### code chunk number 6: MMDiff.Rnw:268-272 ################################################### Cfp1Profiles <- getPeakProfiles(Cfp1, Peaks, bin.length=50, save.files=FALSE,run.parallel=FALSE) setwd(oldwd) ################################################### ### code chunk number 7: MMDiff.Rnw:281-282 ################################################### names(Cfp1Profiles$MD) ################################################### ### code chunk number 8: MMDiff.Rnw:293-295 ################################################### PeakRawHists <- Cfp1Profiles$MD$PeakRawHists names(PeakRawHists)[1:10] ################################################### ### code chunk number 9: MMDiff.Rnw:303-305 ################################################### peak.id <- 5 dim(PeakRawHists[[peak.id]]) ################################################### ### code chunk number 10: MMDiff.Rnw:314-316 ################################################### colnames(PeakRawHists[[peak.id]]) rownames(PeakRawHists[[peak.id]]) ################################################### ### code chunk number 11: MMDiff.Rnw:321-324 ################################################### Peak.id <- 33 Sample.ids <- c("WT.AB2", "Null.AB2", "Resc.AB2") plotPeak(Cfp1Profiles, Peak.id, Sample.ids, NormMethod=NULL) ################################################### ### code chunk number 12: MMDiff.Rnw:335-340 ################################################### 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 13: MMDiff.Rnw:370-372 ################################################### SampleIDs <- c("WT.AB2", "Null.AB2", "Resc.AB2") Cfp1Norm <- getNormFactors(Cfp1Profiles, method = "DESeq", SampleIDs=SampleIDs) ################################################### ### code chunk number 14: MMDiff.Rnw:377-378 ################################################### Cfp1Norm$MD$NormFactors$DESeq ################################################### ### code chunk number 15: MMDiff.Rnw:391-394 ################################################### Peak.id <- 33 Sample.ids <- c("WT.AB2", "Null.AB2", "Resc.AB2") plotPeak(Cfp1Norm, Peak.id, Sample.ids, NormMethod='DESeq') ################################################### ### code chunk number 16: MMDiff.Rnw:411-413 ################################################### Cfp1Dists <- compHistDists(Cfp1Norm, method='GMD', overWrite=FALSE, NormMethod='DESeq') ################################################### ### code chunk number 17: MMDiff.Rnw:419-424 ################################################### 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 18: MMDiff.Rnw:429-431 ################################################### Cfp1Dists <- compHistDists(Cfp1Norm, method='GMD', CompIDs=CompIDs, overWrite=FALSE, NormMethod='DESeq') ################################################### ### code chunk number 19: MMDiff.Rnw:440-441 ################################################### Cfp1Dists$MD$DISTS$GMD[1:10,] ################################################### ### code chunk number 20: MMDiff.Rnw:447-449 ################################################### data(Cfp1Dists) names(Cfp1Dists$MD$DISTS) ################################################### ### code chunk number 21: MMDiff.Rnw:463-474 ################################################### 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 22: MMDiff.Rnw:482-484 ################################################### idxMMD <- which(Cfp1Pvals$MD$Pvals$MMD<0.05) ################################################### ### code chunk number 23: MMDiff.Rnw:489-491 ################################################### idxGMD <- which(Cfp1Pvals$MD$Pvals$GMD<0.05) idxPearson <- which(Cfp1Pvals$MD$Pvals$Pearson<0.05) ################################################### ### code chunk number 24: MMDiff.Rnw:495-496 ################################################### rownames(Cfp1Pvals$MD$Pvals$MMD)[idxMMD[1:10]] ################################################### ### code chunk number 25: MMDiff.Rnw:503-509 ################################################### 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 26: setup ################################################### sessionInfo() ################################################### ### code chunk number 27: MMDiff.Rnw:579-582 ################################################### dev.off() unlink(tmp) setwd(savewd)