系統發生學

系統發生學 是做什麼的？

概述

系統發生分析旨在透過推論物種或序列間的分支模式，重建生物序列（基因、蛋白質、基因體）的演化歷史。本技能涵蓋標準的處理流程：

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 進行序列比對。

    參數：
        input_fasta: 未比對的 FASTA 檔案路徑
        output_fasta: 比對後的輸出檔案路徑
        method: 'auto'（自動選擇）、'einsi'（精確）、'linsi'（精確，較慢）、
                'fftnsi'（中等）、'fftns'（快速）、'retree2'（快速）
        n_threads: CPU 執行緒數量

    回傳：
        比對後 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 修剪比對不佳的欄位。

    方法：
    - '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 建立最大概似法樹。

    參數：
        aligned_fasta: 比對後的 FASTA 檔案
        output_prefix: 輸出檔案的前綴
        model: 'TEST' 表示自動模型選擇，或直接指定（例如 DNA 用 'GTR+G'，
               蛋白質用 'LG+G4' 或 'JTT+G'）
        bootstrap: 超快拔靴法重複次數（建議 1000）
        n_threads: 執行緒數量（使用 'AUTO' 自動偵測）
        extra_args: 其他 IQ-TREE 參數

    回傳：
        包含輸出檔案路徑的字典
    """
    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 — 快速近似最大概似法

對於大型資料集（>1000 條序列），若 IQ-TREE 執行速度過慢時使用：

def run_fasttree(aligned_fasta: str, output_tree: str,
                  sequence_type: str = "nt", model: str = "gtr",
                  n_threads: int = 4) -> str:
    """
    使用 FastTree 建立快速近似最大概似法樹。

    參數：
        sequence_type: 'nt' 代表核苷酸，'aa' 代表胺基酸
        model: 核苷酸用 'gtr'（建議）或 'jc'；胺基酸用 '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:
    """
    將系統發生樹渲染成圖片。

    參數：
        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:
    """
    完整的系統發生學管線：比對 → 修剪 → 建樹 → 視覺化。

    參數：
        input_fasta: 未比對的 FASTA 檔案
        sequence_type: 'nt'（核苷酸）或 'aa'（胺基酸/蛋白質）
        use_fasttree: 使用 FastTree 代替 IQ-TREE（大型資料集時更快）
    """
    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("\n步驟 2：樹狀圖推論")
    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("\n步驟 3：樹狀圖分析")
    t = Tree(tree_file)
    t = midpoint_root(t)

    stats = basic_tree_stats(t)
    print(f"樹狀圖統計資料：{stats}")

    print("\n步驟 4：視覺化")
    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	廣義時間可逆 + Gamma	最靈活的 DNA 模型
HKY+G4	Hasegawa-Kishino-Yano + Gamma	雙速率模型（常見）
TrN+G4	Tamura-Nei	不對等轉換
JC	Jukes-Cantor	最簡單；所有速率相同

蛋白質模型

模型	說明	使用情境
LG+G4	Le-Gascuel + Gamma	最佳平均蛋白質模型
WAG+G4	Whelan-Goldman	廣泛使用
JTT+G4	Jones-Taylor-Thornton	經典模型
Q.pfam+G4	以 pfam 訓練的模型	適用於 Pfam 類別的蛋白質家族
Q.bird+G4	鳥類專用模型	脊椎動物蛋白質

提示： 使用 -m TEST 讓 IQ-TREE 自動選取最佳模型。

最佳作法

• 比對品質優先：比對不良會導致不可靠的樹；請手動檢查比對結果
• 少量序列 (<200) 使用 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（圖形化介面視覺化）：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/