742. 单细胞RNA速度分析scVelo¶

一句话概述：通过比较未剪接RNA（刚转录的）和已剪接RNA（成熟的）的比例，推断每个细胞的发育方向——就像看每个细胞正在"往哪个方向走"，从而预测细胞命运。

核心知识点速查表¶

概念	白话解释	关键工具
RNA速度(RNA velocity)	细胞基因表达的变化速率和方向	scVelo, velocyto
未剪接RNA(Unspliced)	刚从DNA转录出来，还带着内含子	代表"将来的表达"
已剪接RNA(Spliced)	去掉内含子的成熟mRNA	代表"当前的表达"
潜在时间(Latent time)	每个细胞在分化轨迹上的"内部时钟"	scVelo动态模型
稳态模型	假设大多数基因处于稳态	velocyto原始方法
动态模型	允许基因在不同细胞中处于不同阶段	scVelo EM模型

一、原理（白话版）¶

1.1 RNA速度的核心思想¶

基因表达的过程：
DNA → [转录] → pre-mRNA(未剪接) → [剪接] → mRNA(已剪接) → [降解] → 消失

关键洞察：
- 如果一个基因正在"上调"：未剪接RNA多，已剪接RNA还没跟上
  → u/s比值高 → 这个基因在"加速"

- 如果一个基因正在"下调"：未剪接RNA少，已剪接RNA还没降解完
  → u/s比值低 → 这个基因在"减速"

综合所有基因的速度 → 推断细胞整体的发展方向

1.2 稳态模型 vs 动态模型¶

特性	稳态模型(Velocyto)	动态模型(scVelo)
假设	大多数细胞处于稳态	每个细胞可以在任何阶段
方程	ds/dt = β·u - γ·s	完整ODE方程组
参数估计	线性回归	EM算法
适用范围	简单分化过程	复杂分化/多分支
计算速度	快	慢

二、分析流程¶

2.1 准备数据（从BAM文件提取剪接信息）¶

# ===== 用velocyto提取剪接/未剪接计数 =====
# velocyto从10x Chromium的BAM文件中提取spliced和unspliced RNA

# 安装
pip install velocyto  # 安装velocyto

# 运行（10x Chromium数据）
velocyto run10x \
  /path/to/cellranger/output/ \  # Cell Ranger输出目录
  /path/to/genes.gtf \  # 基因组注释GTF文件
  -m /path/to/repeat_msk.gtf \  # 重复区域mask文件（可选但推荐）
  --samtools-threads 8  # samtools线程数

# 输出：sample.loom文件
# 包含三层矩阵：spliced, unspliced, ambiguous

2.2 scVelo标准分析¶

# ===== scVelo完整分析流程 =====
import scvelo as scv  # 导入scVelo
import scanpy as sc   # 导入Scanpy
import matplotlib.pyplot as plt  # 导入绑图

# 设置scVelo参数
scv.settings.verbosity = 3  # 显示详细信息
scv.settings.presenter_view = True  # 演示模式（大字体）
scv.set_figure_params("scvelo")  # 使用scVelo默认绘图参数

# ===== 第一步：加载数据 =====
# 方法一：直接加载loom文件
adata = scv.read("sample.loom", cache=True)  # 读取loom文件

# 方法二：将loom与已有的h5ad合并
adata_rna = sc.read_h5ad("processed.h5ad")  # 已处理的scRNA-seq数据
ldata = scv.read("sample.loom")  # velocyto输出的loom文件

# 合并（确保细胞barcode匹配）
adata = scv.utils.merge(adata_rna, ldata)  # 合并spliced/unspliced矩阵到adata

# ===== 第二步：预处理 =====
# 过滤和归一化
scv.pp.filter_and_normalize(
    adata,
    min_shared_counts=20,  # 至少20个细胞共享的基因
    n_top_genes=2000,      # 使用top 2000高变基因
    retain_genes=None      # 不额外保留基因
)

# 计算矩（moments）— 基于邻近细胞平滑表达值
scv.pp.moments(
    adata,
    n_pcs=30,   # 使用30个主成分
    n_neighbors=30  # 每个细胞的30个最近邻
)

# ===== 第三步：估计RNA速度 =====
# 方法一：稳态模型（快速但简单）
scv.tl.velocity(adata, mode="stochastic")  # 随机模型（稳态的改进版）

# 方法二：动态模型（推荐，更准确）
scv.tl.recover_dynamics(adata, n_jobs=8)  # EM算法估计反应速率参数
scv.tl.velocity(adata, mode="dynamical")  # 动态模型

# 构建速度图
scv.tl.velocity_graph(adata)  # 计算细胞间的速度转换概率

# ===== 第四步：可视化 =====
# 在UMAP上显示速度流向
scv.pl.velocity_embedding_stream(
    adata,
    basis="umap",  # 使用UMAP降维
    color="cell_type",  # 按细胞类型着色
    legend_loc="right",
    title="RNA Velocity Stream",
    figsize=(10, 8),
    save="velocity_stream.png",
    dpi=300
)

# 速度箭头图（每个细胞一个箭头）
scv.pl.velocity_embedding(
    adata,
    basis="umap",
    color="cell_type",
    arrow_length=3,    # 箭头长度
    arrow_size=2,      # 箭头大小
    title="RNA Velocity Arrows",
    save="velocity_arrows.png"
)

2.3 潜在时间分析¶

# ===== 潜在时间（Latent Time）=====
# 只有动态模型才能计算

# 计算潜在时间
scv.tl.latent_time(adata)  # 估计每个细胞的内部时钟

# 可视化潜在时间
scv.pl.scatter(
    adata,
    color="latent_time",  # 按潜在时间着色
    color_map="gnuplot",  # 紫→黄色阶
    size=50,
    title="Latent Time",
    save="latent_time.png"
)

# 查看驱动基因（速度贡献最大的基因）
# 这些基因最能体现分化方向
top_genes = adata.var["fit_likelihood"].sort_values(ascending=False).index[:10]
scv.pl.scatter(
    adata,
    basis=top_genes[:4],  # 显示前4个驱动基因的phase plot
    ncols=2,
    frameon=False,
    save="top_driver_genes.png"
)

# 驱动基因热图
scv.pl.heatmap(
    adata,
    var_names=top_genes,  # 显示top基因
    sortby="latent_time",  # 按潜在时间排序
    col_color="cell_type",  # 列按细胞类型着色
    n_convolve=100,  # 平滑窗口
    save="driver_genes_heatmap.png"
)

2.4 速度置信度评估¶

# ===== 评估速度结果的可靠性 =====
# 速度置信度
scv.tl.velocity_confidence(adata)  # 计算速度置信度

# 可视化置信度
scv.pl.scatter(
    adata,
    c=["velocity_length", "velocity_confidence"],  # 速度大小和置信度
    cmap="coolwarm",
    perc=[5, 95],
    save="velocity_confidence.png"
)

# PAGA速度（结合PAGA轨迹推断）
scv.tl.paga(
    adata,
    groups="cell_type",  # 按细胞类型分组
    minimum_spanning_tree=False
)

scv.pl.paga(
    adata,
    basis="umap",
    size=50,
    alpha=0.1,
    min_edge_width=2,
    node_size_scale=1.5,
    title="PAGA Velocity",
    save="paga_velocity.png"
)

三、常见报错与解决¶

报错信息	原因	解决方案
`No spliced/unspliced layers`	loom文件未正确合并	检查`adata.layers`中是否有spliced和unspliced
`Too few cells`	过滤太严格	降低min_shared_counts阈值
`Negative velocities everywhere`	数据质量差或预处理问题	检查unspliced比例(应占5-25%)
`recover_dynamics: ConvergenceWarning`	EM未收敛	增加max_iter或换用stochastic模式
`velocity_graph: memory error`	细胞数太多	下采样到50k细胞以内
`Inconsistent velocity direction`	可能是生物学真实或方法局限	对比不同方法(scVelo vs CellRank)

四、面试高频问题¶

Q1: RNA速度的基本假设是什么？¶

A: ①转录→剪接→降解的单向过程；②剪接和降解速率在短时间内恒定；③未剪接RNA可以从BAM文件中准确识别；④邻近细胞代表连续的分化状态。

Q2: scVelo的动态模型相比稳态模型好在哪？¶

A: 稳态模型假设所有细胞都在稳态平衡点附近，只能推断平衡时的速度。动态模型允许每个细胞处于诱导或抑制的不同阶段，能捕捉完整的转录动力学，还能估计"潜在时间"。

Q3: RNA速度有什么局限性？¶

A: ①依赖未剪接RNA的准确计量（10x数据中unspliced比例本身就低）；②不能处理非稳态的剪接速率（如温度变化时）；③不能区分"转录增加"和"降解减慢"；④在复杂分支中可能给出错误方向。

五、速查表¶

# ===== scVelo分析速查 =====

# 安装
pip install scvelo

# 提取spliced/unspliced（从Cell Ranger BAM）
velocyto run10x cellranger_output/ genes.gtf

# 标准流程
import scvelo as scv
adata = scv.read("sample.loom")
scv.pp.filter_and_normalize(adata, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

# 动态模型（推荐）
scv.tl.recover_dynamics(adata, n_jobs=8)
scv.tl.velocity(adata, mode="dynamical")
scv.tl.velocity_graph(adata)
scv.tl.latent_time(adata)

# 可视化
scv.pl.velocity_embedding_stream(adata, basis="umap")
scv.pl.scatter(adata, color="latent_time")

# 关键检查
# unspliced比例：5-25%正常
# velocity_confidence > 0.5较可靠