Hello,

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

locale:
[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

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).