DiffBind - question about showing significant difference on peaks with 0 fold change
1
0
Entering edit mode
Ti • 0
@8b5051c7
Last seen 8 months ago
Japan

Hi, I've used DiffBind to analyze differential binding analysis for cut&tag data.

The result is showed in MAplot:

enter image description here

I noticed some peaks with 0 log fold change are differentially biding sites, while it could be a little bit confusing.

I see this result should be caused by the overall gaining of binding affinity, as the red curve indicates.

However, this data is spike-in normalized, and the overall gaining of binding affinity might be reasonable in our case.

So how can I get differential binding sites without considering the overall gaining of signal?

I can manually set a cutoff (e.g. abs(log FC) > 1) to determine, however the p-value and FDR should be re-caculated, which I don't know how to do in DiffBind.

Here is my code, just regular steps but with a spike-in normalization:

sample = read.csv("samplesheet.csv")
CTCF <- dba(sampleSheet=sample)
CTCF <- dba.count(CTCF)
scale_factors = read.delim("scale_factors.txt", header=T, row.names=1)
CTCF <- dba.normalize(CTCF, normalize=scale_factors$scale_factor)
MUT2_vs_WT2 <- dba.contrast(CTCF, design=FALSE, group1 = CTCF$masks$M2, group2 = CTCF$masks$W2, name1 = "MUT2", name2 = "WT2",  minMembers = 2)
MUT2_vs_WT2 <- dba.analyze(MUT2_vs_WT2)

Thanks for your kindly help in advance!

DiffBind • 1.3k views
ADD COMMENT
0
Entering edit mode

Rory Stark Could you please help take a look? Really appreciate it!

ADD REPLY
0
Entering edit mode
Rory Stark ★ 5.1k
@rory-stark-5741
Last seen 8 weeks ago
Cambridge, UK

I'm wondering about how you calculated the scaling factors, they really need to be correct!

It would be good to compare the MA plot using bNormalized=FALSE.

Also, is there a reason you're using design=FALSE? I would think you'd want to set the contrast, ie

MUT2_vs_WT2 <- dba.contrast(CTCF, contrast=c("Condition", "M2","W2"))

Have a look at the data returned by dba.report(), especially the Fold values. Do these values look right gievn the Conc_ values for the two groups?

ADD COMMENT
0
Entering edit mode

Thanks for your kindly reply Mr. Stark!

As our samples are spiked-in with E.coli DNA, the scaling factors should be proportional to the ratio of E.coli DNA, e.g a sample with 0.61% reads mapped to E.coli DNA has scale factor of 0.61. I've read about your previous answer (CUT&Tag differential peak analysis using Diffbind, with pre-calculated normalization factor.) and triple-checked my formula so it should be correct.

After using bNormalized=FALSE, the MAplot looks like:

enter image description here

Does that mean I shouldn't use normalization factors computed by differential analysis method? However I didn't see argument like bNormalized in dba.analyze().

There is no specific reason I use design=FALSE, just because I'm not very familiar with DiffBind.

I have checked about result of dba.report(), and the FC values are right given the conc values. We also have other two groups and they don't show the overall shift of signals. The other MA plots are like below and they seems more reasonable:

enter image description here

Anyway, do you have any suggestions about how to eliminate the overall signal shifting and estimate the true differential binding sites? Please forgive my ignorance, I don't really understand the details behind the differential analysis methods. Thanks again!

ADD REPLY
0
Entering edit mode

Comparing the normalized and non-normalize MA plots, it looks like the normalization is shifting everything up. I'm not sure if this is an issue with how the MA plot is being calculated or if it is really how it is being modeled by the GLM.

Could you make your CTCF DBA object available to me, along with the scale_factors? I'd like to get to the bottom of this.

ADD REPLY
0
Entering edit mode

Sure, I have saved it using

save(MUT2_vs_WT2, file='M2_vs_W2.Rdata')

And the scale factors used can be found in

MUT2_vs_WT2$norm$DESeq2$norm.facs

I have uploaded the Rdata to:

https://github.com/Devalock/tmp

Thanks for your kind help!

ADD REPLY
0
Entering edit mode

I've been away on leave but am back looking at this if you can send an updated link.

ADD REPLY
0
Entering edit mode

Sure, I've updated the link above. Thank you!

Could you please take a quick look when you are available Mr.Strak and really appreciate it! Rory Stark

ADD REPLY

Login before adding your answer.

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