本来搞差异分析的工具和包就一大堆了，而且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 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)
看图好像没什么区别
另外几个函数，我就不一一介绍了，大家可以自行探索。
* samr.plot(samr.obj, del, min.foldchange=0)
* samr.plot(samr.obj, del=.3)
* samr.assess.samplesize.obj
* 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对象的。
原文来自：http://www.bio-info-trainee.com/1608.html