用samr包对芯片数据做差异分析

Posted on by

本来搞差异分析的工具和包就一大堆了，而且limma那个包已经非常完善了，我是不准备再讲这个的，正好有个同学问了一下这个包，我就随手测试了一下，顺便看看它跟limma有什么差异没有！手痒了就记录了测试流程！

学习一个包其实非常简单，就是找到包的官网看看说明书即可！说明书链接

 

samr这个包更简单，就一个函数SAM,但是根据分析数据的不同被包装成了两个函数，分别是处理高通量测序数据的SAMseq和处理芯片数据的samr,本次我只讲解芯片数据的处理，然后跟limma这个包做一个简单比较~

所以，我们只需要制作好数据，然后学会用samr这个函数即可！

我们还是利用CLL这个包的测试数据来讲解这个包的用法，首先也是制作表达矩阵和分组信息。

suppressPackageStartupMessages(library(CLL))
data(sCLLex)
exprSet=exprs(sCLLex)   ##sCLLex是依赖于CLL这个package的一个对象
samples=sampleNames(sCLLex)
pdata=pData(sCLLex)
group_list=as.character(pdata[,2])
group_list
##  [1] "progres." "stable"   "progres." "progres." "progres." "progres."
##  [7] "stable"   "stable"   "progres." "stable"   "progres." "stable"  
## [13] "progres." "stable"   "stable"   "progres." "progres." "progres."
## [19] "progres." "progres." "progres." "stable"
as.numeric(as.factor(group_list))
##  [1] 1 2 1 1 1 1 2 2 1 2 1 2 1 2 2 1 1 1 1 1 1 2

这个表达矩阵exprSet和分组信息group_list就可以直接用来做差异分析啦~！ 它的分组信息要求比较读取，需要1,1,1,2,2,2这样的向量，所以我用了as.numeric(as.factor(group_list))，具体见下面的代码！

suppressPackageStartupMessages(library(samr))
data=list(x=exprSet,y=as.numeric(as.factor(group_list)), 
          geneid=as.character(1:nrow(exprSet)),
          genenames=rownames(exprSet), 
          logged2=TRUE
)
samr.obj<-samr(data, resp.type="Two class unpaired", nperms=100)

这样其实已经OK啦，重点是如何调整这个函数的参数，以及如何理解这个函数返回的结果(samr.obj这个对象非常重要，关乎你能否真正用好samr)~

我这里的genenames其实是探针名，如果真正要做分析，可以修改，而且我的nperms次数为100，也可以修改，一般是1000.

除了直接应用它找差异基因外，它还有几个单独的函数

首先是对表达矩阵进行normalization

x.norm <- samr.norm.data(data$x)
par(mfrow=c(1,2))
boxplot(exprSet, col = rainbow(exprSet),main="before normalization",las=2)
boxplot(x.norm,  col = rainbow(exprSet),main="after normalization",las=2)
QQ截图20160505194154

 

看图好像没什么区别

另外几个函数，我就不一一介绍了，大家可以自行探索。

* samr.plot(samr.obj, del, min.foldchange=0)

* samr.plot(samr.obj, del=.3)

* samr.assess.samplesize.obj<- samr.assess.samplesize(samr.obj, data, log2(1.5))

* samr.assess.samplesize.plot(samr.assess.samplesize.obj)

我们重点看看这个samr得到的差异与limma的差异区别在哪里

## 首先提取samr做差异分析检验的p值
pv=samr.pvalues.from.perms(samr.obj$tt, samr.obj$ttstar)
## 然后提取limma包做差异分析检验的p值
library(limma) 
design=model.matrix(~factor(sCLLex$Disease))
fit=lmFit(sCLLex,design)
fit=eBayes(fit)
options(digits = 4)
DEG_limma=topTable(fit,coef=2,adjust='BH',n=Inf) 
pv_limma=DEG_limma$P.Value
names(pv_limma)=rownames(DEG_limma)
head(pv[sort(names(pv))])
##  100_g_at   1000_at   1001_at 1002_f_at 1003_s_at   1004_at 
##    0.2531    0.4144    0.5671    0.5686    0.4687    0.6340
head(pv_limma[sort(names(pv_limma))])
##  100_g_at   1000_at   1001_at 1002_f_at 1003_s_at   1004_at 
##    0.2497    0.4312    0.5349    0.5498    0.4361    0.6473
cor(pv[sort(names(pv))],pv_limma[sort(names(pv_limma))])
## [1] 0.9976

从数据上来看，没什么本质区别,而且相关系数高达0.9978.

所以结论是，没必要搞那么多的包，用limma就好了，甚至直接用t检验也是OK的

还有plot和summary也是可以直接作用于samr的结果samr.obj对象的

Comments are closed.