limma paired t test

Question

[xaimi0813]

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)

[Ram]

Please use the formatting bar (especially the code option) to present your post better. You can use backticks for inline code (`text` becomes text), 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.

Answer 1

[Kevin Blighe]

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

[xaimi0813]

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

[Kevin Blighe]

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