I'm hoping to use DEseq2 to look for regions of open chromatin that differ between cell types. It seems like the general statistical framework should be highly analogous between RNA-seq and DNase-seq (or, in my case, ATAC-seq). However, I wanted to check that some of the assumptions underlying DEseq2 do in fact translate, as some of my diagnostic plots differ from what I expected from reading the (excellent) DEseq2 documentation.

Data source/background: ATAC-seq was performed on three different cell types within a single tissue (two biological replicates per cell type). Open chromatin peaks were called on independent samples (~30,000 peaks per sample), and I then calculated the per region coverage in the regions defined by the union of all peaks (~50,000 peaks) for each sample. These counts (50,00 regions x 6 samples) constitute the count table I load into DEseq2. I then ran DEseq2 using design = ~ cell_type to identify regions of differentially open chromatin across the three cell types.

Of the 50,000 peaks I examined, 23,000 pass padj < 0.1 (comparing the first two cell types). I get a similar number of differences between the second cell-type comparison, and ~2,000 differences between the third pair (these are qualitatively what I expect for these cell types, but seem quantitatively high). Here's the MA plot for the first cell-type comparison (23,000 differences):

https://www.dropbox.com/s/14nkp0bqtd318jz/Screen%20Shot%202014-09-22%20at%2011.17.48%20AM.png?dl=0

Many of these differences are small or identified in regions with low mean read counts--I can easily filter on LFC to get a more manageable list (setting lfcThreshold=1 and altHypothesis="greaterAbs" reduces this to 1034 regions). When I look at top hits, I get interesting and tissue-relevant results, so I think this is generally working. But I'm worried that getting back so many significant differences might mean my data aren't appropriate for DEseq2 (e.g. too noisy). Additional specific questions:

1) Is the DEseq2 normalization framework appropriate for open chromatin data? Specifically, the normalization by scale factors derived from the median of "pseudo-references" seems to assume similar overall distributions of reads between conditions ("most genes don't change" for RNA-seq, or "most open chromatin doesn't change" if looking at DHS peaks). Given that I get so many significant differences, is this normalization not appropriate for my data, and are there specific alternatives I should consider (total read depth seems inadequate given different signal to noise ratios between samples)? Is part of the problem that RNA-seq considers the entire transcriptome (lots of low-level background), whereas my approach to feature selection has enriched for regions that have higher "expression."

2) I expected some of the lower read depth regions to be removed by independent filtering. But the independent filtering threshold is set at 0%.

https://www.dropbox.com/s/l88z0ajku1ol5k5/Screen%20Shot%202014-09-22%20at%2011.56.21%20AM.png?dl=0

When I test for larger differences (setting lfcThreshold=1 and altHypothesis="greaterAbs"), the independent filtering threshold moves to 40%.

3) Dispersion fit: this isn't related to my concern about too many hits, but parametric fitting of my dispersion curve failed. Local fitting was successful, but I don't see any asymptotic shift at higher read depths as I've seen in other example dispersion fit plots. Is this an acceptable fit or--again--does this suggest a problem with applying DEseq2 to a different type of dataset (do my data have a weaker mean-dispersion relationship and will this undermine my results)?

https://www.dropbox.com/s/8m4xs4s1bjuo4r1/Screen%20Shot%202014-09-22%20at%2011.59.48%20AM.png?dl=0

Please let me know if I can clarify question(s) provide additional information/if I've posted this in the wrong format/forum. I consulted the posting guide--I haven't provided sample code because nothing in DEseq2 is reporting any errors, but I am happy to provide sample data/code if it will help better illustrate my question. I absolutely appreciate any input.

>sessionInfo()

R version 3.1.1 (2014-07-10)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods
[8] base

other attached packages:
[1] DESeq2_1.4.5            RcppArmadillo_0.4.400.0 Rcpp_0.11.2
[4] GenomicRanges_1.16.4    GenomeInfoDb_1.0.2      IRanges_1.22.10
[7] BiocGenerics_0.10.0

loaded via a namespace (and not attached):
 [1] annotate_1.42.1      AnnotationDbi_1.26.0 Biobase_2.24.0
 [4] DBI_0.3.0            genefilter_1.46.1    geneplotter_1.42.0
 [7] grid_3.1.1           lattice_0.20-29      locfit_1.5-9.1
[10] RColorBrewer_1.0-5   RSQLite_0.11.4       splines_3.1.1
[13] stats4_3.1.1         survival_2.37-7      tools_3.1.1
[16] XML_3.98-1.1         xtable_1.7-4         XVector_0.4.0

hi,

These are good questions. 

"1) Is the DEseq2 normalization framework appropriate for open chromatin data? Specifically, the normalization by scale factors derived from the median of "pseudo-references" seems to assume similar overall distributions of reads between conditions ("most genes don't change" for RNA-seq, or "most open chromatin doesn't change" if looking at DHS peaks)...."

It doesn't necessarily assume similar distributions, but that the median ratio will capture the size relationship. If you imagine a scenario with two samples (instead of a psuedoreference):

sample1,sample2
1,8
1,8
2,4
2,4
2,4
8,2
8,2
8,1

In this example suppose the "true" size ratio is 1:2 and the genes 3-5 are not differentially expressed (for sake of simplicity I'll just use DE instead of differentially open chromatin). Even though the majority (5/8) of the rows are DE, the median ratio method will still accurately capture the 1:2 size ratio. So many rows can be DE as long as there is some balancing of up/down, or if there are only one kind (say, "up-regulation"), that this be less than an majority of rows.

Looking at the right side of an MA plot can sometimes help to see if the size factors are well-estimated, in that the log fold changes for the rows with large counts tend to have lower standard error, and if you presume that some of these are non-DE, they should fall on the x-axis. In your MA plot the points are falling generally on the x-axis, with maybe a bit more below than above on the right side. If you have a method for finding those rows which you think are non-DE, you can use the following code:

sf = estimateSizeFactorsForMatrix( counts(dds)[nonDERows,] )
sizeFactors(dds) = sf
dds = DESeq(dds) 

Then DESeq() will use your supplied size factors instead of estimating them. I'll think a bit on some better advice on how to find these rows. RUVSeq might be useful here.

But I think increasing the lfcThreshold makes a lot of sense for this dataset. Note that the dispersion is quite low, dispersion estimates ~ 1e-2, meaning coefficient of variation 10% within cell-type for reasonably large counts. 

"2) I expected some of the lower read depth regions to be removed by independent filtering. But the independent filtering threshold is set at 0%."

If you have large fold changes and high power, then filtering won't increase the number of rejections. By increasing the threshold for LFC you reduce the power and then independent filtering kicks in.

"3) Dispersion fit: this isn't related to my concern about too many hits, but parametric fitting of my dispersion curve failed. Local fitting was successful, but I don't see any asymptotic shift at higher read depths as I've seen in other example dispersion fit plots. Is this an acceptable fit or--again--does this suggest a problem with applying DEseq2 to a different type of dataset (do my data have a weaker mean-dispersion relationship and will this undermine my results)?"

This fit looks good to me, it seems you do have a mean-dispersion relationship which is accurately captured by the local fit. The shape of the curve is definitely dependent on the different protocol. For instance, I've seen chromosome conformation capture sequencing datasets where there was not much dependence, and where the "mean" fitType could be used.