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. This post will step
through a typical pipeline and supply this missing conceptual
background.
The "hello world" of scRNA-seq datasets is a sample of ~3k
mononuclear white blood cells (Peripheral
Blood Mononuclear
Cells) from a healthy donor called "pbmc3k."
Both Seurat and Scanpy have tutorials built around this dataset, so we
can compare our results at the end.
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, the raw dataset is obtained through the following
series of steps:
- 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.
- 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.
- 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.
- 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.
- 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
a 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.
- 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.
- 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,
deduplicate UMIs then write 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 maybe 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.
Load the data
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 stastical algorithms would be prohibitively
slow on a matrix with millions of non-zero entires. 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(). We also want
to transpose the matrix with .T because for cell typing we
want cells x genes instead of what we get from Cell Ranger which is
genes x cells. Replacing the gene indexes with gene names is helpful as
well for later analysis. The loaded matrix conceptually :
$$
\begin{array}{r|ccc}
& \text{cell}_1 & \text{cell}_2 & \text{cell}_3 &
\text{...} \\
\hline
\text{CD3D} & 1 & 0 & 0 \\
\text{MS4A1} & 5 & 2 & 192 \\
\text{LYZ} & 0 & 0 & 97 \\
\text{...} &
\end{array}
$$
Quality Control (QC)
The gene counts we get may still have several sources of experimental
error:
- Doublets, as discussed in step 4.
- The droplets are not perfectly sealed. RNA can leak into empty
droplets from other broken cells.
- The sequencer isn't perfect. It can occasionally report a wrong
base, resulting in a misclassification of the gene by Cell Ranger.
- The read from the sequencer may be a short sequence that could align
with multiple places in the reference genome, known as ambiguous
alignment.
- Some cells may be dying and their gene expression is dominated by
programmed cell death (apoptosis).
To filter out these spurious readings, we do the following
(thresholds just experimentally determined / recommended by the Seurat
folks): - Drop genes seen in fewer than 3 cells. Most likely these are
sequencing/alignment errors or ambient debris. - Drop cells with fewer
than 200 genes. Likely empty cells that just happened to have some
debris in them.
- Drop cells with more than 2500 genes. These are likely doublets. -
Drop cells where there is significant (>5%) mitochondrial gene
expression (MT-*). These cells are likely in the process of dying.
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 &
\text{...} \\
\hline
\text{CD3D} & 1 & 0 & 0 \\
\text{MS4A1} & 5 & 2 & 192 \\
\text{LYZ} & 0 & 0 & 97 \\
\text{...} &
\end{array}
$$
you can see how cell 3 counts 10-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 make graphs more palatable (whole numbers instead of decimals) we
scale up 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->20 can be drowned out by a count
changing from 100->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. Its in almost
every cell, because is the gene the 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...)$$
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: say you had a dataset of the heights and weights
of a group of people. 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. Running PCA on that dataset would yield an axis
(principle component) which is a weighted combination of height and
weight, conceptually corresponding to something like "size".
Importantly, each 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 PC2 is orthogonal to PC1, PC3 is orthogonal to PC1 and PC2 and
so on. In the height/weight analogy, PC1 could represent "size" and PC2
could represent "lanky/stockiness". After you account for someones
overall size, people can still vary by how lanky/stocky they are: its an
orthogonal axis of variation.
After running PCA to determine the PCs, the cell data points can be
redefined from cell1 = (gene1_count, gene2_count...) to cell1 =
(pc1_value, pc2_value, ...). PCA itself can generate up to min(n_cells,
n_genes) dimensions, but you'd only need at most 50 of them to capture
all of the variance.
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
neighbour 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 weight edges between cells are
build 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 suprisingly 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 is actually defined in a relatively complex way: for
each pair of cells in a proposed cluster, subtract the edge weight
between them with the expected weight if cells were randomly connected
to each other using a weight corresponding to their number of actual
connections. The sum of those differences represents the total "above
random clumpiness" in the cluster. Do that for all the clusters to get a
modularity score for the full set of cluster assignments. The goal then
is to maximize modularity or "clumpiness." 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_clusters ^ N_cells - with 2700 cells that 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 adopting the label of the each
connected 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 appart 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.
With a set of clusters, we can now examine their gene expression to
assign a cell type based on known biology. In the pmbc3k dataset, the
marker genes are:
CELL_TYPE_MARKERS = {
"Naive CD4 T": ["IL7R", "CCR7", "CD3D"],
"CD14+ Mono": ["CD14", "LYZ"],
"Memory CD4 T": ["IL7R", "S100A4", "CD3D"],
"B": ["MS4A1", "CD79A"],
"CD8 T": ["CD8A", "CD3D"],
"FCGR3A+ Mono": ["FCGR3A", "MS4A7"],
"NK": ["NKG7", "GNLY"],
"DC": ["FCER1A", "CST3"],
"Platelet": ["PPBP"],
}
It is still hard to visualize the cell type clusters, with ~10 or
more PC dimesnions. Algorithms like UMAP try to compress the data
further into 2 dimensions for visualization while still retaining as
much of the structure as possible.
Questions or comments? Email feedback@cstein.xyz
Intro to Single Cell RNA Sequencing
2026 Jul 23 See all postsSingle 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. This post will step through a typical pipeline and supply this missing conceptual background.
The "hello world" of scRNA-seq datasets is a sample of ~3k mononuclear white blood cells (Peripheral Blood Mononuclear Cells) from a healthy donor called "pbmc3k." Both Seurat and Scanpy have tutorials built around this dataset, so we can compare our results at the end.
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, the raw dataset is obtained through the following series of steps:
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 maybe doublets present, though usually not a significant quantity of them to mess with the end result. It looks like this:
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.
Load the data
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 stastical algorithms would be prohibitively slow on a matrix with millions of non-zero entires. So we use scipy, which is optimized for this kind of work.
scipy.io.mmreadwill 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 withtocsr(). We also want to transpose the matrix with.Tbecause for cell typing we want cells x genes instead of what we get from Cell Ranger which is genes x cells. Replacing the gene indexes with gene names is helpful as well for later analysis. The loaded matrix conceptually :$$ \begin{array}{r|ccc} & \text{cell}_1 & \text{cell}_2 & \text{cell}_3 & \text{...} \\ \hline \text{CD3D} & 1 & 0 & 0 \\ \text{MS4A1} & 5 & 2 & 192 \\ \text{LYZ} & 0 & 0 & 97 \\ \text{...} & \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): - Drop genes seen in fewer than 3 cells. Most likely these are sequencing/alignment errors or ambient debris. - Drop cells with fewer than 200 genes. Likely empty cells that just happened to have some debris in them.
- Drop cells with more than 2500 genes. These are likely doublets. - Drop cells where there is significant (>5%) mitochondrial gene expression (MT-*). These cells are likely in the process of dying.
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 & \text{...} \\ \hline \text{CD3D} & 1 & 0 & 0 \\ \text{MS4A1} & 5 & 2 & 192 \\ \text{LYZ} & 0 & 0 & 97 \\ \text{...} & \end{array} $$
you can see how cell 3 counts 10-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 make graphs more palatable (whole numbers instead of decimals) we scale up 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->20 can be drowned out by a count changing from 100->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. Its in almost every cell, because is the gene the 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...)$$
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: say you had a dataset of the heights and weights of a group of people. 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. Running PCA on that dataset would yield an axis (principle component) which is a weighted combination of height and weight, conceptually corresponding to something like "size".
Importantly, each 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 PC2 is orthogonal to PC1, PC3 is orthogonal to PC1 and PC2 and so on. In the height/weight analogy, PC1 could represent "size" and PC2 could represent "lanky/stockiness". After you account for someones overall size, people can still vary by how lanky/stocky they are: its an orthogonal axis of variation.
After running PCA to determine the PCs, the cell data points can be redefined from cell1 = (gene1_count, gene2_count...) to cell1 = (pc1_value, pc2_value, ...). PCA itself can generate up to min(n_cells, n_genes) dimensions, but you'd only need at most 50 of them to capture all of the variance.
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 neighbour 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 weight edges between cells are build 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 suprisingly 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 is actually defined in a relatively complex way: for each pair of cells in a proposed cluster, subtract the edge weight between them with the expected weight if cells were randomly connected to each other using a weight corresponding to their number of actual connections. The sum of those differences represents the total "above random clumpiness" in the cluster. Do that for all the clusters to get a modularity score for the full set of cluster assignments. The goal then is to maximize modularity or "clumpiness." 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_clusters ^ N_cells - with 2700 cells that 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 adopting the label of the each connected 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 appart 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.
With a set of clusters, we can now examine their gene expression to assign a cell type based on known biology. In the pmbc3k dataset, the marker genes are:
Uniform Manifold Approximation and Projection (UMAP)
It is still hard to visualize the cell type clusters, with ~10 or more PC dimesnions. Algorithms like UMAP try to compress the data further into 2 dimensions for visualization while still retaining as much of the structure as possible.
Questions or comments? Email feedback@cstein.xyz