集群分析

集群分析是根據項目彼此之間的相似性將項目分組到集群中。在生物資訊學中，集群廣泛應用於基因表現資料分析，以尋找具有相似基因表現模式的基因群組。這可以識別功能相關的基因，並提示目前未知基因的功能。

Biopython 模組 Bio.Cluster 提供了常用的集群演算法，其設計初衷是應用於基因表現資料。然而，此模組也可用於其他類型資料的集群分析。Bio.Cluster 和底層的 C 集群函式庫由 De Hoon et al. 描述 [DeHoon2004]。

以下四種集群方法在 Bio.Cluster 中實作

• 階層式集群（成對質心、單一、完全和平均連結）；

• \(k\)-均值、\(k\)-中位數和 \(k\)-中心點集群；

• 自組織映射；

• 主成分分析。

資料表示

要進行集群的資料由 \(n \times m\) Numerical Python 陣列 data 表示。在基因表現資料集群的背景下，通常列對應於不同的基因，而行對應於不同的實驗條件。Bio.Cluster 中的集群演算法既可應用於列（基因）又可應用於行（實驗）。

遺失值

\(n \times m\) Numerical Python 整數陣列 mask 表示 data 中是否有任何遺失值。如果 mask[i, j] == 0，則 data[i, j] 遺失，並且在分析中會被忽略。

隨機數產生器

\(k\)-均值/中位數/中心點集群演算法和自組織映射 (SOM) 包含隨機數產生器的使用。Bio.Cluster 中的均勻隨機數產生器基於 L’Ecuyer 的演算法 [Lecuyer1988]，而服從二項分布的隨機數是使用 Kachitvichyanukul 和 Schmeiser 的 BTPE 演算法產生的 [Kachitvichyanukul1988]。隨機數產生器會在第一次呼叫期間自動初始化。由於此隨機數產生器使用兩個乘法線性同餘產生器的組合，因此需要兩個（整數）種子進行初始化，我們為此使用系統提供的隨機數產生器 rand (在 C 標準函式庫中)。我們通過使用以秒為單位的 epoch 時間呼叫 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\) 由 \(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'，歐幾里得距離）

  定義要使用的距離函數（請參閱距離函數）。

為了節省記憶體，距離矩陣以 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}\]

計算群集屬性

計算群集質心

群集的質心可以定義為所有群集項目在每個維度上的平均值或中位數。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)。質心數據儲存在 2D Numerical Python 陣列 cdata 中，遺失的數據由 2D Numerical Python 整數陣列 cmask 指示。如果 transpose 為 0，則這些陣列的維度為 \(\left(\textrm{群集數}, \textrm{列數}\right)\)，如果 transpose 為 1，則這些陣列的維度為 \(\left(\textrm{列數}, \textrm{群集數}\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\)-中心點集群

這些演算法的不同之處在於如何定義群集中心。在 \(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\)-中心點群集

kmedoids 常式使用距離矩陣和使用者傳遞的群集數量，對一組給定的項目執行 \(k\)-中心點群集

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

其中定義了以下引數：, nclusters=2, npass=1, initialid=None)|

• 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\)-中心點群集演算法的次數，每次都使用不同的（隨機）初始條件。如果給定 initialid，則會忽略 npass 的值，因為在這種情況下群集演算法的行為是確定性的。
• initialid （預設值：None）

  指定要用於 EM 演算法的初始群集。如果 initialid 為 None，則 EM 演算法的每次 npass 執行都會使用不同的隨機初始群集。如果 initialid 不為 None，則它應該等於一個 1D 陣列，其中包含每個項目的群集編號（介於 0 和 nclusters-1 之間）。每個群集應至少包含一個項目。指定初始群集後，EM 演算法是確定性的。

此函式傳回一個元組 (clusterid, error, nfound)，其中 clusterid 是一個陣列，其中包含分配給每個項目的群集編號，error 是最佳 \(k\)-中心點群集解決方案的群集內距離總和，nfound 是找到最佳解決方案的次數。請注意，clusterid 中的群集編號定義為代表群集質心的項目編號。

階層式群集

階層式分群方法與 \(k\)-平均分群方法本質上有所不同。在階層式分群中，基因或實驗條件之間表現譜的相似性以樹狀結構表示。這種樹狀結構可透過 Treeview 和 Java Treeview 等程式以圖形方式顯示，這促使階層式分群在基因表現數據分析中廣為採用。

階層式分群的第一步是計算距離矩陣，指定要分群的項目之間的所有距離。接下來，我們透過將兩個最接近的項目連接起來建立一個節點。後續節點透過成對連接項目或節點來建立，連接的依據是它們之間的距離，直到所有項目都屬於同一個節點。然後，可以透過追溯哪些項目和節點被合併來建立樹狀結構。與 \(k\)-平均分群中使用的 EM 演算法不同，階層式分群的完整過程是確定性的。

階層式分群有多種變體，它們的差異在於如何根據子節點的成員定義子節點之間的距離。在 Bio.Cluster 中，可以使用成對單一、最大、平均和質心連動。

• 在成對單一連動分群中，兩個節點之間的距離定義為兩個節點成員之間成對距離的最短距離。

• 在成對最大連動分群（也稱為成對完整連動分群）中，兩個節點之間的距離定義為兩個節點成員之間成對距離的最長距離。

• 在成對平均連動分群中，兩個節點之間的距離定義為兩個節點項目之間所有成對距離的平均值。

• 在成對質心連動分群中，兩個節點之間的距離定義為它們質心之間的距離。質心是透過計算叢集中所有項目的平均值來計算的。由於在每個步驟都需要計算每個新形成的節點與現有節點和項目之間的距離，因此成對質心連動分群的計算時間可能比其他階層式分群方法長得多。另一個特點是（對於基於皮爾森相關性的距離度量），當在分群樹中向上移動時，距離不一定會增加，甚至可能減少。這是由於使用皮爾森相關性時，質心計算和距離計算之間的不一致造成的：皮爾森相關性在距離計算時有效地將數據正規化，但在質心計算時則不會發生這種正規化。

對於成對單一、完整和平均連動分群，可以直接從各個項目之間的距離找到兩個節點之間的距離。因此，一旦知道距離矩陣，分群演算法就不需要存取原始基因表現數據。然而，對於成對質心連動分群，新形成的子節點的質心只能從原始數據計算，而不能從距離矩陣計算。

成對單一連動階層式分群的實現基於 SLINK 演算法 [Sibson1973]，它比成對單一連動分群的直接實現更快、更節省記憶體。此演算法產生的分群結果與傳統單一連動演算法找到的分群解相同。此函式庫中實現的單一連動階層式分群演算法可用於叢集大型基因表現數據集，而傳統階層式分群演算法由於過多的記憶體需求和運行時間而失敗。

表示階層式分群解

階層式分群的結果由一個節點樹組成，其中每個節點連接兩個項目或子節點。通常，我們不僅對每個節點連接哪些項目或子節點感興趣，還對它們連接時的相似性（或距離）感興趣。為了在階層式分群樹中儲存一個節點，我們使用在 Bio.Cluster 中定義的 Node 類別。 Node 的實例具有三個屬性

• left

• right

• distance

這裡，left 和 right 是指在此節點連接的兩個項目或子節點的整數，而 distance 是它們之間的距離。要分群的項目編號從 0 到 \(\left(\textrm{項目數} - 1\right)\)，而叢集的編號從 -1 到 \(-\left(\textrm{項目數}-1\right)\)。請注意，節點數比項目數少 1。

要建立新的 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 表示完整的階層式分群解。 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

可以使用方括號存取 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 等視覺化程式顯示階層式分群解，最好縮放所有節點距離，使其介於零和一之間。這可以透過在現有的 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)

其中定義了以下參數

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

在此情況下，定義了以下引數

• 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 時，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]。

自組織映射將項目組織成位於某些拓撲結構中的叢集。通常選擇矩形拓撲。自組織映射產生的叢集，其在拓撲結構中相鄰的叢集彼此之間的相似性，高於拓撲結構中彼此相距較遠的叢集。

計算自組織映射的第一步，是將一個資料向量隨機指派給拓撲結構中的每個叢集。如果要對列進行分群，則每個資料向量中的元素數量等於欄數。

然後，透過一次取一列，並找出拓撲結構中哪個叢集的資料向量最接近，來產生自組織映射。該叢集的資料向量，以及相鄰叢集的資料向量，會使用正在考慮的列的資料向量進行調整。調整由下式給定

\[\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 實作完整的演算法，以計算矩形網格上的自組織映射。首先，它會初始化亂數產生器。然後，使用亂數產生器初始化節點資料。用於修改自組織映射的基因或樣本的順序也是隨機的。自組織映射演算法中的迭代總數由使用者指定。

要執行 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)

  在自組織映射演算法中使用的參數 \(\tau\) 的初始值。 inittau 的預設值為 0.02，此值用於 Michael Eisen 的 Cluster/TreeView 程式。
• niter (預設值：1)

  要執行的迭代次數。
• dist（預設值：'e'，歐幾里得距離）

  定義要使用的距離函數（請參閱距離函數）。

此函數會傳回元組 (clusterid, celldata)

• clusterid:

  具有兩欄的陣列，其中列數等於已分群的項目數。每列都包含矩形自組織映射網格中，項目被指派到的格子的 \(x\) 和 \(y\) 座標。
• celldata:

  如果對列進行分群，則為維度 \(\left(\verb|nxgrid|, \verb|nygrid|, \textrm{欄數}\right)\) 的陣列；如果對欄進行分群，則為維度 \(\left(\verb|nxgrid|, \verb|nygrid|, \textrm{列數}\right)\) 的陣列。此陣列的每個元素 [ix][iy] 都是一個一維向量，其中包含座標為 [ix][iy] 的網格格子的叢集質心的基因表現資料。

主成分分析

主成分分析 (PCA) 是一種廣泛用於分析多變數資料的技術。 Yeung 和 Ruzzo (2001) [Yeung2001] 提供了一個將主成分分析應用於基因表現資料的實例。

本質上，主成分分析是一種座標轉換，其中資料矩陣中的每一列都寫成基向量（稱為主成分）的線性總和，這些基向量經過排序和選擇，以便每個基向量都能最大程度地解釋資料向量中的剩餘變異數。例如，一個 \(n \times 3\) 資料矩陣可以表示為三維空間中 \(n\) 個點的橢圓體雲。第一個主成分是橢圓體的最長軸，第二個主成分是橢圓體的第二長軸，第三個主成分是最短軸。資料矩陣中的每一列都可以重建為主成分的適當線性組合。但是，為了降低資料的維度，通常只保留最重要的主成分。資料中存在的剩餘變異數則被視為未解釋的變異數。

主成分可以透過計算資料共變異數矩陣的特徵向量來找到。對應的特徵值決定每個主成分解釋資料中存在的多少變異數。

在應用主成分分析之前，通常會從資料矩陣中的每欄減去平均值。在上面的範例中，這有效地將橢圓體雲圍繞其在 3D 空間中的質心居中，主成分描述了橢圓體雲中點相對於其質心的變異。

以下 pca 函式首先使用奇異值分解來計算資料矩陣的特徵值和特徵向量。奇異值分解是以 C 語言翻譯 Algol 程序 svd [Golub1971] 來實現，該程序使用 Householder 雙對角化和 QR 演算法的變體。然後，將主成分、每個資料向量沿主成分的座標，以及對應於主成分的特徵值進行評估，並以特徵值大小遞減的順序返回。如果需要資料中心化，則在呼叫 pca 常式之前，應從資料矩陣的每一列中減去平均值。

若要對矩形矩陣 data 應用主成分分析，請使用

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

此函式會回傳一個元組 columnmean, coordinates, components, eigenvalues

• columnmean

  包含 data 中每一列平均值的陣列。
• coordinates

  在主成分方面，data 中每一列的座標。
• components

  主成分。
• eigenvalues

  對應於每個主成分的特徵值。

原始矩陣 data 可以透過計算 columnmean +  dot(coordinates, components) 來重建。

處理 Cluster/TreeView 類型的檔案

Cluster/TreeView 是用於叢集基因表現資料的基於 GUI 的程式碼。它們最初由 Michael Eisen 在史丹佛大學時撰寫 [Eisen1998]。Bio.Cluster 包含用於讀取和寫入對應於 Cluster/TreeView 指定格式的資料檔案的函式。特別是，透過以該格式儲存叢集結果，可以使用 TreeView 來視覺化叢集結果。我們建議使用 Alok Saldanha 的 http://jtreeview.sourceforge.net/Java TreeView 程式 [Saldanha2004]，它可以顯示階層式以及 \(k\)-means 叢集結果。

類別 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 命令會讀取以 Tab 分隔的文字檔案 mydatafile.txt，其中包含以 Michael Eisen 的 Cluster/TreeView 程式所指定的格式呈現的基因表現資料。在這種檔案格式中，列代表基因，而欄代表樣本或觀察值。對於簡單的時間序列，最小的輸入檔案看起來會像這樣

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	NAME	GWEIGHT	GORDER	0	30	1	2	4
EWEIGHT				1	1	1	1	0
EORDER				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

  此陣列顯示 data 陣列中哪些元素遺失（若有的話）。如果 mask[i, j] == 0，則表示 data[i, j] 遺失。如果未發現遺失任何資料，則 mask 會設定為 None。
• geneid

  這是一個包含每個基因唯一描述的清單（即 ORF 編號）。
• 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）

  確定是要計算 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

  判斷正在叢集的是基因還是樣本。如果 transpose 為 False，則表示正在叢集基因（列）。如果 transpose 為 True，則表示叢集樣本（欄）。

此函數會傳回一個 Tree 物件。此物件包含 \(\left(\textrm{項目數量} - 1\right)\) 個節點，其中項目數量為如果對列進行分群則為列數，如果對欄進行分群則為欄數。每個節點描述一個成對連結事件，其中節點屬性 left 和 right 各自包含一個項目或子節點的編號，而 distance 則為它們之間的距離。項目的編號從 0 到 \(\left(\textrm{項目數量} - 1\right)\)，而叢集的編號從 -1 到 \(-\left(\textrm{項目數量}-1\right)\)。

執行 \(k\)-means 或 \(k\)-medians 叢集

若要對記錄中儲存的項目執行 \(k\)-means 或 \(k\)-medians 叢集，請使用

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

  在自組織映射演算法中使用的參數 \(\tau\) 的初始值。 inittau 的預設值為 0.02，此值用於 Michael Eisen 的 Cluster/TreeView 程式。
• niter (預設值：1)

  要執行的迭代次數。
• dist（預設值：'e'，歐幾里得距離）

  定義要使用的距離函數（請參閱距離函數）。

此函數會傳回元組 (clusterid, celldata)

• clusterid:

  具有兩欄的陣列，其中列數等於已分群的項目數。每列都包含矩形自組織映射網格中，項目被指派到的格子的 \(x\) 和 \(y\) 座標。
• celldata:

  如果對列進行分群，則為維度 \(\left(\verb|nxgrid|, \verb|nygrid|, \textrm{欄數}\right)\) 的陣列；如果對欄進行分群，則為維度 \(\left(\verb|nxgrid|, \verb|nygrid|, \textrm{列數}\right)\) 的陣列。此陣列的每個元素 [ix][iy] 都是一個一維向量，其中包含座標為 [ix][iy] 的網格格子的叢集質心的基因表現資料。

儲存叢集結果

若要儲存叢集結果，請使用

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

其中定義了以下參數

• jobname

  字串 jobname 用作要儲存檔案名稱的基底名稱。
• geneclusters

  此引數描述基因（以列為單位）叢集結果。若是 \(k\)-means 叢集，這是一個 1D 陣列，其中包含每個基因所屬的叢集編號。它可以使用 kcluster 計算。若是階層式叢集，geneclusters 是一個 Tree 物件。
• expclusters

  此引數描述實驗條件的（以欄為單位）叢集結果。若是 \(k\)-means 叢集，這是一個 1D 陣列，其中包含每個實驗條件所屬的叢集編號。它可以使用 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 *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\)-means 分群解。

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