先说我为什么会被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来模拟敲低之后沿着网络传播信号预测每个细胞转录状态会发生多大程度的偏移。我整理了一个对比表方便你按需选型。工具输入数据显式使用ATAC峰扰动模拟输出侧重WGCNARNA表达否否表达共模块SCENICRNA表达否否regulon模块Pando多组学是否带回归系数的GRNFigR多组学是弱/间接DORCTF活性CellOracle多组学是有base GRN扰动流场需要说清楚的是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扫描的底层库它依赖rdkitrdkit的二进制包用conda装远比pip装省事。我的推荐是混合安装先创建干净的conda环境再用pip补上CellOracle本身和Scanpy。conda create -n celloracle python3.9 -y conda activate celloracle conda install -c conda-forge gimmemotifs rdkit pip install celloracle0.10.13 scanpy python-igraph pip install numpy1.27 pandas2.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峰计数矩阵peak×cell的AnnDatavar_names为chr_起始_终止提供可及性峰和调控元件坐标RNA表达矩阵gene×cell的AnnData基因名用gene symbol提供基因表达信息和细胞状态坐标基因组fasta与peak坐标对应的参考基因组motif扫描时提取峰序列基因注释GTF匹配fasta版本确定TSS和peak-基因对应关系TF motif数据库HOCOMOCO v11等人源/鼠源库判断峰上有哪些TF的结合位点很多人卡在第一步ATAC的AnnData必须用raw count不要先normalize或者做过log变换。你关心的是这个峰在哪些细胞里有可及性信号后续回归步骤自己会处理尺度问题提前标准化反而可能把计数分布压坏。peak的ID格式也注意一下。CellOracle内部会把var_names解析成染色体坐标常见格式是chr1_1000_2000这种下划线分隔。如果你的数据是chr1:1000-2000这种冒号格式先用正则转一下再导入不然到了扫描阶段会报坐标解析错误。这个转换本身不难但很容易被忽略。2.3 基因组注释和motif数据库的准备两个容易搞混的坑先说基因组版本。你的peak来自哪个参考基因组就只能用哪个版本的fasta和gtfhg38的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_comps50) sc.pp.neighbors(adata_rna) sc.tl.umap(adata_rna) sc.tl.leiden(adata_rna, resolution0.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(adataadata_atac, genomehg38) # 检查peak注释是否被正确解析 oracle.get_peak_info() # 导入TF motif信息这里用HOCOMOCO人源v11 TFinfo co.data.load_hocomoco_data(speciesHuman) oracle.import_TF_data(TFinfoTFinfo) # 把scanpy算出的细胞状态同步给Oracle # 不同小版本这一步的函数名可能有差异以官方notebook为准 oracle.get_scanpy_information(adataadata_atac, rna_adataadata_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_conditionperturb_condition, n_simulation50, overwriteTrue, ) # 在UMAP上画流场图 oracle.plot_perturb_umap( plot_graphFalse, plot_arrowTrue, arrow_density0.5, saveoutput/IKZF1_knockdown_stream.png, )n_simulation50表示重复50次模拟然后汇总delta值太小时噪声很大太大则运行时间成倍增加。arrow_density0.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_THREADS85.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敲掉之后细胞群会不会发生状态转换这类问题这个教程的代码应该能直接帮你把流程启动起来。跑通之后我强烈建议再沿着官方文档把网络子图可视化、多种扰动对比这些进阶功能补上那些对文章的提升非常明显。