alterlab-phylogenetics

系統発生学

概要

系統発生解析では、分岐パターンを推定することで、生物学的配列（遺伝子、タンパク質、ゲノム）の進化の歴史を再構築します。このスキルでは、標準的なパイプラインを扱います。

1. MAFFT — 多重配列アラインメント
2. IQ-TREE 3 — モデル選択を伴う最尤系統樹推定
3. FastTree — 高速な近似最尤法（大規模データセット向け）
4. ETE3 — 系統樹操作と可視化のための Python ライブラリ

インストール： アラインメントツールと系統樹構築ツールはコンパイル済みの CLI ツール（PyPI では提供されていません）です。bioconda または Homebrew でバイナリをインストールし、Python の可視化レイヤーは uv でインストールしてください。

# CLI binaries — bioconda (cross-platform) ...
conda install -c bioconda mafft iqtree fasttree trimal   # iqtree is now IQ-TREE 3
# ... or Homebrew on macOS (Apple Silicon): IQ-TREE/TrimAl live in the brewsci/bio tap
brew install mafft fasttree
brew tap brewsci/bio && brew install brewsci/bio/iqtree brewsci/bio/trimal

# Python visualization layer
uv pip install "ete3==3.1.3"   # also needs numpy<2 and PyQt5 for rendering

バイナリ名： bioconda の iqtree パッケージは現在 IQ-TREE 3（3.1.x）であり、実行ファイル iqtree と iqtree3 がインストールされます。このパッケージに iqtree2 はありません。2.x と 3.x の両方のパッケージで動作する iqtree を呼び出してください。PATH 上で複数のバージョンが共存する場合に限り、バージョン付きの名前を使用してください。IQ-TREE 3 では、MixtureFinder、複数系統樹の混合モデル（MAST）、遺伝子／部位の一致係数、大規模な病原体データセット向けの CMAPLE が追加されています。以下のオプションは 2.x から変更されていません。

ETE3（3.1.3）は ete3 の最終リリースであり、古い numpy／PyQt5 にバージョンを固定しているため、Python ≥3.12 ではインストールが不安定になる場合があります。t.render() が失敗する場合は、Newick 系統樹を書き出して FigTree／iTOL で表示するか、メンテナンスされている後継版の ete4 を使用してください（注：ete4 では TreeStyle／render API が変更されているため、以下のコード例は ete3 専用です）。

このスキルを使用する場合

次のような場合に系統発生学を使用します。

• 進化的関係： 自分の配列に最も近縁な生物／遺伝子はどれですか？
• ウイルスの系統動態： アウトブレイクの拡大を追跡し、感染時期を推定する
• タンパク質ファミリー解析： 遺伝子ファミリー内の進化的関係を推定する
• 水平遺伝子伝播の検出： 種系統樹／遺伝子系統樹と一致しない遺伝子を特定する
• 祖先配列の再構築： 祖先のタンパク質配列を推定する
• 分子時計解析： 時系列サンプリングを用いて分岐年代を推定する
• GWAS の補助： 変異を進化の文脈に位置付ける（例：SARS-CoV-2 の変異株）
• 微生物学： 16S rRNA またはコアゲノム系統発生に基づく種の系統発生

該当しないケース

ケース	代わりに使用
既存の Newick 系統樹の操作または比較（枝刈り、ルート設定、Robinson-Foulds、重複／種分化）	alterlab-etetoolkit
通常の配列解析、翻訳、またはフォーマット変換	alterlab-biopython
アラインメント前に配列を検索して見つけるための相同性検索	alterlab-blast
コミュニティの特徴量テーブルに対する多様性／序列化（UniFrac、PCoA、PERMANOVA）	alterlab-scikit-bio
ASV と系統樹を生成する 16S アンプリコンパイプラインの実行	alterlab-qiime2-amplicon

標準的なワークフロー

1. MAFFT による多重配列アラインメント

import subprocess
import os

def run_mafft(input_fasta: str, output_fasta: str, method: str = "auto",
               n_threads: int = 4) -> str:
    """
    Align sequences with MAFFT.

    Args:
        input_fasta: Path to unaligned FASTA file
        output_fasta: Path for aligned output
        method: 'auto' (auto-select), 'einsi' (accurate), 'linsi' (accurate, slow),
                'fftnsi' (medium), 'fftns' (fast), 'retree2' (fast)
        n_threads: Number of CPU threads

    Returns:
        Path to aligned FASTA file
    """
    methods = {
        "auto": ["mafft", "--auto"],
        "einsi": ["mafft", "--genafpair", "--maxiterate", "1000"],
        "linsi": ["mafft", "--localpair", "--maxiterate", "1000"],
        "fftnsi": ["mafft", "--retree", "2", "--maxiterate", "2"],
        "fftns": ["mafft", "--retree", "2", "--maxiterate", "0"],
        "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 failed:\n{result.stderr}")

    # Count aligned sequences
    with open(output_fasta) as f:
        n_seqs = sum(1 for line in f if line.startswith('>'))
    print(f"MAFFT: aligned {n_seqs} sequences → {output_fasta}")

    return output_fasta

# MAFFT method selection guide:
# Few sequences (<200), accurate: linsi or einsi
# Many sequences (<1000), moderate: fftnsi
# Large datasets (>1000): fftns or auto
# Ultra-fast (>10000): mafft --retree 1

2. アラインメントのトリミング（任意ですが推奨）

def trim_alignment_trimal(aligned_fasta: str, output_fasta: str,
                            method: str = "automated1") -> str:
    """
    Trim poorly aligned columns with TrimAl.

    Methods:
    - 'automated1': Automatic heuristic (recommended)
    - 'gappyout': Remove gappy columns
    - 'strict': Strict gap threshold
    """
    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 warning: {result.stderr}")
        # Fall back to using the untrimmed alignment
        import shutil
        shutil.copy(aligned_fasta, output_fasta)
    return output_fasta

3. IQ-TREE 3 — 最尤系統樹

def run_iqtree(aligned_fasta: str, output_prefix: str,
                model: str = "TEST", bootstrap: int = 1000,
                n_threads: int = 4, extra_args: list = None) -> dict:
    """
    Build a maximum likelihood tree with IQ-TREE 3.

    Args:
        aligned_fasta: Aligned FASTA file
        output_prefix: Prefix for output files
        model: 'TEST' for automatic model selection, or specify (e.g., 'GTR+G' for DNA,
               'LG+G4' for proteins, 'JTT+G' for proteins)
        bootstrap: Number of ultrafast bootstrap replicates (1000 recommended)
        n_threads: Number of threads ('AUTO' to auto-detect)
        extra_args: Additional IQ-TREE arguments

    Returns:
        Dict with paths to output files
    """
    cmd = [
        "iqtree",
        "-s", aligned_fasta,
        "--prefix", output_prefix,
        "-m", model,
        "-B", str(bootstrap),   # Ultrafast bootstrap
        "-T", str(n_threads),
        "--redo"                # Overwrite existing results
    ]

    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 failed:\n{result.stderr}")

    # Print model selection result
    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",  # Full report
        "model": f"{output_prefix}.model.gz",
    }

    print(f"IQ-TREE: Tree saved to {output_files['tree']}")
    return output_files

# Model selection: ModelFinder Plus (-m MFP) is the DEFAULT — omitting -m runs it.
# -m TEST is the narrower jModelTest/ProtTest-style search (no FreeRate models);
# use it only when you need to match those tools.
# Typical winners — DNA: GTR+F+I+G4, HKY+F+G4; Protein: LG+F+G4, Q.pfam+G4;
# Codon: MG+F3X4.

# For temporal (molecular clock) analysis, add:
# 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:
    """
    Build a fast approximate ML tree with FastTree.

    Args:
        sequence_type: 'nt' for nucleotide or 'aa' for amino acid
        model: For nt: 'gtr' (recommended) or 'jc'; for 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 failed:\n{result.stderr}")

    print(f"FastTree: Tree saved to {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:
    """Load a Newick tree file."""
    t = Tree(tree_file)
    print(f"Tree: {len(t)} leaves, {len(list(t.traverse()))} nodes")
    return t

def basic_tree_stats(t: Tree) -> dict:
    """Compute basic tree statistics."""
    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:
    """Find the most recent common ancestor of a set of leaves."""
    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:
    """
    Render phylogenetic tree to image.

    Args:
        t: ETE3 Tree object
        color_groups: Dict mapping leaf_name → color (for coloring taxa)
        show_branch_support: Show bootstrap values
    """
    ts = TreeStyle()
    ts.show_leaf_name = True
    ts.show_branch_support = show_branch_support
    ts.mode = "r"  # 'r' = rectangular, 'c' = circular

    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"Tree saved to: {output_file}")

def midpoint_root(t: Tree) -> Tree:
    """Root tree at midpoint (use when outgroup unknown)."""
    t.set_outgroup(t.get_midpoint_outgroup())
    return t

def prune_tree(t: Tree, keep_leaves: list) -> Tree:
    """Prune tree to keep only specified leaves."""
    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:
    """
    Complete phylogenetic pipeline: align → trim → tree → visualize.

    Args:
        input_fasta: Unaligned FASTA
        sequence_type: 'nt' (nucleotide) or 'aa' (amino acid/protein)
        use_fasttree: Use FastTree instead of IQ-TREE (faster for large datasets)
    """
    os.makedirs(output_dir, exist_ok=True)
    prefix = os.path.join(output_dir, "phylo")

    print("=" * 50)
    print("Step 1: Multiple Sequence Alignment (MAFFT)")
    aligned = run_mafft(input_fasta, f"{prefix}_aligned.fasta",
                         method="auto", n_threads=n_threads)

    print("\nStep 2: Tree Inference")
    if use_fasttree:
        tree_file = run_fasttree(
            aligned, f"{prefix}.tree",
            sequence_type=sequence_type,
            model="gtr" if sequence_type == "nt" else "lg"
        )
    else:
        # -m TEST auto-detects the alphabet (nt vs aa) and selects the best model.
        iqtree_files = run_iqtree(
            aligned, prefix,
            model="TEST",
            bootstrap=bootstrap,
            n_threads=n_threads
        )
        tree_file = iqtree_files["tree"]

    print("\nStep 3: Tree Analysis")
    t = Tree(tree_file)
    t = midpoint_root(t)

    stats = basic_tree_stats(t)
    print(f"Tree statistics: {stats}")

    print("\nStep 4: Visualization")
    visualize_tree(t, f"{prefix}_tree.png", show_branch_support=True)

    # Save rooted tree
    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("Phylogenetic analysis complete!")
    print(f"Results in: {output_dir}/")
    return results

IQ-TREE モデルガイド

DNA モデル

モデル	説明	使用例
GTR+G4	一般時間可逆モデル + ガンマ分布	最も柔軟な DNA モデル
HKY+G4	Hasegawa-Kishino-Yano + ガンマ分布	2 速度モデル（一般的）
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で学習（QMaker）	一般的なタンパク質ファミリー
Q.bird+G4	鳥類クレード特異的（QMaker）	鳥類のタンパク質。関連モデル：Q.mammal、Q.insect、Q.yeast、Q.plant

ヒント： IQ-TREEはデフォルトでモデルを自動選択します（-m MFP）。モデルを固定する必要がある場合にのみ、明示的に指定してください。

ベストプラクティス

• まずアラインメントの品質を確認します：アラインメントが不適切だと、信頼性の低い系統樹になります。アラインメントを手動で確認してください
• 小規模（200配列未満）にはlinsi、大規模なアラインメントにはfftnsまたはautoを使用します
• モデル選択：ModelFinder PlusはIQ-TREEのデフォルトです。-mを省略するか、-m MFPを指定してください。 -m TESTはjModelTest／ProtTest形式のモデルに検索を限定するため、これらのツールに合わせる場合にのみ使用してください
• ブートストラップ：枝の支持値を得るには、ウルトラファスト・ブートストラップを1000回以上（-B 1000）実行してください
• 系統樹に根を付けます：根なし系統樹は誤解を招く可能性があります。外群または中点を使って根付けしてください
• 5000配列を超える場合はFastTreeを使用します：IQ-TREEは遅くなるため、FastTreeなら10〜100倍高速です
• 長いアラインメントをトリミングします：TrimAlは信頼性の低い列を除去し、系統樹の精度を向上させます
• 系統樹の構築前に、ウイルス／細菌配列の組換えを確認します（RDP4、GARD）

追加リソース

• MAFFT：https://mafft.cbrc.jp/alignment/software/
• IQ-TREE 3：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/