csaw_example.R

grouping <- c("A","A","B","B")
design<-model.matrix(~factor(grouping))


#############################################################
### csaw, with its combined window methodology.
############################################################

xparam <- readParam(dedup=FALSE)

if (is.tf) { 
  test.widths <- 10
  names <- c("csaw")
} else {
  test.widths <- c(50, 150, 250)
  names <- c("csaw.short", "csaw", "csaw.long")
}

for (w in seq_along(test.widths)) { 
  data <- windowCounts(bam.files, width=test.widths[w], ext=fraglen, param=xparam, filter=20)
  binned <- windowCounts(bam.files, width=2000, bin=TRUE, param=xparam)

  bin.ab <- scaledAverage(asDGEList(binned), scale=median(getWidths(binned))/median(getWidths(data)))
  threshold <- median(bin.ab) + log2(2)
  keep <- aveLogCPM(asDGEList(data)) >  threshold

  data <- data[keep,]
  tabres <- analyzeQLCounts(assay(data), design, totals=data$totals)
  merged <- mergeWindows(rowRanges(data), tol=100, max.width=5000)
  tabneg <- combineTests(merged$id, tabres)

  for (cutoff in fdr.thresholds) { 
    resultDump(merged$region, tabneg, cutoff, out=ofile)
    result<-assessChIP(ofile, lfile, tol=NA, checkfc=FALSE) 
    dump2file(names[w], cutoff, result)
  }
}