I am trying fit a linear model with limma, where I do adjust for a batch variable. Then I want o calculate statistics for all pair-wise comparisons. Here, I'm wondering whether the way I create the model and its contrast-fits make sense and looks correct.
Could anyone consider the example below and give me a feedback on that?
# Fit many adjusted statistics with limma. Here we adjust for the day effect (a surrogate for a batch effect):
##---------------------
#Covariate of interest:
Condition <- as.factor(pData(ExpressionSet$Condition)
#Adjusting variable:
Batch <- (pData(ExpressionSet)$Scan_date)
##---------------------
# The adjusted model
(design_adj <- model.matrix(~0+ Condition + as.factor(Batch))) ##Zero: do not include the intercept
##---------------------
# Fitting the model
fit_adj = lmFit (ExpressionSet, design_adj)
##---------------------
# Making all contrasts
(cont.matrix_adj = t(makeAllContrasts(design_adj, Condition))) ###Question: is this a correct matrix?
# Compute Contrasts from Linear Model Fit:
fit2_adj= contrasts.fit(fit_adj, cont.matrix_adj)
##---------------------
# call eBayes
fit2_adj= eBayes(fit2_adj)
##---------------------
# Write all results:
write.fit(fit2_adj, results = NULL, file="adj_Limma_global.txt", digits = 5,
adjust = "BH", method = "global", F.adjust = "none",
quote = FALSE, sep = "\t", row.names = TRUE,col.names=NA)
Dear Aaron,
I compared the approach you suggested here to that of the
makeAllContrasts
from EMA package, and can say that both method generate same outputs. However, the advantage ofmakeAllContrasts
is that it creates the contrasts while keeping the original name of annotations.