Question: limma paired t test
0
gravatar for xaimi0813
11 months ago by
xaimi08130
xaimi08130 wrote:

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)
R • 1.1k views
ADD COMMENTlink modified 11 months ago by Kevin Blighe63k • written 11 months ago by xaimi08130

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

ADD REPLYlink written 11 months ago by RamRS28k
1
gravatar for Kevin Blighe
11 months ago by
Kevin Blighe63k
Kevin Blighe63k wrote:

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

ADD COMMENTlink written 11 months ago by Kevin Blighe63k

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

ADD REPLYlink written 11 months ago by xaimi08130

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

ADD REPLYlink written 11 months ago by Kevin Blighe63k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1360 users visited in the last hour