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

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
ADD COMMENT
0
Entering edit mode
@mikelove
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().

ADD COMMENT

Login before adding your answer.

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

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

Powered by the version 2.3.6