계통수 만들기 — 재귀와 이진 트리로 진화 지도 그리기

이 토픽을 마치면

교과서에서 배운 이진 트리 와 재귀 를 엮어서, Newick 포맷의 계통수를 파싱하고 순회하며 clade별 분석을 수행하는 도구를 직접 만들 수 있습니다. MEGA나 iTOL 같은 실전 계통 도구가 다루는 자료구조를 코드로 이해합니다.

이 글은 교육용 일반 예제 입니다. 실제 계통수 추정은 UPGMA, Neighbor-Joining, Maximum Likelihood, Bayesian 등 여러 정교한 알고리즘을 사용합니다.

"(((A:0.1,B:0.2):0.05,C:0.3):0.1,D:0.4);" — 이게 뭐지?

여러분이 10종의 미생물 16S rRNA 서열을 분석해 계통수를 만들었습니다. 결과 파일에 이런 문자열이 있습니다.

text

(((Ecoli:0.05,Salmonella:0.06):0.02,Klebsiella:0.08):0.03,(Bacillus:0.15,Staph:0.14):0.10);

이것이 Newick 포맷 입니다. 계통수를 텍스트로 표현하는 표준. 규칙:

소괄호 () 는 하나의 clade(공통 조상에서 갈라진 그룹)
쉼표 , 는 형제 노드
콜론 뒤의 숫자 :0.05 는 부모까지의 진화 거리 (가지 길이)
세미콜론 ; 은 문자열의 끝

여러분이 하고 싶은 것:

파싱: 문자열을 파이썬 자료구조로 변환
순회: 각 clade의 잎(leaf) 종을 나열
거리 계산: 두 종 사이의 총 진화 거리 계산
시각화: 트리를 그림으로 표시

진짜 접근은 재귀입니다. Newick 문자열이 자기 자신을 자식으로 포함하는 재귀적 구조이므로, 파싱과 순회 모두 재귀로 자연스럽게 표현됩니다.

블랙박스에서 부품으로

부품 1: 트리 노드 정의

python

from dataclasses import dataclass, field
from typing import Optional


@dataclass
class TreeNode:
    name: Optional[str] = None
    branch_length: float = 0.0
    children: list["TreeNode"] = field(default_factory=list)
    
    @property
    def is_leaf(self) -> bool:
        return len(self.children) == 0

핵심 관찰: children 이 같은 타입 TreeNode 의 리스트입니다. 이 자기 참조 구조가 재귀의 자연스러운 무대입니다.

부품 2: Newick 파서 (재귀)

파서를 손으로 짜면 다음 재귀 구조가 나옵니다.

python

class NewickParser:
    def __init__(self, s: str) -> None:
        self.s = s.rstrip(";").strip()
        self.pos = 0
    
    def parse(self) -> TreeNode:
        return self._parse_node()
    
    def _parse_node(self) -> TreeNode:
        node = TreeNode()
        
        if self._peek() == "(":
            self._consume("(")
            node.children.append(self._parse_node())
            while self._peek() == ",":
                self._consume(",")
                node.children.append(self._parse_node())
            self._consume(")")
        
        node.name = self._read_name()
        
        if self._peek() == ":":
            self._consume(":")
            node.branch_length = self._read_number()
        
        return node
    
    def _peek(self) -> Optional[str]:
        return self.s[self.pos] if self.pos < len(self.s) else None
    
    def _consume(self, expected: str) -> None:
        assert self._peek() == expected, f"Expected {expected} at pos {self.pos}"
        self.pos += 1
    
    def _read_name(self) -> str:
        start = self.pos
        while self.pos < len(self.s) and self.s[self.pos] not in ",():;":
            self.pos += 1
        return self.s[start:self.pos]
    
    def _read_number(self) -> float:
        start = self.pos
        while self.pos < len(self.s) and self.s[self.pos] not in ",():;":
            self.pos += 1
        return float(self.s[start:self.pos])

핵심은 _parse_node 가 자기 자신을 재귀 호출한다는 점입니다. Newick의 재귀 문법이 재귀 함수로 그대로 매핑됩니다.

사용:

python

tree = NewickParser(
    "(((Ecoli:0.05,Salmonella:0.06):0.02,Klebsiella:0.08):0.03,(Bacillus:0.15,Staph:0.14):0.10);"
).parse()

부품 3: 재귀 순회

이제 트리를 다양한 방식으로 순회합니다. 모두 재귀입니다.

모든 잎 나열:

python

def get_leaves(node: TreeNode) -> list[str]:
    if node.is_leaf:
        return [node.name] if node.name else []
    result = []
    for child in node.children:
        result.extend(get_leaves(child))
    return result

트리 높이 계산:

python

def tree_height(node: TreeNode) -> int:
    if node.is_leaf:
        return 0
    return 1 + max(tree_height(child) for child in node.children)

들여쓰기로 표시:

python

def pretty_print(node: TreeNode, depth: int = 0) -> None:
    label = f"{node.name or '(internal)'}"
    if node.branch_length:
        label += f" [len={node.branch_length}]"
    print("  " * depth + label)
    for child in node.children:
        pretty_print(child, depth + 1)

출력:

text

(internal)
  (internal)
    (internal)
      Ecoli [len=0.05]
      Salmonella [len=0.06]
    Klebsiella [len=0.08]
  (internal)
    Bacillus [len=0.15]
    Staph [len=0.14]

부품 4: 두 잎 사이의 거리

두 종 사이의 진화 거리는 공통 조상까지의 거리 합입니다. 이것을 계산하려면 우선 각 잎의 조상 경로를 찾아야 합니다.

python

def find_path_to_leaf(node: TreeNode, target: str) -> Optional[list[TreeNode]]:
    if node.is_leaf:
        return [node] if node.name == target else None
    
    for child in node.children:
        subpath = find_path_to_leaf(child, target)
        if subpath is not None:
            return [node] + subpath
    
    return None


def evolutionary_distance(root: TreeNode, leaf_a: str, leaf_b: str) -> float:
    path_a = find_path_to_leaf(root, leaf_a)
    path_b = find_path_to_leaf(root, leaf_b)
    
    if path_a is None or path_b is None:
        raise ValueError("Leaf not found")
    
    # 공통 조상 찾기 (경로의 마지막 공통 노드)
    lca_idx = 0
    while lca_idx < min(len(path_a), len(path_b)) and path_a[lca_idx] is path_b[lca_idx]:
        lca_idx += 1
    lca_idx -= 1
    
    # 각 경로에서 LCA 이후 노드들의 branch_length 합
    distance = sum(node.branch_length for node in path_a[lca_idx + 1:])
    distance += sum(node.branch_length for node in path_b[lca_idx + 1:])
    
    return distance


d = evolutionary_distance(tree, "Ecoli", "Bacillus")
print(f"Ecoli ↔ Bacillus: {d}")
# 0.02 + 0.05 (Ecoli 경로) + 0.10 + 0.15 (Bacillus 경로) + 0.03 (LCA) = 0.35

페이딩 — 여러분이 채워야 할 두 빈칸

빈칸 1: clade 통계

특정 clade(내부 노드) 아래의 잎 수와 평균 가지 길이를 계산합니다.

python

def clade_stats(node: TreeNode) -> dict:
    """
    반환: {
        "leaf_count": 아래 잎의 수,
        "total_branch_length": 아래 모든 가지 길이 합,
        "mean_branch_length": 평균
    }
    """
    if node.is_leaf:
        # TODO: 잎 노드일 때의 기본값
        pass
    
    # TODO: 재귀 호출로 각 자식의 통계를 얻고 통합
    pass

힌트: leaf일 때 {"leaf_count": 1, "total_branch_length": node.branch_length, ...}. 내부 노드는 자식들의 결과 합산.

빈칸 2: 잎 이름 필터로 서브트리 추출

관심 있는 잎들만 남겨 축약된 트리를 만듭니다.

python

def prune_tree(node: TreeNode, keep_leaves: set[str]) -> Optional[TreeNode]:
    """
    keep_leaves에 있는 잎들만 포함하는 축약 트리 반환.
    필요 없는 내부 노드는 제거하되 가지 길이는 병합.
    """
    if node.is_leaf:
        # TODO: 잎이 keep_leaves에 있으면 반환, 아니면 None
        pass
    
    # TODO: 자식들을 재귀 pruning, None이 아닌 것만 유지
    # 자식이 0개면 None, 1개면 그 자식을 반환 (branch_length 합침)
    pass

힌트: 자식이 1개만 남으면 이 내부 노드는 불필요, 자식의 branch_length 에 이 노드의 branch_length 를 더해 반환.

성찰 — 실전 계통수 도구와의 차이

계통수 추정 알고리즘: 여러분이 다룬 것은 이미 만들어진 트리를 파싱하고 분석하는 것입니다. 실제 계통수를 서열에서 처음 만드는 것은 별개 문제입니다. UPGMA(가장 단순), Neighbor-Joining(중간 복잡도), Maximum Likelihood (RAxML, IQ-TREE), Bayesian (MrBayes, BEAST) 등이 있습니다.

비이진 트리: Newick은 이진 트리가 아닐 수 있습니다 — 한 노드가 3개 이상 자식을 가질 수 있습니다(polytomy). 여러분의 파서는 이미 이 경우를 처리합니다.

부트스트랩 값: 실전 트리는 각 내부 노드에 부트스트랩 지지도(0~100)를 표시합니다. Newick에서는 대개 내부 노드 이름 위치에 씁니다 — 여러분의 _read_name 이 이것을 잡습니다.

시각화: 실전은 ETE Toolkit, Biopython Phylo, iTOL로 그립니다. matplotlib 로 간단한 방사형·직교 트리를 그리는 것도 가능합니다.

주석 확장: 실전은 Nexus 포맷(계통수 + 서열 + 메타데이터)이나 phyloXML(XML 기반)을 씁니다. 여러분의 Newick 파서는 최소 형식.

확장 프로젝트

1. 시각화: matplotlib 로 여러분의 트리를 직교 dendrogram 형태로 그리기. 잎 이름과 가지 길이 표시.

2. 부트스트랩 필터링: 부트스트랩 지지도가 낮은 내부 노드를 pruning으로 접기.

3. UPGMA 구현: 거리 행렬에서 UPGMA로 트리를 처음부터 만들기. 여러분이 만든 TreeNode 자료구조 재사용.

4. Biopython 통합: Biopython Phylo 모듈로 파일 입출력하고, 여러분의 분석 함수 적용.

이 편의 부품 지도

[F] 이진 탐색 트리(확장): TreeNode 자기 참조 구조. 자식이 여러 개일 수 있는 일반 트리.
[F] 재귀: 파싱, 순회, 거리 계산 모두 재귀로 구현. 재귀-반복 대비.
[W] 파일 I/O: Newick 파일 읽기 등 (완성 스크립트로 제공).

[F] = 여러분이 직접 구현 / [W] = 완성 코드로 제공.