跳至章节信息跳至正文
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

高级整合

🧠 关键要点

当 RNA 测序（RNA Sequencing, RNA-seq）、转座酶可及染色质测序（assay for transposase-accessible chromatin with high-throughput sequencing, ATAC-seq）和抗体衍生标签（Antibody-Derived Tag, ADT）数据来自不同实验、没有完整逐细胞配对时，可以利用共有模态、配对桥接数据或先验关系进行整合。

查询到参考映射（query-to-reference mapping）可将新数据接入已有参考，用于细胞类型预测及缺失模态插补（imputation），例如根据 RNA 数据预测蛋白信号。插补结果是依赖参考和模型的估计，不能替代实际测量。

⚙️ 环境设置
步骤
yml

  1. 安装 conda：

    • 在创建环境之前，请确保 conda 已安装在你的系统中。

  2. 保存 yml 内容：

    • 将 yml 选项卡中的内容保存为文件 environment.yml。

  3. 创建环境：

    • 打开终端或命令提示符。

    • 运行以下命令：

      conda env create -f environment.yml

  4. 激活环境：

    • 创建好环境后，使用以下命令激活它：

      conda activate <environment_name>

    • 请将 <environment_name> 替换为 environment.yml 文件中指定的环境名称。该名称在 yml 文件中如下所示：

      name: <environment_name>

  5. 验证安装：

    • 通过运行以下命令，检查环境是否创建成功：

      conda env list

研究动机

本章介绍未配对整合（unpaired integration）、部分模态重叠数据的镶嵌式整合（mosaic integration），以及查询到参考映射（query-to-reference mapping）。当不同模态来自独立实验、并非在同一个细胞中联合测量时，可根据共有模态、配对桥接数据或先验知识建立联系。

下面分别说明这些情形及对应方法。

环境设置

import logging

import anndata
import anndata2ri
import multigrate as mtg
import networkx as nx
import numpy as np
import pandas as pd
import rpy2.rinterface_lib.callbacks
import scanpy as sc
import scglue
from rpy2.robjects import pandas2ri, r
from scvi.model import TOTALVI

rpy2.rinterface_lib.callbacks.logger.setLevel(logging.ERROR)

pandas2ri.activate()
anndata2ri.activate()

%load_ext rpy2.ipython

import warnings

warnings.filterwarnings("ignore")
输出
Global seed set to 0
/lustre/groups/ml01/workspace/anastasia.litinetskaya/miniconda3/envs/advanced-integration/lib/python3.9/site-packages/pytorch_lightning/utilities/warnings.py:53: LightningDeprecationWarning: pytorch_lightning.utilities.warnings.rank_zero_deprecation has been deprecated in v1.6 and will be removed in v1.8. Use the equivalent function from the pytorch_lightning.utilities.rank_zero module instead.
  new_rank_zero_deprecation(
/lustre/groups/ml01/workspace/anastasia.litinetskaya/miniconda3/envs/advanced-integration/lib/python3.9/site-packages/pytorch_lightning/utilities/warnings.py:58: LightningDeprecationWarning: The `pytorch_lightning.loggers.base.rank_zero_experiment` is deprecated in v1.7 and will be removed in v1.9. Please use `pytorch_lightning.loggers.logger.rank_zero_experiment` instead.
  return new_rank_zero_deprecation(*args, **kwargs)
During startup - Warning message:
Setting LC_CTYPE failed, using "C"

原教程为 Seurat 桥接整合（bridge integration）另列了开发分支安装命令：remotes::install_github("satijalab/seurat", "feat/dictionary", quiet = TRUE) 或 devtools::install_github("https://github.com/satijalab/seurat/tree/feat/dictionary")。这些命令对应当时的 feat/dictionary 分支；本章环境指定的 Seurat 5.1 已包含桥接接口，应按所用版本的文档配置，而不是默认必须安装这个旧开发分支。

%%R
suppressPackageStartupMessages({
    library(SingleCellExperiment)
    library(Seurat)
})
set.seed(123)

    WARNING: The R package "reticulate" only fixed recently
    an issue that caused a segfault when used with rpy2:
    https://github.com/rstudio/reticulate/pull/1188
    Make sure that you use a version of that package that includes
    the fix.

数据准备

继续使用 NeurIPS 2021 单细胞竞赛的数据：Multiome 配对获取 RNA 测序（RNA Sequencing, RNA-seq）与转座酶可及染色质测序（assay for transposase-accessible chromatin with high-throughput sequencing, ATAC-seq）信息；转录组与表位联合测序（Cellular Indexing of Transcriptomes and Epitopes by Sequencing, CITE-seq）则配对获取 RNA 与蛋白信号 Luecken et al., 2021。

各示例需要以下输入：

  • 将 Multiome 的 RNA 与 CITE-seq 的抗体衍生标签（Antibody-Derived Tag, ADT）作为未配对数据，使用 图链接统一嵌入（graph linked unified embedding, GLUE）

  • 使用 CITE-seq 的配对 RNA 与蛋白数据构建参考，并完成查询映射；方法为 totalVI

  • 将 CITE-seq 的 RNA 部分作为查询，演示自 Seurat v4 引入的映射流程；对应教程为 Seurat v4

  • 使用 CITE-seq 配对数据作为 桥接 数据，联系 RNA 参考与 ADT 查询

  • 联合 NeurIPS Multiome 和 CITE-seq，构建三模态参考并映射查询；方法为 multigrate

为演示查询映射，先将数据分为参考集与查询集。仅演示整合的方法也使用这里选出的参考批次。下方文件路径属于原作者的共享文件系统，运行时需替换为对应的本地预处理文件。

Multiome 与 CITE-seq 可使用相同的批次名称，但来自不同实验，并不是同一批细胞。批次名称相同不能作为逐细胞配对的依据。

cite_reference_batches = [
    "s1d1",
    "s1d2",
    "s1d3",
]  # need for totalVI, multigrate and bridge (RNA part)
multiome_reference_batches = ["s1d1", "s1d2", "s1d3"]  # need for GLUE and multigrate
cite_query_batches = ["s2d1", "s2d4"]  # need for totalVI, multigrate and bridge
multiome_query_batches = ["s2d1", "s2d4"]  # need for multigrate
rna_multiome = sc.read(
    "/lustre/groups/ml01/workspace/anastasia.litinetskaya/data/neurips-multiome/rna_hvg.h5ad"
)
rna_multiome
AnnData object with n_obs × n_vars = 69249 × 4000 obs: 'GEX_pct_counts_mt', 'GEX_n_counts', 'GEX_n_genes', 'GEX_size_factors', 'GEX_phase', 'ATAC_nCount_peaks', 'ATAC_atac_fragments', 'ATAC_reads_in_peaks_frac', 'ATAC_blacklist_fraction', 'ATAC_nucleosome_signal', 'cell_type', 'batch', 'ATAC_pseudotime_order', 'GEX_pseudotime_order', 'Samplename', 'Site', 'DonorNumber', 'Modality', 'VendorLot', 'DonorID', 'DonorAge', 'DonorBMI', 'DonorBloodType', 'DonorRace', 'Ethnicity', 'DonorGender', 'QCMeds', 'DonorSmoker' var: 'feature_types', 'gene_id', 'highly_variable', 'means', 'dispersions', 'dispersions_norm' uns: 'ATAC_gene_activity_var_names', 'Site_colors', 'batch_colors', 'cell_type_colors', 'dataset_id', 'genome', 'hvg', 'organism' obsm: 'ATAC_gene_activity', 'ATAC_lsi_full', 'ATAC_lsi_red', 'ATAC_umap', 'GEX_X_pca', 'GEX_X_umap', 'X_umap' layers: 'counts'
atac_multiome = sc.read(
    "/lustre/groups/ml01/workspace/anastasia.litinetskaya/data/trimodal_neurips/atac_hvf_muon.h5ad"
)
atac_multiome
AnnData object with n_obs × n_vars = 69249 × 20000 obs: 'GEX_pct_counts_mt', 'GEX_n_counts', 'GEX_n_genes', 'GEX_size_factors', 'GEX_phase', 'ATAC_nCount_peaks', 'ATAC_atac_fragments', 'ATAC_reads_in_peaks_frac', 'ATAC_blacklist_fraction', 'ATAC_nucleosome_signal', 'cell_type', 'batch', 'ATAC_pseudotime_order', 'GEX_pseudotime_order', 'Samplename', 'Site', 'DonorNumber', 'Modality', 'VendorLot', 'DonorID', 'DonorAge', 'DonorBMI', 'DonorBloodType', 'DonorRace', 'Ethnicity', 'DonorGender', 'QCMeds', 'DonorSmoker', 'technology', 'cell_type_l2', 'cell_type_l1', 'cell_type_l3', 'assay', 'split' var: 'feature_types', 'gene_id', 'highly_variable', 'means', 'dispersions', 'dispersions_norm' uns: 'ATAC_gene_activity_var_names', 'dataset_id', 'genome', 'hvg', 'log1p', 'organism' obsm: 'ATAC_gene_activity', 'ATAC_lsi_full', 'ATAC_lsi_red', 'ATAC_umap', 'GEX_X_pca', 'GEX_X_umap' layers: 'binary', 'counts', 'cpm', 'tf-idf-binary', 'tf-idf-counts'
rna_cite = sc.read(
    "/lustre/groups/ml01/workspace/anastasia.litinetskaya/data/neurips-cite/rna_hvg.h5ad"
)
rna_cite
AnnData object with n_obs × n_vars = 90261 × 4000 obs: 'GEX_n_genes_by_counts', 'GEX_pct_counts_mt', 'GEX_size_factors', 'GEX_phase', 'ADT_n_antibodies_by_counts', 'ADT_total_counts', 'ADT_iso_count', 'cell_type', 'batch', 'ADT_pseudotime_order', 'GEX_pseudotime_order', 'Samplename', 'Site', 'DonorNumber', 'Modality', 'VendorLot', 'DonorID', 'DonorAge', 'DonorBMI', 'DonorBloodType', 'DonorRace', 'Ethnicity', 'DonorGender', 'QCMeds', 'DonorSmoker', 'is_train' var: 'feature_types', 'gene_id', 'highly_variable', 'means', 'dispersions', 'dispersions_norm' uns: 'dataset_id', 'genome', 'hvg', 'organism' obsm: 'ADT_X_pca', 'ADT_X_umap', 'ADT_isotype_controls', 'GEX_X_pca', 'GEX_X_umap' layers: 'counts'
adt_cite = sc.read("/lustre/groups/ml01/workspace/daniel.strobl/neurips_cite_pp.h5ad")
adt_cite
AnnData object with n_obs × n_vars = 120502 × 136 obs: 'donor', 'batch', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'n_counts', 'outliers' var: 'gene_ids', 'feature_types', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts' uns: 'batch_colors', 'donor_colors', 'neighbors', 'pca', 'umap' obsm: 'X_isotypes', 'X_pca', 'X_pcahm', 'X_umap' varm: 'PCs' layers: 'counts' obsp: 'connectivities', 'distances'

CITE-seq 的 RNA 与 ADT 对象在 .obs_names 中使用了不同的名称格式。先统一细胞条形码（cell barcode, CB）及批次后缀，并确认名称唯一，使两种模态能按同一细胞标识对齐。

rna_cite.obs_names
Index(['GCATTAGCATAAGCGG-1-s1d1', 'TACAGGTGTTAGAGTA-1-s1d1', 'AGGATCTAGGTCTACT-1-s1d1', 'GTAGAAAGTGACACAG-1-s1d1', 'TCCGAAAAGGATCATA-1-s1d1', 'CTCCCAATCCATTGGA-1-s1d1', 'GACCAATCAATTTCGG-1-s1d1', 'TTCCGGTAGTTGTAAG-1-s1d1', 'ACCTGTCAGGACTGGT-1-s1d1', 'TTCGATTTCAGGACAG-1-s1d1', ... 'TCTTCCTAGCCAACCC-1-s4d9', 'TTCCACGGTTGAGGAC-1-s4d9', 'ATTCCTAGTCCAAGAG-1-s4d9', 'GCGGAAAGTACGCGTC-1-s4d9', 'TAACTTCAGATACAGT-1-s4d9', 'GAATCACCACGGAAGT-1-s4d9', 'GCTGGGTGTACGGATG-1-s4d9', 'TCGAAGTGTGACAGGT-1-s4d9', 'GCAGGCTGTTGCATAC-1-s4d9', 'ACGTAACAGGTCTACT-1-s4d9'], dtype='object', length=90261)
adt_cite.obs_names
Index(['AAACCCAAGGATGGCT-1-0-0-0-0-0-0-0-0-0-0-0', 'AAACCCAAGGCCTAGA-1-0-0-0-0-0-0-0-0-0-0-0', 'AAACCCAAGTGAGTGC-1-0-0-0-0-0-0-0-0-0-0-0', 'AAACCCACAAGAGGCT-1-0-0-0-0-0-0-0-0-0-0-0', 'AAACCCACATCGTGGC-1-0-0-0-0-0-0-0-0-0-0-0', 'AAACCCACATTCTCTA-1-0-0-0-0-0-0-0-0-0-0-0', 'AAACCCAGTCCGCAGT-1-0-0-0-0-0-0-0-0-0-0-0', 'AAACCCAGTGCATACT-1-0-0-0-0-0-0-0-0-0-0-0', 'AAACCCAGTTGACGGA-1-0-0-0-0-0-0-0-0-0-0-0', 'AAACCCATCGATACTG-1-0-0-0-0-0-0-0-0-0-0-0', ... 'TTTGGTTGTCGTACTA-1-1', 'TTTGGTTTCACCTCGT-1-1', 'TTTGGTTTCGAGAACG-1-1', 'TTTGGTTTCGATACGT-1-1', 'TTTGTTGGTAGCTTTG-1-1', 'TTTGTTGGTCCAATCA-1-1', 'TTTGTTGGTGACTAAA-1-1', 'TTTGTTGTCCCTCTCC-1-1', 'TTTGTTGTCTAGAGCT-1-1', 'TTTGTTGTCTCTGCCA-1-1'], dtype='object', length=120502)
adt_cite.obs_names = [
    name.split("-")[0] + "-" + name.split("-")[1] + "-" + batch
    for batch, name in zip(adt_cite.obs["donor"], adt_cite.obs_names, strict=False)
]

用同一份共有细胞名称列表索引 RNA 与 ADT，既保留交集，也统一行序。

common_idx = list(set(rna_cite.obs_names).intersection(set(adt_cite.obs_names)))
rna_cite = rna_cite[common_idx].copy()
adt_cite = adt_cite[common_idx].copy()

按细胞索引，将 RNA 中的 Samplename 和 cell_type 注释加入 ADT。

adt_cite.obs = adt_cite.obs.join(rna_cite.obs[["Samplename", "cell_type"]])

对 Multiome 也检查 RNA 与 ATAC 的细胞名称及顺序完全一致。下方断言用于验证，不会自动重排数据。

assert np.sum(rna_multiome.obs_names != atac_multiome.obs_names) == 0

为区分不同数据集与模态中重复的批次名称，下面定义辅助函数，在 .obs 中添加带前缀或后缀的标签列。默认新列名为 new_ 加原列名。注意第四个位置参数是 how：ADT 调用中的 "batch" 会选择前缀分支，生成 new_donor，而不是把新列命名为 batch。实际使用时需按函数签名核对。

def update_obs_column(
    adata, obs_column_name, suffix, how="right", new_obs_column_name=None
):
    if new_obs_column_name is None:
        new_obs_column_name = f"new_{obs_column_name}"
    # otherwise can't use + operator with categorical columns
    adata.obs[obs_column_name] = adata.obs[obs_column_name].astype("str").copy()
    # create new column in .obs
    if how == "right":
        adata.obs[new_obs_column_name] = f"_{suffix}"
        adata.obs[new_obs_column_name] = (
            adata.obs[obs_column_name] + adata.obs[new_obs_column_name]
        )
    else:
        adata.obs[new_obs_column_name] = f"{suffix}_"
        adata.obs[new_obs_column_name] = (
            adata.obs[new_obs_column_name] + adata.obs[obs_column_name]
        )
    return adata
rna_multiome = update_obs_column(rna_multiome, "batch", "_rna_multiome")
atac_multiome = update_obs_column(atac_multiome, "batch", "_atac_multiome")
rna_cite = update_obs_column(rna_cite, "batch", "_rna_cite")
adt_cite = update_obs_column(adt_cite, "donor", "_adt_cite", "batch")

根据前面列出的批次，分别提取参考集和查询集。

# query
rna_multiome_query = rna_multiome[
    rna_multiome.obs["batch"].isin(multiome_query_batches)
].copy()
atac_multiome_query = atac_multiome[
    atac_multiome.obs["batch"].isin(multiome_query_batches)
].copy()
rna_cite_query = rna_cite[rna_cite.obs["batch"].isin(cite_query_batches)].copy()
adt_cite_query = adt_cite[adt_cite.obs["donor"].isin(cite_query_batches)].copy()
# reference
rna_multiome = rna_multiome[
    rna_multiome.obs["batch"].isin(multiome_reference_batches)
].copy()
atac_multiome = atac_multiome[
    atac_multiome.obs["batch"].isin(multiome_reference_batches)
].copy()
rna_cite = rna_cite[rna_cite.obs["batch"].isin(cite_reference_batches)].copy()
adt_cite = adt_cite[adt_cite.obs["donor"].isin(cite_reference_batches)].copy()

使用 GLUE 进行未配对整合

未配对整合的关键是没有已知的逐细胞对应关系，而不一定是所有特征（feature）都没有交集。本例 RNA 与 ADT 来自不同细胞，测量的分子类型也不同，因此用编码基因与蛋白之间的先验关系联系两种模态。

GLUE Cao & Gao, 2022 是一种未配对整合模型：各模态使用变分自编码器（Variational Autoencoder, VAE）学习细胞表示，并以批次等信息为条件；先验图连接不同模态的 feature，为 feature 的嵌入（Embedding）提供约束，再结合跨模态对齐学习共享的细胞空间。本例使用 NeurIPS 竞赛的 Multiome RNA 与 CITE-seq ADT 数据（https://openproblems.bio/neurips_2021/），按预先指定的编码基因–蛋白对应关系构图。模型输出每个细胞在共享潜在空间（latent space）中的表示；这种统计对齐不意味着恢复了真实的一一配对细胞。

未配对 RNA 与 ATAC 的整合示例及模型细节见 GLUE 教程：https://scglue.readthedocs.io/en/latest/tutorials.html。

按照 GLUE 教程准备 RNA 编码器输入：从 counts 层取回原始计数（Count），依次按细胞总量归一化（normalization）、log1p 和标准化，再计算 100 维主成分分析（Principal Component Analysis, PCA）表示。原始 Count 层保留给后续重建目标使用。

rna_multiome.X = rna_multiome.layers["counts"].copy()
sc.pp.normalize_total(rna_multiome)
sc.pp.log1p(rna_multiome)
sc.pp.scale(rna_multiome)
sc.tl.pca(rna_multiome, n_comps=100, svd_solver="auto")

ADT 编码器使用 .X 中预先完成中心化对数比转换（Centered Log-Ratio Transformation, CLR）的数据，再计算 PCA。查看最大值只能辅助检查尺度，不能单独证明数据已经按预期归一化。

np.max(adt_cite.X)
69.27085
sc.tl.pca(adt_cite, n_comps=100, svd_solver="auto")

先把蛋白名称映射到对应的基因符号，以便构图。下方字典是简化映射：CD3、TCR 或 HLA-DR 等复合体/标记不一定与单个基因一一对应，实际分析应结合抗体靶标核对先验关系；重命名也不代表 RNA 与蛋白信号等价。

rename_proteins = {
    "CD103": "ITGAE",
    "CD11b": "ITGAM",
    "CD11c": "ITGAX",
    "CD122": "IL2RB",
    "CD124": "IL4R",
    "CD127": "IL7R",
    "CD13": "ANPEP",
    "CD137": "TNFRSF9",
    "CD152": "CTLA4",
    "CD154": "CD40LG",
    "CD16": "FCGR3A",
    "CD161": "KLRB1",
    "CD185": "CXCR5",
    "CD194": "CCR4",
    "CD196": "CCR6",
    "CD20": "MS4A1",
    "CD21": "CR2",
    "CD23": "FCER2",
    "CD25": "IL2RA",
    "CD26": "DPP4",
    "CD268": "TNFRSF13C",
    "CD272": "BTLA",
    "CD278": "ICOS",
    "CD29": "ITGB1",
    "CD3": "CD3G",
    "CD303": "CLEC4C",
    "CD304": "NRP1",
    "CD314": "KLRK1",
    "CD319": "SLAMF7",
    "CD335": "NCR1",
    "CD35": "CR1",
    "CD352": "SLAMF6",
    "CD39": "ENTPD1",
    "CD49a": "ITGA1",
    "CD49f": "ITGA6",
    "CD54": "ICAM1",
    "CD56": "NCAM1",
    "CD62L": "SELL",
    "CD71": "TFRC",
    "CD73": "NT5E",
    "CD79b": "CD79B",
    "CD8": "CD8A",
    "CD85j": "LILRB1",
    "CD88": "C5AR1",
    "CD94": "KLRD1",
    "CD95": "FAS",
    "HLA-DR": "HLA-DRA",
    "IgD": "IGHD",
    "IgM": "IGHM",
    "TCR": "TRAC",
}
adt_cite.var_names = [
    rename_proteins[name] if name in rename_proteins else name
    for name in adt_cite.var_names
]

构建先验图

图的节点覆盖两种模态的全部 feature，先验关系的权重在 0 到 1 之间，零表示不连接。本例先比较重命名后的蛋白与基因名称，以布尔矩阵标记对应关系；只有匹配项进入图，边权为 1，其余项不添加边。GLUE 可使用带权、带符号的先验关系，此处仅演示正向的基因–蛋白连接。

p = np.array(adt_cite.var_names)
r = np.array(rna_multiome.var_names)  # noqa F811
# mask entries are set to 1 where protein name is the same as gene name
mask = np.repeat(p.reshape(-1, 1), r.shape[0], axis=1) == r
mask = np.array(mask)

匹配关系确定后，给基因和蛋白名称分别添加 _rna 与 _prot 后缀，确保它们是图中不同的节点。

rna_vars = [v + "_rna" for v in rna_multiome.var_names]
prot_vars = [v + "_prot" for v in adt_cite.var_names]
rna_multiome.var_names = rna_vars
adt_cite.var_names = prot_vars

再为每个节点添加权重为 1 的自环（self-loop）。本例所有边的 sign 均为 1。

adj = pd.DataFrame(mask, index=prot_vars, columns=rna_vars)
diag_edges = adj[adj > 0].stack().index.tolist()
diag_edges = [(n1, n2, {"weight": 1.0, "sign": 1}) for n1, n2 in diag_edges]
self_loop_rna = [(g, g, {"weight": 1.0, "sign": 1}) for g in rna_vars]
self_loop_prot = [(g, g, {"weight": 1.0, "sign": 1}) for g in prot_vars]

使用 NetworkX 创建无向图，加入两种模态的节点、跨模态连接和自环。

graph = nx.Graph()
graph.add_nodes_from(rna_vars)
graph.add_nodes_from(prot_vars)
graph.add_edges_from(diag_edges)
graph.add_edges_from(self_loop_prot)
graph.add_edges_from(self_loop_rna)

原教程中，该图有 4,136 个节点：4,000 个基因和 136 个蛋白；共有 4,186 条边，其中 4,136 条是自环，50 条连接匹配的基因与蛋白。这些数量取决于实际输入及名称映射，不能看作 GLUE 的固定要求。

graph.number_of_nodes(), graph.number_of_edges()
(4136, 4186)

配置数据

配置 RNA 的编码器与解码器：编码器读取 X_pca 中的低维表示，解码器以 counts 层中的原始 Count 为重建目标，并采用负二项分布（Negative Binomial Distribution, NB）。Samplename 用于指定批次；use_highly_variable=False 表示使用当前对象中的全部 feature。

scglue.models.configure_dataset(
    rna_multiome,
    "NB",
    use_highly_variable=False,
    use_batch="Samplename",
    use_layer="counts",
    use_rep="X_pca",
)

ADT 同样以 PCA 表示作为编码器输入，但解码器重建的是 .X 中的归一化值，采用高斯分布（Gaussian distribution）。这里没有把 CLR 值当作原始 Count 使用。

scglue.models.configure_dataset(
    adt_cite,
    "Normal",
    use_highly_variable=False,
    use_batch="Samplename",
    use_rep="X_pca",
)

传入两种模态及先验图，初始化并训练 GLUE。

glue = scglue.models.fit_SCGLUE(
    {"rna": rna_multiome, "adt": adt_cite},
    graph,
)
输出
[INFO] fit_SCGLUE: Pretraining SCGLUE model...
[INFO] autodevice: Using GPU 0 as computation device.
[INFO] check_graph: Checking variable coverage...
[INFO] check_graph: Checking edge attributes...
[INFO] check_graph: Checking self-loops...
[INFO] check_graph: Checking graph symmetry...
[INFO] check_graph: All checks passed!
[INFO] SCGLUEModel: Setting `graph_batch_size` = 1457
[INFO] SCGLUEModel: Setting `max_epochs` = 102
[INFO] SCGLUEModel: Setting `patience` = 9
[INFO] SCGLUEModel: Setting `reduce_lr_patience` = 5
[INFO] SCGLUETrainer: Using training directory: "/tmp/GLUETMPx6ds9f5b"
[INFO] SCGLUETrainer: [Epoch 10] train={'g_nll': 0.199, 'g_kl': 0.039, 'g_elbo': 0.238, 'x_rna_nll': 0.241, 'x_rna_kl': 0.006, 'x_rna_elbo': 0.247, 'x_adt_nll': 1.934, 'x_adt_kl': 0.21, 'x_adt_elbo': 2.144, 'dsc_loss': 0.636, 'vae_loss': 2.401, 'gen_loss': 2.369}, val={'g_nll': 0.193, 'g_kl': 0.04, 'g_elbo': 0.233, 'x_rna_nll': 0.243, 'x_rna_kl': 0.006, 'x_rna_elbo': 0.249, 'x_adt_nll': 1.862, 'x_adt_kl': 0.199, 'x_adt_elbo': 2.061, 'dsc_loss': 0.63, 'vae_loss': 2.32, 'gen_loss': 2.288}, 4.8s elapsed
[INFO] SCGLUETrainer: [Epoch 20] train={'g_nll': 0.196, 'g_kl': 0.041, 'g_elbo': 0.237, 'x_rna_nll': 0.237, 'x_rna_kl': 0.006, 'x_rna_elbo': 0.242, 'x_adt_nll': 1.845, 'x_adt_kl': 0.175, 'x_adt_elbo': 2.02, 'dsc_loss': 0.674, 'vae_loss': 2.272, 'gen_loss': 2.238}, val={'g_nll': 0.194, 'g_kl': 0.041, 'g_elbo': 0.235, 'x_rna_nll': 0.236, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.242, 'x_adt_nll': 1.798, 'x_adt_kl': 0.173, 'x_adt_elbo': 1.97, 'dsc_loss': 0.67, 'vae_loss': 2.222, 'gen_loss': 2.188}, 5.9s elapsed
[INFO] SCGLUETrainer: [Epoch 30] train={'g_nll': 0.182, 'g_kl': 0.041, 'g_elbo': 0.224, 'x_rna_nll': 0.236, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.241, 'x_adt_nll': 1.835, 'x_adt_kl': 0.172, 'x_adt_elbo': 2.007, 'dsc_loss': 0.676, 'vae_loss': 2.257, 'gen_loss': 2.223}, val={'g_nll': 0.155, 'g_kl': 0.041, 'g_elbo': 0.196, 'x_rna_nll': 0.237, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.242, 'x_adt_nll': 1.782, 'x_adt_kl': 0.166, 'x_adt_elbo': 1.948, 'dsc_loss': 0.669, 'vae_loss': 2.198, 'gen_loss': 2.165}, 12.0s elapsed
[INFO] SCGLUETrainer: [Epoch 40] train={'g_nll': 0.147, 'g_kl': 0.042, 'g_elbo': 0.189, 'x_rna_nll': 0.236, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.241, 'x_adt_nll': 1.833, 'x_adt_kl': 0.169, 'x_adt_elbo': 2.002, 'dsc_loss': 0.677, 'vae_loss': 2.251, 'gen_loss': 2.217}, val={'g_nll': 0.159, 'g_kl': 0.042, 'g_elbo': 0.2, 'x_rna_nll': 0.238, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.243, 'x_adt_nll': 1.781, 'x_adt_kl': 0.166, 'x_adt_elbo': 1.947, 'dsc_loss': 0.674, 'vae_loss': 2.199, 'gen_loss': 2.165}, 12.1s elapsed
Epoch 00044: reducing learning rate of group 0 to 2.0000e-04.
Epoch 00044: reducing learning rate of group 0 to 2.0000e-04.
[INFO] LRScheduler: Learning rate reduction: step 1
[INFO] SCGLUETrainer: [Epoch 50] train={'g_nll': 0.146, 'g_kl': 0.042, 'g_elbo': 0.188, 'x_rna_nll': 0.235, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.24, 'x_adt_nll': 1.824, 'x_adt_kl': 0.17, 'x_adt_elbo': 1.993, 'dsc_loss': 0.676, 'vae_loss': 2.241, 'gen_loss': 2.207}, val={'g_nll': 0.142, 'g_kl': 0.042, 'g_elbo': 0.183, 'x_rna_nll': 0.235, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.241, 'x_adt_nll': 1.771, 'x_adt_kl': 0.166, 'x_adt_elbo': 1.938, 'dsc_loss': 0.669, 'vae_loss': 2.185, 'gen_loss': 2.152}, 14.4s elapsed
Epoch 00058: reducing learning rate of group 0 to 2.0000e-05.
Epoch 00058: reducing learning rate of group 0 to 2.0000e-05.
[INFO] LRScheduler: Learning rate reduction: step 2
[INFO] SCGLUETrainer: [Epoch 60] train={'g_nll': 0.147, 'g_kl': 0.042, 'g_elbo': 0.188, 'x_rna_nll': 0.235, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.24, 'x_adt_nll': 1.818, 'x_adt_kl': 0.17, 'x_adt_elbo': 1.989, 'dsc_loss': 0.677, 'vae_loss': 2.237, 'gen_loss': 2.203}, val={'g_nll': 0.159, 'g_kl': 0.042, 'g_elbo': 0.201, 'x_rna_nll': 0.239, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.244, 'x_adt_nll': 1.77, 'x_adt_kl': 0.167, 'x_adt_elbo': 1.936, 'dsc_loss': 0.668, 'vae_loss': 2.188, 'gen_loss': 2.155}, 6.4s elapsed

2023-01-04 16:22:27,927 ignite.handlers.early_stopping.EarlyStopping INFO: EarlyStopping: Stop training

[INFO] EarlyStopping: Restoring checkpoint "59"...
[INFO] EarlyStopping: Restoring checkpoint "59"...
[INFO] fit_SCGLUE: Estimating balancing weight...
[INFO] estimate_balancing_weight: Clustering cells...
[INFO] estimate_balancing_weight: Matching clusters...
[INFO] estimate_balancing_weight: Matching array shape = (24, 25)...
[INFO] estimate_balancing_weight: Estimating balancing weight...
[INFO] fit_SCGLUE: Fine-tuning SCGLUE model...
[INFO] check_graph: Checking variable coverage...
[INFO] check_graph: Checking edge attributes...
[INFO] check_graph: Checking self-loops...
[INFO] check_graph: Checking graph symmetry...
[INFO] check_graph: All checks passed!
[INFO] SCGLUEModel: Setting `graph_batch_size` = 1457
[INFO] SCGLUEModel: Setting `align_burnin` = 17
[INFO] SCGLUEModel: Setting `max_epochs` = 102
[INFO] SCGLUEModel: Setting `patience` = 9
[INFO] SCGLUEModel: Setting `reduce_lr_patience` = 5
[INFO] SCGLUETrainer: Using training directory: "/tmp/GLUETMPe8ok569r"
[INFO] SCGLUETrainer: [Epoch 10] train={'g_nll': 0.136, 'g_kl': 0.042, 'g_elbo': 0.178, 'x_rna_nll': 0.237, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.242, 'x_adt_nll': 1.831, 'x_adt_kl': 0.168, 'x_adt_elbo': 2.0, 'dsc_loss': 0.662, 'vae_loss': 2.249, 'gen_loss': 2.216}, val={'g_nll': 0.131, 'g_kl': 0.042, 'g_elbo': 0.173, 'x_rna_nll': 0.237, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.243, 'x_adt_nll': 1.787, 'x_adt_kl': 0.166, 'x_adt_elbo': 1.952, 'dsc_loss': 0.659, 'vae_loss': 2.202, 'gen_loss': 2.169}, 8.9s elapsed
[INFO] SCGLUETrainer: [Epoch 20] train={'g_nll': 0.126, 'g_kl': 0.042, 'g_elbo': 0.167, 'x_rna_nll': 0.238, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.244, 'x_adt_nll': 1.829, 'x_adt_kl': 0.168, 'x_adt_elbo': 1.997, 'dsc_loss': 0.641, 'vae_loss': 2.247, 'gen_loss': 2.215}, val={'g_nll': 0.135, 'g_kl': 0.042, 'g_elbo': 0.177, 'x_rna_nll': 0.239, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.244, 'x_adt_nll': 1.786, 'x_adt_kl': 0.164, 'x_adt_elbo': 1.95, 'dsc_loss': 0.635, 'vae_loss': 2.201, 'gen_loss': 2.169}, 5.9s elapsed
[INFO] SCGLUETrainer: [Epoch 30] train={'g_nll': 0.126, 'g_kl': 0.041, 'g_elbo': 0.168, 'x_rna_nll': 0.238, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.243, 'x_adt_nll': 1.828, 'x_adt_kl': 0.168, 'x_adt_elbo': 1.995, 'dsc_loss': 0.647, 'vae_loss': 2.245, 'gen_loss': 2.213}, val={'g_nll': 0.129, 'g_kl': 0.042, 'g_elbo': 0.171, 'x_rna_nll': 0.239, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.244, 'x_adt_nll': 1.779, 'x_adt_kl': 0.164, 'x_adt_elbo': 1.943, 'dsc_loss': 0.634, 'vae_loss': 2.194, 'gen_loss': 2.163}, 4.4s elapsed
Epoch 00036: reducing learning rate of group 0 to 2.0000e-04.
Epoch 00036: reducing learning rate of group 0 to 2.0000e-04.
[INFO] LRScheduler: Learning rate reduction: step 1
[INFO] SCGLUETrainer: [Epoch 40] train={'g_nll': 0.114, 'g_kl': 0.041, 'g_elbo': 0.155, 'x_rna_nll': 0.238, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.243, 'x_adt_nll': 1.815, 'x_adt_kl': 0.168, 'x_adt_elbo': 1.982, 'dsc_loss': 0.667, 'vae_loss': 2.232, 'gen_loss': 2.198}, val={'g_nll': 0.132, 'g_kl': 0.041, 'g_elbo': 0.173, 'x_rna_nll': 0.239, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.244, 'x_adt_nll': 1.777, 'x_adt_kl': 0.165, 'x_adt_elbo': 1.942, 'dsc_loss': 0.632, 'vae_loss': 2.192, 'gen_loss': 2.161}, 10.3s elapsed
[INFO] SCGLUETrainer: [Epoch 50] train={'g_nll': 0.117, 'g_kl': 0.041, 'g_elbo': 0.158, 'x_rna_nll': 0.238, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.243, 'x_adt_nll': 1.814, 'x_adt_kl': 0.169, 'x_adt_elbo': 1.983, 'dsc_loss': 0.671, 'vae_loss': 2.232, 'gen_loss': 2.198}, val={'g_nll': 0.115, 'g_kl': 0.041, 'g_elbo': 0.156, 'x_rna_nll': 0.24, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.245, 'x_adt_nll': 1.77, 'x_adt_kl': 0.166, 'x_adt_elbo': 1.936, 'dsc_loss': 0.632, 'vae_loss': 2.188, 'gen_loss': 2.156}, 9.2s elapsed
Epoch 00050: reducing learning rate of group 0 to 2.0000e-05.
Epoch 00050: reducing learning rate of group 0 to 2.0000e-05.
[INFO] LRScheduler: Learning rate reduction: step 2
Epoch 00058: reducing learning rate of group 0 to 2.0000e-06.
Epoch 00058: reducing learning rate of group 0 to 2.0000e-06.
[INFO] LRScheduler: Learning rate reduction: step 3
[INFO] SCGLUETrainer: [Epoch 60] train={'g_nll': 0.118, 'g_kl': 0.041, 'g_elbo': 0.159, 'x_rna_nll': 0.237, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.242, 'x_adt_nll': 1.815, 'x_adt_kl': 0.169, 'x_adt_elbo': 1.983, 'dsc_loss': 0.668, 'vae_loss': 2.232, 'gen_loss': 2.199}, val={'g_nll': 0.113, 'g_kl': 0.041, 'g_elbo': 0.154, 'x_rna_nll': 0.239, 'x_rna_kl': 0.005, 'x_rna_elbo': 0.244, 'x_adt_nll': 1.771, 'x_adt_kl': 0.166, 'x_adt_elbo': 1.937, 'dsc_loss': 0.643, 'vae_loss': 2.187, 'gen_loss': 2.155}, 14.5s elapsed
Epoch 00065: reducing learning rate of group 0 to 2.0000e-07.
Epoch 00065: reducing learning rate of group 0 to 2.0000e-07.
[INFO] LRScheduler: Learning rate reduction: step 4

2023-01-04 16:32:13,598 ignite.handlers.early_stopping.EarlyStopping INFO: EarlyStopping: Stop training

[INFO] EarlyStopping: Restoring checkpoint "59"...
[INFO] EarlyStopping: Restoring checkpoint "59"...

分别编码 RNA 与 ADT 细胞，得到同维度的潜在表示，再按细胞拼接到一个 AnnData 中供可视化使用。下方拼接不是合并配对测量；后续图和二维可视化将使用 X_glue，而不是直接比较 RNA 与蛋白的原始矩阵。

rna_multiome.obsm["X_glue"] = glue.encode_data("rna", rna_multiome)
adt_cite.obsm["X_glue"] = glue.encode_data("adt", adt_cite)
adt_cite.obs["modality"] = "ADT"
rna_multiome.obs["modality"] = "RNA"

combined = anndata.concat([rna_multiome, adt_cite])

在共享表示中以余弦距离构建近邻图（nearest-neighbor graph），再计算统一流形近似与投影（Uniform Manifold Approximation and Projection, UMAP）。

sc.pp.neighbors(combined, use_rep="X_glue", metric="cosine")
sc.tl.umap(combined)
sc.pl.umap(
    combined, color=["cell_type", "modality", "Samplename"], ncols=1, frameon=False
)
<Figure size 727.8x1440 with 3 Axes>

本例可视化中，各批次和模态有一定重叠。是否达到合适整合，还需同时评估批次混合与生物学信息保留，不能只凭 UMAP 判断。可参考单细胞整合基准评估（Single-Cell Integration Benchmarking, scIB）的相关指标，但这些指标并未全面覆盖多模态整合质量。

使用 totalVI 整合部分重叠数据并映射查询

另一种情形是同时拥有配对与未配对数据，例如 CITE-seq 与仅测 RNA 的数据集。可以利用共有 RNA 模态将它们接入同一潜在空间，也可以由模型估计缺失的蛋白信号。本节演示部分重叠数据的整合与查询映射。原文还提到蛋白插补（imputation），但下面没有展示提取插补蛋白值的调用；映射得到潜在表示本身不等于已经输出蛋白预测。本例使用的方法为 totalVI Gayoso et al., 2021。

将参考中 s1d3 批次的全部蛋白 Count 置零，以模拟未测蛋白。totalVI 默认把某批次中全为零的蛋白视为缺失测量；这里的零是缺失标记，不是已测得蛋白丰度为零。该判定也可能与真实生物学零值混淆，使用自己的数据时需核对缺失设置。

adata = rna_cite.copy()
adata.obsm["protein_counts"] = adt_cite.layers["counts"].A.copy()
adata.obsm["protein_counts"][adata.obs["batch"] == "s1d3"] = 0.0
adata
AnnData object with n_obs × n_vars = 16294 × 4000 obs: 'GEX_n_genes_by_counts', 'GEX_pct_counts_mt', 'GEX_size_factors', 'GEX_phase', 'ADT_n_antibodies_by_counts', 'ADT_total_counts', 'ADT_iso_count', 'cell_type', 'batch', 'ADT_pseudotime_order', 'GEX_pseudotime_order', 'Samplename', 'Site', 'DonorNumber', 'Modality', 'VendorLot', 'DonorID', 'DonorAge', 'DonorBMI', 'DonorBloodType', 'DonorRace', 'Ethnicity', 'DonorGender', 'QCMeds', 'DonorSmoker', 'is_train', 'new_batch' var: 'feature_types', 'gene_id', 'highly_variable', 'means', 'dispersions', 'dispersions_norm' uns: 'dataset_id', 'genome', 'hvg', 'organism' obsm: 'ADT_X_pca', 'ADT_X_umap', 'ADT_isotype_controls', 'GEX_X_pca', 'GEX_X_umap', 'protein_counts' layers: 'counts'

注册 RNA 原始 Count 层、蛋白矩阵与批次标签，再按示例设置编码器和解码器的归一化方式及层数，训练参考模型。

TOTALVI.setup_anndata(
    adata,
    layer="counts",
    batch_key="batch",
    protein_expression_obsm_key="protein_counts",
)
No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)

INFO     Generating sequential column names                                                                        
INFO     Found batches with missing protein expression                                                             

arches_params = {
    "use_layer_norm": "both",
    "use_batch_norm": "none",
    "n_layers_decoder": 2,
    "n_layers_encoder": 2,
}

vae = TOTALVI(adata, **arches_params)
vae.train()
输出
INFO     Computing empirical prior initialization for protein background.                                          

GPU available: True (cuda), used: True
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
HPU available: False, using: 0 HPUs
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]

Epoch 258/400:  64%|██████▍   | 258/400 [06:23<03:23,  1.43s/it, loss=709, v_num=1]Epoch 00258: reducing learning rate of group 0 to 2.4000e-03.
Epoch 381/400:  95%|█████████▌| 381/400 [09:23<00:27,  1.44s/it, loss=703, v_num=1]Epoch 00381: reducing learning rate of group 0 to 1.4400e-03.
Epoch 400/400: 100%|██████████| 400/400 [09:51<00:00,  1.44s/it, loss=710, v_num=1]
`Trainer.fit` stopped: `max_epochs=400` reached.

Epoch 400/400: 100%|██████████| 400/400 [09:51<00:00,  1.48s/it, loss=710, v_num=1]

提取参考的潜在表示，并计算邻居图与 UMAP。

adata.obsm["X_totalvi"] = vae.get_latent_representation()
sc.pp.neighbors(adata, use_rep="X_totalvi")
sc.tl.umap(adata)
sc.pl.umap(adata, color=["cell_type", "batch"], ncols=1, frameon=False)
<Figure size 727.8x960 with 2 Axes>

图中配对数据与仅 RNA 数据在相关细胞类型附近重叠。还需检查局部邻域、批次结构和生物学差异，才能评价整合质量。

查询到参考映射

将新的仅 RNA 查询与 CITE-seq 查询接入上述 totalVI 参考。

把查询中 s2d4 批次的全部蛋白 Count 置零，模拟仅 RNA 的查询；另一个批次保留配对测量。

query = rna_cite_query.copy()
query.obsm["protein_counts"] = adt_cite_query.layers["counts"].A.copy()
query.obsm["protein_counts"][query.obs["batch"] == "s2d4"] = 0.0

在 .obs 中加入参考/查询及具体数据类型的标签，供后续绘图使用。

adata.obs["dataset_name"] = "Reference"
query.obs["dataset_name"] = "Query"
adata.obs["dataset_name_fine"] = "CITE reference"
adata.obs["dataset_name_fine"][adata.obs["batch"] == "s1d3"] = "RNA reference"
query.obs["dataset_name_fine"] = "CITE query"
query.obs["dataset_name_fine"][query.obs["batch"] == "s2d4"] = "RNA query"

用 load_query_data 从参考模型初始化查询模型，并按其查询适配机制训练；本例将 weight_decay 与 scale_adversarial_loss 都设为 0。

vae_q = TOTALVI.load_query_data(
    query,
    vae,
)
vae_q.train(
    plan_kwargs={"weight_decay": 0.0, "scale_adversarial_loss": 0.0},
)
输出
INFO     Found batches with missing protein expression                                                             
INFO     Computing empirical prior initialization for protein background.                                          

GPU available: True (cuda), used: True
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
HPU available: False, using: 0 HPUs
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]

Epoch 233/400:  58%|█████▊    | 233/400 [05:32<03:54,  1.40s/it, loss=546, v_num=1]Epoch 00233: reducing learning rate of group 0 to 2.4000e-03.
Epoch 273/400:  68%|██████▊   | 273/400 [06:29<02:58,  1.40s/it, loss=548, v_num=1]Epoch 00273: reducing learning rate of group 0 to 1.4400e-03.
Epoch 287/400:  72%|███████▏  | 287/400 [06:49<02:41,  1.43s/it, loss=542, v_num=1]
Monitored metric elbo_validation did not improve in the last 45 records. Best score: 1102.225. Signaling Trainer to stop.

提取查询的潜在表示，与保留的参考表示合并，再重新计算邻居图和 UMAP。这里保留的是参考潜在表示；由于重新拟合了 UMAP，参考的二维坐标不会因此自动保持原样。

query.obsm["X_totalvi_scarches"] = vae_q.get_latent_representation(query)
adata.obsm["X_totalvi_scarches"] = adata.obsm["X_totalvi"]
full_data = adata.concatenate(query, batch_key="concat_batch")
full_data
AnnData object with n_obs × n_vars = 32340 × 4000 obs: 'GEX_n_genes_by_counts', 'GEX_pct_counts_mt', 'GEX_size_factors', 'GEX_phase', 'ADT_n_antibodies_by_counts', 'ADT_total_counts', 'ADT_iso_count', 'cell_type', 'batch', 'ADT_pseudotime_order', 'GEX_pseudotime_order', 'Samplename', 'Site', 'DonorNumber', 'Modality', 'VendorLot', 'DonorID', 'DonorAge', 'DonorBMI', 'DonorBloodType', 'DonorRace', 'Ethnicity', 'DonorGender', 'QCMeds', 'DonorSmoker', 'is_train', 'new_batch', '_scvi_labels', '_scvi_batch', 'dataset_name', 'dataset_name_fine', 'concat_batch' var: 'feature_types', 'gene_id', 'highly_variable', 'means', 'dispersions', 'dispersions_norm' obsm: 'ADT_X_pca', 'ADT_X_umap', 'ADT_isotype_controls', 'GEX_X_pca', 'GEX_X_umap', 'protein_counts', 'X_totalvi_scarches' layers: 'counts'
sc.pp.neighbors(full_data, use_rep="X_totalvi_scarches")
sc.tl.umap(full_data)
sc.pl.umap(
    full_data,
    color=["cell_type", "batch", "dataset_name", "dataset_name_fine"],
    ncols=1,
    frameon=False,
)
<Figure size 727.8x1920 with 4 Axes>

使用 Seurat WNN 进行查询映射与插补

Seurat 自 v4 起支持将新的 RNA 查询映射到加权最近邻（Weighted Nearest Neighbor, WNN）多模态参考 Hao et al., 2021。这里使用前一章配对整合得到的参考对象。

该方法在查询与参考之间寻找锚点（anchor），据此转移细胞类型标签，并估计查询中缺失的蛋白信号。结果依赖参考覆盖范围与匹配质量；预测标签和插补值都需要验证。

Seurat 是 R 软件包，因此先借助 anndata2ri 软件包（https://github.com/theislab/anndata2ri）将 Python AnnData 转换为 R 对象。

先读取前一章保存的参考并查看 RNA 查询。原文的 readRDS 调用额外传入 cite；读取该文件通常只需 readRDS(file="wnn_ref.rds")，执行前应核对 R 参数匹配。

%%R
ref <- readRDS(cite, file = "wnn_ref.rds")
ref
An object of class Seurat 
4136 features across 16294 samples within 2 assays 
Active assay: RNA (4000 features, 4000 variable features)
 1 other assay present: ADT
 4 dimensional reductions calculated: pca, harmony_pca, spca, wnn.umap

rna_cite_query
AnnData object with n_obs × n_vars = 16046 × 4000 obs: 'GEX_n_genes_by_counts', 'GEX_pct_counts_mt', 'GEX_size_factors', 'GEX_phase', 'ADT_n_antibodies_by_counts', 'ADT_total_counts', 'ADT_iso_count', 'cell_type', 'batch', 'ADT_pseudotime_order', 'GEX_pseudotime_order', 'Samplename', 'Site', 'DonorNumber', 'Modality', 'VendorLot', 'DonorID', 'DonorAge', 'DonorBMI', 'DonorBloodType', 'DonorRace', 'Ethnicity', 'DonorGender', 'QCMeds', 'DonorSmoker', 'is_train', 'new_batch' var: 'feature_types', 'gene_id', 'highly_variable', 'means', 'dispersions', 'dispersions_norm' uns: 'dataset_id', 'genome', 'hvg', 'organism' obsm: 'ADT_X_pca', 'ADT_X_umap', 'ADT_isotype_controls', 'GEX_X_pca', 'GEX_X_umap' layers: 'counts'

先对 RNA 查询当前 .X 中的数据计算 PCA。原文称其为 CLR 归一化 Count，但此处使用的是 RNA 数据，代码也没有执行 CLR；应确认文件中保存的 RNA 归一化结果与参考处理一致，不能套用前面 ADT 的尺度说明。

sc.pp.pca(rna_cite_query)
rna_cite_query
AnnData object with n_obs × n_vars = 16046 × 4000 obs: 'GEX_n_genes_by_counts', 'GEX_pct_counts_mt', 'GEX_size_factors', 'GEX_phase', 'ADT_n_antibodies_by_counts', 'ADT_total_counts', 'ADT_iso_count', 'cell_type', 'batch', 'ADT_pseudotime_order', 'GEX_pseudotime_order', 'Samplename', 'Site', 'DonorNumber', 'Modality', 'VendorLot', 'DonorID', 'DonorAge', 'DonorBMI', 'DonorBloodType', 'DonorRace', 'Ethnicity', 'DonorGender', 'QCMeds', 'DonorSmoker', 'is_train', 'new_batch' var: 'feature_types', 'gene_id', 'highly_variable', 'means', 'dispersions', 'dispersions_norm' uns: 'dataset_id', 'genome', 'hvg', 'organism', 'pca' obsm: 'ADT_X_pca', 'ADT_X_umap', 'ADT_isotype_controls', 'GEX_X_pca', 'GEX_X_umap', 'X_pca' varm: 'PCs' layers: 'counts'
np.max(rna_cite_query.X)
9.061238

复制查询所需的矩阵、细胞注释和 PCA 表示，转换为 R 中的 Seurat 对象。

adata_ = sc.AnnData(rna_cite_query.X.copy())
adata_.obs_names = rna_cite_query.obs_names.copy()
adata_.var_names = rna_cite_query.var_names.copy()
adata_.obs["cell_type"] = rna_cite_query.obs["cell_type"].copy()
adata_.obs["batch"] = rna_cite_query.obs["batch"].copy()
adata_.obsm["X_pca"] = rna_cite_query.obsm["X_pca"].copy()
%%R -i adata_
query = as.Seurat(adata_, data='X', counts=NULL)
query
An object of class Seurat 
4000 features across 16046 samples within 1 assay 
Active assay: originalexp (4000 features, 0 variable features)
 1 dimensional reduction calculated: PCA

寻找查询与参考之间的锚点，指定使用参考的有监督主成分分析（Supervised Principal Component Analysis, sPCA）结果及前 20 个维度。

%%R
anchors <- FindTransferAnchors(
  reference = ref,
  query = query,
  reference.reduction = "spca",
  dims = 1:20
)

调用 MapQuery，使用参考中保存的 wnn.umap 模型投影查询，以保持参考的二维坐标。通过 refdata 指定转移 cell_type 和 ADT 信息。预测的 ADT 继承参考测定（Assay）对象的数据尺度；前一章保存的是处理后的 ADT 数据，不能直接称为预测原始蛋白 Count。

%%R
query <- MapQuery(
  anchorset = anchors,
  query = query,
  reference = ref,
  refdata = list(
    cell_type = "cell_type",
    predicted_ADT = "ADT"
  ),
  reference.reduction = "spca",
  reduction.model = "wnn.umap",
  verbose=FALSE
)

合并参考与查询的 UMAP 坐标，比较预测标签和已有细胞类型注释。已有注释是比较基准，不一定是无误的生物学真值。

%%R
ref$id <- 'reference'
query$id <- 'query'
refquery <-  merge(ref, query)
refquery[["umap"]] <- merge(ref[["wnn.umap"]], query[["ref.umap"]])
%%R
p1 <- DimPlot(refquery, reduction = "umap", group.by = "cell_type", label = TRUE, label.size = 3, repel = TRUE) + NoLegend()
p2 <- DimPlot(refquery, reduction = "umap", group.by = "predicted.cell_type", label = TRUE, label.size = 3, repel = TRUE) + NoLegend()
p1 + p2
Fontconfig warning: ignoring UTF-8: not a valid region tag
Image produced in Jupyter
%%R
DimPlot(refquery, reduction = "umap", group.by = "id", label = TRUE, label.size = 3, repel = TRUE) + NoLegend()
Image produced in Jupyter

通过桥接整合将 ADT 映射到 RNA 图谱

桥接整合 Hao et al., 2022 利用包含配对 RNA 与目标模态的桥接数据，将仅有目标模态的查询映射到 RNA 参考。本例以配对 CITE-seq 作为桥，连接 RNA 参考与 ADT 查询。

为演示桥接，RNA 参考只取一个批次。若参考包含多个批次，可先按 Seurat 整合教程处理（https://satijalab.org/seurat/articles/integration_introduction.html）。本例将原 CITE-seq 查询中的一个批次作为配对桥接数据，另一个批次的 ADT 作为最终查询；三者的角色与前面的 totalVI 示例不同。

分别取 s1d1 的 RNA 作为参考、s2d1 的配对 RNA/ADT 作为桥、s2d4 的 ADT 作为查询。

rna_cite_ref = rna_cite[rna_cite.obs["batch"] == "s1d1"].copy()
rna_cite_bridge = rna_cite_query[rna_cite_query.obs["batch"] == "s2d1"].copy()
adt_cite_bridge = adt_cite_query[adt_cite_query.obs["donor"] == "s2d1"].copy()
adt_cite_query_bridge = adt_cite_query[adt_cite_query.obs["donor"] == "s2d4"].copy()

将参考 RNA 的原始 Count 转换到 R，使用 SCTransform 归一化，再计算 PCA 和 UMAP，并保存 UMAP 模型以投影后续查询。

adata_ = sc.AnnData(rna_cite_ref.layers["counts"].A.copy())
adata_.obs_names = rna_cite_ref.obs_names.copy()
adata_.var_names = rna_cite_ref.var_names.copy()
adata_.obs["cell_type"] = rna_cite_ref.obs["cell_type"].copy()
%%R -i adata_
ref = as.Seurat(adata_, data=NULL, counts='X')
ref
An object of class Seurat 
4000 features across 5219 samples within 1 assay 
Active assay: originalexp (4000 features, 0 variable features)

%%R
ref <- RenameAssays(object = ref, originalexp = "RNA", verbose=FALSE) 
ref <- SCTransform(ref, verbose=FALSE)
ref <- RunPCA(ref, verbose=FALSE)
ref <- RunUMAP(ref, dims = 1:50, return.model=TRUE, verbose=FALSE)

将桥接数据的 RNA 与 ADT 原始 Count 分别转换到 R，再合入一个对象。RNA 使用 SCTransform；ADT 进行 CLR 归一化和标准化。原文还提到此处计算 ADT PCA，但所列单元没有显式调用 RunPCA。另需核对 CreateAssayObject(counts=...) 所取的 adt@assays$originalexp@data 确实保存原始 Count，不能仅凭槽位名称假定其尺度。

adata_ = sc.AnnData(adt_cite_bridge.layers["counts"].A.copy())
adata_.obs_names = adt_cite_bridge.obs_names.copy()
adata_.var_names = adt_cite_bridge.var_names.copy()
adata_.obs["cell_type"] = adt_cite_bridge.obs["cell_type"].copy()
adata_
AnnData object with n_obs × n_vars = 10465 × 136 obs: 'cell_type'
%%R -i adata_
adt <- as.Seurat(adata_, data=NULL, counts='X')
adata_ = sc.AnnData(rna_cite_bridge.layers["counts"].A.copy())
adata_.obs_names = rna_cite_bridge.obs_names.copy()
adata_.var_names = rna_cite_bridge.var_names.copy()
adata_.obs["cell_type"] = rna_cite_bridge.obs["cell_type"].copy()
%%R -i adata_
rna = as.Seurat(adata_, data=NULL, counts='X')
cite <- rna
cite <- RenameAssays(object = cite, originalexp = "RNA") 
cite[["ADT"]] <- CreateAssayObject(counts = adt@assays$originalexp@data)
%%R
DefaultAssay(cite) <- "RNA"
VariableFeatures(cite) <- rownames(cite)
cite <- SCTransform(cite, verbose = FALSE) 
%%R
DefaultAssay(cite) <- "ADT"
VariableFeatures(cite) <- rownames(cite)
cite <- NormalizeData(cite, normalization.method = 'CLR', margin = 2, verbose=FALSE)
cite <- ScaleData(cite, verbose=FALSE)

对 ADT 查询采用同样的 CLR 归一化与标准化，再运行 PCA，将其命名为 apca。

adata_ = sc.AnnData(adt_cite_query_bridge.layers["counts"].A.copy())
adata_.obs_names = adt_cite_query_bridge.obs_names.copy()
adata_.var_names = adt_cite_query_bridge.var_names.copy()
adata_.obs["cell_type"] = adt_cite_query_bridge.obs["cell_type"].copy()
%%R -i adata_
query <- as.Seurat(adata_, data=NULL, counts='X')
query <- RenameAssays(object = query, originalexp = "ADT", verbose=FALSE) 
VariableFeatures(query) <- rownames(query)
query <- NormalizeData(query, normalization.method = 'CLR', margin = 2, verbose=FALSE)
query <- ScaleData(query, verbose=FALSE)
query <- RunPCA(query, verbose=FALSE, reduction.name = 'apca')

先构建扩展参考，再寻找桥接转移锚点，最后用 MapQuery 将 ADT 查询接入 RNA 参考，并预测细胞类型。各接口的参数与对象结构需对应安装的 Seurat 版本。

%%R
dims.adt <- 1:50
dims.rna <- 1:50

DefaultAssay(cite) <-  "RNA"
DefaultAssay(ref) <- "SCT"
obj.rna.ext <- PrepareBridgeReference(
    reference = ref,
    bridge = cite,
    bridge.query.assay = "ADT",
    reference.reduction = "pca",
    reference.dims = dims.rna,
    normalization.method = "SCT",
    verbose = FALSE
)
  |                                                  | 0 % ~calculating   |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s  
  |                                                  | 0 % ~calculating   |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s  

%%R
bridge.anchor <- FindBridgeTransferAnchors(
    extended.reference = obj.rna.ext, 
    query = query,
    reduction = "pcaproject",
    dims = dims.adt,
    verbose = FALSE
)
  |                                                  | 0 % ~calculating   |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s  

%%R
query <- MapQuery(
    anchorset = bridge.anchor, 
    reference = ref, 
    query = query, 
    refdata = list(
        cell_type = "cell_type"
    ),
    reduction.model = "umap",
    verbose = FALSE
)

比较查询的预测标签与已有注释，并将参考和查询显示在同一坐标系中。

%%R
p1 <- DimPlot(query, group.by = "predicted.cell_type", reduction = "ref.umap", label = TRUE) + NoLegend()
p2 <- DimPlot(query, group.by = "cell_type", reduction = "ref.umap", label = TRUE) + NoLegend()
p1 + p2
Image produced in Jupyter
%%R
ref$id <- 'reference'
query$id <- 'query'
refquery <-  merge(ref, query)
refquery[["umap"]] <- merge(ref[["umap"]], query[["ref.umap"]])
DimPlot(refquery, group.by = 'id', shuffle = TRUE)
Image produced in Jupyter

使用 multigrate 构建三模态参考并映射查询

最后用 RNA、ATAC 和 ADT 构建三模态参考图谱。multigrate Lotfollahi et al., 2022 通过各模态的编码器得到潜在表示的分布参数，再用专家乘积（Product of Experts, PoE）组合各模态的贡献 Lee & Schaar, 2021。条件解码器重建各自输入，并结合整合目标处理技术差异。缺失模态通过掩码排除其专家贡献，而不是把缺失测量当作真实零值重建；得到的是共享潜变量的近似后验，不能把它简单理解为原始模态的边缘分布直接相乘。multigrate 还采用单细胞架构手术（Single-Cell Architectural Surgery, scArches）Lotfollahi et al., 2022 的思路，适配新的单模态或多模态查询。本章环境从 main 安装 multigrate，未锁定提交；示例参数可能与后来版本不同，运行前需匹配接口。

构建参考

重新加载 RNA 数据：CITE-seq 与 Multiome 使用相同的 4,000 个高变基因（Highly Variable Genes, HVGs）。分别按对应蛋白或 ATAC 对象的细胞名称提取参考和查询，并确认两套 RNA 的基因名称及列顺序一致。

rna1 = sc.read(
    "/lustre/groups/ml01/workspace/anastasia.litinetskaya/data/trimodal_neurips/rna_hvg_cite.h5ad"
)
rna1_ref = rna1[adt_cite.obs_names].copy()
rna1_query = rna1[adt_cite_query.obs_names].copy()
rna1_ref.shape, rna1_query.shape
((16294, 4000), (16046, 4000))
rna2 = sc.read(
    "/lustre/groups/ml01/workspace/anastasia.litinetskaya/data/trimodal_neurips/rna_hvg_multiome.h5ad"
)
rna2_ref = rna2[atac_multiome.obs_names].copy()
rna2_query = rna2[atac_multiome_query.obs_names].copy()
rna2_ref.shape, rna2_query.shape
((17243, 4000), (10331, 4000))

按“模态×数据集”组织输入：第一组细胞有 RNA 与 ADT，第二组有 RNA 与 ATAC，None 表示该数据集缺失此模态。RNA 取 counts 层、ATAC 取 cpm 层、ADT 取 .X。这里跨数据集共有的是 RNA，没有要求所有细胞都同时测量三种模态。

adata = mtg.data.organize_multiome_anndatas(
    adatas=[[rna1_ref, rna2_ref], [None, atac_multiome], [adt_cite, None]],
    layers=[["counts", "counts"], [None, "cpm"], [None, None]],
)
adata
AnnData object with n_obs × n_vars = 33537 × 24136 obs: 'GEX_n_genes_by_counts', 'GEX_pct_counts_mt', 'GEX_size_factors', 'GEX_phase', 'ADT_n_antibodies_by_counts', 'ADT_total_counts', 'ADT_iso_count', 'cell_type', 'batch', 'ADT_pseudotime_order', 'GEX_pseudotime_order', 'Samplename', 'Site', 'DonorNumber', 'Modality', 'VendorLot', 'DonorID', 'DonorAge', 'DonorBMI', 'DonorBloodType', 'DonorRace', 'Ethnicity', 'DonorGender', 'QCMeds', 'DonorSmoker', 'is_train', 'GEX_n_counts', 'GEX_n_genes', 'ATAC_nCount_peaks', 'ATAC_atac_fragments', 'ATAC_reads_in_peaks_frac', 'ATAC_blacklist_fraction', 'ATAC_nucleosome_signal', 'ATAC_pseudotime_order', 'technology', 'cell_type_l2', 'cell_type_l1', 'cell_type_l3', 'assay', 'split', 'group', 'balancing_weight', 'donor', 'log1p_n_genes_by_counts', 'log1p_total_counts', 'modality', 'n_counts', 'n_genes_by_counts', 'new_donor', 'outliers', 'total_counts', 'new_batch' var: 'modality' uns: 'modality_lengths' layers: 'counts'

注册 Modality 与 Samplename 两个分类协变量（covariate），并指定前 4,000 列为 RNA。协变量用于模型条件化及相关校正设置，不代表这些标签所关联的所有生物学差异都应被消除。

mtg.model.MultiVAE.setup_anndata(
    adata,
    categorical_covariate_keys=["Modality", "Samplename"],
    rna_indices_end=4000,
)

初始化并训练模型。RNA 使用 NB 重建损失，ATAC 与 ADT 使用均方误差（Mean Squared Error, MSE）；示例给 Kullback–Leibler 散度（Kullback–Leibler divergence, KL divergence）项与整合项指定权重，并以 Modality 为整合分组。原代码的 mmd="marginal" 属于旧接口；较新的 multigrate 接口使用 alignment_type 等参数，不能直接把旧关键字原样传入该版本。应按对应版本的文档核对整合损失设置。

model = mtg.model.MultiVAE(
    adata,
    losses=["nb", "mse", "mse"],
    loss_coefs={
        "kl": 1e-3,
        "integ": 10000,
    },
    integrate_on="Modality",
    mmd="marginal",
)
model.train()
输出
GPU available: True (cuda), used: True
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
HPU available: False, using: 0 HPUs
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]

Epoch 200/200: 100%|██████████| 200/200 [28:53<00:00,  8.43s/it, loss=2.34e+03, v_num=1]
`Trainer.fit` stopped: `max_epochs=200` reached.

Epoch 200/200: 100%|██████████| 200/200 [28:53<00:00,  8.67s/it, loss=2.34e+03, v_num=1]

提取潜在表示，并另存于 .obsm['latent_ref']，供后续对照；查询适配后再次编码时会覆盖 .obsm['latent']。随后计算参考的邻居图和 UMAP。

model.get_latent_representation()
adata.obsm["latent_ref"] = adata.obsm["latent"].copy()
adata
AnnData object with n_obs × n_vars = 33537 × 24136 obs: 'GEX_n_genes_by_counts', 'GEX_pct_counts_mt', 'GEX_size_factors', 'GEX_phase', 'ADT_n_antibodies_by_counts', 'ADT_total_counts', 'ADT_iso_count', 'cell_type', 'batch', 'ADT_pseudotime_order', 'GEX_pseudotime_order', 'Samplename', 'Site', 'DonorNumber', 'Modality', 'VendorLot', 'DonorID', 'DonorAge', 'DonorBMI', 'DonorBloodType', 'DonorRace', 'Ethnicity', 'DonorGender', 'QCMeds', 'DonorSmoker', 'is_train', 'GEX_n_counts', 'GEX_n_genes', 'ATAC_nCount_peaks', 'ATAC_atac_fragments', 'ATAC_reads_in_peaks_frac', 'ATAC_blacklist_fraction', 'ATAC_nucleosome_signal', 'ATAC_pseudotime_order', 'technology', 'cell_type_l2', 'cell_type_l1', 'cell_type_l3', 'assay', 'split', 'group', 'balancing_weight', 'donor', 'log1p_n_genes_by_counts', 'log1p_total_counts', 'modality', 'n_counts', 'n_genes_by_counts', 'new_donor', 'outliers', 'total_counts', 'new_batch', 'size_factors', '_scvi_batch' var: 'modality' uns: 'modality_lengths', '_scvi_uuid', '_scvi_manager_uuid' obsm: '_scvi_extra_categorical_covs', 'latent', 'latent_ref' layers: 'counts'
sc.pp.neighbors(adata, use_rep="latent")
sc.tl.umap(adata)
sc.pl.umap(adata, color=["cell_type", "Modality", "Samplename"], ncols=1, frameon=False)
<Figure size 727.8x1440 with 3 Axes>

映射到三模态参考

按与参考相同的模态顺序、基因集合、数据层和协变量设置组织查询。

query = mtg.data.organize_multiome_anndatas(
    adatas=[
        [rna1_query, rna2_query],
        [None, atac_multiome_query],
        [adt_cite_query, None],
    ],
    layers=[["counts", "counts"], [None, "cpm"], [None, None]],
)
mtg.model.MultiVAE.setup_anndata(
    query,
    categorical_covariate_keys=["Modality", "Samplename"],
    rna_indices_end=4000,
)

模拟单模态查询：将 site2_donor4_multiome 的前 4,000 列 RNA 置零，保留 ATAC；将 site2_donor1_cite 的 RNA 之后全部列置零，其 ATAC 本来就缺失，因此得到仅 RNA 查询。其余两个批次保留配对数据。这些零用于标记缺失模态；实际运行后应检查修改确已写回 query，且与所用版本的缺失掩码约定一致。

idx_atac_query = query.obs["Samplename"] == "site2_donor4_multiome"
idx_scrna_query = query.obs["Samplename"] == "site2_donor1_cite"

idx_mutiome_query = query.obs["Samplename"] == "site2_donor1_multiome"
idx_cite_query = query.obs["Samplename"] == "site2_donor4_cite"

(
    np.sum(idx_atac_query),
    np.sum(idx_scrna_query),
    np.sum(idx_mutiome_query),
    np.sum(idx_cite_query),
)
(6111, 10465, 4220, 5581)
query[idx_atac_query, :4000].X = 0
query[idx_scrna_query, 4000:].X = 0

从参考模型加载查询，加入新批次类别所需的参数，并按查询适配时的冻结规则训练。本例设 weight_decay=0；具体哪些参数保持固定取决于版本实现。

q_model = mtg.model.MultiVAE.load_query_data(query, model)
q_model.train(weight_decay=0)
输出
GPU available: True (cuda), used: True
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
HPU available: False, using: 0 HPUs
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]

Epoch 200/200: 100%|██████████| 200/200 [16:13<00:00,  4.78s/it, loss=2.1e+03, v_num=1] 
`Trainer.fit` stopped: `max_epochs=200` reached.

Epoch 200/200: 100%|██████████| 200/200 [16:13<00:00,  4.87s/it, loss=2.1e+03, v_num=1]

用查询适配后的模型分别编码查询与参考。原教程称参考表示除采样波动外保持不变；实际应与 latent_ref 比较，核对参数冻结与随机采样设置，不能仅凭调用了 load_query_data 就认定结果必然完全一致。

q_model.get_latent_representation(adata=query)
q_model.get_latent_representation(adata=adata)
INFO     Input AnnData not setup with scvi-tools. attempting to transfer AnnData setup

合并参考与查询的潜在表示，重新构建邻居图并拟合 UMAP，查看不同查询类型的映射。重新计算 UMAP 会改变二维布局，不等于固定参考 UMAP 的投影。

adata.obs["reference"] = "reference"
query.obs["reference"] = "query"

adata.obs["type_of_query"] = "reference"
query.obs.loc[idx_atac_query, "type_of_query"] = "ATAC query"
query.obs.loc[idx_scrna_query, "type_of_query"] = "scRNA query"
query.obs.loc[idx_mutiome_query, "type_of_query"] = "multiome query"
query.obs.loc[idx_cite_query, "type_of_query"] = "CITE-seq query"
adata_both = anndata.concat([adata, query])
sc.pp.neighbors(adata_both, use_rep="latent")
sc.tl.umap(adata_both)
sc.pl.umap(
    adata_both,
    color=["cell_type", "Modality", "Samplename", "reference"],
    ncols=1,
    frameon=False,
)
<Figure size 727.8x1920 with 4 Axes>
sc.pl.umap(
    adata_both, color="type_of_query", ncols=1, frameon=False, groups=["ATAC query"]
)
<Figure size 640x480 with 1 Axes>
sc.pl.umap(
    adata_both, color="type_of_query", ncols=1, frameon=False, groups=["CITE-seq query"]
)
<Figure size 640x480 with 1 Axes>
sc.pl.umap(
    adata_both, color="type_of_query", ncols=1, frameon=False, groups=["multiome query"]
)
<Figure size 640x480 with 1 Axes>
sc.pl.umap(
    adata_both, color="type_of_query", ncols=1, frameon=False, groups=["scRNA query"]
)
<Figure size 640x480 with 1 Axes>

原教程图中，不同查询类型落在参考中相关细胞类型附近；即使参考没有仅含单一模态的细胞，模型仍可利用已有模态联系处理单模态查询。这是当前示例的结果，仍需评估新状态、参考缺失类型和跨批次泛化，不能将可视化邻近视为标签预测的充分验证。

会话信息

%%R
sessionInfo()
输出
R version 4.1.3 (2022-03-10)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /lustre/groups/ml01/workspace/anastasia.litinetskaya/miniconda3/envs/advanced-integration/lib/libopenblasp-r0.3.21.so

locale:
 [1] LC_CTYPE=C                 LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats4    tools     stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] Matrix_1.5-3                SeuratObject_4.1.3         
 [3] Seurat_4.1.1.9001           SingleCellExperiment_1.16.0
 [5] SummarizedExperiment_1.24.0 Biobase_2.54.0             
 [7] GenomicRanges_1.46.1        GenomeInfoDb_1.30.1        
 [9] IRanges_2.28.0              S4Vectors_0.32.4           
[11] BiocGenerics_0.40.0         MatrixGenerics_1.6.0       
[13] matrixStats_0.63.0         

loaded via a namespace (and not attached):
  [1] Rtsne_0.16             colorspace_2.0-3       deldir_1.0-6          
  [4] ellipsis_0.3.2         ggridges_0.5.4         XVector_0.34.0        
  [7] RcppHNSW_0.4.1         spatstat.data_3.0-0    farver_2.1.1          
 [10] leiden_0.4.3           listenv_0.9.0          ggrepel_0.9.2         
 [13] RSpectra_0.16-1        fansi_1.0.3            codetools_0.2-18      
 [16] splines_4.1.3          polyclip_1.10-4        jsonlite_1.8.4        
 [19] ica_1.0-3              cluster_2.1.4          png_0.1-8             
 [22] uwot_0.1.14            spatstat.sparse_3.0-0  shiny_1.7.4           
 [25] sctransform_0.3.5      compiler_4.1.3         httr_1.4.4            
 [28] fastmap_1.1.0          lazyeval_0.2.2         cli_3.5.0             
 [31] later_1.3.0            htmltools_0.5.4        igraph_1.3.5          
 [34] gtable_0.3.1           glue_1.6.2             GenomeInfoDbData_1.2.7
 [37] RANN_2.6.1             reshape2_1.4.4         dplyr_1.0.10          
 [40] Rcpp_1.0.9             scattermore_0.8        vctrs_0.5.1           
 [43] nlme_3.1-161           progressr_0.12.0       lmtest_0.9-40         
 [46] spatstat.random_3.0-1  stringr_1.5.0          globals_0.16.2        
 [49] mime_0.12              miniUI_0.1.1.1         lifecycle_1.0.3       
 [52] irlba_2.3.5.1          goftest_1.2-3          future_1.30.0         
 [55] zlibbioc_1.40.0        MASS_7.3-58.1          zoo_1.8-11            
 [58] scales_1.2.1           spatstat.core_2.4-4    promises_1.2.0.1      
 [61] spatstat.utils_3.0-1   parallel_4.1.3         RColorBrewer_1.1-3    
 [64] reticulate_1.26        pbapply_1.6-0          gridExtra_2.3         
 [67] ggplot2_3.4.0          rpart_4.1.19           stringi_1.7.8         
 [70] fastDummies_1.6.3      rlang_1.0.6            pkgconfig_2.0.3       
 [73] bitops_1.0-7           lattice_0.20-45        tensor_1.5            
 [76] ROCR_1.0-11            purrr_1.0.0            labeling_0.4.2        
 [79] patchwork_1.1.2        htmlwidgets_1.6.0      cowplot_1.1.1         
 [82] tidyselect_1.2.0       parallelly_1.33.0      RcppAnnoy_0.0.20      
 [85] plyr_1.8.8             magrittr_2.0.3         R6_2.5.1              
 [88] generics_0.1.3         DelayedArray_0.20.0    withr_2.5.0           
 [91] mgcv_1.8-41            pillar_1.8.1           fitdistrplus_1.1-8    
 [94] abind_1.4-5            survival_3.4-0         RCurl_1.98-1.9        
 [97] sp_1.5-1               tibble_3.1.8           future.apply_1.10.0   
[100] crayon_1.5.2           KernSmooth_2.23-20     utf8_1.2.2            
[103] spatstat.geom_3.0-3    plotly_4.10.1          grid_4.1.3            
[106] data.table_1.14.6      digest_0.6.31          xtable_1.8-4          
[109] tidyr_1.2.1            httpuv_1.6.7           munsell_0.5.0         
[112] viridisLite_0.4.1

贡献者

我们衷心感谢以下人员的贡献：

作者

  • Anastasia Litinetskaya

审阅者

  • Lukas Heumos

References
  1. Luecken, M. D., Burkhardt, D. B., Cannoodt, R., Lance, C., Agrawal, A., Aliee, H., Chen, A. T., Deconinck, L., Detweiler, A. M., Granados, A. A., Huynh, S., Isacco, L., Kim, Y. J., Klein, D., KUMAR, B. D., Kuppasani, S., Lickert, H., McGeever, A., Mekonen, H., … Bloom, J. M. (2021). A sandbox for prediction and integration of DNA, RNA, and proteins in single cells. Thirty-Fifth Conference on Neural Information Processing Systems Datasets and Benchmarks Track (Round 2). https://openreview.net/forum?id=gN35BGa1Rt
  2. Cao, Z.-J., & Gao, G. (2022). Multi-omics single-cell data integration and regulatory inference with graph-linked embedding. Nature Biotechnology. 10.1038/s41587-022-01284-4
  3. Gayoso, A., Steier, Z., Lopez, R., Regier, J., Nazor, K. L., Streets, A., & Yosef, N. (2021). Joint probabilistic modeling of single-cell multi-omic data with totalVI. Nature Methods, 18(3), 272–282. 10.1038/s41592-020-01050-x
  4. Hao, Y., Hao, S., Andersen-Nissen, E., Mauck, W. M., Zheng, S., Butler, A., Lee, M. J., Wilk, A. J., Darby, C., Zager, M., Hoffman, P., Stoeckius, M., Papalexi, E., Mimitou, E. P., Jain, J., Srivastava, A., Stuart, T., Fleming, L. M., Yeung, B., … Satija, R. (2021). Integrated analysis of multimodal single-cell data. Cell, 184(13), 3573-3587.e29. https://doi.org/10.1016/j.cell.2021.04.048
  5. Hao, Y., Stuart, T., Kowalski, M., Choudhary, S., Hoffman, P., Hartman, A., Srivastava, A., Molla, G., Madad, S., Fernandez-Granda, C., & Satija, R. (2022). Dictionary learning for integrative, multimodal, and scalable single-cell analysis. bioRxiv. 10.1101/2022.02.24.481684
  6. Lotfollahi, M., Litinetskaya, A., & Theis, F. J. (2022). Multigrate: single-cell multi-omic data integration. bioRxiv. 10.1101/2022.03.16.484643
  7. Lee, C., & van der Schaar, M. (2021). A Variational Information Bottleneck Approach to Multi-Omics Data Integration. Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, 130, 1513–1521. http://proceedings.mlr.press/v130/lee21a.html
  8. Lotfollahi, M., Naghipourfar, M., Luecken, M. D., Khajavi, M., Büttner, M., Wagenstetter, M., Avsec, Ž., Gayoso, A., Yosef, N., Interlandi, M., Rybakov, S., Misharin, A. V., & Theis, F. J. (2022). Mapping single-cell data to reference atlases by transfer learning. Nature Biotechnology, 40(1), 121–130. 10.1038/s41587-021-01001-7
多模态整合
配对整合
展望
展望