A new unsupervised gene clustering algorithm based on the integration of biological knowledge into expression data

Marie Verbanck; Sébastien Lê; Jérôme Pagès

doi:10.1186/1471-2105-14-42

BMC生物信息学。2013年；14时42分。

2013年2月7日在线发布。数字对象标识：10.1186/1471-2105-14-42

预防性维修识别码：PMC3635920型

PMID：23387364

一种新的基于生物知识与表达数据集成的无监督基因聚类算法

玛丽·韦班克,¹ 塞巴斯蒂安·莱伊,¹和杰罗姆·帕格斯¹

作者信息文章注释版权和许可信息 PMC免责声明

摘要

背景

生物学家在分析组学数据时大量使用基因聚类算法。经典的基因聚类策略仅基于表达数据的使用，直接如在热图中，或间接如在基于共表达网络的聚类中。然而，经典策略可能不足以揭示基因之间的所有潜在关系。

结果

我们提出了一种新的无监督基因聚类算法，该算法将外部生物知识（如基因本体注释）集成到表达数据中。我们引入了一种新的基因间距离，即将生物学知识整合到表达数据分析中。因此，如果两个基因同时具有相似的表达谱和相似的功能谱，则它们是相近的。然后使用经典算法（如K-means）获得基因簇。此外，我们提出了一种基因簇的自动评估程序。该程序基于两个指标，这两个指标衡量基因簇的全球共表达和生物同质性。它们与假设检验相关联，假设检验允许用p值来补充每个指标。

在模拟数据和实际数据上，将我们的聚类算法与热图聚类和基于基因共表达网络的聚类进行了比较。在这两种情况下，它都优于其他方法，因为它提供了最高比例的显著共表达和生物同源的基因簇，这是很好的解释候选者。

结论

我们的新聚类算法为解释提供了更高比例的好候选。因此，我们期望对这些簇的解释能帮助生物学家对基因之间的关系形成新的假设。

背景

由于转录组分析数据等组学数据提供了对大量基因的测量，因此通过将基因聚类到模块中，可以将数据经典地分解到更容易理解的水平。在无监督的聚类策略中，我们可以回忆起主要使用的两种技术：热图[1]这包括主题和基因表达的层次分类，以及基于共表达网络的聚类[2]. 基因聚类不仅实用，因为它减少了研究对象的数量，而且有望传递某种生物现实。事实上，我们期望基因表达之间的相似性反映基因功能之间的相似。然后对基因簇进行解释，以产生关于基因功能作用及其关系的新假设。

在实践中，为了解释基因簇，外部生物学知识，如基因本体论（GO）信息[三]使用。最经典的程序包括基因集富集分析，目的是通过一组生物功能表征每个簇。已经提出了改进基因集富集分析的尝试，例如Bauer等人[4]提出了贝叶斯富集分析。后者包括将GO项表示为贝叶斯网络，每个基因的响应在表达方面被建模为GO项激活的函数。在多元分析（MVA）中，存在一些将生物学知识直接叠加在MVA输出上的尝试[5,6]. 其目的是促进基因表达或基因簇的解释，因为MVA提供了可用于聚类的距离矩阵。

在这些方法中，基因簇仅基于表达数据获得，生物知识是用来充分利用簇的后验知识。这种方法的局限性是明确的：基于表达数据的基因聚类只允许分离共表达单元，但不一定是生物学上的连贯单元[7,8]. 事实上，聚类结构只能与它所基于的距离/相似性矩阵一样好。因此，积极将生物知识集成到表达式数据中，以分离更有意义的生物实体的想法。

在其他情况下，积极将生物知识整合到表达数据中的这一问题也得到了解决。为了进行生物网络推断，Kashima等人[9]提出了一种半监督学习方法。如果在近亲物种中观察到相同的相似性，则特定物种中的表达谱和氨基酸序列之间的相似性会得到加强。为了预测基因功能类别，例如基因和GO术语之间的关联，Azuaje等人[10]结合两种类型的信息：基因表达谱相似性和基于GO的相似性。两个相似性指数的平均值用于基因聚类。与预测基因功能类别的目标相同，Li等人[11]，通过考虑与同一功能注释相关的基因子集，将表达数据与生物学知识相结合。然后根据基因表达谱的相似性对基因子集进行聚类。

本文的目标是提出一种新的无监督聚类算法，该算法基于基因之间的新距离，主动将外部生物知识整合到表达数据中。根据生物学知识，如果一个簇收集了与类似生物功能相关的共表达基因，则认为该簇是令人满意的。这样的星团有望在生物学上引起人们的兴趣，并成为生物学解释的良好候选者。

在实践中，我们引入了共表达生物功能的概念，它允许在功能注释中集成共表达信息。将表达数据与GO注释相结合，定义了基因之间的新距离。如果两个基因同时被表达并牵涉到同一组生物功能中，那么它们是相近的。然后使用经典的聚类算法（K-means或分层升序分类）来获得基因聚类。在本文中，我们将强调支持该方法的生物学原理，并讨论我们提出的距离。

为了补充聚类程序，我们提出了一种基因聚类的自动验证程序，以便于对其进行解释。该程序的目的是突出解释的良好候选者，这些候选者是显著共表达和显著生物学相关的基因簇。它基于与假设检验相关的两个指标。一个指标衡量集群内基因的共表达，而另一个指标量化其生物同质性。

用于执行所有分析的R代码以R包的形式提供，网址为http://marie.veranck.free.fr/packages/.

方法

将生物学知识集成到表达数据中：生物学原理

让我们回顾一下，大多数经典的基因聚类策略仅基于表达数据。表达式数据可以直接用于热图，也可以间接用于基于共表达网络的聚类。由此获得的簇是用于解释的候选者，并且仍有待生物学表征。生物特征化是使用外部生物知识完成的，例如基因本体注释。这些是根据文献中报告的实验或生物信息学推导得出的。这一经典方法依赖于两个隐含假设。首先，共表达簇的生物学特性表明，共表达基因之间存在系统的生物学联系。其次，生物表征纯粹基于外部生物知识，因此，部分外部生物知识预计与研究中的实验相关。

第一个假设可能有问题[7,8]在本文中，我们考虑了一个新的观点，即共存与生物联系之间的联系。广义地说，两个基因之间的共表达可能由两种现象引起，要么是真正的生物联系（例如来自真正的基因调控网络），要么是对同一实验条件的不同生物反应的平行和独立激活。为了区分这两种情况，我们建议更多地相信第二种假设，然后积极依赖外部生物知识。因此，我们认为，如果两个共表达基因在现有生物学知识中已经被描述为生物学相关，那么它们的共表达更有可能反映出真正的生物学联系。

在实践中，我们使用与GO注释的“生物过程”相关的本体，该本体为每个基因提供了该基因所涉及的生物功能列表：此后，该列表将被称为功能配置文件基因。因此，如果两个共表达基因具有相似的功能特征，则推测它们的共表达是由真正的生物学联系引起的。相反，如果两个共表达基因具有完全不同的功能谱，则它们的共表达可能是由于不同生物反应的并行激活所致。

无监督基因聚类算法

在本节中，我们提出了一个新的基因间距离，它符合暴露的生物学原理，并将用于聚类的角度。这个距离可以量化两个基因之间的共表达和功能相似性。

生物知识的编码

让我们考虑一下K（K）基因和J型GO条款。基因和GO注释之间的关联编码在二进制矩阵中T型 ∈ ℳ(K（K）, J型)，其中每行k个表示其中一个K（K）基因和每一列j个其中一个J型GO术语：通用术语T型_k个j个等于1，如果基因k个与GO项关联j个否则为0（图（图1）。1). 一排k个基质的可解释为一个基因功能谱，它是与基因相关的一组生物功能。A列j个基质代表一种生物功能，可以被同化为与该功能相关的基因子集。让K（K）^j个={k个|T型_k个j个＝1}是与该功能相关的基因的子集j个.

在单独的窗口中打开

图1

矩阵T型：编码基因和生物功能之间的关联。基因和生物功能之间的关联在基质中合成T型每行代表一个基因功能概况，而每列代表生物功能和基因之间的关联。通用术语T型_k个j个如果基因等于1k个与生物功能有关j个，否则为0。行边距T型_k个.是基因的生物功能数量k个与关联。列边距T型_.j个是功能基因的数量j个与关联。最后，T型_..等于基因和生物功能之间关联的总数。

基因间的新距离：共表达生物功能

为了符合先前暴露的生物学原理，我们定义了一个距离来量化功能轮廓的相似性{T型_k个j个;j个∈J型}共表达基因。为此，我们通过定义一个共存生物功能作为功能对唯一共存基因的限制。换句话说，如果K（K）^j个可以分为L（左）_j个共存簇，这将导致尽可能多的共存生物功能被考虑。为了获得这些共表达的生物功能，我们提出了以下基于层次聚类的算法。

对于每个生物功能j个:

1.以下基因之间的距离矩阵K（K）^j个基于皮尔逊相关系数进行计算。两个基因之间的距离k个和k个^′可以表示为：

{d日}_{G公司} (k个, {k个}^{'}) = 1 - \frac{1}{我} \sum_{我 = 1}^{我} (\frac{{G公司}_{伊克} - {G公司}_{.k（千）}}{{S公司}_{k个}}) (\frac{{G公司}_{我 {k个}^{'}} - {G公司}_{. {k个}^{'}}}{{S公司}_{{k个}^{'}}})

(1)

哪里我是样本数，G公司_我k个和G公司_我k个^′分别是基因的表达k个和k个^′用于样品我,G公司_.k个和G公司_.k个^′分别是我基因表达值k个和k个^′,S公司_k个和S公司_k个^′分别是我基因表达值k个和k个^′.

2.对前面定义的距离矩阵（1）执行层次聚类过程：let ${P（P）}^{j个} = {{K（K）}_{1}^{j个}; ...; {K（K）}_{我}^{j个}; ...; {K（K）}_{{L（左）}_{j个}}^{j个}}$ 是上的分区K（K）^j个在里面L（左）_j个集群。对于所有人我= 1, …, L（左）_j个, ${K（K）}_{我}^{j个}$ 由共表达基因组成。

3.我们建立了一个矩阵T型^j个 ∈ ℳ(K（K）, L（左）_j个)通过拆分j个^t吨小时第列，共列T型进入之内L（左）_j个柱。在T型^j个每条线k个表示其中一个K（K）基因和每列都是一个虚拟变量，例如 ${T型}_{肯尼亚}^{j个}$ 如果基因等于1k个属于 ${K（K）}_{我}^{j个}$ 否则为0：一列T型^j个可以解释为一种共存的生物功能。

我们定义T型_{c（c）o（o）e（电子）x个第页}作为所有事物的并列J型矩阵T型^j个（图（图2）。2)。T型_{c（c）o（o）e（电子）x个第页}这两种信息的组合结果。分析T型_{c（c）o（o）e（电子）x个第页}允许研究在共表达条件下基因功能谱的相似程度。因此，可以通过以下公式计算基因之间的新距离T型_{c（c）o（o）e（电子）x个第页}:

在单独的窗口中打开

图2

矩阵T型_共存：矩阵的分解T型.将生物功能分解为共存的生物功能导致构建矩阵T型_{c（c）o（o）e（电子）x个第页}其中一行表示基因，一列表示共表达的生物功能。通用术语T型_{c（c）o（o）e（电子）x个第页}, ${T型}_{千焦} 1_{k个 \in {K（K）}_{我}^{j个}}$ 如果基因等于1k个与函数关联j个如果它属于集群 ${K（K）}_{我}^{j个}$ ，否则为0。共存生物功能的列边界我等于相应簇中的基因数，即 $卡片 ({K（K）}_{我}^{j个})$ 此外，对于每个功能j个，与衍生的共表达生物功能相关的列边距之和j个等于与函数关联的列边距j个: $\sum_{我 = 1}^{{L（左）}_{j个}} 卡片 ({K（K）}_{我}^{j个}) = {T型}_{.j英寸}$ 最后，我们可以注意到，行边距和关联总数等于T型.

{d日}_{{T型}_{共存}} (k个, {k个}^{'}) = \sum_{j个 = 1}^{J型} \sum_{我 = 1}^{{L（左）}_{j个}} \frac{{T型}_{..}}{卡片 ({K（K）}_{我}^{j个})} (\frac{{T型}_{千焦}}{{T型}_{k、。}} 1_{k个 \in {K（K）}_{我}^{j个}} - \frac{{T型}_{{k个}^{'} j个}}{{T型}_{{k个}^{'} .}} 1_{{k个}^{'} \in {K（K）}_{我}^{j个}})^{2}

(2)

哪里T型_k个.和T型_k个^′.分别是与基因相关的行边距k个和k个^′,T型_..是基因和生物功能之间的关联总数 $1_{k个 \in {K（K）}_{我}^{j个}}$ 等于1的虚拟变量，如果 $k个 \in {K（K）}_{我}^{j个}$ ，否则为0。基因k个和k个^′都与j个：如果它们不共存，则它们不属于同一共存簇P（P）^j个。在这种情况下j个^t吨小时距离计算项（2）较高。因此，具有相似表达谱和相似功能谱的基因是相近的。此距离对应于对应分析中的基因之间的距离T型_{c（c）o（o）e（电子）x个第页}.

技术说明1：在步骤2中，P（P）^j个分区在L（左）_j个与生物功能相关的共表达基因簇j。P（P）^j个通过切割分类树来确定。切割分类树提供了一个分区，并允许计算所讨论分区的簇内惯量之和。计算中隔板之间的相对惯性损失L（左）簇和中的分区L（左）+1群集为 $\frac{\sum_{我 = 1}^{L（左） + 1} 惯性 (我)}{\sum_{我 = 1}^{L（左）} 惯性 (我)}$ .P（P）^j个通过切割分类树获得相对惯性损失较高的分区。

技术说明2：在所有与j个，是共表达的，j个被认为是一种共存的生物功能。我们添加了一个步骤0。包括过滤生物功能：它允许定义生物功能j个可以认为是共存的。因此，与j个根据下一节介绍的程序，通过计算共表达指示剂的p值进行测试。如果该p值低于所选阈值（例如10%），则所讨论的函数被视为一个共存函数，不会在T型_共存，但它是保守的。

注：在完全不同的背景下，为了预测基因功能类别，Li等人[[11]]基于检测异质功能类中同质共表达基因亚群的思想，提出了一种模糊近聚类算法。这种检测使他们能够更好地预测基因功能类别。

获取基因簇

为了获得基因簇，将聚类算法（如K-means或层次升序分类）应用于距离矩阵。我们希望通过此过程获得共表达和生物相关基因的簇。

基因簇评估

为了使一个簇成为一个很好的解释候选，它必须收集共表达和生物相关的基因。经典评估程序侧重于所谓的生物同质性簇及其生物功能特征。然而，在我们的聚类过程中，共表达必然与生物同质性竞争，因为这两种类型的信息是积极结合的。因此，我们提出了一种基于两个指标的基因簇评估程序：共表达指标和与假设检验相关的生物同质性指标。

共表达指示器

共表达被定义为两个基因之间的正相关。事实上，如果两个基因正相关，那么它们在相同的实验条件下会过度表达或表达不足。我们希望找到一种能综合簇内相关性的共表达指标（CI）。我们考虑一个经验性但方便的指标，即同一簇基因之间相关性的平均值K（K）_我该指标计算如下：

\begin{array}{l} CI公司 ({K（K）}_{我}) = & \frac{1}{\frac{卡片 ({K（K）}_{我}) (卡片 ({K（K）}_{我}) - 1)}{2}} \sum_{k个 | k个 \in {K（K）}_{我}} \\ \times (\sum_{{k个}^{'} | {k个}^{'} \in {K（K）}_{我}, {k个}^{'} > k个} \frac{1}{我} \sum_{我 = 1}^{我} (\frac{{G公司}_{伊克} - {G公司}_{.k（千）}}{{S公司}_{k个}}) (\frac{{G公司}_{我 {k个}^{'}} - {G公司}_{. {k个}^{'}}}{{S公司}_{{k个}^{'}}})) \end{array}

(3)

哪里我是样本数，G公司_我k个和G公司_我k个^′分别是样本的表达式我基因的k个和k个^′,G公司_.k个和G公司_.k个^′分别是我基因的表达值k个和k个^′,S公司_k个和S公司_k个^′分别是我基因的表达值k个和k个^′.

共表达指标确实提供了一种衡量全球基因簇共表达情况的方法。其范围为−⅓ 至1（见附录附录1：共表达指标的变化范围）。如果所有基因完全共存，则指标等于1。相反，让我们考虑一个基因没有共表达的集群，以至于区分了两个子集群：在每个子集群内，基因是正相关的，在子集群之间，它们是负相关的。在这种情况下，指示器接近于0，并且可能小于0。

生物均匀性指示剂

我们的目标是根据基因功能谱的相似性定义一个生物同质性指标。经典地，基因簇的生物同质性是通过与其相关的丰富生物功能的数量和性质来评估的，通过富集试验确定簇的特征并不能保证功能谱的相似性，因为富集试验是针对每个生物功能分别进行的。Datta和Datta[12]提出了一种多维生物同质性指标，目的是评价整个聚类过程，而不是聚类本身。我们采用这个想法来测量基因簇的生物同质性。我们认为，作为生物同质性指标（BHI），一个系数是从克雷默公式得出的V（V）系数[13]它提供了来自K（K）_我该指标计算如下：

\begin{array}{lcr} BHI公司 ({K（K）}_{我}) & = & 1 - \sqrt{\frac{\sum_{k个 \in {K（K）}_{我}} (\sum_{j个 = 1}^{J型} \frac{{({T型}_{千焦} - \frac{{T型}_{k、。} {T型}_{.j英寸}}{{T型}_{..}})}^{2}}{\frac{{T型}_{k、。} {T型}_{.j英寸}}{{T型}_{..}}})}{{T型}_{..} (卡片 ({K（K）}_{我}) - 1)}} \end{array}

（4）

哪里T型_k个j个如果基因等于1k个与生物功能有关j个，否则为0，T型_k个.是与基因相关的行边距k个.

生物均匀性指标在0和1之间变化（见附录2：生物均匀性指数的变化范围）。因此，如果一个簇中的所有基因都具有完全相似的功能，那么生物同质性指标等于1。相反，如果没有一个基因具有相似的功能，以至于没有一个生物功能与来自K（K）_我则生物均匀性指标为0。

尽管这一指标有其局限性，因为生物的同质性应主要依赖于生物学解释，然而，自动评估基因簇的生物利益碰巧是有用的。

假设检验程序

我们用假设检验程序来补充这些指标，因为这两个指标在很大程度上取决于集群的规模，所以更为合理：

•共表达指示剂：在其计算中（3）除以 $\frac{卡片 ({K（K）}_{我}) (卡片 ({K（K）}_{我}) - 1)}{2}$ 执行时，CI的值会随着集群的大小而机械地减少

•生物同质性指标：按c（c）一第页d日(K（K）_我)−1在其计算的第二项（4）中执行，由于第二项在0和1之间变化，BHI的值会随着集群的大小而机械地增加

目的是评估一种方法在多大程度上提供了共存和生物同质性高于随机聚类的聚类。因此，随机聚类对应于检验的零假设，随机聚类的指标值被用作参考情况。实际上，将p值与集群关联K（K）_我对于一个指标，相同大小的簇是由简单地绘制基因而不进行替换构成的。然后计算每个簇的指标，从而获得零假设下指标值的分布。通常，与待测试聚类的指标值相对应的观测值被定位在零假设下的相应分布中。最后，通过指标值优于观测值的随机组成簇的比例来估计p值。

注1：该过程的兴趣在于在零假设下获得分布的方式。由于指标的计算仍然基于实际数据，因此在零假设下的分布与数据的分布有关。

注2：很明显，由单个基因组成的簇未进行测试。

结果

当我们提出一种新的无监督聚类算法并对聚类进行自动评估时，我们通过将其与两种最经典的基因聚类策略进行比较，在模拟和实际数据集上验证了整个方法。一方面，我们将其与来自表达式数据的热图的聚类进行比较。另一方面，我们选择使用加权基因共表达网络（WGCNA）从表达数据生成共表达网络[2]. 共表达网络允许根据网络节点的拓扑重叠计算基因之间的差异矩阵。最后，根据差异矩阵计算层次聚类算法，并提供基因聚类。

模拟研究

模拟数据集

在本节中，我们将解释如何模拟表达式和GO数据集。

为了模拟表达式数据，我们使用与中相同的过程[14]. 表达式数据矩阵G公司_秒我米，由K（K）基因和我样本是在具有特定相关结构的多元高斯分布中随机抽取的，因此我们有共同表达基因的潜在簇。由于这种模拟数值数据的方法相当经典，因此我们坚持模拟GO注释数据，这在文献中并不常见。

为了模拟GO注释数据，我们采用了之前公开的生物学原理：GO注释由与研究中的实验相关的信息和与实验无关的信息组成。换句话说，模拟GO注释的一部分必须具有与表达式数据结构相似的结构，另一部分必须有随机结构。因此，模拟GO矩阵T型_秒我米通过并列两种类型的矩阵获得：

•

{T型}_{模拟}^{e（电子）}

：其基因功能谱模拟基因表达谱，因此当两个基因在G公司_秒我米，它们在

{T型}_{模拟}^{e（电子）}

•

{T型}_{模拟}^{第页}

：其基因功能谱与基因表达谱无关

实际上，为了获得 ${T型}_{模拟}^{e（电子）}$ 首先，我们仅基于它们的表达谱之间的相关性构建一个基因分类树。然后我们考虑每个节点j个分类树的生物功能。如果基因k个与节点关联j个分类树的， ${T型}_{模拟}^{e（电子）} (k个, j个) = 1$ ，否则为0。因此，具有相似表达谱的基因在机械上共享紧密的功能谱。为了获得 ${T型}_{模拟}^{第页}$ ，我们并排第页乘以矩阵 ${T型}_{模拟}^{e（电子）}$ 并独立排列每列中的行，其中第页是一个整数，表示T型_秒我米：具体来说，有第页随机生物功能是结构化生物功能的数倍T型_秒我米.

这种生成类似矩阵的方法 ${T型}_{模拟}^{e（电子）}$ 选择它是因为它模拟了GO信息的层次结构。这种生成随机矩阵的方法 ${T型}_{模拟}^{第页}$ 允许保留生物功能的边缘，重要的是这些边缘代表了与功能相关的基因数量，可以解释为功能的特定程度。

结果

在实践中，我们将这三种方法应用于模拟数据集。我们考虑两种大小的模拟表达式数据。第一种类型由10个个体和300个基因组成，对于每种方法，我们将其划分为20个簇。第二种类型由25个个体和1000个基因组成，对于每种方法，我们将其划分为100个簇。对于这两种类型的模拟表达数据集，我们将随机性强度在1到3之间的模拟GO注释关联起来。对于每个配置，将生成100个数据集。

无论采用何种聚类方法，我们都会将每个聚类关联到两个p值，分别对应于共表达指标和生物同质性指标。对于给定的分区，我们测量以下簇的比例：

•显著共存：与CI相关的p值低于所选阈值

•生物均匀性显著：与BHI相关的p值低于所选阈值

•两者显著共存且生物均一：与CI和BHI相关的p值均低于选定阈值

结果汇总在表中表1。1平均而言，这三种方法都为分区提供了大量显著共存的簇。该比例并不取决于Heatmap和WGCNA的随机性强度。然而，对于我们的聚类算法，当随机性强度增加时，我们观察到显著共存的簇的比例略有下降。这是意料之中的，因为当随机性强度较高时，共表达与生物同质性的竞争更为激烈。

表1

模拟研究结果

			共表达指示器			生物均匀性指示剂			两者都有
我	*K（K）*	第页	热图	WGCNA公司	集成	热图	WGCNA公司	集成	热图	WGCNA公司	集成
10	300	1	92.15	94.90	98.65	65.50	81.5	89.5	64.60	78.95	88.80
10	300	2	92.31	94.80	96.55	50.40	60.15	67.25	49.75	58.30	66.25
10	300	三	92	95.32	94.52	36.77	45.81	54.03	36.61	45	53.39
25	1000	1	88.70	99.12	91.33	7.67	28	45.44	7.35	27.09	44.72
25	1000	2	90.25	99.12	90.55	3.79	11.89	29.62	3.54	11.17	28.95
25	1000	三	89	98.99	85.67	1.94	3.55	18.66	1.80	3.34	18.06

在单独的窗口中打开

三种聚类算法的仿真研究结果：热图分类（Heatmap classification）、基于共表达网络的聚类（WGCNA）和我们的聚类算法（Integration）。模拟数据集因样本数而异(我)，基因数量(K（K）)和随机性的强度(第页). 我们给出了给定分区中簇的平均比例（%），这些簇是显著共存（CI）、生物均匀（BHI）或同时共存和生物均匀（both）。让我们以包含10个个体和300个变量的模拟表达式数据集为例，与随机度为1的模拟GO注释相关联。这些数据集的热图平均为分区提供92.15%的显著共存簇。

平均而言，来自热图的分区具有较低比例的簇，这些簇在生物学上具有明显的均匀性。当随机性的强度增加时，这个比例会严重降低。考虑到基因表达背后的网络结构是有益的，因为它提供了比热图更大比例的显著生物同质簇。然而，当随机性强度非常高时，WGCNA提供的生物均匀簇的比例实际上会下降。我们的聚类算法提供了相当高比例的生物均质聚类，即使随机性强度等于3。

如果我们关注显著共存且生物同质的簇的比例，我们的聚类算法将优于其他两种方法。

鸡肉数据集分析

该方法应用于与已发布数据集相关的转录组数据集示例[15]. 通过这项实验，目的是了解鸡绝食后的遗传机制。因此，在27只处于4种营养状态的鸡中收集到约12000个肝脏基因的表达：16小时禁食“F16”，16小时禁餐后5小时营养期“F16R5”，16个小时禁食后16小时营养期（F16R16），最后是持续喂养状态“F”。在我们的例子中，我们选择了一些基因，这些基因的表达因实验因素而异，这使得我们通过多重测试的因子分析方法保留了大约3600个基因[16].

此外，类似于Busold等人[5]，我们使用GO信息，其中考虑了GO项之间的层次结构：当一个基因与一个项关联时，它会自动与其父项关联。

在模拟研究中，我们进行了三个与热图对应的基因聚类，一个基于共表达网络（WGCNA）的聚类和我们自己的聚类过程。我们选择将每个过程获得的簇数设置为200。对于给定的分区，我们将共表达指示剂和生物同质性指示剂的两个p值与每个簇关联，这两个值在联合图中可视化。在图中图3，三，一个点表示一个簇，其在x轴上的值等于共表达指示剂p值，在y轴上的数值等于生物同质性指示剂p的值。此外，表表22提供了三个分区中每个分区的簇的比例，如模拟研究中所示，这些分区是显著共存（CI）、生物均质（BHI）或同时共存和生物均质的（both）。

在单独的窗口中打开

图3

三种聚类程序中与共表达指标和生物均匀性指标相关的p值的表示。应用于鸡表达数据的三种聚类程序的结果如下：热图（Heatmap）、基于共表达网络的聚类（WGCNA）和我们的聚类算法（Integration）。无论采用何种聚类方法，对于每个聚类，都会关联一个与共表达指标对应的p值和一个与生物同质性指标对应的p-值。P值以联合表示法表示，其中每个点代表一个簇，与共表达指示剂相关的P值在x轴上表示，而与生物同质性指示剂相关联的P值则在y轴上表示。

表2

案例研究结果

	CI公司	BHI公司	两者都有
热图	91.50	13.50	13.50
WGCNA公司	63	68	46
集成	53.50	79.50	53.50

在单独的窗口中打开

针对鸡数据集的三种聚类算法的结果：热图分类（Heatmap）、基于共表达网络的聚类（WGCNA）和我们的聚类算法（Integration）。我们给出了在给定分区中，显著共表达（CI）、生物同质（BHI）或同时表达和生物同质的簇的百分比（%）。

首先，热图提供的分区由大多数显著共存的簇组成（91.50%）。然而，一小部分集群在生物学上具有明显的同质性，以至于与BHI相关的p值似乎是按照均匀分布分布的。QQ图（图（图4）4)实际上证实了与生物均匀性指标相关联的p值分布可以被认为是均匀的，它对应于零假设下p值后面的分布。因此，热图聚类可以归结为独立于任何生物同质性的基因聚类。

在单独的窗口中打开

图4

热图聚类中与生物同质性指标相关的p值QQ图。这里，我们重点关注与从鸡表达数据的热图聚类中获得的聚类生物同质性指标相关的p值。QQ图与均匀分布（x轴）内预期的概率相交，p值来自生物均匀性指标（y轴）。

其次，与热图相比，考虑共表达网络大大改善了结果。因此，WGCNA提供了更高比例的生物均质集群（68%）。然而，共存簇的比例下降。最终，WGCNA提供了合理比例的优秀口译候选人（46%）。

第三，与其他两种方法相比，使用我们自己的聚类算法，显著共存的簇所占比例降低。这是意料之中的，因为共表达与生物同质性竞争。然而，具有明显生物均质性的集群的比例显著增加（79.50%）。这导致优秀口译候选人的比例更高（53.50%）。

注：由一个基因组成的簇自动被视为坏候选。因此，由于我们的聚类策略提供了这些聚类中不可忽略的一部分，因此优秀候选的百分比在机械上较低。

总之，通过将生物学知识集成到表达数据中，我们成功地获得了合理比例的簇，这些簇收集了显著共存的生物相关基因。这些簇是很好的候选，它们的解释可能会揭示基因之间的新关系。

集群解释

通过将生物学知识集成到表达数据中而获得的聚类，以及呈现有趣特性的聚类，是很好的解释候选。为了将具有代表性的GO注释与聚类关联起来，我们选择应用经典的富集测试程序，该程序包括与多重测试校正相关的fisher精确测试（Benjamini-Hochberg，截止值为5%）。对浓缩过程结果的总体印象是与聚类相关的GO注释的一致性。与一个集群关联的丰富GO注释在GO层次结构中非常接近。这直接传达了我们的程序所保证的基因簇的生物同质性。

与Désert等人的论文相比[15]，通过对集群的丰富注释，还强调了针对禁食实施的一般和众所周知的机制。此外，我们的程序还揭示了新的轨迹。例如，一些簇与磷脂和鞘氨醇脂质机制有关，其基因在禁食鸡中表达，Désert等人没有对此进行描述。这些簇聚集了几种与这些脂质水解有关的酶，从而释放脂肪酸。然后，我们认为鸡在禁食一段时间后，可以从质膜中摄取脂肪酸。

讨论和结论

我们提出了一种新的无监督基因聚类算法，该算法通过将生物知识集成到表达数据中，依赖于基因之间的新距离。为此，我们提出了一种基于共存生物功能概念的明智编码。由于一种生物功能可以被同化为一组参与该功能的基因，我们可以将一种共表达的生物功能同化为对该组共表达基因的限制。自然，这个距离是用来聚类基因的。

然后通过我们提出的两个指标来评估基因簇的属性，这两个指标允许量化共表达和生物同质性。一方面，共表达是通过基于基因间相关性的指标来评估的。该指标纯粹是经验指标，但非常方便且易于解释。另一方面，生物同质性是通过一个基于Cramér V系数的指标来测量的，该系数是从编码GO注释的矩阵中计算出来的。尽管这一指标有其局限性，因为生物的同质性主要依赖于生物学解释，但自动了解基因簇的生物学意义恰好是有用的。此外，我们提出假设检验，用p值来增强这些指标，以验证簇是否显著共存和生物均匀。

为了测试我们的聚类算法和评估过程，我们将其应用于模拟数据集和实际数据集。此外，为了定位我们的方法，我们将其与生物学家常用的两种基因聚类策略进行了比较：热图和基于共表达网络的聚类。

具体来说，我们的方法显示出一些局限性，因为它提供了由单个基因构成的相对重要的聚类比例。然而，它优于其他方法：积极将生物知识集成到表达式数据中，可以提供最佳候选分区比例最高的分区。这些簇似乎确实是很好的解释候选，因为可以证明与磷脂和鞘磷脂机制相关的簇。然而，最终的外部生物验证仍有待完成，包括进行更高级的生物解释。

附录

附录1：共表达指示剂的变化范围

共表达指标包括计算簇内基因相关性的平均值K（K）_我让我们回顾一下共表达指示剂的计算（方程式（3））：

\begin{array}{l} CI公司 ({K（K）}_{我}) = & \frac{1}{\frac{卡片 ({K（K）}_{我}) (卡片 ({K（K）}_{我}) - 1)}{2}} \sum_{k个 | k个 \in {K（K）}_{我}} \\ \times (\sum_{{k个}^{'} | {k个}^{'} \in {K（K）}_{我}, {k个}^{'} > k个} \frac{1}{我} \sum_{我 = 1}^{我} (\frac{{G公司}_{伊克} - {G公司}_{.k（千）}}{{S公司}_{k个}}) (\frac{{G公司}_{我 {k个}^{'}} - {G公司}_{. {k个}^{'}}}{{S公司}_{{k个}^{'}}})) \end{array}

CI的最小值根据c（c）一第页d日(K（K）_我). 为了在K（K）_我，我们考虑两个子组，如组内相关等于1和组间相关等于-1。的所有基因K（K）_我在两个子组之间平均分布。

如果卡片(K（K）_我)是偶数

在这种情况下，每个子组由 $\frac{卡片 ({K（K）}_{我})}{2}$ 基因。负相关的最大数量等于 $\frac{卡片 ({K（K）}_{我})}{2} \times \frac{卡片 ({K（K）}_{我})}{2}$ .

CI公司 ({K（K）}_{我}) = \frac{[\frac{卡片 ({K（K）}_{我}) (卡片 ({K（K）}_{我}) - 1)}{2} - {(\frac{卡片 ({K（K）}_{我})}{2})}^{2}] - {(\frac{卡片 ({K（K）}_{我})}{2})}^{2}}{\frac{卡片 ({K（K）}_{我}) (卡片 ({K（K）}_{我}) - 1)}{2}}

CI公司 ({K（K）}_{我}) = - \frac{1}{卡片 ({K（K）}_{我}) - 1}

如果卡片(K（K）_我)很奇怪

在这种情况下，其中一个子组由 $\frac{卡片 ({K（K）}_{我}) - 1}{2}$ 基因，另一个由 $\frac{卡片 ({K（K）}_{我}) + 1}{2}$ .最大负相关数等于 $\frac{卡片 ({K（K）}_{我}) - 1}{2} \times \frac{卡片 ({K（K）}_{我}) + 1}{2}$ .

\begin{array}{l} CI公司 ({K（K）}_{我}) & = & \frac{[\frac{卡片 ({K（K）}_{我}) (卡片 ({K（K）}_{我}) - 1)}{2} - \frac{卡片 ({K（K）}_{我}) - 1}{2} \times \frac{卡片 ({K（K）}_{我}) + 1}{2}] - \frac{卡片 ({K（K）}_{我}) - 1}{2} \times \frac{卡片 ({K（K）}_{我}) + 1}{2}}{\frac{卡片 ({K（K）}_{我}) (卡片 ({K（K）}_{我}) - 1)}{2}} \\ CI公司 ({K（K）}_{我}) & = & - \frac{1}{卡片 ({K（K）}_{我})} \end{array}

当所有基因都存在时，CI最大，等于1K（K）_我完全正相关。

附录2：生物同质性指标的变化范围

让我们回顾一下生物均匀性指标的计算（方程式（4））：

制动马力 ({K（K）}_{我}) = 1 - \sqrt{\frac{\sum_{k个 \in {K（K）}_{我}} (\sum_{j个 = 1}^{J型} \frac{{({T型}_{千焦} - \frac{{T型}_{k、。} {T型}_{.j英寸}}{{T型}_{..}})}^{2}}{\frac{{T型}_{k、。} {T型}_{.j英寸}}{{T型}_{..}}})}{{T型}_{..} (卡片 ({K（K）}_{我}) - 1)}}

哪里T型_k个j个如果基因等于1k个与生物功能有关j个，否则为0，T型_k个.是与基因相关的行边距k、。

BHI最小，当K（K）_我具有相似的功能特征，以至于没有任何生物功能与以下两个基因相关K（K）_我:

BHI公司 ({K（K）}_{我}) = 1 - \sqrt{\frac{\sum_{k个 \in {K（K）}_{我}} (\sum_{j个 = 1}^{J型} \frac{{({T型}_{千焦} - \frac{{T型}_{k、。} {T型}_{.j英寸}}{{T型}_{..}})}^{2}}{\frac{{T型}_{k、。} {T型}_{.j英寸}}{{T型}_{..}}})}{{T型}_{..} (卡片 ({K（K）}_{我}) - 1)}}

∀j个|T型_k个j个=1,T型_.j个=1

∀j个|T型_k个j个=0,T型_.j个=0

BHI公司 ({K（K）}_{我}) = 1 - \sqrt{\frac{\sum_{k个 \in {K（K）}_{我}} ({T型}_{k、。} \frac{{(1 - \frac{{T型}_{k、。}}{{T型}_{..}})}^{2}}{\frac{{T型}_{k、。}}{{T型}_{..}}} + ({T型}_{..} - {T型}_{k、。}) \frac{{T型}_{k、。}}{{T型}_{..}})}{{T型}_{..} (卡片 ({K（K）}_{我}) - 1)}}

BHI公司 ({K（K）}_{我}) = 1 - \sqrt{\frac{\sum_{k个 \in {K（K）}_{我}} ({T型}_{..} {(1 - \frac{{T型}_{k、。}}{{T型}_{..}})}^{2} + {T型}_{k、。} - \frac{{T型}_{k、。}^{2}}{{T型}_{..}})}{{T型}_{..} (卡片 ({K（K）}_{我}) - 1)}}

BHI公司 ({K（K）}_{我}) = 1 - \sqrt{\frac{\sum_{k个 \in {K（K）}_{我}} (\frac{{T型}_{..}^{2} - 2 {T型}_{..} {T型}_{k、。} + {T型}_{k、。}^{2} + {T型}_{..} {T型}_{k、。} - {T型}_{k、。}^{2}}{{T型}_{..}})}{{T型}_{..} (卡片 ({K（K）}_{我}) - 1)}}

BHI公司 ({K（K）}_{我}) = 1 - \sqrt{\frac{\sum_{k个 \in {K（K）}_{我}} {T型}_{..} - \sum_{k个 \in {K（K）}_{我}} {T型}_{k、。}}{{T型}_{..} (卡片 ({K（K）}_{我}) - 1)}}

BHI公司 ({K（K）}_{我}) = 1 - \sqrt{\frac{卡片 ({k个}_{我}) {T型}_{..} - {T型}_{..}}{{T型}_{..} (卡片 ({K（K）}_{我}) - 1)}}

制动马力(K（K）_我) = 0

当所有基因K（K）_我具有完全相似的功能配置文件：

BHI公司 ({K（K）}_{我}) = 1 - \sqrt{\frac{\sum_{k个 \in {K（K）}_{我}} (\sum_{j个 = 1}^{J型} \frac{{({T型}_{千焦} - \frac{{T型}_{k、。} {T型}_{.j英寸}}{{T型}_{..}})}^{2}}{\frac{{T型}_{k、。} {T型}_{.j英寸}}{{T型}_{..}}})}{{T型}_{..} (卡片 ({K（K）}_{我}) - 1)}}

∀j个|T型_k个j个=1,T型_.j个=c（c）一第页d日(K（K）_我) &

{T型}_{.k（千）} = \frac{T.公司。。}{卡片 {K（K）}_{我}}

∀j个|T型_k个j个＝0时，T型_.j个=0

因此：

BHI公司 ({K（K）}_{我}) = 1 - \sqrt{\frac{\sum_{k个 \in {K（K）}_{我}} (\sum_{j个 = 1}^{J型} \frac{{(1 - \frac{\frac{{T型}_{..}}{卡片 ({K（K）}_{我})} 卡片 ({K（K）}_{我})}{T.公司。。})}^{2}}{\frac{\frac{{T型}_{..}}{卡片 ({K（K）}_{我})} 卡片 ({K（K）}_{我})}{T.公司。。}})}{{T型}_{..} (卡片 ({K（K）}_{我}) - 1)}}

BHI公司(K（K）_我) = 1

竞争性利益

作者声明，他们没有相互竞争的利益。

作者的贡献

MV、SL和JP制定了方法并起草了手稿。MV实现了该算法。所有作者都批准了最后的手稿。

致谢

作者们感谢来自Ouest农业大学遗传系的Sandrine Lagarrigue，感谢她能够提供并让他们使用她的数据。作者感谢审稿人的宝贵意见。

工具书类

Eisen MB、Spellman PT、Brown PO、Botstein D.全基因组表达模式的聚类分析和显示。国家科学院院刊。1998;95(25):14863–14868. doi:10.1073/pnas.95.25.14863。 [PMC免费文章][公共医学] [交叉参考][谷歌学者]
Zhang B，Horvath S.加权基因共表达网络分析的一般框架。统计应用基因分子生物学。2005;4：第17条。[公共医学][谷歌学者]
Ashburner M、Ball CA、Blake JA、Botstein D、Butler H、Cherry JM、Davis AP、Dolinski K、Dwight SS、Eppig JT、Harris MA、Hill DP、Issel-Tarver L、Kasarskis A、Lewis S、Matese JC、Richardson JE、Ringwald M、Rubin GM、Sherlock G.基因本体论：生物学统一的工具。基因本体联盟。自然遗传学。2000;25:25–29. doi:10.1038/75556。 [PMC免费文章][公共医学] [交叉参考][谷歌学者]
Bauer S，Gagneur J，Robinson PN。GOing Bayesian：基因组尺度数据的基于模型的基因集分析。核酸研究。2010;38:3523–3532. doi:10.1093/nar/gkq045。 [PMC免费文章][公共医学] [交叉参考][谷歌学者]
Busold CH、Winter S、Hauser N、Bauer A、Dippon J、Hoheisel JD、Fellenberg K。GO注释在对应分析中的整合：促进微阵列数据的解释。生物信息学。2005;21（10）：2424–2429。doi:10.1093/bioinformatics/bti367。[公共医学] [交叉参考][谷歌学者]
Fagan A、Culhane AC、Higgins DG。蛋白质组和基因表达数据集成的多元分析方法。蛋白质组学。2007;7(13):2162–2171. doi:10.1002/pmic.200600898。[公共医学] [交叉参考][谷歌学者]
Yaung MKS，Tegnér J，Collins JJ。使用奇异值分解和稳健回归对基因网络进行反向工程。美国国家科学院院刊。2002;99(9):6163–6168. doi:10.1073/pnas.092576199。 [PMC免费文章][公共医学] [交叉参考][谷歌学者]
Bryan J.基于基因表达数据的基因聚类问题。《多元分析杂志》。2004年；90:44–66. doi:10.1016/j.jmva.2004.02.011。[交叉参考][谷歌学者]
Kashima H、Yamanishi Y、Kato T、Sugiyama M、Tsuda K。从全基因组数据和进化信息同时推断多个物种的生物网络：半监督方法。生物信息学。2009;25(22):2962–2968. doi:10.1093/bioinformatics/btp494。[公共医学] [交叉参考][谷歌学者]
Azuaje F，Wang H，Zheng H，Léonard F，Rolland-Turner M，Zhang L，Devaux Y，Wagner D。基因功能相似性和共同表达的预测整合定义了内皮祖细胞的治疗反应。BMC系统生物。2011;5：46.网址：10.1186/1752-0509-5-46。 [PMC免费文章][公共医学] [交叉参考][谷歌学者]
Li XL，Tan YC，Ng SK.使用模糊最近聚类法从基因表达数据进行系统基因功能预测。BMC生物信息学。2006年；7（补充4）：S23。doi:10.1186/1471-2105-7-S4-S23。 [PMC免费文章][公共医学] [交叉参考][谷歌学者]
Datta S，Datta S.使用功能类参考集评估基因表达数据聚类算法的方法。BMC生物信息学。2006年；7：397.网址：10.1186/1471-2105-7-397。 [PMC免费文章][公共医学] [交叉参考][谷歌学者]
克拉梅·H。统计数学方法（PMS-9）新泽西：普林斯顿大学出版社；1945[谷歌学者]
Dray S.关于主成分的数量：基于矩阵之间相似性度量的维度测试。计算统计数据分析。2008;52（4）：2228–2237。doi:10.1016/j.csda.2007.07.015。[交叉参考][谷歌学者]
Désert C、Duclos M、Blavy P、Lecerf F、Moreews F、Klopp C、Aubry M、Herault F、Le RoyP、Berri C、Douaire M、Diot C、Lagarrigue S。鸡肝脏从喂养到禁食转变的转录组分析。BMC基因组学。2008;9：611.网址：10.1186/1471-2164-9-611。 [PMC免费文章][公共医学] [交叉参考][谷歌学者]
Friguet C，Kloareg M，Causeur D.依赖性下多重测试的因子模型方法。美国统计协会。2009;104(488):1406–1415. doi:10.1198/jasa.2009.tm08332。[交叉参考][谷歌学者]

文章来自BMC生物信息学由以下人员提供BMC公司

一种新的基于生物知识与表达数据集成的无监督基因聚类算法

玛丽·韦班克

塞巴斯蒂安·莱伊

杰罗姆·帕格斯

摘要

背景

结果

结论

背景

方法

将生物学知识集成到表达数据中：生物学原理

无监督基因聚类算法

生物知识的编码

基因间的新距离：共表达生物功能

获取基因簇

基因簇评估

共表达指示器

生物均匀性指示剂

假设检验程序

结果

模拟研究

模拟数据集

结果

表1

鸡肉数据集分析

表2

集群解释

讨论和结论

附录

附录1：共表达指示剂的变化范围

如果卡片(K（K）我)是偶数

如果卡片(K（K）我)很奇怪

附录2：生物同质性指标的变化范围

竞争性利益

作者的贡献

致谢

工具书类

如果卡片(K（K）_我)是偶数

如果卡片(K（K）_我)很奇怪