モチーフ検索：トライブを用いてプロモーター配列中の転写因子結合部位を特定する
このトピックを終えると
教科書で学んだトライ（Trie）と二分探索を組み合わせて、プロモーター配列から複数の転写因子の結合部位を一度に検索するツールを自分で作成できるようになります。Aho-Corasickのような実用的な多重パターンマッチングアルゴリズムの基礎となるデータ構造をコードで理解します。
この記事は教育用の一般的な例です。実際のモチーフ検索は、位置重み行列（PWM）、背景モデルなど、より高度な手法を用います。
「500個のモチーフ × 100kbの配列」— 単純な検索の落とし穴
マウスプロモーター100kbの配列から、既知の転写因子の結合部位である500個のモチーフ（それぞれ6〜12bp）をすべて見つけたいとします。
単純なアプローチ：
def naive_search(sequence: str, motifs: list[str]) -> list[tuple[int, str]]: hits = [] for motif in motifs: for i in range(len(sequence) - len(motif) + 1): if sequence[i:i+len(motif)] == motif: hits.append((i, motif)) return hitsこのアプローチの問題は、モチーフごとに繰り返し探索を行うことです。時間計算量はO（モチーフ数 × 配列の長さ）= O（500 × 100,000）= 5 × 10⁷です。Pythonで数秒で済みます。
もし配列が30億bpのヒトゲノムであればどうでしょうか？ 5 × 10⁵ × 3 × 10⁹ = 1.5 × 10¹⁵です。宇宙が終わるまでに終わらないでしょう。
現実的なアプローチはトライです。500個のモチーフを1つの木構造に統合し、配列を1回だけ走査しながら、木構造をたどってマッチを検出します。時間計算量はO（配列の長さ + マッチの数）に短縮されます。
ブラックボックスからパーツへ — トライを覗いてみよう
パーツ 1：トライデータ構造
トライは、各ノードが 1 つの文字を表し、ルートからリーフへのパスが 1 つの文字列を表すツリーです。
例として、モチーフのセット {"TATA", "TATT", "GCGC"} をトライで表現すると：
root
/ \
T G
| |
A C
| |
T G
/ \ |
A* T* C** は、完結したモチーフの終端を示します。3 つのモチーフが共通の接頭辞（TA）を共有する場合、その部分はツリーに 1 回だけ保存されます。
Python による実装：
class TrieNode: def __init__(self) -> None: self.children: dict[str, TrieNode] = {} self.matches: list[str] = []
class Trie: def __init__(self) -> None: self.root = TrieNode() def insert(self, pattern: str) -> None: node = self.root for char in pattern: if char not in node.children: node.children[char] = TrieNode() node = node.children[char] node.matches.append(pattern)パーツ 2：トライを巡回して検索
次に、配列の各位置でトライをたどりながらマッチングを確認します。
def search(trie: Trie, sequence: str) -> list[tuple[int, str]]: hits: list[tuple[int, str]] = [] for i in range(len(sequence)): node = trie.root j = i while j < len(sequence) and sequence[j] in node.children: node = node.children[sequence[j]] for match in node.matches: hits.append((i, match)) j += 1 return hits各位置 i から開始し、ツリーの最大深さ（最長のモチーフの長さ）までたどります。500 個のモチーフの最大長が 12 の場合、各位置で最大 12 回の検索が行われます。全体の計算量は O（配列の長さ × 最大モチーフの長さ）です。
実践的な最適化（Aho-Corasick）では、ここで フェイルリンク を追加して、再開位置を事前に計算します。このチュートリアルでは、基本的なトライを使って概念を確立します。
二分探索との比較
トライの仕組みを理解するために、代替的な方法である「ソートされた配列 + 二分探索」と比較すると、より明確になります。
ソートされたモチーフ配列によるアプローチ：
def sorted_search(sequence: str, sorted_motifs: list[str], max_motif_len: int) -> list[tuple[int, str]]: from bisect import bisect_left hits = [] for i in range(len(sequence)): for length in range(1, max_motif_len + 1): substring = sequence[i:i+length] idx = bisect_left(sorted_motifs, substring) if idx < len(sorted_motifs) and sorted_motifs[idx] == substring: hits.append((i, substring)) return hitsこのアプローチの時間計算量：O(シーケンスの長さ × 最大モチーフ長 × log(モチーフの数))。
トライと比較すると、logの項はありますが、実際にはトライの方が高速です。その理由はキャッシュの局所性です。トライの探索は、関連するノードが連続したメモリにアクセスされるのに対し、二分探索は配列のあちこちをランダムにアクセスします。
二分探索は、モチーフの集合が頻繁に変化せず、ロード時間が重要な場合に、トライよりも有利です。トライの構築は遅いですが、検索は高速です。二分探索はソートするだけでよいため、初期ロードは高速です。これはトレードオフです。
実践的な使用例
簡単なシナリオで確認してみましょう。
# 既知の転写因子結合モチーフ（例）motifs = [ "TATAAA", # TATA box "CAAT", # CAAT box "GGGCGG", # GC box "CACGTG", # E-box "TGACTCA" # AP-1]
# 仮想プロモーター配列promoter = "GCTATAAACCAATGGGCGGATGCACGTGCCCTGACTCAAG"
# トライを構築trie = Trie()for motif in motifs: trie.insert(motif)
# 検索hits = search(trie, promoter)for pos, motif in sorted(hits): print(f"位置 {pos}: {motif}")
# 出力：# 位置 2: TATAAA# 位置 9: CAAT# 位置 12: GGGCGG# 位置 21: CACGTG# 位置 28: TGACTCA一度の配列走査で、5つのモチーフがすべて検出されました。
フェーディング — 埋めるべき2つの空白
空白1：IUPACコードのサポート
実際の転写因子結合部位は可変性を持ちます。IUPACコード（R = A/G、Y = C/T、W = A/Tなど）で表現されます。
IUPAC = { "A": {"A"}, "C": {"C"}, "G": {"G"}, "T": {"T"}, "R": {"A", "G"}, "Y": {"C", "T"}, "S": {"C", "G"}, "W": {"A", "T"}, "K": {"G", "T"}, "M": {"A", "C"}, "B": {"C", "G", "T"}, "D": {"A", "G", "T"}, "H": {"A", "C", "T"}, "V": {"A", "C", "G"}, "N": {"A", "C", "G", "T"}}
def search_with_iupac(trie: Trie, sequence: str) -> list[tuple[int, str]]: """ モチーフにIUPACコードが含まれている可能性がある。 例："TATA**W**A" はTATAAAまたはTATATAのどちらにもマッチする。 """ # TODO：トライを探索するとき、各ノードのchildrenキーがIUPACコードの場合、拡張セットでマッチさせる。 passヒント: ツリーの構築はそのままにして、検索時にsequence[j]と各子ノードのキーをIUPAC拡張セットでマッチさせる。
空白2：逆相補配列の自動検索
DNAは二重鎖であるため、モチーフは逆相補配列としても現れる可能性があります。
def reverse_complement(seq: str) -> str: complement = {"A": "T", "T": "A", "G": "C", "C": "G"} return "".join(complement.get(b, b) for b in reversed(seq))
def search_both_strands(motifs: list[str], sequence: str) -> list[tuple[int, str, str]]: """ 各モチーフとその逆相補配列をすべてトライに入れて検索する。 返り値：（位置、モチーフ、“forward”または“reverse”） """ # TODO：トライを構築するとき、各モチーフの正方向と逆相補配列の両方を挿入する。 # 各エントリで元のモチーフと方向を記録する。 passヒント: TrieNode.matches を list[tuple[str, str]] （元のモチーフ、方向）に拡張する。
考察：実用的なモチーフ検索との違い
Position Weight Matrix (PWM)：実際には、ほとんどの転写因子の結合部位は、厳密な文字列マッチングではなく確率的なマッチングである。各位置の各塩基に確率が割り当てられ、候補配列のスコアを計算する。JASPARのようなデータベースのモチーフは、この形式である。
背景モデル：実際には、ランダムなマッチング確率を背景モデルとして計算し、統計的に有意なマッチングのみを報告する。FIMO（MEMEスイート）が標準的なツールである。
Aho-Corasick：トライに失敗リンクを追加すると、O（配列長 + マッチ数）となり、真に線形時間になる。実用的な複数パターンマッチングの標準。
全ゲノム規模：ヒトゲノム（30億bp）の検索には、サフィックス配列、FMインデックスなどのインデックスデータ構造が代替手段となる。BWA、Bowtieなどのアライメントツールの中心部。
拡張プロジェクト
1. JASPARデータのロード: 実際のJASPARデータベースから、ヒトの転写因子100個分のPWM（位置特異的スコア行列）をダウンロードし、作成したツールでプロモーター領域のスキャンを行う。
2. 可視化: スキャン結果をプロモーターマップとして表示する。matplotlibやIGVスタイルのトラックを使用する。
3. Aho-Corasickの拡張: フェイルリンクを追加し、実際に線形時間でマッチングを行うツールを完成させる。
このセクションの構成要素
- [F] トライ（Trie）: 複数の文字列を木構造に統合。ノード＝文字、パス＝文字列。
- [F] 二分探索: ソートされた配列＋
bisectによる代替アプローチ。トライとのトレードオフを理解する。 - [W] ファイル入出力: FASTAファイルの解析など（完成したスクリプトとして提供）。
[F] = 自分で実装 / [W] = 完成したコードとして提供。