I need to manually define contrast in a GLM model implemented in DESeq2 model which performs differential expression analysis on read count data of genes using GLM.
I would like to examine which genes have significantly higher expression than the grand mean at each level of a categorical variable let’s say condition2. But while I am doing this, I want to control the effect of condition1.
Here is an example how I define design matrix and contrast matrix
condition1 <- factor(rep(c("A", "B", "C"), each=6))
a <- sample(rep(1:6, each=3))
condition2 <- factor(a, levels=1:6)
mm <- model.matrix(~ condition1 + condition2,
contrasts=list(condition1 = contr.treatment,
condition2 = contr.sum))
First question: does this design matrix and contrast do what I want: control the effect of condition1 while calculates differences between the levels of condition2 and the overall/grand mean of condition2.
Next, I compile an expression data and perform differential expression analysis where I need to give a contrast vector.
library(DESeq2)
count_data <- matrix(sample(1:100, 50 * 18, replace=TRUE), ncol=18)
dds <- DESeqDataSetFromMatrix(count_data, colData=data.frame(condition1, condition2),
design=mm)
dds2 <- DESeq(dds)
res <- results(dds2, contrast=c(0, 1, 0, 1, 0, 0, 0, 0))
Second question: in the last line I need to give a contrast vector which means one numeric value for each coefficient. In this design matrix I have 8 coefficients: Intercept, 2 coefficients for condition1 and 5 coefficients for condition2. The two contrast matrices look like:
condition1:
2 3
A 0 0
B 1 0
C 0 1
condition2:
[,1] [,2] [,3] [,4] [,5]
1 1 0 0 0 0
2 0 1 0 0 0
3 0 0 1 0 0
4 0 0 0 1 0
5 0 0 0 0 1
6 -1 -1 -1 -1 -1
Am I right if I put together the two contrast matrices I can give 0, 1, 0 for the intercept and condition1 and 1, 0, 0, 0, 0 for condition2, so c(0, 1, 0, 1, 0, 0, 0, 0)? Am I right that this contrast and design would give me the difference between condition2 level1 and the grand mean of condition2 while the effect of condition1 is controlled?
Thank you for your help!
Torda
