Newick 树表示到 scipy.cluster.hierarchy 链接矩阵格式

Question

我有一组基于 DNA 序列对齐和聚类的基因，我在 Newick 树表示中拥有这组基因 ( https://en.wikipedia.org/wiki/Newick_format )。有谁知道如何将此格式转换为 scipy.cluster.hierarchy.linkage 矩阵格式？来自链接矩阵的 scipy 文档：

返回一个 (n-1) × 4 的矩阵 Z。在第 i 次迭代中，索引为 Z[i, 0] 和 Z[i, 1] 的簇被组合成簇 n+i。索引小于 n 的聚类对应于 n 个原始观测值之一。簇 Z[i, 0] 和 Z[i, 1] 之间的距离由 Z[i, 2] 给出。第四个值 Z[i, 3] 表示新形成的聚类中原始观察的数量。

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

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

Answer 1

我知道如何从树表示中生成链接矩阵，感谢@cel 的澄清。让我们以 Newick wiki 页面 ( 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 格式中，我们给出了每个叶子和分支之间的距离，从 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 矩阵实际上不会像我在这里展示的那样将字符串放在适当的位置，我们将使用索引方案，因为我们首先AB将A 和 B 组合在一起，索引为 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)