处理单细胞不可避免的一个问题就是样本整合问题。
那如何将不同器官，不同测序平台，不同物种之间的单细胞数据进行整合分析呢？
Scanpy使用python语言构建了一套完整的单细胞分析流程，其中就包括使用ingest和BBKNN整合方法。
单细胞转录数据分析之Scanpy：https://www.jianshu.com/p/e22a947e6c60
单细胞转录组之Scanpy - 轨迹推断/拟时序分析：https://www.jianshu.com/p/0b2ca0e0b544
单细胞转录组之Scanpy - 样本整合分析：https://www.jianshu.com/p/beef8a8be360
单细胞空间转录分析之Scanpy：https://www.jianshu.com/p/8dc231c06932
单细胞空间转录分析之Scanpy-多样本整合：https://www.jianshu.com/p/caa98aeac191
单细胞空间转录分析之Seurat：https://www.jianshu.com/p/c9a601ced91f
单细胞空间转录分析之Seurat-多样本整合（浅谈空间批次）：https://www.jianshu.com/p/609b04096b79
如果仅仅是通过相同基因把不同的矩阵拼接起来，我们通常得到的结果是: Fig1 Raw。
细胞通常是以不同样本，不同平台，或者不同时间点等分开，也就是说两组数据中相同的细胞类别是并没有能合在一起。
当然，如果按基因直接整合后，细胞能均匀混合，我们是可以直接进行矩阵合并后分析的，但是当整合后混合结果不理想时，我们可能需要在整合时去除批次，主要是样本批次，平台批次，不同时间批次等，保留生物学意义。
因此，一系列单细胞整合分析工具（去批次）应运而生，Tran等人在Genome Biology上发表了一篇比较14种单细胞批次矫正的方法的文章：A benchmark of batch-effect correction methods for single-cell RNA sequencing data。 https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1850-9
这儿我们不讲这些整合方法的好坏，主要还是学习Scanpy提供的单细胞整合方法：包括ingest和BBKNN（Batch balanced kNN）。其中ingest主要用于含有reference的数据，而BBKNN主要用于整合（去批次），当然整合的方法也可以用来做reference。
ingest注释：非常简单，过程清晰，工作流程是透明且快速的。主要是在参考数据上拟合模型并使用它来投影新数据。该函数使用knn分类器来映射标签，使用UMAP来映射嵌入。
导入相关包
import scanpy as sc
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as pt
import bbknn
sc.settings.verbosity = 1 # verbosity: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_versions()
sc.settings.set_figure_params(dpi=80, frameon=False, figsize=(3, 3), facecolor='white')
import os
os.getcwd() ##查看当前路径
os.chdir('./Integrate/ingest') ##修改路径
os.getcwd()
results_file = 'ingest.h5ad'
读取数据 （Scanpy自带的两个数据集，一个是pbmc3k的，另一个是pbmc68k的部分细胞，都已经将细胞类别注释好了）
adata_ref = sc.datasets.pbmc3k_processed() # this is an earlier version of the dataset from the pbmc3k tutorial
adata = sc.datasets.pbmc68k_reduced()
我们可以查看一下数据集的内容：
adata_ref
AnnData object with n_obs × n_vars = 2638 × 1838
obs: 'n_genes', 'percent_mito', 'n_counts', 'louvain'
var: 'n_cells'
uns: 'draw_graph', 'louvain', 'louvain_colors', 'neighbors', 'pca', 'rank_genes_groups'
obsm: 'X_pca', 'X_tsne', 'X_umap', 'X_draw_graph_fr'
varm: 'PCs'
obsp: 'distances', 'connectivities'
adata_ref.obs.head()
n_genes percent_mito n_counts louvain
index
AAACATACAACCAC-1 781 0.030178 2419.0 CD4 T cells
AAACATTGAGCTAC-1 1352 0.037936 4903.0 B cells
AAACATTGATCAGC-1 1131 0.008897 3147.0 CD4 T cells
AAACCGTGCTTCCG-1 960 0.017431 2639.0 CD14+ Monocytes
AAACCGTGTATGCG-1 522 0.012245 980.0 NK cells
adata_ref.obs.louvain.value_counts()
CD4 T cells 1144
CD14+ Monocytes 480
B cells 342
CD8 T cells 316
NK cells 154
FCGR3A+ Monocytes 150
Dendritic cells 37
Megakaryocytes 15
Name: louvain, dtype: int64
data.obs.head()
bulk_labels n_genes percent_mito n_counts S_score G2M_score phase louvain
index
AAAGCCTGGCTAAC-1 CD14+ Monocyte 1003 0.023856 2557.0 -0.119160 -0.816889 G1 1
AAATTCGATGCACA-1 Dendritic 1080 0.027458 2695.0 0.067026 -0.889498 S 1
AACACGTGGTCTTT-1 CD56+ NK 1228 0.016819 3389.0 -0.147977 -0.941749 G1 3
AAGTGCACGTGCTA-1 CD4+/CD25 T Reg 1007 0.011797 2204.0 0.065216 1.469291 G2M 9
ACACGAACGGAGTG-1 Dendritic 1178 0.017277 3878.0 -0.122974 -0.868185 G1 2
可以看到参考数据集adata_ref已经注释好了细胞类别信息，这儿作为reference数据集。
可以看到refence细胞类别在umap上的分布。
var_names = adata_ref.var_names.intersection(adata.var_names) #提取共有基因
adata_ref = adata_ref[:, var_names]
adata = adata[:, var_names]
sc.pp.pca(adata_ref) ##对参考数据集进行降维
sc.pp.neighbors(adata_ref)
sc.tl.umap(adata_ref)
sc.pl.umap(adata_ref, color='louvain',legend_loc='on data',legend_fontsize='xx-small')
pt.savefig('pbmc3k_ref_umap.pdf')
使用ingest预测pbmc68k数据集（这儿pbmc68k假设是不知道细胞类别的）
sc.tl.ingest(adata, adata_ref, obs='louvain')
adata.uns['louvain_colors'] = adata_ref.uns['louvain_colors'] # fix colors
sc.pl.umap(adata, color=['louvain', 'bulk_labels'], wspace=0.5,legend_fontsize='xx-small')
pt.savefig('pbmc_data_umap.pdf')
左边显示了预测数据pbmc68k预测到的细胞类别，右边是pbmc68k原来的细胞标签。通过将“ bulk_labels”注释与“ louvain”进行比较，我们看到数据已被合理地映射，只有树突状细胞的注释似乎是模棱两可的。
对数据集合并，查看是否存在批次效应。
adata_concat = adata_ref.concatenate(adata, batch_categories=['ref', 'new']) #合并数据集
adata_concat.obs.louvain = adata_concat.obs.louvain.astype('category')
adata_concat.obs.louvain.cat.reorder_categories(adata_ref.obs.louvain.cat.categories, inplace=True) # fix category ordering
adata_concat.uns['louvain_colors'] = adata_ref.uns['louvain_colors'] # fix category colors
sc.pl.umap(adata_concat, color=['batch', 'louvain'],legend_fontsize='xx-small')
pt.savefig('pbmc_batch_umap.pdf')
尽管在单核细胞和树突状细胞簇中似乎有一些批处理效应，但新数据在其他方面相对均匀地映射。
使用BBKNN整合
import bbknn
sc.tl.pca(adata_concat)
sc.external.pp.bbknn(adata_concat, batch_key='batch')
sc.pl.umap(adata_concat, color=['batch', 'louvain'],legend_fontsize='xx-small')
pt.savefig('pbmc_batch_umap_bbknn.pdf')
我们可以看到使用bbknn整合后，也不能维护巨核细胞簇。但是，似乎可以更均匀地混合细胞。
使用BBKNN整合，Pancreas胰腺数据
BBKNN，一个快速和直观的批处理效果去除工具，可以直接在scanpy工作流中使用。
它包含来自4个不同研究（Segerstolpe，Baron，Wang，Muraro）的人类胰腺数据，这些数据已在有关单细胞数据集集成的开创性论文。数据可以直接在这网址下载：https://www.dropbox.com/s/qj1jlm9w10wmt0u/pancreas.h5ad?dl=1
导入包，设置输出路径
import scanpy as sc
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as pt
import bbknn
sc.settings.set_figure_params(dpi=80, frameon=False, figsize=(3, 3), facecolor='white')
import os
os.getcwd() ##查看当前路径
os.chdir(./Integrate/BBKNN') ##修改路径
os.getcwd()
results_file = 'bbknn.h5ad'
导入数据
adata_all = sc.read("./Integrate/BBKNN/pancreas.h5ad")
adata_all.shape
counts = adata_all.obs.celltype.value_counts()
counts
counts
alpha 4214
beta 3354
ductal 1804
acinar 1368
not applicable 1154
delta 917
gamma 571
endothelial 289
activated_stellate 284
dropped 178
quiescent_stellate 173
mesenchymal 80
macrophage 55
PSC 54
unclassified endocrine 41
co-expression 39
mast 32
epsilon 28
mesenchyme 27
schwann 13
t_cell 7
MHC class II 5
unclear 4
unclassified 2
Name: celltype, dtype: int64
minority_classes = counts.index[-5:].tolist() # 获得细胞数目少的类别
adata_all = adata_all[
~adata_all.obs.celltype.isin(minority_classes)]
adata_all.obs.celltype.cat.reorder_categories(counts.index[:-5].tolist(), inplace=True) # 去除细胞少的类别
降维
sc.pp.pca(adata_all)
sc.pp.neighbors(adata_all)
sc.tl.umap(adata_all)
sc.pl.umap(adata_all, color=['batch', 'celltype'], palette=sc.pl.palettes.vega_20_scanpy,legend_fontsize='xx-small')
pt.savefig('pancreas_batch_umap.pdf')
这儿我们可以看到明显的批次效应存在，样本之间的批次很严重。
使用BBKNN来整合数据
sc.external.pp.bbknn(adata_all, batch_key='batch')
sc.tl.umap(adata_all)
sc.pl.umap(adata_all, color=['batch', 'celltype'],legend_fontsize='xx-small')
pt.savefig('pancreas_batch_umap_bbknn.pdf')
我们通过BBKNN整合之后，可以明显发现批次效应去除非常明显，右边相同细胞类别能较好的合并在一起。
如果知道其中一套数据集中的细胞类别，这个数据集可以作为reference，就可以用上面提到的ingest。
我们这儿试着把batch 0作为reference，给其他数据集进行类别注释，并评估注释类别的准确性。
adata_ref = adata_all[adata_all.obs.batch == '0']
sc.pp.pca(adata_ref)
sc.pp.neighbors(adata_ref)
sc.tl.umap(adata_ref)
sc.pl.umap(adata_ref, color='celltype',legend_loc='on data',legend_fontsize='xx-small')
pt.savefig('pancreas_batch0_umap.pdf')
可以看到单一reference数据集细胞类别清晰。
使用ingest预测其它3个数据集的细胞类别
adatas = [adata_all[adata_all.obs.batch == i].copy() for i in ['1', '2', '3']]
sc.settings.verbosity = 2 # a bit more logging
for iadata, adata in enumerate(adatas):
print(f'... integrating batch {iadata+1}')
adata.obs['celltype_orig'] = adata.obs.celltype # save the original cell type
sc.tl.ingest(adata, adata_ref, obs='celltype')
adata_concat = adata_ref.concatenate(adatas)
adata_concat.obs.celltype = adata_concat.obs.celltype.astype('category')
adata_concat.obs.celltype.cat.reorder_categories(adata_ref.obs.celltype.cat.categories, inplace=True) # fix category ordering
adata_concat.uns['celltype_colors'] = adata_ref.uns['celltype_colors'] # fix category coloring
sc.pl.umap(adata_concat, color=['batch', 'celltype','celltype_orig'],legend_fontsize='xx-small')
pt.savefig('pancreas_batch_umap_ingest.pdf')
与BBKNN的结果相比，这是以一种更加明显的方式保持分群。
预测数据细胞类别一致性评估
adata_query = adata_concat[adata_concat.obs.batch.isin(['1', '2', '3'])] #只提取batch 1，2，3
sc.pl.umap(adata_query, color=['batch', 'celltype', 'celltype_orig'], wspace=0.4,legend_fontsize='xx-small')
pt.savefig('pancreas_batch_umap_ingest_query.pdf')
这个结果依然不能很好的反映一致性，让我们首先关注与参考保守的细胞类型，以简化混淆矩阵的reads。
obs_query = adata_query.obs
conserved_categories = obs_query.celltype.cat.categories.intersection(obs_query.celltype_orig.cat.categories) # intersected categories
obs_query_conserved = obs_query.loc[obs_query.celltype.isin(conserved_categories) & obs_query.celltype_orig.isin(conserved_categories)] # intersect categories
obs_query_conserved.celltype.cat.remove_unused_categories(inplace=True) # remove unused categoriyes
obs_query_conserved.celltype_orig.cat.remove_unused_categories(inplace=True) # remove unused categoriyes
obs_query_conserved.celltype_orig.cat.reorder_categories(obs_query_conserved.celltype.cat.categories, inplace=True) # fix category ordering
pd.crosstab(obs_query_conserved.celltype, obs_query_conserved.celltype_orig)
pd.crosstab(adata_query.obs.celltype, adata_query.obs.celltype_orig)
pd.crosstab(adata_query.obs.celltype, adata_query.obs.celltype_orig)
celltype_orig PSC acinar alpha beta ... mesenchymal mesenchyme not applicable unclassified endocrine
celltype ...
alpha 0 0 1819 3 ... 0 0 310 10
beta 0 1 49 804 ... 0 1 511 24
ductal 0 240 7 5 ... 1 0 102 1
acinar 0 168 2 3 ... 0 0 89 0
delta 0 0 5 4 ... 0 0 102 6
gamma 0 0 1 5 ... 0 0 14 0
endothelial 1 0 2 0 ... 0 6 7 0
activated_stellate 49 1 1 3 ... 79 20 17 0
quiescent_stellate 4 0 1 1 ... 0 0 1 0
macrophage 0 0 1 1 ... 0 0 1 0
mast 0 0 0 0 ... 0 0 0 0
[11 rows x 16 columns]
sc.tl.embedding_density(adata_concat, groupby='batch')
sc.pl.embedding_density(adata_concat, groupby='batch')
pt.savefig('pancreas_batch_umap_density.pdf')
可视化分布的批次
总体而言，保守细胞类型也按预期进行映射。主要的例外是原始注释中的某些腺泡细胞显示为腺泡细胞。然而，已经观察到参考数据具有腺泡和导管细胞的簇，这解释了差异，并指示了初始注释中的潜在不一致。
