scVelo RNA 速率进阶 — 动力学模型推断细胞分化方向¶

# 激活单细胞分析环境
conda activate scanpy_env

# 安装 scVelo（最新版 0.3.4）
pip install scvelo==0.3.4

# 安装可选加速依赖
pip install numba             # 加速核心计算（强烈推荐）
pip install loompy            # 读取 loom 格式文件

# 验证安装
python -c "import scvelo; print(scvelo.__version__)"

# 如果 scVelo 依赖冲突，用 conda 安装
conda install -c bioconda scvelo -y

import scvelo as scv           # 导入 scVelo
import scanpy as sc            # 配合 Scanpy 使用
import matplotlib.pyplot as plt

# ── 设置 ──────────────────────────────────────────────
scv.settings.verbosity = 3     # 输出详细日志（0-4，3最合适调试）
scv.settings.set_figure_params('scvelo', facecolor='white', dpi=120)

# ── 读取数据 ──────────────────────────────────────────
# 方式1：从 velocyto 生成的 loom 文件读取
adata = scv.read('sample.loom', cache=True)

# 方式2：从已处理的 h5ad 读取（已合并剪接/未剪接信息）
adata = sc.read_h5ad('velocity_ready.h5ad')

# ── 预处理 ────────────────────────────────────────────
# 过滤：每个基因至少在 20 个细胞中有剪接和未剪接信号
scv.pp.filter_and_normalize(
    adata,
    min_shared_counts=20,     # 最低共享计数阈值
    n_top_genes=2000           # 选择高变基因数量
)
# 计算 PCA 和近邻图
scv.pp.moments(
    adata,
    n_pcs=30,                 # PCA 主成分数
    n_neighbors=30            # 近邻数（影响平滑程度）
)

import scvelo as scv
import scanpy as sc

# ── 完整 scVelo 动力学分析流程 ───────────────────────────

# 1. 加载数据（假设已有 velocyto loom + Scanpy 注释）
adata_loom = scv.read('sample.loom', cache=True)
adata_annot = sc.read_h5ad('annotated.h5ad')

# 合并 loom 的 velocity 信息和 Scanpy 的注释
# 确保 barcode 格式一致
adata_loom.var_names_make_unique()
adata = scv.utils.merge(adata_loom, adata_annot)

# 2. 过滤和标准化
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

# 3. 估计 RNA velocity（动力学模型，最准确但较慢）
scv.tl.recover_dynamics(
    adata,
    n_jobs=4                  # 使用 4 个 CPU 核心并行计算
)
scv.tl.velocity(
    adata,
    mode='dynamical'          # 动力学模型（比 stochastic 更准确）
)
scv.tl.velocity_graph(
    adata,
    n_jobs=4                  # 计算细胞间转换概率图
)

# 4. 可视化 velocity 箭头（最核心的输出图）
scv.pl.velocity_embedding_stream(
    adata,
    basis='umap',             # 在 UMAP 上画速率流场图
    color='cell_type',        # 按细胞类型着色
    title='RNA Velocity',
    figsize=(8, 6)
)

# 单箭头图（每个细胞一个箭头）
scv.pl.velocity_embedding(
    adata,
    basis='umap',
    arrow_length=3,           # 箭头长度（越大越明显）
    arrow_size=2
)

# 5. 计算伪时间（基于 velocity 的潜伏时间）
scv.tl.latent_time(adata)     # 计算 latent time
scv.pl.scatter(
    adata,
    color='latent_time',      # 按潜伏时间着色（代替 Monocle 的 pseudotime）
    color_map='gnuplot'
)

# 6. 分析 velocity 置信度（检验结果可靠性）
scv.tl.velocity_confidence(adata)
scv.pl.scatter(
    adata,
    c=['velocity_length', 'velocity_confidence'],  # 速率大小和置信度
    cmap='coolwarm',
    perc=[5, 95]              # 过滤掉极端值
)

# 7. 找速率驱动基因（贡献最大的 marker 基因）
scv.tl.rank_velocity_genes(
    adata,
    groupby='cell_type',      # 按细胞类型分组
    min_corr=0.3              # 最低相关性过滤
)
# 显示每个类型的 top3 速率驱动基因
df = scv.DataFrame(adata.uns['rank_velocity_genes']['names'])
print(df.head(3))

# 8. 基因相位图（展示单个基因的 spliced vs unspliced 动力学）
scv.pl.velocity(
    adata,
    ['MYC', 'CDK4'],          # 关注的基因
    ncols=2                   # 每行显示 2 个
)

# 9. 保存结果
adata.write_h5ad('velocity_final.h5ad')
print("RNA Velocity 分析完成！")

# 安装
pip install scvelo==0.3.4 numba

# 标准流程（5步）
scv.pp.filter_and_normalize(adata, min_shared_counts=20)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
scv.tl.velocity(adata, mode='stochastic')   # 快速版（用于探索）
scv.tl.velocity_graph(adata)
scv.pl.velocity_embedding_stream(adata, basis='umap')

# 精准版（发论文用）
scv.tl.recover_dynamics(adata, n_jobs=4)   # 先恢复动力学参数
scv.tl.velocity(adata, mode='dynamical')   # 再用动力学模型