整合 AIR 与转录组学

43. 整合 AIR 与转录组学#

43.1. 动机#

随着新型单细胞技术的发展，AIRR 测序常常与转录组学等其他组学层结合。对于 B 细胞和 T 细胞，这使得我们能够在单细胞层面分析多种特征：转录组提供了细胞当前状态的信息，而 AIR 则提示细胞的特异性，从而解释细胞在感染或接种疫苗后的命运。

43.2. 数据准备#

首先，我们把 GEX 与 AIR 数据合并。注意，这一步也可以在单细胞分析的更早阶段进行：例如，可以在预处理之前就融合两种模态，以过滤 GEX 双细胞。在融合 GEX 与 AIR 数据时，通常执行左连接（left join），即只保留带有 GEX 的细胞用于分析。因此，融合与过滤的先后顺序并不太重要。不过，在早期阶段把 AIR 信息可视化到 GEX 的 UMAP 中会很方便。例如，聚类中缺失 AIR 信息可以在细胞类型指派时提供帮助。

警告

Scirpy 在 v0.13 中更改了其数据结构的格式。虽然整体分析流程没有改变，但本章中展示的一些输出可能不再准确。

详见 scirpy 发布说明了解此变化的更多细节。在本章更新之前，请参考 scirpy 官方文档。

import warnings

warnings.filterwarnings(
    "ignore",
    ".*IProgress not found*",
)
warnings.simplefilter(action="ignore", category=FutureWarning)
warnings.simplefilter(action="ignore", category=UserWarning)

import numpy as np
import pandas as pd
import scanpy as sc
import scirpy as ir

warnings.simplefilter(action="ignore", category=pd.errors.DtypeWarning)

path_data = "./data"
path_gex_tcr = f"{path_data}/TCR_00_GEX.h5ad"
path_gex_bcr = f"{path_data}/BCR_00_GEX.h5ad"
path_tcr = f"{path_data}/TCR_01_preprocessed.h5ad"
path_bcr = f"{path_data}/BCR_01_preprocessed.h5ad"

path_tmp = f"{path_data}/tmp"
path_res = f"{path_data}/res"

首先，我们加载由作者处理过的基因表达数据。你可以从以下位置下载完整数据集： https://www.ebi.ac.uk/arrayexpress/files/E-MTAB-10026/E-MTAB-10026.processed.4.zip。由于我们只使用选定供体的 B 细胞和 T 细胞，我们为你提供了一个下采样的数据集：

! wget -O $path_bcr_input -nc https://figshare.com/ndownloader/files/35574338

下面我们加载 GEX 和此前已注释的 TCR 数据。这两个 AnnData 对象随后可以通过 Scirpy 轻松合并为一个共享的 AnnData 对象，它会把 AIR 信息存储在 adata.obs 中、把 GEX 存储在 adata.X 中。这里执行的是左连接：只有来自带 scRNA 数据细胞的 AIR 才会被添加到 AnnData 对象中。带 scRNA 但没有 AIR 的细胞在分析中通常会保留，因为其中可能包含 B 细胞和 T 细胞以外的细胞。然而，带 AIR 但没有 scRNA 的细胞显然缺少测量数据；由于大多数下游分析是在转录组层面进行的，这些细胞会被排除。

adata_tc = sc.read(path_gex_tcr)
adata_tcr = sc.read(path_tcr)

ir.pp.merge_with_ir(adata_tc, adata_tcr)

下面我们在 T 细胞转录组上计算 Leiden 聚类。

sc.pp.neighbors(adata_tc)
sc.tl.leiden(adata_tc)

为了获得总体概览，我们将把带有聚类分配的数据绘制成 UMAP 可视化。

sc.tl.umap(adata_tc)
sc.pl.umap(adata_tc, color=["full_clustering", "has_ir", "leiden"], ncols=1)

../_images/b5990455a37baae69a548b32bf1da602b971b2ba65209bdc0349e6830c0bbab0.png

可以看到，大多数细胞属于 CD8+ 效应 T 细胞，并且大多数细胞都表达 TCR。AIR 的存在可以指导细胞类型注释，把其他免疫细胞的聚类与 T 细胞、B 细胞区分开来。此外，AIR 的缺失也可用于质量控制：例如，如果数据是按 T 细胞分选的，但某些大聚类却缺少 TCR，就应当从转录组层面进一步调查，看这些细胞是否被错误分选了。

类似地，我们针对 B 细胞合并 BCR 和 GEX 数据。

adata_bc = sc.read(path_gex_bcr)
adata_bcr = sc.read(path_bcr)

ir.pp.merge_with_ir(adata_bc, adata_bcr)

sc.pp.neighbors(adata_bc)
sc.tl.leiden(adata_bc)
sc.tl.umap(adata_bc)
sc.pl.umap(adata_bc, color=["full_clustering", "has_ir", "leiden"], ncols=1)

../_images/c677a00066909238335abc362c507ec9b0885a5bc6fcb52114b93ba44ac64e62.png

在这个可视化中，我们可以看到浆母细胞（Plasmablast）与 B 细胞之间有清晰的分离。

43.4. 多模态整合#

在上一节中，我们展示了如何在由其他模态得到的聚类上进行特定模态的分析。然而，这并未充分利用配对数据所提供的额外信息，因为 AIR 与基因表达是相互关联的：识别相同表位的适应性免疫细胞在激活后会经历相似的发育过程。与其他组学组合一样，近来人们对整合 AIR 与 GEX 以得到聚类或共享表示也产生了兴趣。

这些模型仍然较新，还不是标准“最佳实践”分析流程的一部分。我们认为这些方法或类似方法今后会被越来越多地使用，因此想把这些整合模型作为一种展望来介绍。不过，在应用这些模型时应当谨慎，因为它们仍处于实验阶段，其优缺点尚未经过独立评估。

43.4.1. 整合 TCR 与 GEX#

直到最近，才开发出三种联合利用 TCR 序列和转录组的方法，它们各自针对分析的不同方面：

TESSA：在 TCR 和 GEX 嵌入上应用贝叶斯模型，对转录组和序列相似的 T 细胞进行聚类。
Conga：通过图论方法，依据 GEX 和 TCR 的距离对 T 细胞克隆进行聚类
mvTCR：使用深度变分自编码器（Deep Variational Autoencoder）来得到 GEX 与 TCR 的共享嵌入。

所有多模态方法都需要配对的 GEX 和 AIR 数据。因此，我们先过滤掉所有没有 CDR3β 的细胞。

adata_tc = adata_tc[~adata_tc.obs["IR_VDJ_1_junction_aa"].isna()].copy()

43.4.1.1. TESSA#

Zhang 等人开发了 TESSA（TCR functional landscape estimation supervised with scRNA-seq analysis，即在 scRNA-seq 分析监督下的 TCR 功能景观估计） [Zhang et al., 2021]，它旨在通过贝叶斯建模，根据 T 细胞克隆的 TCR 序列和转录组对其进行嵌入和聚类。首先用一个预训练的自编码器把 CDR3β 序列压缩为 30 维的数值表示。随后，对各维度上调权重，使 TCR 表示与相似 TCR 组的基因表达相关联，从而量化各 TCR 位置对解释细胞基因表达的重要性。在一个迭代过程中，不断更新权重和分组，直到收敛，以在两种模态之间达到最大程度的一致。

在嵌入具有已知表位特异性的 T 细胞时，TESSA 产生了高纯度的聚类，数据来自 [10x Genomics, 2019]，优于常用于聚类 TCR 序列的单模态模型 GLIPH [Glanville et al., 2017]。此外，聚类中心性可指示亲合力更高的克隆——这类克隆表现为克隆扩增和较高的 ADT Count。把 TESSA 用于来自 [Yost et al., 2019]的数据时，作者在接受 PD-1 阻断治疗的患者中检测到新的应答 T 细胞聚类。

安装所需的代码和说明见 TESSA 仓库。

43.4.1.1.1. 数据预处理#

TESSA 需要若干特定格式的 TCR 和 GEX 文件。

首先，我们保存一个包含以下各列的 .csv 文件：

contig_id：以 Cell Barcode 作为索引列
cdr3：去掉起始 C 和结尾 F 之后的 CDR3β 氨基酸序列

df_tcr = adata_tc.obs
df_tcr.index.name = "contig_id"

# trimm cdr3 sequence
df_tcr["cdr3"] = [seq[1:-1] for seq in df_tcr["IR_VDJ_1_junction_aa"]]

# select only columns needed
df_tcr = df_tcr[["cdr3"]]

df_tcr.to_csv(f"{path_tmp}/TESSA_tcrs.csv")
df_tcr.head(5)

	cdr3
contig_id
S11_AAACCTGAGATTACCC-1	ASSLDARDRGRVTEAF
S11_AAACCTGGTAATTGGA-1	ASSPGTGTYGYT
S11_AAACCTGGTAGCACGA-1	ASSIPGAVHEQY
S11_AAACCTGGTCTCAACA-1	ASSLDARDRGRVTEAF
S11_AAACCTGTCGCCGTGA-1	ASSPQTGVARYGYT

我们将像 TESSA 论文中那样，筛选变异程度最高的 10% 的基因。

n_genes = adata_tc.shape[1] // 10
sc.pp.highly_variable_genes(adata_tc, n_top_genes=n_genes)
adata_tessa = adata_tc[:, adata_tc.var["highly_variable"]].copy()

现在我们把 GEX 矩阵存储为一个 “.csv” 文件。TESSA 需要 adata 对象的转置矩阵，即行代表不同的基因、列代表特定的细胞。为此，我们先把 Count 数据存储到一个带有正确索引和列名的 DataFrame 中，然后对该 DataFrame 进行转置。

count_mat = adata_tessa.X.A
df_counts = pd.DataFrame(count_mat)

df_counts.index = adata_tessa.obs.index
df_counts.index.name = ""
df_counts.columns = adata_tessa.var.index

df_counts = df_counts.transpose()

df_counts.to_csv(f"{path_tmp}/TESSA_gex.csv")
df_counts.head()

	S11_AAACCTGAGATTACCC-1	S11_AAACCTGGTAGCACGA-1	S11_AAACCTGGTCTCAACA-1	S11_AAACCTGTCGCCGTGA-1	S11_AAACCTGTCTTATCTG-1	S11_AAACGGGGTCGGCACT-1	S11_AAACGGGGTTACGGAG-1	...	S12_TTTGGTTGTCAGAGGT-1	S12_TTTGGTTTCTACTATC-1	S12_TTTGGTTTCTCGCTTG-1	S12_TTTGTCACACGGTAGA-1	S12_TTTGTCAGTAGAAAGG-1	S12_TTTGTCAGTGAGTGAC-1	S12_TTTGTCATCCACTCCA-1
KLHL17	0.000000	0.00000	0.000000	0.000000	0.0000	0.00000	0.000000	...	0.000000	0.000000	0.000000	1.662003	0.000000	0.000000	0.000000
HES4	0.000000	0.00000	0.000000	0.000000	0.0000	0.00000	0.000000	...	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
ISG15	1.845827	1.54457	2.441434	2.430966	0.9284	1.64661	1.565189	...	1.210719	1.581859	1.511785	0.000000	2.015501	2.088418	1.611682
AGRN	0.000000	0.00000	0.000000	0.000000	0.0000	0.00000	0.000000	...	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
TTLL10	0.000000	0.00000	0.000000	0.000000	0.0000	0.00000	0.000000	...	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000

5 rows × 5281 columns

43.4.1.1.2. 运行模型#

我们需要提供不同的设置选项，其中大多数用于指定输入或输出目录。我们先将它们汇总到一个字典中。

settings_full = {
    # Input files
    "tcr": f"{path_tmp}/TESSA_tcrs.csv",
    "exp": f"{path_tmp}/TESSA_gex.csv",
    # TESSA models
    "model": "TESSA/BriseisEncoder/TrainedEncoder.h5",
    "embeding_vectors": "TESSA/BriseisEncoder/Atchley_factors.csv",
    # Output files
    "output_TCR": f"{path_res}/TESSA_tcr_embedding.csv",
    "output_log": f"{path_res}/TESSA_log.log",
    "output_tessa": f"{path_res}",
    "within_sample_networks": "FALSE",
}

下一步，我们将通过指定环境并添加设置来创建运行 TESSA 的命令。

cmd_tessa = "source ~/.bashrc &&"
cmd_tessa += "conda activate TESSA && "

cmd_tessa += "python TESSA/Tessa_main.py"
for key, value in settings_full.items():
    cmd_tessa += f" -{key} {value}"
cmd_tessa

'source ~/.bashrc &&conda activate TESSA && python TESSA/Tessa_main.py -tcr ./data/tmp/TESSA_tcrs.csv -exp ./data/tmp/TESSA_gex.csv -model TESSA/BriseisEncoder/TrainedEncoder.h5 -embeding_vectors TESSA/BriseisEncoder/Atchley_factors.csv -output_TCR ./data/res/TESSA_tcr_embedding.csv -output_log ./data/res/TESSA_log.log -output_tessa ./data/res -within_sample_networks FALSE'

最后，我们通过调用该命令在指定设置下运行模型。视计算资源而定，这可能需要几分钟。由于命令行输出非常冗长，我们将通过 capture 行魔术命令（line magic）静默该单元格。调试时应删除第一行。

%%capture
!$cmd_tessa

43.4.1.1.3. 输出#

现在我们会得到三个输出文件，可进一步用于下游任务：

TESSA_tcr_embedding.csv：TCR 的嵌入表示（不含 GEX 信息）
result_meta.csv：逐细胞的聚类分配
tessa_final.RData：包含推断得到的模型参数（如加权向量 b）

这里，我们将前两个文件用作嵌入和聚类，并以 AnnData 格式存储以便处理。为避免聚类被高度扩增的克隆主导，我们去除重复的 TCR，将细胞级数据归并到克隆级。

tessa_embedding = pd.read_csv(f"{path_res}/TESSA_tcr_embedding.csv", index_col=0)
tessa_embedding = tessa_embedding.drop_duplicates()
tessa_obs = adata_tc[adata_tc.obs.index.isin(tessa_embedding.index)].obs
tessa_embedding = sc.AnnData(X=tessa_embedding.values, obs=tessa_obs)

接下来，我们将 TESSA 的聚类分配添加到 AnnData 对象中。由于 TESSA 在内部会略微改动 Barcode，我们需要为此做出调整。此外，我们会去除每个聚类内重复的 TCR。

clustering = pd.read_csv(f"{path_res}/result_meta.csv", index_col=0)
clustering.index = clustering.index.str.replace(".", "-")
clustering = clustering[clustering.index.isin(tessa_embedding.obs.index)]
tessa_embedding.obs["TESSA_cluster"] = clustering["cluster_number"].values

print(f"Unique clones: {tessa_embedding.obs['cdr3'].nunique()}")
print(f"Unique clusters: {tessa_embedding.obs['TESSA_cluster'].nunique()}")

Unique clones: 751
Unique clusters: 134

可以看到，751 个克隆被归入 144 个聚类，这些聚类具有相似的 TCR 序列，很可能在功能上相关。在 t-SNE 可视化中，我们将突出显示 TCR 嵌入空间中最大的十个聚类。

sc.pp.neighbors(tessa_embedding)
sc.tl.tsne(tessa_embedding)

top_10_clusters = tessa_embedding.obs["TESSA_cluster"].value_counts().head(10).index
tessa_embedding.obs["TESSA_cluster_top10"] = tessa_embedding.obs["TESSA_cluster"].apply(
    lambda x: x if x in top_10_clusters else np.nan
)

sc.pl.tsne(tessa_embedding, color="TESSA_cluster_top10")

WARNING: Consider installing the package MulticoreTSNE (https://github.com/DmitryUlyanov/Multicore-TSNE). Even for n_jobs=1 this speeds up the computation considerably and might yield better converged results.
... storing 'TESSA_cluster' as categorical
... storing 'TESSA_cluster_top10' as categorical

../_images/a809bf3731309b7fc85567f294b7ced064ef68bcc20f7d1d5c4117e39ba40abe.png

为了在 GEX 空间中绘制聚类分配，我们将该分配添加到原始 adata 对象上。绘图后可以看到，TESSA 聚类在 GEX 层面也共享相似的表型。

mapping_dict = dict(
    tessa_embedding.obs[["IR_VDJ_1_junction_aa", "TESSA_cluster_top10"]].values
)
adata_tc.obs["TESSA_cluster_top10"] = adata_tc.obs["IR_VDJ_1_junction_aa"].map(
    mapping_dict
)
sc.pl.umap(adata_tc, color="TESSA_cluster_top10")

... storing 'cdr3' as categorical
... storing 'TESSA_cluster_top10' as categorical

../_images/92d7b055af9c6b9efd1746d801b16d5938ec4174557cdd7582b1563da6b633af.png

在两种可视化中，我们都能看到聚类在 TCR 和 GEX 层面是相关的。然而，完整的聚类结构无法在单一模态中直接观察到。由于 TCR 与 GEX 图谱相似，这些细胞可能针对相同的表位。例如，我们可以将这一注释用于上文所述的细胞网络之间的 DEG 分析。

43.4.1.2. CoNGA#

克隆型邻接图分析（CoNGA） [Schattgen et al., 2022] 在克隆型层面利用 GEX 和 TCR 上的相似度图，其中使用了著名的距离度量 TCRDist [Dash et al., 2017]。再依据图的邻域，根据共享的 GEX 与 TCR 归属来形成聚类。这些所谓的“CoNGA 聚类”因而具有相似的受体序列以及相似的基因表达。CoNGA 聚类中包含的克隆型，其 TCR 在长度和理化性质等方面高度相似，并被证明能够捕捉特异性。下面的说明将运行主要的 CoNGA 流程。由于流程的不同方面已在 Colab 教程中解释，以下描述主要侧重于如何转换数据并运行分析流程。

43.4.1.2.1. 数据准备#

CoNGA 方便地支持多种基因表达矩阵输入类型，包括 h5ad 文件。对于克隆型信息，我们需要先在 Cell Ranger 输出上运行一个预处理脚本来创建 clone_file。

df_tcr = pd.read_csv(
    f"{path_data}/TCR_00_read_aligned.csv",
    # index_col='barcode',
)
df_tcr = df_tcr.drop("Unnamed: 0", axis=1)
df_tcr = df_tcr[df_tcr["barcode"].isin(adata_tc.obs.index)]

path_reduced_tcr = f"{path_tmp}/conga_tcrs.csv"
df_tcr.to_csv(path_reduced_tcr)

path_conga_clones = f"{path_tmp}/conga_clone_file.tsv"

cmd_conga_pp = (
    "python conga/scripts/setup_10x_for_conga.py "
    f"--filtered_contig_annotations_csvfile {path_reduced_tcr} "
    f"--output_clones_file {path_conga_clones} "
    "--organism human"
)

%%capture
!$cmd_conga_pp

CoNGA 流程需要未取对数的基因表达 Count。由于我们的 SARS-CoV-2 数据集未提供该数据，我们将用其逆变换还原 log1p 变换，并将结果保存到单独的文件中。如果在自己的数据集上运行分析，可以直接链接到原始 Count 数据。

path_conga_gex = f"{path_tmp}/conga_gex.h5ad"

adata_conga = adata_tc.copy()
adata_conga.X = np.expm1(adata_tc.X)
sc.write(adata=adata_conga, filename=path_conga_gex)

43.4.1.2.2. 执行 CoNGA#

我们可以通过以下命令行调用来执行主要的 CoNGA 流程，该调用指定了输入和输出目录。

cmd_conga = f"cd {path_res} && "

cmd_conga += (
    "python ../../conga/scripts/run_conga.py "
    f"--graph_vs_graph --graph_vs_graph_stats --graph_vs_features --tcr_clumping --find_hotspot_features "
    f"--gex_data ../../{path_conga_gex} --gex_data_type h5ad "
    f"--clones_file ../../{path_conga_clones} "
    f"--organism human --outfile_prefix conga"
)

%%capture
!$cmd_conga

43.4.1.2.3. 输出#

分析流程的结果存储在 ./data/res/conga_results_summary.html 中，并附有对每个结果的解释。其他注释和图表见 ./data/res_conga/。解释所有结果超出了本教程的范围。因此，我们建议参阅 web summary 和 colab notebook 中非常详细的解释。

43.4.1.3. mvTCR#

An 等人提出的 mvTCR 是一种多视图变分自编码器，可将 TCR 序列和基因表达压缩为低维表示 [An et al., 2021]。两种深度学习架构——Transformer 和多层感知机（MLP）——分别从 TCR 和 GEX 中提取信息，然后融合以得到联合空间。随后，训练好的模型可用于嵌入相似的数据。

作者表明，相比单模态嵌入，多模态模型能更好地捕捉抗原特异性；用于预测和聚类的数据来自 [10x Genomics, 2019]。此外，他们还表明，在 SARS-CoV-2 数据集上，细胞类型和细胞功能在嵌入空间中得以保留；该数据集来自 [Fischer et al., 2021]。代码和进一步说明见 mvTCR 仓库。

mvTCR 依赖 AnnData 格式，因此我们需要为训练准备好 AnnData 对象。为此，mvTCR 提供了若干函数。

import sys

sys.path.append("mvTCR")
sys.path.append(".")

import tcr_embedding.utils_training as utils
from tcr_embedding.models.model_selection import run_model_selection
from tcr_embedding.utils_preprocessing import encode_tcr, group_shuffle_split

AnnData 对象需要克隆型分配（见克隆型分析）。由于 mvTCR 使用主受体的两条链，我们也将利用这些信息来定义克隆型。# todo => 使用已提供的

adata_mvtcr = adata_tc[~adata_tc.obs["IR_VJ_1_junction_aa"].isna()].copy()

# todo delete once data is unified
ir.tl.define_clonotypes(
    adata_mvtcr, key_added="clonotype", receptor_arms="all", dual_ir="primary_only"
)

100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 728/728 [00:00<00:00, 795.39it/s]

该模型需要 CDR3α 和 CDR3β 链的数值编码。为此，我们需要提供包含氨基酸序列的列。此外，pad 属性表示序列的最大长度。如果之后不打算嵌入更多数据，可以将其设为最大序列长度。

pad = max(
    adata_mvtcr.obs["IR_VJ_1_junction_aa"].str.len().max(),
    adata_mvtcr.obs["IR_VDJ_1_junction_aa"].str.len().max(),
)
encode_tcr(
    adata_mvtcr,
    column_cdr3a="IR_VJ_1_junction_aa",
    column_cdr3b="IR_VDJ_1_junction_aa",
    pad=pad,
)

数据需要划分为训练集和验证集。具体做法是设置值 train 或 val ，用于 adata.obs['set']中的每个细胞。通过选择 clonotype 和 0.2，用于评估模型性能的数据集将包含约 20% 的克隆型，且这些克隆型不与训练集重叠。此设置一般可用于数据分析。

train, val = group_shuffle_split(adata_mvtcr, group_col="clonotype", val_split=0.2)
adata_mvtcr.obs["set"] = "train"
adata_mvtcr.obs.loc[val.obs.index, "set"] = "val"

我们需要提供训练过程中使用的各种信息：

params_experiment = {
    "study_name": "test_haniffa",  # Name that identifies the study
    "model_name": "moe",  # Type of mixture model used during training, authors suggest using Mixture of Export (moe)
    "balanced_sampling": "clonotype",  # Column containing the id for clonotypes
    "comet_workspace": None,  # Can be used to log experiments via comet-ml
    "save_path": path_res,  # Output path, were the selected models are stored
}

mvTCR 会运行超参数优化来选择最佳模型架构。对于一般分析，优化模式 pseudo_metric 最为合适。此处，以能否保留数据集的多种特征来决定最佳模型。下面，我们选择相同的 TCR 序列（clonotype）与细胞类型（functional.cluster），两者权重相同（均为 1）。如果训练后潜空间被某一种模态主导，我们可以改用不同的权重重新运行模型。

params_optimization = {
    "name": "pseudo_metric",
    "prediction_labels": {"clonotype": 1, "full_clustering": 1},
}

最后，我们可以运行模型选择。尤其对于较大的数据集，作者使用了广泛的模型搜索。作为示例，这里我们只训练一个模型。建议使用带 GPU 的计算机。即便如此，模型仍需训练数分钟。

n_runs = 1
run_model_selection(adata_mvtcr, params_experiment, params_optimization, n_runs)

[I 2022-11-17 13:46:07,144] A new study created in RDB with name: test_haniffa
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 200/200 [08:12<00:00,  2.46s/it]
[I 2022-11-17 13:54:41,499] Trial 0 finished with value: 1.681318287535559 and parameters: {'dropout': 0.1, 'activation': 'linear', 'rna_hidden': 1500, 'hdim': 200, 'shared_hidden': 100, 'rna_num_layers': 1, 'tfmr_encoding_layers': 4, 'loss_weights_kl': 4.0428727350273357e-07, 'loss_weights_tcr': 0.034702669886504146, 'lr': 1.0994335574766187e-05, 'zdim': 50, 'tfmr_embedding_size': 16, 'tfmr_num_heads': 8, 'tfmr_dropout': 0.15000000000000002}. Best is trial 0 with value: 1.681318287535559.

Study statistics:
  Number of finished trials: 1
  Number of pruned trials: 0
  Number of complete trials: 1
Best trial: 
  trial_0
  Value: 1.681318287535559

现在我们有了一个训练好的模型，可用它将数据嵌入到同时捕捉 TCR 和 GEX 的表示中。为此，我们将从 studies 文件夹中加载最佳模型。

best_trial = 0
path_model = f"../data/res/trial_{best_trial}/best_model_by_metric.pt"
model = utils.load_model(adata_mvtcr, path_model)

下一步，我们现在可以嵌入这些数据，并将注释分配给新的、已嵌入的 AnnData 对象。

mvtcr_embedding = model.get_latent(adata_mvtcr, metadata=[])
mvtcr_embedding.obs = adata_mvtcr.obs.copy()

现在，这一嵌入可用于可视化等数据分析。为此，我们将选取样本中丰度最高的前 10 个克隆型。

top_10_clones = mvtcr_embedding.obs["clonotype"].value_counts().head(10).index
mvtcr_embedding.obs["clonotypes_top10"] = mvtcr_embedding.obs["clonotype"].apply(
    lambda x: x if x in top_10_clones else np.nan
)

现在，我们将使用 UMAP 来展示细胞类型分配和前 10 个克隆型。

sc.pp.neighbors(mvtcr_embedding)
sc.tl.umap(mvtcr_embedding)
sc.pl.umap(mvtcr_embedding, color=["full_clustering", "clonotypes_top10"], ncols=1)

... storing 'clonotype' as categorical
... storing 'set' as categorical
... storing 'clonotypes_top10' as categorical

../_images/0f0ef6bd2f55863515e93b9abf9ad157ad653e6c60ea627c032283ba468cf742.png

这里我们对 TCR-GEX 嵌入进行可视化。可以清楚地看到，数据按克隆型分离。不过，数据在克隆内部还形成了额外的亚聚类和结构，而这受 GEX 的影响。根据具体分析，你可能需要调整 GEX 与 TCR 之间的权重并运行更长时间的参数搜索。此外，现在也可以使用 Scanpy 函数在该嵌入上执行聚类等不同的下游任务。

43.4.2. 整合 BCR 与 GEX#

上述三种模型都是为 TCR 开发的。不过，由于 TCR 与 BCR 结构相似，CoNGA 和 mvTCR 或许也适用于 B 细胞。但这两个模型在其论文中均未针对 B 细胞进行评估。TESSA 依赖于预训练的 TCR 自编码器，因此不应用于 BCR-GEX 整合。作者最近发表了 Benisse，它与 TESSA 类似，但应用于 B 细胞 [Zhang et al., 2022].

与 TCR 类似，多模态 BCR 整合工具仍较为新颖，尚未用于标准分析流程。因此，应用这些模型时我们建议保持谨慎。

43.4.2.1. Benisse#

Benisse 使用预训练模型将 BCR 的重链嵌入为数值表示。该模型在克隆型层面运作：即对具有相同 BCR 的细胞的 GEX 取平均，得到每个克隆的代表性表示。随后，Benisse 通过构建稀疏图结构，检测在 BCR 和 GEX 层面相似的克隆。

由于该模型的设置与 TESSA 相似，我们将遵循相同的处理步骤，仅在细节与上文不同处加以说明。

由于运行时间较长，我们首先将 B 细胞数据集下采样为某一供体的 1000 个细胞。

adata_benisse = adata_bc[~adata_bc.obs["IR_VDJ_1_junction_aa"].isna()]
adata_benisse = adata_benisse[adata_benisse.obs["sample_id"] == "MH9143275"]
sc.pp.subsample(adata_benisse, n_obs=1000, random_state=0)
print(f"Amount of cells: {len(adata_benisse)}")

Amount of cells: 1000

df_bcr = adata_benisse.obs
df_bcr.index.name = "contigs"

df_bcr = df_bcr.rename(columns={"IR_VDJ_1_junction_aa": "cdr3"})

df_bcr = df_bcr[["cdr3"]]

df_bcr.to_csv(f"{path_tmp}/BENISSE_bcrs.csv")
df_bcr.head(5)

	cdr3
contigs
TTCTCAAAGAATTCCC-MH9143275	CAGGSQWEVKLDYW
TATCAGGAGTGAACAT-MH9143275	CARPRSLIAAAGAFDIW
CATCCACCAACGCACC-MH9143275	CASEEVTMVRGVMFPYGMDVW
GTCGTAATCACCGTAA-MH9143275	CAREGMVYVDYW
AACTCTTAGCACCGCT-MH9143275	CATQNPGYSSSWDDRGAFDIW

我们将筛选变异程度最高的 1000 个基因（与 GitHub 仓库样本数据文件中提供的一致）。

sc.pp.highly_variable_genes(adata_benisse, n_top_genes=1000)
adata_benisse = adata_benisse[:, adata_benisse.var["highly_variable"]].copy()

count_mat = adata_benisse.X.A
df_counts = pd.DataFrame(count_mat)

df_counts.index = adata_benisse.obs.index
df_counts.index.name = ""
df_counts.columns = adata_benisse.var.index

df_counts = df_counts.transpose()

df_counts.to_csv(f"{path_tmp}/BENISSE_gex.csv")
df_counts.head()

	TTCTCAAAGAATTCCC-MH9143275	GTCGTAATCACCGTAA-MH9143275	ACTATCTCACAGGAGT-MH9143275	CAGTAACCAGGAATCG-MH9143275	ATTGGACAGCCATCGC-MH9143275	...	TACACGATCGGGAGTA-MH9143275	TGCGCAGTCTCAAGTG-MH9143275
HES4	0.000000	0.000000	0.00000	0.000000	0.00000	...	0.000000	0.000000
CPTP	0.334078	0.000000	0.00000	0.000000	0.44491	...	0.565442	0.000000
THAP3	0.000000	0.225355	0.45654	0.000000	0.00000	...	0.000000	0.000000
VAMP3	0.334078	0.000000	0.00000	0.000000	0.00000	...	0.565442	0.000000
CASP9	0.000000	0.225355	0.00000	0.355897	0.00000	...	0.000000	0.291166

5 rows × 1000 columns

与 TESSA 不同，Benisse 还需要一个根据 10x Cell Ranger 输出整理的 contigs 文件。这里，我们将过滤原始 contigs 文件，使其仅包含所选供体的 Barcode。

path_bcr_contigs = f"{path_data}/BCR_00_read_aligned.csv"
df_bcr_contigs = pd.read_csv(path_bcr_contigs, index_col=0, low_memory=False)
df_bcr_contigs = df_bcr_contigs[df_bcr_contigs["barcode"].isin(adata_benisse.obs.index)]
df_bcr_contigs.to_csv(f"{path_tmp}/BENISSE_bcr_contigs.csv")

43.4.2.1.1. 运行模型#

首先，我们需要通过预训练的编码器嵌入 BCR 序列：

settings_embedding = {
    "input_data": f"{path_tmp}/BENISSE_bcrs.csv",
    "output_data": f"{path_res}/BENISSE_encoded_bcrs.csv",
    "cuda": "True",
}

cmd_bcr_embedding = "source ~/.bashrc && "
cmd_bcr_embedding += "conda activate BENISSE && "

cmd_bcr_embedding += "python Benisse/AchillesEncoder.py"
for key, value in settings_embedding.items():
    cmd_bcr_embedding += f" --{key} {value}"
cmd_bcr_embedding

'source ~/.bashrc && conda activate BENISSE && python Benisse/AchillesEncoder.py --input_data ./data/tmp/BENISSE_bcrs.csv --output_data ./data/res/BENISSE_encoded_bcrs.csv --cuda True'

%%capture
!$cmd_bcr_embedding

接下来，我们需要创建图表示。这里我们将使用作者提供的默认值。

max_iter = 100
cmd_benisse = [
    "source ~/.bashrc &&",
    "conda activate BENISSE &&",
    "export R_LIBS=$CONDA_PREFIX/lib/R/library &&",
    "Rscript Benisse/Benisse.R",
    f"{path_tmp}/BENISSE_gex.csv",
    f"{path_tmp}/BENISSE_bcr_contigs.csv",
    f"{path_res}/BENISSE_encoded_bcrs.csv",
    f"{path_res}",
    f"1610 1 {max_iter} 1 1 10 1e-10",
]
cmd_benisse = " ".join(cmd_benisse)
cmd_benisse

'source ~/.bashrc && conda activate BENISSE && export R_LIBS=$CONDA_PREFIX/lib/R/library && Rscript Benisse/Benisse.R ./data/tmp/BENISSE_gex.csv ./data/tmp/BENISSE_bcr_contigs.csv ./data/res/BENISSE_encoded_bcrs.csv ./data/res 1610 1 100 1 1 10 1e-10'

即使减少了迭代次数，运行这段代码仍可能需要几分钟。

%%capture
!$cmd_benisse

43.4.2.1.2. 输出#

现在我们会得到几个输出文件，作者的 GitHub 仓库中对此有更详细的描述。下面我们将考察：

connectionplot.pdf：所得聚类的可视化
clone_annotation.csv：逐细胞的聚类分配

from IPython.display import IFrame

IFrame(f"{path_res}/connectionplot.pdf", width=600, height=600)

连接图（connection plot）将 Benisse 所构建的图可视化呈现出来。每个克隆型由学习到的嵌入空间中的一个节点表示。该聚类同时依据 BCR 与基因表达的相似性。

benisse_clusters = pd.read_csv(f"{path_res}/clone_annotation.csv", index_col=0)
benisse_clusters.index = [el.replace(".", "-") for el in benisse_clusters.index]
adata_benisse_out = adata_benisse[
    adata_benisse.obs.index.isin(benisse_clusters.index)
].copy()
adata_benisse_out.obs["benisse_clusters"] = benisse_clusters.loc[
    adata_benisse_out.obs.index
]["graph_label"]
sc.pp.neighbors(adata_benisse_out)
sc.tl.umap(adata_benisse_out)

top_10_clusters = (
    adata_benisse_out.obs["benisse_clusters"].value_counts().head(10).index
)
adata_benisse_out.obs["clusters_top10"] = adata_benisse_out.obs[
    "benisse_clusters"
].apply(lambda x: x if x in top_10_clusters else np.nan)
sc.pl.umap(adata_benisse_out, color="clusters_top10")

... storing 'benisse_clusters' as categorical
... storing 'clusters_top10' as categorical

../_images/61bcb31d35d1c46b5b455fdcb1446f5a2c8a981d96cdfac99deb171b667f4750.png

这里，我们将最大的十个聚类投影到转录组空间上。虽然有些聚类表现出相同的表型，但另一些聚类在 RNA 空间中相互分离。该聚类注释可用于 BCR 或转录组层面的下游分析（见上文），或用于推断细胞之间的祖先关系（见论文）。

43.5. Quiz#

为什么把 AIR 序列信息与基因表达整合起来会有用？

GEX 可用于改进 AIR 序列的 Read；两种模态各自从不同角度揭示细胞的信息，同时又彼此关联，等等。

问：为什么把 AIR 序列信息与基因表达整合起来会有用？答： GEX 可用于改进 AIR 序列的 Read；两种模态各自从不同角度揭示细胞的信息，同时又彼此关联，等等。

AIR 序列能为我们提供哪些 GEX 中没有直接捕获的信息？

细胞与抗体标记表位结合之间的 Count 矩阵；IR 序列可用于在不同供体之间进行解复用（demultiplexing），等等。

问： AIR 序列能为我们提供哪些 GEX 中没有直接捕获的信息？答：细胞与抗体标记表位结合之间的 Count 矩阵；IR 序列可用于在不同供体之间进行解复用（demultiplexing），等等。

多模态整合方法依赖于什么前提？

具有相同或相似 AIR 的细胞往往具有相似的表型；由于 AIR 与 GEX 相互独立，二者提供彼此正交的信息，等等。

问：多模态整合方法依赖于什么前提？答：具有相同或相似 AIR 的细胞往往具有相似的表型；由于 AIR 与 GEX 相互独立，二者提供彼此正交的信息，等等。

43.6. 参考文献#

[10xG19] (1,2)

10x Genomics. A new way of exploring immunity–linking highly multiplexed antigen recognition to immune repertoire and phenotype. Tech. rep, 2019.

[ADT+21]

Yang An, Felix Drost, Fabian Theis, Benjamin Schubert, and Mohammad Lotfollahi. Jointly learning t-cell receptor and transcriptomic information to decipher the immune response. bioRxiv, 2021.

[DFGH+17]

Pradyot Dash, Andrew J Fiore-Gartland, Tomer Hertz, George C Wang, Shalini Sharma, Aisha Souquette, Jeremy Chase Crawford, E Bridie Clemens, Thi HO Nguyen, Katherine Kedzierska, and others. Quantifiable predictive features define epitope-specific t cell receptor repertoires. Nature, 547(7661):89–93, 2017.

[FAW+21]

David S Fischer, Meshal Ansari, Karolin I Wagner, Sebastian Jarosch, Yiqi Huang, Christoph H Mayr, Maximilian Strunz, Niklas J Lang, Elvira D’Ippolito, Monika Hammel, and others. Single-cell rna sequencing reveals ex vivo signatures of sars-cov-2-reactive t cells through ‘reverse phenotyping’. Nature communications, 12(1):1–14, 2021.

[GHN+17]

Jacob Glanville, Huang Huang, Allison Nau, Olivia Hatton, Lisa E Wagar, Florian Rubelt, Xuhuai Ji, Arnold Han, Sheri M Krams, Christina Pettus, and others. Identifying specificity groups in the t cell receptor repertoire. Nature, 547(7661):94–98, 2017.

[SGC+22]

Stefan A Schattgen, Kate Guion, Jeremy Chase Crawford, Aisha Souquette, Alvaro Martinez Barrio, Michael JT Stubbington, Paul G Thomas, and Philip Bradley. Integrating t cell receptor sequences and transcriptional profiles by clonotype neighbor graph analysis (conga). Nature Biotechnology, 40(1):54–63, 2022.

[YSW+19]

Kathryn E Yost, Ansuman T Satpathy, Daniel K Wells, Yanyan Qi, Chunlin Wang, Robin Kageyama, Katherine L McNamara, Jeffrey M Granja, Kavita Y Sarin, Ryanne A Brown, and others. Clonal replacement of tumor-specific t cells following pd-1 blockade. Nature medicine, 25(8):1251–1259, 2019.

[ZCW+22]

Ze Zhang, Woo Yong Chang, Kaiwen Wang, Yuqiu Yang, Xinlei Wang, Chen Yao, Tuoqi Wu, Li Wang, and Tao Wang. Interpreting the b-cell receptor repertoire with single-cell gene expression using benisse. Nature Machine Intelligence, pages 1–9, 2022.

[ZXW+21]

Ze Zhang, Danyi Xiong, Xinlei Wang, Hongyu Liu, and Tao Wang. Mapping the functional landscape of t cell receptor repertoires by single-t cell transcriptomics. Nature methods, 18(1):92–99, 2021.