Digital Cell Sorter (DCS): a cell type identification, anomaly detection, and Hopfield landscapes toolkit for single-cell transcriptomics

Functionality overview and toolkit structure

DCS functionalities include: (i) pre-processing (handling of missing values, removing all-zero genes and cells, converting gene index to a desired convention, normalization, log-transforming); (ii) quality control and batch effects correction; (iii) cells anomaly score quantification; (iv) dimensionality reduction and clustering; (v) cell type annotation; (vi) visualization, and (vii) post-processing analysis.

We classify our tools into three categories: primary processing tools, data query tools, and visualization tools. Primary processing tools consist of functions transforming input files into valid tables and translating gene names to a desired convention, general data pre-processing, dimensionality reduction, and clustering. Cell type identification and anomaly score calculations are also included in the processing tools. Data query Application Programming Interface (API) tools include functions for retrieving gene expression of one or many genes across all cells or in a subset of cells, extraction of new marker genes characteristic of a cell cluster, and other query-type functionalities. Visualization tools include two-dimensional projection, quality control histograms, marker expression projections, marker expression summaries, gene expression heatmaps, individual gene t-texts, cell types assignment matrices, cell types stacked barplots, anomaly scores projections, pDCS null distribution histograms, new markers plots, Sankey diagrams (a.k.a. river plot), and cell type markers summary diagrams.

The DCS package is implemented in Python 3 and compatible with a variety of modern operating systems. The source code is deposited at https://github.com/sdomanskyi/DigitalCellSorter. Each new DCS release is archived on Zenodo, and can be installed via PyPI or GitHub, for example, by running command pip install DigitalCellSorter in user’s terminal. The Sphinx-build interactive documentation is available at https://digital-cell-sorter.readthedocs.io and as appendix of the Supplemental Materials. Tutorials and detailed description of each module and function of DCS as well as installation instructions for specific platforms are included in the documentation.

Cell type annotation

After unsupervised clustering, DCS automatically assigns clusters to cell types without relying on an expert curator to interpret the data, or a training data with cells already labeled. DCS uses all the information available in a knowledgebase of characteristic markers for many cell types. While cell type identification by manual interpretation generally provides good results, DCS assures that all the available information, including the presence and absence of markers, is taken into account, and can automatically identify cell types in very large datasets. Moreover, DCS provides statistical significance of the assignments, labeling as unknown clusters for which there is not sufficient information to make an assignment, and provides multiple labels with their associated scores when the expression pattern is consistent with more than one cell type. In DCS, automatic cell type annotation can be obtained using two methods: a voting algorithm and an Hopfield recurrent network classifier.

Hopfield landscapes visualization

The Hopfield model introduced in the previous section can be used to define a quasi-energy function (7) $$E=-{\displaystyle \frac{1}{2}\sum _{i,j}{\sigma }_i{J}_{ij}{\sigma }_j},$$where σ_i is the Boolean variable describing gene up or down regulation, and the matrix J is defined using Eq. (2). Eq. (7) defines a Lyapunov function for the signaling dynamics in the gene network and allows us to build a Hopfield landscape in which the attractor states, that is, the different cell types, are the minima of a complex multi-dimensional phenotypic landscape. The quasi-energy in Eq. (7) can be defined with or without the matrices A_ij and Q_μν, and using different normalizations for the ${\xi }_i^{\mu }$. Therefore, our DCS software allows for different options for the landscape representation. This landscape can be explicitly visualized by starting at the equilibrium points corresponding to the attractor configurations and adding noise to sample their basins of attraction. The points sampled can then be represented in a 2D plot using principal components projections (Szedlak et al., 2017; Maetschke & Ragan, 2014; Fard et al., 2016; Taherian Fard & Ragan, 2017). An example of output from this DCS visualization tool is in Fig. 1, where the A_ij was used (and not the Q_μν) and the ${\xi }_i^{\mu }$ were normalized as in the table $\tilde{M}$ (See Fig. S2 for a landscape with Q_μν). The figure shows the visualization of the phenotypic landscape of 14 different hematopoietic cell types, where the coordinates are the first two largest principal components of the point attractors, and the colors and contour lines reflect the Hopfield quasi-energy function of Eq. (7). These attractors were built using markers from Newman et al. (2015), modified to merge subtypes of B cells, CD4 T cells, NK cells, Macrophages, Dendritic cells and Mast cells, and their position in the landscape are indicated by stars. This visualization effectively uses the Hopfield quasi-energy to represent the matching/mismatching with the attractors in all dimensions. Note how for some cell types the basins of attraction are closer to each other than for others, and sometimes they overlap. For instance, Mast Cells, Eosinophils, and Plasma cells have overlapping basins, well separated from other cell types such as T cells or Macrophages. Additional considerations and analyses on Hopfield attractors and the role of the Q matrix are discussed in Supplemental Materials (Figs. S2 and S3).

Cell anomaly quantification

This module implements an algorithm that quantifies cell anomaly, that is, how much cells are different from other cells within the same dataset or within one or more clusters. Anomaly detection and can be interpreted as the opposite of clustering, which is based on similarity measures. The module for the quantification of anomaly is based on the Isolation Forest anomaly detection algorithm (Liu, Ting & Zhou, 2008, 2012).

The Isolation Forest algorithm isolates cells by randomly selecting a gene and a split value for its expression value. A schematic illustration of the algorithm is shown in Fig. 2. For the sake of illustration, consider a genome with two genes A and B only (along the two axes in Fig. 2A). Random partitioning is illustrated by the vertical and horizontal lines labeled with s_1,…,9. At each step the cells are partitioned in two sets by choosing one of the two genes and a random threshold for the chosen gene’s expression. For a typical (blue) cell (Fig. 2A), the number of steps necessary to isolate the cell is larger than for an anomalous cell (red). This recursive partitioning can be described by a tree, as shown in Fig. 2B, where each internal node contains a gene and a threshold value that were used in that partitioning, whereas leaf nodes are the isolated cells. The partitioning is carried out until all cells are isolated, that is, no further partition is possible. The cell anomaly algorithm has a training and an evaluation stage. In the training stage, an ensemble of n = 100 isolation trees, called isolation forest, is generated using random sub-samples of ψ = 256 cells from the original dataset. Next, in the evaluation stage, for each cell i in the dataset, or in a subset of cells of interest, and for each tree j in the forest, a path length h_ij is derived by counting the number of edges from the root node of the tree to the node where the cell is isolated by following genes and thresholds stored in tree j. To compute anomaly score of cell i, path lengths h_ij obtained from all n trees in the isolation forest are averaged and normalized according to (8) $$s_i=-{\displaystyle \frac{\sum _j{h}_{ij}}{c(\psi )n}},$$where c(ψ) is a normalization factor that accounts for the sub-sampling size ψ, giving a single score for the cell i, or a measure of anomaly. The algorithm’s pseudocode, and details about normalization factor, special conditions in the anomaly score computation and optimal selection of n and ψ are discussed in the original publication of the anomaly detection algorithm (Liu, Ting & Zhou, 2008). When a forest of random trees collectively results in shorter branches for particular cells, they are likely anomalous cells. Cell anomaly can be used to rank cells from more anomalous to less anomalous. Then the top-ranked cells can be further analyzed to investigate the biological difference from the other cells in the dataset. This module can be useful to identify small numbers of cells that may not be otherwise separated into a distinct cluster.

Quality control

Figure S3 shows an example of output from the quality control module. Count depth, that is, number of reads per gene, and gene counts, that is, number of non-zero genes per cell, are evaluated for each cell. Cell with count depths and gene counts that are below 50% of the median of their respective distributions are tagged as low quality. The quality control parameters defaults can be overwritten by a software user, as different datasets may require adjustment of the cutoffs.

The DCS quality control procedure also accounts for the fraction of mitochondrial genes expressed in each cell, as shown in Fig. S3C. The fraction of mitochondrial genes cutoff is chosen at the abscissa interception with a line that passes through the median of the distribution and median plus 1.5 standard deviations. The list of human mitochondrial genes is taken directly from MitoCarta2.0: an updated inventory of mammalian mitochondrial proteins (Calvo, Clauser & Mootha, 2016).

The integration of DCS with other algorithms for batch correction, clustering, dimensionality reduction and other visualization tools is detailed in the Supplemental Materials (Figs. S4–S7).

Cell identification in mixtures

We compared the performance of our previous pDCS method (Domanskyi et al., 2019b), with the new Hopfield classifier and the consensus annotation classifier described in Voting algorithm. To evaluate the algorithms’ performance, we randomly generated 100 mixtures of pure cell types, representing a gold standard, and evaluated the algorithms based on their automatic annotation. As gold standard, we used data from CD14 monocytes, CD19 B cells, CD34 cells, CD4 memory T cells, CD4 naive T cells, CD4 regulatory T cells, CD4 helper T cells, CD8 cytotoxic T cells, CD8 naive cytotoxic T cells and CD56 NK cells from a FACS-sorted PBMC dataset (Zheng et al., 2017), in addition to endothelial (SRS2397417, SRS4181127, SRS4181128, SRS4181129, SRS4181130) and epithelial cells (SRS2769050, SRS2769051) annotated in the PanglaoDB database (Franzén, Gan & Björkegren, 2019). The total number of cells used was 98,752, 53% of which are T cells. Each independent synthetic mixture contains on average 5,000 cells, chosen randomly from the 16 cell types listed above in random fractions. As the algorithms uses prior knowledge in the form of marker genes, we tested the algorithms using two marker genes tables: (i) CD Marker Handbook containing 11 main cell types and covering all 16 cell-types and sub-types present in the gold standard set (BD Biosciences, 2020), and (ii) CIBERSORT LM22 (Newman et al., 2015), modified to merge subsets of B cells, T cells, NK cells, Macrophages, Dendritic cells and Granulocytes. Marker table (ii) does not have marker information for CD34 cells, Endothelial cells or Epithelial cells, which are 25.6% of all the cells. To quantify each method’s performance we calculated a multi-class weighted F1 score, excluding unassigned cells, and took the median value over 100 independent random mixtures. We also calculated the median of the fraction of cells that had cell type unassigned for the same random mixtures. The results are shown in Table 1. Note how, for both marker tables, the new Hopfield classifier performed better than the original pDCS method and that the consensus annotation outperformed all other methods.

A particular case of cell identification in a mixture is shown in Fig. 3. This figure shows how the platform can be used to analyze complex mixtures involving different tissues. Here we combined scRNA-seq data from: (1) CD138+/CD38+ cells from bone marrow samples obtained from hip replacement surgery (Ledergor et al., 2018) and (2) PBMC (Zheng et al., 2017). The samples from (1) are dominated by plasma cells due to the CD138+/CD38+ pre-selection by flow cytometry, while plasma cells should be relatively rare in the PBMC sample. We therefore expect our algorithm to label most of the cells from (1) as plasma cells, and identify other cell types in the cells coming from (2). Figure 3A shows the t-SNE layout of the mixture of 8,448 CD138/+CD38+ cells and 1,000 PBMCs randomly selected cells from the PBMC dataset, after QC and batch effect correction. Clusters annotated with cell type information are shown in Fig. 3B and the relative size of the clusters are shown in Fig. 3C, including cells that did not pass QC. Annotation was obtained using the CIBERSORT LM22 list as knowledge-base. Note how in this particular case the algorithm correctly identified 96% of the cells from (1) that passed QC as “Plasma cells”, while most of the remaining was assigned to “Monocytes”. Cells from (2) were mostly identified as T cells, which is the dominant cell type expected in dataset (2). The Sankey plot in panel (D) shows on the left side the clusters with their label and on the right side the batches from datasets (1) (labeled by “hip”) and (2). The thickness of the lines is proportional to the number of cells in the corresponding cluster-batch pair. Note how most of the cells from dataset (1) are linked to clusters labeled either as plasma cells or failed QC.

Performance with different marker lists

Since our voting and Hopfield-based algorithms depend on a list of marker genes, we have studied the dependance of cell type assignment on different lists of markers by comparing the labels assigned to a gold standard. In addition to the CIBERSORT and the CD Marker Handbook introduced above, we created a new list of markers obtained by differential expression analysis of data from a large recently-published single-cell study (Han et al., 2020), and we also used a manually-curated list of cell-type markers from PanglaoDB (Franzén, Gan & Björkegren, 2019). These cell marker lists differ in the number of markers per cell type, presence/absence of negative markers, and overlap of markers across cell types. More importantly, these lists are based on different labels, often corresponding to different degrees of type/subtype refinement. As a gold standard we used a set of 68,579 PBMC cells from Zheng et al. (2017). Figure 4 shows Sankey plots connecting cells annotated using the consensus annotation and one of the markers list (on the left side of each plot) with the annotation in the gold standard (right side of each plot). This visualization allows for a direct comparison of cells labeled using different degrees of type/sub-type refinement. The plot also explicitly indicates cells that were labeled as “Unassigned” and cells that did not pass QC. Overall, a comparison among these four plots indicates that the annotation is quite consistent.

Validation of anomaly detection

We have designed an in silico experiment to validate and demonstrate the utility of the cell anomaly score introduced in Cell anomaly quantification. The experiment mimics a scenario in which a small number of circulating endothelial cells are present in peripheral blood, and evaluates the algorithm on its ability to detect these cells using their anomaly score. We selected 1,637 cells (T cells and Monocytes of PBMC) from SRS3363004 (GSM3169075) and 70 endothelial cells from SRS3822686 (GSM3402081) to mimic the presence of cells that are rare in blood. Both datasets were annotated in PanglaoDB and have been sequenced using Illumina HiSeq 2500 and 10× Chromium. The two datasets have similar sequencing depth, median number of expressed genes per cell and fraction of good quality cells, as detailed in Table S1.

We combined 1,637 PBMC and 20 endothelial cells randomly chosen out of the total 70. We then calculated the anomaly score of each cell after DCS normalization. The normalization and anomaly score calculation was repeated 100 times, each time with new randomly chosen 20 cells out of the 70 available cells. Receiver operating characteristic (ROC) curves were calculated for each of the 100 realizations with the cell anomaly rank as threshold. The average ROC curve is shown in Fig. 5A. The AUC is 0.929, indicating that the anomaly score can efficiently identify cells of a phenotype different than the majority of the other cells in the dataset. Next, instead of normal endothelial cells we used Merkel cell carcinoma cells (from SRA749327, SRS3693909) resulting in AUC of 100-fold average ROC 0.951, Fig. 5B, and Kaposi’s sarcoma cells (from SRA843432, SRS4322341) resulting in AUC of 100-fold average ROC 0.844, Fig. 5C. As a negative control, we replaced the 70 endothelial cells with 928 T cells from a bone marrow sample (from SRS3805245, GSM3396161) sequenced using Illumina HiSeq 3000 and 10× chromium, but otherwise very similar to the endothelial cells in sequencing depth and median number of expressed genes per cell. As before, we chose 20 out of the 928 T cells and combined them with 1637 blood cells. The ROC curve averaged over 100 random realizations average is shown in Fig. 5, curve (D). The AUC is 0.594, indicating that cells of similar phenotype get similar score despite being sequenced in different experiments. We have also used 5 and 10 cells instead of 20 in a set of additional in silico experiments to verify that the resulting AUC values are similar to the experiments with 20 cells detailed above.

To emphasize the difference between anomaly detection and clustering, we have quantified how close anomalous cells are to clusters of other cells. First we analyzed the normal PBMC cells to identify two clusters, L and M, with 1,231 lymphoid cells and 407 myeloid cells, respectively. Then, we combined the PBMC cells with 20 randomly chosen endothelial cells, hereon denoted as cluster O, and we normalized and projected the expression data on the principal components. Next, we calculated the inter- and intra-cluster euclidean distance for the clusters L, M and O, and derived the Silhouette scores for each of the clusters (defined as the difference between the nearest inter-cluster distance and intra-cluster distance, divided by the larger of the two) and the average Silhouette score across all cell in clusters L, M and O. Repeating this procedure 100 times, each time selecting a new set of cells O, we calculated the average measures, summarized in Table S2. The table has four sections where the small cluster O corresponds to the four cases in Fig. 5. The analysis shows that the O clusters always have to lowest silhouette score, indicating that clustering algorithms based on an euclidean distance will not separate the cells in O from the L and M cells.