Intro to Single Cell RNA Sequencing

2026 Jul 23 See all posts
Intro to Single Cell RNA Sequencing

Single cell RNA sequencing, known as scRNA-seq, is the process of determining the level of gene expression across individual cells. Before it, we could only measure average gene expression across many cell types in a sample, hiding cell-type-specific responses. In the mid-2010s, special "droplet" machines made scRNA-seq practical by isolating individual cells from a sample into droplets that contain chemical machinery to tag both the cell and individual RNA molecules in the cell. That is a big upgrade. For example, bulk sequencing had made clear in the 2000s that autoimmune diseases like lupus create chronic "alarm" markers in their white blood cells. scRNA-seq showed us exactly which type of white blood cells were responsible for the expression, giving us a much deeper understanding of the pathology.

Although these "droplet" machines are advanced, the raw data they output still requires a significant amount of statistical analysis to classify cells into cell types and build a picture of gene expression per cell type. That downstream analysis is mostly packaged in libraries like Seurat or Scanpy, but they contain many knobs that are hard to know how to use without a deeper understanding of what the analysis is doing or even the broader sample collection process. The "hello world" of scRNA-seq datasets is a sample of ~3k mononuclear white blood cells (Peripheral Blood Mononuclear Cells, PBMCs) from a healthy donor called pbmc3k. This post supplies the missing conceptual background to really understand the Seurat tutorial on that dataset.

If you are coming from a general computing background (like me) one relevant bit of biology that you need to know is the relationship between all the "NAs" (nucleic acids) like DNA, mRNA, cDNA etc. Roughly speaking, the nucleus contains DNA which is a full blueprint of every gene the cell might use. The epigenome, also inside the nucleus, is an on/off filter which can block access to parts of this blueprint. Parts of the blueprint are read (transcription) by nucleus internal machinery and exported as instructions via messenger RNA (mRNA). Builders outside the nucleus (ribosomes) use those instructions to build proteins and other things (translation). It's generally possible to deterministically reverse mRNA (reverse transcription) back into a mirrored copy of DNA called complementary DNA (cDNA). It's mirrored in the sense that DNA itself is constructed of pairs of base molecules (ACTG) which only pair in unique, deterministic ways. So knowing one element of the pair allows you to infer the other one. With that biology in mind, we can start with how the raw dataset is obtained in a wet lab.

Wet Lab

The wet lab steps are particularly important to understand where errors can occur and why we have to take various steps for quality control. The process is roughly as follows.

  1. Extract PBMCs. Blood is drawn from the donor, diluted and layered on top of a special solution called a Ficoll-Paque. It's then spun in a centrifuge where substances denser than ficoll (like red blood cells and non-PBMC white blood cells) slip below it and substances less dense stay above (like PBMCs). PBMCs can then be pipetted out visually from above the ficoll separator.
  2. Cryopreserve. Usually blood sampling and the droplet machine are done by different people/in different places so they have to be decoupled. So the PBMCs are frozen before transport, notably using dimethyl sulfoxide to stop ice crystals from forming and breaking cell membranes.
  3. Prepare for droplet machine. Samples are thawed quickly around 37C. The cells can stay alive for a bit if kept in the right culture (nutrients/temperature), allowing you to for example apply some treatment to them. You wash the cells with a dye that dead/leaky cells will grab on to (trypan blue), so you can identify and remove them.
  4. Run through droplet machine. The sample is then fed through a machine like 10x Genomics' Chromium. The machine will force the sample through one extremely small tunnel and force gel beads containing reactants and chemical barcode machinery through another tunnel. The tunnels meet and create a single droplet with a cell, reactants and barcode information to identify the cell and its molecules. This process is not perfect, the machine can be tuned to either have more doublets (2 cells in a droplet) or more empty droplets. Doublets are basically data loss, so usually it's tuned for a very low doublet rate in exchange for a fair amount of cellless droplets. The output is millions of tiny droplets suspended in solution.
  5. Run the droplet reactions. Heat the solution (using a "thermal cycler" machine) on a schedule to trigger the reactions in each droplet. The reaction that happens is reverse transcription (RT) which converts some of the expressed mRNA floating around the cell back into complementary DNA (cDNA) and simultaneously tags it with both a cellID and specific molecule ID (unique molecular identifier or UMI). Now each piece of cDNA that was captured inside the droplet can be identified globally. Importantly, some reactions might just execute better than others, resulting in a "deeper" capture where the UMI counts are higher across all genes in that cell.
  6. Sequence everything. Since each piece of cDNA has already been tagged with a cell and molecule ID, we can break all the droplets into a single large solution. Use a special enzyme (polymerase) to multiply the copies of the cDNA so it's easier to detect in a sequencer. Notably, not all cDNA gets scaled by the same factor, but since we have the UMIs on each molecule we can deduplicate. We then take the amplified cDNA (different number of duplicates per UMI) and sequence it in a machine like Illumina. There is further sampling error with the sequencing machine: you get a total number of reads which are used on unevenly amplified, cell-tagged cDNA. Illumina outputs a raw sequence file containing the cDNA reads and the quality/confidence of those reads.
  7. Convert sequence file to cell count matrix. Software like Cell Ranger understands the chemistry of a 10x Genomics droplet machine as well as sequence files. It takes a sequence file, peels off the cellID + UMI, matches the remaining cDNA to a known gene, deduplicates UMIs then writes a gene x cell count matrix into a .mtx file (a standardized Matrix Market file).

The "Gene / cell matrix filtered" file from the pbmc3k dataset contains this final .mtx file. Filtered just means it has removed empty droplets that may have arisen in step 4. There still may be doublets present, though usually not a significant quantity of them to mess with the end result. It looks like this:

%%MatrixMarket matrix coordinate real general
%
32738 2700 2286884
32709 1 4
32707 1 1
32706 1 10
32704 1 1
...

The header tells you that the matrix is "coordinate" oriented, meaning the non-zeros are specified instead of a full gene x cell matrix. It's more efficient that way since the matrix is sparse, many genes are only expressed in a few cells. The first row gives you the dimensions: 32,738 genes, 2700 cells and 2,286,884 count records or non-zero matrix entries. The schema for each subsequent row is <gene index in genes.tsv> <cell index in barcodes.tsv> <UMI count>, where the indices are 1-based, not 0-based. You can think of it like a list of (gene, cell, count) records, which we'll represent as a genes x cells matrix which is mostly zeroes.

Given this sparse matrix of gene x cell counts, the statistical task is to first group the cells based on their gene expression and then label those groups with a specific cell type based on a biological understanding of the relationship between gene expression and cell type.

Data Loading

The first step is to load the matrix into an appropriate python library. While it would be possible to use just the standard library (dicts/lists), some of the statistical algorithms would be prohibitively slow on a matrix with millions of non-zero entries. So we use SciPy, which is optimized for this kind of work. scipy.io.mmread will read a matrix market file into a coordinate format (COO) sparse matrix by default (notably changing to a different default soon). For some reason, this COO matrix format doesn't support all the operations we need and the "compressed sparse row" (CSR) format is just superior in all ways so we convert to that with tocsr(). Cell Ranger stores genes × cells; for cell typing we often transpose with .T so each row is one cell, but the counts themselves are easier to read as genes × cells. Replacing the gene indexes with gene names is helpful as well for later analysis. The loaded matrix conceptually looks like:

\[ \begin{array}{r|ccc} & \text{cell}_1 & \text{cell}_2 & \text{cell}_3 & \cdots \\ \hline \text{CD3D} & 1 & 0 & 0 \\ \text{MS4A1} & 5 & 2 & 192 \\ \text{LYZ} & 0 & 0 & 97 \\ \vdots & \vdots & \vdots & \vdots \end{array} \]

Quality Control (QC)

The gene counts we get may still have several sources of experimental error:

To filter out these spurious readings, we do the following (thresholds just experimentally determined / recommended by the Seurat folks):

Normalize

After QC, we're left with a genes x cells matrix containing (we hope) only valid cells. The raw counts of expressed genes per cell in the matrix depend on how well the droplet reaction executes (Step 5), sequencing unevenly amplified cDNA (Step 6) as well as the actual expressed differences. In this example:

\[ \begin{array}{r|ccc} & \text{cell}_1 & \text{cell}_2 & \text{cell}_3 & \cdots \\ \hline \text{CD3D} & 1 & 0 & 0 \\ \text{MS4A1} & 5 & 2 & 192 \\ \text{LYZ} & 0 & 0 & 97 \\ \vdots & \vdots & \vdots & \vdots \end{array} \]

you can see how cell 3 counts 20x+ that of cell 1. Cell 3 may actually have much more total RNA or it just had a strong reaction / was sequenced "deeply." To remove this experimental effect, we can divide each count by the total UMI count in the cell. We lose the potential signal from the absolute expression, but the relative expression is still useful biologically.

\[ \begin{array}{r|ccc} & \text{cell}_1 & \text{cell}_2 & \text{cell}_3 \\ \hline \text{CD3D} & 0.17 & 0 & 0 \\ \text{MS4A1} & 0.83 & 1 & 0.66 \\ \text{LYZ} & 0 & 0 & 0.34 \end{array} \]

To put the values on a more convenient scale (instead of tiny fractions) we multiply by a constant. 10k is the default in Seurat and that works fine.

\[ \begin{array}{r|ccc} & \text{cell}_1 & \text{cell}_2 & \text{cell}_3 \\ \hline \text{CD3D} & 1700 & 0 & 0 \\ \text{MS4A1} & 8300 & 10000 & 6600 \\ \text{LYZ} & 0 & 0 & 3400 \end{array} \]

The remaining problem is that the ranges of 0 to 10000 are still too large in absolute terms for further statistical processing. In absolute terms, a count changing from \(10 \to 20\) can be drowned out by a count changing from \(100 \to 200\) in the same experiment. To map the values into a smaller absolute range we can take \(\log(1 + x)\) (the \(+1\) prevents \(\log(0)\)) to prevent this drowning out effect:

\[ \begin{array}{r|ccc} & \text{cell}_1 & \text{cell}_2 & \text{cell}_3 \\ \hline \text{CD3D} & 7.44 & 0 & 0 \\ \text{MS4A1} & 9.02 & 9.21 & 8.79 \\ \text{LYZ} & 0 & 0 & 8.13 \end{array} \]

Highly Variable Genes (HVG)

The pbmc3k dataset starts with around 32k genes, as per the matrix market file header. QC filters that down to roughly ~13k. However, a relatively large chunk of the remaining genes that pass QC don't actually vary much between cells, so they won't help us in distinguishing cell types. For example, a gene like ACTB (index 11888 in the genes.tsv of pbmc3k) has 2677 hits in the matrix file. It's in almost every cell, because it is the gene that produces actin which is needed for basic cytoskeleton construction. Before we run more computationally intensive steps like PCA, we can simply remove these low variable genes.

Principal Component Analysis (PCA)

So far we've filtered out erroneous data, normalized and removed genes that are roughly constant across cell types. What remains is ~3k cells, each with ~2k count values per variable gene. One can consider the gene counts per cell as coordinates in a ~2k dimensional space, each point in that space a cell. Could we then just compute the "distance" between each cell in that ~2k dimensional space and look for clusters? In principle yes, but the clusters would not be biologically meaningful because of the curse of dimensionality. As the number of dimensions grows, a definition of distance like the Euclidean distance

\[ \sqrt{(x_1-y_1)^2 + (x_2-y_2)^2 + \cdots} \]

will cease to be useful because if the cluster is really defined by low variance in a small handful of dimensions (e.g. T cells defined by a few T cell specific genes), that will be drowned out by high variance in all the other dimensions.

To counteract this problem, we want to construct dimensions upon which the data spreads maximally. You can think of these new dimensions/axes as weighted combinations of genes called "principal components" (PCs), as opposed to the axes just simply being the genes themselves. They represent some kind of biological mechanism that leads to correlated gene expression, for example certain cell types are defined by genes ABC with increased expression and DEF with decreased.

A helpful analogy is imagining a dataset of people's weights and heights. That's a 2 dimensional dataset where the 2 variables are correlated: if you are taller you are more likely to weigh more and vice versa. A principal component of that dataset would be an axis representing something like "size." A weighted combination of height and weight, such that bigger people have higher values and smaller people have lower values.

Importantly, each additional axis should be orthogonal to each other otherwise we'd be adding a dimension which is at least partially capturing the same variance as another axis, wasting a costly dimension. Principal Component Analysis or PCA, is the process by which we build these axes such that \(\mathrm{PC}_2\) is orthogonal to \(\mathrm{PC}_1\), \(\mathrm{PC}_3\) is orthogonal to \(\mathrm{PC}_1\) and \(\mathrm{PC}_2\) and so on. In the height/weight analogy, \(\mathrm{PC}_1\) could represent "size" and \(\mathrm{PC}_2\) could represent "lanky/stockiness." After you account for someone's overall size, people can still vary by how lanky/stocky they are: it's an orthogonal axis of variation.

Before running PCA we also scale each gene across cells (zero mean / unit variance) so genes that just happen to have bigger numeric ranges don't dominate the components. After that, the cell data points can be redefined from \(\text{cell}_1 = (\text{gene}_1, \text{gene}_2, \ldots)\) to \(\text{cell}_1 = (\mathrm{PC}_1, \mathrm{PC}_2, \ldots)\). PCA itself can generate up to \(\min(n_{\text{cells}}, n_{\text{genes}})\) dimensions; in practice you keep a small prefix (often on the order of 10–50) that captures most of the useful variance for clustering.

Clustering

With the data now living in a much smaller number of optimally chosen dimensions ("PC space"), we can cluster the cells. There are other ways to measure the distance between cells, but the simplest and default in Seurat is the Euclidean distance as defined earlier. The k-nearest neighbours approach (KNN) says to pick a value \(k\) and look for the \(k\) closest cells to each cell. \(k\) can be tuned to see if more biologically sound (i.e. the genes expressed in a cluster corresponds to a known cell type) clusters emerge or not. Empirically a starting point is 20.

A cluster of cells can then be defined roughly as cells who share a lot of KNN defined neighbours. Given 2 cells and their \(k\) nearest neighbours, one can define the "Jaccard similarity" as the number of neighbours common to both divided by the total number of unique neighbours across the 2 cells. One then weights edges between cells by this similarity metric to build a weighted graph.

Given this weighted graph, how do we label each cell with a cluster ID? Turns out this is a surprisingly difficult problem. First we must define a score that represents how good the clustering assignment is, known as "modularity" or conceptually "clumpiness." Can we just define it as the total weight of all the edges connecting cells within a cluster? No, because then a single cluster would always have the highest score. Modularity instead compares each within-cluster edge to a random baseline: for every pair of cells \(i, j\) in the same proposed cluster, take the actual edge weight \(A_{ij}\) minus the weight you'd expect if connections were random while preserving each cell's total connection strength \(k_i\). Sum those differences over same-cluster pairs (and scale by the total edge weight \(2m\)):

\[ Q = \frac{1}{2m}\sum_{ij}\left(A_{ij} - \frac{k_i k_j}{2m}\right)\delta(c_i, c_j) \]

where \(\delta(c_i, c_j)\) is 1 only when \(i\) and \(j\) share a cluster label. Positive contributions mean "more linked than chance"; the goal is to maximize \(Q\). The randomly connected weight can be scaled by a "resolution" factor; the higher the resolution, the higher above the random baseline the proposed cluster needs to be in order to move the modularity up.

A brute force approach then could be to just try all possible assignments for increasing cluster sizes and compute the modularity for each. That process grows exponentially though because for say 4 cells and 2 clusters, there are 16 or \(2^4\) different assignments. In general, \(N_{\text{clusters}}^{N_{\text{cells}}}\), which with 2700 cells, is obviously intractable.

The Leiden algorithm's approach to the computational complexity is:

  1. Start with every cell in its own cluster.
  2. Take a random cell, look at its edges and compute whether moving it into the community of a neighbouring cell would improve modularity or not. This is where the resolution parameter comes in. If the resolution is high, that means the random baseline weight is scaled up, so the edge has a higher bar to clear before triggering a merge. The effect is "merge pickiness", and probably more distinct clusters in the end.
  3. Repeat for another random cell, until local moves stop helping.
  4. Repeat the process inside each cluster. This should break apart any loosely connected clusters.
  5. Treat each cluster as a single super-cell with aggregate weighted edges to other super-cells and repeat the whole algorithm.

Leiden is not guaranteed to find the global optimum assignment, it instead finds a local optimum which is usually good enough. The randomness of selecting cells also means that different random seeds can give you different runs.

Note the Seurat actually defaults to Louvain, which is essentially Leiden without the refinement step (4) that avoids leaving poorly connected clusters.

With a set of clusters, we can now examine their gene expression to assign a cell type based on known biology. The pbmc3k dataset had already filtered out many types of other cells that would normally be in the blood like red blood cells, neutrophils, etc. The remaining set we'd expect to see are adaptive immune cells (T cells and B cells) and some innate immune cells (NK, DC, monocytes). Platelets are hard to fully isolate as well because they stick to other cells. Each of those cell populations have subtypes, such as whether the cell has fought a pathogen before or not (naive vs memory) or whether the cell has a specific protein on its surface known as its cluster of differentiation (CD). We may or may not be able to assign a cluster to a specific subtype. The pbmc tutorial, for example, ends up with 9 clusters:

In that particular dataset it just so happens that distinct clusters broke out for the subtypes of T cells and not for B cells. It's possible that by modifying the parameters we used in the analysis (\(k\), resolution, number of PCs etc.) or with a deeper biological understanding of what genes are expressed in which cell types that we'd be able to cluster in a more fine grained manner.

Assigning the labels is non-deterministic and somewhat of an art. For example, in the Seurat tutorial, how do they distinguish between naive CD4 T cells, memory CD4 T cells and CD8 T cells (a broader grouping of both)? They use a set of genes that are known to be expressed as part of specific processes:

Given these known gene expression associations, you can then reason what type of a cell a cluster is. High expression of {CD3D,IL7R,CCR7} is very likely a T-cell. If S100A4 is also present, that's evidence of a memory T-cell. If CD8 is present, that is evidence of a CD8 T cell type, but notably CD8 can also be expressed in NK cells.

Over time, researchers have externalized this type of reasoning into common sets of markers, but again most of it isn't deterministic because of the overlapping expression across cell types.

Summary

To recap, a scRNA-seq pipeline runs a sample of cells through a specialized machine to capture the gene expression in each cell. Determining the set of cell types in the sample from these raw gene counts per cell involves filtering out erroneous data (QC), normalizing for capture/sequencing depth, focusing on highly variable genes, reducing dimensions with PCA, clustering cells that share expression neighbourhoods (e.g. Leiden on a KNN/SNN graph), and finally biological reasoning to assign those patterns of gene expression to a known cell type.

Questions or comments? Email feedback@cstein.xyz