Newick 树表示为 scipy.cluster.hierarchy 链接矩阵格式
Newick tree representation to scipy.cluster.hierarchy linkage matrix format
我有一组基因，这些基因已根据 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 文档来看，他们对这个链接矩阵的结构的描述相当混乱。 "iteration" 是什么意思？此外，这种表示如何跟踪哪些原始观测值在哪个集群中？
我想弄清楚如何进行这种转换，因为我的项目中其他聚类分析的结果已经使用 scipy 表示完成，并且我一直将其用于绘图目的。
我知道了如何从树表示中生成链接矩阵，感谢@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 格式中，我们得到了每片叶子和树枝之间的距离，所以从 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)
我找到了这个解决方案：
import numpy as np
import pandas as pd
from ete3 import ClusterTree
from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage
import logging
def newick_to_linkage(newick: str, label_order: [str] = None) -> (np.ndarray, [str]):
"""
Convert newick tree into scipy linkage matrix
:param newick: newick string, e.g. '(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);'
:param label_order: list of labels, e.g. ['A', 'B', 'C']
:returns: linkage matrix and list of labels
"""
# newick string -> cophenetic_matrix
tree = ClusterTree(newick)
cophenetic_matrix, newick_labels = tree.cophenetic_matrix()
cophenetic_matrix = pd.DataFrame(cophenetic_matrix, columns=newick_labels, index=newick_labels)
if label_order is not None:
# sanity checks
missing_labels = set(label_order).difference(set(newick_labels))
superfluous_labels = set(newick_labels).difference(set(label_order))
assert len(missing_labels) == 0, f'Some labels are not in the newick string: {missing_labels}'
if len(superfluous_labels) > 0:
logging.warning(f'Newick string contains unused labels: {superfluous_labels}')
# reorder the cophenetic_matrix
cophenetic_matrix = cophenetic_matrix.reindex(index=label_order, columns=label_order)
# reduce square distance matrix to condensed distance matrices
pairwise_distances = pdist(cophenetic_matrix)
# return linkage matrix and labels
return linkage(pairwise_distances), list(cophenetic_matrix.columns)
基本用法：
>>> linkage_matrix, labels = newick_to_linkage(
... newick='(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);'
... )
>>> print(linkage_matrix)
[[0. 1. 0.4472136 2. ]
[2. 3. 1. 2. ]
[4. 5. 1.4832397 4. ]]
>>> print(labels)
['A', 'B', 'C', 'D']
同字矩阵是什么样的：
>>> print(cophenetic_matrix)
A B C D
A 0.0 0.3 0.9 1.0
B 0.3 0.0 1.0 1.1
C 0.9 1.0 0.0 0.7
D 1.0 1.1 0.7 0.0
高级用法：
>>> linkage_matrix, labels = newick_to_linkage(
... newick='(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);',
... label_order=['C', 'B', 'A']
... )
WARNING:root:Newick string contains unused labels: {'D'}
>>> print(linkage_matrix)
[[1. 2. 0.43588989 2. ]
[0. 3. 1.4525839 3. ]]
>>> print(labels)
['C', 'B', 'A']
我有一组基因，这些基因已根据 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 文档来看，他们对这个链接矩阵的结构的描述相当混乱。 "iteration" 是什么意思？此外，这种表示如何跟踪哪些原始观测值在哪个集群中？
我想弄清楚如何进行这种转换，因为我的项目中其他聚类分析的结果已经使用 scipy 表示完成，并且我一直将其用于绘图目的。
我知道了如何从树表示中生成链接矩阵，感谢@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 格式中，我们得到了每片叶子和树枝之间的距离，所以从 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)
我找到了这个解决方案：
import numpy as np
import pandas as pd
from ete3 import ClusterTree
from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage
import logging
def newick_to_linkage(newick: str, label_order: [str] = None) -> (np.ndarray, [str]):
"""
Convert newick tree into scipy linkage matrix
:param newick: newick string, e.g. '(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);'
:param label_order: list of labels, e.g. ['A', 'B', 'C']
:returns: linkage matrix and list of labels
"""
# newick string -> cophenetic_matrix
tree = ClusterTree(newick)
cophenetic_matrix, newick_labels = tree.cophenetic_matrix()
cophenetic_matrix = pd.DataFrame(cophenetic_matrix, columns=newick_labels, index=newick_labels)
if label_order is not None:
# sanity checks
missing_labels = set(label_order).difference(set(newick_labels))
superfluous_labels = set(newick_labels).difference(set(label_order))
assert len(missing_labels) == 0, f'Some labels are not in the newick string: {missing_labels}'
if len(superfluous_labels) > 0:
logging.warning(f'Newick string contains unused labels: {superfluous_labels}')
# reorder the cophenetic_matrix
cophenetic_matrix = cophenetic_matrix.reindex(index=label_order, columns=label_order)
# reduce square distance matrix to condensed distance matrices
pairwise_distances = pdist(cophenetic_matrix)
# return linkage matrix and labels
return linkage(pairwise_distances), list(cophenetic_matrix.columns)
基本用法：
>>> linkage_matrix, labels = newick_to_linkage(
... newick='(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);'
... )
>>> print(linkage_matrix)
[[0. 1. 0.4472136 2. ]
[2. 3. 1. 2. ]
[4. 5. 1.4832397 4. ]]
>>> print(labels)
['A', 'B', 'C', 'D']
同字矩阵是什么样的：
>>> print(cophenetic_matrix)
A B C D
A 0.0 0.3 0.9 1.0
B 0.3 0.0 1.0 1.1
C 0.9 1.0 0.0 0.7
D 1.0 1.1 0.7 0.0
高级用法：
>>> linkage_matrix, labels = newick_to_linkage(
... newick='(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);',
... label_order=['C', 'B', 'A']
... )
WARNING:root:Newick string contains unused labels: {'D'}
>>> print(linkage_matrix)
[[1. 2. 0.43588989 2. ]
[0. 3. 1.4525839 3. ]]
>>> print(labels)
['C', 'B', 'A']