Hi,
I am sorry if this is a repeated question, but I couldn't find any example that fits my exact case.
I have a dataset of individuals with alternative genotypes (North or Heterozygous). Each genotype was reared under different temperatures (Seven or Ten). I am interested in understanding the effect of temperature by genotype and of genotype by temperature. Initially I constructed the following model matrix in edgeR following section 3.3.1:
design$group<-factor(paste0(design$TEMP, ".", design$GENO))
design_model <- model.matrix(~ 0 + group)
colnames(design_model) <- levels(group)
SEVEN.Heterozygous SEVEN.North TEN.Heterozygous TEN.North
1 1 0 0 0
2 1 0 0 0
3 1 0 0 0
4 1 0 0 0
5 1 0 0 0
6 1 0 0 0
7 1 0 0 0
8 1 0 0 0
9 1 0 0 0
10 1 0 0 0
11 1 0 0 0
12 1 0 0 0
13 0 1 0 0
14 0 1 0 0
15 0 1 0 0
16 0 1 0 0
17 0 1 0 0
18 0 1 0 0
19 0 0 1 0
20 0 0 1 0
21 0 0 1 0
22 0 0 1 0
23 0 0 1 0
24 0 0 1 0
25 0 0 1 0
26 0 0 1 0
27 0 0 1 0
28 0 0 1 0
29 0 0 1 0
30 0 0 0 1
31 0 0 0 1
32 0 0 0 1
33 0 0 0 1
34 0 0 0 1
35 0 0 0 1
36 0 0 0 1
37 0 0 0 1
38 0 0 0 1
attr(,"assign")
[1] 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"
And testing differential expression for the effect of, let's say, genotype for a given temperature:
# Calculate dispersion
df <- estimateDisp(df, model_design, robust=TRUE)
df <- estimateGLMCommonDisp(df, model_design)
# Fit the model
fit <- glmQLFit(df, model_design, robust=TRUE)
# Make a contrast
con.SEVEN.North_vs_Het <- makeContrasts(groupSEVEN.North - groupSEVEN.Heterozygous, levels=design_model)
qlf.SEVEN.North_vs_Het <- glmQLFTest(fit, contrast=con.SEVEN.North_vs_Het)
Now I wanted to make this a bit more complicated:
1. Testing for overall effect of temperature or genotype
With the model matrix above, is it correct to test for the overall effect of temperature like this?
qlf.Ten_vs_Seven <- glmQLFTest(fit, contrast=c(-1,-1,1,1)
Or should I instead construct a model matrix that includes temperature and genotype explicitly (see code bellow)? But then how do I test for the effect of genotype within temperature, or temperature within genotype, like above; doo I need to have a separate model matrix for this interaction?
TEMP<-factor(paste(design$TEMP))
GENO<-factor(paste(design$GENO))
design_model <- model.matrix(~ GENO+TEMP)
2. How to add a blocking factor?
Sex doesn't seem to have a strong effect on gene expression in this dataset (all groups of GENO x TEMP have equivalent numbers of males and females), but I thought I should still control for this effect. What would be the best way to do this in the model matrix? From the first model matrix, should I just add SEX as a blocking factor in the design matrix above like so:
design$group<-factor(paste0(design$TEMP, ".", design$GENO))
SEX<-factor(paste(design$SEX))
design_model <- model.matrix(~ 0 + group + SEX)
And then proceeding as above to make the contrasts that I want? I am able to fit the model with this model matrix, but I am not sure it is the correct way of adding the blocking factor.
So I also tried this code following the section 3.3.2 of the manual and inspired by this question, but I get a warning that my Design matrix is not of full rank
SEX<-factor(paste(design$SEX))
TEMP<-factor(paste(design$TEMP))
GENO<-factor(paste(design$GENO))
design_model <- model.matrix(~ 0 + SEX + GENO:TEMP)
SEXF SEXM GENOHeterozygous:TEMPSEVEN GENONorth:TEMPSEVEN GENOHeterozygous:TEMPTEN GENONorth:TEMPTEN
1 0 1 1 0 0 0
2 1 0 1 0 0 0
3 0 1 1 0 0 0
4 1 0 1 0 0 0
5 0 1 1 0 0 0
6 0 1 1 0 0 0
7 1 0 1 0 0 0
8 1 0 1 0 0 0
9 0 1 1 0 0 0
10 1 0 1 0 0 0
11 0 1 1 0 0 0
12 0 1 1 0 0 0
13 1 0 0 1 0 0
14 0 1 0 1 0 0
15 1 0 0 1 0 0
16 1 0 0 1 0 0
17 1 0 0 1 0 0
18 0 1 0 1 0 0
19 1 0 0 0 1 0
20 0 1 0 0 1 0
21 0 1 0 0 1 0
22 0 1 0 0 1 0
23 1 0 0 0 1 0
24 1 0 0 0 1 0
25 0 1 0 0 1 0
26 1 0 0 0 1 0
27 1 0 0 0 1 0
28 1 0 0 0 1 0
29 0 1 0 0 1 0
30 1 0 0 0 0 1
31 0 1 0 0 0 1
32 1 0 0 0 0 1
33 0 1 0 0 0 1
34 1 0 0 0 0 1
35 0 1 0 0 0 1
36 0 1 0 0 0 1
37 0 1 0 0 0 1
38 0 1 0 0 0 1
attr(,"assign")
[1] 1 1 2 2 2 2
attr(,"contrasts")
attr(,"contrasts")$SEX
[1] "contr.treatment"
attr(,"contrasts")$GENO
[1] "contr.treatment"
attr(,"contrasts")$TEMP
[1] "contr.treatment"
Then I tried explicitly adding GENO and TEMP:
SEX<-factor(paste(design$SEX))
TEMP<-factor(paste(design$TEMP))
GENO<-factor(paste(design$GENO))
design_model <- model.matrix(~ SEX + GENO + TEMP + GENO:TEMP)
(Intercept) GENONorth TEMPTEN SEXM GENONorth:TEMPTEN
1 1 0 0 1 0
2 1 0 0 0 0
3 1 0 0 1 0
4 1 0 0 0 0
5 1 0 0 1 0
6 1 0 0 1 0
7 1 0 0 0 0
8 1 0 0 0 0
9 1 0 0 1 0
10 1 0 0 0 0
11 1 0 0 1 0
12 1 0 0 1 0
13 1 1 0 0 0
14 1 1 0 1 0
15 1 1 0 0 0
16 1 1 0 0 0
17 1 1 0 0 0
18 1 1 0 1 0
19 1 0 1 0 0
20 1 0 1 1 0
21 1 0 1 1 0
22 1 0 1 1 0
23 1 0 1 0 0
24 1 0 1 0 0
25 1 0 1 1 0
26 1 0 1 0 0
27 1 0 1 0 0
28 1 0 1 0 0
29 1 0 1 1 0
30 1 1 1 0 1
31 1 1 1 1 1
32 1 1 1 0 1
33 1 1 1 1 1
34 1 1 1 0 1
35 1 1 1 1 1
36 1 1 1 1 1
37 1 1 1 1 1
38 1 1 1 1 1
attr(,"assign")
[1] 0 1 2 3 4
attr(,"contrasts")
attr(,"contrasts")$GENO
[1] "contr.treatment"
attr(,"contrasts")$TEMP
[1] "contr.treatment"
attr(,"contrasts")$SEX
[1] "contr.treatment"
And this works with fitting the model, but then I am not sure how to test for differences between temperatures within genotypes, or differences between genotypes within temperatures because I only have the GENONorth:TEMPTEN, and not GENONorth:TEMPSEVEN, GENOHeterozygous:TEMPTEN, or GENOHeterozygous:TEMPSEVEN...
Any advice?
Thank you so much!
