I'm fairly new to limma and just want to double check if I am doing something correctly. I have some microarray data and am looking for genes differentially expressed between states DIS and TE. I have several continuous covariates (V0, V1, V2 etc. in the phenotype data cov.tb) I am only interested in DIS/TE and would like to adjust for these covariates, excluding their effect.
I have looked into it and examined several examples but I would like a more experienced pair of eyes to confirm whether or not I have done this correctly. Is it really as simple as:
> mod1 <- model.matrix(~ 0 + mRNA_labels + V0 + V1 + V2 + V3 + V4, data=cov.tb)
Below is my entire set of commands:
> mRNA_labels <- factor(
c(rep('CNT', 17), rep('DIS', 21), rep('TE', 28)),
levels = c('CNT', 'DIS', 'TE')
)
> mod1 <- model.matrix(~ 0 + mRNA_labels + V0 + V1 + V2 + V3 + V4, data=cov.tb)
> colnames(mod1)[1] <- "CNT"
colnames(mod1)[2] <- "DIS"
colnames(mod1)[3] <- "TE"
colnames(mod1)[4] <- "V0"
colnames(mod1)[4] <- "V1"
colnames(mod1)[5] <- "V2"
colnames(mod1)[6] <- "V3"
colnames(mod1)[7] <- "V4"
> mod1
CNT DIS TE V0 V1 V2 V3 V4
VV51.CEL 1 0 0 -1.16253985 -0.112174410 0.0095394904 -0.159846427 0.124077809
VV50.CEL 1 0 0 -0.18154195 0.148821751 0.2191634922 0.009425092 -0.089006798
VV38.CEL 1 0 0 -6.24420703 -0.274974896 0.0397183693 -0.069572763 0.213245438
...
> contrast_mod1 <- makeContrasts(
DISvTE = DIS - TE,
levels=mod1)
> contrast_mod1
Contrasts
Levels DISvTE
CNT 0
DIS 1
TE -1
V0 0
V1 0
V2 0
V3 0
V4 0
>fit <- lmFit(expression.dat, mod1)
>fit.cont <- contrasts.fit(fit, contrast_mod1)
>fit.eb <- eBayes(fit.cont)