第四章：RNA velocity 需要什么样的单细胞数据？
做 RNA velocity 之前，得先搞清楚手里拿的是什么数据。这玩意儿不是随便一个 scRNA-seq 数据集就能跑的。我在项目中见过不少同学，拿了个 10X 的普通 3' 端数据就直接跑 velocyto，结果出来的箭头乱七八糟。为什么？因为数据本身就不对路。
说白了，RNA velocity 要的是能区分未剪接和已剪接的 mRNA 分子。标准的 scRNA-seq 文库通常只抓 poly-A 尾巴，测出来的都是成熟 mRNA，看不到前体。你得换个思路——要找那些能同时捕获内含子信息的建库方式。
4.1 核心数据类型总览
我先列个清单，你看完就明白自己需要什么了。
| 数据类型 | 建库方式 | 是否支持 velocity | 典型平台 |
|---|---|---|---|
| 标准 3' 端 scRNA-seq | poly-A 捕获 + 逆转录 | ❌ 不支持 | 10X v2/v3, Drop-seq |
| 全长 scRNA-seq | Smart-seq2/3 | ✅ 支持 | Smart-seq2, Smart-seq3 |
| 内参与剪接捕获 | velocyto 协议 / 10X 改进版 | ✅ 原生支持 | 10X + velocyto, 10X 5' + intron |
| sci-Plex / 组合索引 | 组合条形码 + 全长扩增 | ✅ 支持 | sci-RNA-seq3, SPLiT-seq |
| 单细胞核 RNA-seq | snRNA-seq，捕获核内未剪接 | ✅ 支持 | 10X snRNA-seq, DroNc-seq |
关键结论：只要你的实验方案能测到内含子区域的 reads，就能算 RNA velocity。做不到？那箭头就是瞎画的。
4.2 为什么标准 3' 端数据不行？
10X 的 v2 和 v3 试剂盒，用的是 poly-T 引物抓 poly-A 尾巴。mRNA 前体还没剪接好，内含子还在，poly-A 尾巴也就还没加工完整，所以大部分前体分子压根没被捕获。我做过一个对比实验：同一样本分别用 10X v3 和 Smart-seq3 测，结果 10X 那边内含子比例不到 2%，Smart-seq3 这边有 15-20%。这差距，基本上决定了你能不能跑 velocity。
我的小技巧：拿到一个新数据集，先跑个 samtools idxstats 看看哪个染色体区域的 reads 多。如果发现 intron 区域几乎为 0，趁早换数据源，不用浪费时间调参数。
4.3 做 velocity 的数据必须满足的条件
嗯，这里要做个 checklist，你对照着看。
- reads 能跨外显子-内含子边界——建库时不能只截取 poly-A 端，要测足够长的片段，最好覆盖到内含子区域。
- 有足够的未剪接 reads 信号——一般要求未剪接占比 > 5%，太低的话箭头信噪比太差。
- 有准确的基因注释 GTF——内含子/外显子边界定义错了，velocity 的计数矩阵就是错的。
- 细胞数别太少——我个人建议低于 500 个细胞就不太合适了，统计模型跑不稳。
⚠️ 曾经踩过的坑：我用过一个 10X 5' 端数据（基因表达+VDJ），以为能抓内含子，但实际上 5' 端试剂盒的逆转录引物也是 poly-dT，效果和 3' 端一样差。别被“5'”两个字骗了，要看建库原理。
4.4 不同平台的数据格式差异
不管来自哪个平台，最终你都要转成同一个样子：一个跟基因×细胞对应的 spliced / unspliced / ambiguous 计数矩阵。velocyto 和 scVelo 都认这个格式。但不同平台产出的原始数据不一样。
| 平台 | 原始格式 | 需要什么预处理 |
|---|---|---|
| Smart-seq2 | FASTQ 或 BAM（每个细胞一个文件） | 需要自己从 BAM 提取跨剪接信息，或者用 alevin-fry |
| 10X + velocyto 协议 | FASTQ + 10X 白名单 | 用 velocyto run10x 命令直接出 loom 文件 |
| 10X 普通数据（用其它工具） | FASTQ + 10X 白名单 | 需要先做伪比对（kallisto bustools）再转 loom |
| snRNA-seq（Drop-seq 类型） | FASTQ + 条形码处理 | 同样用 velocyto 的 run_dropstar 或自行构建计数矩阵 |
你看，不同平台对技术要求不一样。我个人习惯先在笔记本上跑个 1000 个细胞的小样本，验证一下数据中有没有内含子信号，再上集群跑全量。别一上来就全量跑，万一数据不行，浪费几天时间。
4.5 怎么看数据是否合适？——快速检查命令
拿到 BAM 文件后，我会先跑一行命令看看内含子比例：
# 用 samtools 统计 intron 区域的 reads 数
# 假设你的基因注释文件为 genes.gtf，先提取 intron 区间
# 然后用 bedtools 或 IO 脚本统计
samtools view -b data.bam chr1:10000-20000 | samtools flagstat -更直接一点，用 velocyto run 之后生成的 .loom 文件里含三个矩阵（spliced / unspliced / ambiguous），你可以用 Python 看一眼比例：
import loompy
import numpy as np
loom_path = "your_data.loom"
with loompy.connect(loom_path) as ds:
spliced = ds[:, :][0] # 注意索引方式
unspliced = ds[:, :][1]
ratio = unspliced.sum() / (spliced.sum() + 1e-6)
print(f"未剪接 / 已剪接 比例: {ratio:.3f}")如果 ratio 小于 0.03，趁早放弃这个数据集。
4.6 数据类型的框架图
下面这张图总结了整个数据类型的判断流程。你拿到一个新数据集，顺着箭头走一圈，就能确定能不能做 velocity。
4.7 补充数据：单细胞核 RNA-seq 的特殊性
刚才提到了 snRNA-seq。这个很多人不知道，其实核内 RNA 天然包含大量未剪接前体。我做神经退行性研究的时候，细胞不好解离，就用了单细胞核 RNA-seq。结果 velocity 跑出来效果出奇好——因为核内富含 pre-mRNA，未剪接比例轻轻松松上 20%。
推荐场景：如果你的样本解离后细胞状态差、或者想研究神经元/心肌细胞这类大细胞，不妨直接用 snRNA-seq。省钱省力，velocity 效果反而更好。
4.8 总结一下数据准备清单
我不喜欢记太多理论，直接看列表对照就行。
- ✅ 优先选：Smart-seq3、10X + velocyto 改造协议、snRNA-seq
- ⚠️ 可以试：Smart-seq2、sci-RNA-seq（只要全长）
- ❌ 避免用：标准 10X v2/v3 / Drop-seq / 5' 端基因表达数据
再啰嗦一句：检查 GTF 版本。人类用 GENCODE v32+，小鼠用 GENCODE M25+。版本不对，内含子坐标全错，算出来的矩阵全是垃圾。我见过有人拿 Ensembl 的旧版注释跑，一跑 velocities 箭头全指向同一个方向——笑死，全是系统误差。