Can I model technical replicates in DESeq2?

Question

I’d normally use collapseReplicates (or do the collapsing upstream) to handle technical replicates.

However, in my current RNA-seq experimental design, samples were sequenced twice using different library preparation protocols, leading to marked differences in the resulting count estimates (in fact, the choice of different protocol explains most of the variance in the bulk data according to PCA). I would therefore like to model these differences by adding them as a covariate to the DESeq2 model.

I thought I could just add another factor to the design formula, equivalent to the “pasilla” example in the DESeq2 vignette, which adds the factor type for the library type (single-end vs paired-end). However, the pasilla vignette states that these different libraries are actually independent biological replicates, not technical replicates as I had previously assumed.

I’m now concerned that having technical replicates in the design might artificially inflate the power. Can the design be salvaged?

It’s worth noting that we also have biological replicates, and a full rank design; i.e. every combination of treatment and library prep, with biological replicates.

Answer 1

If they're truly technical replicates, then there's no way to model them using DESeq2*, as you've alluded to with the collapseReplicates function. DESeq2/ Mike Love's general recommendation with collapseReplicates is to just add the reads together for technical replicates.

If you want to model them instead of collapsing them down, you can voom transform your data and follow an example similar to section 11.3 in the Limma users guide.

(Code example from 11.3):

Design

| FileName  | Cy3   | Cy5  |
| --------- |:-----:| ----:|
| File1     | wt1   | mu1  |
| File2     | wt1   | mu1  |
| File3     | wt2   | mu2  |
| File4     | wt2   | mu2  |

Example

> biolrep <- c(1, 1, 2, 2)
> corfit <- duplicateCorrelation(MA, ndups = 1, block = biolrep)
> fit <- lmFit(MA, block = biolrep, cor = corfit$consensus)
> fit <- eBayes(fit)
> topTable(fit, adjust = "BH")

*Edit: As @Devon Ryan mentions, you can design a paired model in DESeq2, but beware of making this full rank.

Answer 2

Practically speaking, there's no way to include the technical replicates in that design (in DESeq2 at least). Your concern regarding inflating the power is exactly correct and the only way to combat that would be to add a pairing factor like one might do with case-control or tumor-normal studies. That is, something like:

  group libraryPrep sample
1    WT           A      1
2    WT           B      1
3   MUT           A      2
4   MUT           B      2
5    WT           A      3
6    WT           B      3
7   MUT           A      4
8   MUT           B      4

Here sample pairs the technical replicates together. However this ends up either being rank deficient or in the end not more informative.