Hi,

I have one question about gain and loss of enrichment with CSAW. Our bam files are:

Our bam files are:

bam.files <- c("sample1_rep1.bam", "sample1_rep2.bam", "sample2_rep1.bam")

I would like to know at which files are corresponding “up” regions detected by CSAW : it is sample1 or sample2 in this case ?

The whole R code I run is :

bam.files <- c("sample1_rep1.bam", "sample1_rep2.bam", "sample2_rep1.bam")
design <- model.matrix(~factor(c("sample1", "sample1", "sample2")))
colnames(design) <- c("intercept", "cell.type")
param <- readParam(minq=50, pe="both")
data <- windowCounts(bam.files, ext=150, width=400, param=param)
require(edgeR)
keep <- aveLogCPM(asDGEList(data)) >= -1
data2 <- data[keep,]
binned <- windowCounts(bam.files, bin=TRUE, width=10000, param=param)
normfacs <- normOffsets(binned)
y <- asDGEList(data2, norm.factors=normfacs)
y <- estimateDisp(y, design)
fit <- glmQLFit(y, design, robust=TRUE)
results <- glmQLFTest(fit)
merged <- mergeWindows(rowRanges(data2), tol=1000L)
tabcom <- combineTests(merged$id, results$table)
tab.best <- getBestTest(merged$id, results$table)
is.sig <- tabcom$FDR <= 0.05
up<-tab.best$logFC>0

require(rtracklayer)
test_up<-merged$region[is.sig & up]
test_up$score<-tab.best$logFC[is.sig & up]
export(test_up,"csaw_file_up.bed")
down<-tab.best$logFC<0
test_down<-merged$region[is.sig & down]
test_down$score<-tab.best$logFC[is.sig & down]
export(test_down,"csaw_file_down.bed")

Thanks a lot !

Joanna

Looking at your design matrix will reveal the answer:

  (Intercept) factor(c("sample1", "sample1", "sample2"))sample2
1           1                                                 0
2           1                                                 0
3           1                                                 1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$`factor(c("sample1", "sample1", "sample2"))`
[1] "contr.treatment"

By default, the GLM-based methods in edgeR will drop the last coefficient. In this case, the last coefficient represents the log-fold change of sample 2 over sample/group 1. So, positive log-fold changes (i.e., "up") will represent an increase in binding in sample 2.