Creating a correct limma model with respect to contrast matrix and design matrix
1
0
Entering edit mode
@payam_delfani1980-17107
Last seen 4.2 years ago

## 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

##---------------------

# Making all contrasts

(cont.matrix_adj = t(makeAllContrasts(design_adj, Condition))) ###Question: is this a correct matrix?

# Compute Contrasts from Linear Model Fit:

##---------------------

# call eBayes

##---------------------

# 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)

3
Entering edit mode
Aaron Lun ★ 28k
@alun
Last seen 2 hours ago
The city by the bay

I don't know where you got makeAllContrasts from, that's not a limma function. I assume that you're using the function from the EMA package, but I can't say that it's giving you what you want. I would do it manually:

# Making up an example design.
Cond <- gl(4, 4)
Batch <- factor(rep(1:4, 4))
design <- model.matrix(~0 + Cond + Batch)

# Making a contrast for each pairwise comparison.
cont.set <- list()
for (i in seq_len(nlevels(Cond))) {
for (j in seq_len(i-1L)) {
new.cont <- integer(ncol(design))
new.cont[i] <- 1
new.cont[j] <- -1
cont.set <- c(cont.set, list(new.cont))
}
}
cont.mat <- do.call(cbind, cont.set)


... and use that in contrasts.fit.

0
Entering edit mode

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 of makeAllContrasts is that it creates the contrasts while keeping the original name of annotations.