使用ete3与曼哈顿距离字典实现邻接法遇TreeError问题求助

问题背景 #

基于曼哈顿距离矩阵转换的字典，使用ete3的Tree数据结构实现邻接法（Neighbor Joining），涉及84个样本个体。实现代码如下：

while len(t.get_leaves()) > 2:
    min_dist = float('inf')
    for node1 in t:
        for node2 in t:
            if node1 != node2:
                dist = sum(distance_dict[(leaf1.name, leaf2.name)] for leaf1 in node1.get_leaves() for leaf2 in node2.get_leaves()) / (len(node1.get_leaves()) * len(node2.get_leaves()))
                if dist < min_dist and dist > 0:
                    min_dist = dist
                    min_node1 = node1
                    min_node2 = node2
    new_node = t.add_child()
    new_node.add_child(min_node1)
    new_node.add_child(min_node2)
    new_node.dist = min_dist / 2.0
    if min_node1.is_leaf() and min_node1 in t:
        t.remove_child(min_node1)
    if min_node2.is_leaf() and min_node2 in t:
        t.remove_child(min_node2)

运行时出现错误：

TreeError: 'child not found'

错误原因与修复方案 #

1. 核心错误：移除节点的逻辑冗余 #

当执行new_node.add_child(min_node1)时，ete3会自动将min_node1从原父节点（t）中移除——因为树结构中一个节点只能有一个父节点。此时再调用t.remove_child(min_node1)，min_node1已不属于t的子节点，自然会抛出"child not found"错误。

修复：直接删除最后两行的remove_child代码，无需手动移除节点。

2. 遍历逻辑错误：错误遍历所有节点 #

for node1 in t会遍历树中所有节点（包括内部节点和叶节点），但邻接法每次应合并的是当前树的顶层聚类单元（即t的直接子节点），否则会重复计算内部节点的组合，导致逻辑混乱。

修复：将遍历对象改为t.children，仅处理当前树的直接子节点：

for node1 in t.children:
    for node2 in t.children:

3. 额外问题：距离计算不符合邻接法逻辑 #

当前代码的距离计算是两个聚类单元的平均距离，这是UPGMA算法的逻辑，而非邻接法（NJ）的校正距离公式。邻接法需使用校正后的距离：
[
D_{ij} = d_{ij} - \frac{r_i + r_j}{N-2}
]
其中(r_i = \sum_{k=1}^N d_{ik})，(N)为当前聚类单元的数量。若目标是实现标准邻接法，需修正距离计算逻辑。

修复后的完整代码示例 #

while len(t.get_leaves()) > 2:
    min_dist = float('inf')
    # 仅遍历当前树的直接子节点（聚类单元）
    cluster_nodes = list(t.children)
    n_clusters = len(cluster_nodes)
    if n_clusters < 2:
        break
    
    # 计算每个聚类单元的r值（用于NJ距离校正）
    r_values = {}
    for node in cluster_nodes:
        total = 0
        node_leaves = node.get_leaves()
        for other_node in cluster_nodes:
            if node == other_node:
                continue
            other_leaves = other_node.get_leaves()
            dist_sum = sum(distance_dict[(l1.name, l2.name)] for l1 in node_leaves for l2 in other_leaves)
            avg_dist = dist_sum / (len(node_leaves) * len(other_leaves))
            total += avg_dist
        r_values[node] = total
    
    # 查找校正后距离最小的聚类对
    avg_dist = 0
    for i in range(n_clusters):
        node1 = cluster_nodes[i]
        for j in range(i+1, n_clusters):
            node2 = cluster_nodes[j]
            # 计算原始平均距离
            node1_leaves = node1.get_leaves()
            node2_leaves = node2.get_leaves()
            dist_sum = sum(distance_dict[(l1.name, l2.name)] for l1 in node1_leaves for l2 in node2_leaves)
            avg_dist = dist_sum / (len(node1_leaves) * len(node2_leaves))
            # NJ校正距离
            nj_dist = avg_dist - (r_values[node1] + r_values[node2]) / (n_clusters - 2)
            if nj_dist < min_dist and nj_dist > 0:
                min_dist = nj_dist
                min_node1 = node1
                min_node2 = node2
    
    # 创建新节点并设置子节点距离（NJ公式）
    new_node = t.add_child()
    node1_size = len(min_node1.get_leaves())
    node2_size = len(min_node2.get_leaves())
    dist_to_node1 = (avg_dist / 2) + (r_values[min_node1] - r_values[min_node2]) / (2 * (n_clusters - 2))
    dist_to_node2 = avg_dist - dist_to_node1
    
    min_node1.dist = dist_to_node1
    min_node2.dist = dist_to_node2
    
    new_node.add_child(min_node1)
    new_node.add_child(min_node2)

内容的提问来源于stack exchange，提问作者Silastander