Hi all, 

I have what may seem like as really simple questions, but being new to using R in general and very new to DEseq2 I just can not figure it out reading the manual (which is great and very informative for a biologist, thanks!)

I have a simple design of 3 donors (A, B and C) and 3 conditions (unit,eth and vitd). 

My basic question is how can I do a contrast with a multifactorial design? (IU may be getting these terms mixed up and therein may lie the problem.

If I simply import the "eth" and "vitd" data and follow the manual for MF design I can run through and get great results (with 85 DE genes):

> dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ condition)

> dds <- dds[ rowSums(counts(dds)) > 1, ]

> ddsMF <- dds

> design(ddsMF) <- formula(~ donor + condition)

> ddsMF <- DESeq(ddsMF)

etcetc

However when I try and do the same thing but also importing the "uns" sample (so now I have "uns", and the two previous "eth" + "vitd") and try and follow the contrast I get 75 DE genes. So same procedure as before but this time 

> res.vitd.eth <- vitdresults(ddsMF, contrast=c("condition","vitd","eth"))

with the ddsMF run as before.

I'm sure I am doing something so basic wrongs here.

Thank you in advance,

Reuben

There are a couple of reasons I can think of that may explain this, but I assume your 75-gene list is fairly similar in content to the 85-gene list (if not, then the most likely explanation is that the data-importing has got sample labels misaligned).

 The most likely cause is due to the way replicate variability is estimated.  When you include an extra experimental group in the dataset, regardless of whether it is involved in your contrast, the residual noise is calculated within each group and then pooled - so if you're 'uns' data is more or less noisy than your existing conditions (for genes of interest), then you'll see a different amount of noise, and that will feed into the hypothesis testing.

The other reason could be (depending on whether you're using the latest version of DESeq2), a 'shrinkage' is applied to your log-fold-changes, so that each gene's lfc is moderated towards a cross-gene average.  Again, the amount of shrinkage is calculated across all pairwise comparisons (I'm not being particularly precise here, but hopefully this gives a picture) - so if the 'uns vs X' comparisons are substantially different in magnitude than the 'eth vs vitd' comparison, then you may see differences if you include a 'phantom' experimental group in your dataset. You can avoid this by having 'results(...., betaPrior=FALSE)' in your calculation, though if you're using the latest version of DESeq2, this is the default, in which case this reason is moot.

As a side issue, you should think carefully (or get a local statistician to) about your design: you may want to analyse it differently depending on whether you have every condition for each donor, or whether each donor only presents one condition, etc.