Question: Creating a correct limma model with respect to contrast matrix and design matrix
0
gravatar for payam_delfani1980
7 months ago by
payam_delfani19800 wrote:

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)

ADD COMMENTlink modified 7 months ago by Aaron Lun23k • written 7 months ago by payam_delfani19800
Answer: Creating a correct limma model with respect to contrast matrix and design matrix
3
gravatar for Aaron Lun
7 months ago by
Aaron Lun23k
Cambridge, United Kingdom
Aaron Lun23k wrote:

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 COMMENTlink modified 7 months ago • written 7 months ago by Aaron Lun23k

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 REPLYlink written 7 months ago by payam_delfani19800
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 138 users visited in the last hour