DiffBind fold parameter vs. IfcThreshold in DESeq2
Entering edit mode
Shuwan ▴ 10
Last seen 5 weeks ago
United States


I am currently using DiffBind to analyze ATAC-seq data. I noticed that when I called DESeq2::results to analyze the DESeq2 object retrieved from dba.analyze using the parameters alpha = 0.05, lfcThreshold = log2(1.3) , it gave me a different number of DB sites from dba.report using th = 0.05, fold = log2(1.3). I figured that DiffBind did ad hoc filtering in which the fold parameter was only used to filter peaks after the DESeq2::results using alpha = 0.05 is done. But what I would like to do is to include the fold change parameter in the statistical testing. The output of DESeq2::results of the retrieved DESeq2 object though lack the peak information (chr, start, end), so I can not directly use it for any further analysis. The rownames of the output dataframe are simply 1, 2, 3, 4, 5...

I was wondering if there is a way to recover the peak info so I can directly use the results of DESeq2::results? Or how can I include the fold value in the statistical testing when calling dba.report? Thank you so much!!!

sessionInfo( )
R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

other attached packages:
 [1] edgeR_3.34.0                              limma_3.48.0                             
 [3] DESeq2_1.32.0                             org.Mm.eg.db_3.13.0                      
 [5] TxDb.Mmusculus.UCSC.mm10.knownGene_3.10.0 GenomicFeatures_1.44.0                   
 [7] AnnotationDbi_1.54.1                      ChIPseeker_1.28.3                        
 [9] data.table_1.14.0                         DiffBind_3.2.2                           
[11] SummarizedExperiment_1.22.0               Biobase_2.52.0                           
[13] MatrixGenerics_1.4.0                      matrixStats_0.59.0                       
[15] GenomicRanges_1.44.0                      GenomeInfoDb_1.28.0                      
[17] IRanges_2.26.0                            S4Vectors_0.30.0                         
[19] BiocGenerics_0.38.0                      

loaded via a namespace (and not attached):
  [1] utf8_1.2.1                              tidyselect_1.1.1                       
  [3] RSQLite_2.2.7                           grid_4.1.0                             
  [5] BiocParallel_1.26.0                     scatterpie_0.1.6                       
  [7] munsell_0.5.0                           base64url_1.4                          
  [9] systemPipeR_1.26.2                      withr_2.4.2                            
 [11] colorspace_2.0-1                        GOSemSim_2.18.0                        
 [13] Category_2.58.0                         filelock_1.0.2                         
 [15] knitr_1.33                              rstudioapi_0.13                        
 [17] DOSE_3.18.0                             bbmle_1.0.23.1                         
 [19] GenomeInfoDbData_1.2.6                  mixsqp_0.3-43                          
 [21] hwriter_1.3.2                           polyclip_1.10-0                        
 [23] bit64_4.0.5                             farver_2.1.0                           
 [25] pheatmap_1.0.12                         treeio_1.16.1                          
 [27] coda_0.19-4                             batchtools_0.9.15                      
 [29] vctrs_0.3.8                             generics_0.1.0                         
 [31] xfun_0.23                               csaw_1.26.0                            
 [33] BiocFileCache_2.0.0                     R6_2.5.0                               
 [35] apeglm_1.14.0                           graphlayouts_0.7.1                     
 [37] invgamma_1.1                            locfit_1.5-9.4                         
 [39] rsvg_2.1.2                              bitops_1.0-7                           
 [41] cachem_1.0.5                            fgsea_1.18.0                           
 [43] DelayedArray_0.18.0                     assertthat_0.2.1                       
 [45] BiocIO_1.2.0                            scales_1.1.1                           
 [47] ggraph_2.0.5                            enrichplot_1.12.0                      
 [49] gtable_0.3.0                            tidygraph_1.2.0                        
 [51] rlang_0.4.11                            genefilter_1.74.0                      
 [53] splines_4.1.0                           lazyeval_0.2.2                         
 [55] rtracklayer_1.52.0                      brew_1.0-6                             
 [57] checkmate_2.0.0                         BiocManager_1.30.15                    
 [59] yaml_2.2.1                              reshape2_1.4.4                         
 [61] backports_1.2.1                         qvalue_2.24.0                          
 [63] RBGL_1.68.0                             tools_4.1.0                            
 [65] ggplot2_3.3.3                           ellipsis_0.3.2                         
 [67] gplots_3.1.1                            RColorBrewer_1.1-2                     
 [69] Rcpp_1.0.6                              plyr_1.8.6                             
 [71] progress_1.2.2                          zlibbioc_1.38.0                        
 [73] purrr_0.3.4                             RCurl_1.98-1.3                         
 [75] prettyunits_1.1.1                       viridis_0.6.1                          
 [77] cowplot_1.1.1                           ashr_2.2-47                            
 [79] ggrepel_0.9.1                           tinytex_0.32                           
 [81] magrittr_2.0.1                          DO.db_2.9                              
 [83] truncnorm_1.0-8                         mvtnorm_1.1-2                          
 [85] SQUAREM_2021.1                          amap_0.8-18                            
 [87] patchwork_1.1.1                         hms_1.1.0                              
 [89] evaluate_0.14                           xtable_1.8-4                           
 [91] XML_3.99-0.6                            emdbook_1.3.12                         
 [93] jpeg_0.1-8.1                            gridExtra_2.3                          
 [95] compiler_4.1.0                          biomaRt_2.48.0                         
 [97] bdsmatrix_1.3-4                         tibble_3.1.2                           
 [99] shadowtext_0.0.8                        KernSmooth_2.23-20                     
[101] V8_3.4.2                                crayon_1.4.1                           
[103] htmltools_0.5.1.1                       GOstats_2.58.0                         
[105] aplot_0.0.6                             tidyr_1.1.3                            
[107] geneplotter_1.70.0                      DBI_1.1.1                              
[109] tweenr_1.0.2                            dbplyr_2.1.1                           
[111] MASS_7.3-54                             rappdirs_0.3.3                         
[113] boot_1.3-28                             ShortRead_1.50.0                       
[115] Matrix_1.3-4                            cli_2.5.0                              
[117] metapod_1.0.0                           igraph_1.2.6                           
[119] pkgconfig_2.0.3                         TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
[121] rvcheck_0.1.8                           GenomicAlignments_1.28.0               
[123] numDeriv_2016.8-1.1                     ggtree_3.0.2                           
[125] annotate_1.70.0                         XVector_0.32.0                         
[127] AnnotationForge_1.34.0                  stringr_1.4.0                          
[129] VariantAnnotation_1.38.0                digest_0.6.27                          
[131] graph_1.70.0                            Biostrings_2.60.1                      
[133] rmarkdown_2.8                           fastmatch_1.1-0                        
[135] tidytree_0.3.4                          GSEABase_1.54.0                        
[137] restfulr_0.0.13                         GreyListChIP_1.24.0                    
[139] curl_4.3.1                              Rsamtools_2.8.0                        
[141] gtools_3.9.2                            rjson_0.2.20                           
[143] nlme_3.1-152                            lifecycle_1.0.0                        
[145] jsonlite_1.7.2                          viridisLite_0.4.0                      
[147] BSgenome_1.60.0                         fansi_0.5.0                            
[149] pillar_1.6.1                            lattice_0.20-44                        
[151] plotrix_3.8-1                           KEGGREST_1.32.0                        
[153] fastmap_1.1.0                           httr_1.4.2                             
[155] survival_3.2-11                         GO.db_3.13.0                           
[157] glue_1.4.2                              png_0.1-7                              
[159] bit_4.0.4                               Rgraphviz_2.36.0                       
[161] ggforce_0.3.3                           stringi_1.6.2                          
[163] blob_1.2.1                              latticeExtra_0.6-29                    
[165] caTools_1.18.2                          memoise_2.0.0                          
[167] DOT_0.1                                 dplyr_1.0.6                            
[169] ape_5.5                                 irlba_2.3.3
DiffBind • 137 views
Entering edit mode
Rory Stark ★ 4.1k
Last seen 19 days ago
CRUK, Cambridge, UK

The good news is that this feature has already been built into the development version of DiffBind, which will be generally available in the next Bioconductor release (usually October). If, when you call dba.report(), you specify a minimum absolute fold change other than 0, the tests will include the lfcThreshold (this also works for an edgeR analysis via edgeR::glmTreat()).

That probably doesn't help you much now!

In the meantime you can map the peak "numbers" you get back from DESeq2::results() as you suggest. Retrieve the full consensus set by calling dba.peakset() with bRetrieve=TRUE. Then the rownames returned by DESeq2::results() can be used as indices into the consensus peakset (returned by default as a GRanges object).

Entering edit mode

Thank you so much for the quick response, Rory!!! Happy to hear that this feature will be included in the next release. I just mapped the consensus peak back to the DESeq2 results and it worked perfectly! Thank you!


Login before adding your answer.

Traffic: 190 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6