單細(xì)胞RNA-seq生信分析全流程——第十篇：數(shù)據(jù)整合（二）

10.7 基于圖的整合模型Graph-based integration

我們要介紹的下一個方法是BBKNNBatch-Balanced k-Nearest Neighbor或“批量平衡KNN”匾乓。這是與scVI截然不同的方法脸候，它不是使用神經(jīng)網(wǎng)絡(luò)將細(xì)胞嵌入批量校正的空間中秦效，而是修改用于聚類和嵌入的k最近鄰(KNN)圖的構(gòu)建方式捷凄。正如我們在前面的章節(jié)中所看到的厘擂，正常的KNN過程將細(xì)胞連接到整個數(shù)據(jù)集中最相似的細(xì)胞程腹。BBKNN所做的更改是強(qiáng)制將細(xì)胞連接到其他批次的細(xì)胞。雖然這是一個簡單的修改舅逸，但它可能非常有效桌肴，特別是當(dāng)批次效應(yīng)非常強(qiáng)時。然而琉历，由于輸出是一個集成圖坠七，它的下游用途可能有限，因為很少有包會接受它作為輸入旗笔。

BBKNN的一個重要參數(shù)是每批次的鄰居數(shù)量彪置。建議的啟發(fā)方法是，如果細(xì)胞數(shù)量超過100,000個蝇恶，則使用25拳魁；如果細(xì)胞數(shù)量少于100,000個，則使用默認(rèn)值3撮弧。

neighbors_within_batch = 25 if adata_hvg.n_obs > 100000 else 3
neighbors_within_batch

在使用BBKNN之前潘懊，我們首先執(zhí)行PCA，就像在構(gòu)建普通KNN圖之前一樣贿衍。與此處模擬原始計數(shù)的scVI不同授舟，我們從對數(shù)歸一化表達(dá)矩陣開始。

adata_bbknn = adata_hvg.copy()
adata_bbknn.X = adata_bbknn.layers["logcounts"].copy()
sc.pp.pca(adata_bbknn)

我們現(xiàn)在可以運行BBKNN贸辈，替換標(biāo)準(zhǔn)工作流程中對scanpy neighbors()函數(shù)的調(diào)用释树。一個重要的區(qū)別是確保設(shè)置了batch_key參數(shù)，該參數(shù)指定adata_hvg.obs中包含批次標(biāo)簽的列。

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'

與默認(rèn)的scanpy函數(shù)不同躏哩，BBKNN不允許指定存儲結(jié)果的key署浩，因此它們始終存儲在默認(rèn)的“neighbors”key下。

我們可以像使用普通KNN圖一樣使用這個新的集成圖來構(gòu)建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(

與將細(xì)胞身份分組在一起的未集成數(shù)據(jù)相比，這種集成也得到了改進(jìn)炊汤，但我們?nèi)匀豢吹脚沃g存在一些變化正驻。

10.8 使用相互最近鄰(MNN)的線性嵌入整合

一些下游應(yīng)用程序無法接受集成的嵌入或鄰域圖，并且需要校正的表達(dá)矩陣抢腐」檬铮可以產(chǎn)生此輸出的一種方法是Seurat中的積分方法。Seurat積分方法屬于一類線性嵌入模型迈倍，它利用相互最近鄰（Seurat 稱之為錨）的思想來糾正批量效應(yīng)伤靠。相互最近鄰是來自兩個不同數(shù)據(jù)集的細(xì)胞對，當(dāng)數(shù)據(jù)集放置在相同（潛在）空間中時啼染，它們彼此相鄰宴合。找到這些細(xì)胞后，可以使用它們來對齊兩個數(shù)據(jù)集并糾正它們之間的差異迹鹅。在一些評估中卦洽，Seurat也被認(rèn)為是最佳混合方法之一。
由于Seurat是一個 R 包斜棚，我們必須將數(shù)據(jù)從Python傳輸?shù)絉阀蒂。這里我們準(zhǔn)備要轉(zhuǎn)換的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蚤霞，準(zhǔn)備好的AnnData現(xiàn)在可以在R中作為SingleCellExperiment對象使用。請注意义钉，與AnnData對象相比昧绣，這是轉(zhuǎn)置的，因此我們的觀察結(jié)果（細(xì)胞）現(xiàn)在是列断医，我們的變量（基因）現(xiàn)在是行滞乙。

%%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使用自己的對象來存儲數(shù)據(jù)。作者提供了一個從SingleCellExperiment轉(zhuǎn)換的函數(shù)鉴嗤，很有幫助斩启。我們只需提供SingleCellExperiment對象并告訴Seurat哪些分析（AnnData對象中的層）包含原始計數(shù)和標(biāo)準(zhǔn)化表達(dá)（Seurat 存儲在名為“數(shù)據(jù)”的槽中）。

%%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集成函數(shù)需要一個對象列表兔簇。我們使用SplitObject()函數(shù)創(chuàng)建它。

%%R -i batch_key
batch_list <- SplitObject(seurat, split.by = batch_key)
batch_list

輸出結(jié)果：

$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

我們現(xiàn)在可以使用此列表來查找每對數(shù)據(jù)集的錨點。通常垄琐，您會首先識別批量感知的高度可變基因（使用FindVariableFeatures()和SelectIntegrationFeatures()函數(shù)）边酒，但正如我們已經(jīng)完成的那樣，我們告訴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可以使用錨點來計算將一個數(shù)據(jù)集映射到另一個數(shù)據(jù)集的轉(zhuǎn)換。這是以成對的方式完成的翻擒，直到合并所有數(shù)據(jù)集氓涣。默認(rèn)情況下，Seurat將確定合并順序陋气，以便首先將更多相似的數(shù)據(jù)集合并在一起劳吠，但也可以定義此順序。

%%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

結(jié)果是另一個Seurat對象巩趁，但現(xiàn)在請注意痒玩，主動測定被稱為“集成”。這包含校正后的表達(dá)矩陣议慰，它是積分的最終輸出蠢古。

在這里，我們提取該矩陣并準(zhǔn)備將其傳輸回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])

輸出結(jié)果：

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

現(xiàn)在便瑟，我們將校正后的表達(dá)矩陣作為圖層存儲在AnnData對象中。我們還設(shè)置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>

現(xiàn)在我們有了積分的結(jié)果到涂，我們可以計算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)

正如我們之前所看到的践啄，批次是混合的，而標(biāo)簽是分開的沉御。人們很容易選擇基于UMAP的集成屿讽，但這并不能完全代表數(shù)據(jù)整合的質(zhì)量。在后面吠裆，我們將提出一些更嚴(yán)格評估集成整合方法的方法伐谈。

10.9 對集成整合進(jìn)行基準(zhǔn)測試

這里演示的方法是根據(jù)基準(zhǔn)測試實驗的結(jié)果選擇的，包括單細(xì)胞集成基準(zhǔn)測試項目试疙。該項目還制作了一個名為scib的軟件包诵棵，可用于運行一系列集成方法以及用于評估的指標(biāo)。在本節(jié)中祝旷，我們將展示如何使用此包來評估集成的質(zhì)量履澳。
scib指標(biāo)可以單獨運行嘶窄，但也有用于一次運行多個指標(biāo)的包裝器。在這里距贷，我們運行指標(biāo)的子集柄冲，使用metrics_fast()函數(shù)可以快速計算這些指標(biāo)。該函數(shù)需要幾個參數(shù)：原始未集成數(shù)據(jù)集忠蝗、集成數(shù)據(jù)集现横、批次鍵和標(biāo)簽鍵。根據(jù)集成方法的輸出阁最，我們可能還需要提供其他參數(shù)长赞，例如，這里我們使用embed參數(shù)指定用于scVI和scANVI的嵌入闽撤。您還可以使用其他參數(shù)控制某些指標(biāo)的運行方式。另請注意脯颜，您可能需要檢查對象的格式是否正確哟旗，以便scIB可以找到所需的信息。

讓我們?yōu)樯厦鎴?zhí)行的每個集成以及未集成的數(shù)據(jù)（在高度可變的基因選擇之后）運行指標(biāo)栋操。

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...

以下是單個集成整合的指標(biāo)結(jié)果之一的示例：

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

每行都是一個不同的指標(biāo)闸餐，值顯示該指標(biāo)的分?jǐn)?shù)。分?jǐn)?shù)介于0和1之間矾芙，其中1表示性能良好舍沙，0表示性能較差（如果需要，scib還可以返回某些指標(biāo)的未縮放分?jǐn)?shù)）剔宪。因為我們在這里只運行快速指標(biāo)拂铡，所以一些指標(biāo)具有NaN分?jǐn)?shù)。另請注意葱绒，某些指標(biāo)無法與某些輸出格式一起使用感帅，這也可能是返回NaN 值的原因。

為了比較這些方法地淀，將所有指標(biāo)結(jié)果放在一張表中非常有用失球。該代碼將它們組合起來并將它們整理成更方便的格式。

# 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

現(xiàn)在帮毁，我們將所有分?jǐn)?shù)放在一張表中实苞，其中指標(biāo)為列，方法為行烈疚。

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

對于某些指標(biāo)黔牵，分?jǐn)?shù)往往處于相對較小的范圍內(nèi)。為了強(qiáng)調(diào)方法之間的差異并將每個指標(biāo)放在相同的等級上胞得，我們對它們進(jìn)行了縮放荧止，以便表現(xiàn)最差的人得到0分屹电，表現(xiàn)最好的人得到1分，其他人的分?jǐn)?shù)介于兩者之間跃巡。

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

這些值現(xiàn)在可以更好地表示方法之間的差異危号。然而，值得注意的是素邪，縮放分?jǐn)?shù)只能用于比較這組特定集成的相對性能外莲。如果我們想添加另一種方法，我們需要再次執(zhí)行縮放兔朦。我們也不能說集成絕對是“好”偷线，只是說它比我們嘗試過的其他方法更好。這種縮放強(qiáng)調(diào)了方法之間的差異沽甥。例如声邦，如果我們的指標(biāo)分?jǐn)?shù)為0.92、0.94和0.96摆舟，則這些分?jǐn)?shù)將縮放為0亥曹、0.5和1.0。這使得第一種方法的得分看起來要差得多恨诱，盡管它只比其他兩種方法低一點點媳瞪，但仍然得到了很高的分?jǐn)?shù)。當(dāng)比較幾種方法并且當(dāng)它們獲得相似的原始分?jǐn)?shù)時照宝，這種效果會更大蛇受。您查看原始分?jǐn)?shù)還是縮放分?jǐn)?shù)取決于您是想關(guān)注絕對性能還是方法之間的性能差異。

評估指標(biāo)可以分為兩類厕鹃，一類衡量批次效應(yīng)的消除兢仰，另一類衡量生物變異的保存。我們可以通過取每組縮放值的平均值來計算每個類別的匯總分?jǐn)?shù)熊响。這種匯總分?jǐn)?shù)對于原始值沒有意義旨别，因為某些指標(biāo)始終比其他指標(biāo)產(chǎn)生更高的分?jǐn)?shù)（因此對平均值的影響更大）。

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

將兩個匯總分?jǐn)?shù)相互繪制可以表明每種方法的優(yōu)先級汗茄。有些會偏向批量校正秸弛，而另一些則傾向于保留生物變異。

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顯然是表現(xiàn)最差的递览，在批次去除和生物保護(hù)方面得分最低。其他三種方法具有相似的批次校正分?jǐn)?shù)瞳腌，其中scANVI在生物保護(hù)方面得分最高绞铃，其次是scVI和Seurat。

為了獲得每種方法的總體分?jǐn)?shù)嫂侍，我們可以將兩個匯總分?jǐn)?shù)結(jié)合起來儿捧。scIB論文建議權(quán)重為40%批量校正和 60% 生物保護(hù)荚坞，但您可能更愿意根據(jù)數(shù)據(jù)集的優(yōu)先級對事物進(jìn)行不同的權(quán)重。

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")

正如我們已經(jīng)看到的颓影，scVI和scANVI是表現(xiàn)最好的，其中scANVI得分略高懒鉴。值得注意的是诡挂，這只是如何針對該特定數(shù)據(jù)集運行這些指標(biāo)的示例，而不是對這些方法的正確評估临谱。為此璃俗，您應(yīng)該參考現(xiàn)有的基準(zhǔn)測試論文。特別是悉默，我們只運行了一小部分高性能方法和一部分指標(biāo)城豁。另外，請記住抄课，分?jǐn)?shù)與所使用的方法相關(guān)钮蛛，因此即使這些方法的性能幾乎相同，微小的差異也會被夸大剖膳。

現(xiàn)有的基準(zhǔn)測試建議的方法通常表現(xiàn)良好，但不同場景下的性能也可能存在很大差異岭辣。對于某些分析吱晒，您自己進(jìn)行集成評估可能是更好的選擇。scib包使這個過程變得更容易沦童，但它仍然是一項艱巨的任務(wù)仑濒，依賴于對基本事實的充分了解和對指標(biāo)的解釋。

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)

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者

人面猴
序言：七十年代末偷遗，一起剝皮案震驚了整個濱河市墩瞳，隨后出現(xiàn)的幾起案子，更是在濱河造成了極大的恐慌氏豌，老刑警劉巖喉酌，帶你破解...
沈念sama閱讀 206,013評論 6贊 481
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件，死亡現(xiàn)場離奇詭異泵喘，居然都是意外死亡泪电，警方通過查閱死者的電腦和手機(jī)，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 88,205評論 2贊 382
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進(jìn)店門纪铺，熙熙樓的掌柜王于貴愁眉苦臉地迎上來相速，“玉大人，你說我怎么就攤上這事鲜锚⊥晃埽” “怎么了苫拍？”我有些...
開封第一講書人閱讀 152,370評論 0贊 342
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵，是天一觀的道長旺隙。經(jīng)常有香客問我绒极，道長，這世上最難降的妖魔是什么催束？我笑而不...
開封第一講書人閱讀 55,168評論 1贊 278
?港島之戀（遺憾婚禮）
正文為了忘掉前任集峦，我火速辦了婚禮，結(jié)果婚禮上抠刺，老公的妹妹穿的比我還像新娘塔淤。我一直安慰自己，他們只是感情好速妖，可當(dāng)我...
茶點故事閱讀 64,153評論 5贊 371
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布高蜂。她就那樣靜靜地躺著，像睡著了一般罕容。火紅的嫁衣襯著肌膚如雪备恤。梳的紋絲不亂的頭發(fā)上，一...
開封第一講書人閱讀 48,954評論 1贊 283
城市分裂傳說
那天锦秒，我揣著相機(jī)與錄音露泊，去河邊找鬼。笑死旅择，一個胖子當(dāng)著我的面吹牛惭笑，可吹牛的內(nèi)容都是我干的。我是一名探鬼主播生真，決...
沈念sama閱讀 38,271評論 3贊 399
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼沉噩，長吁一口氣：“原來是場噩夢啊……” “哼！你這毒婦竟也來了柱蟀？” 一聲冷哼從身側(cè)響起川蒙，我...
開封第一講書人閱讀 36,916評論 0贊 259
萬榮殺人案實錄
序言：老撾萬榮一對情侶失蹤，失蹤者是張志新（化名）和其女友劉穎长已，沒想到半個月后畜眨，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體，經(jīng)...
沈念sama閱讀 43,382評論 1贊 300
?護(hù)林員之死
正文獨居荒郊野嶺守林人離奇死亡术瓮，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點故事閱讀 35,877評論 2贊 323
?白月光啟示錄
正文我和宋清朗相戀三年胶果，在試婚紗的時候發(fā)現(xiàn)自己被綠了。大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片斤斧。...
茶點故事閱讀 37,989評論 1贊 333
活死人
序言：一個原本活蹦亂跳的男人離奇死亡早抠，死狀恐怖，靈堂內(nèi)的尸體忽然破棺而出撬讽，到底是詐尸還是另有隱情蕊连，我是刑警寧澤悬垃，帶...
沈念sama閱讀 33,624評論 4贊 322
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布，位于F島的核電站甘苍，受9級特大地震影響尝蠕，放射性物質(zhì)發(fā)生泄漏。R本人自食惡果不足惜载庭，卻給世界環(huán)境...
茶點故事閱讀 39,209評論 3贊 307
男人毒藥：我在死后第九天來索命
文/蒙蒙一看彼、第九天我趴在偏房一處隱蔽的房頂上張望。院中可真熱鬧囚聚，春花似錦靖榕、人聲如沸。這莊子的主人今日做“春日...
開封第一講書人閱讀 30,199評論 0贊 19
一樁弒父案茁计，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽。三九已至谓松，卻和暖如春星压，著一層夾襖步出監(jiān)牢的瞬間，已是汗流浹背鬼譬。一陣腳步聲響...
開封第一講書人閱讀 31,418評論 1贊 260
情欲美人皮
我被黑心中介騙來泰國打工娜膘，沒想到剛下飛機(jī)就差點兒被人妖公主榨干…… 1. 我叫王不留，地道東北人优质。一個月前我還...
沈念sama閱讀 45,401評論 2贊 352
代替公主和親
正文我出身青樓劲绪，卻偏偏與公主長得像，于是被迫代替她去往敵國和親盆赤。傳聞我的和親對象是個殘疾皇子，可洞房花燭夜當(dāng)晚...
茶點故事閱讀 42,700評論 2贊 345