聚类分析

集群分析是根据项目彼此的相似性将项目分组到集群中。在生物信息学中,集群广泛用于基因表达数据分析,以寻找具有相似基因表达谱的基因组。这可能会识别功能相关基因,并暗示目前未知基因的功能。

Biopython模块 Bio.Cluster 提供了常用的聚类算法,设计时考虑了基因表达数据的应用。但是,该模块也可以用于其他类型数据的聚类分析。 Bio.Cluster 底层的C集群库由De Hoon描述 et al. [DeHoon2004].

以下四种集群方法在中实现 Bio.Cluster :

  • 分层聚集(成对的重心、单一、完全和平均链接);

  • \(k\) - 意味着, \(k\) - 中位数,以及 \(k\) - 中密度聚集;

  • 自组织地图;

  • 主成分分析。

数据表示

要聚集的数据由 \(n \times m\) Python数组 data .在基因表达数据集群的背景下,通常行对应于不同的基因,而列对应于不同的实验条件。中的集群算法 Bio.Cluster 可以应用于行(基因)和列(实验)。

缺失值

\(n \times m\) 数字Python整数数组 mask 指示中是否有任何值 data 不见了如果 mask[i, j] == 0 那么 data[i, j] 在分析中缺失并被忽略。

随机数生成器

\(k\) -means/medoids集群算法和自组织地图(SOM)包括使用随机数生成器。中的均匀随机数生成器 Bio.Cluster 基于L ' Ecuyer的算法 [Lecuyer1988], 而遵循二项分布的随机数是由Kachitvichyanukul和Schmeiser使用BTPE算法生成的 [Kachitvichyanukul1988]. 随机数生成器在第一次调用期间自动初始化。由于此随机数生成器使用两个乘线性全等生成器的组合,因此需要两个(整)种子进行初始化,为此我们使用系统提供的随机数生成器 rand (in C标准库)。我们通过调用初始化此生成器 srand 历元时间以秒为单位,并使用由 rand 作为中均匀随机数生成器的种子 Bio.Cluster .

距离函数

为了根据项目的相似性将其分组,我们应该首先定义我们的确切含义 similar . Bio.Cluster 提供八个距离函数,由单个字符指示,来测量相似性,或相反,测量距离:

  • 'e' 欧氏距离;

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

  • 'c' Pearson相关系数;

  • 'a' :皮尔逊相关系数的绝对值;

  • 'u' :非中心皮尔森相关性(相当于两个数据载体之间的角度的cos);

  • '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},\]

因此被称为 metrics .在日常语言中,这意味着两点之间的最短距离是一条直线。

其余六个距离测量与相关系数有关,其中距离 \(d\) 是根据相关性来定义的 \(r\) 通过 \(d=1-r\) .请注意,这些距离函数是 semi-metrics 不满足三角不等式。例如用于

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

we find a Pearson distance \(d\left(\underline{u},\underline{w}\right) = 1.8660\), while \(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|.\]

这等于你必须在城市中两点之间步行的距离,你必须沿着城市街区步行。至于欧几里得距离,表达数据是直接相互减去的,因此我们应该确保它们被正确规范化。

Pearson相关系数

皮尔逊相关系数定义为

\[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),\]

in which \(\bar{x}, \bar{y}\) are the sample mean of \(x\) and \(y\) respectively, and \(\sigma_x, \sigma_y\) are the sample standard deviation of \(x\) and \(y\). The Pearson correlation coefficient is a measure for how well a straight line can be fitted to a scatterplot of \(x\) and \(y\). If all the points in the scatterplot lie on a straight line, the Pearson correlation coefficient is either +1 or -1, depending on whether the slope of line is positive or negative. If the Pearson correlation coefficient is equal to zero, there is no correlation between \(x\) and \(y\).

Pearson distance 然后被定义为

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

由于Pearson相关系数介于-1和1之间,因此Pearson距离介于0和2之间。

绝对皮尔森相关性

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

相应的距离定义为

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

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

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

非中心相关性(角度的cos)

在某些情况下,最好使用 uncentered correlation 而不是常规的皮尔逊相关系数。非中心相关性定义为

\[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}\]

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

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

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

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

非中心相关性等于中的两个数据载体的角度的cos \(n\) - 维度空间,通常被称为维度空间。

绝对非中心相关

对于常规Pearson相关性,我们可以使用非中心相关性的绝对值定义距离测量:

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

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

从几何上讲,非中心相关性的绝对值等于两个数据载体的支持线之间的cos(即,角度,而不考虑载体的方向)。

Spearman秩相关

Spearman等级相关性是非参数相似性衡量的一个例子,并且对于异常值往往比Pearson相关性更稳健。

为了计算Spearman等级相关性,如果我们按照每个载体中的数据的值对它们进行排序,我们就会用它们的等级替换每个数据值。然后我们计算两个等级载体(而不是数据载体)之间的Pearson相关性。

与Pearson相关性的情况一样,我们可以将对应于Spearman等级相关性的距离测量定义为

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

哪里 \(r_{\mbox{S}}\) 是斯皮尔曼等级相关性。

肯德尔的 \(\tau\)

肯德尔的 \(\tau\) 是非参数相似性度量的另一个示例。它类似于Spearman等级相关性,但仅使用相对等级来计算,而不是等级本身 \(\tau\) (see Snedecor & Cochran [Snedecor1989]) .

我们可以定义与肯德尔对应的距离测量 \(\tau\) 作为

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

作为肯德尔的 \(\tau\) 始终在-1和1之间,相应的距离将在0和2之间。

加权

对于中可用的大多数距离功能 Bio.Cluster ,可以应用权重载体。权重载体包含数据载体中项目的权重。如果物品重量 \(i\)\(w_i\) ,那么该项目将被视为发生了 \(w_i\) 数据中的时间。权重不一定是整数。

计算距离矩阵

距离矩阵是一个平方矩阵,其中的项目之间的所有成对距离 data ,并可以通过函数计算 distancematrixBio.Cluster 模块:

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

其中定义了以下参数:

  • data (必需)
    包含项数据的数组。
  • mask (默认: None )
    显示缺失哪些数据的一个整数数组。如果 mask[i, j] == 0 那么 data[i, j] 失踪了如果 maskNone ,那么所有数据都存在。
  • weight (默认: None )
    计算距离时使用的权重。如果 weightNone ,那么假设相等的权重。
  • transpose (默认: 0 )
    Determines if the distances between the rows of data are to be calculated (transpose is False), or between the columns of data (transpose is True).
  • dist (默认: 'e' 、欧几里得距离)
    定义要使用的距离函数(请参阅 距离函数 ).

为了节省内存,距离矩阵作为1D阵列列表返回。各行中的列数等于行号。因此,第一行有零个元素。例如,

>>> 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}\]

计算集群属性

计算集群重心

集群的重心可以定义为所有集群项上每个维度的平均值或中位数。功能 clustercentroidsBio.Cluster 可用于计算:

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

其中定义了以下参数:

  • data (必需)
    包含项数据的数组。
  • mask (默认: None )
    显示缺失哪些数据的一个整数数组。如果 mask[i, j] == 0 那么 data[i, j] 失踪了如果 maskNone ,那么所有数据都存在。
  • clusterid (默认: None )
    显示每个项目所属的集群的整值载体。如果 clusteridNone ,则假设所有项目属于同一集群。
  • method (默认: 'a' )
    指定算术平均值是否 (method=='a' )或中位数 (method=='m' )用于计算集群中心。
  • transpose (默认: 0 )
    Determines if the centroids of the rows of data are to be calculated (transpose is False), or the centroids of the columns of data (transpose is True).

This function returns the tuple (cdata, cmask). The centroid data are stored in the 2D Numerical Python array cdata, with missing data indicated by the 2D Numerical Python integer array cmask. The dimensions of these arrays are \(\left(\textrm{number of clusters}, \textrm{number of columns}\right)\) if transpose is 0, or \(\left(\textrm{number of rows}, \textrm{number of clusters}\right)\) if transpose is 1. Each row (if transpose is 0) or column (if transpose is 1) contains the averaged data corresponding to the centroid of each cluster.

计算集群之间的距离

给定之间的距离函数 items ,我们可以定义两个之间的距离 clusters 以多种方式。两个聚类的算术平均值之间的距离用于成对质心连锁聚类, \(k\) - 意味着集群。在 \(k\) -medoids集群,而是使用两个集群的中位数之间的距离。两个集群的项目之间的最短成对距离用于成对单连锁集群,而最长成对距离用于成对最大连锁集群。在成对平均连锁集群中,两个集群之间的距离被定义为成对距离的平均值。

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

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

其中定义了以下参数:

  • data (必需)
    包含项数据的数组。
  • mask (默认: None )
    显示缺失哪些数据的一个整数数组。如果 mask[i, j] == 0 那么 data[i, j] 失踪了如果 maskNone ,那么所有数据都存在。
  • weight (默认: None )
    计算距离时使用的权重。如果 weightNone ,那么假设相等的权重。
  • index1 (默认: 0 )
    包含属于第一个集群的项目的索引的列表。仅包含一个项目的集群 \(i\) 可以表示为列表 [i] ,或作为一个整体 i .
  • index2 (默认: 0 )
    包含属于第二个集群的项目的索引的列表。仅包含一个项的集群 \(i\) 可以表示为列表 [i] ,或作为一个整体 i .
  • method (默认: 'a' )
    指定如何定义集群之间的距离:
    • 'a' :两个集群中心之间的距离(算术平均值);

    • 'm' :两个集群重心之间的距离(中位数);

    • 's' :两个聚类中的项目之间的最短成对距离;

    • 'x' :两个集群中项目之间的最长成对距离;

    • 'v' :两个集群中项目之间成对距离的平均值。

  • dist (默认: 'e' 、欧几里得距离)
    定义要使用的距离函数(请参阅 距离函数 ).
  • transpose (默认: 0 )
    如果 transposeFalse ,计算行之间的距离 data .如果 transposeTrue ,计算柱之间的距离 data .

分区算法

分区算法将项目分为 \(k\) 集群,使项目到其集群中心的距离总和最小。簇的数量 \(k\) 由用户指定。有三种分区算法 Bio.Cluster :

  • \(k\) - 意味着集群

  • \(k\) - 中位数聚集

  • \(k\) -medoids集群

这些算法在定义集群中心的方式上有所不同。在 \(k\) - 意味着集群,集群中心被定义为集群中所有项目的平均数据载体。而不是卑鄙,在 \(k\) - 中位数对数据载体中的每个维度计算中位数进行聚集。最后在 \(k\) 对集群中心进行集群的-medoids被定义为与集群中其他项目的距离总和最小的项目。该集群算法适用于距离矩阵已知但原始数据矩阵不可用的情况,例如基于蛋白质的结构相似性对蛋白质进行集群时。

期望最大化(EM)算法用于找到该划分 \(k\) 组在EM算法的初始化中,我们随机将项目分配到聚类。为了确保不产生空簇,我们使用二项分布来随机选择每个簇中的项目数为一个或多个。然后,我们随机地将聚类分配置换到项目,使得每个项目在任何聚类中具有相等的概率。因此,每个集群都保证包含至少一个项目。

然后我们卸载:

  • 计算每个集群的重心,定义为集群的平均值、中位数或中间值;

  • 计算每个物品到集群中心的距离;

  • 对于每个项目,确定哪个集群重心最接近;

  • 将每个项重新分配到其最近的集群,或者如果没有进一步重新分配项,则停止迭代。

为了避免集群在迭代期间变空,在 \(k\) - 意味着和 \(k\) -medians集群该算法会跟踪每个集群中的项目数量,并禁止将集群中最后剩余的项目重新分配到不同的集群。为 \(k\) -medoids集群,不需要这样的检查,因为充当集群重心的项目与自身的距离为零,因此永远不会更接近不同的集群。

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

EM算法的运行频率取决于正在聚集的项目数量。根据经验,我们可以考虑找到最佳解决方案的频率;这个数字由该库中实现的分区算法返回。如果最佳解决方案被发现多次,则不太可能存在比找到的更好的解决方案。但是,如果最佳解决方案只被发现一次,则很可能存在群内距离总和较小的其他解决方案。如果物品数量很大(超过几百个),可能很难找到全球最优的解决方案。

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

\(k\) - 意味着和 \(k\) - 中位数

\(k\) - 意味着和 \(k\) - 中位数算法作为函数实现 kclusterBio.Cluster :

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

其中定义了以下参数:

  • data (必需)
    包含项数据的数组。
  • nclusters (默认: 2 )
    簇的数量 \(k\) .
  • mask (默认: None )
    显示缺失哪些数据的一个整数数组。如果 mask[i, j] == 0 那么 data[i, j] 失踪了如果 maskNone ,那么所有数据都存在。
  • weight (默认: None )
    计算距离时使用的权重。如果 weightNone ,那么假设相等的权重。
  • transpose (默认: 0 )
    确定行是否 (transpose0 )或列 (transpose1 )将被聚集。
  • npass (默认: 1 )
    的次数 \(k\) 执行-means/-medians集群算法,每次都使用不同的(随机)初始条件。如果 initialid 给出的值 npass 被忽略,并且集群算法只运行一次,因为在这种情况下它的行为是确定性的。
  • method (默认: a )
    描述如何找到集群的中心:
    • method=='a' :算术平均数 (\(k\) - 意味着集群);

    • method=='m' :中位数 (\(k\) - 中位数集群)。

    对于其他价值观 method ,使用算术平均值。

  • dist (默认: 'e' 、欧几里得距离)
    定义要使用的距离函数(请参阅 距离函数 ).而所有八项距离测量均由 kcluster ,从理论角度来看,最好使用欧几里得距离来计算 \(k\) - 均值算法,以及 \(k\) - 中位数。
  • initialid (默认: None )
    指定用于EM算法的初始集群。如果 initialidNone ,然后对每个 npass EM算法的运行。如果 initialidNone ,那么它应该等于包含集群号(介于 0nclusters-1 )针对每个项目。每个集群应至少包含一个项。通过指定初始集群,EM算法是确定性的。

该函数返回一个tuple (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=无)|

  • distance (必需)
    包含项目之间距离的矩阵;该矩阵可以通过三种方式指定:
    • 作为2D Numerical Python数组(其中仅访问数组的左下部分):

      distance = array([[0.0, 1.1, 2.3], [1.1, 0.0, 4.5], [2.3, 4.5, 0.0]])
      
    • 作为1D Numerical 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算法的初始集群。如果 initialidNone ,然后对每个 npass EM算法的运行。如果 initialidNone ,那么它应该等于包含集群号(介于 0nclusters-1 )针对每个项目。每个集群应至少包含一个项。通过指定初始集群,EM算法是确定性的。

该函数返回一个tuple (clusterid, error, nfound) ,在哪里 clusterid 是一个数组,包含每个项目被分配到的集群的编号, error 是最优的群内距离总和 \(k\) -medoids集群解决方案,以及 nfound 是找到最佳解决方案的次数。请注意,中的集群号 clusterid 定义为表示集群重心的项的项号。

层次聚类

分层集群方法本质上与 \(k\) - 意味着集群方法。在分层集群中,基因或实验条件之间表达谱的相似性以树结构的形式表示。这种树结构可以通过Treeview和Java Treeview等程序以图形方式显示,这有助于分层集群在基因表达数据分析中的流行。

分层集群的第一步是计算距离矩阵,指定要集群的项目之间的所有距离。接下来,我们通过连接两个最近的项来创建一个节点。后续节点是通过根据项目或节点之间的距离成对连接项目或节点来创建的,直到所有项目属于同一个节点。然后可以通过追溯合并的项和节点来创建树结构。与EM算法不同,EM算法用于 \(k\) - 意味着集群,分层集群的整个过程是确定性的。

存在多种形式的分层集群,其不同之处在于如何根据其成员定义子节点之间的距离。在 Bio.Cluster 、两两单个、最大、平均和重心联动可用。

  • 在成对单链接集群中,两个节点之间的距离被定义为两个节点成员之间成对距离中的最短距离。

  • 在成对最大链接集群(也称为成对完全链接集群)中,两个节点之间的距离被定义为两个节点成员之间的成对距离中的最长距离。

  • 在成对平均链接集群中,两个节点之间的距离被定义为两个节点的项目之间所有成对距离的平均值。

  • 在成对中心链接集群中,两个节点之间的距离被定义为其中心之间的距离。通过取集群中所有项目的平均值来计算重心。由于每一步都需要计算每个新形成的节点到现有节点和项的距离,因此成对中心链接集群的计算时间可能比其他分层集群方法明显长。另一个特点是(对于基于Pearson相关性的距离测量),距离在集群树中上升时不一定增加,甚至可能减少。这是由于使用Pearson相关性时,重心计算和距离计算之间的不一致造成的:虽然Pearson相关性有效地对距离计算的数据进行了标准化,但对于重心计算,没有发生这样的标准化。

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

成对单链接分层集群的实现基于SLINK算法 [Sibson1973], 这比成对单链接集群的直接实现要快得多,内存效率更高。该算法产生的集群结果与传统单链接算法找到的集群解相同。该库中实现的单链接分层集群算法可用于对大型基因表达数据集进行集群,而传统的分层集群算法由于内存需求和运行时间过多而失败。

代表分层集群解决方案

分层集群的结果由一棵节点树组成,其中每个节点连接两个项或子节点。通常,我们不仅对每个节点连接哪些项或子节点感兴趣,而且还对它们连接时的相似性(或距离)感兴趣。为了在分层集群树中存储一个节点,我们使用类 Node ,定义于 Bio.Cluster .的实例 Node 具有三个属性:

  • left

  • right

  • distance

Here, left and right are integers referring to the two items or subnodes that are joined at this node, and distance is the distance between them. The items being clustered are numbered from 0 to \(\left(\textrm{number of items} - 1\right)\), while clusters are numbered from -1 to \(-\left(\textrm{number of items}-1\right)\). Note that the number of nodes is one less than the number of items.

创建一个新 Node 对象,我们需要指定 leftright ; distance is optional.

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

的属性 left , right ,而且 distance 现有 Node 对象可以直接修改:

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

如果出现以下情况,则会引发错误 leftright 不是整,或者如果 distance 无法转换为浮点值。

Python类 Tree 代表了完整的分层集群解决方案。一 Tree 对象可以从列表中创建 Node 对象:

>>> 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

a中的各个节点 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 object.但是,我们可以将树转换为节点列表,修改此列表,并从此列表创建新树:

>>> 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之间。这可以通过调用 scale 现有方法 Tree 对象:

>>> tree.scale()

此方法不接受参数,并返回 None .

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

>>> indices = tree.sort(order)

这样物品 indices[i] 将发生在位置 \(i\) 在树图中。

分层集群后,项目可以分组为 \(k\) 基于存储在 Tree 通过砍树来对象:

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

哪里 nclusters (默认为 1 )是所需的集群数量 \(k\) .这种方法忽略顶部 \(k-1\) 在树结构中链接事件,导致 \(k\) 分开的物品集群。簇的数量 \(k\) 应该是正数,并且小于或等于项目数。此方法返回数组 clusterid 包含每个项目被分配到的集群的编号。集群已编号 \(0\)\(k-1\) 按照树图中从左到右的顺序排列。

执行分层集群

要执行分层集群,请使用 treecluster 功能 Bio.Cluster .

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

其中定义了以下参数:

  • data
    包含项数据的数组。
  • mask (默认: None )
    显示缺失哪些数据的一个整数数组。如果 mask[i, j] == 0 那么 data[i, j] 失踪了如果 maskNone ,那么所有数据都存在。
  • weight (默认: None )
    计算距离时使用的权重。如果 weightNone ,那么假设相等的权重。
  • transpose (默认: 0 )
    确定行是否 (transposeFalse )或列 (transposeTrue )将被聚集。
  • method (默认: 'm' )
    定义要使用的链接方法:
    • method=='s' :成对单连锁聚集

    • method=='m' :成对最大(或完全)连锁聚集

    • method=='c' :成对中心连锁聚集

    • method=='a' :成对平均连锁聚集

  • dist (默认: 'e' 、欧几里得距离)
    定义要使用的距离函数(请参阅 距离函数 ).

要对预先计算的距离矩阵应用分层集群,请指定 distancematrix 呼叫时的争论 treecluster 函数而不是 data 论点:

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

在这种情况下,定义了以下参数:

  • distancematrix
    距离矩阵,可以通过三种方式指定:
    • 作为2D Numerical Python数组(其中仅访问数组的左下部分):

      distance = array([[0.0, 1.1, 2.3], [1.1, 0.0, 4.5], [2.3, 4.5, 0.0]])
      
    • 作为1D Numerical Python数组,连续包含距离矩阵左下半部分的距离:

      distance = array([1.1, 2.3, 4.5])
      
    • 作为包含距离矩阵左下部分行的列表:

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

    这三个公式对应于相同的距离矩阵。作为 treecluster 作为集群算法的一部分,可能会对距离矩阵中的值进行洗牌,请确保在调用之前将此数组保存在不同的变量中 treecluster 如果您以后需要的话。

  • method
    要使用的联动方法:
    • method=='s' :成对单连锁聚集

    • method=='m' :成对最大(或完全)连锁聚集

    • method=='a' :成对平均连锁聚集

    虽然成对的单一连锁、最大连锁和平均连锁聚集可以单独根据距离矩阵计算,但成对的中心连锁却不能。

打电话时 treecluster 或者 datadistancematrixNone .

This function returns a Tree object. This object contains \(\left(\textrm{number of items} - 1\right)\) nodes, where the number of items is the number of rows if rows were clustered, or the number of columns if columns were clustered. Each node describes a pairwise linking event, where the node attributes left and right each contain the number of one item or subnode, and distance the distance between them. Items are numbered from 0 to \(\left(\textrm{number of items} - 1\right)\), while clusters are numbered -1 to \(-\left(\textrm{number of items}-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_{\texttrm {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] 失踪了如果 maskNone ,那么所有数据都存在。
  • weight (默认: None )
    包含计算距离时使用的权重。如果 weightNone ,那么假设相等的权重。
  • transpose (默认: 0 )
    确定行是否 (transpose0 )或列 (transpose1 )将被聚集。
  • nxgrid, nygrid (默认: 2, 1 )
    计算自组织地图的矩形网格中水平和垂直单元的数量。
  • inittau (默认: 0.02 )
    参数的初始值 \(\tau\) 这用于SOM算法。的默认值 inittau 是0.02,用于Michael Thomen的集群/TreeView程序。
  • niter (默认: 1 )
    要执行的迭代次数。
  • dist (默认: 'e' 、欧几里得距离)
    定义要使用的距离函数(请参阅 距离函数 ).

此函数返回该数组 (clusterid, celldata) :

  • clusterid :
    具有两列的数组,其中行的数量等于聚集的项的数量。每一行包含 \(x\)\(y\) 指定该项目的矩形SOM网格中单元格的坐标。
  • celldata :
    An array with dimensions \(\left(\verb|nxgrid|, \verb|nygrid|, \textrm{number of columns}\right)\) if rows are being clustered, or \(\left(\verb|nxgrid|, \verb|nygrid|, \textrm{number of rows}\right)\) if columns are being clustered. Each element [ix][iy] of this array is a 1D vector containing the gene expression data for the centroid of the cluster in the grid cell with coordinates [ix][iy].

主成分分析

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

本质上,PCA是一种坐标变换,其中数据矩阵中的每一行都被写为称为主成分的基向的线性和,这些基向是经过排序和选择的,以便每个最大限度地解释数据载体中的剩余方差。例如一种 \(n \times 3\) 数据矩阵可以表示为椭圆形云 \(n\) 三维空间中的点。第一主成分是椭圆体的最长轴,第二主成分是椭圆体的第二长轴,第三主成分是最短轴。数据矩阵中的每一行都可以重建为主成分的合适线性组合。然而,为了降低数据的维度,通常只保留最重要的主成分。然后,数据中存在的剩余方差被视为无法解释的方差。

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

在应用主成分分析之前,通常从数据矩阵中的每列中减去均值。在上面的示例中,这有效地将椭圆形云集中在3D空间中的质心周围,主成分描述椭圆形云中的点相对于其质心的变化。

功能 pca 下面首先使用奇异值分解来计算数据矩阵的特征值和特征量。奇异值分解作为Algol过程的C翻译实现 svd [Golub1971], 它使用Householder双对角化和QR算法的变体。然后,评估主成分、沿着主成分的每个数据载体的坐标以及与主成分对应的特征值,并以特征值大小的降序返回。如果需要数据居中,则应在调用 pca 套路

将主成分分析应用于矩形矩阵 data ,使用

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

该函数返回一个tuple columnmean, coordinates, components, eigenvalues :

  • columnmean
    包含中每列平均值的数组 data .
  • coordinates
    中每行的坐标 data 关于主要成分。
  • components
    主要组成部分。
  • eigenvalues
    与每个主成分对应的特征值。

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

处理集群/TreeView类型文件

集群/TreeView是用于对基因表达数据进行聚集的基于图形界面的代码。它们最初是由 Michael Eisen 在斯坦福大学的时候 [Eisen1998]. Bio.Cluster 包含用于读写与为集群/TreeView指定的格式相对应的数据文件的函数。特别是,通过以该格式保存集群结果,TreeView可以用于可视化集群结果。我们建议使用Alok SalDanha的http://jreeview.sourceforge.net/Java TreeView程序 [Saldanha2004], 它可以显示分层以及 \(k\) - 表示聚类结果。

类的对象 Record 包含存储在集群/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")

打开gzipped文件,或者

>>> 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 命令读取制表符分隔的文本文件 mydatafile.txt 包含为Michael Thomen的集群/TreeView程序指定的格式的基因表达数据。在此文件格式中,行代表基因,列代表样本或观察结果。对于简单的时间过程,最小输入文件如下所示:

YORF

0分钟

30分钟

1小时

2小时

4小时

YAL 001 C

1

1.3

2.4

5.8

2.4

YAL 002 W

0.9

0.8

0.7

0.5

0.2

YAL 003 W

0.8

2.1

4.2

10.1

10.1

YAL 005 C

1.1

1.3

0.8

0.4

YAL 010 C

1.2

1

1.1

4.5

8.3

每一行(基因)都有一个始终位于第一列的标识符。在这个例子中,我们使用酵母开放阅读框架代码。每列(示例)的第一行都有一个标签。在本示例中,标签描述了采集样本的时间。第一行的第一列包含一个特殊字段,告诉程序每一行中有什么类型的对象。在这种情况下,YORF代表酵母开放阅读框。此字段可以是任何字母数字值。表中的其余单元包含适当基因和样本的数据。第2行第4列中的5.8意味着基因YAL 001 C在2小时时的观察值为5.8。缺失值是可以接受的,并用空单元格指定(例如,2小时时的YAL 004 C)。

输入文件可能包含其他信息。最大输入文件如下:

YORF

NAME

GWEIGHT

GORDER

0

30

1

2

4

EWEIGHT

1

1

1

1

0

EORDER

5

3

2

1

1

YAL 001 C

TFIIIC 138 KD子单位

1

1

1

1.3

2.4

5.8

2.4

YAL 002 W

UNKNOWN

0.4

3

0.9

0.8

0.7

0.5

0.2

YAL 003 W

延伸系数EF 1-Beta

0.4

2

0.8

2.1

4.2

10.1

10.1

YAL 005 C

胞质热休克蛋白70

0.4

5

1.1

1.3

0.8

0.4

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

A Record 对象具有以下属性:

  • data
    包含基因表达数据的数据数组。基因按行存储,而样本按列存储。
  • mask
    该数组显示 data 数组(如果有的话)缺失。如果 mask[i, j] == 0 那么 data[i, j] 失踪了如果没有发现数据丢失, mask 设置为 None .
  • geneid
    这是一个包含每个基因的唯一描述的列表(即,开放阅读框号)。
  • genename
    这是包含每个基因描述的列表(即,基因名称)。如果数据文件中不存在, genename 设置为 None .
  • gweight
    用于计算基因之间表达谱距离的权重。如果数据文件中不存在, gweight 设置为 None .
  • gorder
    基因在输出文件中存储的首选顺序。如果数据文件中不存在, gorder 设置为 None .
  • expid
    这是一个包含每个样本描述的列表,例如实验条件。
  • eweight
    用于计算样本之间表达谱中距离的权重。如果数据文件中不存在, eweight 设置为 None .
  • eorder
    样本存储在输出文件中的首选顺序。如果数据文件中不存在, eorder 设置为 None .
  • uniqid
    在数据文件中用于代替UNIQID的字符串。

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

计算距离矩阵

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

>>> matrix = record.distancematrix()

其中定义了以下参数:

  • transpose (默认: 0 )
    Determines if the distances between the rows of data are to be calculated (transpose is False), or between the columns of data (transpose is True).
  • dist (默认: 'e' 、欧几里得距离)
    定义要使用的距离函数(请参阅 距离函数 ).

此函数以行列表的形式返回距离矩阵,其中每一行的列数等于行号(请参阅部分 计算距离矩阵 ).

计算集群重心

要计算记录中存储的项目集群的重心,请使用

>>> cdata, cmask = record.clustercentroids()
  • clusterid (默认: None )
    显示每个项目所属的集群的整值载体。如果 clusterid 没有给出,则假设所有项目属于同一集群。
  • method (默认: 'a' )
    指定算术平均值是否 (method=='a' )或中位数 (method=='m' )用于计算集群中心。
  • transpose (默认: 0 )
    Determines if the centroids of the rows of data are to be calculated (transpose is False), or the centroids of the columns of data (transpose is 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 )
    如果 transposeFalse ,计算行之间的距离 data .如果 transposeTrue ,计算柱之间的距离 data .

执行分层集群

要对记录中存储的项执行分层集群,请使用

>>> tree = record.treecluster()

其中定义了以下参数:

  • transpose (默认: 0 )
    确定行是否 (transposeFalse )或列 (transposeTrue )将被聚集。
  • method (默认: 'm' )
    定义要使用的链接方法:
    • method=='s' :成对单连锁聚集

    • method=='m' :成对最大(或完全)连锁聚集

    • method=='c' :成对中心连锁聚集

    • method=='a' :成对平均连锁聚集

  • dist (默认: 'e' 、欧几里得距离)
    定义要使用的距离函数(请参阅 距离函数 ).
  • transpose
    确定基因或样本是否正在聚集。如果 transposeFalse ,基因(行)正在聚集。如果 transposeTrue ,样本(列)被聚集。

This function returns a Tree object. This object contains \(\left(\textrm{number of items} - 1\right)\) nodes, where the number of items is the number of rows if rows were clustered, or the number of columns if columns were clustered. Each node describes a pairwise linking event, where the node attributes left and right each contain the number of one item or subnode, and distance the distance between them. Items are numbered from 0 to \(\left(\textrm{number of items} - 1\right)\), while clusters are numbered -1 to \(-\left(\textrm{number of items}-1\right)\).

执行 \(k\) - 意味着或 \(k\) - 中位数聚集

执行 \(k\) - 意味着或 \(k\) - 对存储在记录中的项目进行聚集的中位数,使用

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

其中定义了以下参数:

  • nclusters (默认: 2 )
    簇的数量 \(k\) .
  • transpose (默认: 0 )
    确定行是否 (transpose0 )或列 (transpose1 )将被聚集。
  • npass (默认: 1 )
    的次数 \(k\) 执行-means/-medians集群算法,每次都使用不同的(随机)初始条件。如果 initialid 给出的值 npass 被忽略,并且集群算法只运行一次,因为在这种情况下它的行为是确定性的。
  • method (默认: a )
    描述如何找到集群的中心:
    • method=='a' :算术平均数 (\(k\) - 意味着集群);

    • method=='m' :中位数 (\(k\) - 中位数集群)。

    对于其他价值观 method ,使用算术平均值。

  • dist (默认: 'e' 、欧几里得距离)
    定义要使用的距离函数(请参阅 距离函数 ).

该函数返回一个tuple (clusterid, error, nfound) ,在哪里 clusterid 是一个包含各行或集群被分配到的集群号的整组数组, error 是最佳集群解决方案的集群内距离总和,并且 nfound 是找到该最佳解决方案的次数。

计算自组织地图

要计算记录中存储的项目的自组织地图,请使用

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

其中定义了以下参数:

  • transpose (默认: 0 )
    确定行是否 (transpose0 )或列 (transpose1 )将被聚集。
  • nxgrid, nygrid (默认: 2, 1 )
    计算自组织地图的矩形网格中水平和垂直单元的数量。
  • inittau (默认: 0.02 )
    参数的初始值 \(\tau\) 这用于SOM算法。的默认值 inittau 是0.02,用于Michael Thomen的集群/TreeView程序。
  • niter (默认: 1 )
    要执行的迭代次数。
  • dist (默认: 'e' 、欧几里得距离)
    定义要使用的距离函数(请参阅 距离函数 ).

此函数返回该数组 (clusterid, celldata) :

  • clusterid :
    具有两列的数组,其中行的数量等于聚集的项的数量。每一行包含 \(x\)\(y\) 指定该项目的矩形SOM网格中单元格的坐标。
  • celldata :
    An array with dimensions \(\left(\verb|nxgrid|, \verb|nygrid|, \textrm{number of columns}\right)\) if rows are being clustered, or \(\left(\verb|nxgrid|, \verb|nygrid|, \textrm{number of rows}\right)\) if columns are being clustered. Each element [ix][iy] of this array is a 1D vector containing the gene expression data for the centroid of the cluster in the grid cell with coordinates [ix][iy].

保存聚类结果

要保存集群结果,请使用

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

其中定义了以下参数:

  • jobname
    字符串 jobname 用作要保存的文件名称的基本名称。
  • geneclusters
    此参数描述了基因(行)集群结果。的情况下 \(k\) - 表示聚类,这是一个包含每个基因所属的聚类数的1D数组。它可以用以下公式计算: kcluster .在分层集群的情况下, geneclusters 是一 Tree object.
  • expclusters
    该论点描述了实验条件下的(逐列)聚集结果。的情况下 \(k\) - 意味着集群,这是一个1D数组,包含每个实验条件所属的集群的编号。它可以用以下公式计算: kcluster .在分层集群的情况下, expclusters 是一 Tree object.

此方法写入文本文件 jobname.cdt , jobname.gtr , jobname.atr , jobname*.kgg ,和/或 jobname*.kag 供Java TreeView程序后续阅读。如果 geneclustersexpclusters 都是 None ,该方法只写文本文件 jobname.cdt ;该文件随后可以被读取到新的 Record object.

实例计算

这是分层集群计算的一个例子,对基因使用单连锁集群,对实验条件使用最大连锁集群。由于欧几里得距离被用于基因集群,因此有必要缩放节点距离 genetree 使得它们都在零和一之间。Java TreeView代码正确显示树图需要这一点。为了对实验条件进行聚集,使用了非中心相关性。在这种情况下不需要缩放,因为中的距离 exptree 已经在0到2之间了

示例数据 cyano.txt 可以在Biopython的 Tests/Cluster 这位是日原报社的记者 et al. 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 .