跳转至

基因调控网络推断

一句话概述

利用SCENIC/pySCENIC、GRNBoost2和GENIE3等算法从单细胞或bulk RNA-seq数据推断转录因子(TF)及其靶基因的调控关系，构建基因调控网络(GRN)，是揭示细胞命运决定机制和疾病调控失调的核心分析方法。

核心知识点表格

知识点关键内容常用工具
GRN概念TF→靶基因的有向调控关系-
共表达推断从表达相关性推断调控GENIE3, GRNBoost2
motif验证用TF结合位点验证调控关系RcisTarget, SCENIC
regulon活性每个细胞中TF regulon的活性评分AUCell
SCENIC流程共表达→motif筛→regulon活性SCENIC, pySCENIC
差异调控条件间TF活性变化SCENIC+, CellOracle
因果推断扰动数据验证因果关系Perturb-seq, CRISPRi

各步骤详解

第一步：理解基因调控网络

白话解释： 基因调控网络描述的是"谁控制谁"的关系：哪些转录因子控制哪些基因的表达。就像公司的组织架构图——经理(TF)管理一群员工(靶基因)，这一群人就是一个"regulon"。SCENIC流程通过三步推断这个网络：1）找共表达关系；2）用motif信息验证TF是否真的能结合靶基因；3）计算每个细胞中各TF的活性。

技术细节： - Regulon = TF + 其直接靶基因集合 - SCENIC三步: GRNBoost2(共表达)→RcisTarget(motif)→AUCell(活性) - GRNBoost2: 随机梯度提升树推断TF-gene关联 - RcisTarget: 检查靶基因启动子区是否富集TF的结合motif - AUCell: 用AUC评分量化每个细胞中regulon的活性

第二步：pySCENIC完整流程

代码示例：

# 安装
# pip install pyscenic

import scanpy as sc
import loompy
from pyscenic.rnkdb import FeatherRankingDatabase as RankingDatabase
from pyscenic.utils import modules_from_adjacencies
from pyscenic.prune import prune2df, df2regulons
from pyscenic.aucell import aucell
import pandas as pd
import numpy as np

# 1. 准备数据
adata = sc.read_h5ad("scRNA_processed.h5ad")

# 过滤：保留在至少1%细胞中表达的基因
sc.pp.filter_genes(adata, min_cells=int(0.01 * adata.n_obs))

# 导出为loom格式
row_attrs = {"Gene": np.array(adata.var_names)}
col_attrs = {"CellID": np.array(adata.obs_names),
             "CellType": np.array(adata.obs["cell_type"])}
loompy.create("input.loom", adata.X.T, row_attrs, col_attrs)

# 2. 准备数据库文件（需要预下载）
# TF列表: https://resources.aertslab.org/cistarget/tf_lists/
# Ranking数据库: https://resources.aertslab.org/cistarget/databases/
# pySCENIC 0.12.0+ 要求 Feather v2 格式数据库（需 ctxcore >= 0.2），旧格式 Feather v1 不兼容
# Motif注释: https://resources.aertslab.org/cistarget/motif2tf/

# 3. 命令行运行pySCENIC（推荐，更稳定）


# pySCENIC命令行流程（推荐方式）

# Step 1: GRNBoost2 共表达网络推断
pyscenic grn \
    input.loom \
    /db/allTFs_hg38.txt \           # TF列表
    -o adjacencies.tsv \
    --method grnboost2 \
    --num_workers 16 \
    --seed 42

# Step 2: RcisTarget motif pruning (cisTarget)
pyscenic ctx \
    adjacencies.tsv \
    /db/hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather \
    /db/hg38_500bp_up_100bp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather \
    --annotations_fname /db/motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl \
    --expression_mtx_fname input.loom \
    --output regulons.csv \
    --mask_dropouts \
    --num_workers 16

# Step 3: AUCell regulon活性评分
pyscenic aucell \
    input.loom \
    regulons.csv \
    -o scenic_output.loom \
    --num_workers 16

# 输出: scenic_output.loom 包含原始数据 + regulon AUC scores

# 4. 解析和可视化结果
import loompy

# 读取SCENIC结果
lf = loompy.connect("scenic_output.loom")

# 获取AUCell矩阵
auc_mtx = pd.DataFrame(
    lf.ca.RegulonsAUC,
    index=lf.ca.CellID,
    columns=lf.ra.Regulons)

# 关闭loom文件
lf.close()

# 或使用scanpy读取
adata_scenic = sc.read_loom("scenic_output.loom")

# 5. 可视化TF活性
# 在UMAP上展示特定TF活性
import matplotlib.pyplot as plt

fig, axes = plt.subplots(2, 3, figsize=(15, 10))
tfs_to_show = ["PAX6(+)", "SOX2(+)", "NEUROD1(+)", 
               "GATA3(+)", "FOXP3(+)", "TCF7(+)"]

for ax, tf in zip(axes.flat, tfs_to_show):
    if tf in auc_mtx.columns:
        sc.pl.umap(adata, color=auc_mtx[tf].values, 
            ax=ax, title=tf, show=False)
plt.tight_layout()
plt.savefig("tf_activity_umap.pdf")

# 6. 差异TF活性分析
from scipy.stats import mannwhitneyu

cell_types = adata.obs["cell_type"]
diff_regulons = []

for tf in auc_mtx.columns:
    for ct in cell_types.unique():
        mask = cell_types == ct
        stat, pval = mannwhitneyu(
            auc_mtx.loc[mask, tf],
            auc_mtx.loc[~mask, tf],
            alternative="greater")
        if pval < 0.001:
            diff_regulons.append({"TF": tf, "CellType": ct, "pvalue": pval})

diff_df = pd.DataFrame(diff_regulons)
print(diff_df.sort_values("pvalue").head(20))

第三步：SCENIC R版本

代码示例：

# 安装
# BiocManager::install(c("SCENIC","AUCell","RcisTarget"))
library(SCENIC)
library(AUCell)
library(RcisTarget)
library(Seurat)

# 1. 准备数据
seurat_obj <- readRDS("seurat_processed.rds")
exprMat <- GetAssayData(seurat_obj, slot="counts")
exprMat <- as.matrix(exprMat)
cellInfo <- data.frame(
    CellType = seurat_obj$cell_type,
    row.names = colnames(seurat_obj))

# 2. 初始化SCENIC设置
org <- "hgnc"  # 人类; "mgi"=小鼠
dbDir <- "/db/scenic/"  # 数据库目录
myDatasetTitle <- "My_SCENIC_Analysis"

scenicOptions <- initializeScenic(
    org = org,
    dbDir = dbDir,
    # mc9nr 为旧版数据库，当前推荐使用 mc_v10_clust 系列（motif 数量和注释质量更高）
    dbs = c("hg38__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather",
            "hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather"),
    datasetTitle = myDatasetTitle,
    nCores = 16)

# 3. 基因过滤和共表达模块
genesKept <- geneFiltering(exprMat, scenicOptions,
    minCountsPerGene = 3 * 0.01 * ncol(exprMat),
    minSamples = ncol(exprMat) * 0.01)
exprMat_filtered <- exprMat[genesKept, ]

# 4. 计算相关性矩阵
runCorrelation(exprMat_filtered, scenicOptions)

# 5. GRNBoost2推断共表达网络
# 通常在Python中运行GRNBoost2（更快）
# 或使用GENIE3:
runGenie3(exprMat_filtered, scenicOptions, nParts = 20)

# 6. 构建和修剪regulons (cisTarget)
scenicOptions <- runSCENIC_1_coexNetwork2modules(scenicOptions)
scenicOptions <- runSCENIC_2_createRegulons(scenicOptions)

# 7. AUCell活性评分
scenicOptions <- runSCENIC_3_scoreCells(scenicOptions, exprMat_filtered)

# 8. 二值化（可选：判断regulon在每个细胞中是否"开启"）
scenicOptions <- runSCENIC_4_aucell_binarize(scenicOptions)

# 9. 提取结果
regulonAUC <- getAUC(loadInt(scenicOptions, "aucell_regulonAUC"))
regulons <- loadInt(scenicOptions, "regulons")

# 10. 可视化
# TF活性热图
library(ComplexHeatmap)
avg_auc <- t(apply(regulonAUC, 1, function(x) 
    tapply(x, cellInfo$CellType, mean)))

Heatmap(avg_auc, name="AUC",
    show_row_names=TRUE, row_names_gp=gpar(fontsize=6))

第四步：GENIE3独立使用

白话解释： GENIE3使用随机森林回归来推断每个基因的表达受哪些TF调控。对于每个靶基因，它训练一个随机森林模型用所有TF的表达来预测该基因的表达，根据特征重要性排序得到最可能的调控TF。

代码示例：

library(GENIE3)

# 准备数据
expr_matrix <- as.matrix(GetAssayData(seurat_obj, "data"))

# TF列表
tf_list <- read.table("/db/allTFs_hg38.txt")$V1
tf_list <- intersect(tf_list, rownames(expr_matrix))

# 运行GENIE3
# 注意：计算密集，建议子采样或使用GRNBoost2替代
set.seed(42)
# 子采样细胞（加速）
cells_subset <- sample(colnames(expr_matrix), min(5000, ncol(expr_matrix)))
expr_subset <- expr_matrix[, cells_subset]

weightMat <- GENIE3(
    exprMatrix = expr_subset,
    regulators = tf_list,
    nCores = 16,
    verbose = TRUE)

# 提取Top边
link_list <- getLinkList(weightMat, threshold = 0.01)
# 格式: regulatoryGene, targetGene, weight
head(link_list)

# 筛选高置信度边
top_links <- link_list[link_list$weight > quantile(link_list$weight, 0.99), ]
cat("高置信度调控关系数:", nrow(top_links), "\n")

第五步：GRN可视化和解读

代码示例：

library(igraph)
library(ggraph)

# 构建网络
# 只取top TF及其靶基因
top_tfs <- names(sort(table(top_links$regulatoryGene), decreasing=TRUE))[1:20]
subnet <- top_links[top_links$regulatoryGene %in% top_tfs, ]
subnet <- subnet[1:200, ]  # 限制边数

# 创建igraph对象
g <- graph_from_data_frame(subnet[,1:2], directed=TRUE)
V(g)$type <- ifelse(V(g)$name %in% tf_list, "TF", "Target")
V(g)$size <- ifelse(V(g)$type == "TF", 10, 3)
V(g)$color <- ifelse(V(g)$type == "TF", "red", "lightblue")

# 绘制网络
pdf("GRN_network.pdf", width=12, height=12)
plot(g, layout=layout_with_fr(g),
    vertex.label.cex=0.6,
    edge.arrow.size=0.3,
    edge.width=E(g)$weight * 5,
    main="Gene Regulatory Network")
dev.off()

# 使用ggraph更美观
ggraph(g, layout="fr") +
    geom_edge_link(aes(alpha=weight), arrow=arrow(length=unit(2,"mm"))) +
    geom_node_point(aes(color=type, size=size)) +
    geom_node_text(aes(label=name), size=2, repel=TRUE) +
    theme_void()

第六步：差异调控分析

代码示例：

# 比较不同条件间的TF活性变化
import scanpy as sc
from scipy.stats import ranksums

# 读取SCENIC结果
auc_scores = pd.read_csv("regulon_auc_matrix.csv", index_col=0)

# 比较两个条件
condition = adata.obs["condition"]
cond_A = auc_scores[condition == "disease"]
cond_B = auc_scores[condition == "healthy"]

# 对每个regulon做差异检验
diff_results = []
for tf in auc_scores.columns:
    stat, pval = ranksums(cond_A[tf], cond_B[tf])
    diff_results.append({
        "TF": tf,
        "mean_disease": cond_A[tf].mean(),
        "mean_healthy": cond_B[tf].mean(),
        "diff": cond_A[tf].mean() - cond_B[tf].mean(),
        "pvalue": pval
    })

diff_df = pd.DataFrame(diff_results)
from statsmodels.stats.multitest import multipletests
_, diff_df["padj"], _, _ = multipletests(diff_df["pvalue"], method="fdr_bh")
diff_df = diff_df.sort_values("padj")

# 上调的TF regulons
upregulated_tfs = diff_df[(diff_df["padj"] < 0.01) & (diff_df["diff"] > 0)]
print("疾病中上调的TF:", upregulated_tfs["TF"].tolist()[:10])

实战命令

# === SCENIC/pySCENIC完整流程 ===

# 准备数据库（一次性下载）
wget https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc_v10_clust/gene_based/hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather
wget https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc_v10_clust/gene_based/hg38_500bp_up_100bp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather
wget https://resources.aertslab.org/cistarget/motif2tf/motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl

# 运行pySCENIC
pyscenic grn input.loom allTFs_hg38.txt -o adj.tsv --method grnboost2 --num_workers 16
pyscenic ctx adj.tsv hg38_10kbp.feather hg38_500bp.feather --annotations_fname motifs.tbl --expression_mtx_fname input.loom --output regulons.csv --num_workers 16
pyscenic aucell input.loom regulons.csv -o scenic_output.loom --num_workers 16

面试常问点

Q1: SCENIC流程的三个步骤分别做什么？ A: 1）GRNBoost2/GENIE3：用随机森林推断TF-基因的共表达关系，输出候选调控边；2）RcisTarget：用TF结合motif信息过滤候选边，只保留靶基因启动子区有对应TF motif富集的关系（去除间接关联），形成regulon；3）AUCell：用AUC评分量化每个细胞中每个regulon的活性（regulon基因集是否被该细胞高表达）。

Q2: GENIE3和GRNBoost2的区别？ A: GENIE3用Random Forest，GRNBoost2用Gradient Boosting Trees。GRNBoost2快10-50倍，在大规模单细胞数据上实用。两者的网络推断质量相当，但GRNBoost2更适合>5000细胞的数据。

Q3: 推断的GRN如何验证？ A: 1）ChIP-seq验证TF实际结合靶基因启动子；2）Perturb-seq/CRISPRi验证敲低TF后靶基因表达下降；3）已知调控关系数据库（RegNetwork、TRRUST）比对；4）motif验证（SCENIC内置）；5）跨数据集一致性。

Q4: 为什么单细胞数据特别适合GRN推断？ A: 1）细胞间的表达异质性提供了大量"自然扰动"数据点；2）不同细胞代表不同的调控状态，类似于时间序列采样；3）细胞数量远大于基因数，避免了bulk中的n<p问题；4）可以获得细胞类型特异的调控网络。

易错点

  1. 计算资源不足：SCENIC对大数据集(>50K cells)非常耗时和耗内存。建议子采样(5000-10000 cells per cell type)。

  2. TF列表不准确：使用了过时或不完整的TF列表会遗漏重要调控关系或引入假阳性。

  3. 数据库版本不匹配：ranking database必须与基因组版本和注释版本一致。

  4. 过度解读共表达为调控：即使有motif支持，共表达也可能是共同被第三个TF调控，而非直接调控关系。

  5. 未考虑调控方向：标准SCENIC不区分激活和抑制。需要额外分析（如表达反相关）判断正/负调控。

补充知识

GRN推断工具比较

工具方法输入特点
SCENIC/pySCENICGRNBoost2+cisTarget+AUCellscRNA-seq金标准，motif验证
CellOracleGRN+扰动模拟scRNA+ATAC整合可及性
SCENIC+SCENIC+ATACscRNA+scATAC增强子调控
Dictys动态GRN时间序列scRNA动态网络
GENIE3随机森林bulk/scRNA简单独立
GRNBoost2梯度提升树scRNASCENIC内核

📚 相关文章推荐

单细胞调控网络SCENIC
60% 相关
743. 单细胞基因调控网络SCENIC+
57% 相关
378_单细胞转录因子活性推断
44% 相关
SCENIC 转录因子网络 — 单细胞基因调控网络推断工具
37% 相关
单细胞ATAC-seq分析（scATAC-seq）
26% 相关