arboreto-grn-inference

Arboreto GRN 推論

概要

Arboreto は、並列化されたツリーベースの回帰を使用して、遺伝子発現データから遺伝子制御ネットワーク（GRN）を推論します。各標的遺伝子について、ほかのすべての遺伝子（または指定した TF リスト）を特徴量として回帰モデルを学習し、TF・標的・重要度の三つ組を出力します。相互に置き換え可能な 2 つのアルゴリズム、GRNBoost2（勾配ブースティング、高速）と GENIE3（ランダムフォレスト、従来型）を提供しており、入力形式と出力形式は共通です。計算は Dask で並列化され、ノート PC のコア数から HPC クラスターまで拡張できます。

使用する場面

• バルク RNA-seq 発現データから、転写因子と標的遺伝子の制御関係を推論する
• 単一細胞 RNA-seq カウント行列（行が細胞、列が遺伝子）から遺伝子制御ネットワークを構築する
• pySCENIC 制御解析パイプラインの隣接行列（ステップ 1）を生成する
• 実験条件（例：対照群と処置群）間で制御ネットワークの構造を比較する
• 複数の乱数シードを使って推論を実行し、コンセンサス制御ネットワークを作成する
• 同じデータセットで GRNBoost2 と GENIE3 の出力を比較し、GRN の結果を検証する
• 下流のレギュロン同定と活性スコアリングには、arboreto の出力を pySCENIC で使用します
• GRN 推論前の単一細胞データの前処理（QC、正規化、クラスタリング）には、scanpy-scrna-seq を使用します

前提条件

• Python パッケージ：arboreto、pandas、numpy、dask、distributed、scikit-learn、scipy
• データ要件：TSV/CSV 形式の遺伝子発現行列（列が遺伝子、行が観測値）。TF リストファイル（1 行に遺伝子名を 1 つ）を任意で指定できます
• 環境：Python 3.8 以降。可視化には networkx と matplotlib を任意で使用できます

pip install arboreto distributed networkx matplotlib

事前確認

分析コードを書く前に、以下の点をユーザーと確認してください。

decisions:
  - id: D1
    param: algorithm
    kind: required
    source: user
    ask: "Infer the network with gradient boosting (faster) or random forests (the original GENIE3 formulation)?"
    default: "GRNBoost2"

  - id: D2
    param: expressionMatrix
    kind: required
    source: upstream
    ask: "Which processed expression matrix, and which cells or samples within it, should the network be inferred from?"
    default: null

  - id: D3
    param: regulatorList
    kind: required
    source: literature
    ask: "Restrict candidate regulators to a curated transcription-factor list, or let every gene be a candidate regulator?"
    default: "a species-matched TF list"

  - id: D4
    param: edgeThreshold
    kind: required
    source: user
    ask: "The output ranks every regulator-target pair - how many of the top edges should be kept as the network?"
    default: null

  - id: D5
    param: seed
    kind: never_ask
    source: data
    reason: "Fixes the draw for reproducibility; does not change what the data supports"
    default: 0

  - id: D6
    param: daskClient
    kind: never_ask
    source: data
    reason: "Execution backend; affects runtime only"
    default: "local"

D3 は、実用的なネットワークと毛玉状のネットワークを分ける判断です。制御因子を制限しないと、相関のある遺伝子なら何でも制御因子として現れる可能性があり、その場合でも結果には順位付けされた完全なテーブルが出力されます。D4 には安全なデフォルトがありません。アルゴリズムはすべてのペアをスコアリングするため、リストをどこで切るかは、解釈するネットワークについてのユーザーの判断に委ねられます。

クイックスタート

単一のブロックで GRN 推論を完了できます。Dask はマルチプロセス処理によってワーカープロセスを起動するため、if __name__ == '__main__': ガードが必要です。

import pandas as pd
from arboreto.algo import grnboost2
from arboreto.utils import load_tf_names

if __name__ == '__main__':
    # Load expression matrix (observations x genes)
    expression_matrix = pd.read_csv('expression_data.tsv', sep='\t')
    tf_names = load_tf_names('tf_list.txt')  # optional TF filter

    # Infer GRN (uses all local cores by default)
    network = grnboost2(expression_data=expression_matrix,
                        tf_names=tf_names, seed=777)

    # Filter top links and save
    top_network = network[network['importance'] > 1.0]
    top_network.to_csv('grn_output.tsv', sep='\t', index=False, header=False)
    print(f"Inferred {len(network)} links, kept {len(top_network)} above threshold")
    # Example: Inferred 185432 links, kept 12876 above threshold

ワークフロー

ステップ 1：発現データの読み込み

Arboreto は pandas DataFrame（推奨）または NumPy 配列を受け取ります。行は観測値（細胞またはサンプル）、列は遺伝子です。遺伝子名は列ヘッダーに指定する必要があります。

import pandas as pd

# From TSV (genes as columns, observations as rows)
expression_matrix = pd.read_csv('expression_data.tsv', sep='\t')
print(f"Shape: {expression_matrix.shape}  "
      f"({expression_matrix.shape[0]} observations x {expression_matrix.shape[1]} genes)")
# Example: Shape: (5000, 18654)  (5000 observations x 18654 genes)

# From AnnData (e.g., after scanpy preprocessing)
import anndata as ad
adata = ad.read_h5ad('preprocessed.h5ad')
expression_matrix = pd.DataFrame(
    adata.X.toarray() if hasattr(adata.X, 'toarray') else adata.X,
    columns=adata.var_names.tolist()
)
print(f"Converted AnnData: {expression_matrix.shape}")
# Example: Converted AnnData: (5000, 18654)

ステップ 2：転写因子リストの読み込み

TF リストを指定すると、制御因子が既知の転写因子に限定され、計算時間が短縮されるとともに、生物学的な関連性が高まります。省略した場合は、すべての遺伝子が制御因子候補として扱われます。

from arboreto.utils import load_tf_names

# From file (one TF name per line)
tf_names = load_tf_names('human_tfs.txt')
print(f"Loaded {len(tf_names)} transcription factors")
# Example: Loaded 1639 transcription factors

# Or define directly
tf_names = ['MYC', 'TP53', 'SOX2', 'NANOG', 'POU5F1']

# Verify TFs exist in expression matrix columns
tf_in_data = [tf for tf in tf_names if tf in expression_matrix.columns]
print(f"TFs found in expression data: {len(tf_in_data)}/{len(tf_names)}")
# Example: TFs found in expression data: 1583/1639

ステップ 3：Dask クライアントの設定（任意）

デフォルトでは、arboreto はローカルのすべてのコアを使用する内部 Dask クライアントを作成します。リソース管理、モニタリング、またはクラスターへのデプロイには、明示的にクライアントを作成してください。

from distributed import LocalCluster, Client

# Custom local client with resource limits
local_cluster = LocalCluster(
    n_workers=8,
    threads_per_worker=1,   # avoid GIL contention in scikit-learn
    memory_limit='4GB'       # per worker
)
client = Client(local_cluster)
print(f"Dashboard: {client.dashboard_link}")
# Example: Dashboard: http://127.0.0.1:8787/status

ステップ 4：GRN 推論の実行

grnboost2()（推奨）または genie3() を呼び出します。どちらも同じシグネチャと出力形式を使用します。

from arboreto.algo import grnboost2

if __name__ == '__main__':
    network = grnboost2(
        expression_data=expression_matrix,
        tf_names=tf_names,         # 'all' if no TF list
        client_or_address=client,  # omit to use default local scheduler
        seed=777,                  # for reproducibility
        verbose=True               # print progress
    )
    print(f"Inferred {len(network)} regulatory links")
    print(network.head())
    # Example:
    #       TF  target  importance
    # 0   MYC   CDK4       3.214
    # 1   MYC   CCND1      2.871
    # 2  TP53  CDKN1A      2.654

ステップ 5：結果のフィルタリング

生の出力には、重要度がゼロ以外のすべての TF・標的ペアのリンクが含まれます。信頼度の高い制御相互作用を残すためにフィルタリングしてください。

# Strategy 1: Importance threshold
threshold = 1.0
filtered = network[network['importance'] > threshold]
print(f"Threshold {threshold}: {len(filtered)} links "
      f"({len(filtered)/len(network)*100:.1f}% retained)")
# Example: Threshold 1.0: 12876 links (6.9% retained)

# Strategy 2: Top N links per target gene
top_n = 10
top_per_target = (network.groupby('target')
                  .apply(lambda g: g.nlargest(top_n, 'importance'))
                  .reset_index(drop=True))
print(f"Top {top_n} per target: {len(top_per_target)} links")
# Example: Top 10 per target: 186540 links

ステップ 6：ネットワークの保存

TF、標的、重要度の 3 列を含む TSV ファイルとして保存します。

# Save full network
network.to_csv('full_network.tsv', sep='\t', index=False)
print(f"Saved full network: {len(network)} links")

# Save filtered network (without header for pySCENIC compatibility)
filtered.to_csv('filtered_network.tsv', sep='\t', index=False, header=False)
print(f"Saved filtered network: {len(filtered)} links")

# Clean up Dask client if explicitly created
client.close()
local_cluster.close()

ステップ 7：結果の可視化（任意）

上位の制御相互作用を、有向ネットワークグラフとしてプロットします。

import networkx as nx
import matplotlib.pyplot as plt

# Build directed graph from top links
top_links = network.nlargest(50, 'importance')
G = nx.from_pandas_edgelist(
    top_links, source='TF', target='target',
    edge_attr='importance', create_using=nx.DiGraph()
)

# Draw network
fig, ax = plt.subplots(figsize=(12, 10))
pos = nx.spring_layout(G, k=2, seed=42)
nx.draw_networkx_nodes(G, pos, node_size=300, node_color='lightblue', ax=ax)
nx.draw_networkx_labels(G, pos, font_size=7, ax=ax)
edges = nx.draw_networkx_edges(
    G, pos, edge_color=[G[u][v]['importance'] for u, v in G.edges()],
    edge_cmap=plt.cm.Reds, width=1.5, arrows=True, ax=ax
)
plt.colorbar(edges, ax=ax, label='Importance')
ax.set_title('Top 50 Regulatory Interactions')
plt.tight_layout()
plt.savefig('grn_network.png', dpi=300, bbox_inches='tight')
print("Saved grn_network.png")

主なパラメーター

パラメーター	デフォルト	範囲／オプション	効果
expression_data	(必須)	DataFrame または ndarray	発現行列。観測値 × 遺伝子
tf_names	'all'	文字列のリストまたは 'all'	既知の TF のみに制御因子を制限します。'all' ではすべての遺伝子を使用します
gene_names	None	文字列のリスト	expression_data が NumPy 配列の場合に必須です
client_or_address	'local'	Client、文字列、または 'local'	Dask クライアントのインスタンスまたはスケジューラーのアドレス
seed	None	int	再現性を確保するための乱数シードです。結果を再現可能にするには必ず設定してください
verbose	False	bool	推論中に進行状況のメッセージを表示します

主な概念

アルゴリズムの選択ガイド

どちらのアルゴリズムも、同じ重回帰戦略（ターゲット遺伝子ごとにモデルを学習し、特徴量の重要度を抽出する）を採用していますが、基盤となる回帰器が異なります。

特徴	GRNBoost2	GENIE3
手法	早期停止を伴う確率的勾配ブースティング	ランダムフォレスト（または ExtraTrees）
速度	高速 — 大規模データセット向けに最適化されています	低速 — モデルごとのコストが高くなります
メモリー	少ない（早期停止により木の深さを制限）	多い（ターゲットごとに完全なフォレストを作成）
最適な用途	デフォルトの選択肢。1 万件以上の観測値、シングルセル規模	公開済みの GENIE3 結果との比較、検証
インポート	from arboreto.algo import grnboost2	from arboreto.algo import genie3
出力形式	TF、ターゲット、重要度（同一）	TF、ターゲット、重要度（同一）

選択基準：まず GRNBoost2 を使用してください。公開済みの GENIE3 解析を再現する場合、または GRNBoost2 の結果を独立して検証する場合に限り、GENIE3 を使用してください。

出力形式

ネットワークの DataFrame は 3 列で構成され、重要度の降順に並んでいます。

列	型	説明
TF	str	転写因子（制御因子）の遺伝子名
target	str	ターゲット遺伝子名
importance	float	制御重要度スコア（値が高いほど、予測される制御が強い）

重要度スコアは p 値ではありません。木ベースのモデルから得られる特徴量重要度を表します。後続の解析では、絶対しきい値またはターゲットごとの上位 N 件でフィルタリングしてください。

if __name__ == '__main__': が必要な理由

Dask は Python の multiprocessing モジュールを使用してワーカープロセスを起動します。メインガードがないと、起動した各プロセスがスクリプトを再実行し、プロセスが無限に生成されます。スクリプトとして実行する場合は、必ず arboreto の呼び出しを if __name__ == '__main__': で囲んでください。Jupyter ノートブック内では必要ありません。

よく使われるレシピ

レシピ：クラスターでの分散コンピューティング

使用する場面：単一のマシンでは大きすぎるデータセット（例：5 万個以上の細胞、2 万個以上の遺伝子）。クラスター上で Dask の分散スケジューラーが稼働している必要があります。

# On head node: start scheduler
dask-scheduler
# Output: Scheduler at tcp://10.118.224.134:8786

# On each compute node: start workers
dask-worker tcp://10.118.224.134:8786 --nprocs 4 --nthreads 1 --memory-limit 16GB

from distributed import Client
from arboreto.algo import grnboost2
import pandas as pd

if __name__ == '__main__':
    client = Client('tcp://10.118.224.134:8786')
    print(f"Connected to cluster: {client.scheduler_info()['workers'].__len__()} workers")
    # Example: Connected to cluster: 16 workers

    expression_data = pd.read_csv('large_scrna.tsv', sep='\t')
    tf_names = pd.read_csv('tf_list.txt', header=None)[0].tolist()

    network = grnboost2(
        expression_data=expression_data,
        tf_names=tf_names,
        client_or_address=client,
        seed=42, verbose=True
    )
    network.to_csv('cluster_grn.tsv', sep='\t', index=False)
    print(f"Cluster inference complete: {len(network)} links")
    client.close()

レシピ：複数のシードからコンセンサスネットワークを作成

使用する場面：複数の乱数シードでネットワークを集約し、頑健性を高めます。実行間で一貫して出現するリンクほど信頼性が高くなります。

import pandas as pd
from distributed import LocalCluster, Client
from arboreto.algo import grnboost2

if __name__ == '__main__':
    client = Client(LocalCluster(n_workers=8, threads_per_worker=1))
    expression_data = pd.read_csv('expression_data.tsv', sep='\t')
    tf_names = pd.read_csv('tf_list.txt', header=None)[0].tolist()

    seeds = [42, 123, 456, 789, 1001]
    all_networks = []

    for seed in seeds:
        net = grnboost2(expression_data=expression_data,
                        tf_names=tf_names,
                        client_or_address=client, seed=seed)
        net['seed'] = seed
        all_networks.append(net)
        print(f"Seed {seed}: {len(net)} links")

    combined = pd.concat(all_networks, ignore_index=True)

    # Consensus: mean importance across seeds, keep links found in 3+ runs
    consensus = (combined.groupby(['TF', 'target'])
                 .agg(mean_importance=('importance', 'mean'),
                      n_seeds=('seed', 'nunique'))
                 .reset_index())
    consensus = consensus[consensus['n_seeds'] >= 3]
    consensus = consensus.sort_values('mean_importance', ascending=False)
    consensus.to_csv('consensus_network.tsv', sep='\t', index=False)
    print(f"Consensus network: {len(consensus)} links (found in 3+ of {len(seeds)} runs)")
    # Example: Consensus network: 28453 links (found in 3+ of 5 runs)

    client.close()

レシピ：pySCENIC との統合

使用する場面：arboreto による GRN 推論後の制御解析です。Arboreto は隣接行列を生成します（SCENIC のステップ 1）。pySCENIC はレギュロンの特定（ステップ 2）と活性スコアの算出（ステップ 3）を行います。

from arboreto.algo import grnboost2
from arboreto.utils import load_tf_names
import pandas as pd

if __name__ == '__main__':
    # Step 1 (arboreto): GRN inference
    expression_matrix = pd.read_csv('scrna_expression.tsv', sep='\t')
    tf_names = load_tf_names('allTFs_hg38.txt')

    adjacencies = grnboost2(
        expression_data=expression_matrix,
        tf_names=tf_names, seed=777
    )
    adjacencies.to_csv('adjacencies.tsv', sep='\t', index=False, header=False)
    print(f"SCENIC Step 1 complete: {len(adjacencies)} adjacencies")
    # Example: SCENIC Step 1 complete: 234567 adjacencies

    # Steps 2-3 use pySCENIC CLI (requires pyscenic installation):
    # pyscenic ctx adjacencies.tsv hg38_ranking.feather \
    #     --annotations_fname motifs-v10nr.hgnc-m0.001-o0.0.tbl \
    #     --output regulons.csv
    #
    # pyscenic aucell scrna_expression.tsv regulons.csv \
    #     --output auc_matrix.csv
    print("Next: run pyscenic ctx (Step 2) and pyscenic aucell (Step 3)")

想定される出力

• full_network.tsv -- TF、ターゲット、重要度を含む完全な表。列：TF、target、importance
• filtered_network.tsv -- しきい値または上位 N 件のフィルタリング後に得られる高信頼度のサブセット
• consensus_network.tsv -- 複数のシードによる実行結果を集約したネットワーク。列：TF、target、mean_importance、n_seeds
• adjacencies.tsv -- pySCENIC の入力形式に整形された arboreto の出力（ヘッダーなし）
• grn_network.png -- 上位の制御相互作用を示す有向グラフの可視化

トラブルシューティング

問題	原因	解決策
RuntimeError: freeze_support() またはプロセスが無限に起動する	if __name__ == '__main__': ガードがない	arboreto の呼び出しをすべてメインガードで囲みます。Jupyter では不要です
推論中に MemoryError が発生する	発現行列が利用可能な RAM に対して大きすぎる	まず分散の小さい遺伝子を除外し、上位 5-10k 個の遺伝子に絞るか、分散クラスターを使用します
ネットワークが空、またはほとんど空になる	TF 名が発現行列の列名と一致しない	TF 名が expression_matrix.columns と重複していることを確認し、大文字と小文字の区別も確認します
推論が非常に遅い	大規模データセットで GENIE3 を使用している、または Dask ワーカーが少なすぎる	GRNBoost2 に切り替えます。ワーカーを増やした Dask クライアントを明示的に作成します
ModuleNotFoundError: No module named 'arboreto'	パッケージがインストールされていない	pip install arboreto
重要度スコアがすべて 0、またはすべて同じになる	発現行列に分散がない（例：フィルタリング後にすべて 0 になった）	データの品質を確認し、行列に分散のある 0 以外の発現値が含まれていることを確認します
NumPy 配列の入力で TypeError が発生する	gene_names パラメーターがない	配列の列数に一致する gene_names= リストを指定します

同梱リソース

このエントリーは完全に自己完結しており、references/ ディレクトリはありません。

元の参考ファイルの扱い：

• references/algorithms.md（139 行） -- アルゴリズムの比較、パラメーターの詳細、カスタム回帰器の kwargs。アルゴリズム選択ガイド表と主要パラメーターに統合しました。カスタム回帰器の kwargs（regressor_type、regressor_kwargs）は省略しました。高度な調整が必要になることはまれであり、arboreto API ドキュメントに従って scikit-learn の kwargs を直接渡せるためです。
• references/basic_inference.md（152 行） -- 入力形式、TF の読み込み、基本的なワークフロー、出力形式。ワークフローのステップ 1～2、主要概念（出力形式の表）、クイックスタートにすべて統合しました。
• references/distributed_computing.md（242 行） -- ローカルクライアント、クラスターのセットアップ、監視、パフォーマンスのヒント。ワークフローのステップ 3 とレシピ：分散コンピューティングに統合しました。省略した内容：Dask ダッシュボードのパネル説明（対象範囲を超える監視 UI の詳細）、ネットワーク／ストレージの詳細な調整（クラスター管理者向けの内容）。
• scripts/basic_grn_inference.py（98 行） -- argparse を使用する CLI ラッパー。読み込み、推論、保存の中核ロジックはワークフローのステップ 1～6 とクイックスタートに取り込みました。argparse／CLI の定型コードは除外しました。

ユースケースの説明の扱い（元の一般的なユースケースより）：

• シングルセル RNA-seq 解析 -- 「使用する場面」の箇条書きとワークフローのステップ 1／4 に取り込みました
• TF フィルタリングを伴うバルク RNA-seq -- ワークフローのステップ 2 とクイックスタートに取り込みました
• 比較解析（複数条件） -- 独立したレシピとしては省略しました。条件ごとの単純なループであり、コンセンサスネットワークのレシピでパターンを示しているためです

関連スキル

• scanpy-scrna-seq -- GRN 推論前の上流のシングルセル前処理（QC、正規化、HVG 選択）
• scvi-tools-single-cell -- ネットワーク推論前のバッチ補正と補完に用いる深層生成モデル
• pyscenic (planned) -- arboreto の隣接行列を用いた下流のレギュロン同定と細胞活性スコアリング
• networkx-graph-analysis -- 推論された制御ネットワークのグラフ理論的解析（中心性、コミュニティ）

参考資料

• Arboreto GitHub リポジトリ -- ソースコード、README、サンプル
• Arboreto PyPI -- インストールとバージョン履歴
• Moerman et al.（2019）「GRNBoost2 and Arboreto: efficient and scalable inference of gene regulatory networks」Bioinformatics 35(12):2159-2161 -- アルゴリズム論文
• Aibar et al.（2017）「SCENIC: single-cell regulatory network inference and clustering」Nature Methods 14:1083-1086 -- SCENIC パイプライン論文
• Dask 分散処理のドキュメント -- Dask クライアント、クラスターのセットアップ、ダッシュボード