DiffBind: report only upregulated differential sites
1
0
Entering edit mode
meisan406 • 0
@meisan406-7061
Last seen 3.9 years ago
United States

This is a question on the dba.report function of DiffBind.

I'm trying to retrieve only the upregulated differential sites and have been using the "bUP" and "fold" arguments with the dba.report function, but I've not been able to filter out to obtain only the upregulated differential sites.

My R scripts for this purpose are as follow:

contrast.res <- dba.contrast(count.data, dba.object$masks$IL4c, dba.object$masks$Res, "IL4c", "Res")

diff.analysis.res <- dba.analyze(contrast.res, method=DBA_DESEQ2, bCorPlot=TRUE, bReduceObjects=TRUE)

dba.report(diff.analysis.res, method=DBA_DESEQ2, th=0.1, bUsePval=FALSE, fold=0, file="diff_sites_res_up300")

dba.report(diff.analysis.res, method=DBA_DESEQ2, th=0.1, bUsePval=FALSE, bALL=FASLSE, bUP=TRUE, file="diff_sites_res_up300") 

I'm using R 3.0.2 and DiffBind 1.8.5

I'd appreciate any help on this.

Thanks,

Mei San

 

diffbind • 1.8k views
ADD COMMENT
0
Entering edit mode
Rory Stark ★ 4.1k
@rory-stark-5741
Last seen 23 days ago
CRUK, Cambridge, UK

Hi Mei San-

The parameters you are using don't work quite that way. The fold parameter always refers to absolute values of fold changes (log2), where only results with absolute fold changes greater than or equal to this are included, so it can't be used to distinguish between positive and negative fold changes. The bAll and bUp/bDown   parameters (renamed bGain/bLoss in more recent versions) only have meaning when at least one of the bDB/bNotDB parameters are set to TRUE, indicating that a report-based DBA object should be returned. This is useful if you want the results as peaksets that you can, for example, overlap and use to plot Venn diagrams.

If you want separate reports for the positive and negative fold changes, the most straightforward way is to get a report with both, then subset it:

> db.all  = dba.report(diff.analysis.res, method=DBA_DESEQ2, th=0.1)
> db.gain = db.all[db.all$Fold>0,]
> db.loss = db.all[db.all$Fold<0,]

Hope this helps!

Cheers-

Rory

ADD COMMENT
0
Entering edit mode

Hi Rory,

Thanks, it helped. 

However, I've a couple points for clarification:

1. Could you suggest a way to capture the subset negative and positive fold changes reports as a csv file (equivalent to the "file" parameter in dba.report)? I eventually exported the GRanges object as a bed file using rtracklayer's export.bed function, but I would still like to capture a direct csv file from the subset reports for comparison. The conventional write.table won't work.

2. When I exported the differential sites report as a bed file, I did not have the strand information (i.e. if it's + or -). This is because the count data was used to generate the output and there was no strand information involved, that's why it's lost in the eventual bed file - please correct me if I'm wrong.

Thanks,

Mei San

ADD REPLY
0
Entering edit mode

Hi Mei San-

1, It sounds like you have found a way to do this already. You can;t really do it directly from the interface. I can tell you that the write function used internally is:

> write.csv(db.gain, file="gain.csv", row.names=FALSE, quote=FALSE)

2. You are more or less correct about strand information, Basically peak data from ChIP is not strand-specific -- we can detect that the protein is binding to the DNA at the given location, but not if the binding is strand-specific. Peak intervals are generally specified on the + strand by default.

Cheers-

Rory

ADD REPLY

Login before adding your answer.

Traffic: 299 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