We have a project with a somewhat spectacular design: two cell types in
trans-well culture, each cell population from a human donor.  In each
experimental setup, we use one set of donors, treated with ~10 different
perturbagens performed as a single replicate, plus a vehicle control
performed in duplicate; we repeat this setup independently 5 times with 5
different sets of donors (same perturbagens).  Then, in another batch, we
repeat the entire process with a new set of perturbagens (again with
vehicle controls, though same donors as before).  At the end of the day (or
months, as it may be), we have many hundreds of samples across ~8 batches,
5 sets of cell donors, and ~80 different perturbagens, and we'd of course
like to test for the effects of each individual perturbagen, and
compare/contrast effects across/between perturbagens.

Using edgeR, voom, and DESeq we have analyzed this data two ways

  1) using all of the samples, accounting for a variety of intermingled
batch and donor effects, and modeling all 80 perturbagens simultaneously

  2) using only the small subset of samples relevant to a given perturbagen
(i.e. the 5 treated samples and their 10 matching vehicle controls), and
modeling only 4 donor effects.

When we look for subset of known/expected biological signals across
well-understood perturbagens, we are surprised to find that approach #2
performs much better than #1 (both with respect to sensitivity at any given
FDR, and with respect to statistical ranking, regardless of FDR), despite
our prior expectation that #2 will lack power compared to #1 (which has
many more df than #2, and thus much more accurate tag-level dispersion
estimates).  Our working hypothesis is that the additional batch
corrections required for #1 introduce enough additional uncertainty that
while the same "correct" fold changes are observed by each approach,
statistical significance is rarely achieved with #1.

A second concern is that our dataset as a whole violates the the constant
dispersion assumption inherent to all of the methods, and that by refitting
each subset of the data in #2 we are in effect relaxing that assumption; if
so, I wonder if it would be useful to consider using the whole-dataset
dispersion estimates from #1 as priors when refitting the subsets, to
regain some of the expected power.

A third issue with both approaches is that the same vehicle controls are
reused for multiple perturbagens within a batch, leading to a nested
correlation structure not accounted for in either of the approaches; the
net effect of this is that the estimated perturbagen-induced fold changes
within each batch are correlated with each other, leading to difficulties
when comparing responses within and between batches.  I presume we could
try to address this in the voom pipeline using duplicateCorrelation() to
get a mixed model, but this doesn't help us with edgeR or DEseq.

In sum, I'm sure we're just rediscovering well known modeling pitfalls, but
we'd appreciate any comments or suggestions.  We are currently testing
whether modifying approach #1 to include only donor-level batch
coefficients will recover some accuracy, but any further ideas would be
appreciated.

Thanks,
-Aaron

	[[alternative HTML version deleted]]

