Velocyto/RNA 速度分析

2025年12月19日 · Tutorials

TL;DR: RNA velocity 通过比较新生（含内含子）与成熟（已剪接）mRNA来预测细胞状态转换。使用 Velocyto 从 BAM 文件计数内含子，然后使用 scVelo 进行分析和可视化。这对于轨迹推断和理解动态过程至关重要。

什么是 RNA Velocity？

RNA velocity 估计单个细胞中基因表达的速率和方向。

核心概念

细胞包含两种类型的 mRNA：

未剪接（含内含子）= 新转录的新生 RNA
已剪接（不含内含子）= 成熟的、可翻译的 RNA

关键见解： 如果一个基因的未剪接 mRNA 多于已剪接 mRNA，则它正在被激活。反之，则正在被关闭。

RNA Velocity 概念

为何重要

传统 scRNA-seq	RNA Velocity
静态快照	预测未来状态
"细胞 A 位于位置 X"	"细胞 A 正在向位置 Y 移动"
需要伪时间算法	动态过程的直接测量
分化路径不明确	解析分支

何时使用 RNA Velocity

✅ 理想应用场景：

分化轨迹（干细胞 → 成熟细胞）
细胞周期阶段转换
免疫激活反应
发育时间进程研究
肿瘤进展动态

⚠️ 局限性：

需要高质量的未剪接 reads（并非所有实验方案都能捕获内含子）
假设 RNA 转录/剪接/降解处于平衡状态
在稳态系统中效果较差
计算密集

工作流程概述

BAM 文件 → Velocyto（计数内含子）→ .loom 文件 → scVelo（分析）→ 速度向量

第一部分：Velocyto（计数内含子）

安装

# 创建 Python 3.8 环境
conda create -n velocyto python=3.8
conda activate velocyto

# 安装依赖
conda install numpy scipy cython numba matplotlib scikit-learn h5py click

# 通过 pip 安装
pip install pysam velocyto

测试安装：

velocyto --version

准备所需文件

1. 重复序列屏蔽 GTF

从 UCSC 基因组浏览器下载：

访问 UCSC Table Browser
选择：
- Clade: Mammal
- Genome: Human（或你的生物体）
- Assembly: GRCh38（与你的数据匹配）
- Group: Repeats
- Track: RepeatMasker
- Output format: GTF
点击 "get output"

保存为： GRCh38_rmsk.gtf

2. 基因注释 GTF

使用与你的 Cell Ranger/STAR 比对相同的 GTF：

# Cell Ranger 参考基因组
~/refdata-gex-GRCh38-2024-A/genes/genes.gtf

# 或从 Ensembl 下载
wget ftp://ftp.ensembl.org/pub/release-110/gtf/homo_sapiens/Homo_sapiens.GRCh38.110.gtf.gz
gunzip Homo_sapiens.GRCh38.110.gtf.gz

3. BAM 文件

来自 Cell Ranger 或 STAR 比对：

# Cell Ranger
cellranger_output/outs/possorted_genome_bam.bam

# STAR
STAR_output/Aligned.sortedByCoord.out.bam

4. 细胞条形码

提取有效的细胞条形码：

# 从 Cell Ranger
cut -d',' -f2 filtered_feature_bc_matrix/barcodes.tsv.gz | gunzip > barcodes.txt

# 从 STAR
gunzip -c Solo.out/Gene/filtered/barcodes.tsv.gz > barcodes.txt

运行 Velocyto

单样本

velocyto run \
  -v \  # 详细输出
  -m GRCh38_rmsk.gtf \  # 重复序列屏蔽文件
  -b barcodes.txt \  # 细胞条形码
  -o output_dir \  # 输出目录
  -@ 12 \  # 线程数
  possorted_genome_bam.bam \  # 输入 BAM
  genes.gtf  # 基因注释

输出： output_dir/sample.loom
时间： 约 30 分钟处理 5K 细胞

批量处理

#!/bin/bash
# process_velocyto.sh

# 待处理的样本
samples=(sample1 sample2 sample3)

# 参考文件
REPEAT_MASK=~/RefData/GRCh38_rmsk.gtf
GTF=~/RefData/Homo_sapiens.GRCh38.110.gtf

for sample in "${samples[@]}"; do
  echo "Processing $sample..."
  
  # 提取条形码
  gunzip -c ${sample}/filtered_feature_bc_matrix/barcodes.tsv.gz > ${sample}/barcodes.txt
  
  # 运行 velocyto
  velocyto run \
    -m $REPEAT_MASK \
    -b ${sample}/barcodes.txt \
    -o velocyto_output/${sample} \
    -@ 8 \
    ${sample}/possorted_genome_bam.bam \
    $GTF
    
  echo "$sample complete"
done

第二部分：scVelo（分析与可视化）

安装

# 安装 scVelo 及其依赖
pip install -U scvelo loompy

# 可选：用于增强绘图的 R 集成
conda install r-base
pip install rpy2

基础分析流程

import scvelo as scv
import scanpy as sc

# 加载数据
adata = scv.read("sample.loom", cache=True)

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

# 计算速度
scv.tl.velocity(adata, mode='stochastic')  # 或 'dynamical'
scv.tl.velocity_graph(adata)

# 可视化
scv.pl.velocity_embedding_stream(adata, basis='umap', color='celltype')

与现有 Seurat/Scanpy 分析集成

import scvelo as scv
import scanpy as sc

# 加载你已处理的 Seurat/Scanpy 对象
adata_processed = sc.read_h5ad("processed_data.h5ad")

# 加载 velocyto 计数
adata_velo = scv.read("sample.loom", cache=True)

# 与处理后的数据合并
adata = scv.utils.merge(adata_processed, adata_velo)

# 使用现有的 UMAP/聚类
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
scv.tl.velocity(adata, mode='dynamical')
scv.tl.velocity_graph(adata)

# 在现有 UMAP 上绘图
scv.pl.velocity_embedding_stream(
    adata, 
    basis='umap', 
    color='celltype',
    figsize=(8, 6),
    arrow_size=3,
    arrow_length=3
)

速度模式

模式	描述	使用时机
stochastic	简单稳态模型	快速探索，大多数数据集
dynamical	完整动力学模型	精确计时，复杂动态
deterministic	无噪声模型	大型、高质量数据集

# 随机模式（快速，良好的默认选择）
scv.tl.velocity(adata, mode='stochastic')

# 动力学模式（较慢，更准确）
scv.tl.recover_dynamics(adata, n_jobs=8)
scv.tl.velocity(adata, mode='dynamical')

结果解读

# 速度置信度
scv.tl.velocity_confidence(adata)
scv.pl.scatter(adata, c='velocity_confidence', cmap='coolwarm')

# 识别驱动基因
scv.tl.rank_velocity_genes(adata, groupby='celltype')
scv.pl.velocity_genes(adata, var_names=['TOP1', 'TOP2A', 'MKI67'])

# 速度伪时间
scv.tl.velocity_pseudotime(adata)
scv.pl.scatter(adata, c='velocity_pseudotime', cmap='gnuplot')

高级工作流程

从 STARsolo 获取 RNA Velocity

使用 STARsolo 直接生成已剪接和未剪接计数：

STAR --runThreadN 16 \
  --genomeDir ./STARindex \
  --readFilesIn R2.fastq.gz R1.fastq.gz \
  --readFilesCommand zcat \
  --soloType CB_UMI_Simple \
  --soloCBwhitelist 3M-february-2018.txt \
  --soloUMIlen 12 \
  --soloFeatures Gene Velocyto \  # 启用 velocyto 模式
  --outSAMtype BAM SortedByCoordinate

输出： Solo.out/Velocyto/filtered/ 包含已剪接/未剪接矩阵。

转换为 Loom 格式

import scvelo as scv
import scanpy as sc

# 加载 STARsolo velocyto 输出
spliced = sc.read_mtx("Solo.out/Velocyto/filtered/spliced.mtx.gz").T
unspliced = sc.read_mtx("Solo.out/Velocyto/filtered/unspliced.mtx.gz").T
ambiguous = sc.read_mtx("Solo.out/Velocyto/filtered/ambiguous.mtx.gz").T

# 合并到 AnnData 中
adata = sc.AnnData(spliced)
adata.layers['spliced'] = spliced.X
adata.layers['unspliced'] = unspliced.X
adata.layers['ambiguous'] = ambiguous.X

# 添加细胞/基因名称
barcodes = pd.read_csv("Solo.out/Velocyto/filtered/barcodes.tsv", header=None)
genes = pd.read_csv("Solo.out/Velocyto/filtered/features.tsv", sep="\t", header=None)
adata.obs_names = barcodes[0].values
adata.var_names = genes[0].values

# 保存为 loom 文件
adata.write_loom("sample.loom")

可视化示例

流图（推荐）

scv.pl.velocity_embedding_stream(
    adata, 
    basis='umap',
    color='celltype',
    palette='Set2',
    size=80,
    alpha=0.5,
    arrow_size=2,
    linewidth=1.5,
    density=0.8
)

网格图

scv.pl.velocity_embedding_grid(
    adata,
    basis='umap',
    color='celltype',
    arrow_length=3,
    arrow_size=2,
    density=0.5
)

相图（基因特异性）

# 显示特定基因
scv.pl.velocity(adata, ['MKI67', 'TOP2A', 'PCNA'], ncols=3)

故障排除

速度置信度低

症状： 箭头指向随机方向
原因：

内含子覆盖度不足（测序深度低）
稳态细胞（无活跃转换）
数据质量差

修复：

# 过滤低质量的速度基因
scv.tl.velocity_confidence(adata)
adata = adata[:, adata.var.velocity_confidence > 0.05]

Velocyto 因“无有效细胞”而崩溃

原因： BAM 文件与条形码列表中的条形码不匹配
修复：

# 确保条形码与 BAM 标签完全匹配
samtools view possorted_genome_bam.bam | grep -o "CB:Z:[ACGT]*" | head
# 与 barcodes.txt 格式比较

内存错误

原因： 数据集过大
修复：

# 分批处理
scv.settings.n_jobs = 1  # 减少并行化
scv.pp.filter_genes(adata, min_shared_counts=30)  # 更严格的过滤

比较：Velocyto 与 scVelo

工具	用途	语言	优势
Velocyto	从 BAM 计数内含子	Python/CLI	处理任何比对器输出
scVelo	分析速度	Python	高级模型，美观的绘图

工作流程： 始终使用 Velocyto → scVelo 或 STARsolo → scVelo。

资源

scVelo 教程： scvelo.readthedocs.io
Velocyto 文档： velocyto.org
scVelo 论文： Bergen et al., Nature Biotech 2020
RNA velocity 理论： La Manno et al., Nature 2018

另请参阅：

RNA Velocity