Annotation

Motivation¶

To understand your data better and make use of existing knowledge, it is important to figure out the “cellular identity” of each of the cells in your data. The process of labeling groups of cells in your data based on known (or sometimes unknown) cellular phenotypes is called “cell annotation”. Whereas there are many ways to annotate your cells (e.g. based on batch, disease, sex and more), in this notebook we will focus on the annotation of “cell types”.

A cell type is a cellular phenotype that is robust across datasets, identifiable by specific marker genes or proteins, and tied to a distinct biological function. A classic example is the plasma B cell — a white blood cell that secretes antibodies and can be identified by characteristic markers.

However, like with any categorization the size of categories and the borders drawn between them are partly subjective and can change over time, e.g. because new technologies allow for a higher resolution view of cells, or because specific “sub-phenotypes” that were not considered biologically meaningful are found to have important biological implications (see e.g. Kadur Lakshminarasimha Murthy et al., 2022). Cell types are therefore often further classified into “subtypes” or “cell states” (e.g. activated versus resting) and some researchers use the term “cell identity” to avoid this sometimes arbitrary distinction. For a more detailed discussion of this topic, we recommend the review by Wagner et al. Wagner et al., 2016 and the recently published review by Zeng Zeng, 2022.

Similarly, multiple cell types can be part of a single continuum, where one cell type might transition or differentiate into another. For example, in hematopoiesis cells differentiate from a stem cell into a specific immune cell type. Although hard borders between early and late stages of this differentiation are often drawn, the state of these cells can more accurately be described by the differentiation coordinate between the less and more differentiated cellular phenotypes.

There are multiple ways to annotate cells. We will give an overview of the most widely used approaches below. As we are working with transcriptomic data, each of these methods is ultimately based on the expression of specific genes or gene sets, or general transcriptomic similarity between cells.

Environment setup¶

We will continue working with the scRNA-seq dataset that we earlier preprocessed and will now annotate it.

Set figure parameters:

Load data¶

Let’s read in the toy dataset we will use for this tutorial. It includes a single sample (“site4-donor8”) of the data also used in other parts of the book. Moreover, cells that didn’t pass QC have already been removed.

Manual annotation¶

The classical or oldest way to perform cell type annotation is based on a single or small set of marker genes known to be associated with a particular cell type. This approach dates back to “pre-scRNA-seq times”, when single-cell data was low-dimensional (e.g., FACS data with gene panels consisting of no more than 30-40 genes). However, when no unique markers exist for a specific cell type, this approach can quickly become challenging and less objective, with combinations of markers or expression thresholds necessary for proper annotation. A robust set of marker genes and prior knowledge or annotation experience can help here, but the approach comes with the risk of unclear and subjective decision-making.

For manual annotation, the data is usually clustered before annotation, so that we can annotate groups of cells instead of making a per-cell call. This is not only less laborious but also more robust to noise: a single cell might not have a count for a specific marker even if it was expressed in that cell, simply due to the inherent sparsity of single-cell data. Clustering enables the detection of cells highly similar in overall gene expression and can therefore account for drop-outs at the single-cell level.

Finally, there are two angles from which to approach the marker-gene-based annotation. One option is to work from a table of marker genes for all the cell types you expect in your data and check in which clusters those genes are expressed. The other option is to check which genes are highly expressed in the clusters you defined and then check if they are associated with known cell types or states. If necessary, one can move back and forth between those approaches. We will show examples of both below.

From markers to cluster annotation¶

First, we list a set of markers for cell types in the bone marrow that is based on literature: previous papers that study specific cell types and subtypes and report marker genes for those cell types. Note that markers at the protein level (e.g., used for FACS) sometimes do not work as well in transcriptomic data; hence, using markers from RNA-based papers is often more likely to work. Moreover, sometimes markers in one dataset do not turn out to work as well in other datasets. Ideally, a marker set is therefore validated across multiple datasets. It is often useful to work together with experts: as a bioinformatician, try to team up with a biologist who has more extensive knowledge of the tissue, the biology, the expected cell types and markers, etc.

Subset to only the markers that were detected in our data. We will loop through all cell types and keep only the genes that we find in our adata object as markers for that cell type. This will prevent errors once we start plotting.

To see where these markers are expressed, we can work with a 2-dimensional visualization of the data, such as a UMAP. We calculate such an embedding here based on the scran-normalized count data, using only the highly deviant genes. Note that we first perform a PCA on the normalized counts to reduce the dimensionality of the data before we generate the UMAP.

To start we store our raw counts in .layers['counts'], so that we will still have access to them later if needed. We then set our adata.X to the scran-normalized, log-transformed counts.

We furthermore set our adata.var.highly_variable to the highly deviant genes. Scanpy uses this var column in downstream calculations, such as the PCA below.

Now perform PCA. We use the highly deviant genes (set as “highly variable” above) to reduce noise and strengthen signal in our data and set number of components to the default n=50. 50 is on the high side for data of a single sample, but it will ensure that we don’t ignore important variation in our data.

Calculate the neighbor graph based on the PCs:

And use that neighbor graph to calculate a 2-dimensional UMAP embedding of the data:

Now show expression of the markers using the calculated UMAP. We’ll limit ourselves to B/plasma cell subtypes for this example. Note from the marker dictionary above that there are three negative markers in our list: IGHD and IGHM for B1 B, and PAX5 for plasmablasts, meaning that this cell type is expected not to or to lowly express those markers.

Let’s list the B cell subtypes we want to show the markers for:

And now plot one UMAP per marker for each of the B cell subtypes. Note that we can only plot the markers that are present in our data.

Even markers for a single cell type are often expressed in different subsets of the data, i.e. individual markers are often not uniquely expressed in a single cell type. Rather, it is the intersection of those subsets that will tell you where your cell type of interest is.

Of note is that markers are often sparsely expressed, i.e., it is often only a subset of cells of a cell type in which a marker was detected. This is due to the nature of scRNA-seq data: we only sequence a small subset of the total amount of RNA molecules in the cell, and due to this subsampling, we will sometimes not sample transcripts from specific genes in a cell even if they were expressed in that cell. Therefore, we do not annotate single cells based on a minimum expression threshold of, e.g., a set of markers. Instead, we first subdivide the data into groups of similar cells (i.e., “partition” the data) by clustering, thereby accounting for “missing transcripts” of single genes and rather grouping based on overall transcriptomic similarity. We can then annotate those clusters based on their overall marker expression patterns.

Let us cluster our data now. We will use the Leiden algorithm Traag et al., 2019 as discussed in the Clustering chapter to define a grouping of our data into similar subsets of cells:

In case you want the partitioning to be finer, you can change the resolution to a higher level by changing the resolution parameter of the clustering:

You can also add the cluster numbers to the UMAP.

This clustering is much finer and, in some cases, will help you annotate the data with more detail. You can play around with the resolution parameter to find the setting that best captures the marker expression patterns you observe. In our case, we will stick to the original resolution:

Scrolling back up, you will see that cluster 3 consistently expresses Naive CD20+ B cell markers, while cluster 6 consistently expresses Transitional B cell markers.

We can also visualize this using a dotplot:

Using a combination of visual inspection of the UMAPs and the dotplot above we can now start annotating the clusters:

You might notice that the annotation of B1 B cells is difficult, with none of the clusters expressing all the B1 B markers and several clusters expressing some of the markers. We often see that markers that work for one dataset do not work as well for others. This can be due to differences in sequencing depth, but also due to other sources of variation between datasets or samples.

Let’s visualize our annotations so far:

From cluster differentially expressed genes to cluster annotation¶

Conversely, we can calculate marker genes per cluster and then look up whether we can link those marker genes to any known biology, such as cell types and/or states. For marker gene calculation of clusters, simple methods such as the Wilcoxon rank-sum test are thought to perform best Pullin & McCarthy, 2022. Importantly, as the definition of the clusters is based on the same data as used for these statistical tests, the p-values of these tests will be inflated as also described here Zhang et al., 2019.

Let’s calculate the differentially expressed genes for every cluster, compared to the rest of the cells in our adata:

We can visualize expression of the top differentially expressed genes per cluster with a standard scanpy dotplot:

As you can see above, a lot of the differentially expressed genes are highly expressed in multiple clusters. We can filter the differentially expressed genes to select for more cluster-specific differentially expressed genes:

Visualize the filtered genes :

Let’s take a look at cluster 9, which seems to have a set of relatively unique markers including CD247, MOM2, KLRD1, PRF1, and KLRF1. Some googling tells us that e.g. KLRD1 is a marker for NK cells Isola et al., 2021. In the UMAP we can see that these genes are expressed throughout cluster 9:

However, marker-based annotation can be sensitive to the cluster resolution you choose, the robustness and uniqueness of the marker sets you have, and your knowledge of the cell types to be expected in your data. Therefore, cluster annotations are not guaranteed to be definitive and should be interpreted with caution.

For this reason, the field is partly trying to move away from manual cluster annotation and rather moving towards automated annotation algorithms instead. The rest of this tutorial will focus on those options.

Before we move on, store the final bit of annotation information in our adata:

Automated annotation¶

General remarks¶

The remainder of the discussed methods will be methods for automated, rather than manual annotation of your data. Automated approaches are based on different principles, sometimes requiring pre-defined sets of markers, other times trained on pre-existing full scRNA-seq datasets. As discussed below, the resulting annotations can be of varying quality. It is therefore important to regard these methods as a starting point rather than an end-point of the annotation process. See also several reviews Pasquini et al., 2021, Abdelaal et al., 2019 for a more elaborate discussion of automated annotation methods.

The quality of automatically generated annotations can vary substantially. More specifically, the quality of the annotations depends on:

1. The type of classifier chosen: Previous benchmark studies have shown that different types of classifiers often perform comparably, with neural network-based methods generally not outperforming general-purpose models such as support vector machines or linear regression modelsAbdelaal et al., 2019, Pasquini et al., 2021, Huang & Zhang, 2021.
   

2. The quality of the data that the classifier was trained on. If the training data was not well annotated or annotated at low resolution, the classifier will do the same. Similarly, if the training data and/or its annotation was noisy, the classifier might not perform well.
   

3. The similarity of your own data to the data that the classifier was trained on. For example, if the classifier was trained on a drop-seq single cell dataset and your data is 10X single nucleus rather than single cell drop-seq, this might worsen the quality of the annotation. Classifiers trained on cross-dataset atlases including a diversity of datasets might give more robust and better quality annotations than classifiers trained on a single dataset. An example is the CellTypist (an automated annotation method that will be discussed more extensively below) classifier trained on the Human Lung Cell Atlas Sikkema et al., 2023 which includes 14 different lung datasets. This model is likely to perform better on new lung data than a model that was trained on a single lung dataset.

The aforementioned points highlight possible disadvantages of using classifiers, depending on the training data and model type. Nonetheless, there are several important advantages of using pre-trained classifiers to annotate your data. First, it is a fast and easy way to annotate your data. The annotation does not require the downloading nor preprocessing of the training data and sometimes merely involves the upload of your data to an online webpage. Second, these methods don’t rely on a partitioning of your data into clusters, as the manual annotation does. Third, pre-trained classifiers enable you to directly leverage the knowledge and information from previous studies, such as a high quality annotation. And finally, using such classifiers can help with harmonizing cell-type definitions across a field, thereby clearing the path towards a field-wide consensus on these definitions.

Finally, as these classifiers are often less transparent than e.g. manual marker-based annotation, a good uncertainty measure quantifying annotation uncertainty will improve the quality and usability of the method. We will discuss this more extensively further down.

Marker gene-based classifiers¶

One class of automated cell type annotation methods relies on a predefined set of marker genes. Cells are classified into cell types based on their expression levels of these marker genes. Examples of such methods are Garnett Pliner et al., 2019 and CellAssign Zhang et al., 2019. The more robust and generalizable the set of marker genes these models are based on, the better the model will perform. However, like with other models they are likely to be affected by batch effect-related differences between the data the model was trained on and the data that needs to be labeled. One of the advantages of these methods compared to models based on larger gene sets (see below) is that they are more transparent: we know on the basis of which genes the classification is done.
We will not show an example of marker-based classifiers in this notebook, but encourage you to explore these yourself if you are interested.

Classifiers based on a wider set of genes¶

It is worth noting that the methods discussed so far use only a small subset of the genes detected in the data: often a set of only 1 to ~10 marker genes per cell type is used. An alternative approach is to use a classifier that takes as input a larger set of genes (several thousands or more), thereby making more use of the breadth of scRNA-seq data. Such classifiers are trained on previously annotated datasets or atlases. Examples of these are CellTypist Conde et al., 2022 (see also https://www.celltypist.org, where data can be uploaded to a portal to get automated cell annotations) and Clustifyr Fu et al., 2020.

Let’s try out CellTypist on our data. Based on the CellTypist tutorial (https://www.celltypist.org/tutorials) we know we need to prepare our data so that counts are normalized to 10,000 counts per cell, then log1p-transformed:

We’ll now download the celltypist models for immune cells:

Let’s try out both the Immune_All_Low and Immune_All_High models (these annotate immune cell types finer annotation level (low) and coarser (high)):

For each of these, we can see which cell types it includes to see if bone marrow cell types are included:

Looks like the models include many different immune cell type progenitors!

Now let’s run the models. First the coarse one:

Transform the predictions to adata to get the full output...

...and copy the results to our original AnnData object:

Now the same for the finer annotations:

Now plot:

One way of getting a feeling for the quality of these annotations is by looking if the observed cell type similarities correspond to our expectations:

This dendrogram shows us that the cells are well clustered with each other (e.g., B cells largely clustering together). It is very important to check the automated annotation manually before simply taking them!

Annotation by mapping to a reference¶

A final way to annotate your data is based on mapping your data to an existing, annotated single-cell reference and then performing label transfer using the resulting joint embedding. This reference can for example be a single sample that you annotated manually before, after which you would like to transfer those annotations to the rest of your dataset. Alternatively, it can be a published and ideally well-curated existing reference. In this context we refer to the “new data”, i.e. the data to be mapped and annotated, as the “query”.

There are multiple existing methods that perform such “query-to-reference mapping”, including scArches Lotfollahi et al., 2022, Symphony Kang et al., 2021, and Azimuth (Seurat) Hao et al., 2021. All of these methods enable you to map a new dataset to an existing reference without needing to reintegrate the data from the reference and without needing access to the full reference data.

As query-to-reference mapping involves embedding new data into an existing low-dimensional representation of the reference data, the dimensions and axes of that low-dimensional representation are largely pre-defined before learning from the query. Therefore, learning and incorporating unseen variation that might be present in the query (both new biological variation, e.g. unseen cell types or states and new technical variation, i.e. unseen batch effects that need to be removed) can be a challenge for these models. As a result, integration of the query data with the reference might not always be optimal and batch effects might not be fully removed from the joint query-reference embedding. However, as cell type label transfer does not necessarily require perfect integration but merely close proximity of identical cell types in the embedding, even an imperfect mapping can still be extremely helpful in annotating your data.

We will be using cellmapper Lange, 2025, specifically with scArches as the kernel method for reference-mapping-based label transfer to the new data ()“query”). This takes as its basis an existing (variational autoencoder-based) model that embeds the reference data in a low-dimensional, batch-corrected space. It then slightly extends that model to enable the mapping of an unseen dataset into the same “latent space” (i.e. the low-dimensional embedding). This model extension also enables the learning and removal of batch effects present in the mapped dataset.

Let’s start by preparing our data for the mapping to a reference. scArches, the method that enables us to adapt an existing reference model to new data, requires raw, non-normalized counts. We will therefore keep our counts layer and remove all other layers from our adata to map. We will set our .X to those raw counts as well.

Moreover, it is important that we use the same input features (i.e. genes) as were used for training our reference model and that we put those features in the same order. The reference model’s feature information is stored together with the model. Let’s load the feature table.

The table has both gene names and gene IDs. As gene IDs are usually less subject to change over genome annotation versions, we will use those to subset our data. We will therefore set our row names for both our adata and the reference model features to gene_ids. Importantly, we have to make sure to also store the gene names for later use: these are much easier to understand than the gene IDs.

Now, let’s check if we have all the genes we need in our query data:

We are missing some genes. We will manually add those and set their counts to 0, as it seems like these genes were not detected in our data. Let’s create an AnnData object for those missing genes with only zero values (including our raw counts layer, which will be used for the mapping). We will concatenate that to our own AnnData objects afterwards.

Concatenate our original adata to the missing genes adata. To make sure we can do this concatenation without errors, we’ll remove the PCA matrix from varm.

Now subset to the genes used in the model and order correctly:

Check if our adata gene names correspond exactly to the required gene order:

We can now set the gene indices back to gene names for easy interpretation:

Finally, this reference model used adata.obs['batch'] as our batch variable. We will therefore check that we have this set to one value for our entire sample:

Now let’s talk about our reference model. The better our reference model, the better our label transfer will perform. Using well-annotated reference that integrates many different datasets and that matches your data well (same organ, same single-cell technology etc.) is ideal: such models are trained on a variety of datasets and batches and are therefore expected to be more robust to batch effects. However, such references do not exist yet for all tissues. For this tutorial we will use a reference model trained on the bone marrow samples that we have been using throughout the book, excluding the sample that we will be mapping. The reference model is an scvi-model (used for data integration) that generates a low-dimensional, integrated embedding of the input data, see also the scvi-publication Lopez et al., 2018. Note that this is a toy model generated for this tutorial and it should not be used in other contexts.

We will start by patching things up for compatibility, because this model was built on previous version of pandas.

Now, let’s load the model and pass it the adata which we want to map.

We will now update this reference model so that we can embed our own data (the “query”) in the same latent space as the reference. This requires training on our query data:

Now that we have updated the model, we can calculate the (ideally batch-corrected) latent representation of our query:

We can now use this newly calculated low-dimensional embedding as a basis for visualization and clustering. Let’s calculate the new UMAP using the scVI-based representation of the data.

To see if the mapping-based UMAP makes general sense, let’s look at a few markers and if their expression is localized to specific parts of the UMAP:

Now the essential step is that we can combine the inferred latent space embedding of our query data with the existing reference embedding. Using this joint embedding, we will not only be able to e.g., visualize and cluster the two together, but we can also do label transfer from the query to the reference.

Let’s load the reference embedding: this is often made publicly available with existing atlases.

We’ll store a variable specifying that these cells are from the reference.

Let’s see what’s in this reference object:

As you can see it has only 10 dimensions (in .X) which together represent the latent space embedding of the reference cells. Our query embedding that we calculated for our own data also has 10 dimensions. The 10 dimensions of the reference and query are the same and can be combined!

Moreover, it has cell type labels in .obs['cell_type']. We will use these labels to annotate our own data.

To perform the label transfer, we will first concatenate the reference and query data using the 10-dimensional embedding. To get there, we will create the same type of AnnData object from our query data as we have from the reference (with the embedding under .X) and concatenate the two. With that, we can jointly analyze reference and query including doing transfer from one to the other.

We will set the cell type labels in .obs["cell_type"] as None, so that we can see how we can annotate cell types based on the surroundings.

Let’s visualize the joint embedding with a UMAP.

We can visually get a first impression of whether the reference and query integrated well based on the UMAP:

The (partial) mixing of query and reference in this UMAP is a good sign! When mapping completely fails, you will often see a full separation of query and reference in the UMAP.

Now let’s look at the cell type annotations from the reference. All cells from the query are set to NA here as they don’t have annotations yet and shown in black.

We’ll make this figure a bit bigger so that we can read the legend well:

As you can already tell from the UMAP, we can guess the cell type of each of our own cells (in black) by looking at which cell types from the reference surrounding it. This is exactly what a nearest-neighbor-graph-based label transfer approach does: for each query cell it checks what is the most common cell type among its neighboring reference cells. The higher the fraction of reference cells coming from a single cell type, the more confident the label transfer is.

Let’s perform the KNN-based label transfer.

First we set up the label transfer model:

Now we perform the label transfer:

And store the results in our adata:

Let’s transfer the results to our query adata object which also has our UMAP and gene counts, so that we can visualize all of those together.

We can now visualize the transferred labels in our previously calculated UMAP of our own data:

Let’s set the figure size smaller again:

Based on the neighbors of each of our query cells we can not only guess the cell type these cells belong to, but also generate a measure for certainty of that label: if a cell has neighbors from several different cell types, our guess will be highly uncertain. This is relevant to assess to what extent we can “trust” the transferred labels! Let’s visualize the uncertainty scores:

Let’s check for each cell type label how high the label transfer uncertainty levels were. This gives us a first impression of which annotations are more contentious/need more manual checks.

You’ll notice that e.g. progenitor cells are often more difficult to distinguish than other cell types. Same for the rather unspecific category “Other T” cells in our annotations. We can see that pDC, a cell type that is known to be quite transcriptionally distinct and therefore easier to recognize and label, spreads from 0.0 to 0.5.

To incorporate this uncertainty information in our transferred labels, we can set cells with an uncertainty score above e.g. 0.2 to “unknown”:

Let’s see what our annotations look like after this filtering. Note the Unknown color in the legend and the UMAP.

To ease legibility, we can color only the “unknown” cells. This will make it easier for us to see how many of those there are. You can do the same with any of the other cell type labels.

There are quite many of them! These cells will need particularly careful manual reviewing. However, the low-uncertainty annotations surrounding the “unknown cells” will already give us a first idea of what cell type we can expect each cell to belong to.

Now let’s take a look at our more certain annotations. We will check for a few cell types (chosen at random here) to what extent the reference-transferred annotation matches our known marker genes from above. In reality, this should be done systematically for all annotations!

Conveniently, for each of these cell types we have markers in our dictionary. Let’s plot marker expression for all our newly annotated cell types. You will notice that marker expression generally corresponds to the automated annotations, a good sign!

As you can see, the marker groups are generally most highly expressed in the cells annotated with the matching label. This means these labels are likely (at least partially) correct!

Let’s go back one more time to our UMAP colored by uncertainty:

The uncertainty not only helps us identify regions where the algorithm is unsure which cell type a cell belongs to (e.g., because it falls between two annotated phenotypes), but can also highlight unseen cell types or new cell states. For example, your reference might consist of healthy cells while your query could be from a diseased sample. The uncertainty score can then highlight disease-specific cell states, as they might not have neighbors from the reference that consistently come from a single cell type. Especially when your reference is based on a large dataset, the uncertainty score is useful for flagging parts of the query data that could be worth investigating. Reference-based label transfer not only helps you annotate your data but can also speed up its exploration and interpretation. However, as with any metric, these uncertainty scores are often imperfect and, in some cases, fail to highlight new cell types or states. For a more extensive discussion of uncertainty metrics, see e.g. Engelmann et al., 2022.

Like with any of the methods discussed in this notebook, the quality of the transferred annotations depends on the quality of the “training data” (in this case the reference) and its annotations, the quality of the model, and the match of your own data with the training data!

The quality of the transferred annotations should therefore always be validated with manual inspection using marker gene expression and refinement of the initial annotations might be needed.

Contributors¶

We gratefully acknowledge the contributions of:

Authors¶

• Lisa Sikkema

• Maren Büttner

• Seo H. Kim

Reviewers¶

• Lukas Heumos