Danger
Please access this document in its canonical location as the currently accessed page may not be rendered correctly: Integrating data using ingest and BBKNN
Integrating data using ingest and BBKNN#
The following tutorial describes a simple PCA-based method for integrating data we call ingest and compares it with BBKNN [Polanski19]. BBKNN integrates well with the Scanpy workflow and is accessible through the bbknn function.
The ingest function assumes an annotated reference dataset that captures the biological variability of interest. The rational is to fit a model on the reference data and use it to project new data. For the time being, this model is a PCA combined with a neighbor lookup search tree, for which we use UMAP’s implementation [McInnes18]. Similar PCA-based integrations have been used before, for instance, in [Weinreb18].
As ingest is simple and the procedure clear, the workflow is transparent and fast.
Like BBKNN, ingest leaves the data matrix itself invariant.
Unlike BBKNN, ingest solves the label mapping problem (like scmap) and maintains an embedding that might have desired properties like specific clusters or trajectories.
We refer to this asymmetric dataset integration as ingesting annotations from an annotated reference adata_ref
into an adata
that still lacks this annotation. It is different from learning a joint representation that integrates datasets in a symmetric way as BBKNN, Scanorma, Conos, CCA (e.g. in Seurat) or a conditional VAE (e.g. in scVI, trVAE) would do, but comparable to the initiall MNN implementation in scran. Take a look at tools in the external API or at the ecoystem page to get a start with other tools.
import scanpy as sc
import pandas as pd
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")
scanpy==1.5.0 anndata==0.7.1 umap==0.4.2 numpy==1.18.1 scipy==1.4.1 pandas==1.0.3 scikit-learn==0.22.1 statsmodels==0.11.0
PBMCs#
We consider an annotated reference dataset adata_ref
and a dataset for which you want to query labels and embeddings adata
.
# this is an earlier version of the dataset from the pbmc3k tutorial
adata_ref = sc.datasets.pbmc3k_processed()
adata = sc.datasets.pbmc68k_reduced()
To use sc.tl.ingest
, the datasets need to be defined on the same variables.
var_names = adata_ref.var_names.intersection(adata.var_names)
adata_ref = adata_ref[:, var_names]
adata = adata[:, var_names]
The model and graph (here PCA, neighbors, UMAP) trained on the reference data will explain the biological variation observed within it.
sc.pp.pca(adata_ref)
sc.pp.neighbors(adata_ref)
sc.tl.umap(adata_ref)
The manifold still looks essentially the same as in the clustering tutorial.
Mapping PBMCs using ingest#
Let’s map labels and embeddings from adata_ref
to adata
based on a chosen representation. Here, we use adata_ref.obsm['X_pca']
to map cluster labels and the UMAP coordinates.
sc.tl.ingest(adata, adata_ref, obs="louvain")
adata.uns["louvain_colors"] = adata_ref.uns["louvain_colors"] # fix colors
By comparing the ‘bulk_labels’ annotation with ‘louvain’, we see that the data has been reasonably mapped, only the annotation of dendritic cells seems ambiguous and might have been ambiiguous in adata
already.
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"]
While there seems to be some batch-effect in the monocytes and dendritic cell clusters, the new data is otherwise mapped relatively homogeneously.
The megakaryoctes are only present in adata_ref
and no cells from adata
map onto them. If interchanging reference data and query data, Megakaryocytes do not appear as a separate cluster anymore. This is an extreme case as the reference data is very small; but one should always question if the reference data contain enough biological variation to meaningfully accomodate query data.
Using BBKNN#
sc.tl.pca(adata_concat)
%%time
sc.external.pp.bbknn(adata_concat, batch_key="batch") # running bbknn 1.3.6
CPU times: user 1.67 s, sys: 749 ms, total: 2.42 s
Wall time: 324 ms
sc.tl.umap(adata_concat)
Also BBKNN doesn’t maintain the Megakaryocytes cluster. However, it seems to mix cells more homogeneously.
Pancreas#
The following data has been used in the scGen paper [Lotfollahi2019], has been used here, was curated here and can be downloaded from here (the BBKNN paper).
It contains data for human pancreas from 4 different studies [Segerstolpe2016,Baron2016,Wang2016,Muraro2016], which have been used in the seminal papers on single-cell dataset integration [Butler2018,Haghverdi2018] and many times ever since.
# 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)
Inspect the cell types observed in these studies.
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
To simplify visualization, let’s remove the 5 minority classes.
minority_classes = counts.index[-5:].tolist() # get the minority classes
adata_all = adata_all[~adata_all.obs.celltype.isin(minority_classes)] # actually subset
adata_all.obs.celltype.cat.reorder_categories( # reorder according to abundance
counts.index[:-5].tolist(), inplace=True
)
Seeing the batch effect#
sc.pp.pca(adata_all)
sc.pp.neighbors(adata_all)
sc.tl.umap(adata_all)
We observe a batch effect.
BBKNN#
It can be well-resolved using BBKNN [Polanski19].
%%time
sc.external.pp.bbknn(adata_all, batch_key="batch")
CPU times: user 1.89 s, sys: 810 µs, total: 1.9 s
Wall time: 1.89 s
sc.tl.umap(adata_all)
If one prefers to work more iteratively starting from one reference dataset, one can use ingest.
Mapping onto a reference batch using ingest#
Choose one reference batch for training the model and setting up the neighborhood graph (here, a PCA) and separate out all other batches.
As before, the model trained on the reference batch will explain the biological variation observed within it.
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.
Iteratively map labels (such as ‘celltype’) and embeddings (such as ‘X_pca’ and ‘X_umap’) from the reference data onto the query batches.
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")
... integrating batch 1
running ingest
finished (0:00:06)
... integrating batch 2
running ingest
finished (0:00:07)
... integrating batch 3
running ingest
finished (0:00:03)
Each of the query batches now carries annotation that has been contextualized with adata_ref
. By concatenating, we can view it together.
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"]
Compared to the BBKNN result, this is maintained clusters in a much more pronounced fashion. If one already observed a desired continuous structure (as in the hematopoietic datasets, for instance), ingest
allows to easily maintain this structure.
Evaluating consistency#
Let us subset the data to the query batches.
adata_query = adata_concat[adata_concat.obs.batch.isin(["1", "2", "3"])]
The following plot is a bit hard to read, hence, move on to confusion matrices below.
Cell types conserved across batches#
Let us first focus on cell types that are conserved with the reference, to simplify reading of the confusion matrix.
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 categories
obs_query_conserved.celltype.cat.remove_unused_categories(inplace=True)
# remove unused categories
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)
celltype_orig | alpha | beta | ductal | acinar | delta | gamma | endothelial | mast |
---|---|---|---|---|---|---|---|---|
celltype | ||||||||
alpha | 1819 | 3 | 6 | 0 | 1 | 25 | 0 | 6 |
beta | 49 | 804 | 4 | 1 | 10 | 21 | 0 | 0 |
ductal | 7 | 5 | 692 | 240 | 0 | 0 | 0 | 0 |
acinar | 2 | 3 | 3 | 168 | 0 | 3 | 0 | 0 |
delta | 5 | 4 | 0 | 0 | 305 | 73 | 0 | 0 |
gamma | 1 | 5 | 0 | 0 | 0 | 194 | 0 | 0 |
endothelial | 2 | 0 | 0 | 0 | 0 | 0 | 36 | 0 |
mast | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 1 |
Overall, the conserved cell types are also mapped as expected. The main exception are some acinar cells in the original annotation that appear as acinar cells. However, already the reference data is observed to feature a cluster of both acinar and ductal cells, which explains the discrepancy, and indicates a potential inconsistency in the initial annotation.
All cell types#
Let us now move on to look at all cell types.
pd.crosstab(adata_query.obs.celltype, adata_query.obs.celltype_orig)
celltype_orig | PSC | acinar | alpha | beta | co-expression | delta | dropped | ductal | endothelial | epsilon | gamma | mast | mesenchymal | mesenchyme | not applicable | unclassified endocrine |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
celltype | ||||||||||||||||
alpha | 0 | 0 | 1819 | 3 | 2 | 1 | 36 | 6 | 0 | 4 | 25 | 6 | 0 | 0 | 310 | 10 |
beta | 0 | 1 | 49 | 804 | 35 | 10 | 42 | 4 | 0 | 0 | 21 | 0 | 0 | 1 | 512 | 24 |
ductal | 0 | 240 | 7 | 5 | 0 | 0 | 38 | 692 | 0 | 0 | 0 | 0 | 1 | 0 | 102 | 1 |
acinar | 0 | 168 | 2 | 3 | 0 | 0 | 25 | 3 | 0 | 0 | 3 | 0 | 0 | 0 | 89 | 0 |
delta | 0 | 0 | 5 | 4 | 1 | 305 | 13 | 0 | 0 | 4 | 73 | 0 | 0 | 0 | 101 | 6 |
gamma | 0 | 0 | 1 | 5 | 0 | 0 | 1 | 0 | 0 | 2 | 194 | 0 | 0 | 0 | 14 | 0 |
endothelial | 1 | 0 | 2 | 0 | 1 | 0 | 7 | 0 | 36 | 0 | 0 | 0 | 0 | 6 | 7 | 0 |
activated_stellate | 49 | 1 | 1 | 3 | 0 | 0 | 11 | 8 | 0 | 0 | 0 | 0 | 79 | 20 | 17 | 0 |
quiescent_stellate | 4 | 0 | 1 | 1 | 0 | 0 | 5 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 |
macrophage | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 12 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 |
mast | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 |
We observe that PSC (pancreatic stellate cells) cells are in fact just inconsistently annotated and correctly mapped on ‘activated_stellate’ cells.
Also, it’s nice to see that ‘mesenchyme’ and ‘mesenchymal’ cells both map onto the same category. However, that category is again ‘activated_stellate’ and likely incorrect.
Visualizing distributions across batches#
Often, batches correspond to experiments that one wants to compare. Scanpy offers to convenient visualization possibilities for this.
a density plot
a partial visualization of a subset of categories/groups in an emnbedding
Density plot#
sc.tl.embedding_density(adata_concat, groupby="batch")
computing density on 'umap'