DiffBind.R
In DiffBind: Differential Binding Analysis of ChIP-Seq Peak Data

### R code from vignette source 'DiffBind.Rnw'

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: DiffBind.Rnw:222-226
###################################################
tmp <-  tempfile(as.character(Sys.getpid()))
pdf(tmp)
savewarn <- options("warn")
options(warn=-1)


###################################################
### code chunk number 3: DiffBind.Rnw:229-230
###################################################
library(DiffBind)


###################################################
### code chunk number 4: DiffBind.Rnw:232-233 (eval = FALSE)
###################################################
## setwd(system.file('extra',package='DiffBind'))


###################################################
### code chunk number 5: DiffBind.Rnw:239-240 (eval = FALSE)
###################################################
## tamoxifen <- dba.analyze("tamoxifen.csv")


###################################################
### code chunk number 6: DiffBind.Rnw:246-247 (eval = FALSE)
###################################################
## tamoxifen.DB <- dba.report(tamoxifen)


###################################################
### code chunk number 7: DiffBind.Rnw:258-264 (eval = FALSE)
###################################################
## tamoxifen <- dba(sampleSheet="tamoxifen.csv") %>%
##   dba.blacklist() %>%
##   dba.count()     %>%
##   dba.normalize() %>%
##   dba.contrast()  %>%
##   dba.analyze()


###################################################
### code chunk number 8: sampSheet
###################################################
samples <- read.csv(file.path(system.file("extra", package="DiffBind"),
                              "tamoxifen.csv"))
names(samples)
samples


###################################################
### code chunk number 9: dbaConstruct
###################################################
tamoxifen <- dba(sampleSheet="tamoxifen.csv",
                 dir=system.file("extra", package="DiffBind"))


###################################################
### code chunk number 10: dbaConstructDF (eval = FALSE)
###################################################
## tamoxifen <- dba(sampleSheet=samples)


###################################################
### code chunk number 11: DiffBind.Rnw:305-306
###################################################
tamoxifen


###################################################
### code chunk number 12: tamox_occ_corhm
###################################################
plot(tamoxifen)


###################################################
### code chunk number 13: DiffBind.Rnw:358-359 (eval = FALSE)
###################################################
## tamoxifen <- dba.count(tamoxifen) 


###################################################
### code chunk number 14: DiffBind.Rnw:361-362
###################################################
data(tamoxifen_counts)


###################################################
### code chunk number 15: DiffBind.Rnw:370-371
###################################################
tamoxifen


###################################################
### code chunk number 16: eff_lib_size
###################################################
info <- dba.show(tamoxifen)
libsizes <- cbind(LibReads=info$Reads, FRiP=info$FRiP, 
                  PeakReads=round(info$Reads * info$FRiP))
rownames(libsizes) <- info$ID
libsizes


###################################################
### code chunk number 17: tamox_aff_corhm
###################################################
plot(tamoxifen)


###################################################
### code chunk number 18: normalize
###################################################
tamoxifen <- dba.normalize(tamoxifen)


###################################################
### code chunk number 19: show_norm
###################################################
norm <- dba.normalize(tamoxifen, bRetrieve=TRUE)
norm


###################################################
### code chunk number 20: norm_facs
###################################################
normlibs <- cbind(FullLibSize=norm$lib.sizes, NormFacs=norm$norm.factors,
                  NormLibSize=round(norm$lib.sizes/norm$norm.factors))
rownames(normlibs) <- info$ID
normlibs


###################################################
### code chunk number 21: DiffBind.Rnw:481-484
###################################################
tamoxifen <- dba.contrast(tamoxifen, 
                          reorderMeta=list(Condition="Responsive"))
tamoxifen


###################################################
### code chunk number 22: DiffBind.Rnw:516-518
###################################################
tamoxifen <- dba.analyze(tamoxifen)
dba.show(tamoxifen, bContrasts=TRUE)


###################################################
### code chunk number 23: tamox_sdb_corhm
###################################################
plot(tamoxifen, contrast=1)


###################################################
### code chunk number 24: DiffBind.Rnw:568-569
###################################################
tamoxifen.DB <- dba.report(tamoxifen)


###################################################
### code chunk number 25: DiffBind.Rnw:575-576
###################################################
tamoxifen.DB


###################################################
### code chunk number 26: DiffBind.Rnw:596-598
###################################################
sum(tamoxifen.DB$Fold>0)
sum(tamoxifen.DB$Fold<0)


###################################################
### code chunk number 27: tamox_sdb_venn
###################################################
dba.plotVenn(tamoxifen, contrast=1, bDB=TRUE,
             bGain=TRUE, bLoss=TRUE, bAll=FALSE)


###################################################
### code chunk number 28: tamox_aff_pca
###################################################
dba.plotPCA(tamoxifen,DBA_TISSUE,label=DBA_CONDITION)


###################################################
### code chunk number 29: tamox_sdb_pca
###################################################
dba.plotPCA(tamoxifen, contrast=1, label=DBA_TISSUE)


###################################################
### code chunk number 30: tamox_sdb_ma
###################################################
dba.plotMA(tamoxifen)


###################################################
### code chunk number 31: tamox_sdb_volcano
###################################################
dba.plotVolcano(tamoxifen)


###################################################
### code chunk number 32: DiffBind.Rnw:767-769
###################################################
sum(tamoxifen.DB$Fold<0)
sum(tamoxifen.DB$Fold>0)


###################################################
### code chunk number 33: tamox_sdb_box
###################################################
pvals <- dba.plotBox(tamoxifen)


###################################################
### code chunk number 34: DiffBind.Rnw:808-809
###################################################
pvals


###################################################
### code chunk number 35: DiffBind.Rnw:828-829
###################################################
corvals <- dba.plotHeatmap(tamoxifen)


###################################################
### code chunk number 36: tamox_sdb_hm
###################################################
hmap <- colorRampPalette(c("red", "black", "green"))(n = 13)
readscores <- dba.plotHeatmap(tamoxifen, contrast=1, correlations=FALSE,
                              scale="row", colScheme = hmap)


###################################################
### code chunk number 37: DiffBind.Rnw:908-909
###################################################
tamoxifen <- dba.contrast(tamoxifen,design="~Tissue + Condition")


###################################################
### code chunk number 38: DiffBind.Rnw:915-917
###################################################
tamoxifen <- dba.analyze(tamoxifen)
dba.show(tamoxifen, bContrasts=TRUE)


###################################################
### code chunk number 39: tamox_block_ma
###################################################
dba.plotMA(tamoxifen)


###################################################
### code chunk number 40: tamox_block_vol
###################################################
dba.plotVolcano(tamoxifen)


###################################################
### code chunk number 41: DiffBind.Rnw:966-967
###################################################
multifactor.DB <- dba.report(tamoxifen)


###################################################
### code chunk number 42: DiffBind.Rnw:974-976
###################################################
min(abs(tamoxifen.DB$Fold))
min(abs(multifactor.DB$Fold))


###################################################
### code chunk number 43: DiffBind.Rnw:984-986
###################################################
sum(tamoxifen.DB$Fold > 0) / sum(tamoxifen.DB$Fold < 0)
sum(multifactor.DB$Fold > 0) / sum(multifactor.DB$Fold < 0)


###################################################
### code chunk number 44: DiffBind.Rnw:992-994
###################################################
tamoxifen <- dba.analyze(tamoxifen,method=DBA_ALL_METHODS)
dba.show(tamoxifen,bContrasts=TRUE)


###################################################
### code chunk number 45: tamox_block_venn
###################################################
tamoxifen.OL <- dba.plotVenn(tamoxifen,contrast=1,method=DBA_ALL_METHODS,
                             bDB=TRUE)


###################################################
### code chunk number 46: DiffBind.Rnw:1034-1038
###################################################
tamoxifen <- dba.contrast(tamoxifen,contrast=c("Tissue","MCF7","T47D"),
                          reorderMeta = list(Tissue="MCF7"))
tamoxifen <- dba.analyze(tamoxifen,method=DBA_ALL_METHODS)
dba.show(tamoxifen,bContrasts=TRUE)


###################################################
### code chunk number 47: blacklistPeaks
###################################################
data(tamoxifen_peaks)
tamoxifen
peakdata  <- dba.show(tamoxifen)$Intervals
tamoxifen <- dba.blacklist(tamoxifen, blacklist=DBA_BLACKLIST_HG19, 
                           greylist=FALSE)
tamoxifen
peakdata.BL <- dba.show(tamoxifen)$Intervals
peakdata - peakdata.BL


###################################################
### code chunk number 48: blacklistAnal
###################################################
length(multifactor.DB)
data(tamoxifen_counts)
tamoxifen   <- dba.blacklist(tamoxifen, blacklist=DBA_BLACKLIST_HG19,
                             greylist=FALSE)
blacklisted <- dba.blacklist(tamoxifen, Retrieve=DBA_BLACKLISTED_PEAKS)
tamoxifen   <- dba.contrast(tamoxifen, design="~Tissue + Condition")
tamoxifen   <- dba.analyze(tamoxifen)
blacklisted.DB <- dba.report(tamoxifen)
length(blacklisted.DB)


###################################################
### code chunk number 49: blacklistRes
###################################################
bl_site <- match(blacklisted[[1]], multifactor.DB)
multifactor.DB[bl_site,]
is.na(match(blacklisted[[1]], blacklisted.DB))


###################################################
### code chunk number 50: greylistGet
###################################################
data(tamoxifen_greylist)
names(tamoxifen.greylist)
tamoxifen.greylist$master


###################################################
### code chunk number 51: greylistControls
###################################################
names(tamoxifen.greylist$controls)
tamoxifen.greylist$controls


###################################################
### code chunk number 52: greylistPeaks
###################################################
data(tamoxifen_peaks)
tamoxifen <- dba.blacklist(tamoxifen, blacklist=DBA_BLACKLIST_HG19,
                           greylist=tamoxifen.greylist)


###################################################
### code chunk number 53: greylistCons
###################################################
data(tamoxifen_counts)
cons.peaks <- dba.show(tamoxifen)$Intervals[1]
tamoxifen  <- dba.blacklist(tamoxifen, blacklist=DBA_BLACKLIST_HG19,
                            greylist=tamoxifen.greylist)
cons.peaks.grey <- dba.show(tamoxifen)$Intervals[1]
cons.peaks - cons.peaks.grey


###################################################
### code chunk number 54: greylistMake (eval = FALSE)
###################################################
## tamoxifen <- dba(sampleSheet="tamoxifen.csv")
## tamoxifen <- dba.blacklist(tamoxifen)
## tamoxifen.greylist <- dba.blacklist(tamoxifen, Retrieve=DBA_GREYLIST)


###################################################
### code chunk number 55: greylistP
###################################################
tamoxifen$config$greylist.pval <- 0.999


###################################################
### code chunk number 56: norm0
###################################################
data(tamoxifen_analysis)
dba.plotMA(tamoxifen, contrast=list(Resistant=tamoxifen$masks$Resistant),
           bNormalized=FALSE, sub="Non-Normalized")


###################################################
### code chunk number 57: normDESeq2LibFull
###################################################
tamoxifen <- dba.normalize(tamoxifen, normalize=DBA_NORM_LIB)
tamoxifen <- dba.analyze(tamoxifen)
dba.plotMA(tamoxifen, method=DBA_DESEQ2, sub="DESeq2:lib:full")


###################################################
### code chunk number 58: normRes1
###################################################
dbs <- dba.report(tamoxifen, bDB=TRUE, bGain=TRUE, bLoss=TRUE)
dbs$config$factor <- "normalize"
dbs$class[DBA_ID,] <- colnames(dbs$class)[1] <-  "LIB_Full"
dbs$class[DBA_FACTOR,] <- DBA_NORM_LIB
dbs


###################################################
### code chunk number 59: normDESeq2RLE
###################################################
tamoxifen <- dba.normalize(tamoxifen, normalize=DBA_NORM_NATIVE)
tamoxifen <- dba.analyze(tamoxifen)
dba.plotMA(tamoxifen, method=DBA_DESEQ2, sub="DESeq2:RLE:RiP")


###################################################
### code chunk number 60: normRes2
###################################################
db <- dba.report(tamoxifen, bDB=TRUE, bGain=TRUE, bLoss=TRUE)
db$class[DBA_ID,] <- "RLE_RiP"
db$class[DBA_FACTOR,] <- DBA_NORM_RLE
dbs <- dba.peakset(dbs, db)
db


###################################################
### code chunk number 61: normDESeq2Comparison
###################################################
par(mfrow=c(3,1))
dba.plotVenn(dbs,c(1,4), main="Total DB Sites")
dba.plotVenn(dbs,dbs$masks$Gain,main="Gain in Resistant")
dba.plotVenn(dbs,dbs$masks$Loss,main="Gain in Responsive")
par(mfrow=c(1,1))


###################################################
### code chunk number 62: normDESeq2LibRiP
###################################################
tamoxifen <- dba.normalize(tamoxifen, normalize=DBA_NORM_LIB,
                           library=DBA_LIBSIZE_PEAKREADS, background=FALSE)
tamoxifen <- dba.analyze(tamoxifen)
dba.plotMA(tamoxifen, method=DBA_DESEQ2, sub="DESeq2:lib:RiP")


###################################################
### code chunk number 63: normRes3
###################################################
dbs$class[DBA_CONDITION,1:3] <- DBA_LIBSIZE_FULL
dbs$class[DBA_CONDITION,4:6] <- DBA_LIBSIZE_PEAKREADS
dbs$config$condition <- "lib.size"
db <- dba.report(tamoxifen, bDB=TRUE, bGain=TRUE, bLoss=TRUE)
db$class[DBA_ID,] <- "LIB_RiP"
db$class[DBA_FACTOR,] <- DBA_NORM_LIB
db$class[DBA_CONDITION,] <- DBA_LIBSIZE_PEAKREADS
dbs <- dba.peakset(dbs, db)
db


###################################################
### code chunk number 64: normDESeq2LibsizeVenn
###################################################
dba.plotVenn(dbs,c(1,7,4),main="DB Sites")


###################################################
### code chunk number 65: frip
###################################################
dba.show(tamoxifen,attributes=c(DBA_ID,DBA_FRIP))


###################################################
### code chunk number 66: normBG
###################################################
data(tamoxifen_analysis)
tamoxifen <- dba.normalize(tamoxifen, method=DBA_ALL_METHODS,
                           normalize=DBA_NORM_NATIVE,
                           background=TRUE)
tamoxifen <- dba.analyze(tamoxifen, method=DBA_ALL_METHODS)
dba.show(tamoxifen,bContrasts=TRUE)
par(mfrow=c(2,1))
dba.plotMA(tamoxifen, method=DBA_EDGER, sub="edgeR:TMM:background")
dba.plotMA(tamoxifen, method=DBA_DESEQ2, sub="DESeq2:RLE:background")
par(mfrow=c(1,1))


###################################################
### code chunk number 67: normBgRes
###################################################
db <- dba.report(tamoxifen, bDB=TRUE, bGain=TRUE, bLoss=TRUE)
db$class[DBA_ID,] <- "RLE_BG"
db$class[DBA_FACTOR,] <- DBA_NORM_RLE
db$class[DBA_CONDITION,] <- DBA_LIBSIZE_BACKGROUND
dbs <- dba.peakset(dbs, db)
db


###################################################
### code chunk number 68: normBGVenns
###################################################
par(mfcol=c(3,2))
dba.plotVenn(dbs,c(1,10),   main="All Differentially Bound Sites")
dba.plotVenn(dbs,c(2,11),   main="Gain in Resistant cells")
dba.plotVenn(dbs,c(3,12),   main="Loss in Resistant cells")
dba.plotVenn(dbs,c(1,10,4), main="All Differentially Bound Sies")
dba.plotVenn(dbs,c(2,11,5), main="Gain in Resistant cells")
dba.plotVenn(dbs,c(3,12,6), main="Loss in Resistant cells")


###################################################
### code chunk number 69: MCF7T47D
###################################################
mcf7t47d <- dba(tamoxifen,mask=c(3:7))
dba.plotMA(mcf7t47d,
           contrast=list(MCF7=mcf7t47d$masks$MCF7,
                         T47D=mcf7t47d$masks$T47D), 
           bNormalized=FALSE)


###################################################
### code chunk number 70: loessfit
###################################################
mcf7t47d$config$AnalysisMethod <- DBA_EDGER
mcf7t47d <- dba.normalize(mcf7t47d, offsets=TRUE)
mcf7t47d <- dba.contrast(mcf7t47d, contrast=c("Tissue","MCF7","T47D"))
mcf7t47d <- dba.analyze(mcf7t47d)
dba.plotMA(mcf7t47d)


###################################################
### code chunk number 71: compareLoess
###################################################
mcf7t47d.DB <- dba.report(mcf7t47d)
sum(mcf7t47d.DB$Fold > 0)
sum(mcf7t47d.DB$Fold < 0)


###################################################
### code chunk number 72: normCompare
###################################################
data(tamoxifen_analysis)
dbs.all <- NULL
for(norm in c("lib","RLE","TMM", "loess")) {
  for(libsize in c("full","RiP","background")) {
    tam <- NULL
    background <- offsets <- FALSE
    if(libsize == "full" && norm != "lib") {
      background <- NULL
    }
    if(libsize == DBA_LIBSIZE_BACKGROUND) {
      background <- TRUE
      if(norm == DBA_NORM_LIB) {
        background <- NULL
      }
    } 
    if(norm == "loess" && !is.null(background)) {
      offsets <- TRUE
      if(libsize != "background") {
        background <- FALSE
      } else {
        background <- NULL
      }
    }
    if(!is.null(background)) {
      tam <- dba.normalize(tamoxifen, method=DBA_ALL_METHODS,
                           normalize=norm, library=libsize, 
                           background=background, offsets=offsets)
    }
    if(!is.null(tam)) {
      tam <- dba.analyze(tam, method=DBA_ALL_METHODS)
      for(meth in DBA_ALL_METHODS) {
        db <- dba.report(tam, method=meth, bDB=TRUE)
        if(meth == DBA_EDGER) {
          methstr <- "edgeR"
        } else {
          methstr <- "DESeq2"
        }
        if(libsize == "background") {
          libstr <- "BG"
        } else {
          libstr <- libsize
        }
        id <- paste(norm,libstr,methstr,sep="_")
        if(libsize == "full") {
          libstr <- "BG"
        }
        if(is.null(dbs.all)) {
          dbs.all <- db
          dbs.all$config$factor    <- "Normalization Method"
          dbs.all$config$condition <- "Reference Reads"
          dbs.all$config$treatment <- "Analyis Method"
          dbs.all$class[DBA_ID,]  <- colnames(dbs.all$class)[1] <- id
          
          dbs.all$class[DBA_FACTOR,]    <- norm
          dbs.all$class[DBA_CONDITION,] <- libstr
          dbs.all$class[DBA_TREATMENT,] <- "edgeR"
          dbs.all$class[DBA_TISSUE,]    <- NA
        } else {
          db$class[DBA_ID,]        <- id
          db$class[DBA_FACTOR,]    <- norm
          db$class[DBA_CONDITION,] <- libstr
          db$class[DBA_TREATMENT,] <- methstr
          db$class[DBA_TISSUE,]    <- NA
          dbs.all <- dba.peakset(dbs.all,db)
        }
      }
    }
  }
}
dbs.all <- dba(dbs.all,minOverlap=1)
dbs.all


###################################################
### code chunk number 73: normClusterDESeq2
###################################################
deseq <- dba(dbs.all,mask=dbs.all$masks$DESeq2, minOverlap = 1)
dba.plotHeatmap(deseq, maxSites=nrow(deseq$binding),  bLog=FALSE,
                correlations=FALSE,minval=-5, maxval=5, cexCol=1.3, 
                colScheme = hmap, main = "DESeq2 Differentially Bound Sites", 
                ColAttributes = c(DBA_CONDITION, DBA_FACTOR),
                key.title = "LFC")


###################################################
### code chunk number 74: normCluster1
###################################################
dba.plotHeatmap(dbs.all, maxSites=nrow(dbs.all$binding), bLog=FALSE,
                correlations=FALSE, minval=-5, maxval=5, cexCol=1.3, 
                colScheme = hmap, key.title="LFC",
                ColAttributes = c(DBA_CONDITION, DBA_TREATMENT, DBA_FACTOR),
                main="All Differentially Bound Sites")


###################################################
### code chunk number 75: normCluster2
###################################################
dba.plotHeatmap(dbs.all, cexCol=1.3, main="Correlations of DB Sites",
                ColAttributes = c(DBA_CONDITION, DBA_TREATMENT, DBA_FACTOR))


###################################################
### code chunk number 76: loadSpikes
###################################################
load(system.file('extra/spikes.rda',package='DiffBind'))
spikes


###################################################
### code chunk number 77: spikeins
###################################################
spikes$samples$Spikein


###################################################
### code chunk number 78: spikeMAnone
###################################################
dba.plotMA(spikes, contrast=list(Fulvestrant=spikes$masks$Fulvestrant),
           bNormalized=FALSE, sub="RAW", bSmooth=FALSE, dotSize=1.5)


###################################################
### code chunk number 79: spikeMAstraight
###################################################
par(mfrow=c(3,1))
spikes <- dba.normalize(spikes, normalize=DBA_NORM_LIB, 
                        background=FALSE)
spikes <- dba.analyze(spikes)
dba.plotMA(spikes, sub="LIB full", bSmooth=FALSE, dotSize=1.5)

spikes <- dba.normalize(spikes, normalize="RLE",
                        background=FALSE)
spikes <- dba.analyze(spikes)
dba.plotMA(spikes, sub="RLE RiP", bSmooth=FALSE, dotSize=1.5)

spikes <- dba.normalize(spikes, normalize="RLE", 
                        background=TRUE)
spikes <- dba.analyze(spikes)
dba.plotMA(spikes, sub="RLE BG", bSmooth=FALSE, dotSize=1.5)
par(mfrow=c(1,1))


###################################################
### code chunk number 80: spikeMAspikein
###################################################
par(mfrow=c(2,1))
spikes <- dba.normalize(spikes, spikein = spikes.spikeins)
spikes <- dba.analyze(spikes)
dba.plotMA(spikes, sub="LIB spikein", bSmooth=FALSE, dotSize=1.5)

spikes <- dba.normalize(spikes, normalize=DBA_NORM_RLE, spikein = TRUE)
spikes <- dba.analyze(spikes)
dba.plotMA(spikes, sub="RLE spikein", bSmooth=FALSE, dotSize=1.5)
par(mfrow=c(1,1))


###################################################
### code chunk number 81: loadParfac
###################################################
load(system.file('extra/parallelFactor.rda',package='DiffBind'))
parallelFactor


###################################################
### code chunk number 82: parfacMAnone
###################################################
dba.plotMA(parallelFactor,
           contrast=list(Fulvestrant=parallelFactor$masks$Fulvestrant),
           bNormalized=FALSE, sub="RAW", bSmooth=FALSE, dotSize=1.5)


###################################################
### code chunk number 83: parfacMA
###################################################
par(mfrow=c(2,1))
parallelFactor <- dba.normalize(parallelFactor,spikein = parallelFactor.peaks)
parallelFactor <- dba.analyze(parallelFactor)
dba.plotMA(parallelFactor, sub="LIB CTCF", bSmooth=FALSE, dotSize=1.5)

parallelFactor <- dba.normalize(parallelFactor,norm=DBA_NORM_RLE,
                                spikein = TRUE) 
parallelFactor <- dba.analyze(parallelFactor)
dba.plotMA(parallelFactor, sub="RLE CTCF", bSmooth=FALSE, dotSize=1.5)
par(mfrow=c(1,1))


###################################################
### code chunk number 84: DiffBind.Rnw:2206-2232
###################################################
require(xtable)
collabs <- paste(DBA_NORM_LIB," & ",DBA_NORM_RLE," & ",
                 DBA_NORM_TMM," & ",DBA_OFFSETS_LOESS)

rowlabs <- c(DBA_LIBSIZE_FULL,DBA_LIBSIZE_PEAKREADS,
             DBA_LIBSIZE_BACKGROUND,
             DBA_NORM_SPIKEIN,"parallel factor")

normtab <- matrix("X",5,4)
normtab[1,2:4] <- ""
normtab[c(1,3:5),4] <- ""
normtab <- data.frame(normtab)
rownames(normtab) <- rowlabs
addtorow <- list()
addtorow$pos <- list(0,0)
addtorow$command <- c("& \\multicolumn{4}{c|}{Normalization} \\\\\n",
                      paste("Reference & ",collabs,"\\\\\n"))

captionStr <- paste("Table of allowable normalization schemes.",
                    "Columns are normalization methods set by \\Rcode{normalize} or \\Rcode{offsets}.",
                    "Rows are reference reads, set by \\Rcode{library}, \\Rcode{background}, or \\Rcode{spikein}.")

print(xtable::xtable(normtab,align=c(rep("|c",5),"|"),
                     caption=captionStr, label="DiffBind-normtable"), 
      add.to.row = addtorow,
      include.colnames=FALSE, scalebox=1)


###################################################
### code chunk number 85: DiffBind.Rnw:2273-2274
###################################################
data(tamoxifen_peaks)


###################################################
### code chunk number 86: DiffBind.Rnw:2298-2300
###################################################
olap.rate <- dba.overlap(tamoxifen,mode=DBA_OLAP_RATE)
olap.rate


###################################################
### code chunk number 87: tamox_rate
###################################################
plot(olap.rate,type='b',ylab='# peaks',
     xlab='Overlap at least this many peaksets')


###################################################
### code chunk number 88: DiffBind.Rnw:2357-2358
###################################################
names(tamoxifen$masks)


###################################################
### code chunk number 89: DiffBind.Rnw:2367-2369
###################################################
dba.overlap(tamoxifen,tamoxifen$masks$MCF7 & tamoxifen$masks$Responsive,
            mode=DBA_OLAP_RATE)


###################################################
### code chunk number 90: tamox_mcf7_venn
###################################################
dba.plotVenn(tamoxifen, tamoxifen$masks$MCF7 & tamoxifen$masks$Responsive)


###################################################
### code chunk number 91: DiffBind.Rnw:2400-2403
###################################################
tamoxifen_consensus <- dba.peakset(tamoxifen, 
                                   consensus=c(DBA_TISSUE,DBA_CONDITION),
                                   minOverlap=0.66)


###################################################
### code chunk number 92: DiffBind.Rnw:2418-2422
###################################################
tamoxifen_consensus <- dba(tamoxifen_consensus,
                           mask=tamoxifen_consensus$masks$Consensus,
                           minOverlap=1)
tamoxifen_consensus


###################################################
### code chunk number 93: DiffBind.Rnw:2428-2429
###################################################
consensus_peaks <- dba.peakset(tamoxifen_consensus, bRetrieve=TRUE)


###################################################
### code chunk number 94: tamox_lines_venn
###################################################
data(tamoxifen_peaks)
tamoxifen <- dba.peakset(tamoxifen, consensus=DBA_TISSUE, minOverlap=0.66)
cons.ol <- dba.plotVenn(tamoxifen, tamoxifen$masks$Consensus)


###################################################
### code chunk number 95: DiffBind.Rnw:2469-2470
###################################################
data(tamoxifen_peaks)


###################################################
### code chunk number 96: DiffBind.Rnw:2476-2478
###################################################
dba.overlap(tamoxifen,tamoxifen$masks$Resistant,mode=DBA_OLAP_RATE)
dba.overlap(tamoxifen,tamoxifen$masks$Responsive,mode=DBA_OLAP_RATE)


###################################################
### code chunk number 97: tamox_cons_venn
###################################################
tamoxifen <- dba.peakset(tamoxifen, consensus=DBA_CONDITION, minOverlap=0.33)
dba.plotVenn(tamoxifen,tamoxifen$masks$Consensus)


###################################################
### code chunk number 98: DiffBind.Rnw:2513-2514
###################################################
tamoxifen.OL <- dba.overlap(tamoxifen, tamoxifen$masks$Consensus)


###################################################
### code chunk number 99: DiffBind.Rnw:2521-2523
###################################################
tamoxifen.OL$onlyA
tamoxifen.OL$onlyB


###################################################
### code chunk number 100: tamox_compare_venn
###################################################
tamoxifen <- dba.peakset(tamoxifen,tamoxifen$masks$Consensus,
                         minOverlap=1,sampID="OL Consensus")
tamoxifen <- dba.peakset(tamoxifen,!tamoxifen$masks$Consensus,
                         minOverlap=3,sampID="Consensus_3")
dba.plotVenn(tamoxifen,14:15)


###################################################
### code chunk number 101: DiffBind.Rnw:2559-2560
###################################################
data(tamoxifen_analysis)


###################################################
### code chunk number 102: DiffBind.Rnw:2570-2571
###################################################
tamoxifen.rep <- dba.report(tamoxifen,bCalled=TRUE,th=1)


###################################################
### code chunk number 103: DiffBind.Rnw:2580-2586
###################################################
onlyResistant <- tamoxifen.rep$Called1>=2 & tamoxifen.rep$Called2<3
sum(onlyResistant )
onlyResponsive <- tamoxifen.rep$Called2>=3 &  tamoxifen.rep$Called1<2
sum(onlyResponsive)
bothGroups <- tamoxifen.rep$Called1>= 2 & tamoxifen.rep$Called2>=3
sum(bothGroups)


###################################################
### code chunk number 104: DiffBind.Rnw:2601-2610
###################################################
tamoxifen.DB <- dba.report(tamoxifen,bCalled=TRUE)
onlyResistant.DB  <- (tamoxifen.DB$Called1 >= 2) & (tamoxifen.DB$Called2 < 3)
sum(onlyResistant.DB)
onlyResponsive.DB <- (tamoxifen.DB$Called2 >= 3) & (tamoxifen.DB$Called1 < 2)
sum(onlyResponsive.DB)
bothGroups.DB     <- (tamoxifen.DB$Called1 >= 2) & (tamoxifen.DB$Called2 >= 3)
sum(bothGroups.DB)
neitherGroup.DB   <- (tamoxifen.DB$Called1  < 2) & (tamoxifen.DB$Called2 < 3)
sum(neitherGroup.DB)


###################################################
### code chunk number 105: DiffBind.Rnw:2713-2714 (eval = FALSE)
###################################################
## DBA$config$design <- FALSE


###################################################
### code chunk number 106: sessionInfo
###################################################
toLatex(sessionInfo())