企业官网建设流程全解析

近几年单细胞测序的分析流程已经相当成熟，从质控、降维聚类到细胞类型注释，几乎每个环节都有稳定的工具链。但真正做发育、肿瘤克隆演化、谱系追踪这类课题时，很多人会卡在同一个地方：聚类结果是一张静态的“细胞地图”，而生物学关心的是“细胞怎么从 A 走到 B”。地图画得再漂亮，路径推不出来，故事就讲不完整。最近有团队把强化学习引入单细胞演化分析，把“细胞状态转移”显式建模成一个决策过程，思路确实比单纯堆叠轨迹推断工具更有解释力。这篇文章不复述论文结论，而是把这类工作背后的方法学拆成一套可以自己跑起来的复刻流程，包含环境准备、预处理、演化信号融合、强化学习模块设计、完整代码和排错清单。适合已经做过基础单细胞分析、想进一步做轨迹与克隆演化的读者，也适合想把强化学习落到真实科研场景的算法同学。

1. 从高影响因子工作看方法学价值：为什么要关注“强化学习 + 单细胞演化”

1.1 单细胞演化分析到底在解决什么问题

单细胞转录组测序（scRNA-seq）一次实验可以拿到成千上万个细胞的表达谱，每个细胞相当于一个“快照”。问题是这些快照是同时刻采集的，没有时间标签，我们看到的只是细胞群体在某一时刻的分布状态。

演化分析要做的，就是从这些静态快照里反推出动态过程：一个干细胞如何一步步分化成成熟细胞，一个肿瘤细胞克隆如何积累变异并扩张，一条发育轨迹上哪些细胞处于中间过渡态。通俗地说，聚类告诉我们“这个群体里有几种细胞”，演化分析告诉我们“这几种细胞之间的先后顺序和转换关系”。

常见的下游应用场景包括：造血系统的谱系分化研究、肿瘤耐药克隆的出现顺序、类器官发育成熟度评估、以及药物处理后细胞状态的整体偏移方向。这些问题的共同点是，结论必须依赖“方向性”，而方向性恰恰是普通聚类给不出来的。

1.2 传统生信流程的三处瓶颈

第一处瓶颈是方向性信号弱。常规流程基于表达相似度构建近邻图，相似度是无向的，A 和 B 近就互为邻居，但分不清谁在前谁在后。RNA 速率（RNA velocity）通过引入未剪接与已剪接 reads 的比例来提供方向，是目前最常用的补救方案，但它对测序深度、比对质量和基因覆盖度都比较敏感，很多数据集跑出来的速率场噪声很大。

第二处瓶颈是多信号难以统一。实际课题里往往同时有表达相似度、RNA 速率方向、拷贝数变异（CNV）推断出的克隆结构、以及已知的标记基因先验。这些信号量纲不同、可靠性不同，简单加权或者分步筛选都容易顾此失彼。

第三处瓶颈是路径搜索缺乏全局视角。多数轨迹推断方法本质上是局部贪心：从某个根节点出发，每次选最相似的方向往下走。局部最优不等于全局最优，遇到分支点、遇到细胞密度不均的区域，很容易走偏，甚至出现多条互相矛盾的轨迹。

1.3 把强化学习引入进来的直观理由

强化学习（Reinforcement Learning，RL）研究的是智能体如何在环境中通过试错来最大化累积回报。它有四个核心要素：状态、动作、奖励、策略。把这套框架套到细胞演化上，对应关系其实很自然：

• 状态：当前细胞所处的表达状态（通常用降维后的潜变量表示）。
• 动作：从当前细胞转移到某一个近邻细胞，或者选择“停止”。
• 奖励：这次转移在生物学上是否合理，比如方向是否与 RNA 速率一致、是否走向更高分化程度。
• 策略：学到的“在什么状态下该往哪走”的决策函数。

这样一来，“找一条合理的演化路径”就从局部贪心变成了一个带累积回报的序列决策问题。策略网络在整个细胞图上做全局优化，天然能处理分支和噪声。更重要的是，奖励函数是可解释的——你能明确说出模型为什么选这条路，这对科研场景非常关键。

1.4 本文的复刻范围与读者定位

需要说明的是，顶刊论文里的完整方法体系通常包含大量领域特定的建模细节和专门的验证实验，本文不可能也不应该声称逐行复刻。本文聚焦的是可迁移的方法学骨架：怎样把单细胞数据整理成 RL 能用的环境，怎样设计奖励函数把生物学先验编码进去，怎样训练、评估和解释策略输出。

读者定位上，如果你已经跑通过 scanpy 的标准流程、看得懂 AnnData 的结构，那读起来会比较顺畅；如果你还没接触过单细胞分析，建议先把第 4 节的基础流水线跑通再往下看。强化学习部分我只用最基础的 DQN 实现，不涉及复杂算法，有 PyTorch 基础即可。

2. 核心概念拆解：状态、动作、奖励与细胞命运

2.1 单细胞数据的基本形态

在 Python 生态里，单细胞数据的标准容器是 AnnData。可以把它理解成一张稀疏矩阵加两份注释：X是细胞 × 基因的表达矩阵（行是细胞，列是基因），obs存细胞的元信息（所属样本、聚类标签、QC 指标），var存基因的元信息。

这个结构看着简单，但有几个细节新手容易踩坑。第一，X通常是稀疏矩阵（scipy.sparse），直接转成 numpy 数组在几万个细胞时会吃掉大量内存。第二，obs_names必须唯一，合并多个样本时最容易出现重名。第三，不同来源的数据在拼接前必须先统一基因命名（比如MT-和mt-前缀的差异），否则会出现大量重复基因。

2.2 轨迹推断、RNA 速率与克隆演化三者的关系

这三个概念经常被混着用，其实层次不同。

轨迹推断（trajectory inference）是最宽泛的概念，目标是把细胞排成一条或多条有序的路径，输出通常是伪时间（pseudotime）。它只要求“排序”，不要求物理意义上的方向。

RNA 速率是一种提供方向证据的手段。它利用剪接动力学：基因转录后先产生未剪接的前体 mRNA，再剪接成成熟 mRNA，通过两者的比例可以推断该基因此刻是“上调中”还是“下调中”。把所有基因的速率拼起来，就得到一个高维方向向量，投影到低维空间后形成速率场。

克隆演化关注的是细胞之间的谱系关系，通常依赖 DNA 层面的变异（如 CNV、线粒体变异、CRISPR 条形码）。在 scRNA-seq 里，最常用的是从表达数据反推 CNV（拷贝数变异）来近似克隆结构。

三者的关系可以这样理解：克隆演化提供“骨架”（谁和谁同源），RNA 速率提供“箭头”（谁往谁那边走），轨迹推断负责把骨架和箭头组装成完整路径。RL 模块的作用，是在组装阶段做全局优化。

2.3 马尔可夫决策过程与细胞状态转移的类比

把细胞转移建模成马尔可夫决策过程（MDP），需要定义一个五元组：状态空间、动作空间、转移概率、奖励函数、折扣因子。

这里最容易引起争议的是“马尔可夫性”假设——当前状态包含了预测未来所需的全部信息。显然，一个细胞的表达谱并不能完全决定它的命运，微环境、表观遗传状态、相邻细胞 signaling 都会影响。所以这是一个建模假设，不是生物学事实。它的价值在于：把复杂的连续生物过程离散化成可计算的形式，让全局优化成为可能。理解了这一点，后面解释结果时就不会过度宣称。

折扣因子 γ 的作用是控制“看多远”。γ 接近 1 表示更看重长期回报，路径会走得更长；γ 偏小则偏向短期，容易在局部停住。在细胞演化场景里，通常设置在 0.9 到 0.99 之间比较合理。

2.4 奖励函数设计：把生物学先验写成可计算的分数

奖励函数是整个方法里最需要动脑子的部分，也是最能体现领域知识的地方。我在复刻时参考的思路是把奖励拆成三项：

第一项是方向一致性，计算当前细胞的速率向量与“指向目标细胞的位移向量”之间的余弦相似度。取值在 -1 到 1 之间，越接近 1 说明这次转移越符合速率场指示的方向。这一项直接把 RNA 速率信号编码进来了。

第二项是分化潜能变化。可以用基因集打分（比如干性相关基因集的平均表达）或者基于图结构的中心性指标来近似。走向潜能更低的细胞通常意味着分化推进，给正奖励；反向走则给负奖励。

第三项是回访惩罚。细胞图里存在大量局部环，如果不加约束，智能体很容易在两个细胞之间来回横跳刷分。对已访问节点施加负奖励，能有效抑制这种行为。

三项加权求和时，权重的选择需要根据数据情况调。经验做法是先让方向一致性占主导（比如权重 1.0），潜能变化作为辅助（0.3 到 0.5），回访惩罚保持在能压住死循环的量级（0.2 到 0.3）即可。

3. 环境准备与版本说明

3.1 硬件与操作系统

单细胞分析对内存的需求比 CPU 更敏感。一个 2 万细胞 × 3 万基因的稀疏矩阵，压缩后大约几百 MB 到 1 GB，但在做近邻图构建、速率计算时会临时膨胀数倍。建议至少 32 GB 内存，细胞数超过 5 万时建议 64 GB 起。

RL 训练部分反而是轻量级的，因为状态维度只有几十维（潜变量维度），网络很小，有 GPU 会快一些但不是必须的。CPU 训练一个几万步的 DQN 通常在几分钟到十几分钟量级。

操作系统上，Linux（Ubuntu 20.04 及以上）和 macOS 都可以，Windows 建议用 WSL2，主要是部分生信工具在原生 Windows 上的编译问题较多。

3.2 Python 环境与核心依赖

本文示例基于 Python 3.9 及以上版本，用 conda 管理环境。具体版本号这里不写死，因为 scanpy、scvelo 这类库迭代较快，不同版本间接口有变动，建议按下面的方式安装后在pip list里记录实际版本，方便复现。

# 创建并激活环境 conda create -n scrl python=3.10 -y conda activate scrl # 基础单细胞分析栈 pip install scanpy anndata pip install scvelo pip install leidenalg igraph # 深度学习 pip install torch --index-url https://download.pytorch.org/whl/cpu # 可视化与统计 pip install matplotlib seaborn scipy scikit-learn pandas

如果你要用 CNV 推断，还需要额外装 infercnvpy；如果需要做细胞通讯分析，可以加装其他工具包。本文主线不依赖这些可选组件，先保证主干跑通。

这里提醒一点：torch的安装命令要根据你的硬件调整，CUDA 版本对应的命令请参考官方安装页，不要直接照抄网上博客里的旧命令，很容易装到不匹配的版本。

3.3 数据与目录结构

示例数据用一个公开的 PBMC（外周血单个核细胞）数据集即可，10x Genomics 官方提供的 filtered 矩阵是常见选择，也可以换成自己项目的 h5ad 文件。目录结构建议这样组织：

sc-evolution-rl/ ├── data/ │ ├── raw/ # 原始 10x 矩阵目录 │ └── processed/ # 中间 h5ad 文件 ├── scripts/ │ ├── 01_preprocess.py │ ├── 02_velocity.py │ └── 03_rl_evolution.py ├── outputs/ │ ├── figures/ │ └── metrics/ └── requirements.txt

这种分层的好处是每一步的产物都落盘，中间出问题时不用从头跑。生信流程动辄几十分钟，能断点续跑非常重要。

3.4 版本兼容性提醒

有几个版本相关的坑值得提前说：

一是scanpy从 1.9 到 1.10 之间，部分函数的位置和默认参数有调整，比如高变基因筛选的默认方法。二是scvelo对scanpy版本比较敏感，装完建议先跑一遍官方教程的最小示例验证。三是anndata的AnnData对象在不同版本对obsm里稀疏矩阵的支持不一致，遇到保存后读取报错，先把稀疏转稠密再存。

稳妥做法是：环境搭好后，把pip freeze > requirements.txt存档，写文章或者投稿时附上，复现率会高很多。

4. 标准单细胞预处理流水线：复刻的地基

4.1 读取数据与质控

质控的目标是去掉明显损坏或者空液滴。三个核心指标：每个细胞检测到的基因数（n_genes_by_counts）、总 counts 数、线粒体基因占比（pct_counts_mt）。

线粒体占比高通常意味着细胞膜破损，胞质 RNA 流失而线粒体转录本残留。阈值设置要看组织类型，PBMC 这类细胞一般 10% 到 15% 之间，心肌、肝脏这类高代谢组织可以放宽。不要照抄别人的阈值，正确做法是先画分布图看拐点。

# 文件路径：scripts/01_preprocess.py import scanpy as sc sc.settings.verbosity = 2 sc.settings.n_jobs = 8 def load_and_qc(data_dir, min_genes=200, min_cells=3, max_mito=15.0): """读取 10x 矩阵并完成基础质控。""" adata = sc.read_10x_mtx( data_dir, var_names="gene_symbols", cache=True, ) adata.var_names_make_unique() # 细胞和基因层面的最低过滤 sc.pp.filter_cells(adata, min_genes=min_genes) sc.pp.filter_genes(adata, min_cells=min_cells) # 线粒体基因比例，注意人类是小写 mt-，小鼠是 Mt- adata.var["mt"] = adata.var_names.str.startswith(("MT-", "mt-")) sc.pp.calculate_qc_metrics( adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True, ) before = adata.n_obs adata = adata[adata.obs["pct_counts_mt"] < max_mito].copy() print(f"线粒体过滤: {before} -> {adata.n_obs} 个细胞") return adata

代码里var_names_make_unique()这一步不能省。10x 的基因注释里存在同名基因（比如一些 lncRNA 和蛋白编码基因共用一个 symbol），不去重会在后续按基因名索引时出问题。

4.2 归一化、高变基因与降维聚类

归一化解决的是测序深度差异。normalize_total把每个细胞的总 counts 缩放到同一量级（默认 1e4），log1p做对数变换压缩动态范围。这两步是标准操作，顺序不能反。

高变基因筛选的目的是降维去噪。3 万多个基因里，绝大多数在所有细胞里表达都很稳定，对区分细胞类型没有贡献，反而引入噪声。选 2000 到 3000 个高变基因是常规做法。

def preprocess(adata, n_top_genes=2000, n_pcs=30, n_neighbors=15, resolution=1.0): """归一化 + 降维 + 聚类。""" adata.layers["counts"] = adata.X.copy() # 保留原始 counts，速率分析要用 sc.pp.normalize_total(adata, target_sum=1e4) sc.pp.log1p(adata) adata.raw = adata # 保存对数化全基因矩阵 sc.pp.highly_variable_genes( adata, n_top_genes=n_top_genes, flavor="seurat", subset=True, ) sc.pp.scale(adata, max_value=10) sc.tl.pca(adata, n_comps=n_pcs, svd_solver="arpack") sc.pp.neighbors(adata, n_neighbors=n_neighbors, n_pcs=n_pcs) sc.tl.umap(adata) sc.tl.leiden(adata, resolution=resolution, key_added="leiden") return adata

这里有个细节值得展开：adata.layers["counts"]保存的是归一化之前的原始矩阵。后面的 RNA 速率分析必须用原始 counts，因为剪接动力学建模的是 reads 计数，不是对数化的相对表达。很多人跑到速率那一步发现报错说找不到spliced/unspliced，根源就是前面把原始矩阵覆盖掉了。

另外sc.pp.scale会让矩阵变稠密，细胞数超过 10 万时内存压力很大，可以考虑跳过 scale 直接用未缩放的 PCA，或者改用增量 PCA。

4.3 细胞类型注释与亚群筛选

聚类只是把细胞分成编号（0、1、2……），注释才是给它们赋予生物学身份。常用做法是看每个 cluster 的差异表达基因，再对照已知标记基因。

def annotate(adata, marker_dict): """基于标记基因打分做粗注释。""" for cell_type, markers in marker_dict.items(): valid = [g for g in markers if g in adata.raw.var_names] if len(valid) == 0: continue sc.tl.score_genes(adata, valid, score_name=f"score_{cell_type}") score_cols = [c for c in adata.obs.columns if c.startswith("score_")] adata.obs["cell_type"] = ( adata.obs[score_cols].idxmax(axis=1).str.replace("score_", "", regex=False) ) return adata

自动打分只适合做粗分类，精细亚群还是需要人工核对标记基因的表达模式。这一步不要偷懒，注释错了后面所有演化结论都会跟着错。做演化分析时，通常还要筛掉红细胞、血小板、双细胞（doublet）这些与目标过程无关的群体，否则它们会干扰图结构。

4.4 潜变量表示：为什么要单独做一步

PCA 是线性的，对非线性结构的表达能力有限。做演化分析时，更推荐用变分自编码器（如 scVI）学习一个低维潜变量表示，理由是：它显式建模了 counts 的负二项分布，能更好地处理 dropout（ dropout 是单细胞数据里大量零值的来源），而且潜空间通常更连续，对轨迹推断更友好。

如果你不想引入额外依赖，PCA 前 30 维也够用，只是后续速率场会更嘈杂。这是一个精度和复杂度的权衡，不是必须项。

5. 演化信号的三种来源与融合

5.1 RNA 速率：方向性最强的信号

scVelo 是目前最主流的速率分析工具，它需要每个基因的 spliced（已剪接）和 unspliced（未剪接）计数。这两个矩阵需要在比对阶段用 velocyto 或 STARsolo 单独产出，普通的 CellRanger 输出里没有。

如果你手上只有常规的 filtered 矩阵，是跑不了标准速率的。这一点必须在项目开始前确认清楚，别做到一半才发现数据不支持。

拿到 loom 文件后，流程大致是：

# 文件路径：scripts/02_velocity.py import scvelo as scv def run_velocity(adata, n_top_genes=2000, n_neighbors=30): """计算 RNA 速率并构建速率图。""" scv.pp.filter_and_normalize( adata, min_shared_counts=20, n_top_genes=n_top_genes, ) scv.pp.moments(adata, n_pcs=30, n_neighbors=n_neighbors) scv.tl.velocity(adata, mode="stochastic") scv.tl.velocity_graph(adata) # 潜时间，后续用来和 RL 策略的伪时间做一致性检验 scv.tl.latent_time(adata) return adata

mode="stochastic"是最基础的稳态模型假设，对数据量要求低；mode="dynamical"精度更高但要拟合每个基因的动力学参数，耗时长很多。建议先用 stochastic 跑通看效果，满意了再上 dynamical。

5.2 拷贝数变异推断：克隆结构的骨架

CNV 推断的思路是：用一组正常细胞（比如免疫细胞）作为参考，计算每个细胞在染色体上各基因区域的相对表达，滑动窗口平滑后得到一条“拷贝数谱”。同一克隆的细胞谱型相似，不同克隆差异明显。

这个信号的价值在于它相对稳定，不受转录瞬时波动影响，适合用来定义克隆边界。局限是分辨率有限，只有大片段的增减才能检测到，而且对参考细胞的选择很敏感——参考选错了，整个谱型都会偏。

在 RL 框架里，CNV 更适合作为约束而不是奖励项：如果两个细胞属于不同克隆，转移概率应该被压低。实现上可以在动作选择时做掩码（mask），把跨克隆的邻居屏蔽掉。

5.3 系统发育树与克隆指派

有了 CNV 谱或者线粒体变异，就可以构建细胞间的系统发育树。常见做法是对变异位点做距离矩阵，然后用邻接法（Neighbor-Joining）或者最大简约法建树。

树的用途是给 RL 提供一个“祖先-后代”的先验：树根附近的细胞更可能是祖先状态，树梢是终末状态。这个先验可以直接转成潜能打分的一部分，也可以用来约束起始状态的采样范围。

需要提醒的是，单细胞层面的树推断噪声很大，尤其当细胞数上万时，很多分支的支持度很低。实践中通常先聚成克隆（cluster），再在克隆层面建树，稳定性会好很多。

5.4 多信号一致性打分

三种信号各自有可靠性问题，所以融合时不能简单相加。一个实用做法是计算一致性：对每一对相邻细胞，看速率方向和 CNV 相似度是否指向同一个方向，一致则置信度高，矛盾则降低权重。

这个一致性打分可以直接作为 RL 的环境质量指标：整体一致性高的数据集，训练出来的策略更可信；一致性低的数据集，就应该在结论里明确说明不确定性。

6. 强化学习模块的设计与实现

6.1 把细胞图变成环境

环境类的职责是封装状态转移逻辑。核心数据包括：每个细胞的潜变量表示、速率向量（投影到潜空间）、潜能打分、以及 k 近邻索引。

这里有个实现细节要注意：速率向量原本在基因空间（几万维），必须投影到和潜变量相同的空间才能算余弦相似度。投影方式有两种，一是用 PCA 的载荷矩阵做线性投影，二是训练一个小的回归模型。前者简单可靠，推荐先用它。

# 文件路径：scripts/03_rl_evolution.py import random import numpy as np import torch import torch.nn as nn import torch.optim as optim class CellGraphEnv: """把单细胞近邻图封装成 MDP 环境。""" def __init__(self, z, vel, potency, knn_idx, max_steps=30): self.z = z # (N, d) 潜变量 self.vel = vel # (N, d) 投影后的速率 self.potency = potency # (N,) 分化潜能 0~1 self.knn = knn_idx # (N, k) 近邻索引 self.max_steps = max_steps self.n = z.shape[0] self.k = knn_idx.shape[1] self.cur = 0 self.steps = 0 self.visited = set() @property def obs_dim(self): return self.z.shape[1] * 2 + 1 def reset(self, start=None): self.cur = random.randrange(self.n) if start is None else int(start) self.steps = 0 self.visited = {self.cur} return self._obs() def _obs(self): c = self.cur nb = self.z[self.knn[c]].mean(axis=0) # 邻域聚合，提供局部上下文 return np.concatenate( [self.z[c], nb, [self.potency[c]]] ).astype(np.float32) def step(self, action): c = self.cur if action == self.k: # 动作空间最后一位是“停止” return self._obs(), 0.0, True, {"stop": True} nxt = int(self.knn[c, action]) direction = self.z[nxt] - self.z[c] denom = np.linalg.norm(self.vel[c]) * np.linalg.norm(direction) + 1e-8 cos = float(np.dot(self.vel[c], direction) / denom) dpot = float(self.potency[nxt] - self.potency[c]) reward = 1.0 * cos + 0.5 * dpot if nxt in self.visited: reward -= 0.3 # 回访惩罚 self.visited.add(nxt) self.cur = nxt self.steps += 1 done = self.steps >= self.max_steps return self._obs(), reward, done, {"cos": cos, "dpot": dpot}

动作空间大小是k + 1，多出来的那一位是停止动作。停止动作给 0 奖励，这样智能体在走到终末状态时会自己学会停下，而不是被迫走满步数。

6.2 状态表示与网络结构

状态向量由三部分拼成：当前细胞的潜变量、邻域平均潜变量、当前潜能值。加入邻域信息是为了让智能体感知局部密度和方向，纯靠单点坐标在噪声大的区域容易迷路。

网络用两层全连接就够，输入维度几十维，隐藏层 128 或 256。输出维度等于动作空间大小，每个动作对应一个 Q 值。

class QNet(nn.Module): def __init__(self, obs_dim, n_actions, hidden=128): super().__init__() self.net = nn.Sequential( nn.Linear(obs_dim, hidden), nn.ReLU(), nn.Linear(hidden, hidden), nn.ReLU(), nn.Linear(hidden, n_actions), ) def forward(self, x): return self.net(x)

不要把这个网络做深。状态维度低、样本相关性高，深层网络很容易过拟合到某几条特定路径上，泛化反而变差。

6.3 动作空间与转移约束

动作空间的大小直接由 k 近邻的 k 决定。k 太小，智能体可选的转移方向少，容易卡住；k 太大，动作空间膨胀，探索效率下降。经验值 15 到 30 之间比较合适，和数据规模有关。

转移约束可以通过掩码实现。比如跨克隆转移不允许，那就在计算 Q 值时把对应动作的 Q 值设为极小值（如 -1e9），softmax 或者 argmax 之后自然不会选到。这种方式比在环境里硬性拒绝更平滑，也不会破坏梯度。

6.4 奖励函数与训练循环

训练用标准 DQN：经验回放 + 目标网络。回放池容量设 5 万到 10 万，批大小 128，目标网络每隔几百步同步一次。

from collections import deque def train_dqn(env, total_steps=20000, batch_size=128, gamma=0.95, lr=1e-3, eps_start=1.0, eps_end=0.05, eps_decay=8000, buffer_size=50000, target_sync=500): obs_dim = env.obs_dim n_actions = env.k + 1 policy_net = QNet(obs_dim, n_actions) target_net = QNet(obs_dim, n_actions) target_net.load_state_dict(policy_net.state_dict()) target_net.eval() optimizer = optim.Adam(policy_net.parameters(), lr=lr) buffer = deque(maxlen=buffer_size) obs = env.reset() losses = [] for step in range(total_steps): eps = eps_end + (eps_start - eps_end) * np.exp(-step / eps_decay) if random.random() < eps: action = random.randrange(n_actions) else: with torch.no_grad(): q = policy_net(torch.from_numpy(obs).unsqueeze(0)) action = int(q.argmax(dim=1).item()) next_obs, reward, done, _ = env.step(action) buffer.append((obs, action, reward, next_obs, float(done))) if done: obs = env.reset() else: obs = next_obs if len(buffer) < batch_size: continue batch = random.sample(buffer, batch_size) b_obs, b_act, b_rew, b_next, b_done = zip(*batch) b_obs = torch.from_numpy(np.stack(b_obs)) b_next = torch.from_numpy(np.stack(b_next)) b_act = torch.tensor(b_act, dtype=torch.long) b_rew = torch.tensor(b_rew, dtype=torch.float32) b_done = torch.tensor(b_done, dtype=torch.float32) q_cur = policy_net(b_obs).gather(1, b_act.unsqueeze(1)).squeeze(1) with torch.no_grad(): q_next = target_net(b_next).max(dim=1).values q_target = b_rew + gamma * q_next * (1.0 - b_done) loss = nn.functional.smooth_l1_loss(q_cur, q_target) optimizer.zero_grad() loss.backward() nn.utils.clip_grad_norm_(policy_net.parameters(), 5.0) optimizer.step() losses.append(loss.item()) if step % target_sync == 0: target_net.load_state_dict(policy_net.state_dict()) return policy_net, losses

几个参数值得说明。gamma=0.95是折中值，路径平均长度在 20 步左右。eps_decay=8000控制探索衰减速度，衰减太快会早熟收敛到局部路径。梯度裁剪clip_grad_norm_在奖励尺度不稳定时很有必要，防止某次异常大的奖励把网络带偏。损失用smooth_l1_loss而不是 MSE，对离群奖励更鲁棒。

6.5 策略输出如何回读成生物学结论

训练完成后，策略网络不能直接给出“演化路径”这个结论，需要做后处理。常用做法是：

第一步，从每个细胞出发，用贪心策略（eps=0）rollout 若干步，记录访问序列。第二步，把所有序列汇总成一个有向转移计数矩阵。第三步，对矩阵做归一化，得到细胞间的转移概率。第四步，在转移概率图上跑最长路径或者找吸收态，得到主路径和分支点。

只有这样转成图结构，才能和 scVelo 的 latent time、已知的标记基因表达曲线做对比验证。直接看策略网络的权重是没有生物学意义的。

7. 完整实战：从原始矩阵到演化路径图

7.1 项目结构

按第 3 节给出的目录组织，三个脚本串行执行：01_preprocess.py产出质控后的 h5ad，02_velocity.py补充速率信息，03_rl_evolution.py训练策略并输出路径。

# 文件路径：scripts/run_all.py import subprocess import sys steps = [ "scripts/01_preprocess.py", "scripts/02_velocity.py", "scripts/03_rl_evolution.py", ] for s in steps: print(f"=== 执行 {s} ===") r = subprocess.run([sys.executable, s]) if r.returncode != 0: print(f"{s} 执行失败，终止") break

把执行顺序固化到脚本里，是为了避免手动操作时漏步或者顺序错乱。生信复现出问题，十有八九是步骤顺序或者中间产物版本对不上。

7.2 依赖安装与数据就位

依赖安装见 3.2 节。数据就位后先做一次快速检查，确认矩阵能读、细胞数在合理范围：

# 文件路径：scripts/check_data.py import scanpy as sc adata = sc.read_10x_mtx("data/raw/filtered_feature_bc_matrix", var_names="gene_symbols") print("细胞数:", adata.n_obs) print("基因数:", adata.n_vars) print("矩阵类型:", type(adata.X)) print("前 5 个细胞:", list(adata.obs_names[:5]))

预期输出大致是细胞数几千到上万、基因数两万上下、矩阵类型为scipy.sparse的某个子类。如果细胞数是 0，基本可以确定路径写错了（10x 目录下应该有matrix.mtx.gz、barcodes.tsv.gz、features.tsv.gz三个文件）。

7.3 预处理脚本串起来

把 4.1 和 4.2 的函数组装成完整脚本，最后落盘：

# 文件路径：scripts/01_preprocess.py（完整版下半部分） if __name__ == "__main__": import os marker_dict = { "T": ["CD3D", "CD3E", "IL7R"], "B": ["MS4A1", "CD79A", "CD79B"], "NK": ["GNLY", "NKG7", "KLRD1"], "Monocyte": ["LYZ", "CD14", "FCGR3A"], } adata = load_and_qc("data/raw/filtered_feature_bc_matrix") adata = preprocess(adata) adata = annotate(adata, marker_dict) os.makedirs("data/processed", exist_ok=True) adata.write_h5ad("data/processed/adata_qc.h5ad") print("已保存:", adata.shape) print(adata.obs["cell_type"].value_counts())

跑完之后看cell_type的分布，PBMC 数据里 T 细胞通常占比最高，其次是单核细胞和 B 细胞。如果某个类型的数量异常少（比如 T 细胞只有几十个），说明注释或者质控参数有问题，需要回头调整。

7.4 速率与转移图构建

速率计算完成后，需要把速率向量投影到潜空间，并构建 k 近邻索引。这里用 sklearn 的 NearestNeighbors：

# 文件路径：scripts/02_velocity.py（衔接部分） import numpy as np from sklearn.neighbors import NearestNeighbors def build_graph_inputs(adata, k=20): """准备 RL 环境所需的 z、vel、potency、knn。""" z = adata.obsm["X_pca"][:, :30] # 速率投影：用 PCA 的旋转矩阵把基因空间速率映射到潜空间 vel_gene = adata.layers["velocity"] if hasattr(vel_gene, "toarray"): vel_gene = vel_gene.toarray() rot = adata.varm["PCs"][:, :30] vel_latent = vel_gene @ rot # 简单干性打分：高变基因的平均标准化表达，反向归一化到 0~1 raw = np.asarray(adata.raw[:, adata.var_names].X.mean(axis=1)).ravel() potency = (raw - raw.min()) / (raw.max() - raw.min() + 1e-8) nn = NearestNeighbors(n_neighbors=k + 1, metric="euclidean").fit(z) _, idx = nn.kneighbors(z) knn = idx[:, 1:] # 去掉自身 return z.astype(np.float32), vel_latent.astype(np.float32), \ potency.astype(np.float32), knn.astype(np.int64)

注意idx[:, 1:]去掉自身这一步，kneighbors返回的第一列是细胞自己。

7.5 训练与结果导出

# 文件路径：scripts/03_rl_evolution.py（主流程） def rollout_policy(policy_net, env, n_episodes=2000, max_steps=30): """用贪心策略做 rollout，统计细胞间转移频次。""" transfer = {} paths = [] for _ in range(n_episodes): s = random.randrange(env.n) obs = env.reset(start=s) seq = [s] for _ in range(max_steps): with torch.no_grad(): q = policy_net(torch.from_numpy(obs).unsqueeze(0)) a = int(q.argmax(dim=1).item()) prev = env.cur obs, _, done, _ = env.step(a) if done: break seq.append(env.cur) transfer[(prev, env.cur)] = transfer.get((prev, env.cur), 0) + 1 paths.append(seq) return transfer, paths if __name__ == "__main__": import anndata as ad adata = ad.read_h5ad("data/processed/adata_velocity.h5ad") z, vel, potency, knn = build_graph_inputs(adata, k=20) env = CellGraphEnv(z, vel, potency, knn, max_steps=30) policy_net, losses = train_dqn(env, total_steps=20000) torch.save(policy_net.state_dict(), "outputs/policy.pt") transfer, paths = rollout_policy(policy_net, env) print("采样的路径条数:", len(paths)) print("转移边数:", len(transfer))

7.6 结果验证

结果验证是整篇文章里最容易被忽略、但最不能省的一步。至少做两件事：

一是一致性检验。把 rollout 得到的伪时间（路径上的步序）和 scVelo 的latent_time做 Spearman 相关。相关系数在 0.5 以上说明两者大体一致；如果接近 0 甚至为负，说明 RL 学到的东西和速率场矛盾，需要回查奖励函数权重或者速率质量。

二是生物学合理性检验。挑几个已知的标记基因，看它们沿主路径的表达曲线是否单调。比如 B 细胞分化路径上MS4A1应该逐渐升高，如果曲线上下跳动，说明路径本身有问题。

from scipy.stats import spearmanr def validate(adata, paths, key="latent_time"): n = adata.n_obs step_first = np.full(n, np.nan) for p in paths: for i, c in enumerate(p): if np.isnan(step_first[c]): step_first[c] = i mask = ~np.isnan(step_first) rho, pval = spearmanr(step_first[mask], adata.obs[key].values[mask]) print(f"Spearman rho={rho:.3f}, p={pval:.3e}") return rho

覆盖到的细胞比例也要统计。如果只有一小部分细胞被 rollout 访问过，说明策略探索不足，可以加大训练步数或者提高探索率下限。

8. 常见问题与排查清单

8.1 高频问题对照表

8.2 训练不收敛的三种排查方向

奖励全为常数是最常见的原因。如果速率向量全是 0（投影矩阵维度对不上、或者基因顺序错位），余弦相似度就恒为 0，智能体学不到任何有效信号。排查方法是打印vel的范数分布，正常情况下应该有明显区分度。

第二种是奖励尺度差异过大。cos 在 -1 到 1，potency 差值可能只有 0.001 量级，两项相加时后一项完全被淹没。解决方法是先对 potency 做标准化，或者调大它的权重。

第三种是动作空间过大而样本不足。k 设为 50 而训练步数只有 5000，每个动作平均只被采样到几十次，Q 值估计方差极大。降低 k 或者加大训练量。

8.3 结果不可复现的处理

随机种子必须全程固定，包括random、numpy、torch三处。另外 DataLoader 如果用多进程也要设置 worker 的种子。

import random import numpy as np import torch def set_seed(seed=42): random.seed(seed) np.random.seed(seed) torch.manual_seed(seed) torch.cuda.manual_seed_all(seed) torch.backends.cudnn.deterministic = True

即便如此，GPU 上的浮点运算仍可能有微小差异，在长训练中会被放大。如果对复现要求严格，建议在 CPU 上做最终验证。

9. 最佳实践与工程建议

9.1 数据与实验管理

单细胞项目的中间文件又多又大，必须有清晰的命名规则。我的习惯是{样本}_{处理步骤}_{日期}.h5ad，比如pbmc_qc_20240115.h5ad。每次参数调整都新建文件而不是覆盖，虽然占空间，但回溯时非常省事。

实验记录建议用一个 CSV 或者表格维护：跑了哪个脚本、用了什么参数、产出了什么指标。RL 训练涉及大量超参数组合，靠脑子记一定会乱。

9.2 评估指标要提前定好

不要等训练完了再想“怎么判断结果好不好”。开跑之前就应该确定：主指标是 Spearman 相关系数还是路径覆盖率？可接受的阈值是多少？有没有独立的验证数据集？

如果条件允许，最好留一个独立样本做外部验证。在 PBMC 上训出来的策略换到另一个人的 PBMC 数据上还能用，才说明学到的是普适规律而不是数据集的噪声。

9.3 计算资源与性能

RL 训练本身不慢，瓶颈通常在环境交互。每次step都要做一次 numpy 索引和简单的向量运算，Python 层面开销不小。优化方向：把knn和z预加载到共享内存，或者批量并行多个环境（vectorized env）。

生信部分的时间分布通常是：质控几分钟、PCA 和聚类几分钟、速率计算十几分钟到几小时（取决于基因数和模式）、CNV 推断看数据规模。建议先用 5% 的细胞做小规模试跑，把流程调通再上全量。

9.4 结论边界与生物学解释

这一点必须强调。RL 学到的策略是在一堆建模假设下得到的计算结论，不是直接的生物学事实。至少有三层不确定性：马尔可夫假设是否成立、奖励函数是否完整刻画了分化驱动力、数据本身是否覆盖了关键过渡状态。

写论文或者报告时，建议把结论表述成“在该模型假设下，策略倾向于认为 A 细胞通过 B 中间态过渡到 C”，而不是“A 一定分化成 C”。同时附上敏感性分析：换一组奖励权重，主路径是否稳定？如果不稳定，说明结论不稳健。

另外，涉及临床样本或人源数据时，务必确认数据使用符合伦理审批和授权范围，不要跨项目混用未授权的数据。

10. 学习路线与后续扩展方向

如果这套流程你跑通了，接下来可以往几个方向深入。

第一是强化学习算法本身的替换。DQN 适合离散动作空间，如果你的场景需要连续的动作（比如在潜空间里预测一个连续的方向向量），可以换成 PPO 或者 SAC。代码改动主要在策略网络和训练循环，环境接口基本不用动。

第二是引入更多生物学先验。比如细胞通讯分析得到的配体-受体信号，可以作为奖励项的一部分，让策略倾向于沿着 signaling 活跃的方向转移。再比如把已知的转录因子调控网络编码成约束，抑制不合理的跳跃。

第三是不确定性量化。现在的实现只给出一条最优路径，没有置信度。可以用集成方法（训多个策略取平均）或者贝叶斯神经网络，给出每条转移边的置信区间。这对科研结论的可信度提升很大。

第四是和批次校正结合。多供体、多时间点的数据必然有批次效应，如果批次没有处理好，策略可能学到的是批次差异而不是真实的演化方向。建议在预处理阶段就做好整合（harmony、scVI 等），再进入演化分析。

最后给一个务实的学习顺序：先把第 4 节的预处理跑熟，理解每个参数在干什么；然后花时间啃 scVelo 的文档和论文，把速率这件事搞清楚；最后再回头看 RL 部分。反过来，先学 RL 再学生信，很容易写出“能跑但没意义”的代码。

如果你在复刻过程中遇到策略不收敛、速率场异常、或者和 latent time 对不上的情况，可以先按第 8 节的表格逐条排查，大部分问题都出在数据准备阶段而不是模型本身。