DESeq2 question: complex multifactor experiemnt

[jmgrindheim]

Hello,

I have a complex experimental setup.

What I've done:

countTable=read.table("HTSeq_Table.txt", header=T,row.names=1)

And I have a design like this, where the sample names are the row names.

> design
geno age libType
S1534 Ezh1 P14 single
S1536 Ezh1 P14 single
S8633 Ezh1 P14 single
S1532 Ezh12 P14 single
S8631 Ezh12 P14 single
S1141 Ezh12del P14 single
S1142 Ezh12del P14 single
S1541 Wt P14 single
S1547 Wt P14 single
S8Wrep1 Wt W8 paired
S8Wrep2 Wt W8 paired
SE18rep1 Wt E18 paired
SE18rep2 Wt E18 paired
P0.expt1.bio1 Wt P0 single
P0.expt1.bio2 Wt P0 single
P0.expt2.bio1 Wt P0 single
P0.expt2.bio2 Wt P0 single

So, I have multiple timepoint and single/paired-end sequencing, but only one timepoint has multiple genotypes.

To make a DESeqCountDataSet, I ran

> cds=DESeqDataSetFromMatrix(countData=countTable,colData=designnogroup, design=~geno+age+libType)
Error in DESeqDataSet(se, design = design, ignoreRank) :
the model matrix is not full rank, so the model cannot be fit as specified.
one or more variables or interaction terms in the design formula
are linear combinations of the others and must be removed

1: I'm not sure how to deal with having a timecourse experiment where only one timepoint has multiple genotypes.
2. Do I need to include single/paired-end?

Thank you so much!

Addendum: I found the answer. In the manual of course.

3.12.1 Linear combinations

[dpryan]

Without playing around with the matrix it looks like "libType" is confounding the estimation of the "W8" and "E18" ages. You're going to have to either realign the PE data as SE (just ignore read 2) or accept that W8 and E18 estimations might be confounded by a batch effect (they likely will be regardless, though this will help minimize that). So ~geno+age instead of ~geno+age+libType.

[jmgrindheim]

Thanks Devon, that was really helpful.

Also, do you know what the effect will be if I put ~age+geno vs ~geno+age? It's kinda hard to decide what's more important. I want to see if the genotype causes gene expression to look like a different timepoint, so besides differential expression, I really want normalized count values to create a heatmap and I don't see an easy way to do that in DESeq2. In the first DESeq, I would take baseMeanA or B values from a differential expression test, but I haven't yet found an equivalent with DESeq2.

[dpryan]

The order only affects plotting and what's output by results() by default (it'll default to whichever you specify last), the actual statistics will be the same. For normalized counts, just use counts(dds, normalized=T). That's vastly more meaningful than creating heatmaps with the group means.

[jmgrindheim]

So you think that plotting count values for individual biological replicates is better then plotting them for the replicates grouped?

[dpryan]

Yes, you don't lose the variance that way. If you have a crazy number of groups/samples then that might not be feasible, of course (though then the heatmap likely won't tell you much anyway).