順序整列 — 動的計画法を用いて2つのDNAを正確に重ね合わせる
このトピックを終えると
教科書で学んだ動的計画法、再帰と反復の比較、2次元行列を組み合わせて、2つのDNA配列を正確に重ね合わせ、どこに挿入・欠失・置換があるかを特定するアライメントツールを自分で作成できるようになります。BLASTが魔法のように行うそのアライメントの内部原理をコードで理解できます。
この記事は教育用の一般的な例です。変異解析はバイオインフォマティクスで毎日行う作業なので、題材として選びました。
「なぜ突然変異の位置が一致しないのか？」— 素朴な比較の落とし穴
BRCA1遺伝子の野生型部分配列と、患者からシーケンスされた突然変異配列を比較するとしましょう。
野生型: ATGCTAGCATGCA
突然変異: ATGCAGCATGCA最も単純な比較は、位置ごとに一つずつ対応させることです。
def compare_naive(seq1: str, seq2: str) -> list[tuple[int, str, str]]: diffs = [] for i in range(min(len(seq1), len(seq2))): if seq1[i] != seq2[i]: diffs.append((i, seq1[i], seq2[i])) return diffsこれを上記の2つの配列に適用すると、次のようになります。
位置4: T → A
位置5: A → G
位置6: G → C
位置7: C → A
位置8: A → T
位置9: T → G
位置10: G → C
位置11: C → Aおかしいですね。位置4以降は、ほとんどすべての位置が異なっています。本当にこれほど多くの突然変異があるのでしょうか？
実際に2つの配列を目で見てみると、答えが見えてきます。野生型で位置4のTが削除されただけです。
野生型: ATGC[T]AGCATGCA
突然変異: ATGC AGCATGCA
↑ ここからオフセットが1つずつずれる1文字の削除が、以降の位置をすべてずらして、まるで多くの置換があるかのように見せています。素朴な比較では、この削除・挿入に気づくことができません。
この問題を解決するのが、**配列アライメント（sequence alignment）です。2つの配列の間に、必要なだけギャップ（gap、-文字）**を挿入して、ずれた部分を補正します。
野生型: ATGCTAGCATGCA
突然変異: ATGC-AGCATGCA ← 4番目の位置にギャップこれで、位置5以降は完全に一致します。実際の変異が1つの**欠失（deletion）**であることが明確になります。
問題は、このギャップをどこに挿入するかを決定することです。2つの配列が100文字程度の場合、ギャップを挿入できる位置の組み合わせが爆発的に増加します。どの組み合わせが「最も自然な」アライメントなのか、どのように見つけるのでしょうか？この記事では、その答えを最初から作り上げます。
まずは完成品を見てみましょう（ブラックボックスを先に実行する）
これから作るツールは、2つの配列を受け取り、最適なアラインメントを返します。
aligned1, aligned2, score = align(seq1, seq2)
print(aligned1) # 'ATGCTAGCATGCA'print(aligned2) # 'ATGC-AGCATGCA'print(score) # 8 (例)=== アラインメント結果 ===
Wild: ATGCTAGCATGCA
Mut: ATGC-AGCATGCA
スコア: 8
=== アラインメントなしで単純に比較 ===
Wild: ATGCTAGCATGCA
Mut: ATGCAGCATGCA
不一致の位置: 8箇所（ほぼすべて）
=== アラインメント後の実際の変異 ===
位置4: Tの欠失（deletion） — 実際の変異1箇所同じ2つの配列ですが、アラインメントの有無で解釈が大きく異なります。**アラインメントは、変異の解釈における最初の段階であり、決定的な段階です。**これから、このアラインメントを作成する原理を解説していきます。
このツールは何の部品で構成されているか（部品分解図）
系列整列器（アライナー）
┌──────────────────────────────────────────────────┐
│ [入力] 2つの系列を読み込む ─── 部品：文字列/ファイル入出力 │ ← 完成品として提供（ツール）
│ │ │
│ ▼ │
│ [1段階] スコアマトリックスを埋める │
│ 部品：2D配列 + 動的計画法 │ ← 自分で作成する ★
│ │ │
│ ▼ │
│ [2段階] 後方追跡で最適な経路を復元 │
│ 部品：再帰・反復選択 │ ← 自分で作成する ★
│ │ │
│ ▼ │
│ [出力] 整列された2つの系列 + スコア │
└──────────────────────────────────────────────────┘| 部品 | どこで学んだか | このツールで何をするか |
|---|---|---|
| 文字列入出力 | string-and-file-io | 系列を読み込んで処理する |
| 2D配列 | matrix-2d-array | 部分問題の答えを保存する格子 |
| 動的計画法 | dynamic-programming | 部分問題の答えを再利用して、計算量の爆発を回避する |
| 再帰 vs 反復 | recursion-vs-iteration | 格子を埋める時、結果を追跡する時に、どちらの方法を使うかを選択する |
📌 これらの概念を初めて見る場合は（上部のリンク）
自分で作成する新しい概念は、動的計画法・2Dマトリックス・再帰/反復選択の3つ。 系列の処理はツールなので、完成品として提供します。 たったの3つです。 認知的な限界を超えません。
ステップ1：整列準備（サンプルとして提供）
まず、整列する2つの配列を準備します。実際にはFASTAファイルから読み込みますが、ブラウザでの実習のために文字列で定義します。
# 野生型BRCA1の部分配列（教育用サンプル）wild = "ATGCTAGCATGCA"
# 変異配列（4番目の位置にTの欠失）mut = "ATGCAGCATGCA"
assert len(wild) == 13assert len(mut) == 12assert all(c in "ACGT" for c in wild)assert all(c in "ACGT" for c in mut)この2つの配列を、ギャップ-を適切に挿入して、重ね合わせるのが目標です。
ステップ2：スコアリングルールを定義する（完成版）
まず、アライメントの「品質」をどのように測定するかを定義する必要があります。一般的な方法は、スコアリングシステムを使用することです。
MATCH = 1 # 同じ塩基同士がアライメントされた場合、+1MISMATCH = -1 # 異なる塩基同士がアライメントされた場合、-1GAP = -2 # ギャップが1つ挿入された場合、-2
def score_pair(a: str, b: str) -> int: """2つの塩基（またはギャップ）のアライメントスコア""" if a == b: return MATCH return MISMATCH
assert score_pair("A", "A") == 1assert score_pair("A", "T") == -1このスコアリングルールは、アライメントアルゴリズムの好みを決定します。
- ギャップペナルティが低い場合（例：-0.5）、ギャップが多いアライメントを好む
- ギャップペナルティが高い場合（例：-5）、ギャップよりもミスマッチを許容する
バイオインフォマティクスの実務では、BLOSUM / PAMなどの洗練されたスコアリング行列が使用されます。ここでは、学習のために単純なルールを使用します。
3段階で作成 — まずは再帰から始めて、その爆発を体験する ★ (再帰 vs 繰り返し)
✍️ 自分で埋めるセクション。 部品 = 再帰。まず、再帰を使ってアラインメントスコアを計算し、なぜ実用的ではないのかを実感する。
アラインメントスコアの問題を再帰で表現すると、驚くほど単純になります。2つのシーケンスの末尾から始めて、各位置で3つの選択肢を再帰的に試します。
- 2つのシーケンスの最後の文字をアラインメントする →
score_pair(s1[-1], s2[-1])+align(s1[:-1], s2[:-1]) s1の最後の文字とギャップをアラインメントする →GAP+align(s1[:-1], s2)s2の最後の文字とギャップをアラインメントする →GAP+align(s1, s2[:-1])
3つのオプションのうち、最も高いスコアを選択します。これを再帰で記述すると、次のようになります。
def align_recursive(s1: str, s2: str) -> int: """最適なアラインメントのスコアのみを再帰的に計算（バックトレースなし）。""" if not s1: return len(s2) * GAP if not s2: return len(s1) * GAP return max( align_recursive(s1[:-1], s2[:-1]) + score_pair(s1[-1], s2[-1]), align_recursive(s1[:-1], s2) + GAP, align_recursive(s1, s2[:-1]) + GAP, )
# 検証：短い例でのみassert align_recursive("AA", "AA") == 2 # 2つのマッチassert align_recursive("AT", "A") == 1 + GAP # マッチ + ギャップ = 1 - 2 = -1これは正しい答えを返します。しかし、長いシーケンスには絶対に使用できません。 次のコードを実行すると、その理由がわかります。
import time
# 長さ15程度でも、再帰の地獄short1 = "ATGCTAGCATGCAAG"short2 = "ATGCAGCATGCAAG"
t0 = time.time()score = align_recursive(short1, short2)elapsed = time.time() - t0print(f"長さ15 再帰: {elapsed:.2f}秒")# 数秒かかります。長さ20なら数分、25なら数時間その理由は、再帰が同じサブ問題を何百万回も繰り返し計算するためです。align_recursive("ATGCT", "ATGC") の値は、再帰ツリーのさまざまな場所で何度も計算されます。毎回最初から。
これがO(3^n) の時間複雑度の力です。この時点で、動的計画法が登場します。
🤔 自己説明プロンプト
align_recursive("AT", "A")を再帰ツリーとして描いてみてください。何個のサブ問題が繰り返し登場しますか？（長さ3、4程度に拡張するだけでも、重複が目立ちます。）
4단계：格子に答えを保存する（2D行列 + 動的計画法）
✍️ 手動で埋めるセクション。 コンポーネント = 2D配列 + 動的計画法。部分問題の答えを格子に保存することで、繰り返し計算を回避します。
アイデアはこうです。align(s1[:i], s2[:j])の答えを**格子 dp[i][j]**に保存します。これにより、同じ部分問題を2回計算する必要がなくなります。
🔎 動的計画法とは（図解 — dynamic-programming） 大きな問題の答えが小さな問題の答えから作られる場合、小さな問題の答えを保存しておくと、2回計算する必要がなくなります。再帰の指数関数的な増加を、保存（メモ化）によって多項式時間で解決します。格子に埋める方法は、特に「ボトムアップ動的計画法」と呼ばれます。
格子のサイズは(len(s1)+1) × (len(s2)+1)。余分な1行/1列は、空文字列とのアラインメントのための基本ケースです。
def build_score_matrix(s1: str, s2: str) -> list[list[int]]: """dp[i][j] = s1[:i]とs2[:j]の最適なアラインメントスコア。""" n, m = len(s1), len(s2) dp = [[0] * (m + 1) for _ in range(n + 1)]
# 基本ケース：空文字列とアラインメントすると、ギャップだけ for i in range(1, n + 1): dp[i][0] = dp[i-1][0] + GAP for j in range(1, m + 1): dp[0][j] = dp[0][j-1] + GAP
# 格子を埋める（ボトムアップ） for i in range(1, n + 1): for j in range(1, m + 1): match = dp[i-1][j-1] + score_pair(s1[i-1], s2[j-1]) delete = dp[i-1][j] + GAP # s1の文字とギャップ insert = dp[i][j-1] + GAP # ギャップとs2の文字 dp[i][j] = max(match, delete, insert) return dp
dp = build_score_matrix(wild, mut)
# 検証：最後のセルは再帰的な答えと一致する必要があるassert dp[len(wild)][len(mut)] == align_recursive(wild, mut)# サイズの検証assert len(dp) == len(wild) + 1assert len(dp[0]) == len(mut) + 1# 左側の列は純粋なギャップの累積assert dp[3][0] == 3 * GAP同じ答えを得るのに、速度は全く異なります。
import time
# 先ほど、再帰で数秒かかった長さ15t0 = time.time()score = build_score_matrix(short1, short2)[len(short1)][len(short2)]elapsed = time.time() - t0print(f"長さ15 動的計画法: {elapsed*1000:.2f}ms")# ミリ秒レベル再帰で数秒かかったものが、ミリ秒以内に完了します。これが**O(n×m)**の力です。格子のセルは一つ一つ、一度だけ計算されます。
🔎 なぜ格子が答えなのか（図解 — matrix-2d-array） 2つのインデックス(i, j)に依存する部分問題は、自然に2Dの格子で表現されます。各セルが1つの部分問題を、セル間の矢印が依存関係を表します。この構造がわかると、DP問題へのアプローチが容易になります。
🤔 自己説明プロンプト
dp[i][j]を計算する際に、dp[i-1][j-1]、dp[i-1][j]、dp[i][j-1]のみが必要です。格子全体を保存する必要があるでしょうか？ 実際には、直前の1行だけがあれば、次の行を埋めることができるため、スコアだけはO(min(n,m))の空間で計算できます。いつその最適化ができないでしょうか？（ヒント：逆追跡）
5단계：格子を辿ってアラインメントを復元する ★
✍️ 手動で埋めるセクション。 部品 = 反復的な逆探索。格子の最後のセルから開始し、その値がどのように生成されたかを逆に辿ります。
スコアだけを知っていることと、どのようにアラインメントされたかを正確に知っていること は異なります。実際のアラインメント（ギャップがどこに入るか）を復元するには、格子を最後から最初に向かって逆探索する必要があります。
各セル dp[i][j] で、3つの候補のうちどれを選んだかを辿ります。
dp[i-1][j-1] + score_pair(s1[i-1], s2[j-1])を選んだ場合 → 対角線方向に移動（マッチ/ミスマッチ）dp[i-1][j] + GAPを選んだ場合 → 上方向に移動（s1の文字とギャップ）dp[i][j-1] + GAPを選んだ場合 → 左方向に移動（ギャップとs2の文字）
def traceback(dp: list[list[int]], s1: str, s2: str) -> tuple[str, str]: aligned1: list[str] = [] aligned2: list[str] = [] i, j = len(s1), len(s2)
while i > 0 or j > 0: current = dp[i][j] if i > 0 and j > 0 and current == dp[i-1][j-1] + score_pair(s1[i-1], s2[j-1]): aligned1.append(s1[i-1]) aligned2.append(s2[j-1]) i -= 1 j -= 1 elif i > 0 and current == dp[i-1][j] + GAP: aligned1.append(s1[i-1]) aligned2.append("-") i -= 1 else: aligned1.append("-") aligned2.append(s2[j-1]) j -= 1
return "".join(reversed(aligned1)), "".join(reversed(aligned2))
a1, a2 = traceback(dp, wild, mut)
# 検証：長さが同じ・ギャップを除去すると元のシーケンスが復元されるassert len(a1) == len(a2)assert a1.replace("-", "") == wildassert a2.replace("-", "") == mut# 挿入した欠落箇所が実際に発見される必要があるassert "-" in a2 # 変異側にギャップが含まれている必要があるこの関数は**反復（イテレーション）**で実装されていますが、同じロジックを再帰でも記述できます。格子のサイズが非常に大きい場合は、反復の方がスタックオーバーフローを防ぐのに有利で、コードの可読性は再帰の方が優れている場合もあります。選択は問題のサイズに依存します。
🔎 再帰と反復：いつどちらを使うか（ドロー — recursion-vs-iteration） 再帰は、問題を自然に表現するときに優雅です。反復は、スタックを積み上げないため安全で、多くの場合高速です。小さな問題・優雅さが必要 → 再帰 / 大きな問題・安定性が必要 → 反復。 このツールの3段階目（格子を埋める）は反復で、4段階目（逆探索）も反復で実装しました。どちらも再帰で実装できますが、シーケンスが長くなるとスタックの危険性があるため、反復を選択します。
🤔 自己説明プロンプト
tracebackで、3つの候補のうち1つ を選択します。2つの候補のスコアが同じ場合（引き分け）、どちらを選ぶべきですか？実際のアラインメント結果は異なる可能性がありますか？（ヒント：複数の最適なアラインメントが存在する可能性があります。）
部品を一つに — 完成された整列クラス
2つの関数を1つのツールにまとめます。
class Aligner: def __init__(self, match=1, mismatch=-1, gap=-2): self.match = match self.mismatch = mismatch self.gap = gap
def _score(self, a: str, b: str) -> int: return self.match if a == b else self.mismatch
def align(self, s1: str, s2: str) -> tuple[str, str, int]: n, m = len(s1), len(s2) dp = [[0] * (m + 1) for _ in range(n + 1)] for i in range(1, n + 1): dp[i][0] = dp[i-1][0] + self.gap for j in range(1, m + 1): dp[0][j] = dp[0][j-1] + self.gap for i in range(1, n + 1): for j in range(1, m + 1): dp[i][j] = max( dp[i-1][j-1] + self._score(s1[i-1], s2[j-1]), dp[i-1][j] + self.gap, dp[i][j-1] + self.gap, ) # 逆トレース a1: list[str] = [] a2: list[str] = [] i, j = n, m while i > 0 or j > 0: if i > 0 and j > 0 and dp[i][j] == dp[i-1][j-1] + self._score(s1[i-1], s2[j-1]): a1.append(s1[i-1]); a2.append(s2[j-1]); i -= 1; j -= 1 elif i > 0 and dp[i][j] == dp[i-1][j] + self.gap: a1.append(s1[i-1]); a2.append("-"); i -= 1 else: a1.append("-"); a2.append(s2[j-1]); j -= 1 return "".join(reversed(a1)), "".join(reversed(a2)), dp[n][m]
aligner = Aligner()a1, a2, score = aligner.align(wild, mut)print(a1)print(a2)print(f"score = {score}")
# 検証：クラスの結果が関数結果と一致することdp = build_score_matrix(wild, mut)assert score == dp[len(wild)][len(mut)]このツールは、今日説明する Needleman-Wunschアルゴリズム の縮小版です。 実際のツール（BLAST、EMBOSSのneedleなど）は、これに洗練されたスコア行列、ギャップオープン/拡張の個別のスコア、ローカルアライメントオプションなどを追加しただけです。 基本的な構造は、あなたが今作成したものと変わりません。
パフォーマンス徹底分析 — なぜこれほど高速化されるのか
import time
def bench(n: int, m: int): s1 = "AT" * (n // 2) s2 = "AG" * (m // 2) t0 = time.time() build_score_matrix(s1, s2) return time.time() - t0
for size in [10, 50, 100, 200]: t = bench(size, size) print(f"長さ {size:4d}: {t*1000:.2f}ms ({size*size} セル)")assert bench(100, 100) < bench(200, 200)数値はマシンによって異なりますが、傾向は常に同じです。セルの数（n×m）に正確に比例します。再帰バージョンが 3^n で爆発するのに対し、DPバージョンは n×m で成長するという決定的な違いがあります。
🔎 Big-Oで整理（図解 — big-o-notation）
- 再帰単純バージョン: O(3^n) — 指数関数。長さ30で既に実用的ではない
- 動的計画法バージョン: O(n·m) — 多項式。長さ10,000同士でも数秒
- メモリ（格子）一つが指数関数を多項式に変換しました。 これが動的計画法の魔法です。
別の道もある（マルチパス推論）
- グローバル vs. ローカルアラインメント: 私たちが作ったのは、最初から最後まで全体のアラインメントを合わせる Needleman-Wunsch（グローバル）。実務では、あるシーケンス内で別のシーケンスの一部を見つける Smith-Waterman（ローカル）をよく使用する。アルゴリズムの基本的な構造はほぼ同じで、スコアリングシステムと初期条件だけが異なる。いつ、何を: 進化系統のシーケンス同士 = グローバル / 遺伝子内でドメインを探す = ローカル。
- BLAST系のヒューリスティック: 私たちのDPは、2つのシーケンスの長さn・mに比例するため、ゲノム（数十億bp）を対象とする場合は依然として大きい。BLASTは、正確な解を諦めて近似解を迅速に得るヒューリスティック（シード → 拡張）を使用する。トレードオフ: 100%最適 vs. 実用的な速度。
- スコアリングマトリックスの高度化: 私たちはmatch/mismatchの二値ルールを使用しているが、実務では、BLOSUM62/PAM250のように、残基のペアごとに異なるスコアを割り当てるマトリックスを使用する。これは、タンパク質シーケンスの化学的類似性を反映している。
- アフィンギャップペナルティ: 私たちはギャップ1つに対して-2を固定しているが、実際には、ギャップを最初につけるコストは大きく、続くギャップは小さいという二重ルール（例：オープン -10、拡張 -1）を使用する。これは、短いギャップを複数持つよりも、長いギャップ1つを持つ方が自然であるという進化的な観測を反映している。
要点: 「可能な部分問題をグリッドに保存して、繰り返し計算を回避する」という原則は、シーケンスアラインメントを超えて、パス検索、編集距離、知識検索など、数多くの場所で再利用される。皆さんが今作った基本的な構造が開く可能性は非常に大きい。
次のステップへ（下部のリンク）
- グリッドを埋める際に、なぜこの順序でセルを訪問するのか → 依存関係と訪問順序
- シーケンスが非常に長い場合のメモリ最適化 → 1D DP圧縮技術
- 先ほど作成したインデックスとの組み合わせ → 応用編 シーケンスDBインデックス
実際に試してみる（独立問題）
- 挿入シナリオ: 野生型
ATGCATGCAと変異型ATGCXATGCA(X = 追加された塩基の挿入) をアライメントし、挿入位置が正確に特定されることをアサートで確認してください。 - スコアの調整:
Aligner(gap=-5)とした場合、結果はどう変わりますか？gap=-0.5の場合はどうですか？ギャップペナルティがアライメントの「性質」をどのように決定するかをまとめてください。 - 複数の最適アライメント: 2つの最適アライメントが存在する配列のペアを作成し、
tracebackを拡張して すべての最適アライメントを返す バージョンを作成してください。 - チャレンジ — メモリ O(min(n,m)): スコアだけが必要な場合は、格子全体ではなく 直前の1行だけを保持する ことで済みます。この最適化されたバージョンを実装してください。（逆追跡を諦める代わりに、非常に長い配列に使用できます。）
まとめ
「2つの配列の最適なアラインメントを見つける」という問題を、以下の3つの要素で解決しました。
- 動的計画法により、指数関数的な計算量（O(3^n)）を多項式時間（O(n・m)）に削減しました。
- 2次元マトリックスが、部分問題の答えを格納する格子となりました。
- 再帰と反復の選択により、格子への値の書き込みと、そこから最適なアラインメントを逆追跡する際の安定性を確保しました。
BLASTが魔法のようにアラインメントの結果を提示してくれるとき、その仕組みが理解できるはずです。それは魔法ではなく、教科書で学んだ動的計画法の格子に配列を配置しているだけです。
この記事は、一般的な教育用例です。実用的な配列アラインメントツール（BLAST/EMBOSS/Bowtieなど）には、洗練されたスコアリングマトリックス、ギャップに対する二重のペナルティ、シードアンドエクステンドヒューリスティック、並列化などが追加されています。詳細な実装は、この基本的な枠組みに自分で追加するか、検証済みのツールに任せることができます。