Immune cell subtype discovery from single-cell RNA-seq

Draw a tube of blood, spin it down, and the thin white layer in the middle holds T cells, B cells, NK cells, monocytes and dendritic cells, tangled together. Under a microscope, several of them look nearly identical.

Every chapter before this one told the model the correct answer during training. Logistic regression knew which genomes carried resistance. The network in Chapter 6 knew which peptides provoked a response. The first chapter raised a different possibility and then left it alone for six chapters: biologists already sort cells into types without any of that supervision. This chapter is where that gets tested properly.

The question is this. Given only a cell’s gene expression, nothing else, can an algorithm that has never seen a cell type label recover the immune cell populations a lab already knows are there?

The check is what makes the question answerable rather than just interesting. The cells used here were separated into ten populations by an entirely different method from transcriptomics. Those labels never touch the clustering. They exist only to be compared against afterwards (Figure 7.1).

One limitation needs naming at the outset. Clustering finds structure, not names. It can group T cells together without knowing to call them T cells, and it can just as easily split one real population into two or merge two distinct ones, especially where the biology itself is gradual rather than discrete. Whatever comes out still needs a human, and marker genes, the genes that set one cell type apart from the rest, to say what each group actually is.

7.1 What the labels can and cannot see

The ten populations were pulled out of one donor’s blood with antibody-coated beads, and flow cytometry came in only afterwards, to check how pure each tube was. So every label here, the whole answer key, was made with an antibody against a protein on the cell surface: CD4, CD8, CD14, CD19, CD25, CD34, CD45RA, CD45RO and CD56.

The algorithm will never see a protein. It sees RNA, counted from the tail end of each message. Proteins and messages usually agree. In these cells there are three places where they do not.

One label is close to invisible. Naive and memory T cells were separated using CD45RA versus CD45RO. Those are not two genes. They are two versions of one gene, PTPRC, made by cutting different pieces out of the same message, and those pieces sit near the front, far from the tail this chemistry reads. A 2024 study tested exactly this on tail-end data from blood cells: even with very deep sequencing, fewer than one cell in ten showed either version.

Some labels overlap by definition. The helper T cells were selected on CD4 alone. Naive, memory and regulatory T cells all carry CD4, so the helper tube holds cells that match three other labels. The two CD8 populations, cytotoxic and naive, nest in the same way, with a twist: some of the most experienced killer T cells switch CD45RA back on, so “naive” here is not purely naive.

One label is mostly something else. By the dataset’s own flow check, the CD34+ tube was 45% pure.

In Chapter 5, the gene did not follow the organism. Here the label does not fully follow the RNA.

One decision was made before any code ran, rather than after a disappointing score. The results would be scored against two answer keys. The first is the ten labels, one per tube: B cells, monocytes, NK cells, CD34+ cells and six T-cell labels, four of CD4 T cells (helper, naive, memory and regulatory) and two of CD8 T cells (cytotoxic and naive). The second is six broad lineages: B cells, monocytes, NK cells, CD34+ progenitors, CD4 T cells and CD8 T cells, so that the four CD4 T-cell tubes become one family. Both were declared in advance, so neither could be chosen later to flatter the result.

7.2 What k-means does

Each cell starts as counts for thousands of genes. Many of those genes rise and fall together, so they can be squeezed into 50 combined axes that keep the main patterns. That squeeze is called principal component analysis. Now picture every cell as a point in a space with those 50 directions. Similar cells sit close together.

k-means drops k markers into that space. Every cell joins its nearest marker. Each marker then moves to the middle of the cells that joined it. Cells switch to whichever marker is now nearest, and the process repeats until nothing moves.

That is the whole algorithm. There is no answer key anywhere. The one instruction it gets is k, the number of groups (Figure 7.2).

Chapter 5 set out to use a support vector machine, a model that searches for the widest gap between two outcomes it has been told about. k-means has no outcomes. It decides where the gaps are on its own.

It is not a toy choice for this data. The 2017 study that produced these cells clustered with k-means, on 50 principal components, and for its 68,000-cell blood sample it set k to 10 after reading an error curve. Reading a curve like that is a judgement call made after seeing the data. This chapter replaces it with a rule fixed in advance.

There are two reasons k-means closes the book. It is the plainest version of the idea, the unsupervised counterpart to Chapter 1’s logistic regression. And it drags the most important decision into the open: how many groups are there? Most methods answer that quietly, inside a setting. k-means will not start until someone answers it.

One weakness shows up at the start line. Where the markers begin can change where they finish. So every run starts from 50 different random positions and keeps the tightest result.

k-means is not what single-cell labs usually reach for. Seurat and Scanpy, two of the most widely used toolkits, default to graph-based clustering: link each cell to its nearest neighbours, then find communities more tightly linked inside than out. The version used here is called Leiden. It needs no k, only a resolution dial that does the same job less visibly. Leiden ran alongside k-means, its dial set by the same rule that picks k. If the textbook method lost to the field’s default, that would be reported too.

7.3 The rules

The rules were fixed before any code ran.

The locked drawer is a few lines of code in the script that builds the pooled table (Listing 7.1). Figure 7.3 sets out the protocol as a whole.

    counts = sp.vstack(blocks, format="csr")
    order = rng(cfg["rng_streams"]["shuffle"]).permutation(counts.shape[0])
    counts = counts[order]
    keys = np.asarray(keys)[order]
    barcodes = np.asarray(barcodes)[order]
    ids = [f"c{i:06d}" for i in range(1, counts.shape[0] + 1)]

    obs = pd.DataFrame(index=pd.Index(ids))
    var_out = pd.DataFrame({"gene_symbols": var["gene_symbols"].to_numpy()}, index=pd.Index(var.index.to_numpy()))
    pooled = ad.AnnData(X=counts, obs=obs, var=var_out)
    pooled.uns.clear()
    for column in list(pooled.obs.columns):
        del pooled.obs[column]
    pooled_path = ROOT / qcfg["pooled_file"]
    pooled_path.parent.mkdir(parents=True, exist_ok=True)
    pooled.write_h5ad(pooled_path, compression="gzip")

    labels_path = ROOT / qcfg["labels_file"]
    pd.DataFrame({"cell_id": ids, "population_key": keys, "original_barcode": barcodes}).to_csv(labels_path, index=False, lineterminator="\n")

7.4 Predictions, and how this could go wrong

I wrote down predictions before any code ran. They were not blind: others had clustered these cells before. B cells and monocytes would come back largely intact. NK cells would partly blur into cytotoxic T cells, because both run the same killing programme. The six T-cell labels would not separate cleanly, and the helper label could not separate from the populations nested inside it. The CD34+ label would mostly hold together, as it did in a 2020 re-analysis. That is odd for a tube that is 45% pure, and odd results get checked against the sequencing batch first.

I also set out the ways the chapter could go wrong, each with a check.

The first is that it might sort cells by the tube they came in. Imagine ten piles of holiday photos, each taken with a different camera. Ask someone to sort them by who is in them, and they might sort them by camera instead, because one camera’s photos are brighter and another’s are blurrier. The same risk exists here. Each cell type arrived in its own tube and was measured separately, some far more thoroughly than others. The check is the baseline in the rules. Its three numbers are the “camera settings”: they say how well each cell was measured and nothing about the cell itself. The real data has to clearly beat them. If it cannot, the groups are about the camera, not the people. A second check uses the helper T cells. That tube was collected with a broad net, so it already holds many cells from another tube, the naive T cells. Those two should be hard to separate. If the computer separates them more easily than two tubes that really are different, it is sorting by tube.

The second is that it might count wrong but sort right. The computer decides how many groups there are by looking for the cleanest splits. In blood, the cleanest splits are the big ones, so it may find three or four groups where the lab defined ten. The check is to publish its score for every group count from 2 to 20, plus one extra run told to make ten.

The third is that it never says “I am not sure”. Every cell must go into some group, even cells caught halfway through changing from one type to another. And where the computer starts is random, which can change where it ends. The check is to run it 20 times from different starting points and to compare it with Leiden.

The fourth is that a strange group is not a discovery. A group that matches no label might be a new cell type. It is just as likely to be stray cells, two cells stuck together, or dying cells. Any such group is called “unresolved” and its genes are listed. No new cell types will be claimed.

7.5 Before the key opened

Everything in this section was frozen and committed to the repository before the answer key opened.

94,655 blood cells went in. Each one was just a list of gene counts, with no names attached. The processing recipe, asked for the 1,000 most variable genes, returned 999, and I used its output as it was. The first job was to decide how many groups to make. The rule tried every number from 2 to 20 and kept the one where cells sat most clearly inside their own group (Listing 7.2).

def pick(values, scores):
    """Highest score after rounding to 4 decimal places; ties go to the smallest k or lowest resolution."""
    rounded = {v: round(s, 4) for v, s in zip(values, scores) if s is not None}
    best = max(rounded.values())
    return min(v for v, s in rounded.items() if s == best)

It picked 2. Not ten, and not six. One group of 3,521 cells full of monocyte genes, and one group of 91,134 cells holding everything else. This is what the second risk warned about: in blood, the cleanest split is a coarse one. The silhouette for two groups was 0.63. For ten groups it was 0.21.

So the extra run, the one told to make exactly ten groups, came to matter a great deal. On their marker genes, its groups look like real cell families: B cells, NK cells, monocytes, several kinds of T cells and three separate groups the naming rule called CD34+ progenitors. Leiden, the kind of method most labs use, settled on 23 groups.

Then came the 20 reruns from different random starting points. Leiden gave much the same answer each time. k-means did not. At two groups, 12 of the 20 reruns found a different split from the run fixed in advance, and a slightly tighter one, with 83,404 cells in one group and 11,251 in the other. The rules say the pre-committed run is the one that gets scored, so it is. The other answer is reported next to it.

Finally, the two groups k-means chose for itself were named from their genes, twice. The fixed rule called the big group “CD8 T cells”. The second naming, the one I had said I would do myself, worked from the same marker genes before the key opened, and called the big group “unresolved”: a mix of several cell types. The fixed rule also named every group of the ten-group run and of Leiden.

7.6 Against the answer key

Every grouping beat random easily. I shuffled the labels 1,000 times, and the highest any shuffle scored was 0.003.

That was the warm-up. The real test was the camera-settings baseline: the same method, k-means, given only the three numbers about how well each cell was measured.

At the two groups k-means chose for itself, gene expression lost. On the ten labels it scored 0.01 and the camera settings 0.04. On the six lineages it scored 0.05 and the camera settings 0.13. The reason is almost boring. Gene expression split off a small group, mostly monocytes. The camera settings split off the cells that gave the most RNA. The CD34+ tube gave far more RNA per cell than any other, so 89% of it landed in that group. Against the labels, the camera-settings split happens to score better.

The other two-group answer did better. The split that 12 of the 20 reruns found scored 0.06 on the ten labels and 0.18 on the six lineages, above the camera settings’ 0.04 and 0.13. That comparison was not part of the plan and has no interval, so the verdict still rests on the run fixed in advance.

Then I told k-means to make ten groups, the same number as the lab’s labels. This time gene expression won clearly: 0.44 against 0.08 on the ten labels, and 0.55 against 0.09 on the six lineages.

Leiden, with its 23 groups chosen by the same rule, beat two-group k-means on both keys, with 0.58 on the ten labels and 0.42 on the six lineages. Against ten-group k-means it was a split decision. Leiden scored higher on the ten labels, 0.58 against 0.44, and k-means scored higher on the six lineages, 0.55 against 0.42. That comparison was not planned in advance, so it stays a description rather than a verdict.

The result fits in two sentences. When k-means has to decide how many groups exist, it gets the count wrong, and the camera settings beat it. When it is told the count, it recovers a large part of the lab’s structure from RNA alone. The second risk in my list was “count wrong but sort right”, and that is what happened.

The names were scored too. A name counted as correct if it matched the group’s most common lineage, where that lineage made up at least half of the group, or if it said “unresolved” where no lineage did. The big group k-means chose for itself was only 46% CD4 T cells, its most common lineage, so “unresolved” was the right name and “CD8 T cells” was wrong. The fixed rule got 1 of those 2 groups right, and my second naming got both. My naming covered only those two groups, so its perfect score is a small test. On the ten-group run, the fixed rule named 9 of 10 groups correctly, and on Leiden, 17 of 23. Every group it missed on those two runs was a small one.

7.7 The six predictions

The predictions were scored on the run told to make ten groups, because they are about ten labels. Four passed. Two failed (Figure 7.4).

B cells came back intact, a pass. 99.9% of them landed in one group, and that group was 96% B cells. Monocytes also passed: 93% of them landed in one group.

The NK-cell prediction failed. It needed an overlap with cytotoxic T cells of at least 0.10, and the overlap came out at 0.05. Sharing a killing programme was not enough to mix them: 98% of NK cells sat in one group that was 96% NK cells.

The CD34+ prediction also failed, and it is the most interesting miss. The tube did not scatter into other cell types. Most of it split into three groups of its own, each more than 98% cells from the CD34+ tube. But its biggest single group held only 30% of it, far short of the 80% needed. This was the one prediction shaped by the 2020 re-analysis, where most of these cells stayed in one cluster.

The predictions said odd results get checked against the sequencing batch first, and this result is odd in that way: a tube that was 45% pure kept most of its cells in groups of its own. In the biggest of those groups, 34% of cells had at least one count of the CD34 gene, against 2% of all the cells outside it, so that group does carry the CD34 message. Among cells from the CD34+ tube, the share was 34% inside that group and 32% outside it. By this measure, the tube’s cells outside its biggest group resemble the ones inside it rather than contaminants. What the check cannot do is separate the cells from their run. The camera settings alone kept much of the tube together: 89% of it in one group at two groups, and 54% at ten. Batch and label are the same thing here, so nothing in this chapter can tell how far the tube stands apart because of its cells and how far because of its run.

The six T-cell labels did not separate cleanly, a pass. Not one of them got a group holding 80% of its cells at 80% purity.

The last prediction is the one that matters most: the helper tube could not be pulled apart from the naive T cells inside it. It passed. Helper against naive separated with a pair score of 0.05. Naive against memory, two genuinely different kinds of T cell, separated at 0.49. If the tubes, rather than the cells, had been driving the groups, that order would have flipped. It did not.

Because others had clustered these cells before, B cells and monocytes were safe calls. The two misses are the reason predictions get written down first. Without a record made before any code ran, it would be easy to claim afterwards that these results were expected.

This chapter is rewritten from the posts for Days 62 to 69 of the series.

References

Zheng, G. X. Y., Terry, J. M., Belgrader, P. et al. (2017). Massively parallel digital transcriptional profiling of single cells. Nature Communications, 8, 14049. https://doi.org/10.1038/ncomms14049