一、写在前面

UNAGI是一种深度生成神经网络，专门用于分析时间序列单细胞转录组数据。该工具捕捉了疾病进展背后的复杂细胞动态，增强了药物扰动建模和筛查能力。当应用于特发性肺纤维化患者的数据集时，UNAGI学习到疾病知情的细胞嵌入，加深我们对疾病进展的理解，从而识别潜在的治疗药物候选。作者利用蛋白质组学验证了UNAGI细胞动力学分析的准确性，而使用纤维化混合剂的人体精切肺切片也证实了UNAGI关于硝苯地平（一种降压药物）可能对人体组织具有抗纤维化作用的预测。UNAGI的多功能性延伸至包括新冠在内的其他疾病，展现出适应性，并证实其在解码特发性肺纤维化以外复杂细胞动态中的更广泛应用，进一步扩大其在多样病理领域寻求治疗解决方案的应用。简单来说UNGAI可以对单细胞数据进行时间序列分析，并预测时间序列中的药物与基因干扰。 测试文件与学习手册见文末。

原文链接：

https://www.nature.com/articles/s41551-025-01423-7

代码参考：

https://github.com/mcgilldinglab/UNAGI/blob/main/tutorials/run_UNAGI_using_example_dataset.ipynb

二、实操流程

如果你的计算资源不够用。不妨看看：

足够支持你完成硕博生涯的生信环境

忘记宣传了，独享用户连技术支持都是独享的

RTX5090、4080S、5070显卡上机

访问链接：https://biomamba.xiyoucloud.net/

1.安装

conda create -n unagi python=3.10
conda activate unagi
pip install scUNAGI -i https://pypi.tuna.tsinghua.edu.cn/simple
# 有一个比较恶心的事情，大家需要修改一下这个脚本的源码
# 从github上克隆包相关文件
git clone https://github.com/mcgilldinglab/UNAGI.git
cd UNAGI
# 在本地安装一遍矫正一些包的冲突
pip install ./
locate UNAGI_analyst.py # 定位脚本位置
# 修改定位到的脚本，下面要换成你自己的路径
sed -i '10d' /home/biomamba/miniconda3/lib/python3.9/site-packages/UNAGI/UNAGI_analyst.py # 会用vim的同学也可以手动删除这行内容：from VTK.IO.Xdmf3.Testing.Python.VToXLoop import raiseErrorAndExit
# 如果上面这个问题不解决，你会遇到这个报错：ModuleNotFoundError: No module named 'VTK'
# 并且，numpy的版本需要<2.0，否则会遇到这个报错：AttributeError: `np.string_` was removed in the NumPy 2.0 release. Use `np.bytes_` instead.
pip install numpy==1.26.4 --force-reinstall
# 安装特定版本的torch组件，否则会提示显卡驱动版本不兼容
pip install torch==2.4.1 torchvision==0.19.1 torchaudio==2.4.1 --index-url https://download.pytorch.org/whl/cu128

另外还需要依赖idrem库完成可视化，我已经在本地给大家下载好了，大家可以随手册一起下载

2. 测试数据集演示

本章节将逐步介绍如何使用UNAGI配置模型架构、调整训练参数以及对时间序列单个数据集进行分析。我们通过将UNAGI应用于从snRNA-seq数据，并演示UNAGI的基本功能。

# 以下包都需要能正常导入：
import warnings
import scanpy as sc
import os
warnings.filterwarnings('ignore')
import warnings
from UNAGI import UNAGI

2.1 初始化并加载数据集

加载好UNAGI之后，第一步就是要设置数据分析参数，UNAGI可以直接加载h5ad数据，不明白这部分内容的可以参考一文了解Python中的单细胞通货：AnnData。UNAGi需要一个拥有时间序列作为输入的数据，这个数据保存在annData.obs中，这个数据需要包含KNN连接矩阵(单细胞数据处理可见：告别R依赖！一文搞定Python单细胞基础分析)。

#提交脚本
samples=("Amp-3219" "Amp-3588" "Amp-6762_combined")
count=1
for sample in "${samples[@]}"; do
    id="S${count}"
    echo "Submitting $sample with ID $id ..."
    qsub /data/CRCdata/WES/add_label.bash $sample $id
    count=$((count + 1))
done
import warnings
import scanpy as sc
import os
adata = sc.read_h5ad('./UNAGI-main/UNAGI/data/example/0.h5ad')
adata.obs['stage']

145_94_D__CAACAGTTCCAAGCCG     0 

221_45_Db__TCCATCGAGGCTCACC    0 

248_111_D__GTAACACAGGCTATCT    0 

145_21_Db__ACAGCCGGTGGTCTGC    0 

248_111_D__GGTGTCGCAATACCCA    0                                                            .. 

221_35_Db__GTGACGCGTTTACACG    0 

302_11_D__GAAGCCCGTTCATCTT     0 

145_22_D__CTTCGGTCATGTACGT     0 

145_21_Db__AGATGAAAGGGAGGTG    0 

94_124_D__CCCAACTTCGTCGATA     0 

Name: stage, Length: 13550, dtype: int64

# 看一下这个adata的全部注释信息：
adata.obs

13550 rows × 31 columns

# 看一下包含的细胞类型：
adata.obs['ident'].unique()
# 可以看出这是个肺部的数据

['VEcapillary', 'VEarterial', 'PericyteAlveolar', 'VEaerocyte', 'VEvenous', 'FibroblastAlveolar', 'SMC', 'FibroblastAdventitial', 'VEaerocyte/VEcapillary/FibroblastAlveolar', 'LymphEndo'] Categories (10, object): ['FibroblastAdventitial', 'FibroblastAlveolar', 'LymphEndo', 'PericyteAlveolar', ..., 'VEaerocyte/VEcapillary/FibroblastAlveolar', 'VEarterial', 'VEcapillary', 'VEvenous']

# 实例化unagi
unagi = UNAGI()
# 预处理h5ad文件：
unagi.setup_data('./UNAGI-main/UNAGI/data/example/0.h5ad',# 单细胞的h5ad文件，注意，unagi会在这个路径下生成临时文件夹（样本名文件夹与，model_save文件夹），因此同一时间运行多个unagi，最好放在不同的目录下
                 total_stage=4,# 时序结果的总序列数
                 stage_key='stage',# 生成的时间序列在adata.obs中的名称
                 threads=20# 默认参与计算的线程数
                 )
unagi

Fast mode enabled!

 Cell graphs found, skipping cell graph construction!

<UNAGI.UNAGI_tool.UNAGI at 0x7f2014bd12d0>

2.2 确认模型结构与超参数

run_UNAGI的输出日志会很长，用jupyter学习本教程的同学可以点击左侧的导航栏向下跳转。 我的演示设备中CPU是Intel Xeon Platinum 8383C、GPU是RTX5990，大约花费20min。

# 设置训练参数：
unagi.setup_training(task='example',# 任务名称
                     dist='ziln',# 当前单细胞数据的分布模式，可选的有：'ziln' (zero-inflated log normal), 'zinb' (zero-inflated negative binomial), 'zig' (zero-inflated gamma), and 'nb' (negative binomial)
                     device='cuda:0',# 指定显卡
                     GPU=True,# 是否启用GPU加速
                     epoch_iter=5,# 第二轮循环开始，每轮的训练epoch数
                     epoch_initial=1,# 第一轮的epoch数
                     max_iter=3,# 训练三轮大循环
                     BATCHSIZE=560 # 最小的批次大小
                     )
# 训练：
unagi.run_UNAGI(idrem_dir = './iDREM')# 需要有idrem库，否则可视化内容无法生成

# UNAGI输出了许多文件，存储在0、1、2、model_savve文件夹中
!tree -L 2 UNAGI-main/UNAGI/data/example/

UNAGI-main/UNAGI/data/example/
├── 0
│   ├── averageValues.npy 
│   ├── edges.txt 
│   ├── idremInput
│   ├── idremModel
│   ├── idremResults
│   ├── idremsetting
│   ├── loss.txt 
│   └── stagedata
├── 0.h5ad
 ├── 1
│   ├── averageValues.npy
 │   ├── edges.txt
 │   ├── idremInput
│   ├── idremModel
│   ├── idremResults
│   ├── idremsetting
│   ├── loss.txt 
│   └── stagedata
├── 1.h5ad 
├── 2
│   ├── averageValues.npy 
│   ├── edges.txt 
│   ├── idremInput
│   ├── idremModel
│   ├── idremResults
│   ├── idremsetting
│   ├── loss.txt 
│   └── stagedata
├── 2.h5ad 
├── 3.h5ad
 └── model_save
    ├── example_0.pth
    ├── example_1.pth
    ├── example_2.pth 
    ├── example_dis_0.pth 
    ├── example_dis_1.pth 
    ├── example_dis_2.pth  
    └── training_parameters.json
  19 directories, 20 files

2.3 虚拟干扰分析

在UNAGI模型训练完毕后，时间序列可初步获得四个stage。原输入h5ad的路径下训练完成后会生成一个新的 “升级版单细胞数据”，就是这里的2/stagedata/dataset.h5ad，其中包含：模型学习到的细胞潜变量空间（z）、UMAP、动态基因、细胞轨迹（Tracks）细胞状态转换图、扰动分析所需的全部信息，可用于执行下游任务。

%%time
import warnings
warnings.filterwarnings('ignore')
from UNAGI import UNAGI
unagi = UNAGI()
# 进行药物干扰分析：
unagi.analyse_UNAGI('./UNAGI-main/UNAGI/data/example/2/stagedata/dataset.h5ad',# 输入文件
                    iteration=2,# 选择第几次迭代计算次数
                    random_background_sampling_times=1000,#  寻找动态标记基因时，随机背景基因抽样次数，数值越大，结果越稳定，但跑的会越慢
                    target_dir=None, # 目标文件夹，none则输出
                    customized_drug='./UNAGI-main/UNAGI/data/cmap_drug_target_genes.npy',# 自定义药物数据库，这里是由作者提供的
                    run_perturbation=True# 是否运行基因、药物干扰
                   )

2.4可视化

UNAGI会生成包含.h5ad、.pkl的心文件夹(为指定的情况下输出的是这里的example_2)，可用于后续可视化。

import scanpy as sc
import pickle as pkl
# 读取输出的h5ad单细胞文件
adata = sc.read('./example_2/dataset.h5ad')
# 单细胞对象：
adata

AnnData object with n_obs × n_vars = 27647 × 2484
     obs: 'barcode', 'nUMI', 'nGene', 'percMito', 'percSpliced.v', 'percUnspliced.v', 'percAmbiguous.v', 'Sample.ID', 'Subject.ID', 'Sex', 'Age', 'Surface.Density', 'Percent.Tissue', 'Disease.Ident', 'nCount_RNA', 'final.ident', 'nFeature_RNA', 'nCount_integrated', 'nFeature_integrated', 'name.simple', 'Library.Ident', 'renamedFin', 'integrated_snn_res.2', 'seurat_clusters', 'renamed.2', 'integrated_snn_res.3', 'log10UMI', 'renamed.3', 'stage', 'leiden', 'ident'
     var: 'features'
     obsm: 'X_umap', 'embdedding', 'umap', 'z'
     layers: 'geneWeight'
     obsp: 'gcn_connectivities'

# 读取非结构性注释信息：
adata.uns = pkl.load(open('./example_2/attribute.pkl', 'rb'))
type(adata.uns)

anndata.compat._overloaded_dict.OverloadedDict

2.4.1 细胞降维可视化

from UNAGI.plotting import cell_type_composition, plot_stages_latent_representation
plot_stages_latent_representation(adata, # 读入的单细胞对象
                                  cell_type_key='ident', # 细胞的分类变量
                                  stage_key='stage'# 分组变量
                                 )

ARI:  0.6551458869846446 
NMIs:  0.6503782470882618 
silhouette score:  0.22718687453535993 
ARI:  0.680303236781184 
NMIs:  0.7374966525638296 
silhouette score:  0.30596853452442374 
ARI:  0.7884503852215209 
NMIs:  0.7823241013103959 
silhouette score:  0.20378125738231517 
ARI:  0.7709702766822403 
NMIs:  0.7719353736314634 
silhouette score:  0.22769615347172187

ARIs:  [0.7237174464173974] 
NMI:  [0.7355335936484877] 
silhouette score:  0.2411582049784552

2.4.2 细胞比例可视化

cell_type_composition(adata, cell_type_key='ident', stage_key='stage')

2.4.3 细胞动态发育轨迹

from UNAGI.utils import visualize_dynamic_graphs_by_text
# 展示各个stage中各细胞类型的亚群发育轨迹：
visualize_dynamic_graphs_by_text(edges=adata.uns['edges'],cell_types=adata.uns['clusterType'])

# 动态轨迹图也可以直接保存为txt格式：
# visualize_dynamic_graphs_by_text(edges=adata.uns['edges'],cell_types=adata.uns['clusterType'],write='saved_dynamic_graphs.txt')

2.4.4 药物干扰结果可视化

可以看出ibudilast 的干扰分值非常高，其与纤维化过程中的TGFβ通路相关，其在肺纤维化中具有治疗潜能。

from UNAGI.perturbations import get_top_compounds, get_top_pathways
get_top_compounds(adata,perturb_change=0.5,top_n=10)

2.4.5通路干扰结果

get_top_pathways(adata,perturb_change=0.5,top_n=10)

三、资料领取

更多编程教程可见：

十小时学会Linux

5.5小时掌握R语言

8小时掌握生信Python

本推送的代码文件可以在以下链接中下载：

通过网盘分享的文件链接:

https://pan.baidu.com/s/1No_C678UKW_SWAjS_V13-Q

提取码: 2it8