Extracting interaction contrast when using user supplied model.matrix
Entering edit mode
idetoma • 0
Last seen 19 hours ago

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.

DESeq2 • 80 views
Entering edit mode
Last seen 4 hours ago
United States

It’s probably easiest to just use a numeric contrast, where the values will be applied to coefficients as listed by resultsNames().


Login before adding your answer.

Traffic: 310 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6