我有一组基于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表示完成的,并且我一直在将其用于绘图目的.
解决方法:
我了解了如何从树表示形式生成链接矩阵,感谢@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和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)