Question: CSAW : at which bam files are corresponding “up” regions
0
gravatar for JoannaF
3.7 years ago by
JoannaF0
JoannaF0 wrote:

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

ADD COMMENTlink modified 3.7 years ago by Aaron Lun25k • written 3.7 years ago by JoannaF0
Answer: CSAW : at which bam files are corresponding “up” regions
2
gravatar for Aaron Lun
3.7 years ago by
Aaron Lun25k
Cambridge, United Kingdom
Aaron Lun25k wrote:

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 COMMENTlink written 3.7 years ago by Aaron Lun25k

Thanks for your quick answer !

Best regards,

Joanna

ADD REPLYlink written 3.7 years ago by JoannaF0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 168 users visited in the last hour