Errror in running EdgeR based differential analysis in diffbind
0
0
Entering edit mode
Aswathy • 0
@c8518096
Last seen 9 months ago
United States

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
DiffBind • 309 views
ADD COMMENT

Login before adding your answer.

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