在微阵列实验中，一个常见而重要的任务是识别其表达值在不同组或不同条件下有很大差异的基因。寻找这种差异表达的基因需要能够处理多个测试问题的方法，在这些问题中，数千甚至数万个假设被同时测试。

通常，适合于测试表达水平是否与兴趣的协变量相关的统计数据和相应的p值是按每个基因计算的。之后，针对多重性来调整原始p值，使得I类错误率被严格控制在预先指定的显著性水平。这种错误率的典型例子是家族错误率( family-wise error rate, FWER)，即至少一个假阳性的概率。然而，这个错误率对于检验数千个假说和识别几十个基因的情况来说可能太保守了。因此，在微阵列数据的分析中，另一个错误率变得非常普遍：错误发现率(FDR)，它粗略地说是所有被拒绝的空假设（即已识别的基因）中假阳性的预期比例。

然而，还有其他方法可以调整多重性：例如，可以使用QQ图或贝叶斯框架来实现这一目的。如果将观察到的测试统计量与在零假设下预期的测试统计量的值作图，则大多数点将近似位于对角线上。那些与这条线有很大不同的点对应于最有可能差异表达的基因。Tusher等人提出的微阵列显著性分析(SAM)可用于指定“本质上不同”的含义。Tusher等人(2001)的理论基于适度的t统计量，而Schwender等人(2003)将该方法与基于Wilcoxon秩和的SAM版本进行了比较。

Efron等人(2001)使用经验贝叶斯分析(EBAM)将观察到的测试统计量的分布建模为两个分量的混合，一个用于差异表达的基因，另一个用于非差异表达的基因。根据他们的分析，如果相应的后验概率大于0.9，则称一个基因为差异表达。

SAM和EBAM都在BioConductor软件包siggenes中实现。然而，在本文中，我们将集中讨论SAM。在下面，我们简要描述SAM过程、其在siggenes中的实现(有关更多细节，请参见该软件包中已有的测试统计数据。然后，我们将展示如何为其他测试情况编写您自己的函数。最后，我们将给出一个如何将SAM应用于基因表达数据的例子。

微阵列的显著性分析

在SAM中，针对m个基因中的每一个，计算适合于测试在表达水平和感兴趣的协变量之间是否存在关联的统计量d。这些观察到的测试分数被排序，并与零假设下的预期分数相对照，其中预期测试分数d(i)，i=1，.。。。，m计算如下：如果零分布已知，则d(i)是该零分布的(i−0.5)/m分位数。否则，通过以下方式评估d(i)：

-生成组标签的B个排列；
-计算m个测试统计量并针对B个排列中的每一个对它们进行排序；
-在B的第i个最小分数上求平均。

然后，将距离为∆的两条平行于对角线的直线绘制到这个称为SAM的图中。任何具有d值的基因

-大于或等于基因的d值，称为dup，对应于位于上∆线上方的原点右侧的最左侧点
-小于或等于基因的d值，称为dlow，对应于位于较低∆线以下的原点左侧的最右侧点。

被称为差异表达。然后，FDR通过下面的步骤来估算：

-通过计算有多少mB排列测试分数大于或等于dup或小于或等于dlow，并将该数字除以B，
-将该平均值除以识别的基因的数量，
-将该比率乘以基因没有差异表达的先验概率（默认情况下，sam通过Storey和Tibshiani，2003的过程估计该概率）。

对于∆的几个值重复该过程，并且选择在识别的基因数量和估计的∆之间提供最佳平衡的FDR值。

可以在sam中通过设置参数‘method’来调用不同的统计方法：

d.stat：适度t和F统计信息。如果模糊因子‘s0’被设置为零，则计算“常规的”t或F统计量。(将模糊因子添加到统计的分母中，以防止表达水平非常低的基因变得差异表达。详情见Tusher等人，2001年)。
wilc.stat：单类分析和两类分析的Wilcoxon秩和。
cat.stat：皮尔逊的SNP(单核苷酸多态性)数据等分类数据的χ2统计量(Schwender，2005年)。

示例：ALL数据集

我们采用ALL数据集来演示该包的用法。

<code>library(ALL)
data(ALL)/<code>

尽管在SAM分析之前过滤基因/探针集通常不是一个好主意，因为SAM假设大多数基因没有差异表达，但在这里，我们遵循Scholten和von Heydebreck(2005)的代码，过滤基因并选择样本的子集。

<code>library(genefilter)
subALL /<code>

(filterALL代码见附录)。它最终获得一个包含2391个探针集和79个样本的基因表达数据的exprSet对象。下面我们希望确定表达值在不同分子生物学种类（mol.biol）样本之间("BCR/ABL"和"NEG")有很大差异的探针集合。

<code>mol.biol /<code>

因此，通过指定所需的自变量‘data’和‘cl’来将sam应用于该数据集，其中‘data’可以是矩阵、数据框或包含基因表达数据的exprSet对象，‘cl’是包含样本的分类标签的向量。如果‘data’是一个exprSet对象，那么‘cl’也可以是命名包含类标签的pData(数据)列的字符串。

<code>library(siggenes)
clALL dataALL out1 /<code>

与下面代码等价：

<code>out2             var.equal = TRUE, rand = 123456)/<code>

其中‘var.equal’设置为TRUE，因为我们在这里假设组间方差相等，而‘rand’设置为123456，以使此分析的结果可重现。默认情况下，识别的基因数量和估计的FDR是针对10个等间距分布在0.1和maxi|d(i)−‘d(i)|之间的值来计算的。因此，我们的SAM分析的输出如下所示

其中，P0是估计的基因没有差异表达的先验概率，False是被错误调用的基因的数量，Called是识别的基因的数量，并且FDR=P0*False/Called是估计的FDR。可以使用summary获得更多信息，例如，模糊因子的值。∆的其他值也可以使用summary和print生成上表。例如,

<code>print(out1, seq(1.4, 2, 0.1))/<code>

假设我们选择的∆是1.5%。可以如下方式生成相应的SAM图。

<code>plot(out1, 1.5, sig.col = c(3,2), pch = 16,
     pos.stats = 2, cex = 0.6)/<code>

SAM Plot for ∆ = 1.5

‘sig.col’是指定已识别的下调和上调基因的颜色的数值或矢量，‘pos.stats’表示统计数据在SAM曲线图中显示的位置，‘

cex’指定未识别为差异表达的基因的绘图符号的相对大小。有关SAM特定方法Plot的所有参数，请参见

<code>help.sam(plot)/<code>

有关识别的基因的信息，例如它们的d值、相应的原始p值和q值(参见Storey和Tibshiani，2003)，可以通过以下方式获得

<code>summary(out1,1.5)/<code>

An excerpt from the output of summary(out1, 1.5)

上面信息还可以通过sam2excel存储在CSV文件中，也可以存储在使用sam2html的html文件中。如果‘data’是exprSet对象或在sam2html中指定了‘chipname’，则html文件还将包含基因符号以及指向公共存储库如Entrez、RefSeq和UniGene的链接。如果指定了“cdfname”，则指向已标识探针集合的Affymetrix网页的链接。例如，下面命令生成html格式的报告文件。

<code>sam2html(out1, 1.5, "out1.html", cdfname = "HG-U95Av2")/<code>

sam2html

总结

软件包siggenes包含用于执行微阵列的显著性分析(SAM)和微阵列的经验贝叶斯分析(EBAM)的功能。函数sam不仅为标准检验(如t检验和F检验)提供了一组统计信息，还提供了将用户编写的函数用于其他测试情况的可能性。在确定了一系列基因后，不仅可以获得这些基因的统计数据，如它们的测试分数和p值，还可以链接到包含这些基因的生物信息的公共存储库。EBAM功能目前正在修订中，以提供具有SAM功能已经具有的所有特征的这些功能的对用户更友好和消耗更少内存的版本。

附录

<code>filterALL pdat subsetas.character(pdat$BT)),
which(pdat$mol %in% c("BCR/ABL",
"NEG")))
eset require(genefilter)
f1 f2  0.5
selected filterfun(f1, f2))
esetSub pdat esetSub$mol.biol as.character(esetSub$mol.biol)
esetSub
}
/<code>

分享到:

閱讀更多 濤哥的醫學與生信之旅 的文章

關鍵字: 统计微阵列 FDR