第3章: データ構造とアルゴリズム

3. データ構造とアルゴリズム

学習目標

  • 章の主要概念(データ構造/アルゴリズムの要点)を説明できる
  • 計算量やメモリ効率を意識した実装上の注意を挙げられる
  • サンプルコードや実データに対して適切な手法を選択できる

3.0 アルゴリズム設計の理論的基盤

計算複雑性クラス:

  • P: 多項式時間で解ける決定問題(例: 編集距離がしきい値以下かを判定する問題)
  • NP: 解の候補を多項式時間で検証できる決定問題
  • #P: 解の個数を数える問題(RNA構造数え上げ)
  • PSPACE: 多項式空間(RNA相互作用予測)

近似アルゴリズム理論: 🧪 概念例(擬似コード)

近似比の定義: ρ = max(OPT/ALG, ALG/OPT)
- PTAS (Polynomial-Time Approximation Scheme): 任意のε > 0に対して(1+ε)近似
- FPTAS: さらに計算時間がεに関して多項式
- APX: 定数近似比を持つ問題クラス

確率的アルゴリズム:

  • ラスベガス型: 常に正解、実行時間が確率的
  • モンテカルロ型: 実行時間固定、解の品質が確率的
  • 代表例: MinHashやランダム化探索など、大規模データで近似解を高速に得る手法

文字列アルゴリズムの数学的基礎:

  • 有限オートマトン理論: パターンマッチングの状態機械
  • Z-algorithm: 線形時間での周期性検出
  • KMP法: 失敗関数による効率的検索

3.1 配列データ構造

Suffix Tree:

  • 用途: 高速部分文字列検索
  • 構築時間: O(n)
  • 検索時間: O(m) (mは検索文字列長)
  • メモリ使用量: O(n)

Burrows-Wheeler Transform (BWT):

  • 用途: 配列圧縮と、FM-index などの高速検索構造の基盤
  • 特徴: 後段の圧縮器で高い圧縮効率を得やすい
  • 検索時間: FM-indexと組み合わせると O(m + occ)(occは出現回数)

FM-index: 🧪 概念例(擬似コード)

原配列: ATGCATGC$
BWT変換: CC$GGTTAA
FM-index構築 → 後方検索に利用可能

FM-indexの簡略実装(学習用): 小規模な文字列で後方検索の流れを確認するための簡略版であり、実運用ではより省メモリな実装が用いられる。 🧪 概念例(Python 3.11+)

class FMIndex:
    """FM-index実装(Burrows-Wheeler変換を使用した高速文字列検索)"""

    def __init__(self, text):
        self.text = text + '$'  # 終端記号を追加
        self.sa = self._build_suffix_array()
        self.bwt = self._build_bwt()
        self.alphabet = sorted(set(self.bwt))
        self.occ = self._build_occ_table()
        self.count = self._build_count_table()

    def _build_suffix_array(self):
        """接尾辞配列の構築"""
        return sorted(range(len(self.text)), key=lambda i: self.text[i:])

    def _build_bwt(self):
        """Burrows-Wheeler変換"""
        return ''.join(self.text[i - 1] if i > 0 else self.text[-1] for i in self.sa)

    def _build_occ_table(self):
        """出現回数テーブルの構築"""
        occ = {char: [0] for char in self.alphabet}
        for char in self.bwt:
            for symbol in self.alphabet:
                occ[symbol].append(occ[symbol][-1] + (symbol == char))
        return occ

    def _build_count_table(self):
        """先頭列で各文字が始まる位置を記録"""
        count = {}
        total = 0
        for char in self.alphabet:
            count[char] = total
            total += self.bwt.count(char)
        return count

    def _get_lf_mapping(self, char, index):
        """LF-mappingで前の文字位置へ写像"""
        return self.count[char] + self.occ[char][index]

    def search(self, pattern):
        """検索パターンの出現位置を返す"""
        if not pattern:
            return []

        top = 0
        bottom = len(self.bwt) - 1

        for char in reversed(pattern):
            if char not in self.occ:
                return []

            top = self._get_lf_mapping(char, top)
            bottom = self._get_lf_mapping(char, bottom + 1) - 1

            if top > bottom:
                return []

        return sorted(self.sa[i] for i in range(top, bottom + 1))


fm_index = FMIndex("ATGCATGC")
print(fm_index.search("TGC"))

🧪 概念例(期待される出力)

[1, 5]

実装上の価値: ナイーブな文字列検索はO(nm)の時間計算量を要するが、これらの高度なデータ構造により検索コストを大幅に削減できる。特に参照配列が大きい場合は、インデックスを事前計算しておくことで検索の繰り返しを効率化しやすい。圧縮効果によりメモリ使用量を抑えられる場合もあり、大規模データを扱う実装の基盤として重要である。

3.2 配列アライメントアルゴリズム

動的プログラミングの理論的基盤: 🧪 概念例(擬似コード)

最適性原理 (Bellman): 最適解は最適部分解から構成される
重複部分問題: F(i,j)の値が複数回使用される
メモ化: 計算済み結果の保存による時間短縮

進化距離の数学的モデル:

  • Jukes-Cantor モデル: 全置換が等確率(μt = -3/4 ln(1-4p/3))
  • Kimura 2-parameter モデル: 転移・転換の区別
  • HKY85 モデル: 塩基頻度の偏りを考慮
  • GTR モデル: 6つの置換率パラメータ

統計的有意性評価: 🧪 概念例(擬似コード)

E-value = Kmn exp(-λS)
- K: 検索空間の補正定数
- m,n: 配列長
- λ: Gumbel分布のパラメータ
- S: アライメントスコア

グローバルアライメント(Needleman-Wunsch法): 🧪 概念例(擬似コード)

動的プログラミング:
F(i,j) = max{
    F(i-1,j-1) + s(xi, yj),  // マッチ/ミスマッチ
    F(i-1,j) + gap,          // 挿入
    F(i,j-1) + gap           // 削除
}

初学者向け最小例(Python 3.11+ / 標準ライブラリのみ): Needleman-Wunsch法と同じ動的計画法の骨格を理解するため、まずは2本の短い配列の編集距離を求める。 🧪 概念例(Python 3.11+)

def edit_distance(seq1: str, seq2: str) -> tuple[int, list[list[int]]]:
    rows = len(seq1) + 1
    cols = len(seq2) + 1
    dp = [[0] * cols for _ in range(rows)]

    for i in range(rows):
        dp[i][0] = i
    for j in range(cols):
        dp[0][j] = j

    for i, base1 in enumerate(seq1, start=1):
        for j, base2 in enumerate(seq2, start=1):
            substitution_cost = 0 if base1 == base2 else 1
            dp[i][j] = min(
                dp[i - 1][j] + 1,                      # 削除
                dp[i][j - 1] + 1,                      # 挿入
                dp[i - 1][j - 1] + substitution_cost, # 置換または一致
            )

    return dp[-1][-1], dp


seq1 = "GATTACA"
seq2 = "GACTATA"

distance, matrix = edit_distance(seq1, seq2)
print(f"{seq1}{seq2} の編集距離: {distance}")
print("DP行列の最終行:", matrix[-1])

🧪 概念例(期待される出力)

GATTACA と GACTATA の編集距離: 2
DP行列の最終行: [7, 6, 5, 4, 4, 3, 3, 2]

編集距離が小さいほど、2配列を一致させるための挿入・削除・置換の回数が少ないことを意味する。まずは短い配列でDP行列の更新規則を確認し、その後にスコア行列やギャップペナルティを導入すると、Needleman-Wunsch法やSmith-Waterman法の理解につながる。

ローカルアライメント(Smith-Waterman法):

  • グローバル法を局所最適化に拡張
  • 0以下の値をリセット → 局所的類似領域を検出

高速近似法(BLAST):

  1. シード検索: 短い完全一致領域を特定
  2. 拡張: シード周辺の類似性評価
  3. 統計的有意性検定

スコア行列の理論的基礎: 🧪 概念例(擬似コード)

PAM行列 (Point Accepted Mutation):
- 進化距離1 PAMは1%のアミノ酸変化
- PAM250: 平均2.5回の置換(遠縁関係)

BLOSUM行列 (BLOcks SUbstitution Matrix):
- 保存ブロック内の実際の置換頻度から導出
- BLOSUM62: 62%以上の同一性を持つ配列から構築

ギャップペナルティの最適化:

  • 線形ギャップ: gap penalty = g × length
  • アフィンギャップ: penalty = open + extend × (length-1)
  • 凸ギャップ: 生物学的により現実的なモデル

応用における重要性: 配列アライメントは生物学的類似性の定量化手法である。進化関係の推定、機能推定、疾患原因変異の特定に不可欠である。厳密解法は計算量がO(nm)で大規模データに適用困難だが、BLASTによる近似解法により実用的な解析時間を実現している。これにより、未知配列の機能を既知データベースとの比較により数分で推定可能となる。

3.3 グラフアルゴリズム

グラフ理論の基礎概念: 🧪 概念例(擬似コード)

グラフG = (V, E): 頂点集合Vと辺集合E
- 次数分布: P(k) = ネットワーク中の次数kのノード割合
- クラスタ係数: C = 3×三角形数 / 連結三組数
- 最短経路長: 全ノード対間の最短距離の平均
- 中心性: 重要ノードの定量的指標

生物学的ネットワークの位相的特性:

  • スケールフリー性: P(k) ∝ k^(-γ), γ ≈ 2-3
  • スモールワールド性: 高クラスタ係数 + 短い特性経路長
  • 階層構造: モジュール内高密度、モジュール間低密度
  • 堅牢性: ランダム攻撃に強く、標的攻撃に脆弱

ネットワーク中心性指標: 🧪 概念例(擬似コード)

次数中心性: CD(v) = deg(v) / (n-1)
近接中心性: CC(v) = (n-1) / Σ d(v,u)
媒介中心性: CB(v) = Σ σst(v) / σst
固有ベクトル中心性: Av = λv (主固有ベクトル)
PageRank: PR(v) = (1-d)/n + d Σ PR(u)/deg(u)

コミュニティ検出アルゴリズム:

  • Modularity最適化: Q = 1/2m Σ [Aij - kikj/2m]δ(ci,cj)
  • Louvain法: 階層的クラスタリング
  • Label Propagation: ラベル伝播による高速クラスタリング
  • スペクトラル法: ラプラシアン行列の固有ベクトル

タンパク質相互作用ネットワーク:

  • ノード: タンパク質
  • エッジ: 相互作用関係
  • 問題: 最短経路、クラスタリング、中心性解析

代謝パスウェイ解析:

  • 有向グラフ: 代謝反応の流れ
  • 重み: 反応速度定数
  • 最適化: フラックスバランス解析(FBA)

フラックスバランス解析の数学的定式化: 🧪 概念例(擬似コード)

線形計画問題:
max cᵀv
subject to: Sv = 0 (定常状態制約)
           lb ≤ v ≤ ub (反応速度制約)

S: 化学量論行列
v: フラックスベクトル
c: 目的関数係数

ネットワーク解析の計算複雑性:

  • 最短経路: Dijkstra O((V+E)logV), Floyd-Warshall O(V³)
  • 最大フロー: Ford-Fulkerson O(E f* ), Dinic O(V²E)
  • 最小カット: Karger O(V²logV), Stoer-Wagner O(VE+V²logV)
  • TSP (Traveling Salesman): NP困難、近似比3/2

システム医学への貢献: 生物学的システムはネットワーク構造を持つ。グラフ理論により、疾患は「ネットワークの摂動」として定量的に解析できる。薬剤標的の優先順位付け、副作用予測、多剤併用効果の予測が可能となる。例えば、中心性の高いノードを標的とすることで、より効果的な治療戦略を設計できる。システム全体の堅牢性解析により、耐性機構の理解も深まる。

3.4 シングルセルデータ構造

スパース行列表現: 🧪 概念例(擬似コード)

# 細胞×遺伝子発現行列の効率的表現
import scipy.sparse as sp

class SingleCellDataStructure:
    def __init__(self, expression_matrix):
        # CSR形式でのスパース行列保存
        self.data = sp.csr_matrix(expression_matrix)
        self.n_cells, self.n_genes = self.data.shape

k-NN グラフ構造:

  • 細胞間類似性の効率的計算
  • Approximate Nearest Neighbor (ANN)アルゴリズム
  • LSH (Locality Sensitive Hashing)による高速化

軌跡推定のための構造:

  • 最小全域木(MST)
  • 主曲線(Principal Curves)
  • 拡散マップ(Diffusion Maps)

3.5 集団遺伝学のアルゴリズム

ハプロタイプフェージング:

  • 隠れマルコフモデル(HMM)によるアプローチ
  • Li-Stephensモデル
  • SHAPEIT, Eagle等の実装

連鎖不平衡(LD)計算: 🧪 概念例(擬似コード)

r² = (pAB - pA×pB)² / (pA×pa×pB×pb)
D' = D / Dmax

集団構造解析:

  • 主成分分析(PCA)による祖先成分推定
  • ADMIXTURE/STRUCTUREアルゴリズム
  • f統計量(FST, f3, f4)

最小入出力(期待成果物/期待ログ)

  • 入力: 小規模な配列データ(FASTA等)と、簡単な行列/グラフ(トイデータ)
  • 出力(期待成果物): 代表アルゴリズムの実装(例: 編集距離/アライメント/グラフ探索)と、計算量の見積もり・実測の比較メモ
  • 期待ログ(例): 入力サイズ・実行時間・メモリの計測結果が残る(再現できる形で記録)

前へ: 計算インフラストラクチャ 目次 次へ: ゲノム解析技術

演習

  1. 章内のサンプルコードをベースに、入力規模を10倍にした際の実行時間/メモリを測定し、ボトルネックを特定せよ。
  2. 同じ課題を別手法(データ構造/アルゴリズム)で解き、計算量/実行時間の差を比較して考察せよ。

具体課題例

  • 公開データを用いた再現(SRA/GEO/ArrayExpressから実在アクセッションを選定し)、前処理→主解析→結果要約まで実施。
  • 代替ツールの比較(例: ツールA vs ツールB)。処理時間/メモリ/精度など評価指標を定義し、比較表を作成。
  • 成果物一式(レポート、使用コマンド/パラメータ、ツール/ライブラリのバージョン、入力/出力、実行ログ、MultiQC等のレポート、図表)を添付。