A mismatch between DiffBind's results and the same results retreieved as DESeq2 (from DiffBind)
1
0
Entering edit mode
Sam ▴ 10
@sam-21502
Last seen 23 hours ago
Jerusalem

I retrieve DiffBind results as a DESeq2 object, and then look how many differentially bound regions are within it. When doing this, I get a different number of DB regions compared to the original DiffBind object.

library(DESeq2)
library(DiffBind)
data(tamoxifen_counts)
tamoxifen <- dba.contrast(tamoxifen, design="~Tissue")
tamoxifen <- dba.contrast(tamoxifen, group1=tamoxifen$masks$MCF7, group2=tamoxifen$masks$BT474, name1="MCF7", name2="BT474")
tamoxifen <- dba.analyze(tamoxifen)
dds <- dba.analyze(tamoxifen,bRetrieveAnalysis = DBA_DESEQ2)

res <- results(dds, independentFiltering=F, contrast=c("Tissue","MCF7","BT474"))

dba.show(tamoxifen,bContrasts = T,th=0.1)
summary(res)

sessionInfo( )

Output :

   Group Samples Group2 Samples2 DB.DESeq2
1  MCF7       5  BT474        2       954

out of 2845 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 738, 26%
LFC < 0 (down)     : 437, 15%
outliers [1]       : 0, 0%
low counts [2]     : 0, 0%
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

R version 4.0.4 (2021-02-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=en_IL.UTF-8       LC_NUMERIC=C               LC_TIME=en_IL.UTF-8       
 [4] LC_COLLATE=en_IL.UTF-8     LC_MONETARY=en_IL.UTF-8    LC_MESSAGES=en_IL.UTF-8   
 [7] LC_PAPER=en_IL.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_IL.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
 [1] parallel  stats4    grid      stats     graphics  grDevices utils     datasets  methods  
[10] base     

other attached packages:
 [1] DiffBind_3.0.13             rgl_0.105.13                limma_3.46.0               
 [4] DESeq2_1.30.1               SummarizedExperiment_1.20.0 Biobase_2.50.0             
 [7] MatrixGenerics_1.2.1        matrixStats_0.58.0          GenomicRanges_1.42.0       
[10] GenomeInfoDb_1.26.2         IRanges_2.24.1              S4Vectors_0.28.1           
[13] BiocGenerics_0.36.0         RColorBrewer_1.1-2          pheatmap_1.0.12            
[16] ggrepel_0.9.1               ggplot2_3.3.3              

loaded via a namespace (and not attached):
  [1] GOstats_2.56.0           backports_1.2.1          BiocFileCache_1.14.0    
  [4] plyr_1.8.6               GSEABase_1.52.1          splines_4.0.4           
  [7] BiocParallel_1.24.1      crosstalk_1.1.1          amap_0.8-18             
 [10] digest_0.6.27            invgamma_1.1             htmltools_0.5.1.1       
 [13] GO.db_3.12.1             SQUAREM_2021.1           fansi_0.4.2             
 [16] magrittr_2.0.1           checkmate_2.0.0          memoise_2.0.0           
 [19] BSgenome_1.58.0          base64url_1.4            Biostrings_2.58.0       
 [22] annotate_1.68.0          systemPipeR_1.24.3       bdsmatrix_1.3-4         
 [25] askpass_1.1              prettyunits_1.1.1        jpeg_0.1-8.1            
 [28] colorspace_2.0-0         apeglm_1.12.0            blob_1.2.1              
 [31] rappdirs_0.3.3           xfun_0.21                dplyr_1.0.4             
 [34] crayon_1.4.1             RCurl_1.98-1.2           jsonlite_1.7.2          
 [37] graph_1.68.0             genefilter_1.72.1        brew_1.0-6              
 [40] survival_3.2-7           VariantAnnotation_1.36.0 glue_1.4.2              
 [43] gtable_0.3.0             zlibbioc_1.36.0          XVector_0.30.0          
 [46] webshot_0.5.2            DelayedArray_0.16.1      V8_3.4.0                
 [49] Rgraphviz_2.34.0         scales_1.1.1             mvtnorm_1.1-1           
 [52] DBI_1.1.1                edgeR_3.32.1             miniUI_0.1.1.1          
 [55] Rcpp_1.0.6               emdbook_1.3.12           xtable_1.8-4            
 [58] progress_1.2.2           bit_4.0.4                rsvg_2.1                
 [61] truncnorm_1.0-8          AnnotationForge_1.32.0   htmlwidgets_1.5.3       
 [64] httr_1.4.2               gplots_3.1.1             ellipsis_0.3.1          
 [67] pkgconfig_2.0.3          XML_3.99-0.5             sass_0.3.1              
 [70] dbplyr_2.1.0             locfit_1.5-9.4           utf8_1.1.4              
 [73] tidyselect_1.1.0         rlang_0.4.10             manipulateWidget_0.10.1 
 [76] later_1.1.0.1            AnnotationDbi_1.52.0     munsell_0.5.0           
 [79] tools_4.0.4              cachem_1.0.4             cli_2.3.0               
 [82] generics_0.1.0           RSQLite_2.2.3            evaluate_0.14           
 [85] stringr_1.4.0            fastmap_1.1.0            yaml_2.2.1              
 [88] knitr_1.31               bit64_4.0.5              caTools_1.18.1          
 [91] purrr_0.3.4              RBGL_1.66.0              mime_0.10               
 [94] xml2_1.3.2               biomaRt_2.46.3           compiler_4.0.4          
 [97] rstudioapi_0.13          curl_4.3                 png_0.1-7               
[100] tibble_3.0.6             geneplotter_1.68.0       bslib_0.2.4             
[103] stringi_1.5.3            GenomicFeatures_1.42.1   lattice_0.20-41         
[106] Matrix_1.3-2             vctrs_0.3.6              pillar_1.5.0            
[109] lifecycle_1.0.0          jquerylib_0.1.3          irlba_2.3.3             
[112] data.table_1.14.0        bitops_1.0-6             httpuv_1.5.5            
[115] rtracklayer_1.50.0       R6_2.5.0                 latticeExtra_0.6-29     
[118] hwriter_1.3.2            promises_1.2.0.1         ShortRead_1.48.0        
[121] KernSmooth_2.23-18       MASS_7.3-53.1            gtools_3.8.2            
[124] assertthat_0.2.1         openssl_1.4.3            Category_2.56.0         
[127] rjson_0.2.20             withr_2.4.1              GenomicAlignments_1.26.0
[130] batchtools_0.9.15        Rsamtools_2.6.0          GenomeInfoDbData_1.2.4  
[133] hms_1.0.0                coda_0.19-4              DOT_0.1                 
[136] rmarkdown_2.7            GreyListChIP_1.22.0      ashr_2.2-47             
[139] mixsqp_0.3-43            bbmle_1.0.23.1           numDeriv_2016.8-1.1     
[142] shiny_1.6.0

954 != 738+437 ...

DiffBind • 1.1k views
ADD COMMENT
0
Entering edit mode
ADD REPLY
2
Entering edit mode
Rory Stark ★ 5.2k
@rory-stark-5741
Last seen 14 days ago
Cambridge, UK

The issue is how you set up the contrast. By using the group1 and group2 parameters, you are building a model with only the samples you specify in the groups. This is the way DiffBind used to work by default, prior to version 3.0. In order to use the new way, utilizing all of the samples and respecting the design parameter, you should specify the contrast the same way you would to DESeq2, as follows:

tamoxifen <- dba.contrast(tamoxifen, contrast=c("Tissue","MCF7","BT474"))

yielding the following output:

> tamoxifen <- dba.analyze(tamoxifen)
> dds <- dba.analyze(tamoxifen,bRetrieveAnalysis = DBA_DESEQ2)
> res <- results(dds, independentFiltering=F, contrast=c("Tissue","MCF7","BT474"))
> dba.show(tamoxifen,bContrasts = T,th=0.1)
  Factor Group Samples Group2 Samples2 DB.DESeq2
1 Tissue  MCF7       5  BT474        2      1175
> summary(res)

out of 2845 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 738, 26%
LFC < 0 (down)     : 437, 15%
outliers [1]       : 0, 0%
low counts [2]     : 0, 0%
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

738 + 437 = 1175.

ADD COMMENT

Login before adding your answer.

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