Velocyto RNA 速率 — 计算 RNA 剪接速率推断细胞状态走向¶

一句话说明¶

Velocyto 是 RNA velocity 分析的起点工具，通过区分未剪接（nascent）和已剪接（mature）mRNA，计算每个基因的转录速率，预测细胞未来的状态走向。

安装与配置¶

# 激活生信分析环境
conda activate bioinfo

# 安装 Velocyto（命令行版本，用于从 BAM 文件生成 loom 文件）
pip install velocyto

# 安装必要依赖
pip install numpy scipy matplotlib pandas anndata

# 安装 pysam（处理 BAM 文件）
conda install -c bioconda pysam -y

# 安装 samtools（BAM 文件操作）
conda install -c bioconda samtools -y

# 验证安装
velocyto --version

核心用法¶

# ── 方法1：针对 10X Genomics 数据（最常用） ──────────────
velocyto run10x \
  /path/to/cellranger_output/ \   # Cellranger 输出目录（含 bam/）
  /path/to/refdata-gex-GRCh38/genes/genes.gtf  # 参考基因组 GTF 文件

# ── 方法2：通用 BAM 文件分析 ──────────────────────────────
velocyto run \
  -b barcodes.tsv \               # 有效 barcode 列表
  -o output_loom/ \               # 输出目录
  -m repeat_msk.gtf \             # 重复序列 GTF（可选，用于过滤）
  sample.bam \                    # 输入 BAM 文件（需要排序和索引）
  genes.gtf                       # 参考基因组 GTF 文件

# ── 方法3：Smart-seq2 数据 ────────────────────────────────
velocyto run \
  -o output_loom/ \
  smartseq2_sample.bam \
  genes.gtf

参数详解¶

参数	说明	建议值
`-b`	有效细胞 barcode 列表文件	10X filtered_barcodes.tsv
`-m`	重复序列掩码文件（减少假阳性）	UCSC rmsk.gtf
`-o`	输出 loom 文件保存目录	自定义
`--dtype`	计数矩阵数据类型（uint32 节省内存）	`uint32`
`--samtools-threads`	BAM 解析并行线程数	4-8
`--samtools-memory`	每线程内存限制（MB）	4096

实战案例¶

import velocyto as vcy             # 导入 velocyto Python 接口
import scvelo as scv               # 通常配合 scVelo 使用
import scanpy as sc
import anndata as ad

# ── 在 Python 中加载和分析 Loom 文件 ─────────────────────

# 1. 读取 velocyto 生成的 loom 文件
loom_file = 'sample.loom'
# 使用 scVelo 读取（推荐，功能更全）
adata = scv.read(loom_file, cache=True)      # 读取并缓存，避免重复解析

# 2. 合并多个样本的 loom 文件
import loompy
# 先用命令行合并（推荐）
# loompy merge sample1.loom sample2.loom merged.loom

# 或在 Python 中用 scVelo 合并
loom_files = ['sample1.loom', 'sample2.loom']
adata = scv.read(loom_files[0], cache=True)

# 3. 与已有的 Scanpy/Seurat 注释结合
# 已有注释的 h5ad 文件（含聚类、UMAP 等）
adata_scanpy = sc.read_h5ad('annotated_adata.h5ad')

# 按 barcode 合并（velocyto loom 的细胞 barcode 需要去掉前缀）
adata.obs_names = [x.split(':')[1] for x in adata.obs_names]  # 处理 barcode 格式

# 用 Scanpy 的注释更新 velocyto 数据
adata_merge = scv.utils.merge(adata, adata_scanpy)  # 合并元数据

# 4. 可视化剪接/未剪接比例
scv.pl.proportions(adata_merge)    # 查看每个细胞的 spliced/unspliced 比例

# 5. 后续分析移交 scVelo（见 326 教程）
# scVelo 更快更准确，通常 velocyto 只用于生成 loom 文件

# ── Bash 中的实用操作 ───────────────────────────────────

# 检查 BAM 文件是否已经排序（velocyto 需要排序好的 BAM）
samtools view -H sample.bam | grep "SO:"          # 查看排序方式

# 对 BAM 文件排序（如果未排序）
samtools sort -@ 8 -o sorted.bam unsorted.bam     # 8线程排序

# 给 BAM 建索引（velocyto 必需）
samtools index sorted.bam                          # 生成 .bai 索引文件

# 运行 velocyto（10X 数据，使用 8 个线程）
velocyto run10x \
  --samtools-threads 8 \             # 8 个线程加速
  cellranger_output/ \
  genes.gtf

常见报错与解决¶

报错	原因	解决方法
`BAM file is not sorted`	BAM 未排序	`samtools sort -o sorted.bam input.bam`
`No such file: .bai`	缺少 BAM 索引	`samtools index sorted.bam`
`velocyto not found`	未安装或未激活环境	`pip install velocyto`
loom 文件读取报 UnicodeError	文件损坏或格式问题	重新生成 loom 文件
运行速度很慢（>24h）	默认单线程	增加 `--samtools-threads 8`

速查表¶

# 10X 数据一键运行（最常用）
velocyto run10x cellranger_out/ genes.gtf

# 指定 barcode + 线程数
velocyto run -b barcodes.tsv \
             --samtools-threads 8 \
             sample.bam genes.gtf

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

# Python 中读取
import scvelo as scv
adata = scv.read('sample.loom')