本来搞差异分析的工具和包就一大堆了，而且limma那个包已经非常完善了，我是不准备再讲这个的，正好有个同学问了一下这个包，我就随手测试了一下，顺便看看它跟limma有什么差异没有！手痒了就记录了测试流程！
    学习一个包其实非常简单，就是找到包的官网看看说明书即可！说明书链接
    samr这个包更简单，就一个函数SAM,但是根据分析数据的不同被包装成了两个函数，分别是处理高通量测序数据的SAMseq和处理芯片数据的samr,本次我只讲解芯片数据的处理，然后跟limma这个包做一个简单比较~
    所以，我们只需要制作好数据，然后学会用samr这个函数即可！
    我们还是利用CLL这个包的测试数据来讲解这个包的用法，首先也是制作表达矩阵和分组信息。

    1. suppressPackageStartupMessages(library(CLL))
    3. data(sCLLex)
    4. exprSet=exprs(sCLLex)   ##sCLLex是依赖于CLL这个package的一个对象
    6. samples=sampleNames(sCLLex)
    7. pdata=pData(sCLLex)
    9. group_list=as.character(pdata[,2])
    10. group_list
    11. ##  [1] "progres." "stable"   "progres." "progres." "progres." "progres."
    12. ##  [7] "stable"   "stable"   "progres." "stable"   "progres." "stable"  
    13. ## [13] "progres." "stable"   "stable"   "progres." "progres." "progres."
    14. ## [19] "progres." "progres." "progres." "stable"
    15. as.numeric(as.factor(group_list))
    16. ##  [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))，具体见下面的代码！

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

    这样其实已经OK啦，重点是如何调整这个函数的参数，以及如何理解这个函数返回的结果(samr.obj这个对象非常重要，关乎你能否真正用好samr)~
    我这里的genenames其实是探针名，如果真正要做分析，可以修改，而且我的nperms次数为100，也可以修改，一般是1000.

    除了直接应用它找差异基因外，它还有几个单独的函数
    首先是对表达矩阵进行normalization

    1. x.norm 
    2. par(mfrow=c(1,2))
    3. boxplot(exprSet, col = rainbow(exprSet),main="before normalization",las=2)
    4. boxplot(x.norm,  col = rainbow(exprSet),main="after normalization",las=2)

    image.png
    看图好像没什么区别
    另外几个函数，我就不一一介绍了，大家可以自行探索。
    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的差异区别在哪里

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

    从数据上来看，没什么本质区别,而且相关系数高达0.9978.
    所以结论是，没必要搞那么多的包，用limma就好了，甚至直接用t检验也是OK的
    还有plot和summary也是可以直接作用于samr的结果samr.obj对象的。

    原文来自：http://www.bio-info-trainee.com/1608.html