10.7 基于图的整合模型Graph-based integration
我们要介绍的下一个方法是BBKNNBatch-Balanced k-Nearest Neighbor或“批量平衡KNN”。这是与scVI截然不同的方法，它不是使用神经网络将细胞嵌入批量校正的空间中，而是修改用于聚类和嵌入的k最近邻(KNN)图的构建方式。正如我们在前面的章节中所看到的，正常的KNN过程将细胞连接到整个数据集中最相似的细胞。BBKNN所做的更改是强制将细胞连接到其他批次的细胞。虽然这是一个简单的修改，但它可能非常有效，特别是当批次效应非常强时。然而，由于输出是一个集成图，它的下游用途可能有限，因为很少有包会接受它作为输入。
BBKNN的一个重要参数是每批次的邻居数量。建议的启发方法是，如果细胞数量超过100,000个，则使用25；如果细胞数量少于100,000个，则使用默认值3。
neighbors_within_batch = 25 if adata_hvg.n_obs > 100000 else 3
neighbors_within_batch
3
在使用BBKNN之前，我们首先执行PCA，就像在构建普通KNN图之前一样。与此处模拟原始计数的scVI不同，我们从对数归一化表达矩阵开始。
adata_bbknn = adata_hvg.copy()
adata_bbknn.X = adata_bbknn.layers["logcounts"].copy()
sc.pp.pca(adata_bbknn)
我们现在可以运行BBKNN，替换标准工作流程中对scanpy neighbors()函数的调用。一个重要的区别是确保设置了batch_key参数，该参数指定adata_hvg.obs中包含批次标签的列。
bbknn.bbknn(
adata_bbknn, batch_key=batch_key, neighbors_within_batch=neighbors_within_batch
)
adata_bbknn
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'
与默认的scanpy函数不同，BBKNN不允许指定存储结果的key，因此它们始终存储在默认的“neighbors”key下。
我们可以像使用普通KNN图一样使用这个新的集成图来构建UMAP嵌入。
sc.tl.umap(adata_bbknn)
sc.pl.umap(adata_bbknn, 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.8 使用相互最近邻(MNN)的线性嵌入整合
一些下游应用程序无法接受集成的嵌入或邻域图，并且需要校正的表达矩阵。可以产生此输出的一种方法是Seurat中的积分方法。Seurat积分方法属于一类线性嵌入模型，它利用相互最近邻（Seurat 称之为锚）的思想来纠正批量效应。相互最近邻是来自两个不同数据集的细胞对，当数据集放置在相同（潜在）空间中时，它们彼此相邻。找到这些细胞后，可以使用它们来对齐两个数据集并纠正它们之间的差异。在一些评估中，Seurat也被认为是最佳混合方法之一。
由于Seurat是一个 R 包，我们必须将数据从Python传输到R。这里我们准备要转换的AnnData，以便rpy2和anndata2ri可以处理它。
adata_seurat = adata_hvg.copy()
# Convert categorical columns to strings
adata_seurat.obs[batch_key] = adata_seurat.obs[batch_key].astype(str)
adata_seurat.obs[label_key] = adata_seurat.obs[label_key].astype(str)
# Delete uns as this can contain arbitrary objects which are difficult to convert
del adata_seurat.uns
adata_seurat
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'
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'
基于anndata2ri，准备好的AnnData现在可以在R中作为SingleCellExperiment对象使用。请注意，与AnnData对象相比，这是转置的，因此我们的观察结果（细胞）现在是列，我们的变量（基因）现在是行。
%%R -i adata_seurat
adata_seurat
class: SingleCellExperiment
dim: 2000 10270
metadata(0):
assays(3): X counts logcounts
rownames(2000): GPR153 TNFRSF25 ... TMLHE-AS1 MT-ND3
rowData names(9): feature_types gene_id ... highly_variable_nbatches
highly_variable_intersection
colnames(10270): TCACCTGGTTAGGTTG-3-s1d3 CGTTAACAGGTGTCCA-3-s1d3 ...
AGCAGGTAGGCTATGT-12-s3d7 GCCATGATCCCTTGCG-12-s3d7
colData names(28): GEX_pct_counts_mt GEX_n_counts ... QCMeds
DonorSmoker
reducedDimNames(8): ATAC_gene_activity ATAC_lsi_full ... PCA UMAP
mainExpName: NULL
altExpNames(0):
Seurat使用自己的对象来存储数据。作者提供了一个从SingleCellExperiment转换的函数，很有帮助。我们只需提供SingleCellExperiment对象并告诉Seurat哪些分析（AnnData对象中的层）包含原始计数和标准化表达（Seurat 存储在名为“数据”的槽中）。
%%R -i adata_seurat
seurat <- as.Seurat(adata_seurat, counts = "counts", data = "logcounts")
seurat
An object of class Seurat
2000 features across 10270 samples within 1 assay
Active assay: originalexp (2000 features, 0 variable features)
8 dimensional reductions calculated: ATAC_gene_activity, ATAC_lsi_full, ATAC_lsi_red, ATAC_umap, GEX_X_pca, GEX_X_umap, PCA, UMAP
与我们所看到的采用单个对象和批处理密钥的其他一些方法不同，Seurat集成函数需要一个对象列表。我们使用SplitObject()函数创建它。
%%R -i batch_key
batch_list <- SplitObject(seurat, split.by = batch_key)
batch_list
输出结果：
$s1d3
An object of class Seurat
2000 features across 4279 samples within 1 assay
Active assay: originalexp (2000 features, 0 variable features)
8 dimensional reductions calculated: ATAC_gene_activity, ATAC_lsi_full, ATAC_lsi_red, ATAC_umap, GEX_X_pca, GEX_X_umap, PCA, UMAP
$s2d1
An object of class Seurat
2000 features across 4220 samples within 1 assay
Active assay: originalexp (2000 features, 0 variable features)
8 dimensional reductions calculated: ATAC_gene_activity, ATAC_lsi_full, ATAC_lsi_red, ATAC_umap, GEX_X_pca, GEX_X_umap, PCA, UMAP
$s3d7
An object of class Seurat
2000 features across 1771 samples within 1 assay
Active assay: originalexp (2000 features, 0 variable features)
8 dimensional reductions calculated: ATAC_gene_activity, ATAC_lsi_full, ATAC_lsi_red, ATAC_umap, GEX_X_pca, GEX_X_umap, PCA, UMAP
我们现在可以使用此列表来查找每对数据集的锚点。通常，您会首先识别批量感知的高度可变基因（使用FindVariableFeatures()和SelectIntegrationFeatures()函数），但正如我们已经完成的那样，我们告诉Seurat使用对象中的所有特征。
%%R
anchors <- FindIntegrationAnchors(batch_list, anchor.features = rownames(seurat))
anchors
| | 0 % ~calculating |+++++++++++++++++ | 33% ~03s |++++++++++++++++++++++++++++++++++ | 67% ~01s |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02s
| | 0 % ~calculating |+++++++++++++++++ | 33% ~01m 39s |++++++++++++++++++++++++++++++++++ | 67% ~41s |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01m 56s
An AnchorSet object containing 25352 anchors between 3 Seurat objects
This can be used as input to IntegrateData.
然后，Seurat可以使用锚点来计算将一个数据集映射到另一个数据集的转换。这是以成对的方式完成的，直到合并所有数据集。默认情况下，Seurat将确定合并顺序，以便首先将更多相似的数据集合并在一起，但也可以定义此顺序。
%%R
integrated <- IntegrateData(anchors)
integrated
An object of class Seurat
4000 features across 10270 samples within 2 assays
Active assay: integrated (2000 features, 2000 variable features)
1 other assay present: originalexp
结果是另一个Seurat对象，但现在请注意，主动测定被称为“集成”。这包含校正后的表达矩阵，它是积分的最终输出。
在这里，我们提取该矩阵并准备将其传输回Python。
%%R -o integrated_expr
# Extract the integrated expression matrix
integrated_expr <- GetAssayData(integrated)
# Make sure the rows and columns are in the same order as the original object
integrated_expr <- integrated_expr[rownames(seurat), colnames(seurat)]
# Transpose the matrix to AnnData format
integrated_expr <- t(integrated_expr)
print(integrated_expr[1:10, 1:10])
输出结果：
10 x 10 sparse Matrix of class "dgCMatrix"
TCACCTGGTTAGGTTG-3-s1d3 . -0.0005365199 1.032812e-02 -2.653187e-02
CGTTAACAGGTGTCCA-3-s1d3 0.0001382038 -0.1809919733 -1.454901e-02 3.608087e-03
ATTCGTTTCAGTATTG-3-s1d3 -0.0121073040 -0.0634131457 . 2.144075e-02
GGACCGAAGTGAGGTA-3-s1d3 . . 2.972292e-04 .
ATGAAGCCAGGGAGCT-3-s1d3 -0.0139047071 -0.0313151274 . 2.239855e-02
AGTGCGGAGTAAGGGC-3-s1d3 -0.0004299227 -0.0002657828 . -1.871410e-03
CTACCTCAGACACCGC-3-s1d3 -0.0055208620 -0.0398862297 7.182253e-06 8.240406e-03
CTTCAATTCACGAATC-3-s1d3 . -0.0109928448 . 1.935677e-04
CCATTGTGTAGACAAA-3-s1d3 . 0.0171909615 . 5.711311e-05
CCGTTACTCAATGTGC-3-s1d3 0.0139905515 0.0007981117 2.303346e-03 1.356206e-02
TCACCTGGTTAGGTTG-3-s1d3 -0.023237586 0.031938502 -0.003196878 0.01777767
CGTTAACAGGTGTCCA-3-s1d3 0.114149766 -0.013183394 0.038076743 0.80491293
ATTCGTTTCAGTATTG-3-s1d3 -0.054419895 0.010955783 -0.005951634 0.37223306
GGACCGAAGTGAGGTA-3-s1d3 0.002305526 0.011544715 0.011133475 0.02366670
ATGAAGCCAGGGAGCT-3-s1d3 -0.123505733 -0.009382408 0.002153635 -0.07013584
AGTGCGGAGTAAGGGC-3-s1d3 0.035848768 0.013858992 -0.000379393 0.05617137
CTACCTCAGACACCGC-3-s1d3 -0.003837952 0.082027580 -0.001109391 -0.06307772
CTTCAATTCACGAATC-3-s1d3 0.052970708 0.153601547 0.247920321 -0.01143158
CCATTGTGTAGACAAA-3-s1d3 -0.015445186 0.025763467 -0.003632830 0.02040172
CCGTTACTCAATGTGC-3-s1d3 -0.018025385 0.022560136 0.005755797 0.61496228
TCACCTGGTTAGGTTG-3-s1d3 -6.661643e-03 -0.0183198213
CGTTAACAGGTGTCCA-3-s1d3 5.079864e-02 0.0394717101
ATTCGTTTCAGTATTG-3-s1d3 6.600434e-02 0.0009021682
GGACCGAAGTGAGGTA-3-s1d3 7.172740e-04 0.0095352521
ATGAAGCCAGGGAGCT-3-s1d3 1.226039e-01 0.0063816143
AGTGCGGAGTAAGGGC-3-s1d3 -1.830964e-03 0.0008381943
CTACCTCAGACACCGC-3-s1d3 8.559358e-01 0.0102285085
CTTCAATTCACGAATC-3-s1d3 -5.318167e-06 -0.0341661907
CCATTGTGTAGACAAA-3-s1d3 -6.362298e-03 0.0232357011
CCGTTACTCAATGTGC-3-s1d3 -4.910886e-02 0.0317359576
现在，我们将校正后的表达矩阵作为图层存储在AnnData对象中。 我们还设置adata.X使用这个矩阵。
adata_seurat.X = integrated_expr
adata_seurat.layers["seurat"] = integrated_expr
print(adata_seurat)
adata.X
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'
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', 'seurat'
obsp: 'distances', 'connectivities'
<10270x13431 sparse matrix of type '<class 'numpy.float32'>'
with 14348115 stored elements in Compressed Sparse Row format>
现在我们有了积分的结果，我们可以计算UMAP并将其绘制出来，就像我们对其他方法所做的那样（我们也可以在R中完成此操作）。
# Reset the batch colours because we deleted them earlier
adata_seurat.uns[batch_key + "_colors"] = [
"#1b9e77",
"#d95f02",
"#7570b3",
]
sc.tl.pca(adata_seurat)
sc.pp.neighbors(adata_seurat)
sc.tl.umap(adata_seurat)
sc.pl.umap(adata_seurat, color=[label_key, batch_key], wspace=1)
正如我们之前所看到的，批次是混合的，而标签是分开的。人们很容易选择基于UMAP的集成，但这并不能完全代表数据整合的质量。在后面，我们将提出一些更严格评估集成整合方法的方法。
10.9 对集成整合进行基准测试
这里演示的方法是根据基准测试实验的结果选择的，包括单细胞集成基准测试项目。该项目还制作了一个名为scib的软件包，可用于运行一系列集成方法以及用于评估的指标。在本节中，我们将展示如何使用此包来评估集成的质量。
scib指标可以单独运行，但也有用于一次运行多个指标的包装器。在这里，我们运行指标的子集，使用metrics_fast()函数可以快速计算这些指标。该函数需要几个参数：原始未集成数据集、集成数据集、批次键和标签键。根据集成方法的输出，我们可能还需要提供其他参数，例如，这里我们使用embed参数指定用于scVI和scANVI的嵌入。您还可以使用其他参数控制某些指标的运行方式。另请注意，您可能需要检查对象的格式是否正确，以便scIB可以找到所需的信息。
让我们为上面执行的每个集成以及未集成的数据（在高度可变的基因选择之后）运行指标。
metrics_scvi = scib.metrics.metrics_fast(
adata, adata_scvi, batch_key, label_key, embed="X_scVI"
)
metrics_scanvi = scib.metrics.metrics_fast(
adata, adata_scanvi, batch_key, label_key, embed="X_scANVI"
)
metrics_bbknn = scib.metrics.metrics_fast(adata, adata_bbknn, batch_key, label_key)
metrics_seurat = scib.metrics.metrics_fast(adata, adata_seurat, batch_key, label_key)
metrics_hvg = scib.metrics.metrics_fast(adata, adata_hvg, batch_key, label_key)
Silhouette score...
PC regression...
Isolated labels ASW...
Graph connectivity...
Silhouette score...
PC regression...
Isolated labels ASW...
Graph connectivity...
Silhouette score...
PC regression...
Isolated labels ASW...
Graph connectivity...
Silhouette score...
PC regression...
Isolated labels ASW...
Graph connectivity...
Silhouette score...
PC regression...
Isolated labels ASW...
Graph connectivity...
以下是单个集成整合的指标结果之一的示例：
metrics_hvg
0
NMI_cluster/label NaN
ARI_cluster/label NaN
ASW_label 0.555775
ASW_label/batch 0.833656
PCR_batch 0.537850
cell_cycle_conservation NaN
isolated_label_F1 NaN
isolated_label_silhouette 0.487975
graph_conn 0.985533
kBET NaN
iLISI NaN
cLISI NaN
hvg_overlap 1.000000
trajectory NaN
每行都是一个不同的指标，值显示该指标的分数。分数介于0和1之间，其中1表示性能良好，0表示性能较差（如果需要，scib还可以返回某些指标的未缩放分数）。因为我们在这里只运行快速指标，所以一些指标具有NaN分数。另请注意，某些指标无法与某些输出格式一起使用，这也可能是返回NaN 值的原因。
为了比较这些方法，将所有指标结果放在一张表中非常有用。该代码将它们组合起来并将它们整理成更方便的格式。
# Concatenate metrics results
metrics = pd.concat(
[metrics_scvi, metrics_scanvi, metrics_bbknn, metrics_seurat, metrics_hvg],
axis="columns",
)
# Set methods as column names
metrics = metrics.set_axis(
["scVI", "scANVI", "BBKNN", "Seurat", "Unintegrated"], axis="columns"
)
# Select only the fast metrics
metrics = metrics.loc[
[
"ASW_label",
"ASW_label/batch",
"PCR_batch",
"isolated_label_silhouette",
"graph_conn",
"hvg_overlap",
],
:,
]
# Transpose so that metrics are columns and methods are rows
metrics = metrics.T
# Remove the HVG overlap metric because it's not relevant to embedding outputs
metrics = metrics.drop(columns=["hvg_overlap"])
metrics
ASW_label ASW_label/batch PCR_batch isolated_label_silhouette graph_conn
scVI 0.570953 0.907165 0.913474 0.560875 0.982940
scANVI 0.587934 0.905402 0.887735 0.556232 0.983729
BBKNN 0.555469 0.847772 0.537850 0.481913 0.969514
Seurat 0.571231 0.915237 0.788894 0.488585 0.989366
Unintegrated 0.555775 0.833656 0.537850 0.487975 0.985533
现在，我们将所有分数放在一张表中，其中指标为列，方法为行。
metrics.style.background_gradient(cmap="Blues")
ASW_label ASW_label/batch PCR_batch isolated_label_silhouette graph_conn
scVI 0.570953 0.907165 0.913474 0.560875 0.982940
scANVI 0.587934 0.905402 0.887735 0.556232 0.983729
BBKNN 0.555469 0.847772 0.537850 0.481913 0.969514
Seurat 0.571231 0.915237 0.788894 0.488585 0.989366
Unintegrated 0.555775 0.833656 0.537850 0.487975 0.985533
对于某些指标，分数往往处于相对较小的范围内。为了强调方法之间的差异并将每个指标放在相同的等级上，我们对它们进行了缩放，以便表现最差的人得到0分，表现最好的人得到1分，其他人的分数介于两者之间。
metrics_scaled = (metrics - metrics.min()) / (metrics.max() - metrics.min())
metrics_scaled.style.background_gradient(cmap="Blues")
ASW_label ASW_label/batch PCR_batch isolated_label_silhouette graph_conn
scVI 0.476942 0.901060 1.000000 1.000000 0.676270
scANVI 1.000000 0.879454 0.931475 0.941194 0.716018
BBKNN 0.000000 0.173036 0.000000 0.000000 0.000000
Seurat 0.485521 1.000000 0.668338 0.084500 1.000000
Unintegrated 0.009437 0.000000 0.000000 0.076776 0.806917
这些值现在可以更好地表示方法之间的差异。然而，值得注意的是，缩放分数只能用于比较这组特定集成的相对性能。如果我们想添加另一种方法，我们需要再次执行缩放。我们也不能说集成绝对是“好”，只是说它比我们尝试过的其他方法更好。这种缩放强调了方法之间的差异。例如，如果我们的指标分数为0.92、0.94和0.96，则这些分数将缩放为0、0.5和1.0。这使得第一种方法的得分看起来要差得多，尽管它只比其他两种方法低一点点，但仍然得到了很高的分数。当比较几种方法并且当它们获得相似的原始分数时，这种效果会更大。您查看原始分数还是缩放分数取决于您是想关注绝对性能还是方法之间的性能差异。
评估指标可以分为两类，一类衡量批次效应的消除，另一类衡量生物变异的保存。我们可以通过取每组缩放值的平均值来计算每个类别的汇总分数。这种汇总分数对于原始值没有意义，因为某些指标始终比其他指标产生更高的分数（因此对平均值的影响更大）。
metrics_scaled["Batch"] = metrics_scaled[
["ASW_label/batch", "PCR_batch", "graph_conn"]
].mean(axis=1)
metrics_scaled["Bio"] = metrics_scaled[["ASW_label", "isolated_label_silhouette"]].mean(
axis=1
)
metrics_scaled.style.background_gradient(cmap="Blues")
ASW_label ASW_label/batch PCR_batch isolated_label_silhouette graph_conn Batch Bio
scVI 0.476942 0.901060 1.000000 1.000000 0.676270 0.859110 0.738471
scANVI 1.000000 0.879454 0.931475 0.941194 0.716018 0.842316 0.970597
BBKNN 0.000000 0.173036 0.000000 0.000000 0.000000 0.057679 0.000000
Seurat 0.485521 1.000000 0.668338 0.084500 1.000000 0.889446 0.285010
Unintegrated 0.009437 0.000000 0.000000 0.076776 0.806917 0.268972 0.043106
将两个汇总分数相互绘制可以表明每种方法的优先级。有些会偏向批量校正，而另一些则倾向于保留生物变异。
fig, ax = plt.subplots()
ax.set_xlim(0, 1)
ax.set_ylim(0, 1)
metrics_scaled.plot.scatter(
x="Batch",
y="Bio",
c=range(len(metrics_scaled)),
ax=ax,
)
for k, v in metrics_scaled[["Batch", "Bio"]].iterrows():
ax.annotate(
k,
v,
xytext=(6, -3),
textcoords="offset points",
family="sans-serif",
fontsize=12,
)
在我们的小示例场景中，BBKNN显然是表现最差的，在批次去除和生物保护方面得分最低。其他三种方法具有相似的批次校正分数，其中scANVI在生物保护方面得分最高，其次是scVI和Seurat。
为了获得每种方法的总体分数，我们可以将两个汇总分数结合起来。scIB论文建议权重为40%批量校正和 60% 生物保护，但您可能更愿意根据数据集的优先级对事物进行不同的权重。
metrics_scaled["Overall"] = 0.4 * metrics_scaled["Batch"] + 0.6 * metrics_scaled["Bio"]
metrics_scaled.style.background_gradient(cmap="Blues")
ASW_label ASW_label/batch PCR_batch isolated_label_silhouette graph_conn Batch Bio Overall
scVI 0.476942 0.901060 1.000000 1.000000 0.676270 0.859110 0.738471 0.786726
scANVI 1.000000 0.879454 0.931475 0.941194 0.716018 0.842316 0.970597 0.919285
BBKNN 0.000000 0.173036 0.000000 0.000000 0.000000 0.057679 0.000000 0.023071
Seurat 0.485521 1.000000 0.668338 0.084500 1.000000 0.889446 0.285010 0.526785
Unintegrated 0.009437 0.000000 0.000000 0.076776 0.806917 0.268972 0.043106 0.133453
让我们制作一个快速条形图来可视化整体性能。
metrics_scaled.plot.bar(y="Overall")
正如我们已经看到的，scVI和scANVI是表现最好的，其中scANVI得分略高。值得注意的是，这只是如何针对该特定数据集运行这些指标的示例，而不是对这些方法的正确评估。为此，您应该参考现有的基准测试论文。特别是，我们只运行了一小部分高性能方法和一部分指标。另外，请记住，分数与所使用的方法相关，因此即使这些方法的性能几乎相同，微小的差异也会被夸大。
现有的基准测试建议的方法通常表现良好，但不同场景下的性能也可能存在很大差异。对于某些分析，您自己进行集成评估可能是更好的选择。scib包使这个过程变得更容易，但它仍然是一项艰巨的任务，依赖于对基本事实的充分了解和对指标的解释。
10.10 Session information
10.10.1 Python
import session_info
session_info.show()
-----
anndata 0.8.0
anndata2ri 1.1
bbknn NA
matplotlib 3.6.1
numpy 1.23.4
pandas 1.5.1
rpy2 3.5.1
scanpy 1.9.1
scib 1.0.4
scipy 1.9.3
scvi 0.18.0
session_info 1.0.0
-----
-----
IPython 8.5.0
jupyter_client 7.4.4
jupyter_core 4.11.1
jupyterlab 3.5.0
notebook 6.5.1
-----
Python 3.9.13 | packaged by conda-forge | (main, May 27 2022, 17:01:00) [Clang 13.0.1 ]
macOS-10.15.7-x86_64-i386-64bit
-----
Session information updated at 2022-11-10 12:53
10.10.2 R
─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.1.3 (2022-03-10)
os macOS Catalina 10.15.7
system x86_64, darwin13.4.0
ui unknown
language (EN)
collate C
ctype UTF-8
tz Europe/Berlin
date 2022-11-10
pandoc 2.19.2 @ /Users/luke.zappia/miniconda/envs/bp-integration/bin/pandoc
─ Packages ───────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
abind 1.4-5 2016-07-21 [1] CRAN (R 4.1.3)
Biobase * 2.54.0 2021-10-26 [1] Bioconductor
BiocGenerics * 0.40.0 2021-10-26 [1] Bioconductor
bitops 1.0-7 2021-04-24 [1] CRAN (R 4.1.3)
cli 3.4.1 2022-09-23 [1] CRAN (R 4.1.3)
cluster 2.1.4 2022-08-22 [1] CRAN (R 4.1.3)
codetools 0.2-18 2020-11-04 [1] CRAN (R 4.1.3)
colorspace 2.0-3 2022-02-21 [1] CRAN (R 4.1.3)
cowplot 1.1.1 2020-12-30 [1] CRAN (R 4.1.3)
data.table 1.14.4 2022-10-17 [1] CRAN (R 4.1.3)
DelayedArray 0.20.0 2021-10-26 [1] Bioconductor
deldir 1.0-6 2021-10-23 [1] CRAN (R 4.1.3)
digest 0.6.30 2022-10-18 [1] CRAN (R 4.1.3)
dplyr 1.0.10 2022-09-01 [1] CRAN (R 4.1.3)
ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.1.3)
fansi 1.0.3 2022-03-24 [1] CRAN (R 4.1.3)
fastmap 1.1.0 2021-01-25 [1] CRAN (R 4.1.3)
fitdistrplus 1.1-8 2022-03-10 [1] CRAN (R 4.1.3)
future 1.28.0 2022-09-02 [1] CRAN (R 4.1.3)
future.apply 1.9.1 2022-09-07 [1] CRAN (R 4.1.3)
generics 0.1.3 2022-07-05 [1] CRAN (R 4.1.3)
GenomeInfoDb * 1.30.1 2022-01-30 [1] Bioconductor
GenomeInfoDbData 1.2.7 2022-10-28 [1] Bioconductor
GenomicRanges * 1.46.1 2021-11-18 [1] Bioconductor
ggplot2 3.3.6 2022-05-03 [1] CRAN (R 4.1.3)
ggrepel 0.9.1 2021-01-15 [1] CRAN (R 4.1.3)
ggridges 0.5.4 2022-09-26 [1] CRAN (R 4.1.3)
globals 0.16.1 2022-08-28 [1] CRAN (R 4.1.3)
glue 1.6.2 2022-02-24 [1] CRAN (R 4.1.3)
goftest 1.2-3 2021-10-07 [1] CRAN (R 4.1.3)
gridExtra 2.3 2017-09-09 [1] CRAN (R 4.1.3)
gtable 0.3.1 2022-09-01 [1] CRAN (R 4.1.3)
htmltools 0.5.3 2022-07-18 [1] CRAN (R 4.1.3)
htmlwidgets 1.5.4 2021-09-08 [1] CRAN (R 4.1.3)
httpuv 1.6.6 2022-09-08 [1] CRAN (R 4.1.3)
httr 1.4.4 2022-08-17 [1] CRAN (R 4.1.3)
ica 1.0-3 2022-07-08 [1] CRAN (R 4.1.3)
igraph 1.3.5 2022-09-22 [1] CRAN (R 4.1.3)
IRanges * 2.28.0 2021-10-26 [1] Bioconductor
irlba 2.3.5.1 2022-10-03 [1] CRAN (R 4.1.3)
jsonlite 1.8.3 2022-10-21 [1] CRAN (R 4.1.3)
KernSmooth 2.23-20 2021-05-03 [1] CRAN (R 4.1.3)
later 1.2.0 2021-04-23 [1] CRAN (R 4.1.3)
lattice 0.20-45 2021-09-22 [1] CRAN (R 4.1.3)
lazyeval 0.2.2 2019-03-15 [1] CRAN (R 4.1.3)
leiden 0.4.3 2022-09-10 [1] CRAN (R 4.1.3)
lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.1.3)
listenv 0.8.0 2019-12-05 [1] CRAN (R 4.1.3)
lmtest 0.9-40 2022-03-21 [1] CRAN (R 4.1.3)
magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.1.3)
MASS 7.3-58.1 2022-08-03 [1] CRAN (R 4.1.3)
Matrix * 1.5-1 2022-09-13 [1] CRAN (R 4.1.3)
MatrixGenerics * 1.6.0 2021-10-26 [1] Bioconductor
matrixStats * 0.62.0 2022-04-19 [1] CRAN (R 4.1.3)
mgcv 1.8-41 2022-10-21 [1] CRAN (R 4.1.3)
mime 0.12 2021-09-28 [1] CRAN (R 4.1.3)
miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.1.3)
munsell 0.5.0 2018-06-12 [1] CRAN (R 4.1.3)
nlme 3.1-160 2022-10-10 [1] CRAN (R 4.1.3)
parallelly 1.32.1 2022-07-21 [1] CRAN (R 4.1.3)
patchwork 1.1.2 2022-08-19 [1] CRAN (R 4.1.3)
pbapply 1.5-0 2021-09-16 [1] CRAN (R 4.1.3)
pillar 1.8.1 2022-08-19 [1] CRAN (R 4.1.3)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.1.3)
plotly 4.10.0 2021-10-09 [1] CRAN (R 4.1.3)
plyr 1.8.7 2022-03-24 [1] CRAN (R 4.1.3)
png 0.1-7 2013-12-03 [1] CRAN (R 4.1.3)
polyclip 1.10-4 2022-10-20 [1] CRAN (R 4.1.3)
progressr 0.11.0 2022-09-02 [1] CRAN (R 4.1.3)
promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.1.3)
purrr 0.3.5 2022-10-06 [1] CRAN (R 4.1.3)
R6 2.5.1 2021-08-19 [1] CRAN (R 4.1.3)
RANN 2.6.1 2019-01-08 [1] CRAN (R 4.1.3)
RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.1.3)
Rcpp 1.0.9 2022-07-08 [1] CRAN (R 4.1.3)
RcppAnnoy 0.0.19 2021-07-30 [1] CRAN (R 4.1.3)
RCurl 1.98-1.9 2022-10-03 [1] CRAN (R 4.1.3)
reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.1.3)
reticulate 1.26 2022-08-31 [1] CRAN (R 4.1.3)
rgeos 0.5-9 2021-12-15 [1] CRAN (R 4.1.3)
rlang 1.0.6 2022-09-24 [1] CRAN (R 4.1.3)
ROCR 1.0-11 2020-05-02 [1] CRAN (R 4.1.3)
rpart 4.1.19 2022-10-21 [1] CRAN (R 4.1.3)
Rtsne 0.16 2022-04-17 [1] CRAN (R 4.1.3)
S4Vectors * 0.32.4 2022-03-24 [1] Bioconductor
scales 1.2.1 2022-08-20 [1] CRAN (R 4.1.3)
scattermore 0.8 2022-02-14 [1] CRAN (R 4.1.3)
sctransform 0.3.5 2022-09-21 [1] CRAN (R 4.1.3)
sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.1.3)
Seurat * 4.2.0 2022-09-21 [1] CRAN (R 4.1.3)
SeuratObject * 4.1.2 2022-09-20 [1] CRAN (R 4.1.3)
shiny 1.7.3 2022-10-25 [1] CRAN (R 4.1.3)
SingleCellExperiment * 1.16.0 2021-10-26 [1] Bioconductor
sp * 1.5-0 2022-06-05 [1] CRAN (R 4.1.3)
spatstat.core 2.4-4 2022-05-18 [1] CRAN (R 4.1.3)
spatstat.data 3.0-0 2022-10-21 [1] CRAN (R 4.1.3)
spatstat.geom 3.0-3 2022-10-25 [1] CRAN (R 4.1.3)
spatstat.random 2.2-0 2022-03-30 [1] CRAN (R 4.1.3)
spatstat.sparse 3.0-0 2022-10-21 [1] CRAN (R 4.1.3)
spatstat.utils 3.0-1 2022-10-19 [1] CRAN (R 4.1.3)
stringi 1.7.8 2022-07-11 [1] CRAN (R 4.1.3)
stringr 1.4.1 2022-08-20 [1] CRAN (R 4.1.3)
SummarizedExperiment * 1.24.0 2021-10-26 [1] Bioconductor
survival 3.4-0 2022-08-09 [1] CRAN (R 4.1.3)
tensor 1.5 2012-05-05 [1] CRAN (R 4.1.3)
tibble 3.1.8 2022-07-22 [1] CRAN (R 4.1.3)
tidyr 1.2.1 2022-09-08 [1] CRAN (R 4.1.3)
tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.1.3)
utf8 1.2.2 2021-07-24 [1] CRAN (R 4.1.3)
uwot 0.1.14 2022-08-22 [1] CRAN (R 4.1.3)
vctrs 0.5.0 2022-10-22 [1] CRAN (R 4.1.3)
viridisLite 0.4.1 2022-08-22 [1] CRAN (R 4.1.3)
xtable 1.8-4 2019-04-21 [1] CRAN (R 4.1.3)
XVector 0.34.0 2021-10-26 [1] Bioconductor
zlibbioc 1.40.0 2021-10-26 [1] Bioconductor
zoo 1.8-11 2022-09-17 [1] CRAN (R 4.1.3)
