差异分析的统计学本质探究

  提到表达量数据分析，不管是利用芯片技术还是高通量测序技术得到的表达量矩阵，我们都需要根据样本的分组信息来对所测得的基因或者蛋白分子来做差异分析，从而找到显著性变化的生物大分子。我们通常会用封装好的包来替我们完成这一过程，比如limma, samr, deseq, edgeR等等，它们的本质就是对表达量矩阵做一个归一化，然后利用理想的统计分布检验函数来计算差异的显著性。

所以差异分析的统计学本质，就是如何做归一化，以及使用什么样的统计分布！下面我们来看看Jimmy的分享！

Excel表格做差异分析

 其实主要讲的不是用Excel来做差异分析，只是想讲清楚差异分析的原理，用excel可视化的操作可能会更方便理解。而且想告诉大家，其实生物信息学分析，本来就很简单的，那么多软件，只有你理解了原理，你自己就能写出来的。

首先，还是得到表达矩阵，下面绿色的样本是NASH组，蓝色的样本是normal组。

我们进行差异分析，很简单，就是看两组的表达值，是否差异，而检验的方法就是T检验。

=AVERAGE(D2:L2)    ##求NASH组的平均表达量

=AVERAGE(M2:S2)    ###求normal的平均表达量

=T2-U2     ##计算得到logFOLDchange值

=AVERAGE(D2:S2)    ###得到所有样本的平均表达量

=T.TEST(D2:L2,M2:T2,2,3)    ###用T检验得到两个组的表达量的差异显著程度。

简单检查几个值就可以看到跟limma包得到的结果差不多。

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

下载该R语言包，然后看说明书，需要自己做好三个数据（表达矩阵，分组矩阵，差异比较矩阵），总共三个步骤（lmFit,eBayes,topTable）就可以了。

首先做第一个数据，基因表达矩阵。自己在NCBI里面可以查到下载地址，然后用R语言读取即可。

exprSet=read.table("GSE63067_series_matrix.txt.gz",comment.char = "!",stringsAsFactors=F,header=T)

rownames(exprSet)=exprSet[,1]

exprSet=exprSet[,-1]

然后做好分组矩阵，如下

然后做好，差异比较矩阵，就是说明你想把那些组拿起来做差异分析，如下

最后输出结果即可！

关于limma包差异分析结果的logFC解

首先，我们要明白，limma接受的输入参数就是一个表达矩阵，而且是log后的表达矩阵（以2为底）。那么最后计算得到的logFC这一列的值，其实就是输入的表达矩阵中case一组的平均表达量减去control一组的平均表达量的值，那么就会有正负之分，代表了case相当于control组来说，该基因是上调还是下调。

我之前总是有疑问，明明是case一组的平均表达量和control一组的平均表达量差值呀，跟log foldchange没有什么关系。后来，我终于想通了，因为我们输入的是log后的表达矩阵，那么case一组的平均表达量和control一组的平均表达量都是log了的，那么它们的差值其实就是log的foldchange首先，我们要理解fold change的意义，如果case是平均表达量是8，control是2，那么foldchange就是4，logFC就是2咯那么在limma包里面，输入的时候case的平均表达量被log后是3，control是1，那么差值是2，就是说logFC就是2。这不是巧合，只是一个很简单的数学公式log(x/y)=log(x)-log(y)。

用R语言的DESeq2包来对RNA-seq数据

做差异分析

library(DESeq2)

library(limma)

library(pasilla)

data(pasillaGenes)

exprSet=counts(pasillaGenes)  ##做好表达矩阵

group_list=pasillaGenes$condition##做好分组因子即可

(colData <- data.frame(row.names=colnames(exprSet), group_list=group_list))

dds <- DESeqDataSetFromMatrix(countData = exprSet,

colData = colData,

design = ~ group_list)

##上面是第一步第一步，构建dds这个对象，需要一个表达矩阵和分组矩阵！！！

dds2 <- DESeq(dds)  ##第二步，直接用DESeq函数即可

resultsNames(dds2)

res <-  results(dds2,contrast=c("group_list","treated","untreated"))

## 提取你想要的差异分析结果，我们这里是treated组对untreated组进行比较

resOrdered <- res[order(res$padj),]

resOrdered=as.data.frame(resOrdered)

可以看到程序非常好用！

它只对RNA-seq的基因的reads的counts数进行分析，请不要用RPKM等经过了normlization的表达矩阵来分析。

值得一提的是DESeq2软件独有的normlization方法！

rld <- rlogTransformation(dds2)  ## 得到经过DESeq2软件normlization的表达矩阵！

exprSet_new=assay(rld)

par(cex = 0.7)

n.sample=ncol(exprSet)

if(n.sample>40) par(cex = 0.5)

cols <- rainbow(n.sample*1.2)

par(mfrow=c(2,2))

boxplot(exprSet, col = cols,main=”expression value”,las=2)

boxplot(exprSet_new, col = cols,main=”expression value”,las=2)

hist(exprSet)

hist(exprSet_new)

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

 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)

看图好像没什么区别

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

* 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.Valuenames(pv_limma)=row
names(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。