python - Newick 树表示为 scipy.cluster.hierarchy 链接矩阵格式

Question

我有一组基因，它们已根据 DNA 序列进行比对和聚类，并且我在 Newick 树表示中有这组基因 (https://en.wikipedia.org/wiki/Newick_format)。有谁知道如何将这种格式转换为 scipy.cluster.hierarchy.linkage 矩阵格式？来自链接矩阵的 scipy 文档:

A (n-1) by 4 matrix Z is returned. At the i-th iteration, clusters with indices Z[i, 0] and Z[i, 1] are combined to form cluster n+i. A cluster with an index less than n corresponds to one of the n original observations. The distance between clusters Z[i, 0] and Z[i, 1] is given by Z[i, 2]. The fourth value Z[i, 3] represents the number of original observations in the newly formed cluster.

至少从 scipy 文档来看，他们对这个链接矩阵的结构的描述相当困惑。他们所说的“迭代”是什么意思？此外，这种表示如何跟踪哪些原始观测值位于哪个集群中？

我想弄清楚如何进行这种转换，因为我项目中的其他聚类分析的结果已经用 scipy 表示法完成，并且我一直将其用于绘图目的。

Answer 1

我知道了如何从树表示中生成链接矩阵，感谢@cel 的澄清。让我们以 Newick 维基页面 (https://en.wikipedia.org/wiki/Newick_format) 中的示例为例

字符串格式的树是:

(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);

首先，应该计算所有叶子之间的距离。例如，我们要计算A和B的距离，方法是通过最近的分支从A到B遍历树。因为在 Newick 格式中，我们得到了每片叶子和 Twig 之间的距离，所以从 A 到 B 的距离很简单0.1 + 0.2 = 0.3。对于 A 到 D，我们必须做 0.1 + (0.5 + 0.4) = 1.0，因为从 D 到最近的分支的距离给定为 0.4，而从 D 的分支到 A 的距离是0.5。因此距离矩阵看起来像这样(索引 A=0，B=1，C=2，D=3):

distance_matrix=
 [[0.0, 0.3, 0.9, 1.0],
  [0.3, 0.0, 1.0, 1.1],
  [0.9, 1.0, 0.0, 0.7],
  [1.0, 1.1, 0.1, 0.0]]

从这里，链接矩阵很容易找到。因为我们已经有 n=4 个集群 (A,B,C,D) 作为原始观察结果，我们需要找到树的额外 n-1 簇。每一步只是简单地将两个聚类组合成一个新聚类，我们取彼此最接近的两个聚类。在这种情况下，A 和 B 离得最近，因此链接矩阵的第一行将如下所示:

[A,B,0.3,2]

从现在开始，我们将 A 和 B 视为一个集群，其到最近分支的距离是 A 和 B 之间的距离。

现在我们剩下 3 个簇，AB、C 和 D。我们可以更新距离矩阵以查看哪些集群距离最近。设 AB 在更新后的距离矩阵中的索引为 0。

distance_matrix=
[[0.0, 1.1, 1.2],
 [1.1, 0.0, 0.7],
 [1.2, 0.7, 0.0]]

我们现在可以看到 C 和 D 彼此最接近，所以让我们将它们组合成一个新的集群。链接矩阵中的第二行现在将是

[C,D,0.7,2]

现在，我们只剩下两个集群，AB 和 CD。这些簇到根分支的距离分别为 0.3 和 0.7，因此它们的距离为 1.0。链接矩阵的最后一行将是:

[AB,CD,1.0,4]

现在，scipy 矩阵实际上不会像我在此处显示的那样具有适当的字符串，我们将使用索引方案，因为我们首先组合了 A 和 B，AB 将有索引 4 和 CD 会有索引 5。所以我们应该在 scipy 链接矩阵中看到的实际结果是:

[[0,1,0.3,2],
 [2,3,0.7,2],
 [4,5,1.0,4]]

这是从树表示到 scipy 链接矩阵表示的一般方法。但是，已经存在其他 python 包中的工具可以读取 Newick 格式的树，从这些工具中，我们可以相当容易地找到距离矩阵，然后将其传递给 scipy 的链接函数。下面是一个小脚本，它完全适用于此示例。

from ete2 import ClusterTree, TreeStyle
import scipy.cluster.hierarchy as sch
import scipy.spatial.distance
import matplotlib.pyplot as plt
import numpy as np
from itertools import combinations


tree = ClusterTree('(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);')
leaves = tree.get_leaf_names()
ts = TreeStyle()
ts.show_leaf_name=True
ts.show_branch_length=True
ts.show_branch_support=True

idx_dict = {'A':0,'B':1,'C':2,'D':3}
idx_labels = [idx_dict.keys()[idx_dict.values().index(i)] for i in range(0, len(idx_dict))]

#just going through the construction in my head, this is what we should get in the end
my_link = [[0,1,0.3,2],
        [2,3,0.7,2],
        [4,5,1.0,4]]

my_link = np.array(my_link)


dmat = np.zeros((4,4))

for l1,l2 in combinations(leaves,2):
    d = tree.get_distance(l1,l2)
    dmat[idx_dict[l1],idx_dict[l2]] = dmat[idx_dict[l2],idx_dict[l1]] = d

print 'Distance:'
print dmat


schlink = sch.linkage(scipy.spatial.distance.squareform(dmat),method='average',metric='euclidean')

print 'Linkage from scipy:'
print schlink

print 'My link:'
print my_link

print 'Did it right?: ', schlink == my_link

dendro = sch.dendrogram(my_link,labels=idx_labels)
plt.show()

tree.show(tree_style=ts)