Conc_1 in differnet dba.report result is slightly different
0
Entering edit mode
@shangguandong1996-21805
Last seen 5 hours ago

Hi, Dr Stark

I found when I pool three group in a experiment, like groupA, B, C. The Conc_A in dba.report is slightly different between GroupB_VS_GroupA and GroupC_VS_Group_A.

Just my own data

library(tidyverse)
a <- read.csv("result/Diff_Peak_Anno/H3K27ac/SIM3d_H3K27ac_VS_CIM7d_H3K27ac.csv") %>% 
  select(feature_id, Conc_CIM7d_H3K27ac)
b <- read.csv("result/Diff_Peak_Anno/H3K27ac/SIM8d_H3K27ac_VS_CIM7d_H3K27ac.csv") %>% 
  select(feature_id, Conc_CIM7d_H3K27ac)

# sometimes, the diff will be bigger
> inner_join(a, b, by = "feature_id") %>% 
+   mutate(diff = Conc_CIM7d_H3K27ac.x - Conc_CIM7d_H3K27ac.y) %>% 
+   pull(diff) %>% 
+   table()
.
                  0 0.00999999999999979                0.01 
               3971               13053                  76 
 0.0100000000000002  0.0100000000000007  0.0100000000000016 
                 75                1083                1297

by the way, I find it is hard to open the link posted in the DiffBind manual.

http://https//www.cruk.cam.ac.uk/core-facilities/bioinformatics-core/software/diffbind

Best wishes

Guandong Shang

DiffBind • 53 views
ADD COMMENTlink
0
Entering edit mode
Rory Stark ♦ 3.5k
@rory-stark-5741
Last seen 3 days ago
CRUK, Cambridge, UK

First, thanks for the pointing out the bad link, I've checked in a fix.

On your main issue, I'm not sure how you derived your .csv files. Could you show how you went from a report to the .csv? Did you use dba.report() with file= or some other method?

Here's a check I did to check for this issue:

data(tamoxifen_counts)
tamoxifen <- dba.contrast(tamoxifen,
                          contrast=c("Tissue","BT474","MCF7"))
tamoxifen <- dba.contrast(tamoxifen,
                          contrast=c("Tissue","BT474","T47D"))
tamoxifen <- dba.analyze(tamoxifen)

r1 <- sort(dba.report(tamoxifen, contrast=1, th=1, precision=0))
r2 <- sort(dba.report(tamoxifen, contrast=2, th=1, precision=0))

sum(r1$Conc_BT474 != r2$Conc_BT474)
ADD COMMENTlink
0
Entering edit mode

I use as_tibble to convert Grange objecet, and then use write_csv. The basic command is like below


# I have to admit this example is not reproducibility
# Once the data link is fixed, I will post the test data result here again.

dba_contrast <- dba.contrast(XX)
dba_diff <- dba.analyze(dba_contrast)
dba_report_all <- dba.report(dba_diff,th = 1)

# annotatePeak is the function of ChIPseeker package
peakAnno <- annotatePeak(dba_report_all,
                           TxDb = TxDb,
                           level = "gene",
                           tssRegion = c(-500, 500)
                           )


peakAnno_geneSymbol <- left_join(as_tibble(peakAnno@anno),
                                   gene_alias,
                                   by = c("geneId" = "name"))

peakAnno_geneSymbol %>% 
    readr::write_csv("XX.csv")

By the way, I found the N in Binding Affinity: treat vs. control (N FDR < 0.050)ยท of MA-plot is different from the N result by using filter in R or excel. I hope this phenomenon may be helpful.

For example, the N in MA-plot text may be 6809 FDR < 0.05 while in my result is 6801.

# my MA-plot function code
# Binding Affinity: YE vs. Col_0 (6809 FDR < 0.050)
dba.plotMA(dba_diff,fold = 2,cex.main=0.8)

# neither >= nor > is same as text in MA-plot
read_csv("result/Diff_Peak_Anno/THB_ChIP/YE_VS_Col_0.csv") %>% 
  filter(abs(Fold) >= 2, FDR < 0.05) %>% 
  nrow()

[1] 6801


read_csv("result/Diff_Peak_Anno/THB_ChIP/YE_VS_Col_0.csv") %>% 
  filter(abs(Fold) > 2, FDR < 0.05) %>% 
  nrow()

[1] 6796
ADD REPLYlink

Login before adding your answer.

Similar Posts
Loading Similar Posts
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.3