A clustering method for small scRNA-seq data based on subspace and weighted distance

Datasets

The improved k-means algorithm with EP_dis and relative CH (RCH)

The k-means is a widely used clustering algorithm (MacQueen, 1967; Jain, Murt & Flynn, 1999; Jain, 2008). The algorithm requires the user to provide cluster initialization, distance metric, and the cluster numbers as the parameters (Chiang & Mirkin, 2010). Here we designed an improved k-means algorithm by introducing the EP_dis and RCH, which measure the similarity between two cells more appropriately and can automatically determine the cluster numbers.

Determine the number of clusters

The k-means algorithm needs to be specified the number of clusters. The clustering internal validation (CIV) indices, such as the Calinski-Harabasz (CH) index (Caliński & Harabasz, 1974), Silhouette (Sil) index (Rousseeuw, 1987), and Gap Statistic (Tibshirani & Hastie, 2001), can be used for estimating the cluster numbers. The CH has been proven the best in estimating cluster numbers (Milligan & Cooper, 1985; Chiang & Mirkin, 2010). It is defined as: (2)${\displaystyle CH=\frac{\frac{S{S}_B}{k-1}}{\frac{S{S}_W}{N-k}},k=2,3\dots NC,}$ where N as the sample numbers, NC as the largest cluster numbers. The k with the maximum CH is the suitable cluster numbers. In different k, the CH value is incomparable because the degrees of freedom vary. So, we designed a new index, relative CH (RCH) (Ning et al., 2022), that was relatively comparable under different k: (3)${\displaystyle RC{H}_k=\frac{C{H}_k}{F_{\left(\alpha ,k-1,N-k\right)}}.}$ The workflow of the improved k-means algorithm with EP_dis and RCH is shown in Fig. 1.

The overview of the SSWD

In the scRNA-seq data matrix X_G×N = {x_ij|1 ≤ i ≤ G, 1 ≤ j ≤ N}, rows represent genes, and columns represent cells. x_ij represents the value of gene i in the j th cell. The framework of SSWD is depicted in Fig. 2.

Time complexity of SSWD

The main time-consuming step of SSWD is clustering by the improved k-means with EP_dis and RCH. In step 2 (see The overview of the SSWD), we used the improved k-means algorithm in the density matrix. We denoted n represents the sample numbers, m represents the feature numbers, k represents the cluster numbers, NC represents the range of cluster numbers, l represents the iteration numbers to determine the cluster centers, and w_step represents the search step. Since the k<<n, NC<<n, the step 2 time complexity holds about O(lmn). In step 3, each subspace would be performed the improved k-means algorithm after PCA. We denoted d as the retained dimension after PCA, and s is the number of genes subspace. The SSWD time complexity has roughly O(lmn +lnds). Since d<<m, s<<m, we can simplify the time complexity of SSWD to approximately O(lmn).

Biological insights

We transformed the clustering results of each cell group into “one-against-the-rest”. Then, we executed the Wilcoxon rank-sum test for each gene between the expression value and the binary cluster, adjusting the p-value based on FDR. The gene that adjusted p-value<0.001 was preserved as the differential gene. Next, we used the AUC score to evaluate the performance of genes in distinguishing different cell types. Since AUC was only suitable for dichotomous problems, we constructed a binary classifier based on the mean expression value of each gene and compared the processed values with the binary cluster value. We defined the genes with AUC >0.85 and p-value<0.001 as marker genes.

Evaluation metrics

Two external validation indices, ARI (Adjusted Rand Index) and NMI (Normalized Mutual Information) were used to evaluate the effectiveness of clustering methods.

ARI (Hubert & Arabie, 1985) is a widely used external validation index in clustering, and it is defined as follows: (4)${\displaystyle ARI\left(R,C\right)=\frac{\sum _{ij}\left(\begin{array}{c}\hfill n_{ij}\hfill \\ \hfill 2\hfill \end{array}\right)-\left[\sum _i\left(\begin{array}{c}\hfill a_i\hfill \\ \hfill 2\hfill \end{array}\right)\sum _j\left(\begin{array}{c}\hfill b_j\hfill \\ \hfill 2\hfill \end{array}\right)\right]/\left(\begin{array}{c}\hfill n\hfill \\ \hfill 2\hfill \end{array}\right)}{\frac{1}{2}\left[\sum _i\left(\begin{array}{c}\hfill a_i\hfill \\ \hfill 2\hfill \end{array}\right)+\sum _j\left(\begin{array}{c}\hfill b_j\hfill \\ \hfill 2\hfill \end{array}\right)\right]-\left[\sum _i\left(\begin{array}{c}\hfill a_i\hfill \\ \hfill 2\hfill \end{array}\right)\sum _j\left(\begin{array}{c}\hfill b_j\hfill \\ \hfill 2\hfill \end{array}\right)\right]/\left(\begin{array}{c}\hfill n\hfill \\ \hfill 2\hfill \end{array}\right)},}$ where R and C are published and predicted clusters, respectively. The overlap of samples between R and C can be generalized into a contingency table. n_ij is the times a sample occurs in the ith cluster of R and the jth cluster of C, a_i is the sum of the ith row in the contingency table, b_j is the sum of the jth column in the contingency table, and (.) represents the binomial coefficient.

NMI (Strehl & Ghosh, 2002) is defined as follows: (5)${\displaystyle NMI\left(R,C\right)=\frac{2\ast I\left(R,C\right)}{\left[H\left(R\right)+H\left(C\right)\right]},}$ where $I\left(R,C\right)={\sum }_{i=1}^{\left|R\right|}{\sum }_{j=1}^{\left|C\right|}{p}_{ij}log\left(\frac{p_{ij}}{p_i{p}_j}\right)$ is the mutual information between R and C, $H\left(R\right)=-{\sum }_{i=1}^{\left|R\right|}{p}_{i}log{p}_i$ is the entropy with R, $H\left(C\right)=-{\sum }_{i=1}^{\left|C\right|}{p}_{i}log{p}_i$ is the entropy with C, $p_{ij}=\frac{n_{ij}}{n}$ is the probability that a cell belongs to both the ith cluster in R and the jth cluster in C. The range of ARI and NMI are [0, 1]. The larger ARI (NMI) represent a better performance of clustering.

Reference methods

In this article, seven prevailing clustering algorithms were introduced as reference methods. The SC3 v.1.22.0 (Kiselev et al., 2017), CIDR v.0.1.5 (Lin, Troup & Ho, 2017), Seurat v.4.1.1 (Satija et al., 2015) , SIMLR v.1.20.0 (Wang et al., 2017) were implemented with the original R package in Rstudio4.0. SinNLRR (Zheng et al., 2019) (https://github.com/zrq0123/SinNLRR) and S3C2 (Zhuang et al., 2021) (https://github.com/Cuily-v/S3C2) were implemented in Matlab2017a. The SNN-Cliq (Xu & Su, 2015) (https://github.com/BIOINSu/SNN-Cliq) was run in Matlab2017a and Python3.8. SILMR and SNN-Cliq used the same log transformation as this paper. SinNLRR and S3C2 used the correct cell groups for clustering. Unless specified, the default parameters in the program were used as suggested in the original paper.

Performance evaluation and comparison with reference methods

We compared the performance of SSWD with seven prevailing clustering methods in eight scRNA-seq datasets (Table 4). The SSWD achieved the best clustering performance with an average ARI of 0.791 and was 0.143 higher than the second-ranked SC3, whereas the SNN-Cliq had poor performance (ARI of 0.364). SSWD ranked in the top three for ARI on all other datasets except Yan. SSWD attained the best results for NMI in three datasets (Li, Tian305, Tian307) and the second-best in four datasets (Biase, Yan, Deng, Treutlein). The average NMI of SSWD was the highest (0.850). Seurat had the poorest performance with only 0.579 in NMI because it failed on Biase, and the NMI of Li was only 0.122.

We further demonstrated the SSWD performance by ranking clustering accuracies on eight datasets (Fig. 3). For ARI (Fig. 3A), SSWD was superior to the seven reference methods in rank-wise (median of 2). SNN-Cliq performed the worst, with a median of 7. For NMI (Fig. 3B), SSWD was also superior to others, and the performance of CIDR, Seurat, and SNN-Cliq was all poor. Furthermore, the one-sided Wilcoxon signed-rank test was used to explain the statistical difference between SSWD and the reference methods. Except for SC3 and SinNRLL (in ARI), all the p-value are less than 0.05, which shows SSWD is superior to other methods (Table 5).

The SSWD was also better for estimating the cell groups. Five out of eight datasets (Biase, Patel, Li, Tian307, Tian305) acquired the correct cell groups using SSWD. Especially for Tian307 and Tian305, only the SSWD estimated the correct cell groups and achieved the best ARI of over 0.948. Deng contains seven cell groups, for which all methods failed to identify the correct number of cell groups. For Treutlein, CIDR and Seurat estimated the correct cell groups, but the ARI (0.188 and 0.531) and NMI (0.304 and 0.648) were lower than those of SSWD (0.607 and 0.732). For Yan, the SinNRLL and S3C2 performed very well under the correct cell groups.

Annotate the clusters

We illustrated the effectiveness of the cell annotation using PanglaoDB (Oscar, Li-Ming & Johan, 2019) to clusters taking the Li dataset as an example. Li is a human pancreatic islet cells dataset containing five subtypes (alpha, beta, pp, acinar, and ductal) (Li et al., 2016). According to the AUC score (see Biological insights), we obtained the marker genes for each cluster identified by SSWD. Figure 4 is the expression heatmap of the top 10 marker genes for each cluster, which was divided into five clear modules and indicated that these marker genes could distinguish the clusters well. The keration8 (KRT8) in cluster 1; transthyretin (TTR), glucagon (GEG) in cluster 2; insulin (INS) in cluster 3; pancreatic polypeptide (PPY) in cluster 4; REG1B, REG1A, CTRB2 in cluster 5 were all reported in the original publication. We also annotated the cluster with PanglaoDB. The cluster results annotated with PanglaoDB are consistent with the cell annotations in the original publication (Table 6).

Role of the EP_dis metric

EP_dis was used in SSWD to assess the similarity between cells or genes. According to the EP_dis definition (see Materials & Methods), the optimal w was determined by SS_B/SS_W using a search strategy. When w = 1, the EP_dis equals the Euclidean distance; when w = 0, it is the Pearson distance. We used two simulated datasets, D1 and D2, to display the impact of EP_dis and explain the process of optimizing w by SS_B/SS_W. Figure 5 shows the clustering accuracy of D1 (Fig. 5A) and D2 (Fig. 5B) under different w. It can be seen that the highest scores (CA, Rand, and SS_B/SS_W) in D1 and D2 are not appearing at the endpoints (0.6 in D1 and 0.8 in D2), which indicates that the EP_dis could capture more information between samples.

Role of the relative CH

The RCH in the improved k-means algorithm was used to determine the cluster numbers. In SSWD, we employed RCH to estimate the gene subspace numbers and guide each set of gene subspace grouping. The capability of the RCH directly affects the performance of the SSWD. We utilized simulated datasets D3–D6 with different characteristics, six UCI datasets, and three scRNA-seq datasets to illustrate RCH properties and compare them with CH (Table 7). We can see that CH and RCH were consistent in D3–D6, indicating their good performance in simulated datasets. In the UCI datasets, RCH could estimate the correct cluster numbers except for Dermatology and Yeast, but the corresponding cluster numbers estimated by the RCH was closer to the real value than those of CH. For the scRNA-seq datasets, RCH and CH all failed. Their poor performance may be due to the characteristics of scRNA-seq data. Nonetheless, the RCH result was closer to the true value.

Role of the subspace

After performing steps 1 and 2 of SSWD (see “Materials & Methods”), the Li has been separated into eight sets of gene subspace, and seven participate in consensus clustering (Fig. 6, Fig. S2). The expression heatmaps of the best three sets of genes subspace display clear patterns (Figs. 6A–6C), and their EP_dis heatmap (Figs. 6D–6F) effectively clustered cells with similar expression patterns. Compared with the EP_dis heatmap by 1,000 genes (Fig. 6G), the consensus matrix using sets of genes subspace (Fig. 6H) enhances intercellular signaling. The consensus matrix clustering result was better (ARI of 0.967, NMI of 0.964) than the former (ARI of 0.386, NMI of 0.579) because the former could not distinguish alpha and pp cells well.

Discussion of prevailing methods

We also provided further discussion in Tables 4–5. The average performance of SC3 was only lower than SSWD, and its results were not significantly different from SSWD in the one-sided Wilcoxon signed-rank test. The SC3 combined multiple similarity measures (Euclidean, Pearson, Spearman) in clustering. It used the consistency matrix to integrate the multiple clustering results, and the consistency matrix strengthened the consensus signal between cells. At the same time, we can see that Deng and Treutlein, the best accuracy performers in SC3, could not obtain the correct number of cell groups. Although Biase estimated the correct number of cell groups in SC3, one cell was classified mistakenly, while SinNRLL could classify all cells accurately. Both SinNRLL and S3C2 introduce the idea of subspace clustering. Their average performances were better than other methods except for SSWD and SC3. However, this result was based on the cell group numbers being provided. Evaluating the cluster numbers is an important aspect of clustering methods. Although SinNRLL could estimate the cluster numbers by other methods, its accuracy is still unsatisfactory (Zheng et al., 2019).

SNN-Cliq performed the worst (ARI = 0.364, NMI = 0.624), with none of the seven datasets estimating the correct number of cell groups. SNN-Cliq tended to divide more clusters, probably because the method requires providing three suitable parameters, and the results depend on the graphical representation of the data. CIDR used an implicit imputation approach to reduce the impact of dropout in scRNA-seq and used CH to estimate the cell groups. The method determined the correct number of cell groups in Treutlein and Tian307, but their clustering accuracies were poor. SIMLR adopts a multi-kernel strategy to adaptively select an appropriate distance metric and automatically determine the cell groups. However, this method achieved good performance only in Patel because it used Euclidean distance as the metric to construct a Gaussian kernel function (Taiyun et al., 2018). For Seurat, Biase failed, and the ARI of Li was only 0.084. The results show that Seurat may be unsuitable for small datasets, consistent with the literature (Kiselev, Andrew & Hemberg, 2019).

The SSWD had the best performance in experiments. However, the performance of Patel and Yan were mediocre. Although Patel estimated the correct cell groups, the clustering accuracies were only ranked the third (in ARI) and the fourth (in NMI), probably because there were negative values in Patel datasets. All methods failed to estimate the correct cell group numbers in Yan. The poor performance of Yan in SSWD was because the estimated cell groups was far from the actual numbers.

We can draw the following conclusions from the above observations: (1) Due to the complex structure of scRNA-seq data, developing an optimal clustering method for all situations is impossible. (2) Determining the cluster numbers is difficult, so assigning cells to appropriate types is more important. (3) Selecting suitable similarity measures and using subspace in single-cell clustering help obtain better clustering results.