Hello, 

I am new to ATAC-analysis, and wish to ask for some advice in identifying significant differences between two ATAC-seq samples. 

Firstly, I have two conditions - control and knockdown. 
I have mapped the PE reads to a reference genome, and then subsequently taken reads lower than 100bp (i.e. open chromatin regions). 
I then use MACS2 to turn my two un-normalised BAM files to peaks using this command: 

    macs2 callpeak -t macs2 callpeak -t bamfile --outdir /path/to/ -f BAMPE --keep-dup all --pvalue 1e-2 --call-summits --bdg

I then use feature counts to quantify mapped sequencing reads my 'peak' genomic features. 

    awk 'OFS="\t" {print $1"."$2+1"."$3, $1, $2+1, $3, "."}' peaks.narrowPeak > contol.peaks.saf 

    featureCounts -a peaks.saf -F SAF -T <int> -p -o control.peaks_countMatrix.txt control.lessthan100bp.bam

So I result in two matrix outputs from feature counts:  control.peaks_countMatrix.txt; and knockdown.peaks_countMatrix.txt

I am not used to using DESeq2 for differential analysis (have been using Kallisto/Sleuth for RNA-seq). Am I right in thinking that the two matrix files can be inputted to DESeq2, to peform differential analysis against? Does anyone have a pipeline that they could share in getting this done? 




I would call peaks with the command:

    macs2 callpeak -t bamfile --nomodel --keep-dup=all

The `--nomodel` turns off the ChIPseq-specific shifting model. Do not use the p-value argument, as correction for multiple testing prevents you from excessive numbers of false-positives. This excessive number of peaks will contain lots of regions with small counts (because of the relaxed cutoff) which lack the statistical power to call differential accessability. Better use the default q-value cutoff of 0.05 or even 0.01. As you already have the sub-100bp reads, I think the `--call-summits` option is neither necessary, nor beneficial. Most peaks will already be around 100-200bp. 

There are multiple options how to call the peaks. Either call them one the merged BAM files of one condition, or call on every replicate and then take the intersect for each condition. There is no standard for this. I would probably call peaks with the above command and a q-value of 0.01 on each single file (I hope you have replicates, if not there is no point in using DESeq2), and then take the intersect between the files of one condition. You'll end up with one peak set for each condition, which is then to be merged (`bedtools merge`) and used as input for featureCounts to count the reads. DESeq2 takes **one** count matrix which contains the counts for each single replicate over one peak set (the merge between the conditions). From there on, you can pretty much follow the [DESeq2 manual][1]. It contains all the necessary code and is outstandingly extensive and informative.


  [1]: https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html