请教一个基因芯片分析中的t检验的问题

fengl

对于同一批数据，先用mas5处理，再用t检验分析差异表达基因，一切正常。但用rma处理，再用t检验，死活有问题。而且错误提示还不同。

library(affy)<br />
LM2.PA=ReadAffy()<br />
LM2.PA.mas5=mas5(LM2.PA)<br />
eSet <- exprs(LM2.PA.mas5)<br />
p.value.all.genes = apply(eSet,1,function(x){t.test(x[1:2],x[3:5]) $p.value})<br />

用mas5不出错，但用rma代替mas5就会出错。

错误提示，

错误在t.test.default(x[1:2], x[3:5]) : 数据是恆量

但我用str()查看分别用mas5和rma得到的eSet，看起来都一样啊。请高手指点。谢谢了。

pengchy

错误提示说：数据是恒量：

<br />
> t.test(rep(1,4),rep(3,5))<br />
Error in t.test.default(rep(1, 4), rep(3, 5)) :<br />
  data are essentially constant<br />

这个返回是正常的吧，说明rma处理以后，有一些基因在各个sample的表达值是一样的（可能都为0，或者最小值），你可以在做t检验之前先筛选一下，比如，去掉在所有sample中表达值都低于背景值（或者某个值）的基因。

fengl

谢谢提醒。我用了个笨办法，用subset把数据分别拆出来看，看到哪一行，t.test开始报错。后来发现，在第614行开始报错，原因是rma算法会对原始值进行数据校正（废话），会有背景去除，scaling，log2转换等。我想关键就是log2转换，导致数据变小，可能因为精度的原因，当t.test认为某一行的所有数据都相同的时候（并不是为0或者最小值），就开始报错。而mas5因为使用的是原始值，所以不报错。

pengchy

或者先用倍数的方法筛选，然后再t.test

fengl

[quote]引用第3楼pengchy于2009-08-13 10:14发表的回 2楼(fengl) 的帖子 :

或者先用倍数的方法筛选，然后再t.test

[/quote]

Hi,pengchy,你有兴趣帮我分析下这几个数据的差异表达基因吗。我用R得出的结果居然和作者的不一样。PA (GSM49953.CEL.gz and GSM49954.CEL.gz) and LM2(GSM49955.CEL.gz, GSM49956.CEL.gz and GSM49957.CEL.gz).paper "Genes that mediate breast cancer metastasis to lung",supplementary method 2 in this paper.

我的代码如下

library(affy)<br />
LM2.PA <- ReadAffy()<br />
LM2.PA.mas5 <- mas5(LM2.PA)<br />
eSet <- exprs(LM2.PA.mas5)<br />
LM2.PA.p.value.all.genes = apply(eSet,1,function(x){t.test(x[1:2],x[3:5]) $p.value})<br />
LM2.PA.sub1 = subset(LM2.PA.p.value.all.genes,LM2.PA.p.value.all.genes<0.05)<br />
length(LM2.PA.sub1)<br />
names.sub1 <- names(LM2.PA.sub1)<br />
eSet.sub1 = eSet[names.sub1,]<br />
eSet.sub2 = subset(eSet.sub1,eSet.sub1[,1]>200 | eSet.sub1[,2]>200 | eSet.sub1[,3]>200 | eSet.sub1[,4]>200 | eSet.sub1[,5]>200)<br />
length(row.names(eSet.sub2))<br />
PA.raw.mean=apply(eSet.sub2[,c(1:2)],1,mean)<br />
LM2.raw.mean=apply(eSet.sub2[,c(3:5)],1,mean)<br />
eSet.sub3=rbind(subset(eSet.sub2,(PA.raw.mean/LM2.raw.mean)>3),subset(eSet.sub2,(LM2.raw.mean/PA.raw.mean)>3))<br />
length(row.names(eSet.sub3))

luzifer

做 differently expressed gene，有好多包可以参考，如limma， siggene等，一般不直接使用t-test，而是其它的modified t-test方法，可以参考http://bioconductor.org/packages/2.5/bioc/

pengchy

关于找差异基因，如luzifer战友建议，你可以看看siggene包。如果你是严格按照作者的方案做的，理论上差别不应该太大。你再仔细看看作者的思路吧。

另外，你找差异基因，还可以区分一下上调和下调吧。p值选择0.05是不是太高了？

fengl

最开始，我是想用作者的数据学习R来着。但后来发现结果不对，P值考虑FDR肯定是太高了，但我的重点不在这里，我现在只是想知道如果按照相同的条件，理论上应该是非常接近的。但错误出在哪里呢？

fengl

[quote]引用第6楼pengchy于2009-08-14 16:01发表的回 4楼(fengl) 的帖子 :

关于找差异基因，如luzifer战友建议，你可以看看siggene包。如果你是严格按照作者的方案做的，理论上差别不应该太大。你再仔细看看作者的思路吧。

另外，你找差异基因，还可以区分一下上调和下调吧。p值选择0.05是不是太高了？

[/quote]

嗯。我试下modified 的t test包的分析结果与作者的有什么不同。