10. 数据整合
10.1 总论
大多数scRNA-seq数据分析的一个核心挑战是批次效应。批次效应是测量的表达水平的变化，这是处理不同组或“批次”中的细胞的结果。例如，如果两个实验室从同一队列中采集样本，但这些样本的解离方式不同，则可能会出现批次效应。如果实验室A优化其解离方案以解离样品中的细胞，同时最大限度地减少对细胞的压力，而实验室B没有这样做，那么B组数据中的细胞很可能会表达更多与压力相关的基因（JUN、 JUNB、FOS 等），即使细胞在原始组织中具有相同的profile。一般来说，批次效应的根源多种多样且难以确定。一些批次效应来源可能是技术性的，例如样本处理、实验方案或测序深度的差异，但供体变异、组织或采样位置等生物效应通常也被解释为批次效应。是否应考虑生物因素批次效应取决于实验设计和所提出的问题。消除批次效应对于实现联合分析至关重要，联合分析可以专注于查找跨批次数据中的共同结构，并使我们能够跨数据集执行查询。通常，只有在消除这些影响后，才能识别出以前因批次之间的差异而被掩盖的稀有细胞群。启用跨数据集查询使我们能够提出通过分析单个数据集无法回答的问题。
当从组学数据中消除批次效应时，必须做出两个主要选择：（1）方法和参数化parameterization，以及（2）批次协变量batch covariate。由于不同级别（即样本、供体、数据集等）的细胞分组之间可能会出现批次效应，因此批次协变量的选择表明应保留哪个级别的变异以及应删除哪个级别。批次分辨率越精细，消除的影响就越多。然而，精细批次变化也更有可能与具有生物学意义的信号混淆。例如，样本通常来自不同个体或组织中的不同位置。这些影响可能对检查有意义。因此，批次协变量的选择将取决于整合任务的目标。您想了解个体之间的差异，还是关注常见的细胞类型变异？最近在构建人肺综合图谱的过程中，开创了一种基于定量分析的批量协变量选择方法，其中使用不同技术协变量的方差来做出这种选择。
10.1.1 整合模型的类型
在scRNA-seq中消除批次效应的方法通常由（最多）三个步骤组成：
- 降维
- 建模并消除批次效应
- 投影回高维空间
虽然建模和消除批次效应（步骤 2）是任何批次消除方法的核心部分，但许多方法首先将数据投影到较低维度的空间（步骤 1）以提高信噪比并在该空间中执行批量校正以获得更好的性能。第三步，方法可以在去除拟合的批次效应后将数据投影回原始高维特征空间，以输出批次校正的基因表达矩阵。
批量效应消除方法在这三个步骤中的每一个步骤中都可能有所不同。它们可以使用各种线性或非线性降维方法、线性或非线性批量效应模型，并且它们可以输出不同格式的批量校正数据。总的来说，我们可以将批量效应消除方法分为 4 类。按照发展顺序，它们是全局模型、线性嵌入模型、基于图的方法和深度学习方法（图1）。
全局模型源自bulk转录组学，并将批次效应建模为所有细胞中一致（加法和/或乘法）效应。 一个常见的例子是ComBat。
线性嵌入模型是第一个单细胞特异性批量去除方法。这些方法通常使用奇异值分解 (singular value decomposition，SVD) 的变体来嵌入数据，然后在嵌入中跨批次查找相似单元的局部邻域，并使用它们以局部自适应（非线性）方式校正批次效应。方法通常使用SVD载荷将数据投影回基因表达空间，但也可能仅输出正确的嵌入。这是最常见的一组方法，突出的例子包括开创性的相互最近邻 (MNN) 方法（不执行任何降维）、Seurat整合、Scanorama、FastMNN和Harmony。
基于图的方法通常是运行速度最快的方法。这些方法使用最近邻图来表示每个批次的数据。通过强制不同批次的细胞之间的连接，然后通过修剪强制边缘来允许细胞类型组成的差异，可以纠正批次效应。这些方法最突出的例子是批量平衡k最近邻 (BBKNN) 方法。
深度学习(DL)方法是最新、最复杂的批量效应消除方法，通常需要最多的数据才能获得良好的性能。大多数深度学习集成方法都基于自动编码器网络，并且要么在条件变分自动编码器（conditional variational autoencoder，CVAE）中对批量协变量进行降维，要么在嵌入空间中拟合局部线性校正。深度学习方法的著名例子有scVI、scANVI和scGen。
一些方法可以使用细胞身份标签来为该方法提供参考，以了解哪些生物变异不应作为批次效应而被去除。由于批量效应去除通常是一项预处理任务，因此此类方法可能不适用于许多集成场景，因为现阶段标签通常尚不可用。
批量效应去除方法的更详细概述可以在Argelaguet et al., 2021和Luecken et al., 2021中找到。
10.1.2 批量去除复杂度
scRNA-seq数据中批量效应的消除之前被分为两个子任务：批次校正和数据集成。这些子任务的不同之处在于必须删除的批次效应的复杂性。批次校正方法处理同一实验中样品之间的批次效应，其中细胞身份组成一致，并且效应通常是准线性的。相比之下，数据集成方法处理可能使用不同protocol生成的数据集之间复杂的、通常嵌套的批次效应，并且细胞身份可能无法跨批次共享。虽然我们在这里使用这种区别，但我们不应该认为这些术语在一般用途中经常互换使用。
10.1.3 数据集成方法比较
之前有几个基准测试评估了批次校正和数据集成方法的性能。当消除批次效应时，方法可能会过度校正并消除除批次效应之外的有意义的生物变异。因此，通过考虑批次效应消除和生物变异的保存来评估集成性能非常重要。
k最近邻批量效应测试 (kBET) 是量化scRNA-seq数据批量校正的第一个指标。作者使用kBET发现，在比较主要全局模型时，ComBat优于其他批量校正方法。线性嵌入模型Seurat和Harmony在简单的批量校正任务中表现良好。
由于数据集的大小和数量以及场景的多样性，对复杂的集成任务进行基准测试带来了额外的挑战。最近，一项大型研究使用14个指标对5个RNA任务和两个模拟中跨集成方法类的16种方法进行基准测试 。虽然每个任务的最佳表现方法有所不同，但使用细胞类型标签的方法在不同任务中表现更好。此外，深度学习方法scANVI（带标签）、scVI和scGen（带标签）以及线性嵌入模型Scanorama表现最佳，尤其是在复杂任务上，而Harmony在不太复杂的任务上表现良好。为整合视网膜数据集以构建眼部巨型图集的特定目的而执行的类似基准也发现scVI优于其他方法。
10.1.4 选择集成方法
虽然集成方法现已进行了广泛的基准测试，但并不存在适用于所有场景的最佳方法。集成性能指标和评估pipelines（例如scIB和batchbench）可用于评估您自己的数据的集成性能。然而，许多指标（特别是那些衡量生物变异保守性的指标）需要真实的细胞身份标签。参数优化可以调整许多方法以适用于特定任务，但总的来说，可以说Harmony和Seurat在简单的批量校正任务中始终表现良好，而scVI、scGen、scANVI和Scanorama在更复杂的数据集成任务中表现良好。选择方法时，我们建议首先研究这些选项。此外，集成方法的选择可以根据所需的输出数据格式来指导。在选择一种方法之前，谨慎的做法是测试多种方法并根据成功的定量定义评估结果。
在本章的其余部分中，我们将演示一些性能最佳的方法，并快速演示如何评估集成性能。
# Python packages
import scanpy as sc
import scvi
import bbknn
import scib
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
# R interface
from rpy2.robjects import pandas2ri
from rpy2.robjects import r
import rpy2.rinterface_lib.callbacks
import anndata2ri
pandas2ri.activate()
anndata2ri.activate()
%load_ext rpy2.ipython
输出结果：
Global seed set to 0
/Users/luke.zappia/miniconda/envs/bp-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(
/Users/luke.zappia/miniconda/envs/bp-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)
/Users/luke.zappia/miniconda/envs/bp-integration/lib/python3.9/site-packages/rpy2/robjects/pandas2ri.py:263: DeprecationWarning: The global conversion available with activate() is deprecated and will be removed in the next major release. Use a local converter.
warnings.warn('The global conversion available with activate() '
/Users/luke.zappia/miniconda/envs/bp-integration/lib/python3.9/site-packages/rpy2/robjects/numpy2ri.py:205: DeprecationWarning: The global conversion available with activate() is deprecated and will be removed in the next major release. Use a local converter.
warnings.warn('The global conversion available with activate() '
%%R
# R packages
library(Seurat)
WARNING: The R package "reticulate" does not
consider that it could be called from a Python process. This
results in a quasi-obligatory segfault when rpy2 is evaluating
R code using it. On the hand, rpy2 is accounting for the
fact that it might already be running embedded in a Python
process. This is why:
- Python -> rpy2 -> R -> reticulate: crashes
- R -> reticulate -> Python -> rpy2: works
The issue with reticulate is tracked here:
https://github.com/rstudio/reticulate/issues/208
10.2 数据集
我们将用来演示数据集成的数据集包含多个骨髓单核细胞样本。这些样本最初是为2021年单细胞分析NeurIPS竞赛中的开放问题创建的。使用10x Multiome方案测量同一细胞中的RNA表达(scRNA-seq)和染色质可及性(scATAC-seq)。我们此处使用的数据版本已经过预处理，以去除低质量的细胞。
让我们使用scanpy读取数据集以获取AnnData对象。
adata_raw = sc.read_h5ad(
"../../datasets/openproblems_bmmc_multiome_genes_filtered.h5ad"
)
adata_raw.layers["logcounts"] = adata_raw.X
adata_raw
AnnData object with n_obs × n_vars = 69249 × 129921
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'
uns: 'ATAC_gene_activity_var_names', 'dataset_id', 'genome', 'organism'
obsm: 'ATAC_gene_activity', 'ATAC_lsi_full', 'ATAC_lsi_red', 'ATAC_umap', 'GEX_X_pca', 'GEX_X_umap'
layers: 'counts', 'logcounts'
完整数据集包含69,249个细胞和129,921个特征的测量值。表达矩阵有两个版本，包含原始计数值的counts和包含log标准化计数的logcounts（这些值也存储在adata.X中）。
obs槽包含多个变量，其中一些是在预处理期间计算的（用于质量控制），其他变量包含有关样本的元数据。我们感兴趣的是：
-
cell_type- 每个细胞的带注释的标签 -
batch- 每个细胞的测序批次
对于真正的分析，考虑更多变量很重要，但为了简单起见，我们在这里只考虑这些变量。
我们定义变量来保存这些名称，以便清楚我们如何在代码中使用它们。这也有助于提高可重复性，因为如果我们决定出于某种原因更改其中一个，我们可以确定它在整个笔记本中已更改。
label_key = "cell_type"
batch_key = "batch"
让我们看看不同的批次以及每个批次有多少个细胞。
adata_raw.obs[batch_key].value_counts()
s4d8 9876
s4d1 8023
s3d10 6781
s1d2 6740
s1d1 6224
s2d4 6111
s2d5 4895
s3d3 4325
s4d9 4325
s1d3 4279
s2d1 4220
s3d7 1771
s3d6 1679
Name: batch, dtype: int64
数据集中有13个不同批次。在此实验期间，从一组捐赠者身上采集了多个样本，并在不同的设施进行测序，因此此处的名称是样本编号（例如“s1”）和捐赠者（例如“d2”）的组合。为了简单起见，并减少计算时间，我们将选择三个样本来使用。
keep_batches = ["s1d3", "s2d1", "s3d7"]
adata = adata_raw[adata_raw.obs[batch_key].isin(keep_batches)].copy()
adata
AnnData object with n_obs × n_vars = 10270 × 129921
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'
uns: 'ATAC_gene_activity_var_names', 'dataset_id', 'genome', 'organism'
obsm: 'ATAC_gene_activity', 'ATAC_lsi_full', 'ATAC_lsi_red', 'ATAC_umap', 'GEX_X_pca', 'GEX_X_umap'
layers: 'counts', 'logcounts'
在进行子集化以选择这些批次后，我们剩下10,270个cells。
对于存储在var中的特征，我们有两个注释：
-
feature_types- 每个特征的类型（RNA或ATAC） -
gene_id- 与每个特征相关的基因
我们来看看特征类型。
adata.var["feature_types"].value_counts()
ATAC 116490
GEX 13431
Name: feature_types, dtype: int64
我们可以看到，ATAC特征有超过100,000个，但基因表达（“GEX”）特征只有大约13,000个。多种模态的整合是一个复杂的问题，因此现在我们将仅对基因表达特征进行子集化。我们还执行简单的过滤，以确保我们没有计数为零的特征（这是必要的，因为通过选择样本子集，我们可能已经删除了表达特定特征的所有细胞）。
adata = adata[:, adata.var["feature_types"] == "GEX"].copy()
sc.pp.filter_genes(adata, min_cells=1)
adata
AnnData object with n_obs × n_vars = 10270 × 13431
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', 'n_cells'
uns: 'ATAC_gene_activity_var_names', 'dataset_id', 'genome', 'organism'
obsm: 'ATAC_gene_activity', 'ATAC_lsi_full', 'ATAC_lsi_red', 'ATAC_umap', 'GEX_X_pca', 'GEX_X_umap'
layers: 'counts', 'logcounts'
由于子集化，我们还需要重新规范化数据。在这里，我们只是使用全局缩放通过每个细胞的总计数进行标准化。
adata.X = adata.layers["counts"].copy()
sc.pp.normalize_total(adata)
sc.pp.log1p(adata)
adata.layers["logcounts"] = adata.X.copy()
我们将使用该数据集来演示集成整合。
大多数集成方法需要一个包含所有样本和一个批次变量的单个对象（就像我们这里的那样）。相反，如果每个样本都有单独的对象，则可以使用anndata concat()函数将它们连接起来。
未整合的数据
始终建议在执行任何集成之前查看原始数据。 这可以提供一些关于任何批次效应有多大以及可能导致它们的原因的指示（以及因此应将哪些变量视为批次标签）。 对于某些实验，如果样本已经重叠，它甚至可能表明不需要积分。 例如，对于单个实验室的小鼠或细胞系研究来说，这种情况并不罕见，其中大多数影响批次效应的变量都是可以控制的（即批次校正设置）。
我们将执行高可变基因 (HVG) 选择、PCA 和 UMAP 降维，正如我们在前面的章节中看到的那样。
10.3 未整合的数据
始终建议在执行任何集成之前查看原始数据。这可以提供一些关于任何批次效应有多大以及可能导致它们的原因的指示（以及因此应将哪些变量视为批次标签）。对于某些实验，如果样本已经重叠，它甚至可能表明不需要整合。例如，对于单个实验室的小鼠或细胞系研究来说，这种情况并不罕见，其中大多数影响批次效应的变量都是可以控制的（即批次校正设置）。
我们将执行高可变基因 (HVG) 选择、PCA和UMAP降维，正如我们在前面的章节中看到的那样。
sc.pp.highly_variable_genes(adata)
sc.tl.pca(adata)
sc.pp.neighbors(adata)
sc.tl.umap(adata)
adata
AnnData object with n_obs × n_vars = 10270 × 13431
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', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
uns: 'ATAC_gene_activity_var_names', 'dataset_id', 'genome', 'organism', 'log1p', 'hvg', 'pca', 'neighbors', 'umap'
obsm: 'ATAC_gene_activity', 'ATAC_lsi_full', 'ATAC_lsi_red', 'ATAC_umap', 'GEX_X_pca', 'GEX_X_umap', 'X_pca', 'X_umap'
varm: 'PCs'
layers: 'counts', 'logcounts'
obsp: 'distances', 'connectivities'
这会向我们的AnnData对象添加几个新项目。var槽现在包括平均值、分散度和选定的变量基因。在obsp槽中，我们有KNN图的距离和连接性，在obsm中是PCA和UMAP嵌入。
让我们绘制UMAP，根据细胞身份和批次标签对点进行着色。如果数据集尚未被标记（通常是这种情况），我们将只能考虑批次标签。
adata.uns[batch_key + "_colors"] = [
"#1b9e77",
"#d95f02",
"#7570b3",
] # Set custom colours for batches
sc.pl.umap(adata, color=[label_key, batch_key], wspace=1)
输出结果：
/Users/luke.zappia/miniconda/envs/bp-integration/lib/python3.9/site-packages/scanpy/plotting/_tools/scatterplots.py:163: PendingDeprecationWarning: The get_cmap function will be deprecated in a future version. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
cmap = copy(get_cmap(cmap))
/Users/luke.zappia/miniconda/envs/bp-integration/lib/python3.9/site-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
cax = scatter(
/Users/luke.zappia/miniconda/envs/bp-integration/lib/python3.9/site-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
cax = scatter(
通常，当查看这些图时，您会发现批次之间有明显的分离。在这种情况下，我们看到的情况更加微妙，虽然来自同一标签的细胞通常彼此靠近，但批次之间存在差异。如果我们要使用这些原始数据执行聚类分析，我们可能最终会得到一些包含单个批次的聚类，这在注释阶段很难解释。我们也可能忽略稀有细胞类型，这些细胞类型在任何单个样本中都不够常见，无法产生自己的细胞簇。虽然UMAP通常可以显示批次效果，但在考虑这些2D表示时，一如既往，重要的是不要过度解释它们。对于真正的分析，您应该通过其他方式确认整合，例如通过检查标记基因的分布。
现在我们已经确认存在需要纠正的批次效应，我们可以继续使用不同的集成方法。如果批次彼此完美重叠，或者我们无需校正即可发现有意义的细胞簇，则无需执行一体化。
10.4 批量感知特征选择
如前几章所示，我们经常选择基因的子集用于分析，以减少噪音和处理时间。当我们有多个样本时，我们会做同样的事情，但是，重要的是，以批次感知的方式进行基因选择。这是因为在整个数据集中可变的基因可能会捕获批次效应，而不是我们感兴趣的生物信号。它还有助于选择与稀有细胞身份相关的基因，例如，如果身份仅存在于一个样本中，那么它的标记可能不会在所有样本中变化，但应该在该样本中。
我们可以通过在scanpy height_variable_genes()函数中设置batch_key参数来执行批量感知的高度可变基因选择。然后，scanpy将分别计算每个批次的HVG，并通过选择在最多批次中高度可变的基因来组合结果。我们在这里使用scanpy函数，因为它内置了批次感知功能。对于其他方法，我们必须在每个批次上单独运行它们，然后手动合并结果。
sc.pp.highly_variable_genes(
adata, n_top_genes=2000, flavor="cell_ranger", batch_key=batch_key
)
adata
adata.var
/Users/luke.zappia/miniconda/envs/bp-integration/lib/python3.9/site-packages/scanpy/preprocessing/_highly_variable_genes.py:478: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
hvg = hvg.append(missing_hvg, ignore_index=True)
/Users/luke.zappia/miniconda/envs/bp-integration/lib/python3.9/site-packages/scanpy/preprocessing/_highly_variable_genes.py:478: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
hvg = hvg.append(missing_hvg, ignore_index=True)
/Users/luke.zappia/miniconda/envs/bp-integration/lib/python3.9/site-packages/scanpy/preprocessing/_highly_variable_genes.py:478: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
hvg = hvg.append(missing_hvg, ignore_index=True)
feature_types gene_id n_cells highly_variable means dispersions dispersions_norm highly_variable_nbatches highly_variable_intersection
AL627309.5 GEX ENSG00000241860 112 False 0.006533 0.775289 0.357973 1 False
LINC01409 GEX ENSG00000237491 422 False 0.024462 0.716935 -0.126119 0 False
LINC01128 GEX ENSG00000228794 569 False 0.030714 0.709340 -0.296701 0 False
NOC2L GEX ENSG00000188976 675 False 0.037059 0.704363 -0.494025 0 False
KLHL17 GEX ENSG00000187961 88 False 0.005295 0.721757 -0.028456 0 False
... ... ... ... ... ... ... ... ... ...
MT-ND5 GEX ENSG00000198786 4056 False 0.269375 0.645837 -0.457491 0 False
MT-ND6 GEX ENSG00000198695 1102 False 0.051506 0.697710 -0.248421 0 False
MT-CYB GEX ENSG00000198727 5647 False 0.520368 0.613233 -0.441362 0 False
AL592183.1 GEX ENSG00000273748 732 False 0.047486 0.753417 0.413699 0 False
AC240274.1 GEX ENSG00000271254 43 False 0.001871 0.413772 -0.386635 0 False
13431 rows × 9 columns
我们可以看到var中现在多了一些额外的列：
-
height_variable_nbatches- 发现每个基因高度可变的批次数 -
height_variable_intersection- 每个基因在每批中是否高度可变 -
height_variable- 组合每批结果后是否将每个基因选择为高度可变
让我们检查一下每个基因在多少批次中是可变的：
n_batches = adata.var["highly_variable_nbatches"].value_counts()
ax = n_batches.plot(kind="bar")
n_batches
输出结果：
0 9931
1 1824
2 852
3 824
Name: highly_variable_nbatches, dtype: int64
我们注意到的第一件事是大多数基因的变异性并不高。这是通常的情况，但这可能取决于我们尝试集成的样本的差异程度。随着我们添加更多样本，重叠会减少，在所有三个批次中，相对较少的基因具有高度可变性。通过选择前2000个基因，我们选择了两批或三批中存在的所有HVG以及一批中存在的大部分HVG。
我们将创建一个仅包含选定基因的对象以用于集成。
adata_hvg = adata[:, adata.var["highly_variable"]].copy()
adata_hvg
AnnData object with n_obs × n_vars = 10270 × 2000
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', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'highly_variable_nbatches', 'highly_variable_intersection'
uns: 'ATAC_gene_activity_var_names', 'dataset_id', 'genome', 'organism', 'log1p', 'hvg', 'pca', 'neighbors', 'umap', 'batch_colors', 'cell_type_colors'
obsm: 'ATAC_gene_activity', 'ATAC_lsi_full', 'ATAC_lsi_red', 'ATAC_umap', 'GEX_X_pca', 'GEX_X_umap', 'X_pca', 'X_umap'
varm: 'PCs'
layers: 'counts', 'logcounts'
obsp: 'distances', 'connectivities'
10.5 基于变分自动编码器 (VAE) 的集成
我们将使用的第一个集成方法是scVI（single-cell Variational Inference），这是一种基于scvi-tools包中提供的条件变分自动编码器的方法。变分自动编码器是一种试图降低数据集维数的人工神经网络。条件部分是指在特定协变量（在本例中为批次）上调节此降维过程，以便协变量不会影响低维表示。在基准测试研究中，scVI已被证明在一系列数据集上表现良好，在批次校正方面取得了良好的平衡，同时保留了生物变异性。scVI直接对原始计数进行建模，因此我们为其提供计数矩阵而不是归一化表达矩阵非常重要。
首先，让我们复制数据集以用于此集成。通常没有必要这样做，但由于我们将演示多种集成方法，因此制作副本可以更轻松地显示每种方法添加的内容。
adata_scvi = adata_hvg.copy()
10.5.1 数据准备
使用scVI的第一步是准备我们的AnnData对象。此步骤存储scVI所需的一些信息，例如要使用哪个表达矩阵以及批处理密钥是什么。
scvi.model.SCVI.setup_anndata(adata_scvi, layer="counts", batch_key=batch_key)
adata_scvi
AnnData object with n_obs × n_vars = 10270 × 2000
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', '_scvi_batch', '_scvi_labels'
var: 'feature_types', 'gene_id', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'highly_variable_nbatches', 'highly_variable_intersection'
uns: 'ATAC_gene_activity_var_names', 'dataset_id', 'genome', 'organism', 'log1p', 'hvg', 'pca', 'neighbors', 'umap', 'batch_colors', 'cell_type_colors', '_scvi_uuid', '_scvi_manager_uuid'
obsm: 'ATAC_gene_activity', 'ATAC_lsi_full', 'ATAC_lsi_red', 'ATAC_umap', 'GEX_X_pca', 'GEX_X_umap', 'X_pca', 'X_umap'
varm: 'PCs'
layers: 'counts', 'logcounts'
obsp: 'distances', 'connectivities'
scVI创建的字段以_scvi为前缀。这些是专为内部使用而设计的，不应手动修改。scvi-tools作者的一般建议是在模型训练完成之前不要对我们的对象进行任何更改。在其他数据集上，您可能会看到有关包含非标准化计数数据的输入表达矩阵的警告。这通常意味着您应该检查设置函数提供的层实际上包含计数值，但如果您有对来自全长协议的数据执行基因长度校正的值，或者来自不产生整数计数的另一种量化方法的值，也会发生这种情况。
10.5.2 建立模型
我们现在可以构建一个scVI模型对象。除了我们在这里使用的scVI模型之外，scvi-tools包还包含各种其他模型（我们将在下面使用 scANVI 模型）。
model_scvi = scvi.model.SCVI(adata_scvi)
model_scvi
SCVI Model with the following params:
n_hidden: 128, n_latent: 10, n_layers: 1, dropout_rate: 0.1, dispersion: gene, gene_likelihood: zinb,
latent_distribution: normal
Training status: Not Trained
scVI模型对象包含提供的AnnData对象以及模型本身的神经网络。您可以看到当前模型尚未经过训练。如果我们想修改网络的结构，我们可以为模型构建函数提供额外的参数，但这里我们只使用默认值。
我们还可以print模型的更详细描述，向我们展示相关AnnData对象中事物的存储位置。
model_scvi.view_anndata_setup()
Anndata setup with scvi-tools version 0.18.0.
Setup via `SCVI.setup_anndata` with arguments:
{
│ 'layer': 'counts',
│ 'batch_key': 'batch',
│ 'labels_key': None,
│ 'size_factor_key': None,
│ 'categorical_covariate_keys': None,
│ 'continuous_covariate_keys': None
}
Summary Statistics
┏━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━┓
┃ Summary Stat Key ┃ Value ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━┩
│ n_batch │ 3 │
│ n_cells │ 10270 │
│ n_extra_categorical_covs │ 0 │
│ n_extra_continuous_covs │ 0 │
│ n_labels │ 1 │
│ n_vars │ 2000 │
└──────────────────────────┴───────┘
Data Registry
┏━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Registry Key ┃ scvi-tools Location ┃
┡━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│ X │ adata.layers['counts'] │
│ batch │ adata.obs['_scvi_batch'] │
│ labels │ adata.obs['_scvi_labels'] │
└──────────────┴───────────────────────────┘
batch State Registry
┏━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓
┃ Source Location ┃ Categories ┃ scvi-tools Encoding ┃
┡━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩
│ adata.obs['batch'] │ s1d3 │ 0 │
│ │ s2d1 │ 1 │
│ │ s3d7 │ 2 │
└────────────────────┴────────────┴─────────────────────┘
labels State Registry
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓
┃ Source Location ┃ Categories ┃ scvi-tools Encoding ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩
│ adata.obs['_scvi_labels'] │ 0 │ 0 │
└───────────────────────────┴────────────┴─────────────────────┘
在这里，我们可以准确地看到scVI分配了哪些信息，包括每个不同批次如何在模型中编码等详细信息。
10.5.3 训练模型
该模型将接受给定数量的epoch训练，即每个单元都通过网络的训练迭代。默认情况下，scVI使用以下启发式方法来设置epochs数。对于细胞数少于20,000个的数据集，将使用400个epoch，随着细胞数量增加到20,000以上，epoch的数量会不断减少。这背后的原因是，当网络在每个时期看到更多的单元时，它可以学到与从更多时期和更少的单元中学到的信息相同的信息量。
max_epochs_scvi = np.min([round((20000 / adata.n_obs) * 400), 400])
max_epochs_scvi
400
现在，我们针对选定的epoch数训练模型（这将需要大约20-40 分钟，具体取决于您使用的计算机）。
model_scvi.train()
GPU available: False, used: False
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
HPU available: False, using: 0 HPUs
Epoch 400/400: 100%|██████████| 400/400 [15:49<00:00, 2.27s/it, loss=648, v_num=1]
`Trainer.fit` stopped: `max_epochs=400` reached.
Epoch 400/400: 100%|██████████| 400/400 [15:49<00:00, 2.37s/it, loss=648, v_num=1]
10.5.4 提取嵌入
我们想要从训练模型中提取的主要结果是每个细胞的潜在表示。这是一种多维嵌入，其中批次效应已被移除，其使用方式与我们在分析单个数据集时使用PCA维度的方式类似。我们使用密钥X_scvi将其存储在obsm中。
adata_scvi.obsm["X_scVI"] = model_scvi.get_latent_representation()
10.5.5 计算批量校正的UMAP
我们现在将像集成之前一样可视化数据。我们计算一个新的UMAP嵌入，但我们不是在PCA空间中寻找最近邻，而是从scVI的校正表示开始。
sc.pp.neighbors(adata_scvi, use_rep="X_scVI")
sc.tl.umap(adata_scvi)
adata_scvi
AnnData object with n_obs × n_vars = 10270 × 2000
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', '_scvi_batch', '_scvi_labels'
var: 'feature_types', 'gene_id', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'highly_variable_nbatches', 'highly_variable_intersection'
uns: 'ATAC_gene_activity_var_names', 'dataset_id', 'genome', 'organism', 'log1p', 'hvg', 'pca', 'neighbors', 'umap', 'batch_colors', 'cell_type_colors', '_scvi_uuid', '_scvi_manager_uuid'
obsm: 'ATAC_gene_activity', 'ATAC_lsi_full', 'ATAC_lsi_red', 'ATAC_umap', 'GEX_X_pca', 'GEX_X_umap', 'X_pca', 'X_umap', 'X_scVI'
varm: 'PCs'
layers: 'counts', 'logcounts'
obsp: 'distances', 'connectivities'
一旦我们有了新的UMAP表示，我们就可以像以前一样按批次和身份标签将其着色。
sc.pl.umap(adata_scvi, color=[label_key, batch_key], wspace=1)
/Users/luke.zappia/miniconda/envs/bp-integration/lib/python3.9/site-packages/scanpy/plotting/_tools/scatterplots.py:163: PendingDeprecationWarning: The get_cmap function will be deprecated in a future version. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
cmap = copy(get_cmap(cmap))
/Users/luke.zappia/miniconda/envs/bp-integration/lib/python3.9/site-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
cax = scatter(
/Users/luke.zappia/miniconda/envs/bp-integration/lib/python3.9/site-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
cax = scatter(
这看起来更好！之前，不同批次的产品是相互分开的。现在，批次重叠更多，每个细胞身份标签都有一个斑点。
在许多情况下，我们还没有身份标签，因此从这个阶段开始，我们将继续进行聚类、注释和进一步分析，如其他章节中所述。
10.6 使用细胞标签进行VAE集成
当与scVI进行集成时，我们假装我们还没有任何细胞标签（尽管我们在图中显示了它们）。虽然这种情况很常见，但在某些情况下，我们确实提前了解了有关细胞身份的信息。最常见的是当我们想要将一个或多个公开可用的数据集与一项新研究的数据结合起来时。当我们至少有一些细胞的标签时，我们可以使用scANVI（使用变分推理的单细胞ANnotation）。这是scVI模型的扩展，可以合并细胞身份标签信息以及批次信息。因为它有这些额外的信息，所以它可以尝试保留细胞标签之间的差异，同时消除批次效应。基准测试表明，与scVI相比，scANVI往往能更好地保留生物信号，但有时它在消除批次效应方面并不那么有效。虽然我们在这里为所有细胞提供了标签，但也可以以半监督方式使用scANVI，其中仅为某些细胞提供标签。
我们首先创建一个scANVI模型对象。请注意，由于scANVI改进了已训练的scVI模型，因此我们提供scVI模型而不是AnnData对象。如果我们还没有训练scVI模型，我们需要首先训练。我们还为adata.obs列提供了一个键，其中包含我们的细胞标签以及对应于未标记细胞的标签。在这种情况下，我们所有的细胞都被标记，所以我们只提供一个虚拟值。在大多数情况下，检查设置是否正确非常重要，以便scANVI知道在训练期间要忽略哪个标签。
# Normally we would need to run scVI first but we have already done that here
# model_scvi = scvi.model.SCVI(adata_scvi) etc.
model_scanvi = scvi.model.SCANVI.from_scvi_model(
model_scvi, labels_key=label_key, unlabeled_category="unlabelled"
)
print(model_scanvi)
model_scanvi.view_anndata_setup()
ScanVI Model with the following params:
unlabeled_category: unlabelled, n_hidden: 128, n_latent: 10, n_layers: 1, dropout_rate: 0.1, dispersion: gene,
gene_likelihood: zinb
Training status: Not Trained
Anndata setup with scvi-tools version 0.18.0.
Setup via `SCANVI.setup_anndata` with arguments:
{
│ 'labels_key': 'cell_type',
│ 'unlabeled_category': 'unlabelled',
│ 'layer': 'counts',
│ 'batch_key': 'batch',
│ 'size_factor_key': None,
│ 'categorical_covariate_keys': None,
│ 'continuous_covariate_keys': None
}
Summary Statistics
┏━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━┓
┃ Summary Stat Key ┃ Value ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━┩
│ n_batch │ 3 │
│ n_cells │ 10270 │
│ n_extra_categorical_covs │ 0 │
│ n_extra_continuous_covs │ 0 │
│ n_labels │ 22 │
│ n_vars │ 2000 │
└──────────────────────────┴───────┘
Data Registry
┏━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Registry Key ┃ scvi-tools Location ┃
┡━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│ X │ adata.layers['counts'] │
│ batch │ adata.obs['_scvi_batch'] │
│ labels │ adata.obs['_scvi_labels'] │
└──────────────┴───────────────────────────┘
batch State Registry
┏━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓
┃ Source Location ┃ Categories ┃ scvi-tools Encoding ┃
┡━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩
│ adata.obs['batch'] │ s1d3 │ 0 │
│ │ s2d1 │ 1 │
│ │ s3d7 │ 2 │
└────────────────────┴────────────┴─────────────────────┘
labels State Registry
┏━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓
┃ Source Location ┃ Categories ┃ scvi-tools Encoding ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩
│ adata.obs['cell_type'] │ B1 B │ 0 │
│ │ CD4+ T activated │ 1 │
│ │ CD4+ T naive │ 2 │
│ │ CD8+ T │ 3 │
│ │ CD8+ T naive │ 4 │
│ │ CD14+ Mono │ 5 │
│ │ CD16+ Mono │ 6 │
│ │ Erythroblast │ 7 │
│ │ G/M prog │ 8 │
│ │ HSC │ 9 │
│ │ ILC │ 10 │
│ │ Lymph prog │ 11 │
│ │ MK/E prog │ 12 │
│ │ NK │ 13 │
│ │ Naive CD20+ B │ 14 │
│ │ Normoblast │ 15 │
│ │ Plasma cell │ 16 │
│ │ Proerythroblast │ 17 │
│ │ Transitional B │ 18 │
│ │ cDC2 │ 19 │
│ │ pDC │ 20 │
│ │ unlabelled │ 21 │
└────────────────────────┴──────────────────┴─────────────────────┘
这个scANVI模型对象与我们之前看到的scVI非常相似。如前所述，我们可以修改模型网络的结构，但这里我们只使用默认参数。
同样，我们有一个启发式的方法来选择训练时期的数量。请注意，这比以前要少得多，因为我们只是改进scVI模型，而不是从头开始训练整个网络。
max_epochs_scanvi = int(np.min([10, np.max([2, round(max_epochs_scvi / 3.0)])]))
model_scanvi.train(max_epochs=max_epochs_scanvi)
INFO Training for 10 epochs.
GPU available: False, used: False
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
HPU available: False, using: 0 HPUs
Epoch 10/10: 100%|██████████| 10/10 [00:35<00:00, 3.69s/it, loss=747, v_num=1]
`Trainer.fit` stopped: `max_epochs=10` reached.
Epoch 10/10: 100%|██████████| 10/10 [00:35<00:00, 3.56s/it, loss=747, v_num=1]
我们可以从模型中提取新的潜在表示并创建新的UMAP嵌入，就像我们对scVI所做的那样。
adata_scanvi = adata_scvi.copy()
adata_scanvi.obsm["X_scANVI"] = model_scanvi.get_latent_representation()
sc.pp.neighbors(adata_scanvi, use_rep="X_scANVI")
sc.tl.umap(adata_scanvi)
sc.pl.umap(adata_scanvi, color=[label_key, batch_key], wspace=1)
/Users/luke.zappia/miniconda/envs/bp-integration/lib/python3.9/site-packages/scanpy/plotting/_tools/scatterplots.py:163: PendingDeprecationWarning: The get_cmap function will be deprecated in a future version. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
cmap = copy(get_cmap(cmap))
/Users/luke.zappia/miniconda/envs/bp-integration/lib/python3.9/site-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
cax = scatter(
/Users/luke.zappia/miniconda/envs/bp-integration/lib/python3.9/site-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
cax = scatter(
通过查看UMAP表示，很难区分scANVI和scVI之间的差异，但正如我们将在之后（下一篇）看到的，当量化集成的质量时，度量分数存在差异。这提醒我们不应该过度解释这些二维表示，尤其是在比较方法时。
由于数据整合的内容和方法很多，将在下周继续介绍其他的数据整合方法（基于图的方法和线性嵌入整合模型），以及对整合的结果进行打分。
