EDIT:

Cross-posted on Bioconductor: https://support.bioconductor.org/p/114506/

---------

Hello,

I am trying to figure out how to handle biological replicates across batches in DESeq2's design matrix. They're not quite technical replicates, but trying to include the sample ID as a factor in the design matrix produces a 'Model matrix not full rank' error. What do you think I should do?

# Experiment:

We have three variables of interest in our analysis: Timepoint (A and B), Sex (M or F), and Condition (cond1 or cond2). 

Note that a given sample can only have one timepoint, one sex, and one condition - so it's not like we can split the same tissue and look at it at Timepoint A and B. Each possible Timepoint x Sex x Condition has its own unique samples.

In our first sequencing batch, we collected samples for each possible combination of conditions.  

In our second batch, we took some of the same RNA samples from the first sequencing batch PLUS some new RNA samples, re-generated libraries from all of these, and then sequenced.

In the end, we have a sample table that looks like this:

    SAMPLE | BATCH | CONDITION | SEX | TIMEPOINT
    ------  -------  ---------   ----  ---------
    samp1  |  1    |   cond1   |  F  |    A
    samp1  |  2    |   cond1   |  F  |    A
    samp2  |  1    |   cond1   |  M  |    A
    samp2  |  2    |   cond1   |  M  |    A
    samp3  |  1    |   cond1   |  F  |    A
    samp4  |  2    |   cond1   |  F  |    A

...and so forth.

# Problem

I tried to run DESeq with this:

    design(dds) <- ~Genotype + Sex + Condition + SequencingBatch + Sample

and also with this:

    design(dds) <- ~Genotype + Sex + Condition + SequencingBatch_Sample

...where SequencingBatch_Sample is a combination of SequencingBatch and Sample, like "samp1_1", samp1_2" for sample 1/batch1 and sample1/batch2.
 

These designs, of course, throw a 'Matrix not full rank' error, presumably because it is not possible to have the same sample have every possible combination of Condition/Sex/Timepoint, since each sample can only have one quality each of Condition, Sex, and Timepoint.

However, I'm not quite sure what to do instead.

Using collapseReplicates() doesn't feel appropriate because these aren't really technical replicates - a whole new library was produced for the second batch, just using the same RNA.

Doing nothing seems problematic because then the 'duplicated' signal from a subset of samples threatens to overwhelm the differential expression results.

All suggestions are very welcome - thank you!



Hello,

Back-propagating the answer that Mike Love suggested on Bioconductor forums:

If there is not a lot of technical variation across replicates > consider treating the samples like technical replicates, and collapse them in DESeq using the collapseReplicates() function.

If there is a lot of technical variation across replicates (as in my case) - "even though there is no biological variation that you would expect" -

"The only approach I can think of, to inform the model that you have technical replicates nested within condition which are not biological replicates, and yet also not similar enough to justify collapsing them, is to use the duplicateCorrelation() approach with limma-voom. So I'd take a look at that framework for approaching this experiment. 

It's not possible to use fixed effects to account for the duplicate nature of technical replicates within a condition, and this is the way that DESeq2 controls for variables. You can do within-individual comparisons for various complex designs, but you can't "control" for the correlations among certain related samples within a condition and then compare directly across condition with our implementation. "