scRNA-seq Analysis of Human Skin Samples

Original paper: Solé-Boldo, L., Raddatz, G., Schütz, S. et al. Single-cell transcriptomes of the human skin reveal age-related loss of fibroblast priming. Commun Biol 3, 188 (2020). https://doi.org/10.1038/s42003-020-0922-4

The paper analyses five human skin samples of young and old donors and compares them to understand what changes in skin cells over time. The skin is taken from a region that is not exposed to the sun to avoid the effects of UV radiation and it includes around 16 thousand cells before QC.

Here, we'll focus primarily on one young sample and try to replicate the paper's findings about fibroblast subpopulations. We won't go as deep into it as the study authors do, but we'll also explore sample integration and functional annotation.

Table of contents

  1. Technical preparation
  2. Cleaning and quality control
  3. Dimensionality reduction and clustering
  4. Marker genes
  5. Functional analysis
  6. Sample integration with Harmony
  7. Closing thoughts

Technical preparation

This notebook is intended to run with a kernel located at /mnt/storage/r0978323/system/.local/share/jupyter/kernels/project_1/kernel.json. It should be provided with the submission, but either way, the purpose of it is to set up a python environment running python 3.11.11 with a collection of all the necessary packages. A requirements.txt file with these contents should work to replicate the environment:

scanpy[leiden]
anndata==0.10.9
scrublet
numpy==1.26
gprofiler-official
harmonypy

A copy of the output of running the notebook should be present in /mnt/storage/r0978323/project_data/skin_aging/output, where it could be investigated if there are problems with running it independently.

The relevant modules are imported below and a random state is picked for reproducibility purposes:

Scanpy can automatically save generated plots, though it doesn't provide a lot of fine-grained control, but only a global "settings" object. We'll increase the figure size a bit to improve image resolution:

The original files from the study have been downloaded from GEO: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE130973

Three files have been downloaded to /mnt/storage/r0978323/project_data/skin_aging, unzipped and renamed:

Note that the labels of the data files say "filtered", but that is only with respect to CellRanger. The authors explain their methods like this:

16,062 cells passed the quality control steps performed by Cell Ranger. To remove possible cell doublets, we filtered out cells with more than 7500 expressed genes, and to remove potential apoptotic cells we discarded cells with more than 5% mitochondrial reads. The application of these filters resulted in a final dataset of 15,457 single-cell transcriptomes.

As we'll see below, the "filtered" dataset contains exactly 16,062 cells, so we can be sure that these files will need some additional cleanup. We'll start by creating a local directory that symlinks to the downloaded data:

Scanpy has some trouble reading the matrix and genes with its read_10x_mtx function, so we can't use it as-is. It seems that the function expects version 3 of the CellRanger format, but the study used version 2. This seems to be a problem with gzipping in particular, but unzipping doesn't quite fix it. There are more details in this github issue: https://github.com/scverse/scanpy/issues/1916.

There's a number of proposed fixes, but we've chosen here to construct an AnnData matrix manually:

Sample barcoding is described in the individual sample pages on GEO, for instance this one for sample y1: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM3758115

the cell barcode of this samples for the processed data files is -1

The sample barcodes are ordered from 1 to 5 and they correspond to the following sample names used in the paper:

We've saved the y1 and y2 samples to files, but going forward, we'll mostly work with the 3130 cells of y1.

Cleaning and quality control

Before analyzing cell clusters, we'll perform some cleanup and doublet detection. We'll start by measuring the highest-expressed genes in the downloaded dataset:

The highest-expressed gene seems to be MALAT1, a long non-coding RNA found in the nucleus, which is what we'd expect from a single-cell experiment. (In fact, a recent bioRxiv preprint argues that it can be used as a quality measure: https://www.biorxiv.org/content/10.1101/2024.07.14.603469v2).

The other genes seem to be mostly for ribosomal proteins, and there are also two mitochondrial reads. We'll filter out most of the cells with high MT- expression, but we'll see in a bit that ribosomal genes are a bit trickier to deal with.

For now, let's ask Scrublet to predict doublets in the collection:

The automatic threshold picked by scrublet doesn't seem appropriate. We can clearly see a bimodal distribution in the right image, but the line is chosen to the right of the second peak rather than in the valley between them. A better threshold might be 0.20, so we'll pick that one, save the resulting predictions and re-plot the histogram:

We'll try to visualize doublet locations relative to cell type clusters via UMAP:

Some of the doublets are a bit spread out, but there are also concentrated clusters. We might expect that the locations are between clusters or at the edges of them. We'll remove them later by clustering with a high resolution parameter. For now, we'll proceed with the rest of the cleanup, like ignoring genes with low cell expression and removing cells with low gene expression:

We'll try to locate mitochondrial genes by looking for names starting with MT-. A high percentage of those might indicate dead cells that we should remove. Hemoglobin-heavy cells should also be removed, though as we'll see, this does not seem to be a problem in this particular sample.

We'll note that the HBEGF gene does not seem to be hemoglobin-related, according to both uniprot and the human protein atlas. The "HB" here stands for "heparin-binding", so we'll make sure not to include it in the "hb" set of genes.

We'll also detect ribosomal genes and observe how they distribute.

Plotting the cells by the number of genes they express, we can see that the bulk is located below 3500 with a spread of outliers above. We'll use 3500 as a threshold later.

Mitochondrial gene percentages are not too high, but we'll use the threshold of 5% to avoid potential apoptotic cells. There are essentially no blood cells that we can detect -- hemoglobin percentages go up to 0.4 at the maximum. Other samples do have a few, but here we don't need to do anything. Still, we'll add a filter for hb < 5 just to show that it does not change the cell count.

Ribosomal proteins are a different topic. There's a wide range of their expression, up to around 60% in some cases. The 10x genomics documentation suggests that this is not a signal of poor quality (https://kb.10xgenomics.com/hc/en-us/articles/218169723-What-fraction-of-reads-map-to-ribosomal-proteins):

Single Cell 3’ v2 libraries do contain reads mapping to ribosomal protein transcripts (Rps, Rpl). The fraction of reads varies based on cell type and overall cell health.

They note that ribosomal proteins are not the same as rRNA, these are poly-adenylated (see e.g. https://pmc.ncbi.nlm.nih.gov/articles/PMC1474067/) and do not neccesarily show quality issues:

Although both ribosomal proteins and ribosomal RNA (rRNA) make up the ribosome complex, ribosomal protein transcripts are not equivalent to ribosomal RNA (rRNA). Ribosomal protein transcript detection will not necessarily correlate with either rRNA or mitochondrial transcripts.

As we'll see later, these are distributed fairly evenly across cell types, so they should hopefully not interfere much with our analysis. We'll use highly variable genes, which is what the 10x genomics website recommends as well:

Only include in the PC analysis genes that are "highly variable". We do this by selecting the genes with the highest dispersion across the dataset and performing PCA on those genes only. The Seurat package uses a similar approach.

To try to exclude outliers, however, we'll cut off cells that include more than 50% ribosomal proteins, based on the violin plot above.

We still have the MT-CO1 gene expressed, but it's a fairly low percentage. Many of the ribosomal protein genes are still in this list, but the scanpy documentaiton suggests that this is normal:

We expect to see the “usual suspects”, i.e., mitochondrial genes, actin, ribosomal protein, MALAT1

Let's write the cleaned dataset of 2728 observations to disk, so we can rerun the next sections independently:

Dimensionality reduction and clustering

We'll start by allocating a place for visual inspection of clustering. We'll set a figdir for scanpy that will store the figures for later inspection. Particularly when generating a large amount of figures, we'll opt to store them and attach select results to markdown cells.

Now, we'll read the data from the previous section, normalize counts to a sum of 1,000 and store a record of highly variable genes to the AnnData object:

We'll use PCA to find linear combinations of genes that best express the variability in the data. This will not necessarily correspond to cell types, but it can give us guidance towards particular categories of genes.

We can see some potential marker genes that we should try visualizing. COL6A2 seems to be the biggest source of variability for the first principal component with other collagen-related genes. CST3, FTL, and a number of HLA (Human Leukocyte Antigen) genes in PC2 suggest a grouping of immune cells. PC3 includes PTPRCAP, CD3D and CD3E, and IL32, which also seem immunity-related. According to the scanpy tutorial, CD3D is a well-known marker for T-cells.

Some top marker genes suggested by the paper for different cell types include IL8, LUM, KRT1, so let's keep those in a list for visualisation as well.

For now, let's try to pick a limit of principal components to use for clustering by investigating their explained variability.

The paper picks 28 as their number of principal components, which seems somewhat reasonable, but it's hard to see a sharp reduction in in explained variability for a while longer. We'll pick principal component 45 as a point where it looks as though variability really tapers off.

This is a very subjective choice and it's hard to know what's "correct". Increasing the number of components might amplify minor distinctions, but setting it too low may miss important signals of differentiation. Tweaking the number also shifts umap cluster positions quite a lot, making it hard to visually determine correspondences. We'll leave this constant in a separate cell, so it's easy to experiment with:

With that, let's try calculating a UMAP embedding and plotting it with some of the genes we looked at above.

The separation looks promising, our principal components have been able to find some genes specific to major clusters. The suggested marker genes by the paper also tell us our results are compatible with theirs so far.

Let's plot the QC-related genes and investigate their distributions across clusters:

Mitochondrial genes are evenly spread out and limited to below 5%, since that was the threshold we picked. There are some outliers with higher gene counts, though some may correspond to doublets and will likely go away when we remove them.

Ribosomal genes have a much higher percentage of expression, as we noted before, but they are mostly spread out. There is one broad cluster in the center that has a fairly high expression and one smaller on the upper-right that in a previous chart corresponded to KRT1. As we'll see, we can find non-ribosomal marker genes for both of those.

At this stage, we might also take a look at doublets in the filtered and normalized datasets to see how they're distributed.

We can see some predicted doublets sprinkled between the center and upper-left clusters, which might be an indication of doublets between these two cell types. A few are present inside of islands, where we might expect the clusters to separate with a higher leiden resolution.

Speaking of which, let's try to run the leiden clustering algorithm with different resolutions and decide which ones we should try to focus on. We'll prepare a simple function that runs it on various resolutions and stores the images for later inspection:

Let's pick a range of values with varying granularity. We'll pay particular attention to the 0.5 to 0.9 range:

The three resolutions we'll focus on are 0.1, 0.9, and 5. The last one is only interesting for picking out doublet clusters with some level of precision. The one on the left has 7 clusters that identify the major subtypes of cells and we should be able to make easy inferences about what kinds of cells they represent. Its number will drop to 6 once we remove some doublets and recalculate.

The middle clustering uses a resolution of 0.9 and is the finest resolution by 0.1 increments that separates the island in the top right (what we'll later show is keratinocytes) into two parts. It also produces a separation of fibroblasts into four clusters and one outlier, which mostly corresponds to the organization the paper settles on. (Though removing doublets seems to produce another axis of separation on what is currently cluster 5).

Let's handle potential doublets by filtering out clusters above a certain threshold of predicted doublets. We'll use the function from our course tutorials with some minor stylistic modifications:

Predicted doublets seem to be mostly spread out across clusters, so the removal of entire clusters may not be too useful. Still, we'll drop the top 3 clusters with a predicted doublet rate above 0.1 and filter out the other cells individiually:

Let's rerun clustering on the cleaned dataset:

Let's also save a snapshot of this anndata object so we can rerun the next section without repeating the clustering process, if we need to.

Marker genes

We'll start from the broad clusters defined by the 0.1 resolution and try to understand the functionality of these 7 groups of cells. Let's set up some configuration for this section and then define a class that will let us generate lists and query for particular genes.

We can calculate gene markers by using Wilcoxon tests of differential expression between the clusters. We'll pick a log-fold threshold of 2 for the up-regulated genes and we'll require certain levels of in-group and out-group fractions. We'll also keep the full list of differentially expressed genes for GOrilla.

We'll process the ranking inside of the constructor and then provide a gene_clusters method to look up where certain genes are located. We'll also add a simple gprofiler lookup method to enable us to make simple functional enrichment queries.

Broad clustering

There are 6 clusters at the 0.1 resolution that we can try to understand by using g:Profiler information. For instance, here's a list of biological processes that are upregulated in cluster 0 and a list of possible cell types from the Human Protein Atlas:

Based on these biological processes, it's clear that the cluster is made up of immune cells: lymphocytes, T-cells, or other kinds of leukocytes. We know that this is a skin sample and not one from the tonsils or spleen, but the "lymph node" suggestion by the Human Protein Atlas (HPA) supports the guess of "lymphocites". Since we're exploring this cell map in broad strokes, we can call this cluster "Immunity 1" (since there will be one more).

Let's take a look at cluster 2:

These might be smooth muscle cells, possibly the ones lining blood vessels, or maybe stem cells that can differentiate into them.

With this in mind, we can rerun the gprofile queries on all clusters and attach some very rough annotations with lists of processes that seem relevant. We'll save the results into files and add the annotations manually. In the comments, we'll quote select biological processes from the g:Profiler CSVs.

Narrow clustering

Let's start with a few statements from the paper to try to map our fine-grained clusters to their descriptions.

Keratinocytes were detected in three clusters (#5, #7 and #15) and their diversity was mainly due to their degree of differentiation. While epidermal stem cells (EpSC) and other undifferentiated progenitors (#7 and #15) expressed markers such as KRT5, KRT14, TP63, ITGA6, and ITGB1, differentiated keratinocytes (#5) were defined by KRT1, KRT10, SBSN, and KRTDAP expression

We'll use the 0.9 resolution and see if we can find these specific genes in our rankings. We'll re-plot our umap clustering to try to locate them visually.

Cluster 7 on this chart was part of what we originally classifed as "Development 1 (Epidermis)", but now we see that this is because the larger cluster includes undifferentiated stem cells, whose signal hides the presence of keratinocytes. All genes suggested by the paper are present in that cluster with very high log fold changes of 6-7.

Cluster 15 we flagged as development-related, so it makes sense that it's a group of undifferentiated stem cells. Cluster 18 was previously part of the smooth muscle cell group, but we did notice some developmental functions in our initial g:Profile exploration, which might be explained by this group being absorbed into it due to the broad leiden resolution.

We can see that both clusters 7 and 8 express KRT5 and TP63, which the paper suggests should be marker genes that separate the two. The log fold changes are higher for 8, and so are the adjusted p-values, but it might also be the case that integrating several samples or making the resolution finer might make the thresholds cleaner.

The main focus on the study is fibroblasts, so it's interesting to explore how they are clustered in our UMAP. Fibroblasts take up the entire top left cluster, which in our case is split into 5 clusters (and an extra outlier) instead of 4. We'll use the marker genes listed by the paper:

Fibroblasts were identified by their archetypal markers LUM, DCN, VIM, PDGFRA, and COL1A227

Instead of four clusters, we end up with six, whose gene expression log fold changes do not seem particularly different. More than that, only two of the genes seem to show up in the upregulated list.

We can expand our criteria for log fold changes and look for "all" genes with a log fold change above 1:

Cluster 10 looks like an outlier in log fold expression of LUM and DCN. We might be looking at another undetected doublet cluster or a different kind of issue. We'll exclude it in further analysis, since it's also a comparatively small group of cells.

Let's save the fibroblast clusters' marker genes to files, excluding 10. We'll continue with a more detailed functional analysis on them in the g:Profiler website, and submit the full list of differentially-expressed genes to GOrilla to avoid artefacts of our arbitrary log-fold threshold.

Functional analysis

We'll go over the individual clusters one by one and investigate the biological processes associated with them in both g:Profiler and GOrilla. To guide us in the right direction, we'll consult with the paper's analysis and try to map our clusters to their descriptions and see what we can figure out for the extra cluster.

The way they categorize the fibroblast subpopulations is:

Let's see if we can align the clusters we have to that categorization.

Cluster 5: Pro-inflammatory

We'll start with cluster 5, which seems to have a notable immune response functionality. We can see the relevant GO:BP terms in the g:Profiler screenshot and the relevant cluster in the GOrilla graph:

Cluster g:Profiler GOrilla
fibroblasts_05.png 2025-01-02-141858_1920x1080_scrot.png 2025-01-02-142523_1920x1080_scrot.png

We can also take note of the RelA transcription factor here, otherwise known as TP65, part of the NF-kappa-B complex, which is "Essential for cytokine gene expression in T-cells":

2025-01-02-142844_1920x1080_scrot.png

iRegulon also predicts that RelA binds 45 of the upregulated targets in this group: 2025-01-04-142013_1920x1080_scrot.png

We can expand a few of the GOrilla clusters to see mostly chemotaxis-related proteins, which according to the study has to do with positioning leukocytes:

Significant examples include inflammatory response, cell chemotaxis or negative regulation of cell proliferation, necessary for the final anchoring of leukocytes

2025-01-02-143301_1920x1080_scrot.png 2025-01-02-143436_1920x1080_scrot.png 2025-01-02-143448_1920x1080_scrot.png

It's a very nice match to the paper's cluster 2, down to the fairly low significance values compared to the other clusters. The p-values do cross a significance threshold, and, as the paper notes, this cluster does not seem enriched in "standard" fibroblast processes like building the extracellular matrix.

Cluster 11: Mesenchymal

Let's give a similar treatment to cluster 11. We can see that we have very high-confidence enrichment for "system development" and "developmental process", which strongly suggests that these are mesenchymal cells, able to differentiate into multiple different cell types.

Cluster g:Profiler GOrilla
fibroblasts_11.png 2025-01-02-151512_1920x1080_scrot.png 2025-01-02-151359_1920x1080_scrot.png

While they do express genes related to extracellular matrix organisation, the main functionality of fibroblasts, they differentiate from the others based on this strong developmental component. The g:Profiler results also include the largest number of transcription factors associated with expression in these cell types:

2025-01-02-185603_1920x1080_scrot.png

iRegulon suggests somewhat different factors, though several of them seem development-related like TCF12, involved in neuronal differentiation, EP300, a chromatin-remodeling protein, and NFIC, reportedly "involved in cellular differentiation and stem cell maintenance":

2025-01-04-143655_1920x1080_scrot.png

Cluster 13: Secretory-reticular

Going forward, it'll be more difficult to find associated biological processes that differentiate the clusters. Looking at clusters 13 and 17, they both score highly significant hits into for ECM organisation and collagen binding. We can see that cluster 13 has some mesenchymal functions, which might suggest the "reticular" subtype.

Cluster g:Profiler GOrilla
fibroblasts_13.png 2025-01-02-152808_1920x1080_scrot.png 2025-01-02-153219_1920x1080_scrot.png

At this point, the authors of the paper seem to have focused more on known marker genes rather than on GO annotations. Cluster 13 has several marker genes with a high log fold change compared to the other clusters that indicate its reticular subtype: WISP2, SLPI, CTHRC1.

Cluster 17: Secretory-papillary

Likewise, cluster 17 is mostly focused on the ECM and collagen binding, with some developmental processes mentioned, but also some processes closer to the outer layer of the skin like wound healing.

Cluster g:Profiler GOrilla
fibroblasts_17.png 2025-01-02-154346_1920x1080_scrot.png 2025-01-02-154849_1920x1080_scrot.png

The marker genes upregulated in this cluster include APCDD1, ID1, COL18A1.

Cluster 6: Signaling and/or neuronal differentiation?

So how do we sort cluster 6? So far, we've been able to associate each cluster with one in the paper, either based on the functional annotation or based on specific marker genes. Deciphering the role of this group of cells will be a bit harder. We can start with some of the best-supported biological processes suggested by g:Profiler:

Cluster g:Profiler
fibroblasts_06.png 2025-01-02-160629_1920x1080_scrot.png

What we can see at the top is an immunity regulation function, but less significant than in the pro-inflammatory cluster 5. It's also mixed with ECM organisation regulation and developmental processes.

GOrilla shows several different clusters of functionality with low, but significant p-values:

Regulation Catabolic processes ECM and developmental
2025-01-02-161046_1920x1080_scrot.png 2025-01-02-161102_1920x1080_scrot.png 2025-01-02-161109_1920x1080_scrot.png

In terms of marker genes, APOE is a protein expressed in the central and peripheral nervous system (uniprot), PTGDE is present in the tissues of the blood-brain barrier and in leukocytes (uniprot), and the highest-expressed protein by log fold change, DPT, "May serve as a communication link between the dermal fibroblast cell surface and its extracellular matrix environment" (uniprot).

Given the broad spread of functionality, we might speculate that this is a type of signaling cell or possibly a multipotent cell that is able to differentiate into neurons. A KEGG pathway reference that has relatively high support for this cell type is "complement and coagulation cascades": KEGG. This involves both inflammation and ECM changes and might be coordinated by these kinds of cells. Their catabolic functionality may be coming from a larger ribosomal protein expression:

2025-01-02-163421_1920x1080_scrot.png

iRegulon finds support for TCF12 and EP300, both present in the mesenchymal cluster. It also suggests STAT4 (which is added to the network as "CAT") that binds 41 targets and is categorized as a "Signal transducer" that is "mainly expressed in hematopoietic cells that plays a critical role in cellular growth, differentiation and immune response".

2025-01-04-145251_1920x1080_scrot.png

It's hard to say anything for sure if we're not domain experts. The cell type could also be an artefact of the particular choice of clustering resolution and an abnormal amount of ribosomal protein capture. If we decrease the granularity of clustering, the group splits between clusters 5 (pro-inflammatory) and 11 (mesenchymal). This does lend support to the multipotency of these cells and explains its participation in the aforementioned KEGG pathway.

Sample integration with Harmony

We've done this analysis on a single sample, but the paper merges the two young samples together. Let's do the same and see if cluster 6 turns out to be a feature of the sample and disappears, or if something stranger happens.

We'll read the y1 data we cleaned and marked with doublets and we'll run through the steps to clean up the y2 samples.

Now that we have two cleaned datasets with samples, we can merge them and perform the following normalization, scaling, and clustering steps to the combined dataset:

We've calculated a UMAP embedding before applying harmony, so we can plot how the two samples fit together in terms of embedding density:

There are regions of low density, though the area that we'll later confirm represents the fibroblasts is mostly well-covered by both samples. Let's integrate the principal components with harmony regardless:

The samples' datapoints seem better mixed now. Let's rerun the leiden clustering function on this new sample and visualize a few of the resolutions:

While the 0.1 resolution shapes up well with a few extra groups showing up, separating at 0.9 ends up splitting fibroblasts into a number of additional groups. Let's clean up doublets to avoid one source of noise:

Unfortunately, this doesn't seem to help. Note that cluster 6 is the previous clusters 7 and 8 made up of keratinocytes. We could broaden the resolution to try to maintain fibroblasts into 4-5 clusters, but even at this resolution, keratinocytes can't be separated anymore.

We can confirm that the top cluster consists of fibroblasts by using our MarkerGenes class to search for the same marker genes as before:

Both marker genes and biological processes from GO suggest that the top cluster is, in fact, made up of fibroblasts. So why are they split into so many groups now?

One possibility that has to be acknowledged is some technical error of ours in mixing the samples, either through harmony or through the selection of principal components, QC thresholds, etc. It's hard to see where we've gone wrong, but it's certainly not impossible.

Provided we've performed the steps correctly, it could be that fibroblasts simply vary too much between samples compared to all the other cell types. Fibroblasts are a flexible cell type with the ability to differentiate into mesenchymal cells or into epithelial cells, as described by Thiery and Sleeman in "Complex networks orchestrate epithelial–mesenchymal transitions". Between-sample variability may be detectable by leiden clustering even after adjusting for batch effects with Harmony.

It's hard to say why the paper wouldn't run into similar issues when mixing replicates, however. This is a description of their fibroblast clustering:

For the second-level clustering of the fibroblasts (Supplementary Figs. 4a, b) we subsetted the cells and ran again the functions FindNeighbors() and FindClusters(). In this case, we used 20 PCA dimensions for both functions and a resolution of 0.5 for the FindClusters() function. For visualization we re-calculated the UMAP plot with RunUMAP() function with default parameters and using 20 PCA dimensions.

They reduced their number of principal components to 20, but this produces more or less the same results for us. The Seurat functions use a different algorithm, "shared nearest neighbour modularity optimisation", but we'd have to redo the analysis by using Seurat to make the comparison. No easy answer to be found at this time.

Closing thoughts

This skin sample analysis by Solé-Boldo et al is fascinating, showing the complex variation in cell roles in this outermost layer of the human body, our main contact with the outside world. Fibroblasts show remarkable flexibility and the capacity to take multiple roles: organizing skin tissues, repairing wounds, or differentiating into different types of cells.

While we had some trouble with sample integration, we successfully reproduced the four fibroblast subpopulations that the paper describes and identified them with good support from several different functional enrichment services. Both GO function annotation and transcription factor enrichment were useful in drawing the lines between cell types. And when all else fails, leaning on the paper for hints about specific marker genes worked out well.

The extra cluster that we found could be a mix of cells from the mesenchymal and pro-inflammatory groups, or it could have its own distinct role. In the first place, our taxonomy of types is mostly there to simplify things for our own inferences, and it's unsurprising that cells would rarely fit neatly into each box. Biology is messy, but as long as we accept that and apply our best efforts, we can draw out interesting conclusions both to enhance our understanding, and hopefully to improve our health and well-being.