前言
DeepVelo是一款基于深度学习的推断RNA速率的软件。开发者声称它能更快更精准地推断RNA速率，事实是这样吗？让读者们跟小编一起一探究竟吧。
Github链接：https://github.com/bowang-lab/DeepVelo
学习本教程之前建议读者事先了解一下RNA速率的基础知识，可以点击这里进行学习。
安装
DeepVelo支持的python版本为3.7-3.9，在此环境下使用下面命令即可安装。
pip install deepvelo
这样CPU版的DeepVelo就安装好啦。但是既然这是个基于深度学习的软件，还是用GPU加速才能充分体现它的优势嘛。要安装GPU版本的DeepVelo，要先把CPU版本的pytorch和dgl都卸载。
pip uninstall pytorch dgl
然后使用以下命令查看当前GPU所需的最大cuda版本
nvidia-smi
Tue Aug 27 09:05:27 2024
+---------------------------------------------------------------------------------------+
| NVIDIA-SMI 535.171.04 Driver Version: 535.171.04 CUDA Version: 12.2 |
|-----------------------------------------+----------------------+----------------------+
| GPU Name Persistence-M | Bus-Id Disp.A | Volatile Uncorr. ECC |
| Fan Temp Perf Pwr:Usage/Cap | Memory-Usage | GPU-Util Compute M. |
| | | MIG M. |
|=========================================+======================+======================|
| 0 Quadro RTX 4000 Off | 00000000:65:00.0 Off | N/A |
| 30% 39C P8 2W / 125W | 6217MiB / 8192MiB | 0% Default |
| | | N/A |
+-----------------------------------------+----------------------+----------------------+
然后根据自己的cuda版本安装合适的pytorch（需要根据自己的配置进行调整，不要直接copy下面的命令），因为软件比较老了，所以建议在不报错的情况下尽可能安装旧一点的版本。
pip install torch==1.12.0+cu113 torchvision==0.13.0+cu113 torchaudio==0.12.0 --extra-index-url https://download.pytorch.org/whl/cu113
dgl也不要安装最新版的，需要>=0.4.3,<0.7，cuda版本就只能从pytorch的cuda往低处试了，比如我这是从cu113开始试到了cu111才安装成功。
pip install "dgl-cu111>=0.4.3,<0.7"
成功安装dgl后，可能还会报错缺乏libcusparse、libnvjitlink或torchdata等依赖包，这些就用conda安装就行。
使用
数据预处理
首先导入需要的包，以及随机种子和可视化相关设定。
# %%
# Import top level libraries, including the deepvelo package
import numpy as np
import scvelo as scv
import torch
from deepvelo import Constants, train
from deepvelo.utils import update_dict, velocity
from deepvelo.utils.plot import compare_plot, statplot
from deepvelo.utils.preprocess import autoset_coeff_s
# fix random seeds for reproducibility
SEED = 123
torch.manual_seed(SEED)
torch.backends.cudnn.deterministic = True
torch.backends.cudnn.benchmark = False
np.random.seed(SEED)
# set options for visualization and verbosity
scv.settings.verbosity = 3 # show errors(0), warnings(1), info(2), hints(3)
scv.settings.set_figure_params(
"scvelo", transparent=False
) # for beautified visualization
%load_ext autoreload
%autoreload 2
使用下面命令下载和导入数据集，如果检测到数据集已经下载，则会直接导入。
adata = scv.datasets.dentategyrus_lamanno()
使用下面两行命令进行数据预处理。首先对数据进行过滤和标准化，然后计算一阶和二阶矩便于进行速率估计。
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_neighbors=30, n_pcs=30)
经典RNA速率分析
使用以下代码计算RNA速率。
scv.tl.velocity(adata)
计算的RNA速率储存在adata对象的layers里面，其中velocity是剪接后mRNA的速率，用于后续的分析。
adata.layers
Layers with keys: ambiguous, matrix, spliced, unspliced, Ms, Mu, velocity, variance_velocity, cell_specific_beta, cell_specific_gamma
计算出速率后，我们需要将RNA速率映射到低维空间中，进而推断细胞间的分化发育轨迹。
scv.tl.velocity_graph(adata)
然后不出意外的话应该要出意外了，会出现这个报错。首先是上面的warning，提示没有安装tqdm和ipywidgets，这个按照提示走或者干脆不要进度条都行，问题不大。
主要在于下面的报错，首先去GitHub上搜一下，搜到了这个问题，这个问题是scVelo内部的问题，只需要将scVelo更新一下就能继续使用scVelo进行RNA速率分析了。但问题在于DeepVelo是不支持更高版本的scVelo的，一旦更新scVelo就不能使用DeepVelo了，所以只能硬着头皮去源代码里看看究竟出了什么问题了。
根据报错提示我们首先找到了velocity_graph.py这个模块，发现了第二层报错是出在第363行代码上
vgraph.compute_cosines(n_jobs=n_jobs, backend=backend)
这里可以看出compute_cosines是vgraph对象的成员函数，所以我们在VelocityGraph类中找到了这个函数的定义，并且发现了报错的是当前模块的第175行代码
res = parallelize(
self._compute_cosines,
range(self.X.shape[0]),
n_jobs=n_jobs,
unit="cells",
backend=backend
)()
这行代码报错的是parallelize这个函数，它是从开发者自定义的包中导入的，它的作用是将VelocityGraph类中的另一个成员函数_compute_cosines回调给每一个细胞的表达数据，用于计算它们的余弦相似度。它有一个参数是as_array，默认设置为True，表示将输出转换为数组，False则是按元组输出。但由于_compute_cosines输出的二维元组中第二维长度不一致（即uncertainties, vals, rows, cols长度不相等），因而转换为数组就会产生上面的报错。所以解决方法是不转换为数组，按照元组输出，也就是将上面的代码改为
res = parallelize(
self._compute_cosines,
range(self.X.shape[0]),
n_jobs=n_jobs,
unit="cells",
backend=backend,
as_array=False
)()
这样报错就解决啦。
然后用下面代码绘制RNA速率图。
scv.pl.velocity_embedding_stream(
adata,
basis="tsne",
color="clusters",
legend_fontsize=9,
dpi=150 # increase dpi for higher resolution
)
也能计算细胞的拟时间并且进行可视化，可以看出左下角的细胞是发育起点。
scv.tl.velocity_pseudotime(adata)
scv.pl.scatter(
adata,
color="velocity_pseudotime",
cmap="gnuplot",
dpi=150,
)
DeepVelo分析
那么，我们再使用DeepVelo进行RNA速率分析看看效果如何。首先设定一些基本配置，包括实验的名字和损失函数相关的超参数设置，如批次大小等。
configs = {
"name": "DeepVelo", # name of the experiment
"loss": {"args": {"coeff_s": autoset_coeff_s(adata), # Automatic setting of the spliced correlation objective
"inner_batch_size": 100}}
}
configs = update_dict(Constants.default_configs, configs)
然后训练深度学习模型用于RNA速率的计算，这一步耗时会比经典RNA速率长。
# initial velocity
velocity(adata, mask_zero=False)
trainer = train(adata, configs)
将RNA速率投影到低维空间，这一步运行时间明显比经典RNA速率缩短。
# velocity plot
scv.tl.velocity_graph(adata, n_jobs=8)
然后绘制RNA速率图。
scv.pl.velocity_embedding_stream(
adata,
basis="tsne",
color="clusters",
legend_fontsize=9,
dpi=150 # increase dpi for higher resolution
)
计算细胞的拟时间并进行可视化。
scv.tl.velocity_pseudotime(adata)
scv.pl.scatter(
adata,
color="velocity_pseudotime",
cmap="gnuplot",
dpi=150,
)
总结
通过对经典RNA速率和DeepVelo进行比较，我们发现二者结果差不多。从时间上，经典RNA速率的速率计算+低维投影总时间约为1分25秒，而DeepVelo相应的时间约为1分11秒，比经典RNA速率快了约16.47%。
但是DeepVelo不支持更高版本的python和scVelo，说明这个包开发出来以后就没有维护了，所以小编并不是很推荐用这个包分析RNA速率，直接用scVelo的最新版效果就很好了，而且也不需要像教程里提到的那样修改源代码，开发者已经把bug修复了。当然如果是有志于做算法开发的小伙伴还是建议借着改bug的机会多阅读源代码，对自己的编程能力会有显著提升。
