I am using the Diffbind package to do the differential binding analysis for my ATAC-seq library. I followed the instruction to use blocking methods for my library, because I assumed there were substantial batch effects. Since I didn't know which one of DESeq2_block and Edger_block would generate the best result for me, I tried my library with both of the normalization methods. I characterized log(fold change) > 1 and p-value <0.05 as "positively enriched peaks", and log(fold change) < -1 and p-value<0.05 as "negatively enriched peaks". However, DESeq2_block generated 1530 positively enriched peaks and 556 negatively enriched peaks, while EdgeR_block generated 606 positively enriched peaks and 1811 negatively enriched peaks. The 606 positively enriched peaks from EdgeR_block is almost a subset of 1530 positively enriched peaks generated by DESeq2, and no negatively enriched peaks were characterized as positively enriched peaks in the other methods. Therefore I have the following questions:
1) What is the best way to make a judgement call between these two blocking methods?
2) I was advised to look at the normalized data through the "ReportingTools" library, but this library requires a DESeqDataSet or DESeqResults as input. How do I extract the DESeqDataSet from DBA object properly? I tried to follow the manual by typing: atac.analyze$contrasts[]$DESeq2$block$DEdata, but it does not seem to have $DEdata.
3) If I also want to extract the EdgeR object from DBA what should I do?
Thank you so much!
atac.list <- dba(sampleSheet=filename)
atac.count <- dba.count(atac.list)
atac.contrast <- dba.contrast(atac.count, categories=DBA_CONDITION, block=DBA_REPLICATE, minMembers = 2)
atac.analyze <- dba.analyze(atac.contrast, method=DBA_ALL_METHODS)
###try to see the normalized signal values for a peak
des2Report <- HTMLReport(shortName = "Differential_analysis_screening_library",
title = 'differential analysis of screening library using DESeq2', reportDirectory ="./DESeq2/reports")
factor = colData(dds)$facs, reportDir="./DESeq2/reports")