生命的家谱

((A,B),(C,D))

((A,C),(B,D))

((A,D),(B,C))

A、B | C、D

物种A  A C G T A C G T
物种B  A C G T G C G T
物种C  G C A T A C A T
物种D  G C A T G C A T

A

A → G → C

C

A → G → A

    A
     \
      ●
     /
    B

物种A  A
物种B  A
物种C  G
物种D  G

A A G G

A G A G

A G G A

A A A A

AB | CD

from itertools import combinations


def neighbor_joining(names, matrix):
    """
    根据距离矩阵建立无根Neighbor-Joining树。

    names:
        叶节点名称列表，例如 ["A", "B", "C", "D"]

    matrix:
        对称距离矩阵，对角线为0

    返回：
        graph:
            无根树邻接表
            graph[node] = [(neighbor, branch_length), ...]

        steps:
            每次合并的记录

    说明：
        不对负分支长度进行静默截断。
        若输入不是加性距离，可能出现负长度估计。
    """

    if len(names) < 3:
        raise ValueError("至少需要三个物种")

    if len(matrix) != len(names):
        raise ValueError("距离矩阵维度错误")

    if len(set(names)) != len(names):
        raise ValueError("叶节点名称必须唯一")

    n = len(names)

    for i in range(n):
        if len(matrix[i]) != n:
            raise ValueError("距离矩阵不是方阵")

        if abs(matrix[i][i]) > 1e-9:
            raise ValueError("对角线必须为0")

        for j in range(n):
            if matrix[i][j] < 0:
                raise ValueError("输入距离不能为负")

            if abs(matrix[i][j] - matrix[j][i]) > 1e-9:
                raise ValueError("距离矩阵必须对称")

    d = {
        a: {
            b: float(matrix[i][j])
            for j, b in enumerate(names)
        }
        for i, a in enumerate(names)
    }

    active = list(names)

    graph = {
        name: [] for name in names
    }

    steps = []

    def link(a, b, length):
        graph.setdefault(a, []).append((b, length))
        graph.setdefault(b, []).append((a, length))

    internal_id = 0

    while len(active) > 3:

        size = len(active)

        row_sum = {
            a: sum(d[a][k] for k in active)
            for a in active
        }

        # 选择Q值最小的一对
        def q_value(a, b):
            return (
                (size - 2) * d[a][b]
                - row_sum[a]
                - row_sum[b]
            )

        a, b = min(
            combinations(active, 2),
            key=lambda pair: (
                q_value(pair[0], pair[1]),
                pair
            )
        )

        distance_ab = d[a][b]

        branch_a = (
            distance_ab / 2
            + (row_sum[a] - row_sum[b])
            / (2 * (size - 2))
        )

        branch_b = distance_ab - branch_a

        new_node = f"U{internal_id}"
        internal_id += 1

        link(new_node, a, branch_a)
        link(new_node, b, branch_b)

        steps.append({
            "join": (a, b),
            "q": q_value(a, b),
            "new_node": new_node,
            "branches": (branch_a, branch_b)
        })

        # 为新节点计算到其余对象的距离
        d[new_node] = {}

        for k in active:
            if k in (a, b):
                continue

            new_distance = (
                d[a][k] + d[b][k] - distance_ab
            ) / 2

            d[new_node][k] = new_distance
            d[k][new_node] = new_distance

        # 删除已经合并的两个对象
        for k in active:
            d[k].pop(a, None)
            d[k].pop(b, None)

        d.pop(a)
        d.pop(b)

        active = [
            k for k in active
            if k not in (a, b)
        ]

        active.append(new_node)

    # 最后只剩三个对象
    a, b, c = active

    center = "Ufinal"

    for x, y, z in (
        (a, b, c),
        (b, a, c),
        (c, a, b)
    ):
        length = (
            d[x][y] + d[x][z] - d[y][z]
        ) / 2

        link(center, x, length)

    return graph, steps

NAMES = ["A", "B", "C", "D"]

DISTANCES = [
    [0, 5, 9, 9],
    [5, 0, 10, 10],
    [9, 10, 0, 8],
    [9, 10, 8, 0],
]

graph, steps = neighbor_joining(
    NAMES,
    DISTANCES
)

print("合并过程:")

for step in steps:
    print(step)

print()
print("无根树的全部边:")

visited = set()

for node, neighbors in graph.items():
    for other, length in neighbors:

        edge = tuple(sorted((node, other)))

        if edge in visited:
            continue

        visited.add(edge)

        print(
            node,
            "<->",
            other,
            "长度:",
            round(length, 6)
        )

A      <-> U0       长度 2
B      <-> U0       长度 3
U0     <-> Ufinal   长度 3
C      <-> Ufinal   长度 4
D      <-> Ufinal   长度 4

def tree_distance(graph, start, target):
    """
    在已经构造好的树中，
    计算两个节点之间的唯一简单路径长度。
    """

    stack = [(start, None, 0.0)]

    while stack:
        node, parent, distance = stack.pop()

        if node == target:
            return distance

        for neighbor, length in graph[node]:
            if neighbor == parent:
                continue

            stack.append((
                neighbor,
                node,
                distance + length
            ))

    raise ValueError("图中找不到目标节点")


for i, a in enumerate(NAMES):
    for j in range(i + 1, len(NAMES)):
        b = NAMES[j]

        reconstructed = tree_distance(
            graph, a, b
        )

        expected = DISTANCES[i][j]

        assert abs(
            reconstructed - expected
        ) < 1e-9

        print(
            a, b,
            "原距离:", expected,
            "树上距离:", reconstructed
        )

print("全部六项距离恢复正确")

((A:2,B:3):3,C:4,D:4);

def to_newick(graph, root="Ufinal"):
    """
    从指定的技术根开始输出Newick。

    注意：
        root只是序列化起点。
        不代表已经确定生物学上的真正祖先根。
    """

    def visit(node, parent):
        children = [
            (other, length)
            for other, length in graph[node]
            if other != parent
        ]

        if not children:
            return node

        parts = []

        for child, length in sorted(children):
            subtree = visit(child, node)

            parts.append(
                f"{subtree}:{length:.6f}"
            )

        return "(" + ",".join(parts) + ")"

    return visit(root, None) + ";"


print(to_newick(graph))

A A G G

A G A G

A G G A

A A A A

import math

BASES = "ACGT"


def jc69_probability(a, b, branch_length):
    """
    JC69模型：
    在给定分支长度下，
    祖先碱基a变成后代碱基b的概率。
    """

    factor = math.exp(
        -4.0 * branch_length / 3.0
    )

    if a == b:
        return 0.25 + 0.75 * factor

    return 0.25 - 0.25 * factor


def partial_likelihood(node, column, sequences):
    """
    在树上自底向上计算条件似然。

    叶节点：
        用物种名字符串表示。

    内部节点：
        [(child, branch_length), ...]

    返回：
        {A:概率, C:概率, G:概率, T:概率}
    """

    if isinstance(node, str):
        observed = sequences[node][column]

        return {
            base: float(base == observed)
            for base in BASES
        }

    child_results = [
        (
            partial_likelihood(
                child, column, sequences
            ),
            branch_length
        )
        for child, branch_length in node
    ]

    result = {}

    for ancestral_base in BASES:

        probability = 1.0

        for child_scores, length in child_results:

            subtotal = sum(
                jc69_probability(
                    ancestral_base,
                    descendant_base,
                    length
                ) * child_scores[descendant_base]

                for descendant_base in BASES
            )

            probability *= subtotal

        result[ancestral_base] = probability

    return result


def tree_log_likelihood(
    sequences,
    left_pair,
    right_pair,
):
    """
    对指定四叶拓扑计算观测对数似然。

    固定分支长度：
        末端边0.08
        内部总边长0.20

    这里只比较固定参数下的树。
    尚未进行分支长度优化。
    """

    def clade(pair):
        return [
            (pair[0], 0.08),
            (pair[1], 0.08),
        ]

    # 将内部边在计算表示上拆成两段0.10
    root = [
        (clade(left_pair), 0.10),
        (clade(right_pair), 0.10),
    ]

    length = len(next(iter(sequences.values())))

    total_loglik = 0.0

    for column in range(length):

        root_scores = partial_likelihood(
            root, column, sequences
        )

        # JC69平衡频率：A/C/G/T各为0.25
        site_probability = sum(
            root_scores.values()
        ) / 4.0

        total_loglik += math.log(
            site_probability
        )

    return total_loglik

patterns = (
    [("A", "A", "G", "G")] * 12
    + [("A", "G", "A", "G")] * 5
    + [("A", "G", "G", "A")] * 3
    + [("A", "A", "A", "A")] * 10
)

sequences = {
    taxon: "".join(
        column[i]
        for column in patterns
    )
    for i, taxon in enumerate("ABCD")
}


topologies = [
    (("A", "B"), ("C", "D")),
    (("A", "C"), ("B", "D")),
    (("A", "D"), ("B", "C")),
]

results = []

for left, right in topologies:
    score = tree_log_likelihood(
        sequences, left, right
    )

    label = (
        "".join(left)
        + " | "
        + "".join(right)
    )

    results.append((score, label))

    print(
        label,
        "log-likelihood:",
        round(score, 4)
    )


best = max(results)

print()
print("固定参数下得分最高的拓扑:", best[1])

这是五千万年前真正发生过的全部祖先分化。