SingleR tutorial: easy cell type annotation

In this easy, step-by-step tutorial you will learn how to to do cell type annotation with the R package SingleR. SingleR is a popular reference-based automatic cell type annotation tool used to predict cell identities from gene expression profiles.

We will go step by step through the workflow, including preparing our input data, running SingleR, interpreting the results and some tips and tricks to get the most out of SingleR.

You will learn how to:

• get normalised counts from your Seurat object
• prepare a reference dataset for SingleR
• interpret SingleR’s cell type predictions (deltas and pruned labels)
• plot the cell types
• add your cell types to your original Seurat object

For this tutorial, I’ll be using RStudio, and you’ll need the packages SingleR, tidyverse. I’ll also use a reference dataset from the package celldex().

If you are more of a video-based learner, you can check out my Youtube video on how to perform cell type annotation with SingleR in R. Otherwise, just keep reading! And remember you can just copy the code by clicking on the top right corner of the code snippets in this post.

Let’s dive in!

You can also find the script for this tutorial in Biostatsquid’s GitHub page.

Before you start…

In this tutorial, we will learn how to do a cell type annotation of a single-cell dataset using R package SingleR. We will use sample data from the package celldex so if you don’t have a dataset to play around with, you can totally follow this tutorial!

If you are not that familiar with cell type annotation or automatic methods for cell type annnotation, you might want to check publications like this one, or this one. Or you can check my post on cell type annotation for scRNAseq by clicking here.

Let’s talk briefly about SingleR.

When we’re analysing a scRNAseq dataset, at some point, we need to figure out the cell identities of our cells. Are they macrophages? Tcells? Epithelial cells? As you know, a cell identity is based on gene expression: for example, T cells express CD3 (which is a gene that codes for the T cell receptor. Then, different subtypes of T cells express different genes, for example, T helper cells have a higher expression of CD4, whereas cytotoxic T cells express high levels of CD8. Obviously, it does not depend on one gene, but the main idea here is that a cell type is determined by gene expression profile, and the up and downregulation of many sets of genes define a cell type.

Since a scRNAseq dataset usually has thousands of cells, sometimes tens of thousands, it would be very tedious if we had to look into each of the cells separately, to see what genes they express. That’s where automatic cell type annotation methods come in.

One of the most popular strategies is reference-based annotation. In essence, the algorithm takes a scRNAseq reference dataset that’s already been annotated, and uses the gene expression values from the reference dataset to compare them to your dataset and annotate your own cells. In a way, it transfers the labels from one dataset to the other. This is what SingleR does. Given a reference dataset of samples (single-cell or bulk) with known labels, it assigns those labels to new cells from a test dataset based on similarities in their expression profiles. I won’t go into the method itself, but you can read more about it in the publication.

Great! Ready to start coding?

If you are a beginner in R, don’t be overwhelmed! This tutorial will go step-by-step and I will explain (almost!) every line of code so you know what is happening at each point of the workflow.

For this tutorial you will need R, or Rstudio, and you will need to install the following packages:

Set up your environment

It’s good practice to clear up your environment before you load in any data. This is the general structure I like to keep in all my scripts:

We will also need to load the relevant libraries:

Preparing inputs

We’ll start by loading the data. I’m going to use a dataset from the scRNAseq package of PBMCs, which are immune cells in the blood. We load it as a SingleCellExperiment (sce) object. SingleR can already take this as input. But since in single cell pipelines it’s more common to work with Seurat objects, we will convert this sce into a Seurat object and extract our counts from there.

1. Get normalised or raw counts

SingleR can use raw or normalised counts from our dataset. I’ll show you how to get both.

In the new Seurat version 5, we extract the raw counts using the function LayerData(), and the raw counts are stored in ‘counts’. Let’s check the raw counts for a few genes and the first 5 samples. To extract the normalised counts, we use the layer ‘data’.

In the code below, I assign raw counts to raw_counts, normalised counts to norm_counts, and then I check the counts for a few genes and the first 5 cells to compare raw and normalised counts.

> raw_counts <- LayerData(seu, assay = "RNA", layer = 'counts') #
> raw_counts[c('VIM', 'BCL2', 'TP53', 'CD4'),1:5]
4 x 5 sparse Matrix of class "dgCMatrix"
     AAACCTGAGAGCCCAA_H1B1ln1 AAACCTGAGGCGTACA_H1B1ln1 AAACCTGCAGGTGGAT_H1B1ln1 AAACCTGCAGTATCTG_H1B1ln1
VIM                         2                        .                        1                        2
BCL2                        .                        .                        2                        .
TP53                        .                        .                        .                        .
CD4                         .                        .                        1                        .
     AAACCTGCATCACAAC_H1B1ln1
VIM                         4
BCL2                        .
TP53                        .
CD4                         .
> norm_counts <- LayerData(seu, assay = "RNA", layer = 'data') #
> norm_counts[c('VIM', 'BCL2', 'TP53', 'CD4'),1:5]
4 x 5 sparse Matrix of class "dgCMatrix"
     AAACCTGAGAGCCCAA_H1B1ln1 AAACCTGAGGCGTACA_H1B1ln1 AAACCTGCAGGTGGAT_H1B1ln1 AAACCTGCAGTATCTG_H1B1ln1
VIM                  2.622001                        .                 1.614255                 2.492596
BCL2                 .                               .                 2.202576                 .       
TP53                 .                               .                 .                        .       
CD4                  .                               .                 1.614255                 .       
     AAACCTGCATCACAAC_H1B1ln1
VIM                  3.212371
BCL2                 .       
TP53                 .       
CD4                  .  

Nice! That’s the first thing SingleR will need. If we check the function with ?SingleR, you can see that the argument test, the dataset we want to annotate, is a numeric matrix were the rows are genes and the columns are cells. Or you can also use SummarizedExperiment object as input if you prefer.

2. Get the reference dataset

Next, let’s find a good reference dataset for our dataset. I’m going to use the HumanPrimaryCellAtlasData() from the R package celldex to annotate our cells. We’ll talk more about the choice of reference datasets at the end of the video in the tips & tricks section. Let’s load the data.

This celldex dataset has label.main, which has rougher or more general cell type labels, and label fine. We’ll go with label.main for this example.

Since we have a PBMC dataset, there’s a lot of cell types in the reference dataset, like astrocytes, osteoblasts, hepatocytes.. that are not present in our dataset. This is why I’m going to only keep some blood cell types. However, you need to be careful when removing cell types from your reference, since they may be present in your dataset.

Run SingleR

Now we’re ready to run SingleR. Excited?

SingleR can take many arguments, but the main ones are:

• test, our normalised/raw counts.
• ref, our reference counts, normalised (you can give it just the matrix of log transformed expression values, but since I have the summarizedexperiment format already, I’m just going to give it this).
• labels, the labels of each of the cells in our reference dataset.
• de.method, the differential expression method it uses to predict the cell types. For single-cell data, the authors recommend the Wilcoxon ranked test.

Let’s run SingleR! It will take more or less time depending on the size of your dataset and your reference.

Nice! Let’s check out the output. As you can see, the function returns a dataframe. We have our cell IDs, the predicted label, and then some statistics on that prediction.

> ct_ann %>% head()

DataFrame with 6 rows and 4 columns
                                                 scores           labels              delta.next    pruned.labels
                                               <matrix>      <character>               <numeric>      <character>
AAACCTGAGAGCCCAA_H1B1ln1 0.327776:0.321938:0.262826:...          NK_cell 0.000000000000000000000          NK_cell
AAACCTGAGGCGTACA_H1B1ln1 0.297983:0.273244:0.231393:... Pre-B_cell_CD34- 0.000000000000000000000 Pre-B_cell_CD34-
AAACCTGCAGGTGGAT_H1B1ln1 0.218938:0.185317:0.199509:...          T_cells 0.566511846755758274874          T_cells
AAACCTGCAGTATCTG_H1B1ln1 0.267314:0.226187:0.235404:... Pre-B_cell_CD34- 0.000000000000000000000 Pre-B_cell_CD34-
AAACCTGCATCACAAC_H1B1ln1 0.258897:0.350308:0.218881:...         Monocyte 0.000000000000000111022         Monocyte
AAACCTGCATGGTCAT_H1B1ln1 0.344614:0.419826:0.266799:...         Monocyte 0.000000000000000000000         Monocyte

What are ‘pruned.labels’ in SingleR output?

What about delta and pruned labels?

Well, by itself, the SingleR algorithm will always assign a label to every cell. But what if the cell’s true label isn’t in the reference dataset? It will still assign it a label. For this and other reasons, SingleR can assign incorrect labels.

The developers tried to mitigate this  by removing poor-quality predictions with “low” scores. They basically compute a “delta” value for each cell, which is the gap, or the difference between the score for the assigned label and the median score across all labels. If the delta is small, this indicates that the cell matches all labels with the same confidence, so the assigned label is not very meaningful.

This way, SingleR can discard cells with low delta values caused by (i) ambiguous assignments with closely related reference labels and (ii) incorrect assignments that match poorly to all reference labels – so in the pruned_labels column you will find ‘cleaner’ or more reliable labels.

We can check and compare the pruned labels and the initial labels. As we see the numbers change slightly. And if we check the number of cells without a pruned label, we see that 227 cells were cleaned up and assigned NA because their deltas were too small.

> unique(ct_ann$pruned.labels)

[1] "NK_cell"          "Pre-B_cell_CD34-" "T_cells"          "Monocyte"         "B_cell"          
 [6] "Pro-B_cell_CD34+" "Pro-Myelocyte"    "Platelets"        "Neutrophils"      "Myelocyte"       
[11] "DC"               "Erythroblast"     NA                
> table(ct_ann$pruned.labels)
B_cell               DC     Erythroblast         Monocyte        Myelocyte      Neutrophils 
            5580              714              145             8835             1049              736 
         NK_cell        Platelets Pre-B_cell_CD34- Pro-B_cell_CD34+    Pro-Myelocyte          T_cells 
8764              256            14474              720              313            16841 
> table(ct_ann$labels)
B_cell               DC     Erythroblast         Monocyte        Myelocyte      Neutrophils 
            5580              714              145             9056             1049              736 
         NK_cell        Platelets Pre-B_cell_CD34- Pro-B_cell_CD34+    Pro-Myelocyte          T_cells 
8765              256            14475              720              314            16844 
> summary(is.na(ct_ann$pruned.labels))
Mode   FALSE    TRUE 
logical   58427     227 

Inspect the quality of SingleR predictions

We can inspect the quality of the predictions using several functions.

The plotScoreHeatmap function plots a heatmap with the scores. In the y axis we have the different reference labels, which form these rows. Each column is a cell. So we expect cells with a high score for T cells and low for the other clusters to have the final label to be T cells. If the prediction for a certain cell is ambiguous, it will have high scores (yellow) for many of the cell identities.

We can also plot the delta distribution. It basically shows the distribution of the gap between the score for the assigned label and the score of the remaining labels) across cells assigned to each reference label. Again, remember that we want high deltas, as they give us more confidence in the assigned label.

Add SingleR predictions to Seurat object

Great! So now if you want to add these cell type annotations back to our original Seurat object, you can do it very easily using AddMetadata() from the Seurat object. We just need to make sure that the rownames are the cell IDs that the function will use to match the predictions.

Cell type annotation tips and tricks

• Use the same gene IDs in your reference and your query dataset. You need to make sure the gene annotation in your reference, is the same as the gene annotation in your dataset! Otherwise it won’t be able to use the genes that don’t match to compare gene expression profiles to predict the cell types.
• Choose a good reference dataset. The reference dataset you use is CRUCIAL. It will define which cell types can potentially be identified in your dataset. So choose your reference dataset carefully. Celldex is an R package that has some nice scRNAseq datasets, but in general it is better to look for publications with as similar methods and tissue of origin as you can. My recommendation is to run SingleR several times with different reference datasets and then find a consensus between all of them – this will give you more reliable cell type annotations.
• Use different annotation methods. Don’t trust a single method to do your cell type annotation! You can use other automatic annotation tools, but also check markers to see if things make sense. If you are interested in more cell type annotation videos, leave me a comment down below!

sessionInfo()

Check my sessionInfo() here in case you have trouble reproducting my steps:

And that is the end of this tutorial!

In this SingleR tutorial, I explained how to perform cell type annotation using the R package SingleR. Hope you found it useful and easy to follow!

Before you go, you might want to check:

Squids don't care much for coffee,

but Laura loves a hot cup in the morning!

If you like my content, you might consider buying me a coffee.

You can also leave a comment or a 'like' in my posts or Youtube channel, knowing that they're helpful really motivates me to keep going:)

Cheers and have a 'squidtastic' day!