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