跳转至

151_单细胞RNA速率

一句话概述

RNA velocity(RNA速率)通过区分单细胞中未剪接(unspliced)和已剪接(spliced)mRNA的比例，推断每个细胞的基因表达变化方向和速度，从而预测细胞的未来状态和分化轨迹。

核心知识点总览

知识点说明
RNA velocity原理利用unspliced/spliced mRNA比值推断表达变化方向
稳态模型假设转录-剪接达到平衡，偏离平衡指示变化方向
动态模型允许参数随时间变化，更准确但计算更复杂
velocyto从BAM中提取spliced/unspliced计数的原始工具
scVeloPython包，实现稳态和动态RNA velocity
UniTVelo统一时间RNA velocity，解决多基因时间不一致问题
STARsoloSTAR的单细胞模式，可直接输出velocity矩阵
速率向量可视化在UMAP/tSNE上显示箭头指示细胞转变方向
潜在时间(latent time)从velocity推断的细胞伪时间序列
局限性稳态假设不总成立、测序深度影响、多基因不一致性

各步骤详解

第一步：RNA velocity的生物学原理

白话解释： 基因表达是一个动态过程：DNA先转录出前体mRNA(含intron，即"unspliced")，然后经过剪接变成成熟mRNA(去除intron，即"spliced")。如果一个基因正在被"开启"(转录上调)，它的unspliced mRNA会先增多(因为新转录的还来不及剪接)。反之，如果被"关闭"(转录下调)，unspliced会先减少。通过看unspliced和spliced的比值，就能推断这个基因正在被上调还是下调——就像看河流的水位变化来判断上游是在下雨还是天晴。

技术原理： 基因表达动力学的微分方程：

du/dt = α(t) - β·u       (unspliced mRNA变化)
ds/dt = β·u - γ·s         (spliced mRNA变化)
- α：转录速率 - β：剪接速率 - γ：降解速率 - u：unspliced mRNA量 - s：spliced mRNA量

稳态模型 (scVelo stochastic)： - 当 du/dt = 0 且 ds/dt = 0 时，达到稳态：u = α/β, s = α/γ - 因此稳态时 u = (γ/β)·s (线性关系) - velocity = ds/dt = β·u - γ·s - 正velocity → 基因正在上调；负velocity → 基因正在下调

代码概念：

import numpy as np
import matplotlib.pyplot as plt

# 模拟一个基因的表达动力学
alpha = 5.0   # 转录速率
beta = 0.5    # 剪接速率
gamma = 0.3   # 降解速率

t = np.linspace(0, 30, 1000)
u = np.zeros_like(t)
s = np.zeros_like(t)

dt = t[1] - t[0]
for i in range(1, len(t)):
    # 在t=10时转录上调(alpha翻倍)
    a = alpha * 2 if 10 <= t[i] <= 20 else alpha
    u[i] = u[i-1] + (a - beta * u[i-1]) * dt
    s[i] = s[i-1] + (beta * u[i-1] - gamma * s[i-1]) * dt

# u/s 的相位图
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.plot(t, u, label='unspliced')
plt.plot(t, s, label='spliced')
plt.xlabel('Time'); plt.ylabel('Expression'); plt.legend()
plt.subplot(1, 2, 2)
plt.scatter(s, u, c=t, cmap='viridis', s=1)
plt.xlabel('spliced'); plt.ylabel('unspliced')
plt.colorbar(label='Time')
plt.tight_layout()
plt.savefig('velocity_concept.png', dpi=150)

第二步：从比对文件提取spliced/unspliced计数

白话解释： 普通scRNA-seq分析只看"总表达量"。做RNA velocity需要额外区分哪些reads来自外显子(spliced)、哪些来自内含子(unspliced)、哪些横跨外显子-内含子边界(ambiguous)。

方法一：velocyto

# 安装velocyto
pip install velocyto

# 10x Genomics数据
velocyto run10x \
  /path/to/cellranger/output/ \
  genes.gtf \
  -m repeat_msk.gtf \
  --samtools-threads 16

# 输出: velocyto/sample.loom (包含spliced/unspliced/ambiguous矩阵)

# Smart-seq2数据
velocyto run_smartseq2 \
  -o output_dir/ \
  -e sample_name \
  bam_files/*.bam \
  genes.gtf

# 参数说明:
# -m repeat_msk.gtf : 重复序列mask文件(推荐, 减少假unspliced)
# 下载: UCSC Table Browser → RepeatMasker

方法二：STARsolo (推荐)

# STARsolo可以直接在比对时生成velocity矩阵
STAR --genomeDir star_index/ \
  --readFilesIn read2.fastq.gz read1.fastq.gz \
  --readFilesCommand zcat \
  --soloType CB_UMI_Simple \
  --soloCBwhitelist 3M-february-2018.txt \
  --soloFeatures Gene Velocyto \
  --outSAMtype BAM SortedByCoordinate \
  --outSAMattributes NH HI nM AS CR UR CB UB GX GN sS sQ sM \
  --runThreadN 16

# 输出目录结构:
# Solo.out/Velocyto/
#   ├── filtered/
#   │   ├── spliced.mtx
#   │   ├── unspliced.mtx
#   │   ├── ambiguous.mtx
#   │   ├── barcodes.tsv
#   │   └── features.tsv

方法三：alevin-fry (Salmon生态)

# 使用扩展的参考(包含intron序列)
# 先生成splici参考
pyroe make-splici \
  genome.fa \
  genes.gtf \
  90 \
  splici_ref/

# Salmon alevin量化
salmon alevin -l ISR \
  --index splici_index \
  -1 read1.fastq.gz \
  -2 read2.fastq.gz \
  --chromiumV3 \
  -o alevin_output/ \
  -p 16 \
  --sketch

# 使用pyroe加载

第三步：scVelo 标准分析流程

白话解释： scVelo是做RNA velocity的主力工具。它接收spliced/unspliced计数矩阵，拟合动力学模型，计算每个基因在每个细胞中的velocity，然后在降维空间上展示细胞的"运动方向"。

代码示例：

import scvelo as scv
import scanpy as sc
import anndata as ad

# === 1. 数据加载 ===
# 从loom文件加载
adata_loom = scv.read("velocyto_output.loom", cache=True)

# 或从STARsolo加载
from scipy.io import mmread
import pandas as pd

spliced = mmread("Solo.out/Velocyto/filtered/spliced.mtx").T.tocsr()
unspliced = mmread("Solo.out/Velocyto/filtered/unspliced.mtx").T.tocsr()
barcodes = pd.read_csv("Solo.out/Velocyto/filtered/barcodes.tsv", header=None)[0]
features = pd.read_csv("Solo.out/Velocyto/filtered/features.tsv", header=None, sep="\t")[1]

adata = ad.AnnData(X=spliced)
adata.layers["spliced"] = spliced
adata.layers["unspliced"] = unspliced
adata.obs_names = barcodes.values
adata.var_names = features.values

# === 2. 合并已有的scanpy分析结果(UMAP, clusters) ===
adata_processed = sc.read("processed_scanpy.h5ad")
# 确保barcode匹配
common_cells = adata.obs_names.intersection(adata_processed.obs_names)
adata = adata[common_cells]
adata_processed = adata_processed[common_cells]
adata.obs = adata_processed.obs.copy()
adata.obsm = adata_processed.obsm.copy()

# === 3. 预处理 ===
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

# === 4. 稳态模型 (stochastic) ===
scv.tl.velocity(adata, mode="stochastic")
scv.tl.velocity_graph(adata)

# 可视化
scv.pl.velocity_embedding_stream(adata, basis="umap", color="celltype",
                                   save="velocity_stream.png", dpi=150)
scv.pl.velocity_embedding_grid(adata, basis="umap", color="celltype",
                                 save="velocity_grid.png", dpi=150)

# === 5. 动态模型 (推荐, 更准确) ===
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.scatter(adata, color="latent_time", color_map="gnuplot",
               save="latent_time.png")

# === 6. 关键基因分析 ===
# 查看velocity最高/最低的基因
scv.tl.rank_velocity_genes(adata, groupby="celltype")
top_genes = adata.uns["rank_velocity_genes"]["names"]

# 单基因phase portrait
scv.pl.velocity(adata, var_names=["SOX2", "GATA1", "KLF1"],
                save="gene_velocity.png")

# 速度置信度
scv.tl.velocity_confidence(adata)
scv.pl.scatter(adata, color="velocity_confidence",
               save="confidence.png")

第四步：UniTVelo — 统一时间模型

白话解释： scVelo的一个问题是每个基因独立估计参数，不同基因推断出的"时间"可能不一致(基因A说细胞在早期，基因B说在晚期)。UniTVelo通过一个统一的时间框架解决这个问题——让所有基因共享同一套时间轴。

代码示例：

import unitvelo as utv
import scvelo as scv

# 1. 加载数据 (同scVelo的adata)
adata = scv.read("velocity_data.h5ad")

# 2. 配置UniTVelo
velo_config = utv.config.Configuration()
velo_config.R2_ADJUST = True       # R2校正
velo_config.IROOT = None           # 自动检测root cell
velo_config.FIT_OPTION = "1"       # 拟合选项: "1"=基本模型, "2"=统一时间模型
velo_config.GPU = 0                # GPU编号 (-1表示CPU)

# 3. 运行UniTVelo
adata = utv.run_model(
    adata,
    label="celltype",
    config_file=velo_config
)

# 4. 可视化
scv.pl.velocity_embedding_stream(adata, basis="umap", color="celltype",
                                   save="unitvelo_stream.png")
# UniTVelo输出的latent_time
scv.pl.scatter(adata, color="latent_time",
               save="unitvelo_latent_time.png")

第五步：velocity结果解读与下游分析

白话解释： velocity箭头只是起点——更重要的是据此做下游分析：推断分化路径、识别驱动基因、估计细胞间转变概率。

代码示例：

import scvelo as scv
import cellrank as cr

# === 1. 速度伪时间 ===
scv.tl.latent_time(adata)
scv.pl.scatter(adata, color="latent_time", color_map="gnuplot", size=80)

# === 2. 驱动基因排序 ===
scv.tl.rank_velocity_genes(adata, groupby="celltype", min_corr=0.3)
df = scv.DataFrame(adata.uns["rank_velocity_genes"]["names"])
print(df.head(10))  # 每种细胞类型的top velocity基因

# === 3. PAGA velocity (轨迹图) ===
# 结合PAGA拓扑图和velocity方向
adata.uns["neighbors"]["distances"] = adata.obsp["distances"]
adata.uns["neighbors"]["connectivities"] = adata.obsp["connectivities"]
scv.tl.paga(adata, groups="celltype")
scv.pl.paga(adata, basis="umap", size=50, alpha=0.1,
            min_edge_width=2, node_size_scale=1.5,
            save="paga_velocity.png")

# === 4. CellRank (基于velocity的命运概率) ===
# CellRank将velocity信息转化为细胞命运概率
vk = cr.kernels.VelocityKernel(adata)
vk.compute_transition_matrix()

ck = cr.kernels.ConnectivityKernel(adata)
ck.compute_transition_matrix()

# 组合kernel
combined_kernel = 0.8 * vk + 0.2 * ck

# 估计终末态
estimator = cr.estimators.GPCCA(combined_kernel)
estimator.compute_macrostates(n_states=5, cluster_key="celltype")
estimator.set_terminal_states_from_macrostates()
estimator.compute_absorption_probabilities()

# 可视化命运概率
estimator.plot_absorption_probabilities(same_plot=False)

实战命令

# === 环境搭建 ===
conda create -n velocity python=3.9
conda activate velocity
pip install scvelo scanpy velocyto cellrank unitvelo

# === velocyto提取 (10x数据) ===
velocyto run10x \
  /path/to/cellranger/outs/ \
  genes.gtf \
  -m repeat_msk.gtf \
  --samtools-threads 16

# === STARsolo一步生成velocity矩阵 ===
STAR --genomeDir star_index/ \
  --readFilesIn R2.fq.gz R1.fq.gz \
  --readFilesCommand zcat \
  --soloType CB_UMI_Simple \
  --soloCBwhitelist whitelist.txt \
  --soloFeatures Gene Velocyto \
  --runThreadN 16

# === scVelo完整分析脚本 ===
python -c "
import scvelo as scv
adata = scv.read('data.loom')
scv.pp.filter_and_normalize(adata)
scv.pp.moments(adata)
scv.tl.recover_dynamics(adata, n_jobs=8)
scv.tl.velocity(adata, mode='dynamical')
scv.tl.velocity_graph(adata)
scv.tl.latent_time(adata)
adata.write('velocity_result.h5ad')
scv.pl.velocity_embedding_stream(adata, basis='umap', save='velocity.png')
"

面试常问点

Q1：RNA velocity的基本原理是什么？

A： RNA velocity利用单细胞中unspliced(含内含子)和spliced(成熟)mRNA的比例来推断基因表达变化的方向。当基因转录上调时，新产生的pre-mRNA(unspliced)会先增多；下调时unspliced先减少。通过拟合unspliced和spliced的动力学关系(du/dt和ds/dt微分方程)，计算出每个基因在每个细胞中的velocity(ds/dt)，正值表示上调，负值表示下调。

Q2：稳态模型和动态模型有什么区别？

A： 稳态模型假设每个基因的转录(α)、剪接(β)、降解(γ)速率是恒定的，通过全局拟合u-s的线性关系来估计γ/β比值。动态模型(scVelo dynamical)允许这些参数随时间变化，能够恢复完整的基因表达动力学过程(诱导期→稳态→抑制期)。动态模型更准确但计算量大(约30分钟到数小时)，稳态模型几分钟就能完成。

Q3：velocity分析的主要局限性有哪些？

A： (1) 稳态假设在快速变化的系统中不成立；(2) 10x 3'末端测序对unspliced mRNA的捕获效率低；(3) 不同基因的"时间"可能不一致(UniTVelo试图解决)；(4) 低表达基因的velocity估计噪声大；(5) velocity反映的是"瞬时方向"而非长期轨迹；(6) 循环过程(如细胞周期)可能干扰分化trajectory的解读。

Q4：如何获取unspliced/spliced计数矩阵？

A： 三种主要方法：(1) velocyto从CellRanger的BAM文件中根据reads与外显子/内含子的overlap进行分类；(2) STARsolo在比对时直接通过--soloFeatures Velocyto生成spliced/unspliced/ambiguous矩阵；(3) alevin-fry使用splici参考(外显子+内含子序列)进行量化。STARsolo推荐用于新项目，因为可以一步完成。

Q5：velocity的结果如何与CellRank结合？

A： CellRank将velocity信息转化为细胞间的转移概率矩阵(Markov chain)，然后通过吸收概率分析推断每个细胞到达各终末状态的概率。这比单纯看velocity箭头更robust，因为它结合了velocity方向和细胞间的转录组相似性。CellRank可以识别初始态、终末态、命运决定点，以及驱动命运决定的基因。

易错点

1. 未正确合并velocity数据与已有分析

错误： loom文件中的barcode与h5ad中的barcode格式不匹配(如缺少"-1"后缀)。 正确做法： 检查两边barcode格式，统一后再合并。velocyto产生的barcode可能带有样本前缀(如"sample:")，需要清理。

2. 对所有数据类型盲目使用velocity

错误： 对CITE-seq蛋白数据或ATAC-seq数据做RNA velocity。 正确做法： RNA velocity仅适用于mRNA数据。ATAC-seq有ChromatinVelocity等替代方法。

3. 忽略unspliced reads的比例

错误： 不检查unspliced reads的比例就跑分析。 正确做法： 检查adata.layers["unspliced"].sum() / adata.layers["spliced"].sum()。比值太低(<0.1)说明unspliced捕获不足，velocity结果不可靠。10x 3'数据通常为0.1-0.3，full-length数据为0.3-0.5。

4. 过度解读velocity箭头方向

错误： 将UMAP上的箭头方向等同于真实的分化方向。 正确做法： velocity箭头受UMAP降维影响。建议结合PAGA拓扑图、latent time、已知生物学marker来验证。箭头显示的是局部趋势而非绝对方向。

5. 对低质量细胞不做过滤

错误： 包含doublet、dying cells等低质量细胞做velocity分析。 正确做法： 先做好严格的QC过滤(mito%、gene count、doublet检测)，再进行velocity分析。低质量细胞的unspliced/spliced比例异常，会严重干扰结果。

补充知识

velocity工具生态演化

时期工具特点
2018velocyto开创性工作，稳态模型
2020scVelo动态模型，Python生态集成
2022UniTVelo统一时间框架
2022MultiVelo多组学velocity(RNA+ATAC)
2023VeloVAE变分自编码器建模velocity
2024Dynamo更丰富的动力学建模

多组学velocity

  • MultiVelo: 同时利用RNA和ATAC-seq信息推断velocity
  • chromatin velocity: 从scATAC-seq推断染色质开放性变化方向
  • protein velocity: 利用CITE-seq蛋白数据的velocity

计算资源

  • scVelo稳态模型: ~5分钟 (5万细胞)
  • scVelo动态模型: ~30分钟-2小时 (5万细胞)
  • UniTVelo: ~1小时 (5万细胞, GPU加速)
  • 内存需求: 约8-16GB (5万细胞)

📚 相关文章推荐

742. 单细胞RNA速度分析scVelo
50% 相关
382_单细胞时间序列分析
43% 相关
单细胞轨迹分析（Trajectory Analysis）
42% 相关
scVelo RNA 速率进阶 — 动力学模型推断细胞分化方向
37% 相关
Velocyto RNA 速率 — 计算 RNA 剪接速率推断细胞状态走向
26% 相关