CSAW : at which bam files are corresponding “up” regions
1
0
Entering edit mode
JoannaF ▴ 10
@joannaf-9881
Last seen 12 weeks ago
France

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

csaw up regions bam file gain and loss of enrichment • 918 views
ADD COMMENT
2
Entering edit mode
Aaron Lun ★ 27k
@alun
Last seen 6 hours ago
The city by the bay

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.

ADD COMMENT
0
Entering edit mode

Thanks for your quick answer !

Best regards,

Joanna

ADD REPLY

Login before adding your answer.

Traffic: 218 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6