This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Confusion over what is reference for design matrix in limma

Hi, I need some confirmation to check I am getting this right as i have been reading examples and it is just adding to the confusion!

I am doing a case versus control study for differential expression analysis using Limma package in R. My control group is "0" and my case group is "1". I am wondering what is my reference for comparison between case and control? So when i look at the output topTable FC values how do i know what is under or over expressed?

Here is my class label:

> class
 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0
Levels: 0 1

So the above would imply to me that i am making the comparison between control (0) and case (1). So if i consider the first entry in my top table:

> topTable(resultmulti1,coef = 2, n = 1)

     logFC  AveExpr        t      P.Value    adj.P.Val        B
59285 3.503369 1.645189 87.00804 1.614614e-32 2.340221e-28 63.84269

Where the FC is 3.50033 is this fold change with reference to the control? The reason i ask is that if i look at the mean expression for this gene for the two class labels:

> typeMean <- tapply(expression_iGiS[genename,], phenodata$class, mean)

> typeMean

         0          1 
-0.4548455  3.0452115

This would imply that my class 1 is the reference as if you take away 0-1 this gives - 3.50033 but if you do 1-0 this gives you +3.50033 (which is in the fold change column.

The only thing i can think of is that my design matrix is influencing what is my reference: The deisgn matrix is made:

>phenodata <- as.data.frame(cbind(c(rep(1, 15),rep(0 ,10))

>colnames(phenodata) <- c("class")

>design <- model.matrix(  ~ class  )

This create a matrix which includes an intercept (which is all 1's) and my class variable. Is the intercept being set as 1 effecting my model and what is being used as a reference?

Any comments would be really helpful! I want to know if i can say for example gene 59285 is upregulated in control compared to case. From the lack of clarity in examples i am unsure if this is the right away around. The resource i used was:http://kasperdanielhansen.github.io/genbioconductor/html/limma.html

Thanks Danielle

r limma designmatrix

1 answer

0 is the reference, since it's lexicographically first (and given first in levels(class)). This also makes sense given that the example fold-change is basically mean(1) - mean(0).

Log in to answer this question.