系統樹の作成 — 再帰と二分木を使って進化の道筋を可視化する
このトピックを終えると
教科書で学んだ二分木と再帰を組み合わせて、Newick形式の系統樹を解析し、巡回してクレードごとの分析を行うツールを自分で作成できるようになります。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配列を解析して系統樹を作成しました。結果ファイルに次のような文字列があります。
(((Ecoli:0.05,Salmonella:0.06):0.02,Klebsiella:0.08):0.03,(Bacillus:0.15,Staph:0.14):0.10);これはNewick形式です。系統樹をテキストで表現する標準的な形式です。規則：
- 丸括弧
()は、1つのクレード（共通の祖先から分岐したグループ） - カンマ
,は、兄弟ノード - コロンの後の数字
:0.05は、親ノードまでの進化距離（枝の長さ） - セミコロン
;は、文字列の終わり
あなたがやりたいこと：
- 解析: 文字列をPythonのデータ構造に変換
- 探索: 各クレードのリーフノード（葉ノード）の種を列挙
- 距離計算: 2つの種間の総進化距離を計算
- 可視化: ツリーを図で表示
本質的なアプローチは再帰です。Newick文字列は、自身を子として含む再帰的な構造であるため、解析と探索はどちらも再帰によって自然に表現できます。
ブラックボックスからコンポーネントへ
コンポーネント1：ツリーノードの定義
from dataclasses import dataclass, fieldfrom typing import Optional
@dataclassclass 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パーサー（再帰）
パーサーを手動で実装すると、次の再帰構造が得られます。
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の再帰的な文法が、再帰関数に直接マッピングされます。
使用例：
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：再帰的なトラバース
これで、さまざまな方法でツリーをトラバースします。すべて再帰です。
すべての葉をリスト化する:
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ツリーの高さの計算:
def tree_height(node: TreeNode) -> int: if node.is_leaf: return 0 return 1 + max(tree_height(child) for child in node.children)インデントで表示:
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)出力：
(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：2つの葉の間の距離
2つの種間の進化距離は、共通の祖先までの距離の合計です。これを計算するには、まず各葉の祖先へのパスを見つける必要があります。
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フェーディング — 埋めるべき2つの空白
空白1：クラード統計
特定のクラード（内部ノード）の下にある葉の数と平均分岐長を計算します。
def clade_stats(node: TreeNode) -> dict: """ 返り値: { "leaf_count": 下の葉の数, "total_branch_length": 下のすべての分岐の長さの合計, "mean_branch_length": 平均 } """ if node.is_leaf: # TODO: 葉ノードの場合のデフォルト値 pass # TODO: 再帰呼び出しを使用して、各子ノードの統計を取得し、統合する passヒント: 葉の場合、{"leaf_count": 1, "total_branch_length": node.branch_length, ...}。内部ノードは、子ノードの結果を合計します。
空白2：葉の名前フィルターによるサブツリーの抽出
関心のある葉だけを残して、縮小されたツリーを作成します。
def prune_tree(node: TreeNode, keep_leaves: set[str]) -> Optional[TreeNode]: """ keep_leavesに含まれる葉だけを含む縮小されたツリーを返します。 不要な内部ノードを削除しますが、分岐長は結合します。 """ if node.is_leaf: # TODO: 葉がkeep_leavesに含まれている場合は返します。そうでない場合はNoneを返します。 pass # TODO: 子ノードを再帰的にプルーニングし、Noneでないものを保持します。 # 子ノードが0個の場合はNoneを返します。1個の場合は、その子ノードを返します（分岐長を結合します）。 passヒント: 子ノードが1つだけ残っている場合、この内部ノードは不要です。子ノードのbranch_lengthに、このノードのbranch_lengthを加えて返します。
考察 — 実際の系統樹ツールの違い
系統樹推定アルゴリズム: 皆さんが扱ったのは、すでに作成されたツリーを解析し、分析することです。実際の系統樹を配列から最初から作成することは、別の問題です。UPGMA（最も単純）、Neighbor-Joining（中程度の複雑さ）、Maximum Likelihood（RAxML、IQ-TREE）、Bayesian（MrBayes、BEAST）などがあります。
非二分木: Newick形式は、二分木でない場合があります。つまり、あるノードが3つ以上の子を持つことがあります（多分岐）。皆さんのパーサーは、すでにこの場合を処理しています。
ブートストラップ値: 実際のツリーは、各内部ノードにブートストラップ支持値（0〜100）を表示します。Newick形式では、通常、内部ノードの名前の場所に記述します。皆さんの _read_name 関数がこれを処理します。
視覚化: 実際のツールでは、ETE Toolkit、Biopython Phylo、iTOLを使用して描画します。matplotlibを使用して、単純な放射状または直交のツリーを描画することも可能です。
注釈の拡張: 実際のツールは、Nexus形式（系統樹+配列+メタデータ）またはphyloXML（XMLベース）を使用します。皆さんのNewickパーサーは、最小限の形式を扱います。
拡張プロジェクト
1. 可視化: matplotlib を使用して、作成した木を直交樹形図として描画する。葉の名前と枝の長さを表示する。
2. ブートストラップフィルタリング: ブートストラップのサポート値が低い内部ノードをプルーニングによって削除する。
3. UPGMAの実装: 距離行列からUPGMAを用いて木を最初から作成する。作成した TreeNode データ構造を再利用する。
4. Biopythonとの統合: BiopythonのPhyloモジュールを使用して、ファイルの入出力を行い、作成した解析関数を適用する。
この機能の部品図
- [F] 二分探索木（拡張）: TreeNodeによる自己参照構造。子を複数持つことができる一般的な木。
- [F] 再帰: パース、走査、距離計算はすべて再帰で実装。再帰と反復の比較。
- [W] ファイルI/O: Newickファイルの読み込みなど（完成したスクリプトとして提供）。
[F] = ユーザーが自分で実装 / [W] = 完成したコードとして提供。