DESeq2 with nested design?

[aliceb]

Hi all,

Is is possible to use a nested model in DESeq2 or EdgeR??

I have 3 treatments (evolution lines), but each treatment was performed in 4 replicate lines.

My experimental design is as follows:

TmtA x 4 replicate lines X 3 biological replicates sequenced in each replicate line
TmtB x 4 replicate lines X 3 biological replicates sequenced in each replicate line
TmtC x 4 replicate lines X 3 biological replicates sequenced in each replicate line

This has given me a total of 36 libraries.

My main interest is the difference between treatments, but I also want to take the independence of these 4 replicate lines into account while keeping them nested within their treatment. Is there an easy way to incorporate this in DESeqDataSeqFromMatrix (or similar, this is what I'm using now)?

Thanks.

Alice

[dpryan]

Nesting implies a mixed-effect model, so no, neither DESeq2 nor edgeR will support that. There's some limited support for that in limma through duplicateCorrelation, which you might be able to use.

Having said that, you can simply use the nesting factor (I assume this is the "replicate line") as a blocking factor. In practice, this works quite well (the GLMs allow you to moderate dispersion anyway).

[aliceb]

Thanks for the quick response!

Yes, I think I've tried that using: "design= ~ sel.line + sel.line:block"

However, I'm worried that I'm not correctly setting up my design. I have also tried "design= ~ block + block:sel.line", but that didn't seem right.

The troubling thing is that in both of the designs above, there is virtually no overlap between the significant genes identified in this and an unblocked run (using "design = ~ sel.line"). And when I look at the results of the first two scenarios above, I can see that they are driven by "blocks" more than sel.line (i.e. similar responses in blocks across sel.lines). Does that seem likely? Or am I setting this up wrong?

I've attached my script here in case that's helpful.

Thanks!

[dpryan]

You probably want either "~block+sel.line" or ~block*sel.line". I suspect the former will work better, but perhaps you have a number of block-specific interactions present. The former is what we would use for human experiments, where the patient to patient variability can be quite high and we don't have access to the within-block factor replicates that you have.

[aliceb]

Ah yes, that has much more overlap with the unblocked design than the other two.

Thanks very much!

[Michael Love]

I'd second Devon's suggestion of using "~ block + sel.line". Then use

results(dds, contrast=c("sel.line","GroupE","GroupD"))

etc. to build results tables for comparisons of the three treatments.

[aliceb]

Ah, thank you! I knew I was missing something in the results