R语言实现基因表达模式识别

先验概率:是指根据以往经验和分析得到的概率{P(事件)}.

后验概率:事情已经发生，要求这件事情发生的原因是由某个因素引起的可能性的大小{P（事件|原因）}。

基础的知识有了那么我们看下贝叶斯的理论公式：

对待模型的自我心得就是所有公式都是纸老虎，捅破了其实就那么回事。就怕你不敢去捅。接下来我们就捅一捅，看看在我们的这个基因表达数据中是如何应用的。

首先是提出假设：所有样本之间的表达是可以互换的，那么就可以通过构建每个基因在不同样本表达偏差（相对于每个基因均值大小）进行后期的表达分布模式函数fobs (·|µj )的构建。

接下来，有两个情况：

一个是单条件的两组样本的情况，也就形成在两组样本之间存在差异（DE）和没有差异的基因(EE)。针对这个情况我们得到我们所有样本的整体分布函数：

其中参数π(u)指的系统特定的分布模型。也就是边缘分布函数。

那么，以上的情况是只能对EE基因进行评估，如果是DE的基因需要各自计算均值，那么就得到如下的公式：

那既然既存在差异基因也存在无差异基因，那么我们需要引入一个概率参数P（差异基因），无差异基因就是1- P（差异基因）。至此，我们前面介绍的贝叶斯理论就开始派上用场了：

第二种情况就是多模式混合的情况,计算原理比前面的更加复杂，那就是增加了模式的概率P的参数。

最后，我们需要把上面的这个函数思想应用到指定的分布函数模型中，作者便引入了下面的三个模型：

最后，我们就来看下作者开发的包里面具体函数应用，主要的函数功能如下（就不翻译了）：

包的安装还是通过bioconductor来安装：

BiocManager::install("EBarrays")

首先看下基因表达模式创建函数ebPatterns。具体参数不介绍了，需要注意的是如果设置为0那么就是不列入分析的样本。看下实例：

library(EBarrays)
patterns <- ebPatterns(c("1 1 1 1 11 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1",
                         "1 1 1 1 1 1 1 11 1 1 1 1 2 2 2 2 2 2 2 2 2 2
                          2 2 2"), TRUE)
show(patterns)

上面就得到了我们的样本分布模式，接下来就是基于以上的模式对表达数据进行模型的构建，利用EM算法构建了emfit函数：

其中主要的参数：

Family 设置分布模型，可以选择作者提出的几个模型，也可以自己设置。

Hypotheses 样本的分布模式，多组或者单组。

接下来我们看下实例：

data(sample.ExpressionSet) ## from Biobase
eset <- exprs(sample.ExpressionSet)
patterns <- ebPatterns(c("1 1 1 1 11 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1",
                         "1 1 1 1 1 1 1 11 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2"))
gg.fit <- emfit(data = eset, family ="GG", hypotheses = patterns, verbose = TRUE)
show(gg.fit)

接下来我们可以输出后验概率：

gg.post.out <- postprob(gg.fit,eset)$pattern

然后就是FDR计算所需要的阈值分计算，需要用到函数crit.fun。直接我们看下实例：

gg.threshold <-crit.fun(gg.post.out[,1], 0.05)

这个阈值意义就是某一种表达模式中那些具有此模式表的基因的后验概率必须大于我们的阈值。这也就是我们最终想要的符合某种表达模式基因列表。