The need to reassess single-cell RNA sequencing datasets: more is not always better

Preprocessing of fibroblast sample data

Fibroblast sample data originated from five donors, as described by Reynolds et al.,¹³ and were processed from raw fastq files (E-MATB-8142). The ID numbers are 4820STDY7388991 [S1], 4820STDY7388999 [S2], 4820STDY7389007 [S3], SKN8104899 [S4], SKN8105197 [S5]. Fastq files were processed using the loompy fromfq pipeline described in https://linnarssonlab.org/loompy/kallisto/index.html. Loompy (RRID:SCR_016666) and kallisto (RRID:SCR_016582) versions are 3.0.6 and 0.46.0. Genome fasta index and annotations were based on GRCh38 Gencode v31 (RRID:SCR_014966). Additionally, for other annotations and analysis of other populations, the processed h5ad adata from¹³ was downloaded from the Zenodo repository (ID: 4536165).

Each individual sample (S1–C fibroblast types and subtypes that we had just descS5) data was processed equally using the following scanpy (RRID:SCR_018139, v1.7.0rc1)¹⁴ procedure. To map the clusters from the original publication, cells from the processed data set were extracted and mapped to the samples. Genes with fewer than 30 counts were rejected. The sample was normalised (sc.pp.normalize_per_cell) and log-transformed. Then, Principal Component Analysis (PCA) with 30 components was calculated and feature selection was performed with triku¹⁵ (RRID:SCR_020977, v1.3.1), and kNN with cosine metric were computed. Finally, UMAP (RRID:SCR_018217, v0.4.6)¹⁶ and leiden (v0.8.3)¹⁷ were applied to detect the fibroblast populations.

Most of the cells from the preprocessed adata were mapped to the raw dataset. However, additional unmapped cells appear, some of them related to other cell types (e.g. keratinocytes, immune cells or perivascular cells). To assign unmapped cells to their corresponding cell types a population matching algorithm was applied (described below). This algorithm requires a dictionary of cell types and markers. The markers used were the following:

Once cell types have been assigned, non-fibroblast cells were discarded, and the PCA, triku, kNN, UMAP, leiden cycle was repeated to recalculate the new cell projection.

The sample S5 was discarded from the analysis due to its lack of SFRP2 expression, a well established fibroblast marker that is expressed in the rest of samples.¹²

Then, we separated the Fb2 population from the Fb1 and Fb3 populations for each dataset and applied the population matching algorithm to annotate them with the labels assigned from.¹² The genes used for the population assignation were the following:

Next, all datasets with Fb1 and Fb3, or Fb2 populations were joined. We applied the previous processing routine and, to correct for batch effects, we used bbknn (v1.4.0)¹⁸ with metric=angular and neighbors_within_batch=2 parameters.

To analyse the transcriptomic profile between Fb1 and Fb3, and Fb2 populations, we joined the two datasets and applied the same processing pipeline as before. We first characterised the genes driving the differences by obtaining the DEGs between the two sets of populations, and running GOEA with the first 150 DEGs of each category. The set of ontologies used was GO Biological Process 2018 with the module gseapy (v0.10.4).¹⁹ Then, to assess that the differences were due to cellular stress in the Fb2 population, we downloaded the lists of genes mentioned in the Results section (gene lists are available in the Github repository below), and genes appearing in more than two lists were selected. Then, the population matching algorithm was run against this list, and clusters with scores lower than 0.55 were assigned as ”Non-stress” clusters.

To analyse the differences in transcriptomic profiles within Fb1 and Fb3 populations, we obtained the DEGs between the two sets of A2 populations, which were the easiest to separate in clusters. By using that list of DEGs, we applied the population matching algorithm and divided the Fb1 and Fb3 populations into two halves. We then obtained the DEGs between the two halves and ran GOEA with the first 150 DEGs of each category, which revealed a hypoxia pattern in one of the halves. To assess that the differences were due to hypoxia, we downloaded the lists of hypoxia-related genes, and genes appearing in more than two lists were selected. Since some key genes (some glycolysis genes, or important genes appearing in one list) were missing, they were manually added to obtain a more robust list. Then, the population matching algorithm was run against this list, as well as the list of stress-related genes, and clusters with scores lower than 0.5 were assigned as ”Normal” clusters.

To replicate the analysis on the rest of the cell types, we used the processed h5ad file.

Correction of stress and hypoxia cell states

In order to correct for stress and hypoxia cell states we used the sc.pp.regress_out implementation from scanpy on the stress and hypoxia scores. We first created two sub-datasets, one containing stress and normal cells, and another one with hypoxia and normal cells, and then the scores were regressed out. Finally, the common processing pipeline was applied. Additional correction methods can be seen in the notebooks in the Zenodo repository.²⁰

Population matching algorithm

The aim of this algorithm is to assign a set of clusters to a set of labels, where each label contains a list of representative markers. For each label we extracted the matrix of counts of the genes belonging to the label. Then, we created a new matrix, where we assigned to each cell and gene the sum of the counts of the gene within its kNN, divided by the number of neighbours. This step reduced the noisiness of the expression, and also exacerbated the local expression of a gene and dampened the expression of sparse genes.

Gene expression values were substituted by the ranked index of their expression; and the values were divided by the largest index to sum 1. Therefore, the cell with the highest expression had a value of 1 for that gene, while the lowest expressed cell had a near 0 value. After this normalisation was applied to the rest of genes within the label, the mean of the normalised values across genes was computed, so that each cell had one value for that label.

After the previous steps were computed for the rest of labels, a new matrix with the number of clusters by the number of labels was computed. For each label and each cluster, the percentile of the normalised values within cells of that cluster was computed (percentile 70 by default). This helped reduce noise on normalised values, and assigned a unique number per cluster.

This algorithm allowed choosing of intermediate states, that is cell labels with a high similarity. By default, the label with the highest score per cluster was chosen. With the intermediate state option, labels that had a similar value as the label with the highest value were included. The difference in values was set as a threshold (0.05 by default), and labels with a difference of a value greater than the threshold were not merged.

Reassessment of the main cell populations in a large skin dataset reveals the presence of clusters with stress- and hypoxia-related gene signatures

By using an unsupervised population-matching algorithm (details in processed notebooks available online²⁰) we observed that in each of the healthy donors analysed by Reynolds et al.,¹³ at least two independent fibroblast clusters expressed signature markers of the A1, A2, B1 and B2 populations. One set of cells corresponded to the Fb2 population, and the second set corresponded to the Fb1 and Fb3 populations. A joint analysis of all donors after batch effect correction showed that the cluster duplication observed in each individual donor could be replicated jointly. We therefore assumed that some global effect should be affecting the cells, i.e. Fb2 might be a copy of Fb1+Fb3 cells, although perhaps affected by some alteration. Differential gene expression (DEG) analysis between Fb2 and Fb1+Fb3 revealed an enrichment in ontology terms associated to cell stress (e.g. unfolded protein response, regulation of apoptotic process, mRNA catabolic process). We then designed a signature gene list composed of 50 DEGs commonly associated to stress in very different scRNAseq settings (e.g. ATF3, BTG2, FOS, FOSB, GADD45B, HSPA1A/B, IER2/3, JUN, JUNB, NFKBIA, NR4A1/2, PPP1R15A, RHOB).^21-25 Using this signature, the Fb2 population over-expressed BTG2, EGR1, FOSB, IER2, SOCS3, and ZFP36, among others, indicating that these cells clustered together mainly due to cellular stress.

In a further analysis of the Fb1 and Fb3 cells, we observed that the A1, A2, B1 and B2 populations appeared to duplicate once more. A DEG analysis between each pair of duplicated populations disclosed genes in one of the split populations that were related to glycolysis (ALDOC, ENO2, GAPDH, PGK1, PDK1, PFKFB4, PYGL), cell integrity, hypoxia and apoptosis (BNIP3, BNIP3L, ANGPTL4, LOX, HILPDA); whereas the second split population over-expressed units of the mitochondrial ATPase and complex I, indicating an active oxidative metabolism. It is well known that cells under hypoxic conditions switch from aerobic to anaerobic metabolism to keep energy homeostasis within the cell.^26-28 We therefore generated a curated list of hypoxia-related genes, and managed to separate the non-hypoxic from the hypoxic group with the population-matching algorithm. Once stressed or hypoxic cells were removed on the basis of a set threshold of expression of signature genes, we mapped the main types of fibroblasts in what we termed normal cell subset of Reynolds et al. (Figure 1A). Fibroblast A1, A2 and B2 populations were independently mapped, and we also found clusters which seemingly were mixtures of previously defined populations e.g. B1/B2, A1/A2, or A2/B2. No type C fibroblasts were detected.

Figure 1. A re-analysis of the Reynolds et al. dataset in search of dermal fibroblast subpopulations reveals the presence of substantial proportions of stressed and hypoxic cells.

(A) UMAP projection of normal fibroblasts (after removal of hypoxic and stressed cell subsets) reveals conservation of some, but not all, cell types previously described in independent datasets.¹² (B) UMAP projections of fibroblast, vascular endothelium, pericyte, keratinocyte, lymphoid and APC cell populations from healthy donors, labeled to highlight hypoxic and stressed cell subpopulations as characterized by overexpression of defined gene signatures.

To understand whether the stress and hypoxic signatures were only present in fibroblast subsets or could also be traced to other populations within the Reynolds et al. dataset, we mapped the stress and hypoxia gene signatures to perivascular cell, keratinocyte, vascular endothelial cell, lymphoid cell, and antigen presenting cell (APC) clusters. In our reanalysis of healthy donors, fibroblasts, perivascular cells, keratinocytes, and vascular endothelial cells showed clear hypoxia and stress-related clusters (Figure 1B). For instance, the VE3 population, described by Reynolds et al. as increased in patients suffering from inflammatory conditions, presented a clear stress-related transcriptomic profile. On the other hand, most of the VE2 population over-expressed hypoxia-related genes. On lymphoid cells we did observe a sub-cluster of stressed Tc/Th cells but no clear hypoxic profiles. On APCs, an inflammatory macrophage cluster showed hypoxia, and the M2 and DC2 clusters showed stress-related profiles. Some of these results may be expected in physiological conditions for immune cells , but others could be attributed to sample handling.

Finally, we tested if the aforementioned stress and hypoxia related signatures were present in the previously published scRNAseq datasets of human skin.^8-11 The levels of expression of these genes were clearly higher in the Reynolds et al. dataset as compared to other available resources (Figure 2).

Figure 2. Stress and hypoxia-related signatures in published human dermal fibroblast datasets.

(A) UMAP projection of normal fibroblasts (after removal of hypoxic and stressed cell subsets) reveals conservation of some, but not all, cell types previously described in independent datasets (1). (B) UMAP projections of human dermal fibroblast subsets as defined in¹² are shown here for five published datasets,^8-11,13 and depicted by the average levels of expression of stress and hypoxia gene signatures.

Correction of stress and hypoxia signatures shows that stressed cells show a non-recoverable gene signature

Since the stress and hypoxia related expression profiles are apparent, we were interested in studying the ”reversibility” of the transcriptomic signatures, and creating a normalised dataset where hypoxic and stressed cells could merge with the normal cells, and classifying the whole dataset into the original cell types described in.¹² To this end, we applied two approaches with similar results. On the one hand, we considered cell states as batches, and applied batch effect correction with bbknn and harmony. On the other hand, we applied regression on the stress and hypoxia scores shown in Figure 2 based on the Seurat’s linear regression function implemented in scanpy. Since both approaches showed similar results, we show the results of the latter case in Figure 3.

Figure 3. Dataset merging of stress and hypoxia populations show mixed degrees of integration with the ”normal” dataset.

(A) UMAP projection of merged ”stress” and ”normal” cells. There is a low degree of integration between both cell types. (B) UMAP projection of merged ”hypoxia” and ”normal” cells. There is a high degree of integration between both cell types. (C) Unsupervised assignation of fibroblast types from (B) reveals, similar to results from Figure 1A, major fibroblast types.

To further study if stress and hypoxia transcriptomic profiles are ”recoverable”, we generated two types of datasets, one each with the stress or hypoxia cells, and another one containing normal cells. When applied the correction to the stress + normal dataset we observed that there was no integration between the two states (Figure 3A). On the other hand, there was a good integration between the hypoxia and normal cell states (Figure 3B), and the main fibroblast populations could be correctly mapped (Figure 3C). From these results we infer that the transcriptome from stressed cells is much more altered than the one from hypoxic cells, to the extent that stressed cells are in a computationally non-reversible state.