I using DiffBind to compare WT vs Treatment ATAC-seq samples and I have two questions/problems regarding the DiffBind
dba.plotMA() indicates more differential sites than
dba.report() with the same thresholds.
The input csv:
SampleID,Tissue,Factor,Condition,Replicate,bamReads,bamControl,Peaks WT_1,hs,ATAC,WT,1,WT_1.bam,,WT_1.bed WT_2,hs,ATAC,WT,2,WT_2.bam,,WT_2.bed Treat_1,hs,ATAC,Treat,1,Treat_1.bam,,Treat_1.bed Treat_2,hs,ATAC,Treat,2,Treat_2.bam,,Treat_2.bed
## Load DiffBind sample design. ATAC.dba = dba(sampleSheet = "csv/samples.csv", peakFormat = "bed") ## Count reads in binding site intervals. 4 Samples, 53778 sites in matrix: ATAC.dba.count = dba.count(ATAC.dba, bLog = T, score = DBA_SCORE_TMM_READS_FULL, filter = 10) ## Set up contrasts for differential binding affinity analysis. ATAC.dba.count = dba.contrast(ATAC.dba.count, categories = DBA_CONDITION, minMembers = 2) ## Perform differential binding affinity analysis. Contrast: ## Group1 Members1 Group2 Members2 DB.DESeq2 ## WT 2 Treat 2 2805 ATAC.dba.count = dba.analyze(ATAC.dba.count, bFullLibrarySize = TRUE) ## The MA-plot title: ' 1277 FDR < 0.050' dba.plotMA(ATAC.dba.count, th = 0.05, bUsePval = FALSE, fold = 2, bFlip = TRUE) ## dba.report gives 1274 differential sites however, with the same parameters. ATAC.dba.count.DB = data.frame(dba.report(ATAC.dba.count, th = 0.05, bUsePval = FALSE, fold = 2, bFlip = TRUE)) ## With the approach below, I can reproduce the MA-plot DB count. ATAC.dba.count.MA = data.frame(dba.report(ATAC.dba.count, th = 1, bCounts = FALSE, bFlip = TRUE)) ATAC.dba.count.MA = mutate(ATAC.dba.count.MA, sig = ifelse(abs(Fold) >= 2 & FDR <= 0.05, "sig", "ns")) table(ATAC.dba.count.MA$sig) ## ns (52501) sig (1277) ## With the code below I could find the three missing DB site, and they are true DB sites after manual inspections. tmpdf = merge(ATAC.dba.count.MA, ATAC.dba.count.DB, by = c("seqnames", "start", "end"), all = TRUE) subset(tmpdf, sig == "sig" & is.na(width.y) == T)
dba.report() seems to report not normalised values for a sample. This is also present in the tamoxifen dataset. When I try to report the DB sites MCF71 sample have the same values even is bNormalized is turned on.
Am I missing something?
data("tamoxifen_analysis") dba.report(tamoxifen, th = 0.05, bUsePval = FALSE, fold = 2, bCounts = TRUE, bNormalized = FALSE) dba.report(tamoxifen, th = 0.05, bUsePval = FALSE, fold = 2, bCounts = TRUE, bNormalized = TRUE)
Thank you for the help!