Quantitative Modeling of Stemness in Single-Cell RNA Sequencing Data: A Nonlinear One-Class Support Vector Machine Method

Abstract

Intratumoral heterogeneity and the presence of cancer stem cells are challenging issues in cancer therapy. An appropriate quantification of the stemness of individual cells for assessing the potential for self-renewal and differentiation from the cell of origin can define a measurement for quantifying different cell states, which is important in understanding the dynamics of cancer evolution, and might further provide possible targeted therapies aimed at tumor stem cells. Nevertheless, it is usually difficult to quantify the stemness of a cell based on molecular information associated with the cell. In this study, we proposed a stemness definition method with one-class Hadamard kernel support vector machine (OCHSVM) based on single-cell RNA sequencing (scRNA-seq) data. Applications of the proposed OCHSVM stemness are assessed by various data sets, including preimplantation embryo cells, induced pluripotent stem cells, or tumor cells. We further compared the OCHSVM model with state-of-the-art methods CytoTRACE, one-class logistic regression, or one-class SVM methods with different kernels. The computational results demonstrate that the OCHSVM method is more suitable for stemness identification using scRNA-seq data.

1. INTRODUCTION

Intratumoral heterogeneity is among the greatest challenges in precision cancer therapy. Tumors are extremely heterogeneous and contain various subpopulation cells that response differently to cancer therapy. There is a population of cancer stem cells that are highly resistant to current therapeutic approaches. The cancer stem cells persist during therapy, contribute to relapse, develop to mix type of cancers, and are capable of recapitulating the heterogeneous features in a tumor (Damavandi et al., 2023; Le et al., 2018). The evolution of the stemness features of cancer cells is crucial in predicting the effects of cancer therapy (Bao, 2021). However, the determination of tumor stemness is technically challenging. A major hurdle has been the lack of clear methods to quantify the stemness of the tumor environment or individual cells.

The concept of cancer stem cell (CSC) was proposed four decades ago, which states that tumor growth, analogous to the renewal of healthy tissues, was driven by a small number of dedicated stem cells. Many reports on the identification of CSCs have been seen in many cancer types in the last decades (Batlle and Clevers, 2017; Qin et al., 2023). The model of CSC has been developing in these years, and it has gradually become clear that CSCs do not necessarily have to be quiescent; many examples show that they can be abundant and can proliferate vigorously (Batlle and Clevers, 2017; Walcher et al., 2020). CSCs have important therapeutic implications because of their association with cell plasticity, distant metastasis, and drug resistance (Cao et al., 2021; Shibue and Weinberg, 2017; Yejing et al., 2017).

Cancer cells can stochastically transit to different phenotypes toward a dynamical equilibrium in a cellular population (Friedmann-Morvinski and Verma, 2014; Gupta et al., 2010; Su et al., 2017); the phenotypic transition can be a continuous spectrum revealed by single-cell RNA sequencing (Jiang et al., 2021; Macaulay et al., 2016). In Miranda et al. (2018), cancer stemness was calculated using GSVA package to explore the strong correlation with mutation load, cancer-testis antigen expression, and intratumoral heterogeneity. To better understand cancer development, it is important to quantify the ability of self-renewal and differentiation of cancer cells (Tu et al., 2019; Zheng et al., 2021b), that is, quantification of the stemness (Alqarni et al., 2021).

Advances in single-cell sequencing offer an opportunity to study tumor heterogeneity at the resolution of individual cells (Nawy, 2014; Ren et al., 2018). The single-cell RNA sequencing data enable us to better understand the evolution dynamics, including the identification of diverse cellular subpopulations based on variance in gene expressions (Grün et al., 2015; Jiang et al., 2016), the establishment of the cell lineage relationship (Bendall et al., 2014; Farrell et al., 2018), exploring the diversity of infiltrating T cells (Li et al., 2019; Zheng et al., 2017), and the elucidation of cancer development (Ellsworth et al., 2017; Navin et al., 2011). It is natural to ask whether we can quantify the stemness of cancer cells based on single-cell transcriptome data.

In a recent study of therapy-resistant stem cells in chronic myeloid leukemia (CML) (Giustacchini et al., 2017), the authors applied single-cell RNA sequencing to characterize cells that persist in individuals with CML who were treated with tyrosine kinase inhibitors (TKIs). Results show that cellular differences exist in different patients and different CSC populations, which can assist in a further understanding on the stem cell status evolution and tumor heterogeneity during cancer therapy. These observations motivate our study to develop a method to quantify the stemness of cancer cells, and to investigate the evolution of cancer stemness during cancer therapy.

In the recent years, a lot of transcriptome-based methods have been developed for analyzing the potential of cells for self-renewal and differentiation. A core regulatory network was reconstructed based on single-cell RNA-seq data. Through mathematical formulation describing the relationship among the critical genes “KLF8,” “ZEB1,” and “CDH1” (Tu et al., 2019), authors revealed the molecular mechanism of cell fate decisions. StemSC (Zheng et al., 2021b) proposed a stemness index based on relative expression orderings (REOs) of gene pairs. Through selection of genes significantly related to a differentiation process, reference REOs are constructed for evaluation of stemness in cells. It is a knowledge-based stemness evaluation method. In Zheng et al. (2021a), an absolute REO-based stemness index was proposed and shown to be highly correlated with tumor grade and has prognostic ability.

Attention has been paid on developing machine learning approaches for quantifying stemness. In Malta et al. (2018), a one-class logistic regression (OCLR)-based method was proposed to quantify the stemness for features associated with cancer cell dedifferentiation (Sokolov et al., 2016). The OCLR approach can be applied to different data sets, such as transcription, DNA methylation, or epigenetic regulation, and hence can define different type stemness indices. Application of the stemness indices to single-cell data revealed patterns of intratumor molecular heterogeneity, and the dedifferentiated oncogenic phenotype was generally most prominent in metastatic tumors. The OCLR method, comparable with linear support vector machines (SVMs) (Schölkopf et al., 2000), can only describe the linear relationship between different data samples. When the data are nonlinearly related as in many cellular subpopulations, nonlinear methods can be more appropriate to quantify the stemness at the single-cell level.

CytoTRACE (Gulati et al., 2020) was another recently proposed unsupervised learning framework for cellular differentiation evaluation. However, it can only provide a relative degree of differentiation status for cells in specific single-cell RNA sequencing (scRNA-seq) data sets.

In this study, we proposed a nonlinear machine learning method to define the stemness index based on single-cell transcriptome data: the one-class SVM with Hadamard kernel (OCHSVM) method. The Hadamard kernel function (Jiang et al., 2017) for modeling the nonlinear relationship between cells was effectively incorporated in the one-class SVM learning framework. We applied the proposed method to analyze the single-cell RNA sequencing in CML after TKI therapy (Giustacchini et al., 2017), and explore evolution of tumor stemness in response to the therapy. Results show that OCHSVM is robust in terms of cross-validation, and provides a reasonable stemness measure for modeling the differentiation degree of the considered single cells. Stem cells are also present in embryo and various tissues. To validate the proposed stemness index, we further test the OCHSVM stemness with data sets for mammalian embryo stem cells (ESCs) and induced pluripotent stem cells (iPSCs).

The rest of the article is organized as follows. First, we present the framework of one-class SVM and introduced the Hadamard kernel to define the stemness index. The OCHSVM model is trained with normal hematopoietic stem cells (HSCs), and then validated through quantitative performance measures and biological interpretations. Next, we applied the proposed stemness to CSCs obtained from CML patients, and analyzed the stemness distribution for the CML single cells throughout the disease course. We further test the proposed stemness with data sets from mammalian ESCs, iPSCs, and differentiation of human embryonic stem cells, which show the features of distinct different types of cells. Finally, we give the concluding remarks with comparisons between different state-of-the-art methods.

2. METHODS

2.1. Formulation of the stemness

Given a gene expression vector $x$ of a cell, we try to define a stemness index $S_{i n d}$ ( $0 < S_{i n d} < 1$ ) that gives a likelihood of the cell to be a stem cell. To this end, given a set of training stem cells $X = {x_{1}, x_{2}, \dots, x_{n}},$ we define a support function $f (x)$ that captures the training data set, that is, most of the training cells' gene expression vectors live in the tight region $f (x) > 0$ in the input space.

To this end, we try to find a hyperplane (the “maximum-margin hyperplane”) $w^{T} \cdot x - b = 0$ that separates all the data points from the origin with maximum margin and tightly supports the training data set at the same time. This leads to the following optimization problem: $\{\begin{matrix} m i n i m i z e_{w, b} \frac{1}{2} ∥ w ∥^{2} - b \\ s u b j e c t t o (w^{T} \cdot x_{i} - b) \geq 0, \\ \forall i = 1, \dots, n . \end{matrix}$ (1)

Here $x_{i} \in X$ represents the i-th stem cell in the training data set. The parameters $w$ and b solving this problem determine the support function as $f (x) = w^{T} \cdot x - b,$ which is positive for all examples contained in the training set.

To extend the above hard margin to cases in which the data are not linearly separable, we introduced hinge loss parameters $ξ_{i}$ , and the optimization problem [Eq. (1)] becomes $\{\begin{matrix} m i n i m i z e_{w, ξ, b} \frac{1}{2} | | w {| |}^{2} + \frac{1}{v n} \sum_{i = 1}^{n} ξ_{i} - b \\ s u b j e c t t o w \cdot x_{i} - (b - ξ_{i}) \geq 0, ξ_{i} \geq 0 . \end{matrix}$ (2)

Here v is a parameter that controls the fraction of training data that fall outside of the hypersphere, and $ξ_{i}$ are the slack variables whose values indicate how far these outliers deviate from the surface of the hypersphere.

Through the Lagrangian multiplier method, the problem becomes $\{\begin{matrix} {max}_{α} - \frac{1}{2} \sum_{i, j = 1}^{n} α_{i} α_{j} x_{i}^{T} \cdot x_{j} \\ s u b j e c t t o \frac{1}{v n} \geq α_{i} \geq 0, \sum_{i = 1}^{n} α_{i} = 1 . \end{matrix}$ (3)

In the case of nonlinear data separation, we replaced the bilinear function $x_{i}^{T} \cdot x_{j}$ in Equation (3) with a nonlinear kernel function $K (x_{i}, x_{j})$ (Wang et al., 2014). For example, we can take the linear kernel $K (x_{i}, x_{j}) = x_{i}^{T} \cdot x_{j}$ , polynomial kernel $K (x_{i}, x_{j}) = {(x_{i}^{T} \cdot x_{j} + 1)}^{d}$ , or the Radial Basis Function (RBF) kernel $K (x_{i}, x_{j}) = exp (- \frac{∥ x_{i} - x_{j} ∥^{2}}{σ^{2}})$ (Jiang et al., 2018; Lai et al., 2006). Recently, a Hadamard kernel was proposed and successfully applied to gene expression-based breast cancer outcome predictions (Jiang et al., 2017). The Hadamard kernel is defined as follows:

where $γ > 0$ is an adjustable parameter. From the study in Jiang et al. (2017), the Hadamard kernel is nonlinear and less likely to be stuck into the overfitting problems, and is robust when the optimal parameter locates near the origin.

Here, we applied the Hadamard kernel function, and defined a stemness index-based one-class SVM with Hadamard kernel (OCHSVM) $S_{i n d} (x) = m i n \{\frac{g_{γ} (x)}{{max}_{1 \leq j \leq n} g_{γ} (x_{j})}, 1\},$ (4)

where $g_{γ} (x) = \sum_{i = 1}^{n} α_{i} K_{γ} (x_{i}, x),$ and $α_{i}$ is the solution to the problem

The intuition behind the definition is as follows. $g_{γ} (x) = \sum_{i = 1}^{n} α_{i} K_{γ} (x_{i}, x),$ is largely influenced by the attribute space of $x_{i}, i = 1, 2, \dots, n$ and also the parameter $γ$ involved. Based on the assumptions that the stemness index $S_{i n d}$ ( $0 \leq S_{i n d} \leq 1$ ) gives the likelihood of the cell to be a stem cell, we skated round the problem by considering $S_{i n d} (x) = m i n \{\frac{g_{γ} (x)}{{max}_{1 \leq j \leq n} g_{γ} (x_{j})}, 1\} .$ In the OCHSVM stemness, the parameter $γ$ should be determined beforehand, and the parameters $α_{i}$ are obtained by solving the problem [Eq. (5)] based on the train data set of predetermined stem cells. For each gene expression vector given by scRNA-seq data $x$ , the OCHSVM stemness index satisfies $0 \leq S_{i n d} \leq 1$ , and a larger $S_{i n d}$ means more similar to a stem cell.

2.2. Identification of the parameter

γ

To determine the adjustable parameter $γ$ in the Hadamard kernel, we selected a data set of stem cells for model training, and calculated the overall stemness of the training set with varying $γ$ . Here, we selected a data set (GSE76312) of scRNA-seq data of normal bone marrow cells for model training (Giustacchini et al., 2017). In this data set, the HSC marker CD34 shows obvious bimodal distribution (Fig. 1A). There are two subpopulation cells with either high or low CD34 expression. We selected the subpopulation with a CD34 expression higher than 4.5 as the training data set.

FIG. 1.

Determine the parameter $γ$ in the Hadamard kernel. (A) Distribution of CD34 expression levels in the HSCs from the data set GSE76312. Here the expression levels are given by $log (R P K M + 1)$ . (B) Stemness index calculated from OCHSVM with different values of the parameter $γ$ . We selected CD34-positive bone marrow cells as the training data set for HSCs. Dots show the average stemness index for the testing cells with various parameters $γ$ , and 5 runs for each value. Squares show the mean of 5 runs. The optimal value at $γ = 0.7$ is marked by the arrow. HSC, hematopoietic stem cell.

There are 134 cells in the training data set. For a given parameter $γ$ , we conducted one-class Hadamard SVM training on 80% of cells, and the other 20% of cells were left for testing. We varied $γ$ from 0.1 to 1, and performed 5 independent runs with random permutation of training cells, and evaluated the stemness indexes on the testing cells. Figure 1B shows the dependence of average stemness indexes for different values of $γ$ over different runs. The average stemness is larger than 0.75 with $γ$ varying from 0.1 to 1, and reaches a maximum value of 0.87 at $γ = 0.7$ . Therefore, we can take the optimal value $γ = 0.7$ in the kernel function hereafter. Thus, the parameter $γ = 0.7$ and the training data set of normal HSCs give a definition of OCHSVM stemness through Equation (4).

3. RESULTS

3.1. Stemness of tumor stem cells in CML

We calculated the stemness of cancer cells from patients with CML after TKI treatment. CML is mainly caused by the presence of the BCR-ABL fusion gene and is often treated by TKI (Pasic and Lipton, 2017). The existence of CML stem cells (CML-SCs) is the main reason for the resistance to TKI therapy (Chu et al., 2011; Paolo et al., 2011), and can lead to frequent relapse following treatment discontinuation (Holyoake and Vetrie, 2017; Mahon et al., 2010; Zhang and Li, 2019). Here, we considered a data set of scRNA-seq data for more than 2000 CML-SCs selected with $C D 3 4^{+} C D 3 8^{-}$ surface phenotype from patients with CML throughout the disease course from diagnosis to 5 years of TKI therapy (Giustacchini et al., 2017). Figure 2A shows the stemness index of cells with either BCR-ABL-positive or -negative at different stages.

FIG. 2.

Stemness of cancer stem cells in CML patients. (A) Stemness of cells at stages of TKI treatment. Dark line refers to BCR-ABL⁺ cells, and light line refers to BCR-ABL⁻ cells. (B) Distribution of the stemness over all cells in the entire stage. CML, chronic myeloid leukemia; TKI, tyrosine kinase inhibitor.

The stemnesses remain at relatively high levels (>0.4), which is in agree with the selection condition of $C D 3 4^{+} C D 3 8^{-}$ surface makers. The stemness of CML-SCs decreases at 12 m after TKI therapy and increases back to the original level at diagnosis. Clinically, patients with prolonged TKI therapy showed two distinct clusters of remission of BCR-ABL⁺ SCs; group A cells were enriched for quiescence and HSC-associated gene expression, and hence were more similar to the normal (BCR-ABL⁻) HSCs, whereas group B cells showed enrichment of MYC, E2F, and proliferation-associated gene sets (Giustacchini et al., 2017).

Prolonged TKI treatment may result in the selective persistence of distinct and highly quiescent BCR-ABL⁺ CML-SCs (group A), and the proportion of group B BCR-ABL⁺ SCs decreases. Hence, the stemness of BCR-ABL⁺ SCs evolved similar to that of BCR-ABL⁻ SCs. Increasing of the stemness from 12 to 18 m after TKI therapy may be associated with the selection of this resistant and quiescent CML-SC subset.

We further measure the heterogeneity of CML-SCs through the distribution of stemnesses of all cells (Fig. 2B). The stemness of all sampled cells varies from 0.4 to 1.0, and showed a unimodal distribution for both BCR-ABL⁺ and BCR-ABL⁻ cells. Hence, the stemness cannot distinguish the two cluster cells as shown in Giustacchini et al. (2017). Previous studies have shown that the heterogeneity of cells can be recovered by a subpopulation of homogeneous cells after about 2 weeks of cell culturing, which indicates the transition of cell states during cell regeneration due to transcriptome-wide noise (Chang et al., 2008; Li et al., 2016). The heterogeneity and the transition of cell stemness during cell cycling are essential for cancer development (Lei, 2020).

3.2. Stemness of mammalian ESCs

We applied the proposed stemness to study the progress of human preimplantation, considering that the data set GSE36552 contains 121 individual cells from human preimplantation embryos (Yan et al., 2013). The cells cover six early developmental stages: zygote (3 cells), 2-cell stage (6 cells), 4-cell stage (12 cells), 8-cell stage (20 cells), morula (16 cells), late blastocyst at hatching stage (30 cells), and the seventh stage of primary outgrowth during human embryonic stem cell (hESC) derivation (34 cells). We calculated the OCHSVM stemness of all cells based on the training set with zygote cells, shown in Figure 3.

FIG. 3.

Stemness of human preimplantation embryo cells at different stages. Dots show the OCHSVM stemness of cells at different stages, and the box plots show the spread and centers for the stemnesses at different stages.

From Figure 3 we can see that the stemnesses of cells at different stages decrease from zygote to hESCs, which is in accordance with the formation of ESCs from zygote during preimplantation. During mammalian early embryo development, the cells undergo epigenetic reprogramming so that the epigenetic markers are erased and rebuilt at the early stage from zygote to the ESCs (Liu et al., 2016; Zhang et al., 2016). Figure 3 shows the evolution of stemness values in cells during embryo development, and reach a relatively low level at the state of hESC, in agree with the epigenetic reprogramming during early embryo development.

3.3. Stemness of iPSC-cardiomyocyte differentiation

We can apply the OCHSVM stemness to study the cardiac maturation process through cardiomyocyte-iPSCs. Here, we study a data set (GSE72612) of differentiated iPSC-cardiomyocytes at three stages: days 0, 21, and 30. The HSCs with CD34 expression higher than 4.5 are used as the training data set (Giustacchini et al., 2017). We calculated the OCHSVM stemness over the tested data set (GSE72612) through the above trained model. The obtained stemness shown in Figure 4 obviously decreases along the differentiation process, which is consistent with the biological intuition of the decreasing of stemness during differentiation. At day 0, the undifferentiated iPSC-cardiomyocytes show high stemness values of about 0.8.

FIG. 4.

OCHSVM stemness of the differentiated iPSC-cardiomyocytes at three stages after iPSC differentiation induction. iPSC, induced pluripotent stem cell.

At day 21, followed by an obvious decrement at day 21, the stemness decreases to 0.68 after the induction of cell differentiation and further decreases to about 0.6, with the minimum value of 0.52, at day 30, with further differentiation. Moreover, the variance of the stemness increases over the differentiation process, which suggests the diversity of iPSC-cardiomyocytes during differentiation, in agreement with the experimental observation.

3.4. Stemness in the dynamic interconversion of pluripotent stem cells

We applied the OCHSVM stemness to study the dynamic interconversion of pluripotent stem cells (PSCs). We obtain the scRNA-seq data from the data set GSE60749 that contains 277 mouse embryonic stem cells (mESCs) and 54 nestin-positive neural precursor cells (NPCs) (Kumar et al., 2014). In our analysis, we took 183 mESCs cultured in serum+LIF media, 94 mESCs cultured in 2i+LIF media (ground-state conditions), and 54 nestin-positive NPCs derived from the mESCs. We calculated the OCHSVM stemness through the trained model on mESCs cultured in 2i+LIF media. The stemness in Figure 5 shows a decrease along the differentiation process from mESCs to NPCs. It should be noted that mESCs cultured in 2i+LIF media establish the ground state, and the mESCs cultured in serum+LIF media exhibit a large degree of heterogeneity, meaning that the cells cultured in serum+LIF media are more prone to be differentiated (Tosolini and Jouneau, 2016).

FIG. 5.

OCHSVM stemness of mESCs and nestin-positive NPCs. mESCs, mouse embryonic stem cells; NPCs, neural precursor cells.

Our computational results are consistent with the experimental analysis, where the stemness values of the cells cultured in 2i+LIF are mainly 1with only a few cells whose stemness values are less than 1, which suggests the decreasing of stemness during differentiation from mESCs to NPCs.

3.5. Stemness in the differentiation of human embryonic stem cells

We examined the stemness of cells along the differentiation of human embryonic stem cells (hESCs) from the data set GSE75748 that contains snapshot progenitors from mesendoderm to definitive endoderm (Chu et al., 2016). The data set contains 374 hESCs, 173 NPCs, 138 definite endoderm progenitors (DEPs), 69 trophoblast cells (TBs), 159 human foreskin fibroblasts (HFFs), and 105 endothelial cells (ECs). First, we trained the model to yield high stemness for hESCs, and then applied the model to calculate the OCHSVM stemness of all cells. Results show the decreasing of stemness along the differentiation process, from the highest level for hESCs to the lowest level for ECs (Fig. 6).

FIG. 6.

Stemness index by OCHSVM with human embryonic stem cell entry into endoderm progenitors.

4. DISCUSSION

We have examined the proposed OCHSVM stemness over four different data sets and have shown that they yield reasonable results in accordance with changes in the stemness of tumor stem cells after TKI therapy, embryo development, and iPSC differentiation. Therefore, the OCHSVM stemness can be a measurement of stemness at a single-cell level based on scRNA-seq data. To show the effectiveness of the proposed OCHSVM method, we compare the results with other state-of-the-art methods. As a supervised, one-class nonlinear method, we first compare the results with representative supervised methods, the OCLR stemness in Sokolov et al. (2016), and one-class SVM methods with different kernel functions: the linear kernel $K (x_{i}, x_{j}) = x_{i}^{T} \cdot x_{j} (O C L S V M)$ , the polynomial kernel $K (x_{i}, x_{j}) = {(x_{i}^{T} \cdot x_{j} + 1)}^{d}$ , and the RBF function $K (x_{i}, x_{j}) = exp (- \frac{∥ x_{i} - x_{j} ∥^{2}}{σ^{2}}) (O C R S V M)$ .

There are two parameters d (in polynomial kernel) and $σ$ (in RBF kernel) that need to be specified, where we used the default parameter $d = 2$ in stemness evaluation for all the considered single-cell data sets. Since $σ$ is very sensitive to the data structure, we use leave-out cross-validations (80% for training, 20% for validations in the training data) among candidate values $1, 1 0^{1}, 1 0^{2}, \dots 1 0^{6}$ to select the robust and optimal $σ$ .

First, we compared the stemness of human preimplantation embryo cells at different stages (GSE36552), shown in Figure 7. We use default parameter $d = 2$ in OCPSVM, and the parameter $σ$ in OCRSVM is set to be $1 0^{5}$ after cross-validations. For all the comparison methods, we calculated the stemness of all cells based on the training set with zygote cells. Biologically, we would expect to see the obvious decrease of stemness from zygote to hESCs, and reach relatively low-level stemness for hESCs, as shown in Figure 3. The results from OCRSVM did not show the obvious tendency of decreasing stemness along embryo development. The OCLR stemness showed similar tendency as the OCHSVM stemness, however, with a lower value of the hESCs. The OCLSVM and OCPSVM stemness cannot capture the tendency, and with high value of stemness for all types of cells.

FIG. 7.

Stemness for different stages of human embryonic stem cell differentiation obtained from different methods (A) OCLR, (B) OCLSVM, (C) OCRSVM, (D) OCPSVM, and (E) OCHSVM. OCLR, one-class logistic regression.

We applied the above methods to the iPSC-cardiomyocyte differentiation data set (GSE72612). The HSCs with CD34 expression higher than 4.5 are used as the training data set (Giustacchini et al., 2017) in model comparisons, the same as OCHSVM model training. We use default parameter $d = 2$ in OCPSVM, and the parameter $σ$ in OCRSVM is set to be 10⁵ after cross-validations. Biologically, the stemness decreases along with the differentiation process. The stemnesses obtained from OCPSVM and OCRSVM, however, showed either extremely high or low levels in all the three stages of cells, and failed to distinguish the stemness at days 0, 20, and 30 (Fig. 8). The OCLR and OCLSVM methods showed the tendency of stemness decreasing from day 0 to 30 after differentiation.

FIG. 8.

Stemnesses for iPSC-cardiomyocyte differentiation at different stages obtained from different methods (A) OCLR, (B) OCLSVM, (C) OCRSVM, (D) OCPSVM, and (E) OCHSVM.

We applied the methods to dynamic interconversion of PSCs from mESCs to NPCs (GSE60749). In the data set, there are 94 mESCs cultured in 2i+LIF media (ground-state conditions). We hence calculated the stemness values through the trained model on mESCs cultured in 2i+LIF media. We use the default parameter $d = 2$ in OPSVM, and the parameter $σ$ is set to be 10³ after cross-validations. The OCLR method, OCLSVM, and OCPSVM showed the decrease of stemness in NPCs compared with mESCs. In the OCPSVM method, the stemness values in mESCs with serum+LIF were shown to be lowly expressed, not consistent with the experimental analysis. OCRSVM, however, failed to distinguish the two types of cells and give high-level stemnesses for both mESCs and NPCs (Fig. 9).

FIG. 9.

Stemnesses for dynamic interconversion of pluripotent stem cells from mESCs to NPCs obtained from different methods (A) OCLR, (B) OCLSVM, (C) OCRSVM, (D) OCPSVM, and (E) OCHSVM.

We applied the methods to human embryonic stem cells undergoing differentiation from hESC to EC (GSE75748). The data set contains 374 hESCs, 173 NPCs, 138 DEPs, 69 TBs, 159 HFFs, and 105 ECs. Similar to OCHSVM, we trained the model to yield high stemness for hESCs, and then applied the model to calculate the stemness of all cells in the comparisons. We use default parameter $d = 2$ in OCPSVM, and the parameter $σ$ is set to be 10⁵ after cross-validations. Results obtained from the four methods, shown in Figure 10, showed different patterns, and none of them revealed the decreasing of stemness along the differentiation process, as seen in Figure 6.

FIG. 10.

Stemness for different stages of human embryonic stem cell differentiation obtained from different methods (A) OCLR, (B) OCLSVM, (C) OCRSVM, and (D) OCPSVM, and (E) OCHSVM.

In summary, we compared the results obtained from OCHSVM and the other four methods; only the OCHSVM stemness can give results that biologically make sense for the four data sets. Hence, the proposed OCHSVM method is valuable in analyzing the dynamic processes related to stem cell differentiation and/or embryo development.

Moreover, we further compared our proposed method with a recent method CytoTRACE for evaluating the differentiation state (Gulati et al., 2020). CytoTRACE is a simple but robust method to evaluate the developmental potential through the number of expressed genes per cell. Different from the above methods, CytoTRACE is an unsupervised learning framework to predict differentiation states from scRNA-seq data.

CytoTRACE is an unsupervised learning framework to reconstruct trajectories of cellular differentiation (Gulati et al., 2020). Figure 11 shows the prediction results of differentiation states for the four scRNA-seq data sets in both CytoTRACE and OCHSVM: GSE36552, GSE72612, GSE60749, and GSE75748. In human preimplantation embryo cells (GSE36552), CytoTRACE in general can capture the tendency of differentiation trajectory from zygote to hESCs (Fig. 11a). Different from the results obtained by the proposed OCHSVM model, the differentiation trajectory underwent slight decrement in the 4-cell stage that may contradict our usual understanding.

FIG. 11.

Cellular differentiation states' comparison in CytoTRACE and OCHSVM. (A, E) Human preimplantation embryo cell (Data set: GSE36552), (A) for CytoTRACE and (E) for OCHSVM. (B, F) iPSC-cardiomyocytes (Data set: GSE72612), (B) for CytoTRACE and (F) for OCHSVM. (C, G) Interconversion of pluripotent stem cells (Data set: GSE60749), (C) for CytoTRACE and (G) for OCHSVM. (D, H) Embryonic stem cell differentiation (Data set: GSE75748), (D) for CytoTRACE and (H) for OCHSVM.

In the data set of iPSC-cardiomyocyte differentiation (GSE72612), both OCHSVM and CytoTRACE can predict the right differentiation trajectory (Fig. 11b, e). The difference between OCHSVM and CytoTRACE lies in that OCHSVM is a supervised method, while CytoTRACE is unsupervised. Moreover, CytoTRACE provides ordering of the considered cells, ranging in [0,1], while OCHSVM can be more elegant to show the degree of differentiation. For example, in GSE36552, the stemness values obtained by the OCHSVM model are in the range [0.8,1], indicating that almost all the cells have great potential of pluripotency.

In GSE60749, mESCs cultured in 2i+LIF media establish the ground state, and the mESCs cultured in serum+LIF media are more prone to be differentiated (Tosolini and Jouneau, 2016). It can be revealed by ordering the results of CytoTRACE (Fig. 11c) that mESCs cultured in 2i-LIF media have larger stemness values.

We applied the CytoTRACE to human embryonic stem cells undergoing differentiation from hESC to EC (GSE75748). Results reveal the increment of differentiation degree along the differentiation process, as seen in Figure 11d and h, where the OCHSVM model can capture the same tendency of cellular differentiation trajectories.

In a further step, we checked the stemness distribution for each data set generated by CytoTRACE and OCHSVM. We estimate probability density f for the sample data in the stemness vector generated by cytoTRACE and OCHSVM, respectively. The estimate is based on a normal kernel function that covers the range of the data, as shown in Figure 12, upper figures refer to OCHSVM probability density, and lower figures refer to CytoTRACE probability density. Besides, we conducted statistical tests to quantitatively compare the performances of both methods, shown in Table 1. We performed t-test on the hypothesis H that the stemness values follow normal distribution with mean 0.5 and unknown variance, and returns the testing result in H. Here $H = 0$ indicates that the null hypothesis (“means are equal”) cannot be rejected at the 5% significance level.

FIG. 12.

Stemness distribution comparison in CytoTRACE and OCHSVM. Upper figures refer to OCHSVM probability density, and lower figures refer to CytoTRACE probability density in human preimplantation embryo cell (Data set: GSE36552), iPSC-cardiomyocyte (Data set: GSE72612), interconversion of pluripotent stem cells (Data set: GSE60749), and embryonic stem cell differentiation (Data set: GSE75748).

Table 1.

Performance Comparisons of CytoTRACE and OCHSVM in Stemness Distribution

Algorithm/data set	GSE36552	GSE72612	GSE60749	GSE75748
CytoTRACE	$H = 0, p = 0.6488$	$H = 0, p = 1$	$H = 0, p = 0.2173$	$H = 0, p = 1$
OCHSVM	$H = 1, p = 2.97 e^{- 104}$	$H = 1, p = 1.57 e^{- 35}$	$H = 1, p = 1.19 e^{- 120}$	$H = 1, p = 0$

Null hypothesis: The stemness vector is from normal distribution with mean 0.5 and unknown variance.

It is interesting to see that, for all the data sets, CytoTRACE always estimates the stemness vector following a normal distribution with mean 0.5, also indicated by Figure 12, whose probability density values achieve optimal at 0.5. However, it is not quite reasonable that the stemness values in all the data sets have the same expectation in distribution. In comparison, the stemness values estimated by the OCHSVM model are diversely distributed over different data sets. Although the general tendency of stemness in both methods is the same, from the statistical tests we can conclude that the OCHSVM model is more reasonable.

4.1. Applications

Intratumoral heterogeneity and the presence of CSCs are challenging issues in cancer therapy. In this study, we apply the stemness identification method to check the heterogeneity in CSCs. We extracted the data set from the NCBI GEO database GSE81730 for scRNA-seq data of HSCs from three CML patients. There are 240 HSCs, which were separated into CML-SCs ( $B C R - {A B L}^{+}$ ) or normal HSCs ( $B C R - {A B L}^{-}$ ) based on the presence of the BCR-ABL transcript. The number of CML-SCs in three patient samples was 50, 25, and 51, respectively, and the number of normal HSCs was 46, 23, and 45 correspondingly. Specifically, we trained the model with CML-SCs in patient 3, and evaluated the stemnesses of CML-SCs in other two patients. From the stemness distribution, we discovered heterogeneity in different patients (Fig. 13). Although all the cells were CML stem cells, they behave differently in different patients.

FIG. 13.

Stemness distribution for chronic myeloid leukemia stem cells in three patient samples.

In patient 3, the stem cells were similar to each other in terms of stemness, which showed compact distribution with a small variance. While the stem cells in patient 2 showed more differentiated features and are more scattered in the stemness distributions. These results showed that the proposed stemness can be an indication of cancer heterogeneity.

5. CONCLUSION

In this article, we proposed a nonlinear one-class SVM method, OCHSVM, to calculate the stemness of cells using scRNA-seq data. The OCHSVM method considered the one-class SVM framework for model training, and used the Hadamard kernel to evaluate the nonlinear correlation between different cells. We tested the trained model to analyze CSC development in CML patients after TKI therapy, and explored the evolution of cancer cell stemness in response to the therapy. We also applied the OCHSVM stemness to study the process of human ESC formation during preimplantation, which showed an obvious decrease in the stemness from zygote to hESCs, in agreement with the process of epigenetic reprogramming in the mammalian early embryo development. Moreover, OCHSVM stemness can characterize the differentiation process of iPSC-cardiomyocytes, and the early development of human stem cell differentiations.

The proposed method was compared with a number of state-of-the-art machine learning methods for stemness evaluation, including OCLR and one-class SVM method with several different functions: OCLSVM, OCPSVM, and OCRSVM. Results show that OCHSVM method provides an effective way for modeling the cellular relationships. It was further compared with the most-recent method (an unsupervised learning framework) for differentiation trajectory evaluation, CytoTrace, and can still show superiority and robustness. Besides, the method was also applicable for evaluating tumor heterogeneity. In summary, stemness is an important quantification to measure the ability of self-renewal of stem cells; the proposed OCHSVM stemness offers a method to quantitatively analyze the differentiation status of cells, and to characterize tumor heterogeneity and the patterns in different stages of treatment.

Footnotes

ACKNOWLEDGMENTS

The authors would like to thank the anonymous reviewers for their valuable suggestions.

AUTHOR DISCLOSURE STATEMENT

The authors declare they have no conflicting financial interests.

FUNDING INFORMATION

This research is supported by the National Natural Science Foundation of China Grant Nos. 12271522, 11901575, and 91730301.

References

Alqarni

, Rambely

, Alharbi

, et al. Dynamic behavior and stabilization of brain cell reconstitution after stroke under the proliferation and differentiation processes for stem cells. Math Biosci Eng, 2021; 18(5):6288.

Bao

. An elementary mathematical modeling of drug resistance in cancer. Math Biosci Eng, 2021; 18(1):339.

Batlle

, Clevers

. Cancer stem cells revisited. Nat Med, 2017; 23(10):1124–1134.

Bendall

, Davis

, David Amir

, et al. Single-cell trajectory detection uncovers progression and regulatory coordination in human B cell development. Cell, 2014; 157(3):714–725.

Cao

, Bhatnagar

, Wang

, et al. Cancer stem cells and strategies for targeted drug delivery. Drug Deliv Transl Res, 2021; 11:1779–1805.

Chang

, Hemberg

, Barahona

, et al. Transcriptome-wide noise controls lineage choice in mammalian progenitor cells. Nature, 2008; 453(7194):544–547.

Chu

, Leng

, Zhang

, et al. Single-cell RNA-seq reveals novel regulators of human embryonic stem cell differentiation to definitive endoderm. Genome Biol, 2016; 17(1):173.

Chu

, McDonald

, Lin

, et al. Persistence of leukemia stem cells in chronic myelogenous leukemia patients in prolonged remission with imatinib treatment. Blood, 2011; 118:5565–5572.

Damavandi

, Yi

, Buczacki

. Cancer Stem Cells. In: Encyclopedia of Cell Biology (2nd edition). ( Bradshaw

, Hart

, Stahl

. eds.) Academic Press: Oxford, 2023; pp. 114–123.

10.

Ellsworth

, Blackburn

, Shriver

, et al. Single-cell sequencing and tumorigenesis: Improved understanding of tumor evolution and metastasis. Clin Transl Med, 2017; 6:15.

11.

Farrell

, Yiqun

, Riesenfeld

, et al. Single-cell reconstruction of developmental trajectories during zebrafish embryogenesis. Science, 2018; 360:eaar3131.

12.

Friedmann-Morvinski

, Verma

. Dedifferentiation and reprogramming: Origins of cancer stem cells. EMBO Rep, 2014; 15:244–253.

13.

Giustacchini

, Thongjuea

, Barkas

, et al. Single-cell transcriptomics uncovers distinct molecular signatures of stem cells in chronic myeloid leukemia. Nat Med, 2017; 23(6):692–702.

14.

Grün

, Lyubimova

, Kester

, et al. Single-cell messenger RNA sequencing reveals rare intestinal cell types. Nature, 2015; 525:251–255.

15.

Gulati

, Sikandar

, Wesche

, et al. Single-cell transcriptional diversity is a hallmark of developmental potential. Science, 2020; 367(6476):405.

16.

Gupta

, Fillmore

, Jiang

, et al. Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells. Cell, 2010; 146(4):633–644.

17.

Holyoake

, Vetrie

. The chronic myeloid leukemia stem cell: Stemming the tide of persistence. Blood, 2017; 129(12):1595–1606.

18.

Jiang

, Tu

, Zhou

, et al. Single-cell transcriptomics reveal the heterogeneity and dynamic of cancer stem-like cells during breast tumor progression. Cell Death Dis, 2021; 979(12).

19.

Jiang

, Ching

W-K

, Cheung

W-S

, et al. Hadamard Kernel SVM with applications for breast cancer outcome predictions. BMC Syst Biol, 2017; 11:138.

20.

Jiang

, Ching

W-K

, Yiu

CKF

, et al. Stationary Mahalanobis kernel SVM for credit risk evaluation. Appl Soft Comput, 2018; 71:407–417.

21.

Jiang

, Chen

, Pinello

, et al. GiniClust: Detecting rare cell types from single-cell gene expression data with Gini index. Genome Biol, 2016; 17:144.

22.

Kumar

, Cahan

, Shalek

, et al. Deconstructing transcriptional heterogeneity in pluripotent stem cells. Nature, 2014; 516(7529):56–61.

23.

Lai

, Yu

, Huang

, et al. A Novel Support Vector Machine Metamodel for Business Risk Identification. ( Yang

, Webb

. eds.) PRICAI 2006: Trends in Artificial Intelligence, Volume 4099 of Lectures Notes in Computer Science. Springer: Berlin, Heidelberg; 2006.

24.

, Shen

, Cory

A-S

. Lineage plasticity in cancer progression and treatment. Annu Rev Cancer Biol, 2018; 2:271–289.

25.

Lei

A general mathematical framework for understanding the behavior of heterogeneous stem cell regeneration. J Theor Biol, 2020; 492(110196).

26.

, van der Leun

, Yofe

, et al. Dysfunctional CD8 T cells form a proliferative, dynamically regulated compartment within human melanoma. Cell, 2019; 176(4):775–789.

27.

, Wennborg

, Aurell

, et al. Dynamics inside the cancer cell attractor reveal cell heterogeneity, limits of stability, and escape. Proc Natl Acad Sci USA, 2016; 113(10):2672–2677.

28.

Liu

, Wang

, Liu

, et al. Distinct features of H3K4me3 and H3K27me3 chromatin domains in pre-implantation embryos. Nature, 2016; 537(7621):558–562.

29.

Macaulay

, Svensson

, Labalette

, et al. Single-cell RNA-sequencing reveals a continuous spectrum of differentiation in hematopoietic cells. Cell Rep, 2016; 14(4):1–13.

30.

Mahon

F-X

, Réa

, Guihot

, et al. Discontinuation of imatinib in patients with chronic myeloid leukaemia who have maintained complete molecular remission for at least 2 years: The prospective, multicentre Stop Imatinib (STIM) trial. Lancet Oncol, 2010; 11(11):1029–1035.

31.

Malta

, Sokolov

, Gentles

, et al. Machine learning identifies stemness features associated with oncogenic dedifferentiation. Cell, 2018; 173(2):338–354.e15.

32.

Miranda

, Hamilton

, Zhang

, et al. Cancer stemness, intratumoral heterogeneity, and immune response across cancers. Proc Natl Acad Sci, 2019; 116(18):9020–9029.

33.

Navin

, Kendall

, Troge

, et al. Tumour evolution inferred by single-cell sequencing. Nature, 2011; 472:90–94.

34.

Nawy

. Single-cell sequencing. Nat Methods, 2014; 11(1):18.

35.

Paolo

, Abraham

, Holyoaka

. Hurdles toward a cure for CML: The CML stem cell. Hematol Oncol Clin North Am, 2011; 25(5):951–966.

36.

Pasic

, Lipton

. Current approach to the treatment of chronic myeloid leukaemia. Leuk Res, 2017; 55:65–78.

37.

Qin

, Furutani

, Yonezawa

, et al. Targeting transglutaminase 2 mediated exostosin glycosyltransferase 1 signaling in liver cancer stem cells with acyclic retinoid. Cell Death Dis, 2023; 14(6):385.

38.

Ren

, Kang

, Zhang

. Understanding tumor ecosystems by single-cell sequencing: Promises and limitations. Genome Biol, 2018; 19(1):211.

39.

Schölkopf

, Williamson

, Smola

, et al. Support Vector Method for Novelty Detection. In: Advances in Neural Information Processing Systems. ( Solla

, Leen

, Muller

K-R

, eds.) MIT Press: Denver, Co, USA, 2000; pp. 582–588.

40.

Shibue

, Weinberg

. EMT, CSCs, and drug resistance: The mechanistic link and clinical implications. Nat Rev Clin Oncol, 2017; 14:611–629.

41.

Sokolov

, Paull

, Stuart

. One-class detection of cell states in tumor subtypes. Pac Symp Biocomput, 2016; 21:405–416.

42.

, Wei

, Robert

, et al. Single-cell analysis resolves the cell state transition and signaling dynamics associated with melanoma drug-induced resistance. Proc Natl Acad Sci USA, 2017; 114(52):13679–13684.

43.

Tosolini

, Jouneau

. Acquiring ground state pluripotency: Switching mouse embryonic stem cells from serum/lif medium to 2i/lif medium. Methods Mol Biol, 2016; 1341:41–48.

44.

, Zhang

, et al. Single-cell data-driven mathematical model reveals possible molecular mechanisms of embryonic stem-cell differentiation. Math Biosci Eng, 2019; 16(5):5877–5896.

45.

Walcher

, Kistenmacher

, Suo

, et al. Cancer stem cellsorigins and biomarkers: Perspectives for targeted personalized therapies. Front Immunol, 2020; 11:1280.

46.

Wang

, Xing

, Schaid

. Kernel methods for large-scale genomic data analysis. Brief Bioinform, 2014; 16(2):183–192.

47.

Yan

, Yang

, Guo

, et al. Single-cell RNA-Seq profiling of human preimplantation 115 embryos and embryonic stem cells. Nat Struct Mol Biol, 2013; 20:1131–1139.

48.

Yejing

, Gomez

, Adam

, et al. Stem cell lineage infidelity drives wound repair and cancer. Cell, 2017; 169:636–650.

49.

Zhang

, Zheng

, Huang

, et al. Allelic reprogramming of the histone modification H3K4me3 in early mammalian development. Nature, 2016; 537(7621):553–557.

50.

Zhang

, Li

. Exploiting unique biological features of leukemia stem cells for therapeutic benefit. Stem Cells Transl Med, 2019; 8(8):768–774.

51.

Zheng

, Zheng

, Yoo

J-K

, et al. Landscape of infiltrating T cells in liver cancer revealed by single-cell sequencing. Cell, 2017; 169(7):1342–1356.e16.

52.

Zheng

, Song

, Fu

, et al. An absolute human stemness index associated with oncogenic dedifferentiation. Brief Bioinform, 2021a;22(2):2151–2160.

53.

Zheng

, Xie

, Song

, et al. Stemsc: A cross-dataset human stemness index for single-cell samples. Stem Cell Res Ther, 2021b;13(1):115.