系統解析

系統解析 は何をしますか？

概要

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

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

    # アライメントされた配列数を数える
    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メソッド選択ガイド：
# 少数配列（<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 warning: {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 failed:\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: Tree saved to {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

大規模データセット（>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で高速近似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 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:
    """Newick形式の系統樹ファイルを読み込みます。"""
    t = Tree(tree_file)
    print(f"Tree: {len(t)} leaves, {len(list(t.traverse()))} nodes")
    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:
    """指定したリーフ群の最新共通祖先（MRCA）を見つけます。"""
    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"Tree saved to: {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("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:
        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("\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)

    # ルート付けされた系統樹を保存
    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	General Time Reversible + 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/ | Tutorial: 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（Web可視化）: https://itol.embl.de/
• MUSCLE（代替アライナー）: https://www.drive5.com/muscle/
• TrimAl（アライメントトリミング）: https://vicfero.github.io/trimal/