PythonでIBDの疾患モジュールを解析する：既存薬からリポジショニング候補を絞り込むまで（前半無料）

Disease Module入門 #1では疾患モジュールの考え方を、#2ではシード・ネットワーク・検定の設計を、#3では IBD を例に結果の読み方を解説しました。本記事はその続編です。自分の疾患で同じ解析を走らせるための Python 実装を、入力データの取得からランダム化検定、薬剤近接性、機能解釈まで一気通貫で扱います。

題材は炎症性腸疾患（IBD）です。ただし、ここで作るパイプラインは疾患IDと候補薬剤の表を差し替えれば、ほかの疾患にも使えます。目標は、きれいなネットワーク図を作ることではありません。疾患遺伝子が本当にネットワーク上で凝集しているかを検定し、薬剤候補を仮説として順位づけ、その限界まで説明できる状態になることです。

この解析で得られるもの

出力は、単なる遺伝子リストではありません。どの疾患シードがネットワークに載ったか、シードの最大連結成分がランダム期待よりどの程度大きいか、候補薬の標的が疾患モジュールへどの程度近いか、どの生物学的機能に収束するかを、すべて表と図で残します。解析の各段階に入力データの版と乱数シードを記録するため、あとから疾患ID・閾値・ネットワークを替えた感度分析にもつなげられます。

今回のネットワークには STRING v12.0 のヒト物理相互作用を使います。共発現や文献共起を含む統合スコアではなく、物理相互作用に限定することで、「ネットワーク上の距離」を解釈しやすくします。STRING のスコアは 0〜1000 なので、本記事では高信頼の combined_score >= 700 を採用します。

解析の全体像

Open Targets（疾患関連遺伝子）       STRING physical PPI
             │                               │
             └──── IDをそろえてネットワークへ投影 ────┐
                                                       ↓
                         最大連結成分・最近接距離・次数一致ランダム化
                                                       ↓
                         疾患モジュールの有意性（z-score / p値）
                                                       ↓
Open Targets / 手元の候補薬剤 ── 薬剤標的 ── closest proximity ── 候補順位
                                                       ↓
                              g:Profiler による GO / Reactome / KEGG 解釈

環境を作る

Python 3.11 以上を前提にします。解析コードは、#3の IBD ケーススタディで使った実装を、現行のデータソースで動くように整理したものです。公開データベースは更新されるため、記事中の数値をそのまま再現できるとは限りません が、同じ前提から同じ結論に到達できることは確認しています。まず仮想環境と作業ディレクトリを作ります。

mkdir disease-module && cd disease-module
python -m venv .venv
source .venv/bin/activate

python -m pip install --upgrade pip
pip install pandas numpy networkx matplotlib requests gprofiler-official

mkdir -p data out

次に、STRING のダウンロードページで Homo sapiens に絞った protein.physical.links.v12.0 と protein.info.v12.0 を取得します。下のURLは種別を絞ったファイルです。前者は約9 MB、後者は約2 MBなので、全生物種を含む巨大なファイルを取る必要はありません。

curl -L -o data/9606.protein.physical.links.v12.0.txt.gz \
  https://stringdb-downloads.org/download/protein.physical.links.v12.0/9606.protein.physical.links.v12.0.txt.gz

curl -L -o data/9606.protein.info.v12.0.txt.gz \
  https://stringdb-downloads.org/download/protein.info.v12.0/9606.protein.info.v12.0.txt.gz

Open Targets から IBD シードを取得する

Open Targets の現行 Platform API では、IBD は MONDO_0005265 として取得できます。

post('query($id:String!){ disease(efoId:$id){ id name } }', {"id": DISEASE_ID})

以下のスクリプトは関連ターゲットを全ページ取得し、スコアが 0.5 以上のシードを data/seeds.tsv に保存します。DISEASE_ID を変えれば、ほかの疾患を同じ形式で入力できます。

# 00_fetch_seeds.py
from pathlib import Path
import requests
import pandas as pd

ENDPOINT = "https://api.platform.opentargets.org/api/v4/graphql"
DISEASE_ID = "MONDO_0005265"  # inflammatory bowel disease
SCORE_CUTOFF = 0.50

QUERY = """
query AssociatedTargets($efoId: String!, $index: Int!, $size: Int!) {
  disease(efoId: $efoId) {
    associatedTargets(page: {index: $index, size: $size}) {
      count
      rows { target { approvedSymbol id } score }
    }
  }
}
"""


def post(query, variables):
    response = requests.post(ENDPOINT, json={"query": query, "variables": variables}, timeout=60)
    response.raise_for_status()
    payload = response.json()
    if payload.get("errors"):
        raise RuntimeError(payload["errors"])
    return payload["data"]


def fetch_all_targets(disease_id, page_size=500):
    rows, index, total = [], 0, None
    while total is None or len(rows) < total:
        data = post(QUERY, {"efoId": disease_id, "index": index, "size": page_size})
        page = data["disease"]["associatedTargets"]
        total = page["count"]
        if not page["rows"]:  # 空ページで打ち切る（無限ループ防止）
            break
        rows.extend(page["rows"])
        index += 1
    return rows


rows = fetch_all_targets(DISEASE_ID)
seeds = pd.DataFrame([
    {"symbol": row["target"]["approvedSymbol"], "ensembl_id": row["target"]["id"], "ot_score": row["score"]}
    for row in rows if row["target"]["approvedSymbol"] and row["score"] >= SCORE_CUTOFF
]).sort_values("ot_score", ascending=False)

Path("data").mkdir(exist_ok=True)
seeds.to_csv("data/seeds.tsv", sep="\t", index=False)
print(f"saved {len(seeds)} seeds with Open Targets score >= {SCORE_CUTOFF}")

STRING 物理PPIネットワークを作る

次に STRING ID を遺伝子シンボルに対応づけ、スコア700以上の物理相互作用だけからグラフを作ります。最大連結成分だけを残すのは、到達不能なノードを含めたまま最短距離を計算しないためです。

# 01_build_network.py
import gzip
import pickle
import networkx as nx

SCORE_CUTOFF = 700

string_to_symbol = {}
with gzip.open("data/9606.protein.info.v12.0.txt.gz", "rt") as handle:
    next(handle)
    for line in handle:
        string_id, symbol, *_ = line.rstrip("\n").split("\t")
        if symbol:
            string_to_symbol[string_id] = symbol

graph = nx.Graph()
with gzip.open("data/9606.protein.physical.links.v12.0.txt.gz", "rt") as handle:
    next(handle)
    for line in handle:
        protein_a, protein_b, score = line.split()
        if int(score) < SCORE_CUTOFF:
            continue
        gene_a = string_to_symbol.get(protein_a)
        gene_b = string_to_symbol.get(protein_b)
        if gene_a and gene_b and gene_a != gene_b:
            graph.add_edge(gene_a, gene_b, score=int(score))

lcc = max(nx.connected_components(graph), key=len)
graph = graph.subgraph(lcc).copy()

with open("out/string_physical_v12_score700.pkl", "wb") as handle:
    pickle.dump(graph, handle)

print(f"nodes={graph.number_of_nodes():,}, edges={graph.number_of_edges():,}")

ここまでで、再現可能な入力とネットワークが揃いました。次からは、シードがこのネットワーク上で本当に凝集しているかを検定し、薬剤標的との近さを測り、リポジショニング候補を絞り込んでいきます。

先に概念を押さえたい方へ

本記事は実装に絞っています。疾患モジュールとは何か、シードやネットワークをどう設計するか、結果をどう解釈するか は、無料の3記事にまとめてあります。手を動かす前に理屈を確認したい場合は、こちらからどうぞ。

関連記事Disease Module入門 #1：疾患をネットワークの局所的な破綻として読む →

関連記事Disease Module入門 #2：シード・ネットワーク・検出法と統計的検証 →

関連記事Disease Module入門 #3：創薬応用とIBDで実際に計算する →

シードの取得元である Open Targets の API の使い方は、こちらで解説しています。

関連記事Open Targets入門 #3：Python + GraphQL APIでデータを取得する →