企业官网建设流程全解析

MultiVI 多组学整合分析实战：基于 scvi-tools 的 Multiome RNA+ATAC 联合建模指南

【免费下载链接】knowledge-work-pluginsOpen source repository of plugins primarily intended for knowledge workers to use in Claude Cowork项目地址: https://gitcode.com/GitHub_Trending/kn/knowledge-work-plugins

本篇指南以 scvi-tools 生态中的MultiVI深度生成模型为核心，系统讲解如何在同一批细胞上同时分析 RNA 表达与 ATAC 染色质可及性数据（即 10x Multiome 实验产出）。文中完整覆盖从 MuData 数据准备、HVG/峰过滤、setup_mudata注册、模型训练、联合潜空间聚类，到缺失模态补全（imputation）与差异表达/差异可及性分析的全流程，并配套可复制的 Python 代码与仓库内 CLI 脚本用法，读者完成后可直接在自己的多组学数据集上落地。

MultiVI 是什么：一次理解模型定位

MultiVI 是 scvi-tools 家族中专用于multiome 数据（同一批细胞的 RNA-seq + ATAC-seq）的深度生成模型（deep generative model）。在原仓库的模型选型表中，它的定位非常明确：

模型选型完整说明见 bio-research/skills/scvi-tools/SKILL.md。与 PeakVI、scVI 各自只处理单一模态不同，MultiVI 的能力集中在四个方面：

• 学习跨模态的联合潜空间表征（joint latent representation）：把 RNA 与 ATAC 两套特征映射到同一个低维空间，便于下游统一聚类与可视化；
• 处理缺失模态：可以接受只有 RNA 或只有 ATAC 的细胞，模型在训练中自动利用配对信息学习缺失模态的分布；
• 跨实验批次校正：通过batch_key指定批次列，消除技术差异；
• 缺失模态补全（imputation）：训练完成后，可为 ATAC-only 细胞补全 RNA 表达量，为 RNA-only 细胞补全染色质可及性。

环境准备与依赖安装

MultiVI 依赖 scvi-tools、scanpy、mudata（MuData 支持）与 PyTorch。仓库的 环境搭建参考 给出了多种安装方式，其中支持 Multiome 的安装组合为：

# Conda 环境（推荐） conda create -n scvi-env python=3.10 conda activate scvi-env # 核心依赖 pip install scvi-tools scanpy leidenalg mudata muon # GPU 加速（大样本强烈建议） pip install torch --index-url https://download.pytorch.org/whl/cu118

验证安装是否就绪，并确认 GPU 可用：

import scvi import torch import scanpy as sc import mudata as md print(f"scvi-tools version: {scvi.__version__}") print(f"scanpy version: {sc.__version__}") print(f"PyTorch version: {torch.__version__}") print(f"GPU available: {torch.cuda.is_available()}")

环境参考文档中同时强调：MultiVI 依赖 scvi-tools ≥ 1.0（1.x 之后的 API 中数据注册统一走scvi.model.MULTIVI.setup_mudata，不再是 0.x 时代的scvi.data.setup_anndata），mudata ≥ 0.2.0。如果遇到setup_mudata找不到、MuData 读写报错，优先检查这两个版本。

数据格式：MuData 还是分离的 AnnData

MultiVI 官方推荐以MuData作为输入容器，因为它在内部显式维护rna与atac两个模态，细胞一一对应，天然契合多组学数据的结构。

方式一：MuData（推荐）

import mudata as md # 加载 multiome 数据 mdata = md.read("multiome.h5mu") # 结构约定： # mdata.mod['rna'] - 存放 RNA 计数的 AnnData # mdata.mod['atac'] - 存放 ATAC 计数的 AnnData print(f"RNA: {mdata.mod['rna'].shape}") print(f"ATAC: {mdata.mod['atac'].shape}")

方式二：分离的 AnnData 对象

如果上游分别输出了rna.h5ad与atac.h5ad，需要先手动对齐细胞并保持顺序一致，这是最容易出错的一步：

import scanpy as sc adata_rna = sc.read_h5ad("rna.h5ad") adata_atac = sc.read_h5ad("atac.h5ad") # 确保两个模态使用相同细胞且顺序一致 common_cells = adata_rna.obs_names.intersection(adata_atac.obs_names) adata_rna = adata_rna[common_cells].copy() adata_atac = adata_atac[common_cells].copy()

仓库的 train_model.py 在实现train_multivi时对输入格式做了硬性校验：MultiVI requires MuData format with 'rna' and 'atac' modalities。也就是说，最终交给模型的一定是 MuData；从分离 AnnData 起步时，必须在第 3 步之前把它们合并为 MuData。

Step 1：RNA 模态预处理

RNA 侧沿用标准 scvi-tools 流水线，核心要点是先存原始计数，再做 HVG 选择。仓库的 prepare_data.py 与 model_utils.py 中的prepare_adata()封装了同一套逻辑：QC 过滤（min_genes=200、max_genes=5000、线粒体比例 < 20%）→ 把原始计数存入layers["counts"]→ 批次感知的 HVG 选择 → 子集到 HVG。

import scanpy as sc # RNA 预处理（标准 scvi-tools 流水线） adata_rna = mdata.mod['rna'].copy() # 过滤细胞与基因 sc.pp.filter_cells(adata_rna, min_genes=200) sc.pp.filter_genes(adata_rna, min_cells=3) # 关键：在任何归一化之前保存原始计数 adata_rna.layers["counts"] = adata_rna.X.copy() # HVG 选择 sc.pp.highly_variable_genes( adata_rna, n_top_genes=4000, flavor="seurat_v3", layer="counts", batch_key="batch" # 多批次时必填 ) # 子集到高变基因 adata_rna = adata_rna[:, adata_rna.var['highly_variable']].copy()

关于参数取值：仓库在 data_preparation.md 中说明 scvi-tools 在 1,500–5,000 个 HVG 区间内表现最佳，SKILL.md 建议 2,000–4,000；多批次数据推荐flavor="seurat_v3"+batch_key，单批次数据也可以退化为先normalize_total+log1p再选 HVG（prepare_data.py中正是这样做的），选完 HVG 后必须把X恢复为原始计数，保证layers["counts"]与X一致。

Step 2：ATAC 模态预处理

ATAC 侧是峰（peak）× 细胞的可及性矩阵。与 RNA 的整数计数不同，MultiVI 对 ATAC 输入要求二值化（1 表示可及，0 表示不可及），这也是仓库train_model.py中train_peakvi对 scATAC 数据的处理方式（if adata.X.max() > 1: adata.X = (adata.X > 0).astype(np.float32)）。

import numpy as np # ATAC 预处理 adata_atac = mdata.mod['atac'].copy() # 过滤低质量峰 sc.pp.filter_genes(adata_atac, min_cells=10) # 二值化可及性 adata_atac.X = (adata_atac.X > 0).astype(np.float32) # 峰数量过多时，按总可及性取前 N 个峰 if adata_atac.n_vars > 50000: peak_accessibility = np.array(adata_atac.X.sum(axis=0)).flatten() top_peaks = np.argsort(peak_accessibility)[-50000:] adata_atac = adata_atac[:, top_peaks].copy() # 存入 layer adata_atac.layers["counts"] = adata_atac.X.copy()

峰数量阈值可以依据内存与 GPU 显存调整：n_top_peaks默认 50,000，若训练不稳定可适当减少。注意此时 ATAC 的layers["counts"]存储的是二值化后的结果，这与 RNA 的原始整数计数含义不同，但 MultiVI 会按模态区分处理（RNA 走负二项式似然、ATAC 走伯努利似然）。

Step 3：合并为组合 MuData

把预处理后的两个模态封装回 MuData，并再次确认细胞对齐：

import mudata as md # 再次确保匹配细胞 common_cells = adata_rna.obs_names.intersection(adata_atac.obs_names) adata_rna = adata_rna[common_cells].copy() adata_atac = adata_atac[common_cells].copy() # 创建 MuData mdata = md.MuData({ "rna": adata_rna, "atac": adata_atac }) print(f"Combined multiome: {mdata.n_obs} cells") print(f"RNA features: {mdata.mod['rna'].n_vars}") print(f"ATAC features: {mdata.mod['atac'].n_vars}")

Step 4：Setup MultiVI（数据注册）

MultiVI 使用setup_mudata而非setup_anndata。这里的modalities参数完成模态命名到模型角色的映射：告诉模型 RNA 数据在哪、ATAC 数据在哪、batch 信息挂在哪个模态上。

import scvi scvi.model.MULTIVI.setup_mudata( mdata, rna_layer="counts", atac_layer="counts", batch_key="batch", # 可选，多批次时必填 modalities={ "rna_layer": "rna", "batch_key": "rna", "atac_layer": "atac" } )

该调用与 train_model.py 中train_multivi的实现完全一致（rna_layer="counts"、atac_layer="counts"、modalities={"rna_layer": "rna", "batch_key": "rna", "atac_layer": "atac"}），说明这是仓库内验证过的标准注册范式。batch_key缺省时模型退化为纯多组学整合（不做批次校正）；一旦提供，MultiVI 会在联合潜空间中去掉批次带来的技术偏移。

Step 5：训练 MultiVI

模型架构参数集中在n_latent（潜空间维度）与编码器/解码器层数上，训练流程与仓库 train_model.py 的train_multivi保持一致：

# 创建模型 model = scvi.model.MULTIVI( mdata, n_latent=20, n_layers_encoder=2, n_layers_decoder=2 ) # 训练 model.train( max_epochs=300, early_stopping=True, early_stopping_patience=10, batch_size=128 ) # 检查训练收敛情况（ELBO 曲线） model.history['elbo_train'].plot()

参数建议：n_latent常用 20–50（仓库train_model.py的默认值是 20）；max_epochs300 是文档基线，配合early_stopping_patience=10可避免过拟合；batch_size根据显存调整。训练完成后可以通过model.history观察训练/验证 ELBO 是否收敛——仓库 model_utils.py 中的plot_training_history()还额外绘制了 reconstruction loss 曲线，方便判断两个模态的重建质量。

模型与训练产物的持久化可直接复用仓库 CLI：

# 训练 MultiVI（MuData 输入，自动保存 model/ 与 adata_trained.h5ad） python bio-research/skills/scvi-tools/scripts/train_model.py multiome.h5mu results/ \ --model multivi --batch-key batch --n-latent 20 --max-epochs 300

Step 6：联合潜空间与聚类

训练完成后，get_latent_representation()返回每个细胞的联合表征（形状为 细胞 ×n_latent），将其写入 MuData 的obsm，即可无缝衔接 scanpy 的下游分析：

# 联合潜空间表征 latent = model.get_latent_representation() # 写入 MuData mdata.obsm["X_MultiVI"] = latent # 在联合空间上聚类 sc.pp.neighbors(mdata, use_rep="X_MultiVI") sc.tl.umap(mdata) sc.tl.leiden(mdata, resolution=1.0) # 可视化 sc.pl.umap(mdata, color=['leiden', 'batch'], ncols=2)

注意这里邻居图建立在X_MultiVI上，而不是默认的 PCA——这正是 MultiVI 联合整合的价值所在：同一个 UMAP 上，RNA 与 ATAC 的生物学信号被统一呈现，同时批次被校正混合。

仓库的 cluster_embed.py 会自动检测X_MultiVI（其候选列表["X_scANVI", "X_scVI", "X_totalVI", "X_PeakVI", "X_MultiVI"]中包含该键），并自动输出umap_clusters.png与cluster_counts.csv：

python bio-research/skills/scvi-tools/scripts/cluster_embed.py results/adata_trained.h5ad results/ \ --resolution 1.0 --batch-key batch

Step 7：模态特异分析

补全缺失模态（Imputation）

这是 MultiVI 最具特色的能力：跨模态推断。当要整合 ATAC-only 数据集时，可以用训练好的模型为这些细胞补全 RNA 表达量；反过来，为 RNA-only 细胞补全可及性。

# 为 ATAC-only 细胞补全 RNA 表达量 imputed_rna = model.get_normalized_expression( modality="rna" ) # 为 RNA-only 细胞补全可及性 imputed_atac = model.get_accessibility_estimates()

差异表达与差异可及性

基于联合潜空间的分组（如 Leiden 聚类结果）直接做差异检验，RNA 侧与 ATAC 侧各有一套方法：

# 差异表达（RNA） de_results = model.differential_expression( groupby="leiden", group1="0", group2="1" ) # 差异可及性（ATAC） da_results = model.differential_accessibility( groupby="leiden", group1="0", group2="1" )

仓库的 differential_expression.py 将该能力脚本化，支持一次性对全部分组做 one-vs-rest 比较、按lfc_mean取 top N 基因，并可生成火山图：

# 全部分组 one-vs-rest 的差异表达 python bio-research/skills/scvi-tools/scripts/differential_expression.py model/ adata.h5ad de.csv \ --groupby leiden --n-genes 50 --plot # 指定两组比较 python bio-research/skills/scvi-tools/scripts/differential_expression.py model/ adata.h5ad de.csv \ --groupby cell_type --group1 "T cells" --group2 "B cells"

需要说明的是：该 CLI 当前加载模型的分支（SCVI/SCANVI/TOTALVI.load）尚未覆盖 MultiVI，MultiVI 场景下建议直接以 Python API 调用model.differential_expression()。

处理部分模态数据（Partial Data）

MultiVI 支持把不同模态完整度的数据集合并训练：部分细胞同时有 RNA+ATAC（paired），部分只有 RNA，部分只有 ATAC。做法是在obs中标记缺失模态，将缺失侧的数据置为缺失（NaN），MultiVI 会在训练中自动学习跨模态关联：

# 数据集 1：完整 multiome # 数据集 2：只有 RNA # 数据集 3：只有 ATAC # 标记模态 mdata.obs['modality'] = 'paired' # 双模态细胞 # RNA-only 细胞：ATAC 数据置为缺失/NaN # ATAC-only 细胞：RNA 数据置为缺失/NaN # MultiVI 在训练中自动处理缺失模态

这一特性使得「多组学配对数据 + 历史单模态数据」的联合建模成为可能，是全流程中最能体现 MultiVI 相对朴素拼接方法优势的一步。

完整流水线函数

以下函数整合了上述全部步骤（对齐细胞 → 双模态预处理 → MuData → setup → 训练 → 潜空间 → 聚类），可直接复制使用：

def analyze_multiome( adata_rna, adata_atac, batch_key=None, n_top_genes=4000, n_top_peaks=50000, n_latent=20, max_epochs=300 ): """ Complete multiome analysis with MultiVI. Parameters ---------- adata_rna : AnnData RNA count data adata_atac : AnnData ATAC peak data batch_key : str, optional Batch column name n_top_genes : int Number of HVGs for RNA n_top_peaks : int Number of top peaks for ATAC n_latent : int Latent dimensions max_epochs : int Maximum training epochs Returns ------- MuData with joint representation """ import scvi import scanpy as sc import mudata as md import numpy as np # Get common cells common_cells = adata_rna.obs_names.intersection(adata_atac.obs_names) adata_rna = adata_rna[common_cells].copy() adata_atac = adata_atac[common_cells].copy() # RNA preprocessing sc.pp.filter_genes(adata_rna, min_cells=3) adata_rna.layers["counts"] = adata_rna.X.copy() if batch_key: sc.pp.highly_variable_genes( adata_rna, n_top_genes=n_top_genes, flavor="seurat_v3", layer="counts", batch_key=batch_key ) else: sc.pp.normalize_total(adata_rna, target_sum=1e4) sc.pp.log1p(adata_rna) sc.pp.highly_variable_genes(adata_rna, n_top_genes=n_top_genes) adata_rna.X = adata_rna.layers["counts"].copy() adata_rna = adata_rna[:, adata_rna.var['highly_variable']].copy() # ATAC preprocessing sc.pp.filter_genes(adata_atac, min_cells=10) adata_atac.X = (adata_atac.X > 0).astype(np.float32) if adata_atac.n_vars > n_top_peaks: peak_acc = np.array(adata_atac.X.sum(axis=0)).flatten() top_idx = np.argsort(peak_acc)[-n_top_peaks:] adata_atac = adata_atac[:, top_idx].copy() adata_atac.layers["counts"] = adata_atac.X.copy() # Create MuData mdata = md.MuData({"rna": adata_rna, "atac": adata_atac}) # Setup and train scvi.model.MULTIVI.setup_mudata( mdata, rna_layer="counts", atac_layer="counts", batch_key=batch_key, modalities={"rna_layer": "rna", "batch_key": "rna", "atac_layer": "atac"} ) model = scvi.model.MULTIVI(mdata, n_latent=n_latent) model.train(max_epochs=max_epochs, early_stopping=True) # Add representation mdata.obsm["X_MultiVI"] = model.get_latent_representation() # Cluster sc.pp.neighbors(mdata, use_rep="X_MultiVI") sc.tl.umap(mdata) sc.tl.leiden(mdata) return mdata, model # 使用示例 mdata, model = analyze_multiome( adata_rna, adata_atac, batch_key="sample" ) sc.pl.umap(mdata, color=['leiden', 'sample'])

该函数与 prepare_data.py 的差异在于：后者面向单模态 AnnData（train_model.py流水线），而analyze_multiome专门处理双模态的 MuData 组装与 MultiVI 注册，二者互为补充。

Peak-to-Gene 关联分析

利用训练好的模型可以进一步在潜空间做峰-基因关联（Peak-to-Gene Linking），用于识别候选调控关系。思路是：对每个基因，取附近（如 TSS 上下游 100 kb）的峰，计算补全后的可及性与基因表达之间的相关性：

# 在潜空间基于相关性把 ATAC 峰关联到基因，用于识别调控关系 def link_peaks_to_genes(model, mdata, distance_threshold=100000): """ Link peaks to nearby genes based on correlation. Parameters ---------- model : MULTIVI Trained model mdata : MuData Multiome data distance_threshold : int Maximum distance (bp) to link peak to gene Returns ------- DataFrame of peak-gene links """ # 获取补全后的值 rna_imputed = model.get_normalized_expression() atac_imputed = model.get_accessibility_estimates() # 对基因启动子附近的峰，计算峰可及性与基因表达的相关性 # ... （需要基因组坐标信息） return peak_gene_links

这一分析需要额外的基因组坐标信息（peak 的染色体位置、基因的 TSS 坐标），建议结合 PeakVI 参考（atac_peakvi.md）中的染色质分析流程一起使用。

常见问题排查

针对上表可补充两点仓库证据：其一，「双模态细胞数不一致」同样可能出现在setup_mudata阶段——如果传入了不同 cell 集合的 MuData，模型会报错或静默错位，因此 Step 3 中的intersection对齐是硬性要求；其二，「内存不足」问题在 validate_adata.py 中也有对应检查逻辑——当基因数超过 30,000 时它会建议子集到 2,000–4,000 HVG，并提示adata.X = adata.X.toarray()这类稀疏矩阵转换方案。

另外，在正式训练前建议先跑一遍数据校验脚本，它能自动检查整数计数、NaN/负值、批次列是否存在、HVG 是否已选择等硬性条件：

# 校验 AnnData 与 scvi-tools 的兼容性，并给出模型建议 python bio-research/skills/scvi-tools/scripts/validate_adata.py multiome_rna.h5ad \ --batch-key batch --suggest

关键参考

• Ashuach et al. (2023) "MultiVI: deep generative model for the integration of multimodal data"（MultiVI 原始论文）
• scvi-tools 技能总览：模型选型表与决策树（多模态数据 → MultiVI）
• 数据准备参考：HVG 选择、setup_anndata参数详解
• 环境搭建参考：版本兼容性与 GPU 配置
• 仓库脚本：train_model.py（train_multivi）、prepare_data.py、cluster_embed.py、model_utils.py、validate_adata.py

【免费下载链接】knowledge-work-pluginsOpen source repository of plugins primarily intended for knowledge workers to use in Claude Cowork项目地址: https://gitcode.com/GitHub_Trending/kn/knowledge-work-plugins