arboreto-grn-inference

Arboreto GRN 推断

概述

Arboreto 使用并行化的树模型回归，根据基因表达数据推断基因调控网络（GRN）。对于每个目标基因，它会使用所有其他基因（或指定的 TF 列表）作为特征训练回归模型，并输出 TF-目标基因-重要性三元组。它提供两种可互换的算法——GRNBoost2（梯度提升，速度快）和 GENIE3（随机森林，经典算法）——输入和输出格式完全相同。计算通过 Dask 并行化，可从笔记本电脑的 CPU 核心扩展到 HPC 集群。

适用场景

• 根据 bulk RNA-seq 表达数据推断转录因子到目标基因的调控关系
• 根据单细胞 RNA-seq 计数矩阵构建基因调控网络（行是细胞，列是基因）
• 生成 pySCENIC 调控分析工作流的邻接矩阵（步骤 1）
• 比较不同实验条件下的调控网络结构（例如，对照组与处理组）
• 通过多个随机种子运行推断，生成共识调控网络
• 在同一数据集上比较 GRNBoost2 和 GENIE3 的输出，验证 GRN 结果
• 若要进行下游 regulon 识别和活性评分，请将 arboreto 输出与 pySCENIC 配合使用
• 若要在 GRN 推断前进行单细胞预处理（QC、归一化、聚类），请使用 scanpy-scrna-seq

前置条件

• Python 软件包：arboreto、pandas、numpy、dask、distributed、scikit-learn、scipy
• 数据要求：TSV/CSV 格式的基因表达矩阵（列是基因，行是观测值）；可选的 TF 列表文件（每行一个基因名称）
• 环境：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 推断。if __name__ == '__main__': 保护语句是必需的，因为 Dask 会通过 multiprocessing 启动工作进程。

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 会创建一个使用所有本地 CPU 核心的内部 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：保存网络

保存为包含三列的 TSV 文件：TF、目标基因、重要性。

# 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	整数	用于保证结果可复现的随机种子；要获得可复现结果，请务必设置
verbose	False	布尔值	推断期间打印进度消息

关键概念

算法选择指南

这两种算法都采用相同的多元回归策略（为每个目标基因训练一个模型，并提取特征重要性），但底层回归器不同。

特征	GRNBoost2	GENIE3
方法	带提前停止的随机梯度提升	随机森林（或 ExtraTrees）
速度	快 — 针对大型数据集进行了优化	较慢 — 单个模型的开销更高
内存	较低（提前停止限制了树的深度）	较高（每个目标都构建完整森林）
适用场景	默认选择；1 万个以上观测值，单细胞规模	与已发表的 GENIE3 结果对比；验证
导入	from arboreto.algo import grnboost2	from arboreto.algo import genie3
输出格式	TF、目标、重要性（相同）	TF、目标、重要性（相同）

选择规则：从 GRNBoost2 开始。仅在重现已发表的 GENIE3 分析或独立验证 GRNBoost2 结果时使用 GENIE3。

输出格式

网络 DataFrame 包含三列，按重要性降序排列：

列	类型	说明
TF	字符串	转录因子（调控因子）基因名称
target	字符串	目标基因名称
importance	浮点数	调控重要性分数（越高表示预测的调控作用越强）

重要性分数不是 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 执行 regulon 识别（第 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 调用放在 main 保护中；在 Jupyter 中无需此操作
推断期间出现 MemoryError	表达矩阵对于可用 RAM 来说过大	先筛除低方差基因，将基因数缩减至前 5-10k 个，或使用分布式集群
网络为空或几乎为空	TF 名称与表达矩阵的列名不匹配	确认 TF 名称与 expression_matrix.columns 有重叠；检查大小写是否一致
推断速度非常慢	在大型数据集上使用 GENIE3，或 Dask worker 数量过少	切换到 GRNBoost2；创建包含更多 worker 的显式 Dask 客户端
ModuleNotFoundError: No module named 'arboreto'	未安装软件包	pip install arboreto
所有重要性分数均为 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 筛选的 Bulk RNA-seq——已融入工作流步骤 2 和快速入门
• 比较分析（多个条件）——未作为独立操作指南保留：只需对条件进行简单循环，相关模式已在共识网络操作指南中展示

相关技能

• scanpy-scrna-seq——GRN 推断前的上游单细胞预处理（QC、标准化、HVG 选择）
• scvi-tools-single-cell——网络推断前用于批次校正和插补的深度生成模型
• pyscenic (planned)——使用 arboreto 邻接关系进行下游 regulon 识别和细胞活性评分
• 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 客户端、集群设置、仪表板