编码期刊俱乐部：在骨髓单细胞数据上跑通 ViVAE

2026年2月11日 · 9 min read · Tutorials

一句话结论： ViVAE 是 Saeys 实验室做的变分自编码器（VAE）降维方法，想解决 t-SNE 和 UMAP 丢失全局结构的问题。我们把官方 scRNA-seq 教程在 12,297 个骨髓细胞上完整跑了一遍。ViVAE 的全局结构保留分最高（0.569，t-SNE 0.419，UMAP 0.410），局部结构保留分却垫底（0.325，t-SNE 0.432，UMAP 0.345）。它有用，但不能无脑替换 UMAP。

ViVAE 想解决什么

t-SNE 和 UMAP 都属于邻居嵌入方法，擅长让每个细胞挨着它最近的邻居，却很难让相距较远的细胞群之间的距离保持真实。UMAP 图上两个簇离得多远，大多说明不了什么。

ViVAE（论文：Interpretable models for scRNA-seq data embedding with multi-scale structure preservation）在 VAE 的训练目标里加了一项 MDS 式的距离损失，让二维布局在各个尺度上都尽量尊重原始距离。按 README 的说法，它：

追求多尺度结构保留，对有发育轨迹、离群细胞群或疑似批次效应的数据尤其有用
提供 encoder indicatrices（编码器指示图），标出嵌入在哪里、朝哪个方向被局部拉伸或压缩
是参数化模型，新细胞可以直接投影，不用重新训练
能和 FlowSOM 联用

配套的 ViScore 负责打分：不用标签时给出 Local/Global 结构保留分，有标签时按细胞群算 xNPE。

下面的流程照着仓库里的 tutorials/example_scrnaseq.ipynb 走。每一步我们都实际跑过（打印格式略有整理），文中所有数字都来自我们自己的运行结果，不是从论文里抄的。

测试环境

机器： Linux 工作站，只用 CPU，16 线程
版本： Python 3.11.9，vivae 0.0.2，viscore 0.0.2，flowsom 0.2.3，torch 2.14.0，scanpy 1.11.5，scikit-learn 1.9.1，umap-learn 0.5.9.post2
数据： CELLxGENE 上的 Tabula Sapiens 骨髓数据，12,297 个细胞，15 种注释细胞类型（.h5ad 文件 274 MB）
开销： 整个流程约 3.5 分钟，内存峰值 10.4 GB

第 1 步：安装

conda create --name ViVAE python=3.11.9
conda activate ViVAE

pip install numpy==1.26.3 pandas==2.2.0 matplotlib==3.8.2 scipy==1.12.0 \
  scikit-learn==1.4.0 scanpy==1.9.8 pynndescent==0.5.11 \
  umap-learn==0.5.6 pyemd==1.0.0

pip install git+https://github.com/saeyslab/ViVAE.git \
  git+https://github.com/saeyslab/ViScore.git \
  git+https://github.com/saeyslab/FlowSOM_Python.git

# 实测必须补上这一行，原因见下
pip install "setuptools<81" umap-learn==0.5.9.post2

照 README 的版本锁定来装，我们踩了两个坑，所以加了最后一行：

import viscore 报 ModuleNotFoundError: No module named 'pkg_resources'。ViScore 依赖的 pyemd 还在导入 pkg_resources，新版 setuptools 已经不带这个模块了。装 setuptools<81 就好。
装 FlowSOM 时 scikit-learn 被顺带升到了 1.9.1，结果 umap-learn 0.5.6 直接崩溃：check_array() got an unexpected keyword argument 'force_all_finite'。换成 umap-learn 0.5.9.post2 能正常运行。不装 FlowSOM 的话，原来的版本组合应该没问题。

计算后端与可重复性。 ViVAE 读三个环境变量，必须在 import vivae 之前设置：

import os
os.environ['VIVAE_DETERMINISTIC'] = '1'  # 默认开启：确定性 PyTorch 运算（更慢）
os.environ['VIVAE_CUDA'] = '0'           # 默认 '1'：有 CUDA 就用
os.environ['VIVAE_MPS'] = '0'            # 默认 '0'：不用 Apple GPU
import vivae as vv

确定性默认是开的，但只要跑在 CUDA 上，包会自动把它关掉。README 建议复现教程结果时用 CPU 后端，我们就是这么跑的。

第 2 步：下载与预处理

wget -O scrnaseq.h5ad https://datasets.cellxgene.cziscience.com/f56e0393-c2f3-42f3-93f1-7c8a4ee49d75.h5ad

import scanpy as sc
import numpy as np

hd = sc.read_h5ad('scrnaseq.h5ad')
labs = np.array(hd.obs['cell_type'])

sc.pp.scale(hd, max_value=10.)                      # z 标准化
sc.tl.pca(hd, svd_solver='arpack', n_comps=100)
pc = hd.obsm['X_pca']

这份文件已经标准化过，所以没有 normalize_total/log1p。换成原始计数矩阵的话，得先跑这两步。教程保留 100 个主成分而不是常见的 50 个，算是更保守的选择。

第 3 步：构建 k-NN 图并降噪

knn = vv.make_knn(x=pc, k=1000, fname='scrnaseq_knn.npy')   # 会存盘，下次直接读取
pc_denoised = vv.smooth(x=pc, knn=knn, k=50, coef=1., n_iter=1)

smooth 会把每个点朝它 k 个最近邻的均值移动 coef 比例的距离。coef=1. 时，每个细胞直接被替换成 50 个邻居的均值。教程这样做，是为了让嵌入的骨架更清晰。

建 1,000 近邻图花了 53 秒，平滑本身不到 1 秒。

第 4 步：训练 ViVAE

import torch
np.random.seed(42)
torch.manual_seed(42)

model = vv.ViVAE(input_dim=pc.shape[1], latent_dim=2, random_state=42)
model.fit(X=pc_denoised, n_epochs=100, batch_size=512, verbose=True)
emb = model.transform(pc_denoised)

CPU 上训练用了 40–52 秒。verbose=True 时每个 epoch 打印一行，各项损失分开列出：

Epoch 1/100	recon: 73.6558	kldiv: 0.7139	geom: 0.0000	egeom: 0.0000	mds: 431.6573	imit: 0.0000
Epoch 2/100	recon: 36.1623	kldiv: 0.3566	geom: 0.0000	egeom: 0.0000	mds: 192.0285	imit: 0.0000
...
Epoch 100/100	recon: 0.0802	kldiv: 0.0259	geom: 0.0000	egeom: 0.0000	mds: 0.8671	imit: 0.0000

默认权重下起作用的是三项：重构损失（lam_recon=1）、与潜空间先验的 KL 散度（lam_kldiv=1）、MDS 距离损失（lam_mds=100，文档字符串推荐 scRNA-seq 用 100，流式/质谱细胞术数据用 10）。几何损失和模仿损失默认关闭。

几个默认值值得记一下，教程改了其中两个：n_epochs=50、batch_size=256、learning_rate=1e-3。教程实际训练 100 个 epoch，批大小 512。

第 5 步：生成 t-SNE 和 UMAP 作对照

from sklearn.manifold import TSNE
import umap

emb_tsne = TSNE(n_components=2, perplexity=30., random_state=42).fit_transform(pc)

model_umap = umap.UMAP(random_state=42)
model_umap.fit(pc)
emb_umap = model_umap.transform(pc)

和教程一致，t-SNE 和 UMAP 用的是未降噪的主成分和默认参数。两者在这台机器上都不到一分钟（t-SNE 约 20 秒，UMAP 约 30 秒）。设了 random_state 之后 UMAP 只能单线程跑。

vv.plot_embedding(emb, labels=labs) 可以把三种嵌入按细胞类型上色画出来。

第 6 步：编码器指示图（encoder indicatrices）

ei = model.encoder_indicatrices(X=pc_denoised, radius=1e-3, n_steps=30, n_polygon=200)
fig, ax = vv.plot_indicatrices(indicatrices=ei, scale_factor=4e2, show=False)

encoder_indicatrices 返回的是一个 EncoderIndicatome 对象，不是每个细胞一个数值。它在潜空间上取 30 × 30 的网格，在每个网格点把输入空间里一个极小的球面映射到二维，得到一个多边形。多边形的大小和朝向，就是嵌入在那个位置把数据拉伸或压缩了多少、朝哪个方向。这些多边形叠在嵌入图上看，是一张形变地图，不是颜色刻度。拿它当每个细胞的颜色用，比如 plt.scatter(c=ei)，会直接报 TypeError。

教程对这张图的解读是：单核细胞和粒细胞所在区域被人为压缩，浆细胞则出现了"失控"式的拉伸。浆细胞的问题在第 8 步的 xNPE 分数里也能看到。

第 7 步：用 ViScore 做无监督评估

import viscore as vs

for name, e in [('t-SNE', emb_tsne), ('UMAP', emb_umap), ('ViVAE', emb)]:
    s = vs.score(hd=pc, ld=e)          # 返回字典：Sl、Sg、RNX、QNX
    print(f"{name:<6} Local SP {s['Sl']:.3f}  Global SP {s['Sg']:.3f}")

t-SNE  Local SP 0.432  Global SP 0.419
UMAP   Local SP 0.345  Global SP 0.410
ViVAE  Local SP 0.325  Global SP 0.569

两个分数都是 R_NX 邻域保留曲线下的面积，范围 −1 到 1，0 表示和随机嵌入一样。Local SP 用对数尺度，更看重小邻域；Global SP 对各种邻域大小一视同仁。第一次调用 vs.score 要编译 numba，花了约 20 秒，之后每次不到 1 秒。

结果很清楚：ViVAE 在全局结构上明显领先，在局部结构上排最后。 近邻关系保留得最好的是 t-SNE。换硬件或库版本，具体数字会有小幅变动，但这个格局和教程里报告的一致。

打分时一定要用未降噪的 pc。 ViVAE 是在平滑后的数据上训练的，拿 pc_denoised 当参照会高估它：Global SP 从 0.569 涨到 0.687，Local SP 从 0.325 涨到 0.468。教程里专门提醒过这一点。

第 8 步：哪些细胞群被扭曲了（xNPE）

xNPE（Extended Neighbourhood-Proportion-Error）要用到细胞类型标签。它比较每个细胞群的邻居构成在原始数据和嵌入里的差异，每个细胞群给一个分数：0 表示完美，开启默认的基线校正时，1 相当于随机嵌入的误差水平。

xnpe_umap  = vs.xnpe(hd=pc, ld=emb_umap, annot=labs, knn=knn)   # 字典：细胞群 -> 分数
xnpe_vivae = vs.xnpe(hd=pc, ld=emb,      annot=labs, knn=knn)

for pop in sorted(xnpe_vivae, key=xnpe_vivae.get, reverse=True)[:4]:
    print(f"{pop:<33} {xnpe_vivae[pop]:.3f}")
print("UMAP better than ViVAE for:", [p for p in xnpe_vivae if xnpe_umap[p] < xnpe_vivae[p]])

plasmablast                       0.600
mature NK T cell                  0.363
plasma cell                       0.276
CD8-positive, alpha-beta T cell   0.228
UMAP better than ViVAE for: ['CD8-positive, alpha-beta T cell', 'erythrocyte', 'mature NK T cell', 'memory B cell', 'monocyte', 'neutrophil']

记得传入第 3 步的 knn。不传它（也不指定 k）会报错。两次调用加起来约 1 分钟。vs.plot_xnpe_barplot 可以把结果画成柱状图，vs.neighbourhood_composition_plot 能看出嵌入得不好的细胞群到底混进了哪些邻居。

ViVAE 误差最大的是浆细胞谱系（浆母细胞和浆细胞），和指示图里看到的拉伸对得上。15 个细胞群里有 6 个 UMAP 比 ViVAE 好，其中就有教程也点名的 T/NK 细胞群。

附加：在嵌入上叠加 FlowSOM

FlowSOM 要的是带通道/标记名的 AnnData 对象，所以先把主成分包一层：

import anndata as ad
import pandas as pd
import flowsom as fs

hdf = ad.AnnData(pc)
cols = [f'PC{i+1}' for i in range(pc.shape[1])]
hdf.var = pd.DataFrame({'n': range(pc.shape[1]), 'channel': cols, 'marker': cols})

fsom = fs.FlowSOM(hdf, n_clusters=40, xdim=12, ydim=12, seed=42)
fig, ax = vv.plot_embedding(embedding=emb, labels=labs, s=1., fsom=fsom,
                            dr_model=model, fsom_show_nodes=False, show=False)

聚类是在主成分空间里做的，不是在嵌入上。dr_model=model 让 ViVAE 把 FlowSOM 的节点和最小生成树投影到它的二维布局里，这正是参数化模型的好处。UMAP 也有 .transform，把 dr_model 换成 model_umap 同样能画；t-SNE 就不行。这一步总共约 4 秒。

什么时候用 ViVAE

适合 ViVAE 的情况：

你关心细胞群之间的距离：谱系、轨迹、两群细胞到底隔得多远（这次 Global SP 领先一大截）
你想知道二维图在哪里扭曲了数据（指示图、xNPE）
之后要把新细胞投影进已有布局

继续用 t-SNE/UMAP 的情况：

主要看局部邻域是否紧凑（这次两者的 Local SP 都比 ViVAE 高）
想要更轻的依赖。ViVAE 流程要建 1,000 近邻图、装 PyTorch，还得锁好几个版本才能装干净；12k 细胞的完整运行内存峰值到了 10.4 GB

不管选哪个： 打分时拿未平滑的输入当参照；想根据二维图下结论之前，先看看具体细胞群的误差。

链接

论文： Interpretable models for scRNA-seq data embedding with multi-scale structure preservation（bioRxiv）
代码： ViVAE（教程：tutorials/example_scrnaseq.ipynb）· ViScore · FlowSOM_Python
基准测试： ViScore benchmarking 框架，用于在 HPC 集群上跨数据集比较降维方法

编码期刊俱乐部系列：把已发表的生信工具放到真实数据上跑，如实记录结果。

single-cell dimensionality-reduction machine-learning python coding-journal-club vivae umap