I have two states and two cohorts that create 4 groups (that I call type in the model. Inside each of this group I have individuals nested within group.
I created a variable called nested_subject which distinguishes the individuals nested within a group and I created my model as explained here: http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#group-specific-condition-effects-individuals-nested-within-groups
Since in my case I have an unbalanced design, I supplied a model matrix as explained in the DESeq2 vignette.
dds <- DESeqDataSetFromMatrix(countData = counts,
colData = metadata,
design = ~ 1) ### set a simple model to avoid full rank problem (it is an unbalanced design)
mod.matrix=model.matrix(~ type + type:nested_subject, as.data.frame(colData(dds))) # subject distinguish the individual nested within a group
mod.matrix=mod.matrix[,which(colMeans(mod.matrix)!=0)] #remove column that are all zero (if any)
dds <- DESeq(dds, betaPrior = FALSE, full=mod.matrix)
Type can have 4 values:
> levels(dds$type)
[1] "healthy.G1" "healthy.G2" "disease.G1" "disease.G2"
Therefore this is the resultsName output (i am not printing the coefficient relative to the nested subjects)
> resultsNames(filter)[1:4]
[1] "Intercept" "typehealthy.G2" "typedisease.G1" "typedisease.G2"
Now I would like to calculate the contrasts disease - healthy in each of the two group, and this can be easily done like this:
disease_vs_healthy_G2=results(dds, contrast = list("typedisease.G2", "typehealthy.G2" ))
disease_vs_healthy_G1= results(dds, name = "typedisease.G1")
However, how could I calculate the interaction contrast: (disease.G2 - healthy.G2) - (disease.G1 - healthy.G1)
I tried
interaction <- results(dds, contrast = list(c("type47XXY.Saudi"),c("type46XY.Saudi","type47XXY.European_NorthAmerican")), listValues = c(1,-1/2))
But it does not seem to give the expected results.
