学习资料来源：
- scanpy主页：https://scanpy.readthedocs.io/en/stable/
- 官网：https://scanpy-tutorials.readthedocs.io/en/latest/integrating-data-using-ingest.html【注意教程有两个版本，这里是latest版本的学习笔记】
此教程描述了一种简单的基于pca的方法来整合数据，我们称之为ingest，并将其与BBKNN [Polanski19]进行了比较。BBKNN与Scanpy工作流整合得很好，可以通过BBKNN功能访问。
ingest函数假设有一个带注释的参考数据集，该数据集捕获了感兴趣的生物变异。合理的做法是在参考数据上拟合一个模型，并使用它来投射新的数据。
目前，这个模型是一个结合了邻居查找搜索树的PCA，我们将使用UMAP来实现这个功能。类似的基于pca的整合之前已经被使用过，例如，[Weinreb18]。ingest的优点：
- ingest简单，流程清晰且运行快
- 与BBKNN一样，ingest不改变数据matrix本身
- 不同于BBKNN的是，ingest解决了标签映射问题(如scmap)，并保持了可能具有特定cluster或轨迹等所需属性的嵌入
我们将这种非对称数据集整合称为将注释从带注释的参考adata_ref摄取到没有注释的数据中。它不同于学习以对称方式整合数据集的joint，如BBKNN, Scanorma, Conos, CCA(如Seurat)或有条件的VAE(如scVI, trVAE)，但可与scran中的初始MNN实现相比较。
01 PBMCs数据
用到两个数据集：一个带注释的参考数据集adata_ref和一个需要被注释和嵌入数据的数据集。
import scanpy as sc
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as pl
# verbosity: errors (0), warnings (1), info (2), hints (3)
sc.settings.verbosity = 1
sc.logging.print_versions()
sc.settings.set_figure_params(dpi=80, frameon=False, figsize=(3, 3), facecolor='white')
outdir = '/Pub/Users/zhangjuan/project/scanpy/Intergrating/'
# this is an earlier version of the dataset from the pbmc3k tutorial
adata_ref = sc.datasets.pbmc3k_processed()
adata = sc.datasets.pbmc68k_reduced()
环境中软件版本如下：
-----
anndata 0.8.0
scanpy 1.9.1
-----
PIL 9.2.0
beta_ufunc NA
binom_ufunc NA
colorama 0.4.5
cycler 0.10.0
cython_runtime NA
dateutil 2.8.2
h5py 3.7.0
hypergeom_ufunc NA
igraph 0.9.11
joblib 1.1.0
kiwisolver 1.4.4
leidenalg 0.8.10
llvmlite 0.38.1
louvain 0.7.1
matplotlib 3.4.3
mpl_toolkits NA
natsort 8.1.0
nbinom_ufunc NA
numba 0.55.2
numexpr 2.8.3
numpy 1.22.4
packaging 21.3
pandas 1.4.3
pkg_resources NA
pyparsing 3.0.9
pytz 2022.1
scipy 1.8.1
seaborn 0.11.2
session_info 1.0.0
setuptools 63.2.0
six 1.16.0
sklearn 1.1.1
statsmodels 0.13.2
texttable 1.6.4
threadpoolctl 3.1.0
-----
Python 3.9.13 | packaged by conda-forge | (main, May 27 2022, 16:58:50) [GCC 10.3.0]
Linux-5.10.0-9-amd64-x86_64-with-glibc2.31
-----
Session information updated at 2022-08-03 23:55
参考数据以及查询数据的情况：
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
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
... ... ... ... ...
TTTCGAACTCTCAT-1 1155 0.021104 3459.0 CD14+ Monocytes
TTTCTACTGAGGCA-1 1227 0.009294 3443.0 B cells
TTTCTACTTCCTCG-1 622 0.021971 1684.0 B cells
TTTGCATGAGAGGC-1 454 0.020548 1022.0 B cells
TTTGCATGCCTCAC-1 724 0.008065 1984.0 CD4 T cells
[2638 rows x 4 columns]
adata
AnnData object with n_obs × n_vars = 700 × 765
obs: 'bulk_labels', 'n_genes', 'percent_mito', 'n_counts', 'S_score', 'G2M_score', 'phase', 'louvain'
var: 'n_counts', 'means', 'dispersions', 'dispersions_norm', 'highly_variable'
uns: 'bulk_labels_colors', 'louvain', 'louvain_colors', 'neighbors', 'pca', 'rank_genes_groups'
obsm: 'X_pca', 'X_umap'
varm: 'PCs'
obsp: 'distances', 'connectivities'
使用sc.tl函数，需要在相同的变量（即需要提共有基因数据）上定义数据集
var_names = adata_ref.var_names.intersection(adata.var_names)
adata_ref = adata_ref[:, var_names]
adata = adata[:, var_names]
只有208个共同基因！
在参考数据上训练的模型和图形(这里是PCA，邻居，UMAP)将解释其中观察到的生物变异。
sc.pp.pca(adata_ref)
sc.pp.neighbors(adata_ref)
sc.tl.umap(adata_ref)
该流形看起来与聚类教程中的流形基本相同。
sc.pl.umap(adata_ref, color='louvain')
pl.savefig(outdir + "01-adata_ref-UMAP.png")
结果图：
使用ingest整合pbmc
基于选择的特征，我们使用adata_ref数据集合的labels and embeddings映射到adata。在这里，我们使用adata_ref.obsm['X_pca']来map 细胞亚群标签 and the UMAP 坐标。
sc.tl.ingest函数：
- obs：Labels’ keys in
adata_ref.obswhich need to be mapped toadata.obs(inferred for observation ofadata)
sc.tl.ingest(adata, adata_ref, obs='louvain')
# fix colors
adata.uns['louvain_colors'] = adata_ref.uns['louvain_colors']
sc.pl.umap(adata, color=['louvain', 'bulk_labels'], wspace=0.5)
## 修改图片边距，之前图片保存字体显示不全
pl.margins(0, 0)
pl.savefig(outdir + "02-adata-UMAP.png", bbox_inches='tight', dpi=300, pad_inches=0.1)
结果图：这个图片跟官方的结果是个镜像，180度反的，左边louvain注释是通过ingest整合的结果，右边bulk_labels是数据原有的注释标签：
通过比较louvain注释结果与bulk_labels注释结果，我们发现基本上都能对上，除了DC细胞。但是DC细胞可能在adata中就存在注释问题。
adata_concat = adata_ref.concatenate(adata, batch_categories=['ref', 'new'])
adata_concat.obs.louvain = adata_concat.obs.louvain.astype('category')
# fix category ordering
adata_concat.obs.louvain.cat.reorder_categories(adata_ref.obs.louvain.cat.categories, inplace=True)
# fix category colors
adata_concat.uns['louvain_colors'] = adata_ref.uns['louvain_colors']
sc.pl.umap(adata_concat, color=['batch', 'louvain'])
pl.margins(0, 0)
pl.savefig(outdir + "03-adata_ref-UMAP-new.png", bbox_inches='tight', dpi=300, pad_inches=0.1)
矫正之后的结果：虽然在单核细胞和树突状细胞簇中似乎存在一些批量效应，但是新的数据在其他方面被映射得相对均匀。
巨核细胞megakaryoctes只存在于 adata _ ref 中，adata中没有cell映射到它们
如果交换adata _ ref 和adata，巨核细胞megakaryoctes不再作为一个单独的cluster出现。这是一个极端的情况，因为参考数据非常小; 但是人们应该始终质疑参考数据是否包含足够的生物变异来有意义地容纳查询数据。
使用BBKNN整合数据
adata_concat是adata_ref（2638个细胞）与adata（700个细胞）两个矩阵合并后的数据
参数join默认采用inner方式即交集合并数据 : str (default: 'inner')，Use intersection ('inner') or union ('outer') of variables
合并后，adata_concat包含：3338个细胞和208个共同基因
adata_concat
AnnData object with n_obs × n_vars = 3338 × 208
obs: 'n_genes', 'percent_mito', 'n_counts', 'louvain', 'bulk_labels', 'S_score', 'G2M_score', 'phase', 'batch'
var: 'n_counts-new', 'means-new', 'dispersions-new', 'dispersions_norm-new', 'highly_variable-new', 'n_cells-ref'
uns: 'louvain_colors', 'batch_colors', 'pca'
obsm: 'X_pca', 'X_umap'
varm: 'PCs'
整合数据：
sc.tl.pca(adata_concat)
# running bbknn 1.3.6
sc.external.pp.bbknn(adata_concat, batch_key='batch')
sc.tl.umap(adata_concat)
sc.pl.umap(adata_concat, color=['batch', 'louvain'])
pl.savefig(outdir + "04-adata_concat-UMAP-BBKNN.png",bbox_inches='tight', dpi=300, pad_inches=0.1)
结果图：左图为adata_ref数据与adata数据整合在一起的情况，右边为注释后的结果
02 Pancreas
本数据集在the scGen paper [Lotfollahi19]，here中使用过，在here被核实过，可以在here (the BBKNN paper)这里进行下载。
它包含了来自4个不同研究（Segerstolpe16, Baron16, Wang16, Muraro16）的人类胰腺的数据，它已经在单细胞数据集整合的开创性论文中被使用(Butler18，Haghverdi18) ，并且从那时起被多次使用。
将数据下载下来并存放在data目录下：
# note that this collection of batches is already intersected on the genes
adata_all = sc.read('data/pancreas.h5ad', backup_url='https://www.dropbox.com/s/qj1jlm9w10wmt0u/pancreas.h5ad?dl=1')
adata_all.shape
(14693, 2448)
检查在这些研究中观察到的细胞类型：
counts = adata_all.obs.celltype.value_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
为了简化可视化，让我们删除5个细胞数量比较少的cluster：
# get the minority classes
minority_classes = counts.index[-5:].tolist()
# actually subset
adata_all = adata_all[~adata_all.obs.celltype.isin(minority_classes)]
# reorder according to abundance
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)
pl.savefig(outdir + "05-adata_all-UMAP.png",bbox_inches='tight', dpi=300, pad_inches=0.1)
可以看到很明显的批次效应：这里跟官网的图有点不一样，聚类形状差不多，但是就是有些cluster旋转了某个角度
使用BBKNN整合数据
上面数据中的批次可以很好地使用 BBKNN [Polanski19]来解决。
sc.external.pp.bbknn(adata_all, batch_key='batch')
sc.tl.umap(adata_all)
sc.pl.umap(adata_all, color=['batch', 'celltype'])
pl.savefig(outdir + "06-adata_all-UMAP_batchAdjust.png",bbox_inches='tight', dpi=300, pad_inches=0.1)
结果如下：
使用ingest整合数据
选择一个参考批处理来训练模型和建立邻域图(这里是PCA) ，并分离出所有其他批处理。
和以前一样，对参考批次进行训练的模型将解释在其中观察到的生物学变异。
这里数据有4个batch，选择batch 0 作为adata_ref，其余为adata
adata_ref = adata_all[adata_all.obs.batch == '0']
# Compute the PCA, neighbors and UMAP on the reference data
sc.pp.pca(adata_ref)
sc.pp.neighbors(adata_ref)
sc.tl.umap(adata_ref)
# The reference batch contains 12 of the 19 cell types across all batches.
sc.pl.umap(adata_ref, color='celltype')
pl.savefig(outdir + "07-adata_all-UMAP_ingest.png",bbox_inches='tight', dpi=300, pad_inches=0.1)
结果图如下：参考batch中包含所有19中细胞类型中的12种细胞
迭代地将label(例如‘ 细胞类型’)和embeddings(例如‘ X _ pca’和‘ X _ umap’)从参考数据映射到需要注释的数据中：
adatas = [adata_all[adata_all.obs.batch == i].copy() for i in ['1', '2', '3']]
# a bit more logging
sc.settings.verbosity = 2
for iadata, adata in enumerate(adatas):
print(f'... integrating batch {iadata+1}')
# save the original cell type
adata.obs['celltype_orig'] = adata.obs.celltype
sc.tl.ingest(adata, adata_ref, obs='celltype')
... integrating batch 1
running ingest
finished (0:00:15)
... integrating batch 2
running ingest
finished (0:00:06)
... integrating batch 3
running ingest
finished (0:00:02)
现在，每个查询的batch都带有已经与 adata _ ref 进行上下文关联的注释
使用concatenate，我们可以一起可视化看看结果：
adata_concat = adata_ref.concatenate(adatas)
adata_concat.obs.celltype = adata_concat.obs.celltype.astype('category')
# fix category ordering
adata_concat.obs.celltype.cat.reorder_categories(adata_ref.obs.celltype.cat.categories, inplace=True)
# fix category coloring
adata_concat.uns['celltype_colors'] = adata_ref.uns['celltype_colors']
sc.pl.umap(adata_concat, color=['batch', 'celltype'])
pl.savefig(outdir + "08-adata_concat-UMAP_ingest.png",bbox_inches='tight', dpi=300, pad_inches=0.1)
结果图：与BBKNN的结果相比，ingest 以一种更明显的方式维持cluster。如果已经观察到想要的连续结构(例如在造血数据集中)，ingest 可以轻松地维持这个结构。
评估一致性
提取数据中刚刚被当做查询的部分，即batch1,2,3，包含6113个细胞
## 评估一致性
adata_query = adata_concat[adata_concat.obs.batch.isin(['1', '2', '3'])]
adata_query
View of AnnData object with n_obs × n_vars = 6113 × 2448
# The following plot is a bit hard to read, hence, move on to confusion matrices below
sc.pl.umap(adata_query, color=['batch', 'celltype', 'celltype_orig'], wspace=0.4)
pl.savefig(outdir + "09-adata_query-consistency.png",bbox_inches='tight', dpi=300, pad_inches=0.1)
结果图：整合后数据的batch，整合后数据的注释，数据原有的注释label
不同batches中保守的细胞类型
我们首先来查看在reference batch与query batch中保守的细胞类型：
obs_query = adata_query.obs
# intersected categories
conserved_categories = obs_query.celltype.cat.categories.intersection(obs_query.celltype_orig.cat.categories)
# intersect categories
obs_query_conserved = obs_query.loc[obs_query.celltype.isin(conserved_categories) & obs_query.celltype_orig.isin(conserved_categories)]
# remove unused categoriyes
obs_query_conserved.celltype.cat.remove_unused_categories(inplace=True)
# remove unused categoriyes
obs_query_conserved.celltype_orig.cat.remove_unused_categories(inplace=True)
# fix category ordering
obs_query_conserved.celltype_orig.cat.reorder_categories(obs_query_conserved.celltype.cat.categories, inplace=True)
pd.crosstab(obs_query_conserved.celltype, obs_query_conserved.celltype_orig)
总的来说，保守的细胞类型按预期那样映射。意外的是原注释中一些腺泡细胞表现为腺泡细胞。然而，已经观察到的参考数据显示腺泡和导管细胞聚集，这解释了差异，并表明初始注释中可能存在不一致。
所有的细胞类型
现在，我们来看看所有的细胞类型：
pd.crosstab(adata_query.obs.celltype, adata_query.obs.celltype_orig)
结果：
我们观察到PSC(胰腺星状细胞)细胞实际上只是不一致地注释和正确地映射在“激活星状细胞”上；
此外，很高兴看到“间充质”和“间充质”细胞都属于同一类别。然而，这个类别仍然是'activated_stellate'，而且很可能是错误的
各批次的可视化分布
通常，批量与人们想要比较的实验相对应。Scanpy为此提供了方便的可视化可能性
Density plot
sc.tl.embedding_density(adata_concat, groupby='batch')
sc.pl.embedding_density(adata_concat, groupby='batch')
pl.savefig(outdir + "10-adata_concat-density_plot.png",bbox_inches='tight', dpi=300, pad_inches=0.1)
结果图：
嵌入群子集的部分可视化
for batch in ['1', '2', '3']:
sc.pl.umap(adata_concat, color='batch', groups=[batch])
out = outdir + "11-adata_concat-Partial_visualizaton"+batch+".png"
pl.savefig(out, bbox_inches='tight', dpi=300, pad_inches=0.1)
下回见~
