scVelo — RNA 速率分析
概述
scVelo 是領先的 Python 套件，用於單細胞 RNA-seq 資料的 RNA 速率分析。它透過對 mRNA 剪接動力學建模——利用未剪切（pre-mRNA）與已剪切（成熟 mRNA）豐度的比值，判斷每個細胞中基因是否正在上調或下調——來推估細胞狀態轉換。這使得我們能夠重建發育軌跡並識別細胞命運決策，而無需時間序列資料。
安裝： uv pip install scvelo
主要資源：
- 文件：https://scvelo.readthedocs.io/
- GitHub：https://github.com/theislab/scvelo
- 論文：Bergen et al. (2020) Nature Biotechnology. PMID: 32747759
何時使用此技能
在以下情況使用 scVelo：
- 從快照資料進行軌跡推斷：判定細胞分化方向
- 細胞命運預測：識別前驅細胞及其下游 fate
- 驅動基因識別：找出最能解釋觀察軌跡的基因動態
- 發育生物學：建模造血作用、神經發育、上皮-間質轉換
- 潛在時間估計：根據剪接動力學對細胞進行偽時間排序
- Scanpy 補充：為 UMAP 嵌入增添方向性資訊
先決條件
scVelo 需要同時提供 未剪切 和 已剪切 RNA 的計數矩陣。這些可由以下工具產生：
- STARsolo 或 kallisto|bustools 搭配
lamanno模式 - velocyto CLI：
velocyto run10x/velocyto run - alevin-fry / simpleaf 輸出已剪切/未剪切結果
資料儲存於包含 layers["spliced"] 和 layers["unspliced"] 的 AnnData 物件中。
標準 RNA 速率工作流程
1. 設定與資料載入
import scvelo as scv
import scanpy as sc
import numpy as np
import matplotlib.pyplot as plt
# 配置設定
scv.settings.verbosity = 3 # 顯示計算步驟
scv.settings.presenter_view = True
scv.settings.set_figure_params('scvelo')
# 載入資料（包含已剪切/未剪切層的 AnnData）
# 選項 A：從 loom 檔案載入（velocyto 輸出）
adata = scv.read("cellranger_output.loom", cache=True)
# 選項 B：將 velocyto loom 與 Scanpy 處理過的 AnnData 合併
adata_processed = sc.read_h5ad("processed.h5ad") # 包含 UMAP、分群
adata_velocity = scv.read("velocyto.loom")
adata = scv.utils.merge(adata_processed, adata_velocity)
# 驗證層
print(adata)
# obs × var: N × G
# layers: 'spliced', 'unspliced'（必要）
# obsm['X_umap']（用於視覺化）
2. 預處理
# 篩選與歸一化。從 scVelo 0.3 起，filter_and_normalize() 僅篩選
# 基因並對每個細胞進行歸一化——不再接受 n_top_genes 參數，也不再進行
# 對數轉換，因此對數轉換步驟和 HVG 選擇需由 Scanpy 完成。
scv.pp.filter_and_normalize(
adata,
min_shared_counts=20 # 已剪切+未剪切中的最小計數
)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000, subset=True)
# 計算一階和二階動差（均值與變異數）
# 必須先計算 knn_connectivities
sc.pp.neighbors(adata, n_neighbors=30, n_pcs=30)
scv.pp.moments(
adata,
n_pcs=30,
n_neighbors=30
)
3. 速度估計 — 隨機模型
隨機模型速度快，適合探索性分析：
# 隨機速度（較快，準確度較低）
scv.tl.velocity(adata, mode='stochastic')
scv.tl.velocity_graph(adata)
# 視覺化
scv.pl.velocity_embedding_stream(
adata,
basis='umap',
color='leiden',
title="RNA Velocity (Stochastic)"
)
4. 速度估計 — 動態模型（推薦）
動態模型擬合完整的剪接動力學，準確度更高：
# 恢復動力學（計算密集；10K 細胞約需 10-30 分鐘）
scv.tl.recover_dynamics(adata, n_jobs=4)
# 從動態模型計算速度
scv.tl.velocity(adata, mode='dynamical')
scv.tl.velocity_graph(adata)
5. 潛在時間
動態模型可計算共享的潛在時間（偽時間）：
# 計算潛在時間
scv.tl.latent_time(adata)
# 在 UMAP 上視覺化潛在時間
scv.pl.scatter(
adata,
color='latent_time',
color_map='gnuplot',
size=80,
title='Latent time'
)
# 識別按潛在時間排序的頂尖基因
top_genes = adata.var['fit_likelihood'].sort_values(ascending=False).index[:300]
scv.pl.heatmap(
adata,
var_names=top_genes,
sortby='latent_time',
col_color='leiden',
n_convolve=100
)
6. 驅動基因分析
# 識別擬合度最高的速度基因
scv.tl.rank_velocity_genes(adata, groupby='leiden', min_corr=0.3)
df = scv.DataFrame(adata.uns['rank_velocity_genes']['names'])
print(df.head(10))
# 速度與一致性
scv.tl.velocity_confidence(adata)
scv.pl.scatter(
adata,
c=['velocity_length', 'velocity_confidence'],
cmap='coolwarm',
perc=[5, 95]
)
# 特定基因的相圖
scv.pl.velocity(adata, ['Cpe', 'Gnao1', 'Ins2'],
ncols=3, figsize=(16, 4))
7. 速度箭頭與偽時間
# UMAP 上的箭頭圖
scv.pl.velocity_embedding(
adata,
arrow_length=3,
arrow_size=2,
color='leiden',
basis='umap'
)
# 流線圖（更清晰的視覺化）
scv.pl.velocity_embedding_stream(
adata,
basis='umap',
color='leiden',
smooth=0.8,
min_mass=4
)
# 速度偽時間（潛在時間的替代方案）
scv.tl.velocity_pseudotime(adata)
scv.pl.scatter(adata, color='velocity_pseudotime', cmap='gnuplot')
8. PAGA 軌跡圖
# 結合速度資訊的 PAGA 圖
scv.tl.paga(adata, groups='leiden')
df = scv.get_df(adata, 'paga/transitions_confidence', precision=2).T
df.style.background_gradient(cmap='Blues').format('{:.2g}')
# 繪製含速度的 PAGA 圖
scv.pl.paga(
adata,
basis='umap',
size=50,
alpha=0.1,
min_edge_width=2,
node_size_scale=1.5
)
完整工作流程腳本
import scvelo as scv
import scanpy as sc
def run_rna_velocity(adata, n_top_genes=2000, mode='dynamical', n_jobs=4):
"""
完整 RNA 速率工作流程。
Args:
adata: 包含 'spliced' 和 'unspliced' 層、UMAP 在 obsm 中的 AnnData
n_top_genes: 用於速度的頂尖 HVG 數量
mode: 'stochastic'（快速）或 'dynamical'（準確）
n_jobs: 動態模型的並行作業數
Returns:
包含速度資訊的處理後 AnnData
"""
scv.settings.verbosity = 2
# 1. 預處理（scVelo 0.3 已將對數轉換和 HVG 從 filter_and_normalize 移除）
scv.pp.filter_and_normalize(adata, min_shared_counts=20)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=n_top_genes, subset=True)
if 'neighbors' not in adata.uns:
sc.pp.neighbors(adata, n_neighbors=30)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
# 2. 速度估計
if mode == 'dynamical':
scv.tl.recover_dynamics(adata, n_jobs=n_jobs)
scv.tl.velocity(adata, mode=mode)
scv.tl.velocity_graph(adata)
# 3. 下游分析
if mode == 'dynamical':
scv.tl.latent_time(adata)
scv.tl.rank_velocity_genes(adata, groupby='leiden', min_corr=0.3)
scv.tl.velocity_confidence(adata)
scv.tl.velocity_pseudotime(adata)
return adata
AnnData 中的主要輸出欄位
執行工作流程後，會新增以下欄位：
| 位置 | 鍵 | 說明 |
|---|---|---|
adata.layers | velocity | 每個基因每個細胞的 RNA 速率 |
adata.layers | fit_t | 每個基因每個細胞的擬合潛在時間 |
adata.obsm | velocity_umap | UMAP 上的 2D 速度向量 |
adata.obs | velocity_pseudotime | 從速度推導的偽時間 |
adata.obs | latent_time | 從動態模型推導的潛在時間 |
adata.obs | velocity_length | 每個細胞的速度大小 |
adata.obs | velocity_confidence | 每個細胞的信心分數 |
adata.var | fit_likelihood | 基因層級的模型擬合品質 |
adata.var | fit_alpha | 轉錄速率 |
adata.var | fit_beta | 剪接速率 |
adata.var | fit_gamma | 降解速率 |
adata.uns | velocity_graph | 細胞-細胞轉移機率矩陣 |
速度模型比較
| 模型 | 速度 | 準確度 | 使用時機 |
|---|---|---|---|
stochastic | 快 | 中等 | 探索性分析；大規模資料集 |
deterministic | 中等 | 中等 | 簡單線性動力學 |
dynamical | 慢 | 高 | 發表級分析；識別驅動基因 |
最佳實踐
- 先以隨機模式探索，最終分析再轉為動態模式
- 需要良好的未剪切讀取覆蓋率：短讀取（< 100 bp）可能無法涵蓋內含子
- 至少 2,000 個細胞：細胞數較少時 RNA 速率雜訊較大
- 速度應具有一致性：箭頭應遵循已知生物學；隨機性表示有問題
- k-NN 頻寬很重要：鄰居太少 → 速度雜訊大；太多 → 過度平滑
- 合理性檢查：根細胞（前驅細胞）的標記基因應具有高未剪切/已剪切比率
- 動態模型需要不同的動力學狀態：最適合明確的分化過程
疑難排解
| 問題 | 解決方案 |
|---|---|
| 缺少未剪切層 | 重新執行 velocyto 或使用 STARsolo 搭配 --soloFeatures Gene Velocyto |
| 速度基因極少 | 降低 min_shared_counts；檢查定序深度 |
| 箭頭看似隨機 | 嘗試不同的 n_neighbors 或速度模型 |
| 動態模型記憶體錯誤 | 設定 n_jobs=1；減少 n_top_genes |
| 到處都是負速度 | 檢查已剪切/未剪切層是否未交換 |
其他資源
- scVelo 文件：https://scvelo.readthedocs.io/
- 教學筆記本：https://scvelo.readthedocs.io/tutorials/
- GitHub：https://github.com/theislab/scvelo
- 論文：Bergen V et al. (2020) Nature Biotechnology. PMID: 32747759
- velocyto（預處理）：http://velocyto.org/
- CellRank（命運預測，擴展 scVelo）：https://cellrank.readthedocs.io/
- dynamo（代謝標記替代方案）：https://dynamo-release.readthedocs.io/