聚类分析

聚类分析是将项目分组到聚类中，基于项目之间彼此的相似性。在生物信息学中，聚类广泛用于基因表达数据分析，以找到具有相似基因表达谱的基因组。这可能会识别出功能相关的基因，并提示当前未知基因的功能。

Biopython 模块 Bio.Cluster 提供常用的聚类算法，并针对基因表达数据的应用进行了设计。但是，此模块也可用于对其他类型的数据进行聚类分析。Bio.Cluster 和底层的 C 聚类库由 De Hoon 等人描述 [DeHoon2004].

以下四种聚类方法在 Bio.Cluster 中实现

• 层次聚类（成对质心、单链接、完全链接和平均链接）；

• \(k\)-均值、\(k\)-中位数和 \(k\)-medoids 聚类；

• 自组织映射；

• 主成分分析。

数据表示

要聚类的数据用一个 \(n \times m\) 数值 Python 数组 data 表示。在基因表达数据聚类的情况下，通常行对应于不同的基因，而列对应于不同的实验条件。Bio.Cluster 中的聚类算法既可以应用于行（基因），也可以应用于列（实验）。

缺失值

\(n \times m\) 数值 Python 整数数组 mask 指示 data 中的任何值是否缺失。如果 mask[i, j] == 0，则 data[i, j] 缺失，并在分析中被忽略。

随机数生成器

\(k\)-均值/中位数/medoids 聚类算法和自组织映射 (SOM) 包括使用随机数生成器。Bio.Cluster 中的均匀随机数生成器基于 L’Ecuyer 的算法 [Lecuyer1988]，而遵循二项分布的随机数是使用 Kachitvichyanukul 和 Schmeiser 的 BTPE 算法生成的 [Kachitvichyanukul1988]。随机数生成器在其第一次调用时自动初始化。由于此随机数生成器使用两个乘法线性同余生成器的组合，因此初始化需要两个（整数）种子，为此我们使用系统提供的随机数生成器 rand（在 C 标准库中）。我们通过调用 srand 并使用以秒为单位的纪元时间来初始化此生成器，并将 rand 生成的前两个随机数用作 Bio.Cluster 中的均匀随机数生成器的种子。

距离函数

为了将项目聚集成组，我们应该首先定义我们对“相似”的精确含义。Bio.Cluster 提供八个距离函数，用单个字符表示，来衡量相似性，或者相反，衡量距离

• 'e'：欧几里得距离；

• 'b'：城市街区距离。

• 'c'：皮尔森相关系数；

• 'a'：皮尔森相关系数的绝对值；

• 'u'：非中心化皮尔森相关性（等同于两个数据向量之间的角度的余弦）；

• 'x'：绝对非中心化皮尔森相关性；

• 's'：斯皮尔曼秩相关；

• 'k'：肯德尔的 \(\tau\).

前两个是满足三角不等式的真实距离函数

\[d\left(\underline{u},\underline{v}\right) \leq d\left(\underline{u},\underline{w}\right) + d\left(\underline{w},\underline{v}\right) \textrm{ for all } \underline{u}, \underline{v}, \underline{w},\]

因此被称为度量。在日常语言中，这意味着两点之间的最短距离是直线。

其余六个距离度量与相关系数有关，其中距离 \(d\) 是用相关性 \(r\) 定义的，由 \(d=1-r\) 表示。请注意，这些距离函数是半度量，不满足三角不等式。例如，对于

\[\underline{u}=\left(1,0,-1\right);\]

\[\underline{v}=\left(1,1,0\right);\]

\[\underline{w}=\left(0,1,1\right);\]

我们发现皮尔森距离 \(d\left(\underline{u},\underline{w}\right) = 1.8660\)，而 \(d\left(\underline{u},\underline{v}\right)+d\left(\underline{v},\underline{w}\right) = 1.6340\).

欧几里得距离

在 Bio.Cluster 中，我们将欧几里得距离定义为

\[d = {1 \over n} \sum_{i=1}^{n} \left(x_i-y_i\right)^{2}.\]

只有 \(x_i\) 和 \(y_i\) 都存在时，这些项才包含在求和中，并且分母 \(n\) 应相应地选择。由于表达式数据 \(x_i\) 和 \(y_i\) 直接从彼此中减去，因此在使用欧几里得距离时，我们应该确保表达式数据已正确标准化。

城市街区距离

城市街区距离，也称为曼哈顿距离，与欧几里得距离有关。虽然欧几里得距离对应于两点之间最短路径的长度，但城市街区距离是沿每个维度的距离之和。由于基因表达数据往往具有缺失值，因此在 Bio.Cluster 中，我们将城市街区距离定义为距离之和除以维数

\[d = {1 \over n} \sum_{i=1}^n \left|x_i-y_i\right|.\]

这等于您必须在城市中走过的两点之间的距离，您必须沿着城市街区行走。对于欧几里得距离，表达式数据直接从彼此中减去，因此我们应该确保它们已正确标准化。

皮尔森相关系数

皮尔森相关系数定义为

\[r = \frac{1}{n} \sum_{i=1}^n \left( \frac{x_i -\bar{x}}{\sigma_x} \right) \left(\frac{y_i -\bar{y}}{\sigma_y} \right),\]

其中 \(\bar{x}, \bar{y}\) 分别是 \(x\) 和 \(y\) 的样本均值，\(\sigma_x, \sigma_y\) 分别是 \(x\) 和 \(y\) 的样本标准差。皮尔逊相关系数用于衡量一条直线拟合 \(x\) 和 \(y\) 散点图的程度。如果散点图中的所有点都位于一条直线上，则皮尔逊相关系数为 +1 或 -1，取决于直线的斜率是正还是负。如果皮尔逊相关系数等于零，则 \(x\) 和 \(y\) 之间不存在相关性。

皮尔逊距离 定义为

\[d_{\textrm{P}} \equiv 1 - r.\]

由于皮尔逊相关系数介于 -1 和 1 之间，因此皮尔逊距离介于 0 和 2 之间。

绝对皮尔逊相关

通过取皮尔逊相关的绝对值，我们得到一个介于 0 和 1 之间的数字。如果绝对值为 1，则散点图中的所有点都位于一条具有正或负斜率的直线上。如果绝对值等于零，则 \(x\) 和 \(y\) 之间不存在相关性。

相应的距离定义为

\[d_{\textrm A} \equiv 1 - \left|r\right|,\]

其中 \(r\) 是皮尔逊相关系数。由于皮尔逊相关系数的绝对值介于 0 和 1 之间，因此相应的距离也介于 0 和 1 之间。

在基因表达实验的背景下，如果两个基因的基因表达谱完全相同或完全相反，则绝对相关性等于 1。因此，应谨慎使用绝对相关系数。

非中心相关（角度的余弦）

在某些情况下，可能更喜欢使用非中心相关 而不是常规的皮尔逊相关系数。非中心相关定义为

\[r_{\textrm U} = \frac{1}{n} \sum_{i=1}^{n} \left(\frac{x_i}{\sigma_x^{(0)}} \right) \left(\frac{y_i}{\sigma_y^{(0)}} \right),\]

其中

\[\begin{split}\begin{aligned} \sigma_x^{(0)} & = & \sqrt{{\frac{1}{n}} \sum_{i=1}^{n}x_i^2}; \nonumber \\ \sigma_y^{(0)} & = & \sqrt{{\frac{1}{n}} \sum_{i=1}^{n}y_i^2}. \nonumber \end{aligned}\end{split}\]

这与常规皮尔逊相关系数的表达式相同，只是样本均值 \(\bar{x}, \bar{y}\) 被设置为零。如果存在零参考状态，则非中心相关可能适用。例如，在以对数比率形式给出的基因表达数据的情况下，等于零的对数比率对应于绿色和红色信号相等，这意味着实验操作没有影响基因表达。

对应于非中心相关系数的距离定义为

\[d_{\mbox{U}} \equiv 1 - r_{\mbox{U}},\]

其中 \(r_{\mbox{U}}\) 是非中心相关。由于非中心相关系数介于 -1 和 1 之间，因此相应的距离介于 0 和 2 之间。

非中心相关等于 \(n\) 维空间中两个数据向量角度的余弦，通常被称为余弦相似度。

绝对非中心相关

与常规皮尔逊相关一样，我们可以使用非中心相关的绝对值来定义距离度量

\[d_{\mbox{AU}} \equiv 1 - \left|r_{\mbox{U}}\right|,\]

其中 \(r_{\mbox{U}}\) 是非中心相关系数。由于非中心相关系数的绝对值介于 0 和 1 之间，因此相应的距离也介于 0 和 1 之间。

在几何上，非中心相关的绝对值等于两个数据向量的支撑线之间的余弦（即不考虑向量方向的角度）。

斯皮尔曼等级相关

斯皮尔曼等级相关是无参数相似度度量的一个例子，它比皮尔逊相关更能抵抗异常值。

要计算斯皮尔曼等级相关，我们将每个数据值替换为其等级，如果我们按照每个向量中的值对数据进行排序。然后，我们计算两个等级向量之间的皮尔逊相关，而不是数据向量。

与皮尔逊相关一样，我们可以定义一个对应于斯皮尔曼等级相关的距离度量，如下所示

\[d_{\mbox{S}} \equiv 1 - r_{\mbox{S}},\]

其中 \(r_{\mbox{S}}\) 是斯皮尔曼等级相关。

肯德尔的 \(\tau\)

肯德尔的 \(\tau\) 是无参数相似度度量的另一个例子。它类似于斯皮尔曼等级相关，但它不使用等级本身，而是使用相对等级来计算 \(\tau\)（参见 Snedecor & Cochran [Snedecor1989]）。

我们可以定义一个对应于肯德尔的 \(\tau\) 的距离度量，如下所示

\[d_{\mbox{K}} \equiv 1 - \tau.\]

由于肯德尔的 \(\tau\) 始终介于 -1 和 1 之间，因此相应的距离将在 0 和 2 之间。

加权

对于 Bio.Cluster 中提供的多数距离函数，可以应用权重向量。权重向量包含数据向量中项目的权重。如果项目 \(i\) 的权重为 \(w_i\)，则该项目将被视为在数据中出现 \(w_i\) 次。权重不必是整数。

计算距离矩阵

距离矩阵是一个方阵，其中包含 data 中项目之间所有成对距离，可以通过 Bio.Cluster 模块中的 distancematrix 函数计算

>>> from Bio.Cluster import distancematrix
>>> matrix = distancematrix(data)

其中定义了以下参数

• data（必需）

  包含项目数据的数组。
• mask（默认：None）

  整数数组，显示哪些数据丢失。如果 mask[i, j] == 0，则 data[i, j] 丢失。如果 mask 为 None，则所有数据都存在。
• weight（默认：None）

  计算距离时使用的权重。如果 weight 为 None，则假定权重相等。
• transpose（默认：0）

  确定是否要计算 data 行之间的距离（transpose 为 False），还是计算 data 列之间的距离（transpose 为 True）。
• dist（默认：'e'，欧几里得距离）

  定义要使用的距离函数（参见 距离函数）。

为了节省内存，距离矩阵将作为一维数组列表返回。每一行中的列数等于行号。因此，第一行有零个元素。例如，

>>> from numpy import array
>>> from Bio.Cluster import distancematrix
>>> data = array([[0, 1,  2,  3],
...               [4, 5,  6,  7],
...               [8, 9, 10, 11],
...               [1, 2,  3,  4]])  # fmt: skip
...
>>> distances = distancematrix(data, dist="e")

生成距离矩阵

>>> distances
[array([], dtype=float64), array([ 16.]), array([ 64.,  16.]), array([  1.,   9.,  49.])]

可以改写为

[array([], dtype=float64), array([16.0]), array([64.0, 16.0]), array([1.0, 9.0, 49.0])]

这对应于距离矩阵

\[\begin{split}\left( \begin{array}{cccc} 0 & 16 & 64 & 1 \\ 16 & 0 & 16 & 9 \\ 64 & 16 & 0 & 49 \\ 1 & 9 & 49 & 0 \end{array} \right).\end{split}\]

计算聚类属性

计算聚类中心

聚类的中心可以定义为所有聚类项目在每个维度上的平均值或中位数。Bio.Cluster 中的 clustercentroids 函数可用于计算以下任一值

>>> from Bio.Cluster import clustercentroids
>>> cdata, cmask = clustercentroids(data)

其中定义了以下参数

• data（必需）

  包含项目数据的数组。
• mask（默认：None）

  整数数组，显示哪些数据丢失。如果 mask[i, j] == 0，则 data[i, j] 丢失。如果 mask 为 None，则所有数据都存在。
• clusterid（默认：None）

  整数向量，显示每个项目属于哪个聚类。如果 clusterid 为 None，则假定所有项目都属于同一个聚类。
• method（默认：'a'）

  指定使用算术平均值（method=='a'）还是中位数（method=='m'）来计算聚类中心。
• transpose（默认：0）

  确定是否要计算 data 行的中心（transpose 为 False），还是计算 data 列的中心（transpose 为 True）。

此函数返回元组 (cdata, cmask)。质心数据存储在二维 NumPy 数组 cdata 中，缺失数据由二维 NumPy 整数数组 cmask 表示。如果 transpose 为 0，则这些数组的维度为 \(\left(\textrm{number of clusters}, \textrm{number of columns}\right)\)，如果 transpose 为 1，则维度为 \(\left(\textrm{number of rows}, \textrm{number of clusters}\right)\)。每行（如果 transpose 为 0）或每列（如果 transpose 为 1）包含与每个聚类的质心相对应的平均数据。

计算聚类之间的距离

给定项之间的距离函数，我们可以用几种方法定义两个聚类之间的距离。成对质心连接聚类和 k 均值聚类中使用两个聚类的算术平均值之间的距离。在 k 中值聚类中，使用两个聚类中值之间的距离。在成对单连接聚类中，使用两个聚类中项之间的最短成对距离，而在成对最大连接聚类中，使用最长的成对距离。在成对平均连接聚类中，两个聚类之间的距离定义为成对距离的平均值。

要计算两个聚类之间的距离，请使用

>>> from Bio.Cluster import clusterdistance
>>> distance = clusterdistance(data)

其中定义了以下参数

• data（必需）

  包含项目数据的数组。
• mask（默认：None）

  整数数组，显示哪些数据丢失。如果 mask[i, j] == 0，则 data[i, j] 丢失。如果 mask 为 None，则所有数据都存在。
• weight（默认：None）

  计算距离时使用的权重。如果 weight 为 None，则假定权重相等。
• index1（默认值：0）

  包含属于第一个聚类的项的索引的列表。包含一个项 \(i\) 的聚类可以表示为列表 [i] 或整数 i。
• index2（默认值：0）

  包含属于第二个聚类的项的索引的列表。包含一个项 \(i\) 的聚类可以表示为列表 [i] 或整数 i。
• method（默认：'a'）

  指定如何定义聚类之间的距离

  • 'a'：两个聚类质心之间的距离（算术平均值）；

  • 'm'：两个聚类质心之间的距离（中值）；

  • 's'：两个聚类中项之间的最短成对距离；

  • 'x'：两个聚类中项之间的最长成对距离；

  • 'v'：两个聚类中项之间的成对距离的平均值。

• dist（默认：'e'，欧几里得距离）

  定义要使用的距离函数（参见 距离函数）。
• transpose（默认：0）

  如果 transpose 为 False，则计算 data 行之间的距离。如果 transpose 为 True，则计算 data 列之间的距离。

分区算法

分区算法将项划分为 \(k\) 个聚类，使项到其聚类中心的距离总和最小。聚类数 \(k\) 由用户指定。Bio.Cluster 中提供了三种分区算法

• \(k\) 均值聚类

• \(k\) 中值聚类

• \(k\)-medoids 聚类

这些算法在聚类中心的定义方式上有所不同。在 \(k\) 均值聚类中，聚类中心定义为对聚类中所有项求平均得到的平均数据向量。在 \(k\) 中值聚类中，不是使用均值，而是计算数据向量中每个维度的中值。最后，在 \(k\) 中值聚类中，聚类中心定义为距离聚类中其他项的距离总和最小的项。此聚类算法适用于距离矩阵已知但原始数据矩阵不可用的情况，例如，根据蛋白质的结构相似性对蛋白质进行聚类时。

期望最大化 (EM) 算法用于找到这种将项划分为 \(k\) 个组的分区。在 EM 算法的初始化阶段，我们将项随机分配给聚类。为了确保不会产生空聚类，我们使用二项式分布随机选择每个聚类中的项数，使其至少为一个。然后，我们随机排列项的聚类分配，使每个项在任何聚类中都有相同的概率。因此，每个聚类都保证至少包含一个项。

然后，我们迭代

• 计算每个聚类的质心，定义为聚类的均值、中值或中值；

• 计算每个项到聚类中心的距离；

• 对于每个项，确定哪个聚类质心最接近；

• 将每个项重新分配到其最接近的聚类，或者如果不再进行项重新分配，则停止迭代。

为了避免聚类在迭代期间变为空，在 \(k\) 均值和 \(k\) 中值聚类中，算法跟踪每个聚类中的项数，并禁止聚类中最后一个剩余的项重新分配到另一个聚类。对于 \(k\) 中值聚类，不需要进行这种检查，因为充当聚类质心的项与其自身的距离为零，因此永远不会更接近于另一个聚类。

由于项到聚类的初始分配是随机完成的，因此通常每次执行 EM 算法时都会找到不同的聚类解决方案。为了找到最佳聚类解决方案，\(k\) 均值算法重复执行多次，每次从不同的初始随机聚类开始。对于每次运行，都会保存项到其聚类中心的距离总和，并将该总和值最小的解决方案返回作为整体聚类解决方案。

EM 算法应该运行多少次取决于要聚类的项数。作为一个经验法则，我们可以考虑最佳解决方案找到的频率；此数字由本库中实现的分区算法返回。如果最佳解决方案找到了很多次，那么找到比已找到的解决方案更好的解决方案的可能性很小。但是，如果最佳解决方案只找到了一次，那么可能存在其他解决方案，其聚类内距离总和更小。如果项数很大（超过几百个），则可能难以找到全局最优解。

当不再进行重新分配时，EM 算法结束。我们注意到，对于某些初始聚类分配集，EM 算法无法收敛，因为在执行少量迭代步骤后，相同的聚类解决方案会周期性地重新出现。因此，我们在迭代期间检查这种周期性解决方案的出现。执行一定数量的迭代步骤后，当前聚类结果将保存为参考。通过将每次后续迭代步骤后的聚类结果与参考状态进行比较，我们可以确定是否找到了以前遇到的聚类结果。在这种情况下，迭代将停止。如果在执行一定数量的迭代步骤后仍未遇到参考状态，则当前聚类解决方案将保存以用作新的参考状态。最初，执行十个迭代步骤，然后重新保存参考状态。每次将此迭代步骤数加倍，以确保也可以检测到周期较长的周期性行为。

\(k\) 均值和 \(k\) 中值

\(k\) 均值和 \(k\) 中值算法在 Bio.Cluster 中实现为函数 kcluster

>>> from Bio.Cluster import kcluster
>>> clusterid, error, nfound = kcluster(data)

其中定义了以下参数

• data（必需）

  包含项目数据的数组。
• nclusters（默认值：2）

  聚类数 \(k\)。
• mask（默认：None）

  整数数组，显示哪些数据丢失。如果 mask[i, j] == 0，则 data[i, j] 丢失。如果 mask 为 None，则所有数据都存在。
• weight（默认：None）

  计算距离时使用的权重。如果 weight 为 None，则假定权重相等。
• transpose（默认：0）

  确定要聚类的行（transpose 为 0）还是列（transpose 为 1）。
• npass（默认值：1）

  执行 \(k\) 均值/中值聚类算法的次数，每次使用不同的（随机）初始条件。如果给出了 initialid，则忽略 npass 的值，并且聚类算法只运行一次，因为它在这种情况下表现出确定性。
• method（默认值：a）

  描述如何找到聚类的中心

  • method=='a'：算术平均值（\(k\) 均值聚类）；

  • method=='m'：中值（\(k\) 中值聚类）。

  对于 method 的其他值，使用算术平均值。

• dist（默认：'e'，欧几里得距离）

  定义要使用的距离函数（请参阅 距离函数）。虽然 kcluster 接受所有八种距离度量，但从理论角度来看，最好对 \(k\) 均值算法使用欧几里得距离，对 \(k\) 中值算法使用城市街区距离。
• initialid（默认值：None）

  指定用于 EM 算法的初始聚类。如果 initialid 为 None，则每次 EM 算法的 npass 次运行使用不同的随机初始聚类。如果 initialid 不是 None，则它应该等于一个 1D 数组，其中包含每个项目的聚类编号（介于 0 和 nclusters-1 之间）。每个聚类至少应包含一个项目。在指定初始聚类的情况下，EM 算法是确定性的。

此函数返回一个元组 (clusterid, error, nfound)，其中 clusterid 是一个包含每个行或聚类所分配的聚类编号的整数数组，error 是最佳聚类解决方案的簇内距离之和，而 nfound 是找到此最佳解决方案的次数。

\(k\)-medoids 聚类

kmedoids 例程对给定的一组项目执行 \(k\)-medoids 聚类，使用用户传递的距离矩阵和聚类数量

>>> from Bio.Cluster import kmedoids
>>> clusterid, error, nfound = kmedoids(distance)

其中定义了以下参数：，nclusters=2, npass=1, initialid=None)|

• distance (必需)

  包含项目之间距离的矩阵；此矩阵可以通过三种方式指定

  • 作为 2D 数值 Python 数组（其中只访问数组的左下部分）

    distance = array([[0.0, 1.1, 2.3], [1.1, 0.0, 4.5], [2.3, 4.5, 0.0]])
    

  • 作为 1D 数值 Python 数组，依次包含距离矩阵左下部分的距离

    distance = array([1.1, 2.3, 4.5])
    

  • 作为包含距离矩阵左下部分行的列表

    distance = [array([]), array([1.1]), array([2.3, 4.5])]
    

  这三个表达式对应于同一个距离矩阵。

• nclusters（默认值：2）

  聚类数 \(k\)。
• npass（默认值：1）

  \(k\)-medoids 聚类算法执行的次数，每次使用不同的（随机）初始条件。如果给出 initialid，则忽略 npass 的值，因为在这种情况下，聚类算法的行为是确定性的。
• initialid（默认值：None）

  指定用于 EM 算法的初始聚类。如果 initialid 为 None，则每次 EM 算法的 npass 次运行使用不同的随机初始聚类。如果 initialid 不是 None，则它应该等于一个 1D 数组，其中包含每个项目的聚类编号（介于 0 和 nclusters-1 之间）。每个聚类至少应包含一个项目。在指定初始聚类的情况下，EM 算法是确定性的。

此函数返回一个元组 (clusterid, error, nfound)，其中 clusterid 是一个数组，包含每个项目所分配的聚类编号，error 是最佳 \(k\)-medoids 聚类解决方案的簇内距离之和，而 nfound 是找到最佳解决方案的次数。请注意，clusterid 中的聚类编号定义为代表聚类中心的项目的项目编号。

层次聚类

层次聚类方法从本质上不同于 \(k\)-均值聚类方法。在层次聚类中，基因或实验条件之间表达谱的相似性以树形结构的形式表示。这种树形结构可以通过 Treeview 和 Java Treeview 等程序以图形方式显示，这有助于层次聚类在基因表达数据分析中的普及。

层次聚类的第一步是计算距离矩阵，指定要聚类的所有项目之间的距离。接下来，我们通过连接两个最近的项目来创建一个节点。后续节点通过根据项目或节点之间的距离成对连接项目或节点来创建，直到所有项目都属于同一个节点。然后，可以通过回溯哪些项目和节点已合并来创建树形结构。与用于 \(k\)-均值聚类的 EM 算法不同，层次聚类的整个过程是确定性的。

存在多种层次聚类方法，它们的不同之处在于如何根据成员定义子节点之间的距离。在 Bio.Cluster 中，可以使用成对单链接、最大链接、平均链接和质心链接。

• 在成对单链接聚类中，两个节点之间的距离定义为两个节点成员之间成对距离中最短的距离。

• 在成对最大链接聚类中，也称为成对完全链接聚类，两个节点之间的距离定义为两个节点成员之间成对距离中最长的距离。

• 在成对平均链接聚类中，两个节点之间的距离定义为两个节点中所有项目之间成对距离的平均值。

• 在成对质心链接聚类中，两个节点之间的距离定义为它们质心之间的距离。质心是通过对簇中的所有项目取平均值来计算的。由于在每一步都需要计算新形成的节点与现有节点和项目的距离，因此成对质心链接聚类的计算时间可能比其他层次聚类方法长得多。另一个特点是（对于基于 Pearson 相关性的距离度量），距离在向上进入聚类树时不一定增加，甚至可能减小。这是由于使用 Pearson 相关性时，质心计算和距离计算之间的不一致导致的：虽然 Pearson 相关性有效地对数据进行规范化以进行距离计算，但质心计算没有进行这种规范化。

对于成对单链接、完全链接和平均链接聚类，两个节点之间的距离可以直接从各个项目之间的距离中找到。因此，一旦知道距离矩阵，聚类算法就不需要访问原始基因表达数据。然而，对于成对质心链接聚类，新形成的子节点的质心只能从原始数据而不是从距离矩阵中计算出来。

成对单链接层次聚类的实现基于 SLINK 算法 [Sibson1973]，该算法比成对单链接聚类的直接实现快得多，而且内存效率更高。此算法生成的聚类结果与传统单链接算法找到的聚类解决方案相同。此库中实现的单链接层次聚类算法可用于对大型基因表达数据集进行聚类，而传统层次聚类算法由于内存需求过高和运行时间过长而无法处理这些数据集。

表示层次聚类解决方案

层次聚类的结果由一个节点树组成，其中每个节点连接两个项目或子节点。通常，我们不仅对每个节点连接哪些项目或子节点感兴趣，还对它们的相似性（或距离）感兴趣，因为它们是连接在一起的。为了在层次聚类树中存储一个节点，我们使用 Node 类，该类在 Bio.Cluster 中定义。Node 的实例具有三个属性

• left

• right

• distance

这里，left 和 right 是指在这个节点连接的两个项目或子节点的整数，而 distance 是它们之间的距离。要聚类的项目从 0 到 \(\left(\textrm{项目数} - 1\right)\) 编号，而簇从 -1 到 \(-\left(\textrm{项目数}-1\right)\) 编号。请注意，节点数量比项目数量少一个。

要创建一个新的 Node 对象，我们需要指定 left 和 right；distance 是可选的。

>>> from Bio.Cluster import Node
>>> Node(2, 3)
(2, 3): 0
>>> Node(2, 3, 0.91)
(2, 3): 0.91

现有 Node 对象的属性 left、right 和 distance 可以直接修改

>>> node = Node(4, 5)
>>> node.left = 6
>>> node.right = 2
>>> node.distance = 0.73
>>> node
(6, 2): 0.73

如果 left 和 right 不是整数，或者如果 distance 无法转换为浮点数，则会引发错误。

Python 类 Tree 表示完整的层次聚类解决方案。可以从 Node 对象列表创建 Tree 对象

>>> from Bio.Cluster import Node, Tree
>>> nodes = [Node(1, 2, 0.2), Node(0, 3, 0.5), Node(-2, 4, 0.6), Node(-1, -3, 0.9)]
>>> tree = Tree(nodes)
>>> print(tree)
(1, 2): 0.2
(0, 3): 0.5
(-2, 4): 0.6
(-1, -3): 0.9

Tree 初始化器检查节点列表是否为有效的层次聚类结果

>>> nodes = [Node(1, 2, 0.2), Node(0, 2, 0.5)]
>>> Tree(nodes)
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
ValueError: Inconsistent tree

可以使用方括号访问 Tree 对象中的各个节点

>>> nodes = [Node(1, 2, 0.2), Node(0, -1, 0.5)]
>>> tree = Tree(nodes)
>>> tree[0]
(1, 2): 0.2
>>> tree[1]
(0, -1): 0.5
>>> tree[-1]
(0, -1): 0.5

由于 Tree 对象是不可变的，因此我们无法更改 Tree 对象中的各个节点。但是，我们可以将树转换为节点列表，修改此列表，然后从此列表创建一棵新树

>>> tree = Tree([Node(1, 2, 0.1), Node(0, -1, 0.5), Node(-2, 3, 0.9)])
>>> print(tree)
(1, 2): 0.1
(0, -1): 0.5
(-2, 3): 0.9
>>> nodes = tree[:]
>>> nodes[0] = Node(0, 1, 0.2)
>>> nodes[1].left = 2
>>> tree = Tree(nodes)
>>> print(tree)
(0, 1): 0.2
(2, -1): 0.5
(-2, 3): 0.9

这保证任何 Tree 对象始终格式良好。

为了使用 Java Treeview 等可视化程序显示层次聚类解决方案，最好将所有节点距离缩放到 0 到 1 之间。这可以通过在现有 Tree 对象上调用 scale 方法来实现

>>> tree.scale()

此方法不接受任何参数，并返回 None。

在绘制树之前，您可能还希望重新排序树节点。\(n\) 个项目的层次聚类解决方案可以通过在每个节点上交换左右子节点来绘制为 \(2^{n-1}\) 个不同但等效的树状图。tree.sort(order) 方法访问层次聚类树中的每个节点，并验证左子节点的平均顺序值是否小于或等于右子节点的平均顺序值。如果不是，则交换左右子节点。这里，项目的顺序值由用户给出。在生成的树状图中，从左到右的项目往往会具有递增的顺序值。该方法将返回排序后从左到右的元素索引

>>> indices = tree.sort(order)

以便项目 indices[i] 将出现在树状图中的位置 \(i\)。

在层次聚类之后，可以通过剪切树将项目根据存储在 Tree 对象中的树结构分组到 \(k\) 个簇中

>>> clusterid = tree.cut(nclusters=1)

其中 nclusters（默认值为 1）是所需的聚类数量 \(k\)。此方法忽略树结构中排名前 \(k-1\) 的链接事件，从而导致 \(k\) 个分离的项目聚类。聚类数量 \(k\) 应为正数，且小于或等于项目数量。此方法返回一个数组 clusterid，其中包含分配给每个项目的聚类编号。聚类从 \(0\) 到 \(k-1\) 编号，其在树状图中的左右顺序排列。

执行层次聚类

要执行层次聚类，请使用 Bio.Cluster 中的 treecluster 函数。

>>> from Bio.Cluster import treecluster
>>> tree = treecluster(data)

其中定义了以下参数

• 数据

  包含项目数据的数组。
• mask（默认：None）

  整数数组，显示哪些数据丢失。如果 mask[i, j] == 0，则 data[i, j] 丢失。如果 mask 为 None，则所有数据都存在。
• weight（默认：None）

  计算距离时使用的权重。如果 weight 为 None，则假定权重相等。
• transpose（默认：0）

  确定是否对行（transpose 为 False）或列（transpose 为 True）进行聚类。
• method（默认值：'m'）

  定义要使用的链接方法

  • method=='s'：成对单链接聚类

  • method=='m'：成对最大（或完全）链接聚类

  • method=='c'：成对质心链接聚类

  • method=='a'：成对平均链接聚类

• dist（默认：'e'，欧几里得距离）

  定义要使用的距离函数（参见 距离函数）。

要对预先计算的距离矩阵应用层次聚类，请在调用 treecluster 函数时指定 distancematrix 参数，而不是 data 参数

>>> from Bio.Cluster import treecluster
>>> tree = treecluster(distancematrix=distance)

在这种情况下，将定义以下参数

• 距离矩阵

  距离矩阵，可以通过三种方式指定

  • 作为 2D 数值 Python 数组（其中只访问数组的左下部分）

    distance = array([[0.0, 1.1, 2.3], [1.1, 0.0, 4.5], [2.3, 4.5, 0.0]])
    

  • 作为 1D 数值 Python 数组，依次包含距离矩阵左下部分的距离

    distance = array([1.1, 2.3, 4.5])
    

  • 作为包含距离矩阵左下部分行的列表

    distance = [array([]), array([1.1]), array([2.3, 4.5])]
    

  这三个表达式对应于同一个距离矩阵。由于 treecluster 可能会在聚类算法的过程中对距离矩阵中的值进行洗牌，因此请确保在调用 treecluster 之前将此数组保存在另一个变量中，如果您稍后需要它。

• 方法

  要使用的链接方法

  • method=='s'：成对单链接聚类

  • method=='m'：成对最大（或完全）链接聚类

  • method=='a'：成对平均链接聚类

  虽然可以仅从距离矩阵中计算出成对单链接、最大链接和平均链接聚类，但不能计算出成对质心链接聚类。

在调用 treecluster 时，data 或 distancematrix 应为 None。

此函数返回一个 Tree 对象。此对象包含 \(\left(\textrm{项目数量} - 1\right)\) 个节点，其中项目数量是如果对行进行聚类，则为行数，如果对列进行聚类，则为列数。每个节点描述一个成对链接事件，其中节点属性 left 和 right 分别包含一个项目或子节点的编号，而 distance 则包含它们之间的距离。项目从 0 到 \(\left(\textrm{项目数量} - 1\right)\) 编号，而聚类从 -1 到 \(-\left(\textrm{项目数量}-1\right)\) 编号。

自组织映射

自组织映射 (SOM) 由 Kohonen 发明，用于描述神经网络（例如，参见 Kohonen，1997 [Kohonen1997]）。Tamayo (1999) 首次将自组织映射应用于基因表达数据 [Tamayo1999]。

SOM 将项目组织成位于某个拓扑结构中的聚类。通常选择矩形拓扑结构。SOM 生成的聚类使得拓扑结构中的相邻聚类彼此之间的相似度高于拓扑结构中彼此远离的聚类。

计算 SOM 的第一步是将数据向量随机分配给拓扑结构中的每个聚类。如果对行进行聚类，则每个数据向量中的元素数量等于列数。

然后通过一次取一行并找到拓扑结构中具有最接近数据向量的聚类来生成 SOM。该聚类的数据向量以及相邻聚类的数据向量将使用所考虑行的向量进行调整。调整由下式给出

\[\Delta \underline{x}_{\textrm{cell}} = \tau \cdot \left(\underline{x}_{\textrm{row}} - \underline{x}_{\textrm{cell}} \right).\]

参数 \(\tau\) 是一个在每次迭代步骤中都会减小的参数。我们使用了迭代步骤的简单线性函数

\[\tau = \tau_{\textrm{init}} \cdot \left(1 - {i \over n}\right),\]

\(\tau_{\textrm{init}}\) 是用户指定的 \(\tau\) 的初始值，\(i\) 是当前迭代步骤的编号，而 \(n\) 是要执行的总迭代步骤数。虽然在迭代开始时会快速进行更改，但在迭代结束时只会进行少量更改。

半径为 \(R\) 内的所有聚类都将根据所考虑的基因进行调整。随着计算的进行，该半径也会减小，如

\[R = R_{\textrm{max}} \cdot \left(1 - {i \over n}\right),\]

其中最大半径定义为

\[R_{\textrm{max}} = \sqrt{N_x^2 + N_y^2},\]

其中 \(\left(N_x, N_y\right)\) 是定义拓扑结构的矩形的尺寸。

函数 somcluster 实施了在矩形网格上计算自组织映射的完整算法。首先，它初始化随机数生成器。然后使用随机数生成器初始化节点数据。用于修改 SOM 的基因或样本的顺序也是随机的。SOM 算法中的总迭代次数由用户指定。

要运行 somcluster，请使用

>>> from Bio.Cluster import somcluster
>>> clusterid, celldata = somcluster(data)

其中定义了以下参数

• data（必需）

  包含项目数据的数组。
• mask（默认：None）

  整数数组，显示哪些数据丢失。如果 mask[i, j] == 0，则 data[i, j] 丢失。如果 mask 为 None，则所有数据都存在。
• weight（默认：None）

  包含计算距离时要使用的权重。如果 weight 为 None，则假定权重相等。
• transpose（默认：0）

  确定要聚类的行（transpose 为 0）还是列（transpose 为 1）。
• nxgrid, nygrid（默认值：2, 1）

  计算自组织映射的矩形网格中水平和垂直方向上的单元格数量。
• inittau（默认值：0.02）

  SOM 算法中使用的参数 \(\tau\) 的初始值。inittau 的默认值为 0.02，这是 Michael Eisen 的 Cluster/TreeView 程序中使用的值。
• niter（默认值：1）

  要执行的迭代次数。
• dist（默认：'e'，欧几里得距离）

  定义要使用的距离函数（参见 距离函数）。

此函数返回元组 (clusterid, celldata)

• clusterid:

  一个包含两列的数组，其中行数等于聚类的项目数量。每行包含分配给项目的矩形 SOM 网格中单元格的 \(x\) 和 \(y\) 坐标。
• celldata:

  一个维度为 \(\left(\verb|nxgrid|, \verb|nygrid|, \textrm{列数}\right)\) 的数组，如果对行进行聚类，则为 \(\left(\verb|nxgrid|, \verb|nygrid|, \textrm{行数}\right)\)，如果对列进行聚类。此数组的每个元素 [ix][iy] 是一个 1D 向量，包含坐标为 [ix][iy] 的网格单元格中聚类的质心的基因表达数据。

主成分分析

主成分分析 (PCA) 是一种广泛用于分析多元数据的技术。Yeung 和 Ruzzo (2001) [Yeung2001] 提供了将主成分分析应用于基因表达数据的实际示例。

本质上，PCA 是一种坐标变换，其中数据矩阵中的每一行都写成对称为主成分的基向量的线性总和，这些基向量按顺序排列，并被选择以使每个基向量最大程度地解释数据向量中剩余的方差。例如，一个 \(n \times 3\) 数据矩阵可以表示为三维空间中 \(n\) 个点的椭球云。第一个主成分是椭球的最长轴，第二个主成分是椭球的第二长轴，而第三个主成分是最短轴。数据矩阵中的每一行都可以重建为主成分的适当线性组合。然而，为了降低数据的维度，通常只保留最重要的主成分。然后，数据中存在的剩余方差被视为未解释的方差。

主成分可以通过计算数据协方差矩阵的特征向量来找到。相应的特征值决定了每个主成分解释的数据中存在的方差的多少。

在应用主成分分析之前，通常会从数据矩阵中的每一列中减去均值。在上面的例子中，这实际上将椭球云在 3D 空间中围绕其质心居中，其中主成分描述了椭球云中点相对于其质心的变化。

以下函数 pca 首先使用奇异值分解来计算数据矩阵的特征值和特征向量。奇异值分解作为 Algol 过程 svd [Golub1971] 的 C 语言翻译实现，该过程使用 Householder 双对角化和 QR 算法的一种变体。然后，以特征值大小的降序计算主成分、每个数据向量沿主成分的坐标以及对应于主成分的特征值，并返回结果。如果需要数据居中，则应在调用 pca 例程之前从数据矩阵中的每一列中减去均值。

要对矩形矩阵 data 应用主成分分析，请使用

>>> from Bio.Cluster import pca
>>> columnmean, coordinates, components, eigenvalues = pca(data)

此函数返回一个元组 columnmean, coordinates, components, eigenvalues

• columnmean

  包含 data 中每列平均值的数组。
• 坐标

  每个 data 行相对于主成分的坐标。
• 成分

  主成分。
• 特征值

  对应于每个主成分的特征值。

原始矩阵 data 可以通过计算 columnmean +  dot(coordinates, components) 来重建。

处理 Cluster/TreeView 类型文件

Cluster/TreeView 是用于对基因表达数据进行聚类的基于 GUI 的代码。它们最初由 迈克尔·艾森 在斯坦福大学编写 [Eisen1998]。 Bio.Cluster 包含用于读取和写入对应于 Cluster/TreeView 指定格式的数据文件的函数。特别是，通过将聚类结果保存为该格式，TreeView 可用于可视化聚类结果。我们建议使用 Alok Saldanha 的 http://jtreeview.sourceforge.net/ Java TreeView 程序 [Saldanha2004]，它可以显示层次结构以及 \(k\) 均值聚类结果。

Record 类的对象包含存储在 Cluster/TreeView 类型数据文件中的所有信息。要将数据文件中包含的信息存储在 Record 对象中，我们首先打开文件，然后读取它

>>> from Bio import Cluster
>>> with open("mydatafile.txt") as handle:
...     record = Cluster.read(handle)
...

此两步过程在数据源方面提供了一些灵活性。例如，您可以使用

>>> import gzip  # Python standard library
>>> handle = gzip.open("mydatafile.txt.gz", "rt")

打开 gzip 压缩文件，或

>>> from urllib.request import urlopen
>>> from io import TextIOWrapper
>>> url = "https://raw.githubusercontent.com/biopython/biopython/master/Tests/Cluster/cyano.txt"
>>> handle = TextIOWrapper(urlopen(url))

在调用 read 之前打开存储在互联网上的文件。

read 命令读取包含以 Michael Eisen 的 Cluster/TreeView 程序指定的格式存储的基因表达数据的制表符分隔文本文件 mydatafile.txt。在此文件格式中，行代表基因，列代表样本或观察值。对于一个简单的时序，一个最小的输入文件将如下所示

YORF	0 分钟	30 分钟	1 小时	2 小时	4 小时
YAL001C	1	1.3	2.4	5.8	2.4
YAL002W	0.9	0.8	0.7	0.5	0.2
YAL003W	0.8	2.1	4.2	10.1	10.1
YAL005C	1.1	1.3	0.8		0.4
YAL010C	1.2	1	1.1	4.5	8.3

每行（基因）都有一个标识符，该标识符始终位于第一列。在此示例中，我们使用酵母开放阅读框代码。每列（样本）在第一行都有一个标签。在此示例中，标签描述了采集样本的时间。第一行第一列包含一个特殊字段，用于告知程序每行中是什么类型的对象。在本例中，YORF 代表酵母开放阅读框。此字段可以是任何字母数字值。表中的剩余单元格包含适当基因和样本的数据。第 2 行第 4 列中的 5.8 表示基因 YAL001C 在 2 小时的观察值为 5.8。缺失值是可接受的，并且由空单元格指定（例如，YAL004C 在 2 小时）。

输入文件可能包含其他信息。一个最大的输入文件将如下所示

YORF	名称	基因权重	基因顺序	0	30	1	2	4
实验权重				1	1	1	1	0
实验顺序				5	3	2	1	1
YAL001C	TFIIIC 138 KD 亚基	1	1	1	1.3	2.4	5.8	2.4
YAL002W	未知	0.4	3	0.9	0.8	0.7	0.5	0.2
YAL003W	延长因子 EF1-BETA	0.4	2	0.8	2.1	4.2	10.1	10.1
YAL005C	胞质 HSP70	0.4	5	1.1	1.3	0.8		0.4

添加的列 NAME、GWEIGHT 和 GORDER 以及行 EWEIGHT 和 EORDER 是可选的。NAME 列允许您为每个基因指定一个与第 1 列中的 ID 不同的标签。

Record 对象具有以下属性

• 数据

  包含基因表达数据的数组。基因按行存储，而样本按列存储。
• 掩码

  此数组显示 data 数组中哪些元素（如果有）是缺失的。如果 mask[i, j] == 0，则 data[i, j] 缺失。如果没有找到缺失的数据，则将 mask 设置为 None。
• 基因 ID

  这是一个包含每个基因的唯一描述的列表（即 ORF 编号）。
• 基因名称

  这是一个包含每个基因的描述的列表（即基因名称）。如果数据文件中没有，则将 genename 设置为 None。
• 基因权重

  用于计算基因之间表达谱距离的权重。如果数据文件中没有，则将 gweight 设置为 None。
• 基因顺序

  基因在输出文件中存储的首选顺序。如果数据文件中没有，则将 gorder 设置为 None。
• 实验 ID

  这是一个包含每个样本的描述的列表，例如实验条件。
• 实验权重

  用于计算样本之间表达谱距离的权重。如果数据文件中没有，则将 eweight 设置为 None。
• 实验顺序

  样本在输出文件中存储的首选顺序。如果数据文件中没有，则将 eorder 设置为 None。
• 唯一 ID

  数据文件中使用的字符串代替 UNIQID。

加载 Record 对象后，可以直接访问和修改这些属性中的每一个。例如，可以通过对 record.data 取对数来对数据进行对数转换。

计算距离矩阵

要计算记录中存储的项目之间的距离矩阵，请使用

>>> matrix = record.distancematrix()

其中定义了以下参数

• transpose（默认：0）

  确定是否要计算 data 行之间的距离（transpose 为 False），还是计算 data 列之间的距离（transpose 为 True）。
• dist（默认：'e'，欧几里得距离）

  定义要使用的距离函数（参见 距离函数）。

此函数将距离矩阵作为行列表返回，其中每行的列数等于行号（见部分 计算距离矩阵）。

计算聚类中心

要计算记录中存储的项目聚类的中心，请使用

>>> cdata, cmask = record.clustercentroids()

• clusterid（默认：None）

  显示每个项目属于哪个聚类的整数向量。如果没有给出 clusterid，则假定所有项目都属于同一个聚类。
• method（默认：'a'）

  指定使用算术平均值（method=='a'）还是中位数（method=='m'）来计算聚类中心。
• transpose（默认：0）

  确定是否要计算 data 行的中心（transpose 为 False），还是计算 data 列的中心（transpose 为 True）。

此函数返回元组 cdata, cmask；有关描述，请参阅部分 计算聚类中心。

计算聚类之间的距离

要计算记录中存储的项目聚类之间的距离，请使用

>>> distance = record.clusterdistance()

其中定义了以下参数

• index1（默认值：0）

  包含属于第一个聚类的项的索引的列表。包含一个项 \(i\) 的聚类可以表示为列表 [i] 或整数 i。
• index2（默认值：0）

  包含属于第二个聚类的项目的索引的列表。仅包含一个项目 \(i\) 的聚类可以表示为列表 [i] 或整数 i。
• method（默认：'a'）

  指定如何定义聚类之间的距离

  • 'a'：两个聚类质心之间的距离（算术平均值）；

  • 'm'：两个聚类质心之间的距离（中值）；

  • 's'：两个聚类中项之间的最短成对距离；

  • 'x'：两个聚类中项之间的最长成对距离；

  • 'v'：两个聚类中项之间的成对距离的平均值。

• dist（默认：'e'，欧几里得距离）

  定义要使用的距离函数（参见 距离函数）。
• transpose（默认：0）

  如果 transpose 为 False，则计算 data 行之间的距离。如果 transpose 为 True，则计算 data 列之间的距离。

执行层次聚类

要对记录中存储的项目执行层次聚类，请使用

>>> tree = record.treecluster()

其中定义了以下参数

• transpose（默认：0）

  确定是否对行（transpose 为 False）或列（transpose 为 True）进行聚类。
• method（默认值：'m'）

  定义要使用的链接方法

  • method=='s'：成对单链接聚类

  • method=='m'：成对最大（或完全）链接聚类

  • method=='c'：成对质心链接聚类

  • method=='a'：成对平均链接聚类

• dist（默认：'e'，欧几里得距离）

  定义要使用的距离函数（参见 距离函数）。
• 转置

  确定是否要对基因或样本进行聚类。如果 transpose 为 False，则对基因（行）进行聚类。如果 transpose 为 True，则对样本（列）进行聚类。

此函数返回一个 Tree 对象。此对象包含 \(\left(\textrm{项目数量} - 1\right)\) 个节点，其中项目数量是如果对行进行聚类，则为行数，如果对列进行聚类，则为列数。每个节点描述一个成对链接事件，其中节点属性 left 和 right 分别包含一个项目或子节点的编号，而 distance 则包含它们之间的距离。项目从 0 到 \(\left(\textrm{项目数量} - 1\right)\) 编号，而聚类从 -1 到 \(-\left(\textrm{项目数量}-1\right)\) 编号。

执行 \(k\) 均值或 \(k\) 中位数聚类

要对记录中存储的项目执行 \(k\) 均值或 \(k\) 中位数聚类，请使用

>>> clusterid, error, nfound = record.kcluster()

其中定义了以下参数

• nclusters（默认值：2）

  聚类数 \(k\)。
• transpose（默认：0）

  确定要聚类的行（transpose 为 0）还是列（transpose 为 1）。
• npass（默认值：1）

  执行 \(k\) 均值/中值聚类算法的次数，每次使用不同的（随机）初始条件。如果给出了 initialid，则忽略 npass 的值，并且聚类算法只运行一次，因为它在这种情况下表现出确定性。
• method（默认值：a）

  描述如何找到聚类的中心

  • method=='a'：算术平均值（\(k\) 均值聚类）；

  • method=='m'：中值（\(k\) 中值聚类）。

  对于 method 的其他值，使用算术平均值。

• dist（默认：'e'，欧几里得距离）

  定义要使用的距离函数（参见 距离函数）。

此函数返回一个元组 (clusterid, error, nfound)，其中 clusterid 是一个包含每个行或聚类所分配的聚类编号的整数数组，error 是最佳聚类解决方案的簇内距离之和，而 nfound 是找到此最佳解决方案的次数。

计算自组织映射

要计算记录中存储的项目的自组织映射，请使用

>>> clusterid, celldata = record.somcluster()

其中定义了以下参数

• transpose（默认：0）

  确定要聚类的行（transpose 为 0）还是列（transpose 为 1）。
• nxgrid, nygrid（默认值：2, 1）

  计算自组织映射的矩形网格中水平和垂直方向上的单元格数量。
• inittau（默认值：0.02）

  SOM 算法中使用的参数 \(\tau\) 的初始值。inittau 的默认值为 0.02，这是 Michael Eisen 的 Cluster/TreeView 程序中使用的值。
• niter（默认值：1）

  要执行的迭代次数。
• dist（默认：'e'，欧几里得距离）

  定义要使用的距离函数（参见 距离函数）。

此函数返回元组 (clusterid, celldata)

• clusterid:

  一个包含两列的数组，其中行数等于聚类的项目数量。每行包含分配给项目的矩形 SOM 网格中单元格的 \(x\) 和 \(y\) 坐标。
• celldata:

  一个维度为 \(\left(\verb|nxgrid|, \verb|nygrid|, \textrm{列数}\right)\) 的数组，如果对行进行聚类，则为 \(\left(\verb|nxgrid|, \verb|nygrid|, \textrm{行数}\right)\)，如果对列进行聚类。此数组的每个元素 [ix][iy] 是一个 1D 向量，包含坐标为 [ix][iy] 的网格单元格中聚类的质心的基因表达数据。

保存聚类结果

要保存聚类结果，请使用

>>> record.save(jobname, geneclusters, expclusters)

其中定义了以下参数

• 作业名称

  字符串 jobname 用作要保存的文件的名称的基本名称。
• 基因聚类

  此参数描述了基因（按行）聚类结果。在 \(k\) 均值聚类的情况下，这是一个包含每个基因所属聚类编号的一维数组。它可以使用 kcluster 计算。在层次聚类的情况下，geneclusters 是一个 Tree 对象。
• 实验聚类

  此参数描述了实验条件的（按列）聚类结果。在 \(k\) 均值聚类的情况下，这是一个包含每个实验条件所属聚类编号的一维数组。它可以使用 kcluster 计算。在层次聚类的情况下，expclusters 是一个 Tree 对象。

此方法写入文本文件 jobname.cdt、jobname.gtr、jobname.atr、jobname*.kgg 和/或 jobname*.kag，以便随后由 Java TreeView 程序读取。如果 geneclusters 和 expclusters 都为 None，则此方法仅写入文本文件 jobname.cdt；该文件随后可以读入一个新的 Record 对象。

示例计算

这是一个分层聚类计算的示例，使用单链接聚类方法对基因进行聚类，使用最大链接聚类方法对实验条件进行聚类。由于使用欧几里得距离对基因进行聚类，因此有必要对节点距离 genetree 进行缩放，使其全部介于零和一之间。这是 Java TreeView 代码正确显示树状图所需的。为了对实验条件进行聚类，使用了未中心化的相关性。在这种情况下不需要缩放，因为 exptree 中的距离已经介于零和二之间。

示例数据 cyano.txt 可以在 Biopython 的 Tests/Cluster 子目录中找到，来自 Hihara 等人 2001 年的论文 [Hihara2001]。

>>> from Bio import Cluster
>>> with open("cyano.txt") as handle:
...     record = Cluster.read(handle)
...
>>> genetree = record.treecluster(method="s")
>>> genetree.scale()
>>> exptree = record.treecluster(dist="u", transpose=1)
>>> record.save("cyano_result", genetree, exptree)

这将创建文件 cyano_result.cdt，cyano_result.gtr 和 cyano_result.atr。

类似地，我们可以保存一个 \(k\) 均值聚类解决方案

>>> from Bio import Cluster
>>> with open("cyano.txt") as handle:
...     record = Cluster.read(handle)
...
>>> (geneclusters, error, ifound) = record.kcluster(nclusters=5, npass=1000)
>>> (expclusters, error, ifound) = record.kcluster(nclusters=2, npass=100, transpose=1)
>>> record.save("cyano_result", geneclusters, expclusters)

这将创建文件 cyano_result_K_G2_A2.cdt，cyano_result_K_G2.kgg 和 cyano_result_K_A2.kag。