Hi everyone, I have been analyzing some ATAC Seq data to find differential peaks and the enriched motifs within them. I have 2 replicates for each condition; the comparisons are of the type noTx vs Tx, and I have 5 comparisons in total (20 samples total, 10 conditions, each Tx has a corresponding noTx). I had some preliminary results, where I took on macs2 the parameter -t as the noTx samples and -c as Tx (Treatment has been applied, and I wanted this as reference/control). bams have been processed to remove: mitochondrial genome, duplicates and black listed regions. I then ran HOMER findmotifGenome and found enriched motifs. I realized afterwards that these peaks were actually not differentially enriched peaks, maybe I could categorize them as enriched regions after removing the background, where the background would be Tx.
So, I called peaks for all of the samples independently and defined a set of consensus peaks merging the .narrowPeak files of all of the 20 samples (about 70,000 peaks total). I then annotated the peaks with HOMER to try to find useful labeling, and I used the function multiBamSummary of deepTools with the BED-file option, and with the --outRawCounts the program gave me the count matrix. The thing is that when I took the matrix as input for DESeq2, the results were bad, and no peak had a good padj value, all of them were 1. Literally the volcanos are just a horizontal line.
Do you think voom and limma can make a difference? Maybe I should merge peaks only based on the comparison I will make, so instead of a general peak set, define one for every one of the 5 comparisons, only merging the .narrowPeaks of the invloved samples and not the total list? Is the macs2 approach to just consider the samples you would set as control in DESeq2, as the -c parameter is just wrong? What is the difference between that and diffPeaks?
Any advice is highly appreciated.
Thanks a lot in advance, Bruno
