Hello guys,
One question I have met recently is that when i handle the GEO data(GSE100186) with limma package, the result by t~test I could get was about 1,561 different express genes. However ,when I chose paired t~test, I could get none genes when I set p.value= 0.05. I did not know what wrong with this or how could I deal with this case. The process I operated is as follows: t~test:
group_list<- c(rep("normal",4), rep("Tumor",4))
group_list<- factor(group_list, levels = c("Tumor","normal"))
group_list<- factor(group_list, levels = c("Tumor","normal"), ordered = F)
library("limma")
design<- model.matrix(~0+group_list)
colnames(design)<- levels(group_list)
contrast.matrix<- makeContrasts(Tumor-normal, levels = design)
fit<-lmFit(exprSet, design)
fit2<- contrasts.fit(fit, contrast.matrix)
fit2<- eBayes(fit2)
allDiff<- topTable(fit2, adjust="fdr", coef = 1, number = Inf)
paired~T test:
pairinfo = factor(rep(1:4,2))
group_list<- c(rep("tumor",4), rep("normal",4))
group_list<- factor(group_list,levels = c("tumor","normal"),ordered = F)
design=model.matrix(~ pairinfo+group_list)
fit=lmFit(exprSet,design)
fit=eBayes(fit)
allDiff_pair=topTable(fit,adjust='BH',coef=2,number=Inf,p.value=0.05)
•
written
by
xaimi0813 •
You are taking the incorrect coefficient for the paired test, I believe.
I re-ran your code, creating my own compatible random data matrix with exprSet <- matrix(rexp(200, rate=.1), ncol=8), and have the following design for the paired test:
design
(Intercept) pairinfo2 pairinfo3 pairinfo4 group_listnormal
1 1 0 0 0 0
2 1 1 0 0 0
3 1 0 1 0 0
4 1 0 0 1 0
5 1 0 0 0 1
6 1 1 0 0 1
7 1 0 1 0 1
8 1 0 0 1 1
So, you will want to use the following for the paired test:
allDiff_pair <- topTable(
fit,adjust = 'BH',
coef = 5,
number = Inf)
When deriving the test statistics for the 5th coefficient, an adjustment will be made to control for the tumor-normal pairing. You may still end up with nothing passing FDR Q <= 0.05, though.
Kevin
Thanks. I can not agree you more. Before i run the limma package, i do not know how to choose coef value or i did not fully understand the coef. And in this case, i only have two group~ tumor and normal, and if i chose t-test, the coef.value could be 1. So could you please explain more about coef in limma?Thanks. Aiming XU
Hey, here is what the entry in the function documentation says:
coef: column number or column name specifying which coefficient or
contrast of the linear model is of interest. For ‘topTable’,
can also be a vector of column subscripts, in which case the
gene ranking is by F-statistic for that set of contrasts.
So, it is the column number or column name from design
Log in to answer this question.
Please use the formatting bar (especially the
codeoption) to present your post better. You can use backticks for inline code (`text` becomestext), or select a chunk of text and use the highlighted button to format it as a code block. I've done it for you this time.