Hi, I am running Diffbind differential analysis of ATAC-seq samples. I am testing out the different options in diffbind and I could successfully run Differential analysis with DBA_DESEQ2 method with both library-size normalization and RLE-Rip normalization.
However, after normalizing with TMM and running differential analysis with method = DBA_EDGER, I get the error
DESeq2 analysis has not been run for this contrast
Could you please help me to resolve this issue? Thank you!
My samplesheet is
SampleID,Tissue,Condition,Treatment,Replicate,bamReads,Peaks,PeakCaller WT1_DMSO,HeLa,WT,DMSO,1,WT1_DMSO.bam,WT1_DMSO_peaks.narrowPeak,narrow WT2_DMSO,HeLa,WT,DMSO,2,WT2_DMSO.bam,WT2_DMSO_peaks.narrowPeak,narrow WT3_DMSO,HeLa,WT,DMSO,3,WT3_DMSO.bam,WT3_DMSO_peaks.narrowPeak,narrow WT1_TCDF,HeLa,WT,TCDF,1,WT1_TCDF.bam,WT1_TCDF_peaks.narrowPeak,narrow WT2_TCDF,HeLa,WT,TCDF,2,WT2_TCDF.bam,WT2_TCDF_peaks.narrowPeak,narrow WT3_TCDF,HeLa,WT,TCDF,3,WT3_TCDF.bam,WT3_TCDF_peaks.narrowPeak,narrow KO1_DMSO,HeLa,KO,DMSO,1,KO1_DMSO.bam,KO1_DMSO_peaks.narrowPeak,narrow KO2_DMSO,HeLa,KO,DMSO,2,KO2_DMSO.bam,KO2_DMSO_peaks.narrowPeak,narrow KO3_DMSO,HeLa,KO,DMSO,3,KO3_DMSO.bam,KO3_DMSO_peaks.narrowPeak,narrow KO1_TCDF,HeLa,KO,TCDF,1,KO1_TCDF.bam,KO1_TCDF_peaks.narrowPeak,narrow KO2_TCDF,HeLa,KO,TCDF,2,KO2_TCDF.bam,KO2_TCDF_peaks.narrowPeak,narrow KO3_TCDF,HeLa,KO,TCDF,3,KO3_TCDF.bam,KO3_TCDF_peaks.narrowPeak,narrow
My code is
# mamba install bioconda::bioconductor-diffbind
library(DiffBind)
SAMPLESHEET="samplesheet_PE_q0.05.csv"
# Minimum overlap for the samples to build the consensus peak set
# Default is 2
MIN_OVERLAP = 1
# Summit value; default = 200
SUMMITS=75
# Maximum sites to plot in the profile plot
MAX_SITES = 1000
# Method to perform differential analysis; deseq2 or edger
MOD = "edger"
#if (MOD=="deseq2") { method=DBA_DESEQ2 }
if (MOD=="edger") { method=DBA_EDGER }
# Counts Object; Takes time to compute in single core
# Efficient to run in parallel and save the counts object
# Counts object
COUNTS_OBJ="counts_PE_q0.05_dba_edger_object.rds"
# Read samplesheet
samples <- read.csv(SAMPLESHEET)
# Create DBA object using samplesheet
atac <- dba(sampleSheet=samples)
# Get read counts for consensus peaks
if (file.exists(COUNTS_OBJ)){
counts=readRDS(file = COUNTS_OBJ)
} else {
counts <- dba.count(atac, minOverlap=MIN_OVERLAP, summits=SUMMITS,bParallel=TRUE)
saveRDS(counts, file = COUNTS_OBJ)
}
# Normalize by library size
counts_norm <- dba.normalize(counts)
# Normalize by the DESeq2 or EdgeR
counts_norm <- dba.normalize(counts, method = method, normalize=DBA_NORM_NATIVE)
run_DBA <- function(obj, res_file, heatmap_file) {
print(method)
# Perform diff analysis
obj<- dba.analyze(obj,method=method)
# Differential expression summary
summary = dba.show(obj, bContrasts=TRUE)
# Get result files with all sites
res = dba.report(obj, th=1, bCounts=TRUE)
# Sorting the data frame by FDR in ascending order
res <- res[order(res$FDR),]
print(summary)
summary = as.list(summary)
# write Results to a file
write.csv(res, res_file, row.names = FALSE, quote=FALSE)
}
# KO vs WT
obj = dba.contrast(counts_norm, contrast=c("Condition", "KO", "WT"))
run_DBA(obj, "diff_KO_WT.csv", "diff_KO_WT_heatmpap.pdf")
> sessionInfo()
R version 4.3.2 (2023-10-31)
Platform: x86_64-apple-darwin20 (64-bit)
Running under: macOS Monterey 12.7.1
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/New_York
tzcode source: internal
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] DiffBind_3.12.0 SummarizedExperiment_1.32.0 Biobase_2.62.0
[4] MatrixGenerics_1.14.0 matrixStats_1.1.0 GenomicRanges_1.54.1
[7] GenomeInfoDb_1.38.0 IRanges_2.36.0 S4Vectors_0.40.1
[10] BiocGenerics_0.48.1
loaded via a namespace (and not attached):
[1] amap_0.8-19 tidyselect_1.2.0 dplyr_1.1.3 Biostrings_2.70.1
[5] bitops_1.0-7 fastmap_1.1.1 RCurl_1.98-1.13 apeglm_1.24.0
[9] GenomicAlignments_1.38.0 XML_3.99-0.15 digest_0.6.33 lifecycle_1.0.4
[13] statmod_1.5.0 invgamma_1.1 magrittr_2.0.3 compiler_4.3.2
[17] rlang_1.1.2 tools_4.3.2 utf8_1.2.4 yaml_2.3.7
[21] rtracklayer_1.62.0 S4Arrays_1.2.0 htmlwidgets_1.6.2 interp_1.1-5
[25] DelayedArray_0.28.0 plyr_1.8.9 RColorBrewer_1.1-3 ShortRead_1.60.0
[29] abind_1.4-5 BiocParallel_1.36.0 KernSmooth_2.23-22 numDeriv_2016.8-1.1
[33] hwriter_1.3.2.1 grid_4.3.2 fansi_1.0.5 latticeExtra_0.6-30
[37] caTools_1.18.2 colorspace_2.1-0 edgeR_4.0.2 ggplot2_3.4.4
[41] scales_1.2.1 gtools_3.9.5 MASS_7.3-60 bbmle_1.0.25.1
[45] cli_3.6.1 mvtnorm_1.2-4 crayon_1.5.2 generics_0.1.3
[49] rstudioapi_0.15.0 rjson_0.2.21 bdsmatrix_1.3-6 stringr_1.5.0
[53] splines_4.3.2 zlibbioc_1.48.0 parallel_4.3.2 restfulr_0.0.15
[57] XVector_0.42.0 vctrs_0.6.4 Matrix_1.6-1.1 mixsqp_0.3-54
[61] ggrepel_0.9.4 irlba_2.3.5.1 systemPipeR_2.8.0 jpeg_0.1-10
[65] locfit_1.5-9.8 limma_3.58.1 glue_1.6.2 emdbook_1.3.13
[69] codetools_0.2-19 stringi_1.7.12 gtable_0.3.4 deldir_2.0-2
[73] BiocIO_1.12.0 munsell_0.5.0 tibble_3.2.1 pillar_1.9.0
[77] htmltools_0.5.7 gplots_3.1.3 BSgenome_1.70.1 GenomeInfoDbData_1.2.11
[81] truncnorm_1.0-9 R6_2.5.1 GreyListChIP_1.34.0 lattice_0.22-5
[85] png_0.1-8 Rsamtools_2.18.0 SQUAREM_2021.1 ashr_2.2-63
[89] Rcpp_1.0.11 coda_0.19-4 SparseArray_1.2.2 DESeq2_1.42.0
[93] pkgconfig_2.0.3
