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