NJ·UPGMA 계통수 완전 이해 — 거리 행렬에서 이진 트리로 가는 두 그리디 알고리즘

서열에서 진화의 나무로

MSA를 만들었으면 다음 단계는 계통수(phylogenetic tree) 입니다. 종·유전자 사이의 진화적 관계를 이진 트리로 표현. 이 트리 하나로 종 분화 순서, 공통 조상 추정, 진화 속도 측정 등 무수한 분석이 파생됩니다.

계통수 알고리즘은 크게 세 갈래로 나뉩니다.

거리 기반 (이 편, M26): 서열 사이의 거리 행렬로 트리 구축. 빠르고 직관적
파시모니 (M27): 관측 데이터에 가장 적은 변이 수를 요구하는 트리
최대우도 (M28): 관측 데이터의 확률을 최대화하는 트리

이 편에서는 거리 기반의 두 대표 UPGMA와 NJ를 손으로 유도하고 파이썬으로 구현합니다.

거리 행렬 만들기

MSA에서 서열 쌍 사이의 거리를 계산합니다. 가장 단순한 것은 p-distance: 정렬된 컬럼 중 문자가 다른 컬럼의 비율.

$p_{ij} = \frac{\text{mismatch 컬럼 수}}{\text{총 컬럼 수}}$

실무에서는 Jukes-Cantor 보정 등을 적용한 진화 거리로 확장합니다. 여러 번의 치환이 있었을 가능성을 통계적으로 보정.

$d_{JC}(p) = -\frac{3}{4} \ln\left(1 - \frac{4p}{3}\right)$

이 편의 유도에는 편의상 p-distance 그대로 씁니다. 4개 종의 예제.

	A	B	C	D
A	0	5	9	9
B	5	0	10	10
C	9	10	0	8
D	9	10	8	0

UPGMA — 가장 단순한 접근

UPGMA(Unweighted Pair Group Method with Arithmetic mean)은 계층적 클러스터링(average linkage)의 특수 사례입니다.

알고리즘:

거리 행렬에서 가장 작은 거리 쌍을 찾는다
그 쌍을 하나의 클러스터로 병합
새 클러스터와 나머지 노드 사이 거리를 계산 (기존 두 노드 거리의 평균)
노드가 하나 남을 때까지 반복

위 예제로 손 유도.

Step 1: 최소 거리 = d(A,B) = 5. A와 B 병합. 병합 노드까지의 branch length = 5/2 = 2.5.

Step 2: 새 클러스터 (AB)와 나머지 거리 재계산.

$d((AB), C) = \frac{d(A,C) + d(B,C)}{2} = \frac{9 + 10}{2} = 9.5$

$d((AB), D) = \frac{9 + 10}{2} = 9.5$

새 거리 행렬:

	(AB)	C	D
(AB)	0	9.5	9.5
C	9.5	0	8
D	9.5	8	0

Step 3: 최소 거리 = d(C,D) = 8. C와 D 병합. Branch length = 8/2 = 4.

Step 4: (CD)와 (AB) 병합.

$d((AB), (CD)) = \frac{d((AB),C) + d((AB),D)}{2} = \frac{9.5 + 9.5}{2} = 9.5$

Branch length = 9.5/2 - 4 = 0.75 (C·D 서브트리의 깊이 4를 빼줌).

결과 트리 (branch length 포함):

text

┌── A (2.5)
     ┌─┤
(AB) │ └── B (2.5)
─────┤
(CD) │ ┌── C (4)
     └─┤
       └── D (4)

UPGMA의 결정적 가정 — 분자 시계

UPGMA는 분자 시계(molecular clock) 을 가정합니다. 모든 branch가 root부터 같은 시간을 흘렀다 = 모든 leaf가 root에서 같은 거리에 있다(ultrametric tree).

이 가정은 자주 깨집니다. 어떤 계통은 진화 속도가 더 빠르고, 어떤 계통은 느립니다. 분자 시계가 성립 안 하면 UPGMA는 트리 위상 자체를 틀리게 그립니다.

NJ — 분자 시계를 놓아준 접근

Saitou & Nei(MBE 1987)의 Neighbor Joining은 분자 시계 가정을 놓아버립니다. Branch length가 leaf마다 다를 수 있음. 대신 병합할 쌍을 고르는 기준이 더 정교합니다.

Q 행렬의 정의가 NJ의 심장입니다.

$Q(i, j) = (n - 2) \cdot d(i, j) - \sum_{k} d(i, k) - \sum_{k} d(j, k)$

n은 남은 노드 수, 뒤 두 항은 각 노드의 총 거리 합. Q 값이 가장 작은 쌍을 병합합니다.

왜 Q 행렬인가? 두 노드 i, j가 실제로 진화적으로 가깝다면 두 노드 각각의 "나머지 세계와의 총 거리"가 비슷할 것이고, 그 값들이 Q에서 상쇄되어 순수한 pairwise 거리에 집중하게 됩니다. UPGMA가 단순히 최소 거리를 골랐다면, NJ는 주변 상황을 고려한 최적 병합을 고릅니다.

같은 예제로 손 유도.

Q 계산:

Σ_k d(A,k) = 5 + 9 + 9 = 23
Σ_k d(B,k) = 5 + 10 + 10 = 25
Σ_k d(C,k) = 9 + 10 + 8 = 27
Σ_k d(D,k) = 9 + 10 + 8 = 27

n = 4. 그러면.

$Q(A,B) = 2 \times 5 - 23 - 25 = 10 - 48 = -38$ $Q(A,C) = 2 \times 9 - 23 - 27 = 18 - 50 = -32$ $Q(A,D) = 2 \times 9 - 23 - 27 = 18 - 50 = -32$ $Q(B,C) = 2 \times 10 - 25 - 27 = 20 - 52 = -32$ $Q(B,D) = 2 \times 10 - 25 - 27 = 20 - 52 = -32$ $Q(C,D) = 2 \times 8 - 27 - 27 = 16 - 54 = -38$

최소 Q = -38 (동률: A-B와 C-D). 첫 병합은 A-B 또는 C-D. NJ에서는 이런 동률이 실전에서 자주 발생하고, 결과 트리 위상이 미묘하게 달라질 수 있습니다.

Branch length 계산은 더 정교합니다.

$L(A) = \frac{d(A,B)}{2} + \frac{\sum_k d(A,k) - \sum_k d(B,k)}{2(n-2)} = \frac{5}{2} + \frac{23 - 25}{4} = 2.5 - 0.5 = 2$ $L(B) = d(A,B) - L(A) = 5 - 2 = 3$

UPGMA에서는 두 branch가 대칭(2.5, 2.5)이었지만, NJ에서는 비대칭(2, 3). B가 A보다 살짝 빠르게 진화한 것으로 해석.

파이썬 NJ 구현 (~40줄)

python

def neighbor_joining(D, labels):
    """NJ 알고리즘. D는 거리 행렬(리스트의 리스트), labels는 노드 이름."""
    tree = []
    n = len(labels)
    nodes = list(range(n))
    dist = [row[:] for row in D]
    next_id = n

    while len(nodes) > 2:
        m = len(nodes)
        # Q 행렬 계산
        total = [sum(dist[i][k] for k in range(m)) for i in range(m)]
        Q = [[0]*m for _ in range(m)]
        for i in range(m):
            for j in range(m):
                if i != j:
                    Q[i][j] = (m-2) * dist[i][j] - total[i] - total[j]

        # 최소 Q 쌍
        min_val, mi, mj = float('inf'), 0, 1
        for i in range(m):
            for j in range(i+1, m):
                if Q[i][j] < min_val:
                    min_val, mi, mj = Q[i][j], i, j

        # branch length
        L_i = dist[mi][mj]/2 + (total[mi] - total[mj]) / (2*(m-2))
        L_j = dist[mi][mj] - L_i
        tree.append((labels[nodes[mi]], next_id, L_i))
        tree.append((labels[nodes[mj]], next_id, L_j))
        labels.append(f"N{next_id}")

        # 새 거리 행렬
        new_dist = []
        for k in range(m):
            if k in (mi, mj):
                continue
            row = []
            for l in range(m):
                if l in (mi, mj):
                    continue
                row.append(dist[k][l])
            # 새 노드까지 거리
            new_d = (dist[k][mi] + dist[k][mj] - dist[mi][mj]) / 2
            row.append(new_d)
            new_dist.append(row)
        last_row = [(dist[mi][k] + dist[mj][k] - dist[mi][mj])/2
                    for k in range(m) if k not in (mi, mj)] + [0]
        new_dist.append(last_row)

        dist = new_dist
        nodes = [nodes[k] for k in range(m) if k not in (mi, mj)] + [next_id]
        next_id += 1

    # 남은 두 노드 병합
    if len(nodes) == 2:
        tree.append((labels[nodes[0]], labels[nodes[1]], dist[0][1]))
    return tree

# 위 예제
D = [[0, 5, 9, 9],
     [5, 0, 10, 10],
     [9, 10, 0, 8],
     [9, 10, 8, 0]]
labels = ['A', 'B', 'C', 'D']
tree = neighbor_joining(D, labels[:])
for edge in tree:
    print(edge)

실행하면 위 손 유도와 같은 결과가 나옵니다. NJ는 40줄. 트리 알고리즘은 원리 자체는 어렵지 않고, 오히려 실전에서 부트스트랩·모델 선택 등 통계적 검증이 훨씬 큰 무게를 차지합니다.

UPGMA vs NJ 정리

항목	UPGMA	NJ
분자 시계	가정 O	가정 X
Branch length	leaf 대칭 (ultrametric)	leaf 비대칭
병합 기준	최소 pairwise 거리	최소 Q 값
정확도	분자 시계 성립 시만 정확	훨씬 자주 정확
계산 복잡도	O(N³)	O(N³)
오늘 실무	drug discovery의 짧은 SAR 트리	종·유전자 계통수 표준

의심 없는 default는 NJ. UPGMA는 분자 시계 가정을 명시적으로 원할 때만.

복잡도

두 알고리즘 다 O(N³). N은 노드 수. 실전 서열 수백 개까지는 초 단위. 수천～수만 개는 fast NJ 변종(RapidNJ 등)이 필요합니다.

부트스트랩 — 신뢰도 평가

계통수의 각 브랜치는 얼마나 신뢰할 수 있을까요? 부트스트랩 리샘플링이 표준 도구입니다.

MSA의 컬럼을 무작위로 반복 추출해서 새 MSA 만듬
그 MSA로 새 트리 구축
100～1000회 반복
각 브랜치가 몇 %의 부트스트랩 트리에서 지지되었는지 표시

부트스트랩 값 ≥ 70～80이면 신뢰할 만함. 미만이면 그 브랜치의 위상은 데이터가 부족합니다.

CS 매핑 — 그리디 트리 병합

DryBench의 계층적 클러스터링과 정확히 같은 뼈대입니다.

UPGMA = average linkage 계층적 클러스터링
NJ = Q 기반 가중 병합 (더 정교한 링크 함수)
두 접근 다 greedy bottom-up 병합

이 패턴은 다른 곳에도 나타납니다. 데이터 사이언스의 hierarchical clustering, k-means의 initialization으로 사용되는 Ward's method, 심지어 컴파일러의 함수 인라이닝 결정. 거리 정보에서 트리 구조를 뽑는 것이 데이터 분석의 흔한 패턴입니다.

Rosalind에서 채점받기

Rosalind NJ 문제. 거리 행렬과 노드 수가 주어지고 NJ 트리를 요구합니다. 위 코드의 tree 결과를 Newick 포맷으로 변환하면 통과.

다음 편으로 이어지는 갈래

다음 편 (M27): 파시모니 — Fitch, Sankoff — 관측 MSA에 필요한 최소 변이 수로 트리 평가
두 편 뒤 (M28): 최대우도 — IQ-TREE, RAxML — 통계 모델로 트리 우도 최대화
실전 도구: MEGA, FastTree, IQ-TREE — 오늘의 표준 계통수 도구

더 깊게 파고 싶다면

Saitou & Nei (1987), The neighbor-joining method: a new method for reconstructing phylogenetic trees, MBE — NJ 원 논문.
Felsenstein Inferring Phylogenies — 계통수 알고리즘의 표준 교재.
UC Berkeley BIDS — Yun S. Song 교수의 계통수 강의 (자막 완비).
Compau & Pevzner Bioinformatics Algorithms Chapter 7 — Rosalind 연동 NJ 튜토리얼.

MEGA 또는 FastTree로 실제 단백질 패밀리(예: 헤모글로빈) 계통수를 뽑아봅시다. NJ와 최대우도 트리를 비교하면, 이 편의 감이 훨씬 깊어집니다.

BioPlayground

계통수 거리 기반: NJ와 UPGMA 알고리즘의 손 유도와 파이썬 구현