階層的アセンブリ — k-merグラフを用いて短い断片を繋ぎ合わせる

このトピックを終えると

教科書で学んだ グラフの隣接リスト と DFS/BFS 巡回 を組み合わせて、短いシーケンスリード数百個を1つの長いコンティグに繋ぎ合わせるアセンブラツールを自分で作成できます。SPAdesやVelvetのような実際のアセンブラが魔法のように実行する、そのアセンブリの内部原理、つまりDe Bruijnグラフをコードを通じて理解できるようになります。

この記事は 教育用一般例 です。実際の配列アセンブリは、より複雑な実務パイプラインですが、その中心にあるアルゴリズムの概念を正確に扱っています。

「なぜ3B bpのゲノムが300bpのリードから得られるのか？」—アセンブリの落とし穴

人がヒトゲノムのシーケンスを行うとします。最新のシーケンサー（Illumina）は、一度に約300bpの短いリードを数億個出力します。30億bpのヒトゲノムを、これらの断片から一連の配列に再構築する必要があります。

最も単純なアプローチは、リードを後ろから見て、重なり合う部分を見つけてつなぎ合わせることです。

def naive_assemble(reads: list[str]) -> str:
    result = reads[0]
    for read in reads[1:]:
        overlap = find_overlap(result, read)
        result += read[overlap:]
    return result

このアプローチには、2つの致命的な問題があります。

問題1: リードの順序は意味がありません。シーケンサーから得られたリードは、ゲノム上のランダムな位置から取得されたものであり、reads[0]とreads[1]が実際にゲノム上で隣接している確率は非常に低いです。

問題2: 重なりを計算する処理が指数関数的に増加します。1億個のリードの中から、実際に重なり合う2つのリードをすべて確認するには、最悪の場合、（1億）² = 1兆回の比較が必要です。宇宙が終わるまでに終わることはありません。

実際のアプローチは、グラフを使用することです。リードを短い断片（k-mer）に分割し、各k-merをノードとして、互いに接続する関係をエッジとするド・ブリュイングラフを作成します。次に、このグラフのパスを見つけると、それがアセンブルされた配列になります。

ブラックボックスからパーツへ — De Bruijnグラフを覗いてみる

このフレームを正確に理解するには、3つのパーツが必要です。

パーツ1：k-mer — 文字列を細かく分割する

k-mer は、長さkの短い配列です。ATGCAT を k=3（3-mer）に分割すると、ATG、TGC、GCA、CAT の4つが得られます。

def get_kmers(sequence: str, k: int) -> list[str]:
    return [sequence[i:i+k] for i in range(len(sequence) - k + 1)]

重要な特性は、連続するk-merが(k-1)-merのprefix/suffixを共有するということです。ATG のsuffix TG = TGC のprefix。この共有がグラフのエッジになります。

パーツ2：グラフの隣接リスト

De Bruijnグラフのノードは**(k-1)-merです。エッジはk-mer**です。あるノードから別のノードへエッジがあるということは、その2つの(k-1)-merが、あるk-merのprefix・suffixの関係で繋がっていることを意味します。

このグラフを隣接リスト(辞書)として保存します。

from collections import defaultdict

def build_de_bruijn(reads: list[str], k: int) -> dict[str, list[str]]:
    graph = defaultdict(list)
    for read in reads:
        for kmer in get_kmers(read, k):
            prefix = kmer[:-1]
            suffix = kmer[1:]
            graph[prefix].append(suffix)
    return graph

graph["ATG"] を参照すると、ATGノードから出るエッジの目的地のノードのリストが得られます。

パーツ3：DFS/BFSによる探索

グラフが作成されたら、すべてのエッジを正確に1回ずつ通過するパス(オイラーパス)を見つける必要があります。このパスを繋ぎ合わせると、元の配列を再構成できます。

DFSでアプローチしましょう。開始ノードから再帰的に隣接ノードを訪問し、エッジを消費するたびにグラフから削除します。

def find_eulerian_path(graph: dict[str, list[str]], start: str) -> list[str]:
    stack = [start]
    path = []
    graph = {k: list(v) for k, v in graph.items()}
    
    while stack:
        node = stack[-1]
        if graph.get(node):
            next_node = graph[node].pop()
            stack.append(next_node)
        else:
            path.append(stack.pop())
    
    return path[::-1]

これはHierholzerアルゴリズムです。時間計算量はO(エッジ数)。

3つのパーツを組み合わせて — アセンブリパイプライン

それでは、3つのパーツを1つのパイプラインにまとめましょう。

from collections import defaultdict


def kmerize(sequence: str, k: int) -> list[str]:
    return [sequence[i:i+k] for i in range(len(sequence) - k + 1)]


def build_de_bruijn_graph(reads: list[str], k: int) -> dict[str, list[str]]:
    graph = defaultdict(list)
    for read in reads:
        for kmer in kmerize(read, k):
            prefix, suffix = kmer[:-1], kmer[1:]
            graph[prefix].append(suffix)
    return dict(graph)


def find_start_node(graph: dict[str, list[str]]) -> str:
    out_degree = {node: len(edges) for node, edges in graph.items()}
    in_degree: dict[str, int] = defaultdict(int)
    for edges in graph.values():
        for target in edges:
            in_degree[target] += 1
    
    for node in graph:
        if out_degree[node] > in_degree[node]:
            return node
    return next(iter(graph))


def eulerian_path(graph: dict[str, list[str]], start: str) -> list[str]:
    graph = {k: list(v) for k, v in graph.items()}
    stack, path = [start], []
    while stack:
        node = stack[-1]
        if node in graph and graph[node]:
            stack.append(graph[node].pop())
        else:
            path.append(stack.pop())
    return path[::-1]


def path_to_sequence(path: list[str]) -> str:
    if not path:
        return ""
    return path[0] + "".join(node[-1] for node in path[1:])


def assemble(reads: list[str], k: int = 5) -> str:
    graph = build_de_bruijn_graph(reads, k)
    start = find_start_node(graph)
    path = eulerian_path(graph, start)
    return path_to_sequence(path)

簡単な例で確認してみましょう。

original = "ATGCATGCATGATG"
reads = [original[i:i+8] for i in range(0, len(original) - 7)]
# ['ATGCATGC', 'TGCATGCA', 'GCATGCAT', 'CATGCATG', 'ATGCATGA', 'TGCATGAT', 'GCATGATG']

result = assemble(reads, k=5)
print(result)  # ATGCATGCATGATG (または類似の再構成)

フェーディング — あなたが埋めるべき3つの空白

さあ、あなたの番です。上記のコードには、実際に扱う必要のある3つの空白があります。それぞれにヒントを与えます。

空白1：シークエンシングエラーのフィルタリング

実際のリードには、エラー文字が混ざっています。エラーのあるk-merは、グラフに極端に低い頻度で発生する誤ったノードを作成します。このノードがアセンブリを汚染します。

def filter_low_frequency_kmers(reads: list[str], k: int, threshold: int = 2) -> list[str]:
    from collections import Counter
    kmer_counts: Counter[str] = Counter()
    for read in reads:
        for kmer in kmerize(read, k):
            kmer_counts[kmer] += 1
    
    # TODO: threshold未満の頻度で発生したk-merを含むリードをフィルタリング
    filtered = []
    for read in reads:
        # あなたが埋める場所：このリードのすべてのk-merの頻度がthreshold以上であるか確認
        pass
    return filtered

ヒント: all(kmer_counts[kmer] >= threshold for kmer in kmerize(read, k))

空白2：複数の開始点の処理

実際のシークエンシングデータでは、ゲノムが単一のオイラーパスで完全に繋がらない場合があります。複数のコンティグが生成されることが一般的です。

def assemble_multi_contigs(reads: list[str], k: int) -> list[str]:
    graph = build_de_bruijn_graph(reads, k)
    contigs = []
    
    while graph:
        # TODO: 残りのグラフから開始ノードを選択し、1つのコンティグをアセンブル
        # アセンブル後、使用されたエッジを削除
        # 残りのグラフが空になるまで繰り返す
        pass
    
    return contigs

ヒント: find_start_node を繰り返し呼び出し、各反復で eulerian_path で1つずつアセンブルします。各アセンブルの後、使用されたエッジはすでに削除されているため、自然にグラフが減少します。

空白3：リピート配列への対応

ゲノムにはリピート配列がよくあります。同じ短い配列がゲノムの複数の場所に現れると、De Bruijnグラフでその配列に対応するノードが複数のエッジの交差点となり、アセンブリが曖昧になります。

この問題の完全な解決策は、実際のアセンブラがいくつかの洗練されたヒューリスティックを使用することです（ペアエンドリード、ロングリードなど）。あなたが試すことができるアプローチ：

def detect_repeat_nodes(graph: dict[str, list[str]]) -> set[str]:
    """エッジが複数の方向に分岐するノードは、リピート配列の候補"""
    # TODO: out-degree > 1またはin-degree > 1のノードを返す
    pass

これらのリピートノードがどこにあるかを示すだけで、ユーザーに**「ここにアセンブリが曖昧です」**という警告を与えることができます。

考察 — このコードは実用的なアセンブラとどのように異なるか

あなたが作成したアセンブラは、概念的には SPAdes や Velvet と同じルーツを共有していますが、実際の製品版アセンブラははるかに洗練されています。主な違いをまとめます。

規模への対応: 実際の阿仙ブラは、数億のリードを処理します。あなたの Python 実装は、辞書によるオーバーヘッドのために、この規模を処理できません。実際の阿仙ブラは C++ で記述されており、k-mer を整数としてエンコード (2 ビット/塩基) してメモリを圧縮します。

ペアエンドリード: 実際のシーケンサーは、1 つの断片の両端から 2 つのリードを取得し、その間の距離をほぼ正確に知っています。この距離情報は、反復配列の問題を解決する上で非常に重要です。あなたの実装にはまだこの機能がありません。

長いリード (PacBio/Nanopore): 最近のシーケンサーは、はるかに長いリード (10kb+) を生成します。この場合、アプローチは完全に異なり (OLC — Overlap-Layout-Consensus)、De Bruijn グラフではなく、文字列オーバーラップグラフを使用します。

品質スコア: 実際のリードには、各塩基に品質スコアが関連付けられています。この情報により、エラーフィルタリングをより高度に行うことができます。

拡張プロジェクト

このアセンブラをいくつかの方向に拡張することができます。

1. FASTQファイル入力: 実際のシーケンスデータはFASTQ形式で提供されます。このファイルを解析し、リードと品質スコアを一緒に読み込む関数を追加します。

2. 統計レポート: アセンブリ結果のN50（コンティグ長分布の指標）、総アセンブリ長、コンティグ数などを計算して表示します。

3. 可視化: networkx + matplotlib を使用して、小さなDe Bruijnグラフを描画します。反復ノードがどのようにアセンブリを複雑にするかを目で確認します。

4. 実際のデータで試す: 細菌ゲノム（大腸菌は約4.6M bp）の公開されているシーケンスデータをダウンロードし、作成したアセンブラを実行し、SPAdesの結果と比較します。

この演習における部品の対応表

この演習でどのような個別の CS の概念を組み合わせたかをまとめます。各概念の詳細な説明は DryBench にあります。

• [F] グラフの隣接リスト: graph[node] = [neighbor1, neighbor2, ...] という辞書形式での表現。De Bruijn グラフを格納します。
• [F] DFS/BFS: Hierholzer のスタックベースの DFS を使用してオイラーパスを求めます。実際の組み立て順を巡回します。
• [W] ハッシュテーブル・辞書: defaultdict(list) を使用して O(1) でノードを検索します。
• [W] k-mer スライス: 文字列のスライスとリスト内包表記を使用します。

[F] = 自分で実装する / [W] = 完成したコードとして提供するツール概念。