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

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)

limma limma contrast matrix limma design matrix • 1.3k views
ADD COMMENT
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.

ADD COMMENT
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.

ADD REPLY

Login before adding your answer.

Traffic: 406 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6