Phylogenetics (계통발생학)

계통발생학 (Phylogenetics)

개요

계통발생학 분석은 생물학적 서열(유전자, 단백질, 유전체)의 분기 패턴을 추론하여 진화 역사를 재구성합니다. 이 스킬은 표준 파이프라인을 다룹니다:

1. MAFFT — 다중 서열 정렬
2. IQ-TREE 2 — 모델 선택을 포함한 최대우도 계통수 추론
3. FastTree — 빠른 근사 최대우도 (대규모 데이터셋용)
4. ETE3 — 계통수 조작 및 시각화를 위한 Python 라이브러리

설치:

# Conda (CLI 도구에 권장)
conda install -c bioconda mafft iqtree fasttree
uv pip install ete3

# ete3의 TreeStyle/NodeStyle 렌더링은 Qt 백엔드가 필요하므로, 이미지 출력을 위해
# PyQt5가 필요합니다. 트리 파싱 및 통계는 없이도 작동합니다.
uv pip install PyQt5

이 스킬을 사용해야 하는 경우

다음과 같은 경우 계통발생학을 사용하세요:

• 진화적 관계: 내 서열과 가장 가까운 유연관계를 가진 유기체/유전자는 무엇인가?
• 바이러스 계통역학: 발병 확산을 추적하고 전파 날짜 추정
• 단백질 패밀리 분석: 유전자 패밀리 내 진화 관계 추론
• 수평적 유전자 전이 탐지: 일치하지 않는 종/유전자 트리를 가진 유전자 식별
• 조상 서열 재구성: 조상 단백질 서열 추론
• 분자시계 분석: 시간 샘플링을 이용한 분기 시기 추정
• GWAS 동반: 변이를 진화적 맥락에 배치 (예: SARS-CoV-2 변이)
• 미생물학: 16S rRNA 또는 코어 유전체 계통수를 이용한 종 계통수

표준 워크플로우

1. MAFFT를 이용한 다중 서열 정렬

import subprocess
import os

def run_mafft(input_fasta: str, output_fasta: str, method: str = "auto",
               n_threads: int = 4) -> str:
    """
    MAFFT로 서열을 정렬합니다.

    Args:
        input_fasta: 정렬되지 않은 FASTA 파일 경로
        output_fasta: 정렬된 출력 파일 경로
        method: 'auto' (자동 선택), 'einsi' (정확), 'linsi' (정확, 느림),
                'fftnsi' (중간), 'fftns' (빠름), 'retree2' (빠름)
        n_threads: CPU 스레드 수

    Returns:
        정렬된 FASTA 파일 경로
    """
    methods = {
        "auto": ["mafft", "--auto"],
        "einsi": ["mafft", "--genafpair", "--maxiterate", "1000"],
        "linsi": ["mafft", "--localpair", "--maxiterate", "1000"],
        "fftnsi": ["mafft", "--fftnsi"],
        "fftns": ["mafft", "--fftns"],
        "retree2": ["mafft", "--retree", "2"],
    }

    cmd = methods.get(method, methods["auto"])
    cmd += ["--thread", str(n_threads), "--inputorder", input_fasta]

    with open(output_fasta, 'w') as out:
        result = subprocess.run(cmd, stdout=out, stderr=subprocess.PIPE, text=True)

    if result.returncode != 0:
        raise RuntimeError(f"MAFFT 실패:\n{result.stderr}")

    # 정렬된 서열 수 계산
    with open(output_fasta) as f:
        n_seqs = sum(1 for line in f if line.startswith('>'))
    print(f"MAFFT: {n_seqs}개 서열 정렬 완료 → {output_fasta}")

    return output_fasta

# MAFFT 메서드 선택 가이드:
# 적은 서열(<200), 정확: linsi 또는 einsi
# 많은 서열(<1000), 중간: fftnsi
# 대규모 데이터셋(>1000): fftns 또는 auto
# 초고속(>10000): mafft --retree 1

2. 정렬 다듬기 (선택 사항이지만 권장)

def trim_alignment_trimal(aligned_fasta: str, output_fasta: str,
                            method: str = "automated1") -> str:
    """
    TrimAl로 잘 정렬되지 않은 열을 다듬습니다.

    Methods:
    - 'automated1': 자동 휴리스틱 (권장)
    - 'gappyout': 갭이 많은 열 제거
    - 'strict': 엄격한 갭 임계값
    """
    cmd = ["trimal", f"-{method}", "-in", aligned_fasta, "-out", output_fasta, "-fasta"]
    result = subprocess.run(cmd, capture_output=True, text=True)
    if result.returncode != 0:
        print(f"TrimAl 경고: {result.stderr}")
        # 다듬지 않은 정렬로 대체
        import shutil
        shutil.copy(aligned_fasta, output_fasta)
    return output_fasta

3. IQ-TREE 2 — 최대우도 계통수

def run_iqtree(aligned_fasta: str, output_prefix: str,
                model: str = "TEST", bootstrap: int = 1000,
                n_threads: int = 4, extra_args: list = None) -> dict:
    """
    IQ-TREE 2로 최대우도 계통수를 구축합니다.

    Args:
        aligned_fasta: 정렬된 FASTA 파일
        output_prefix: 출력 파일 접두사
        model: 'TEST'는 자동 모델 선택, 또는 직접 지정 (예: DNA는 'GTR+G',
               단백질은 'LG+G4', 'JTT+G')
        bootstrap: 초고속 부트스트랩 반복 횟수 (1000 권장)
        n_threads: 스레드 수 ('AUTO'로 자동 감지)
        extra_args: 추가 IQ-TREE 인수

    Returns:
        출력 파일 경로를 포함한 딕셔너리
    """
    cmd = [
        "iqtree2",
        "-s", aligned_fasta,
        "--prefix", output_prefix,
        "-m", model,
        "-B", str(bootstrap),   # 초고속 부트스트랩
        "-T", str(n_threads),
        "--redo"                # 기존 결과 덮어쓰기
    ]

    if extra_args:
        cmd.extend(extra_args)

    result = subprocess.run(cmd, capture_output=True, text=True)

    if result.returncode != 0:
        raise RuntimeError(f"IQ-TREE 실패:\n{result.stderr}")

    # 모델 선택 결과 출력
    log_file = f"{output_prefix}.log"
    if os.path.exists(log_file):
        with open(log_file) as f:
            for line in f:
                if "Best-fit model" in line:
                    print(f"IQ-TREE: {line.strip()}")

    output_files = {
        "tree": f"{output_prefix}.treefile",
        "log": f"{output_prefix}.log",
        "iqtree": f"{output_prefix}.iqtree",  # 전체 보고서
        "model": f"{output_prefix}.model.gz",
    }

    print(f"IQ-TREE: 계통수가 {output_files['tree']}에 저장됨")
    return output_files

# IQ-TREE 모델 선택 가이드:
# DNA:     TEST → GTR+G, HKY+G, TrN+G
# 단백질: TEST → LG+G4, WAG+G, JTT+G, Q.pfam+G
# 코돈:   TEST → MG+F3X4

# 시간(분자시계) 분석을 위해 다음 추가:
# extra_args = ["--date", "dates.txt", "--clock-test", "--date-CI", "95"]

4. FastTree — 빠른 근사 ML

IQ-TREE가 너무 느린 대규모 데이터셋(>1000 서열)의 경우:

def run_fasttree(aligned_fasta: str, output_tree: str,
                  sequence_type: str = "nt", model: str = "gtr",
                  n_threads: int = 4) -> str:
    """
    FastTree로 빠른 근사 ML 계통수를 구축합니다.

    Args:
        sequence_type: 'nt'는 뉴클레오타이드, 'aa'는 아미노산
        model: nt의 경우 'gtr' (권장) 또는 'jc'; aa의 경우 'lg', 'wag', 'jtt'
    """
    if sequence_type == "nt":
        cmd = ["FastTree", "-nt", "-gtr"]
    else:
        cmd = ["FastTree", f"-{model}"]

    cmd += [aligned_fasta]

    with open(output_tree, 'w') as out:
        result = subprocess.run(cmd, stdout=out, stderr=subprocess.PIPE, text=True)

    if result.returncode != 0:
        raise RuntimeError(f"FastTree 실패:\n{result.stderr}")

    print(f"FastTree: 계통수가 {output_tree}에 저장됨")
    return output_tree

5. ETE3를 이용한 계통수 분석 및 시각화

from ete3 import Tree, TreeStyle, NodeStyle, TextFace, PhyloTree
import matplotlib.pyplot as plt

def load_tree(tree_file: str) -> Tree:
    """Newick 계통수 파일을 로드합니다."""
    t = Tree(tree_file)
    print(f"계통수: 잎 {len(t)}개, 노드 {len(list(t.traverse()))}개")
    return t

def basic_tree_stats(t: Tree) -> dict:
    """기본 계통수 통계를 계산합니다."""
    leaves = t.get_leaves()
    distances = [t.get_distance(l1, l2) for l1 in leaves[:min(50, len(leaves))]
                 for l2 in leaves[:min(50, len(leaves))] if l1 != l2]

    stats = {
        "n_leaves": len(leaves),
        "n_internal_nodes": len(t) - len(leaves),
        "total_branch_length": sum(n.dist for n in t.traverse()),
        "max_leaf_distance": max(distances) if distances else 0,
        "mean_leaf_distance": sum(distances)/len(distances) if distances else 0,
    }
    return stats

def find_mrca(t: Tree, leaf_names: list) -> Tree:
    """잎 집합의 가장 최근 공통 조상을 찾습니다."""
    return t.get_common_ancestor(*leaf_names)

def visualize_tree(t: Tree, output_file: str = "tree.png",
                    show_branch_support: bool = True,
                    color_groups: dict = None,
                    width: int = 800) -> None:
    """
    계통수를 이미지로 렌더링합니다.

    Args:
        t: ETE3 Tree 객체
        color_groups: 잎 이름 → 색상 매핑 딕셔너리 (분류군 색상 지정용)
        show_branch_support: 부트스트랩 값 표시
    """
    ts = TreeStyle()
    ts.show_leaf_name = True
    ts.show_branch_support = show_branch_support
    ts.mode = "r"  # 'r' = 직사각형, 'c' = 원형

    if color_groups:
        for node in t.traverse():
            if node.is_leaf() and node.name in color_groups:
                nstyle = NodeStyle()
                nstyle["fgcolor"] = color_groups[node.name]
                nstyle["size"] = 8
                node.set_style(nstyle)

    t.render(output_file, tree_style=ts, w=width, units="px")
    print(f"계통수가 다음에 저장됨: {output_file}")

def midpoint_root(t: Tree) -> Tree:
    """중간점에 계통수를 뿌리내립니다 (외부군을 모를 때 사용)."""
    t.set_outgroup(t.get_midpoint_outgroup())
    return t

def prune_tree(t: Tree, keep_leaves: list) -> Tree:
    """지정된 잎만 남기도록 계통수를 가지치기합니다."""
    t.prune(keep_leaves, preserve_branch_length=True)
    return t

6. 전체 분석 스크립트

import subprocess, os
from ete3 import Tree

def full_phylogenetic_analysis(
    input_fasta: str,
    output_dir: str = "phylo_results",
    sequence_type: str = "nt",
    n_threads: int = 4,
    bootstrap: int = 1000,
    use_fasttree: bool = False
) -> dict:
    """
    완전한 계통발생학 파이프라인: 정렬 → 다듬기 → 계통수 → 시각화.

    Args:
        input_fasta: 정렬되지 않은 FASTA
        sequence_type: 'nt' (뉴클레오타이드) 또는 'aa' (아미노산/단백질)
        use_fasttree: IQ-TREE 대신 FastTree 사용 (대규모 데이터셋에서 더 빠름)
    """
    os.makedirs(output_dir, exist_ok=True)
    prefix = os.path.join(output_dir, "phylo")

    print("=" * 50)
    print("1단계: 다중 서열 정렬 (MAFFT)")
    aligned = run_mafft(input_fasta, f"{prefix}_aligned.fasta",
                         method="auto", n_threads=n_threads)

    print("\n2단계: 계통수 추론")
    if use_fasttree:
        tree_file = run_fasttree(
            aligned, f"{prefix}.tree",
            sequence_type=sequence_type,
            model="gtr" if sequence_type == "nt" else "lg"
        )
    else:
        model = "TEST" if sequence_type == "nt" else "TEST"
        iqtree_files = run_iqtree(
            aligned, prefix,
            model=model,
            bootstrap=bootstrap,
            n_threads=n_threads
        )
        tree_file = iqtree_files["tree"]

    print("\n3단계: 계통수 분석")
    t = Tree(tree_file)
    t = midpoint_root(t)

    stats = basic_tree_stats(t)
    print(f"계통수 통계: {stats}")

    print("\n4단계: 시각화")
    visualize_tree(t, f"{prefix}_tree.png", show_branch_support=True)

    # 뿌리내린 계통수 저장
    rooted_tree_file = f"{prefix}_rooted.nwk"
    t.write(format=1, outfile=rooted_tree_file)

    results = {
        "aligned_fasta": aligned,
        "tree_file": tree_file,
        "rooted_tree": rooted_tree_file,
        "visualization": f"{prefix}_tree.png",
        "stats": stats
    }

    print("\n" + "=" * 50)
    print("계통발생학 분석 완료!")
    print(f"결과 위치: {output_dir}/")
    return results

IQ-TREE 모델 가이드

DNA 모델

모델	설명	사용 사례
GTR+G4	일반 시간 가역 모델 + 감마	가장 유연한 DNA 모델
HKY+G4	Hasegawa-Kishino-Yano + 감마	두 비율 모델 (일반적)
TrN+G4	Tamura-Nei	불균등한 전이
JC	Jukes-Cantor	가장 간단함; 모든 비율 동일

단백질 모델

모델	설명	사용 사례
LG+G4	Le-Gascuel + 감마	최고 평균 단백질 모델
WAG+G4	Whelan-Goldman	널리 사용됨
JTT+G4	Jones-Taylor-Thornton	고전 모델
Q.pfam+G4	pfam 훈련	Pfam 유사 단백질 패밀리용
Q.bird+G4	조류 특이	척추동물 단백질

Tip: -m TEST를 사용하여 IQ-TREE가 최적의 모델을 자동으로 선택하도록 하세요.

모범 사례

• 정렬 품질 우선: 잘못된 정렬 → 신뢰할 수 없는 계통수; 정렬을 수동으로 확인
• 작은 정렬(<200 seq)에는 linsi, 큰 정렬에는 fftns 또는 auto 사용
• 모델 선택: 특별한 이유가 없는 한 항상 IQ-TREE에서 -m TEST 사용
• 부트스트랩: 가지 지지도에 ≥1000 초고속 부트스트랩(-B 1000) 사용
• 계통수 뿌리내리기: 뿌리 없는 계통수는 오해의 소지가 있음; 외부군 또는 중간점 뿌리내리기 사용
• >5000 서열은 FastTree: IQ-TREE가 느려짐; FastTree는 10–100배 빠름
• 긴 정렬 다듬기: TrimAl로 신뢰도 낮은 열 제거; 계통수 정확도 향상
• 재조합 확인: 계통수 구축 전 바이러스/박테리아 서열에서 재조합 확인 (RDP4, GARD)

추가 자료

• MAFFT: https://mafft.cbrc.jp/alignment/software/
• IQ-TREE 2: http://www.iqtree.org/ | 튜토리얼: https://www.iqtree.org/workshop/molevol2022
• FastTree: http://www.microbesonline.org/fasttree/
• ETE3: http://etetoolkit.org/
• FigTree (GUI 시각화): https://tree.bio.ed.ac.uk/software/figtree/
• iTOL (웹 시각화): https://itol.embl.de/
• MUSCLE (대체 정렬기): https://www.drive5.com/muscle/
• TrimAl (정렬 다듬기): https://vicfero.github.io/trimal/