Limma Package For Differential Expression Between Control And Stress

Question

[GR]

Hi All,

I started using LIMMA for checking differential expression of genes in control and stress treatments (2 replicates for each condition). I constructed a design matrix where replicates for each one goes to the same group like group1 contains rep1 and rep2 for control and group2 contains rep1 and rep2 for stress treatment. Next, I build a contrast matrix comparing group2 with group1. Finally, I got the results from LIMMA that gives me the fold change values.

My question is how do I know that LIMMA is looking for fold change of control over stress or stress over control? In other words, how to decide whether the genes with more positive log2 fold change values are more highly expressed (overexpressed) in the stress condition or in the control condition.

Hope my query is understandable and any help would be very much appreciated.

Many thanks, RT

[Sean Davis]

You need to visualize your data. I always make a small heatmap of gene expression findings to make sure that the data do, indeed, make sense.

[Damian Kao]

Can't you just pick a few LIMMA results with high significance and compare it to your raw expression data to see what the comparison direction it is?

Answer 1

[seidel]

It depends on how you specify your design matrix. You should look in the limma users guide, as there is a section explaining design matrices. There are two options for how to specify the stress vs. control. In one you can find the difference in a single step by specifying a coefficient for the control and then a coefficient for the difference between control and stress, but the orientation of the result will depend on how you configure those coefficients. The second method takes an extra step, but is more intuitive (to me), where you specify coefficients for each group (i.e one for control, one for stress), and then in a second step you create a contrast between the groups, that you specify like so: myRatio=stress-control. In that case, the resulting ratio depends on which term you put first. The first term is the numerator, and the second term is the denominator of the ratio. For example:

design=cbind(stress=c(1,1,0,0), control=c(0,0,1,1))
fit <- lmFit(eset, design)
cont.matrix <- makeContrasts(STvsCO=stress-control, levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)

The nice thing about this approach is that you get values for the averaged arrays in the fit object, and the ratio in the fit2 object so you can easily examine them separately for various properties. But as DK points out, whatever you do, perform a sanity check by looking at values to see that they make sense. (and read the user guide carefully - it's full of info! http://bioconductor.org/packages/2.9/bioc/vignettes/limma/inst/doc/usersguide.pdf )