Single-cell Curriculum Learning-based Deep Graph Embedding Clustering

Given that the majority of the structure and information inherent in the scRNA-seq data $\mathbf{X}$ is conserved within the latent embedded representation generated by the scCLG encoder. The adjacency matrix decoder of the graph autoencoder can be defined as the inner product between the latent embedding:

where, $f_{E}$ represents the scCLG encoder function; $\mathbf{A}_{rec}$ is the reconstructed adjacency matrix. Therefore, the reconstruction loss of $\mathbf{A}$ and $\mathbf{A}_{rec}$ should be minimized in the learning process as below:

In order to more effectively capture the structure of scRNA-seq data by decoding from the latent embedded representation $\mathbf{Z}$, we integrate the ZINB model into a Chebyshev graph convolutional autoencoder to capture the global probability structure of the scRNA-seq data. Based on this foundation, we propose to apply the ZINB distribution model to simulate the data distribution to capture the characters of scRNA-seq data as follows:

where $\mu$ and $\theta$ are the mean and dispersion in the negative binomial distribution, respectively. $\pi$ is the weight of the point mass at zero. The proportion $\frac{\theta}{\theta+\mu}$ replaces the probability p. After that, to model the ZINB distribution, the decoder network has three output layers to compute the three sets of parameters. The estimated parameters can be defined as follows:

where $f_{D_{zinb}}$ is a three-layer fully connected neural network with hidden layers of 128, 256 and 512 nodes. $\mathbf{W}$ represents the learned weights parameter matrices. $\hat{\pi},\hat{\mu}$ and $\hat{\theta}$ are parameters denoting the estimations of $\pi,\mu$ and $\theta$, respectively. The selection of the activation function depends on the range and definition of the parameters. In terms of the parameter $\pi$, the suitable activation function for it is sigmoid because the interval of $\pi$ is between 0 and 1. Due to the non-negative value of the mean $\mu$ and dispersion $\theta$, we choose the exponential activation function for them. The negative log-likelihood of the ZINB distribution can be used as the reconstruction loss function of the original data $\mathbf{X}$, which can be defined as below:

scRNA-seq clustering clustering as an unsupervised learning task, lacks guidance from labels, which makes it difficult to capture effective optimization signals during the training phase. To overcome this challenge, we apply a clustering module to guide the algorithm to adjust the cluster centers to ensure that the distribution of samples within each cluster is as consistent as possible while minimizing inter-cluster differences. The objective takes the form of Kullback–Leibler (KL) divergence and is formulated as follows:

where $q_{iu}$ is the soft label of the embedding node $z_{i}$ which is defined as the similarity between $z_{i}$ and cluster centre $\mu_{u}$ measured by Student’s t-distribution. This can be described as follows:

Meanwhile, $p_{iu}$ is the auxiliary target distribution, which puts more emphasis on the similar data points assigned with high confidence on the basis of $q_{iu}$, as below:

Since the target distribution $P$ is defined based on $Q$, the embedding learning of $Q$ is supervised in a self-optimizing way to enable it to be close to the target distribution $P$.

Our Hierarchical Difficulty Measurer consists of two difficulty measures from different perspectives. In this section, we present the definition of two difficulty measures and how to calculate them.

Local Difficulty Measurer. We introduce how to identify difficult nodes from a local perspective. Nodes located at the boundaries of multiple classes may reside in transitional regions within the feature space, leading to less distinct or consistent feature representations, thereby increasing the difficulty of classification. The first type of difficult node should have diverse neighbors that belong to multiple classes. Intuitively, features of nodes within the same class tend to be more similar. This is due to the influence of neighboring node features, resulting in nodes with similar connectivity patterns exhibiting comparable feature representations. In order to identify these difficult nodes, we calculate the diversity of the neighborhood’s features:

where $S(u,v)$ denotes the similarity between cell $u$ and cell $v$. A larger $D_{local}(u)$ indicates a more diverse neighborhood. $\mathcal{N}(u)=\{v\in\mathcal{V}|(u,v)\in\mathcal{E}\}$ is the neighborhood of cell $u$. As a result, during neighborhood aggregation, these nodes aggregate neighbors’ features to get an unclear representation, making them difficult for GNNs to learn. By paying less attention to these difficult nodes, scCLG learns more useful information and effectively improves the accuracy of backbone GNNs.

Global Difficulty Measurer. Then we introduce how to identify difficult nodes from a global perspective. Entropy plays a pivotal role in feature selection as a metric from information theory used to quantify uncertainty. In the process of feature selection, we leverage entropy to assess a feature’s contribution to the target variable. When a feature better distinguishes between different categories of the target variable, its entropy value tends to be relatively low, signifying that it provides more information and reduces overall uncertainty. Consequently, in feature selection, lower entropy values indicate features that offer greater discriminatory power, aiding in the differentiation of target variable categories. We assume nodes that have lower entropy have fewer contributions to the graph. Therefore, this type of node is difficult to classify. Inspired by Entropy Variation [23], We consider the node contribution as the variation of network entropy before and after its removal.

For a node $v_{i}$ in graph $\mathcal{G}$, we define $p(v)$ as probabilities:

where $\sum_{v\in\mathcal{V}}p(v)=1$.

The entropy of the graph is as follows:

where $D(v)$ is the degree of node $v$. $Ent(\mathcal{G})$ is the entropy of graph $\mathcal{G}$ with degree matrix.

The global difficulty of the node is as follows:

where $Ent(v)$ is the change if one node and its connections are removed from the network. $\hat{\mathcal{G}}_{v}$ is the modified graph under the removel of $v$. A lower $Ent(v)$ indicates a lower influence on graph structure and is also more difficult. The global difficulty of node $v$ is to subtract the normalized $Ent(v)$ from 1.

Considering two difficulty measurers from local and global perspectives, we finally define the difficulty of $v$ as:

where $\beta$ is the weight coefficient assigned to each difficulty measurer to control the balance of the total difficulty.

With the hierarchical difficulty measurer, we can get a list of nodes sorted in ascending order of nodes based on difficulty. The node at the end of the list is a nuisance for the overall model learning, so should it be retained? The sources of noise in graph neural networks can be varied, firstly, the attribute information of the nodes may contain noise, which affects the representation of the node features and hence the learning of the GNN. Secondly, the presence of anomalous data may cause the spectral energy of the graph to be ”right-shifted”, the distribution of spectral energy shifts from low to high frequencies. These noises will not only reduce the performance of the graph neural network but also propagate through the GNN in the topology, affecting the prediction results of the whole network. In order to solve this problem, we designed a data pruning strategy based on the calculated node difficulty. Specifically, we define a data discarding hyperparameter $\alpha$. The value of $\alpha$ is set while balancing data integrity and model generalization performance. As shown in Fig. 4, the scRNA-seq clustering performance of the scCLG improves after removing the node features with the highest difficulty which prove our hypothesis.

$k$ represents the number of nearest neighbors to consider when constructing cell graph. In order to investigate the impact of $k$, we ran our model with the parameters 5, 10, 15, 25. Fig. 3 (A) shows the NMI and ARI values with different numbers of $k$. As depicted in Fig. 3 (A), we observe that the two metrics first increase rapidly from parameter 5 to 10, reach the best value at $k=20$, and then decrease slowly from parameter 20 to 25. Therefore, we set the neighbor parameter k as 20 in our scCLG model.

In single-cell data analysis, highly variable genes vary significantly among different cells, which helps to reveal the heterogeneity within the cell population and more accurately identify cell subpopulations. To explore the impact of the number of selected highly variable genes, we apply scCLG on real datasets with gene numbers from 300 to 1500. Fig. 3 (B) shows the line plot of the average NMI and ARI on the 7 datasets selecting 300, 500, 1000 and 1500 genes with high variability, respectively. It can be seen that the performance with 500 highly variable genes is better, while the performance with 300 genes is much worse than the others. Therefore, to save computational resources and reduce running time, we set the number of selected high-variance genes in the model to 500.

In single-cell data analysis, data quality can be improved by pruning lower-quality samples thereby affecting the ability to generalize the model. To explore the impact of the selected data, we run our model with pruning rate parameters from 0.06 to 0.21 to drop difficult nodes. We also compared our pruning strategy with two different pruning strategies, namely pruning easy nodes and randomly pruning nodes. Fig. 4 shows the ARI and NMI values with different numbers of $\alpha$ and pruning strategy. As depicted in Fig. 4, we can observe that the best performance is achieved when the $\alpha$ is 0.11 and when difficult nodes are pruned. This indicates that the improvement of data quality can significantly improve the performance of the model. Compared to pruning easy nodes and randomly pruning nodes, pruning difficult nodes brings higher profit because difficult nodes have a negative impact on the representation of the graph. Furthermore, randomly pruning nodes is better than pruning easy nodes, indicating the effectiveness of our difficulty measurer which can assign reasonable difficulty scores to nodes.