Question

DiffBind: report only upregulated differential sites

0

Entering edit mode

meisan406 • 0

@meisan406-7061

Last seen 7.0 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 • 2.5k views

ADD COMMENT • link updated 10.0 years ago by Rory Stark ★ 5.2k • written 10.0 years ago by meisan406 • 0

score 0 · Answer 1 · 2014-11-27

0

Entering edit mode

Rory Stark ★ 5.2k

@rory-stark-5741

Last seen 4 weeks ago

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 • link 10.0 years ago Rory Stark ★ 5.2k

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 • link 10.0 years ago meisan406 • 0

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 • link 10.0 years ago Rory Stark ★ 5.2k