6

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

返回一个 (n-1) x 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 表示完成，并且我一直将其用于绘图目的。

4

2 回答 2

8

我知道了如何从树表示中生成链接矩阵，感谢@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 矩阵实际上并没有我在这里展示的字符串，我们将使用索引方案，因为我们首先将 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)
于 2015-06-24T20:35:21.737 回答
1

我找到了这个解决方案：

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']
于 2021-03-25T14:28:43.377 回答