To identify significant genes associated with spinal cord neuronal cell subpopulations, we used three feature selection methods (Mutual Information Coefficient: MIC, Coefficient of Variation Squared: CV2, and Principal Component Analysis: PCA) to assess the significance of 27,998 genes and ranked them according to their contribution values. Genes with importance scores less than or equal to zero were excluded. Next, the machine learning models were combined with IFS to determine the optimal subset of genes (Figures 2A–C). Machine learning models (SVM, RFC, XGBoost, and KNN) were trained using single-cell gene expression matrices (Normalized of raw read count) as input features, based on five-fold cross-validation.

The results from the training dataset showed that MIC combination with XGBoost (ScnML), achieved the optimal prediction performance by using the first 210 genes, with an accuracy of 94.33% (Table 1). Based on the 210 best genes, ScnML also achieved the best performance on the test dataset, with accuracy, precision, recall, and F1_metrics of 94.08%, 94.24%, 94.26%, and 94.24%, respectively (Table 2). It is notable that the four machine learning models, when combined with PCA, also yielded superior predictive performance. To avoid feature selection methods having the same scoring preferences, we compared the top 100 genes scored by the three feature selection methods. As observed from Figure 2D, there is almost no intersection among the top 100 genes selected by MIC, CV2, and PCA, demonstrating the effectiveness of each feature selection method.

To further validate the robustness of the model, receiver operating characteristic (ROC) curves and confusion matrices were used to evaluate the prediction performance of ScnML. We observe that the AUC of the ScnML model is 0.96 (Figure 3A). The confusion matrix validates the predictive performance of the model for each type of spinal cord neural subpopulation, and the low misclassification rate demonstrates the robustness of the model (Figure 3B). In addition, Uniform Manifold Approximation and Projection (UMAP) of 6,000 single cells revealed that the overall performance of the 210 marker genes was significantly better than that of all genes (Figures 3C, D). In particular, the samples from different categories were almost blended together in the clustering process that exploited all genes (Figure 3C). However, employing the 210 optimal genes generates a distinct distribution of cell subpopulations, demonstrating clear clustering findings (Figure 3D). We also performed heat map clustering analysis on the ScnML gene set and obtained excellent clustering results, demonstrating the advantages of machine learning (Supplementary Figure S1).

We performed functional enrichment analysis of the ScnML gene set to explore biological processes related to sci pathophysiology and potential recovery mechanisms. The analysis revealed significant enrichment in genes associated with axon ensheathment, myelination, and the ensheathment of neurons, highlighting the pivotal role of myelin repair and axonal regeneration post-injury (Franklin and Ffrench-Constant, 2008; Lee et al., 2012) (Supplementary Figure S2). Additionally, processes such as glial cell differentiation and gliogenesis were prominently featured, underscoring the importance of glial responses in scar formation and neural tissue remodeling (Sofroniew, 2009). Importantly, our findings also suggest that the regulation of cell-substrate adhesion and leukocyte migration, including myeloid cells, as key components in the inflammatory response and subsequent healing processes (Supplementary Figures S3, S4). The modulation of cell adhesion dynamics is particularly critical, as it influences axonal growth and neural cell interaction with the extracellular matrix, which are essential for effective nerve repair (Zhu et al., 2015).

In addition, we explored the representation of the 210 marker genes in the biological landscape. We identified potential marker genes such as Atp1a2, which is highly expressed in astrocytes; C1qa and Ly86, which are specifically expressed in microglia; and Vtn, which characterizes a subpopulation of endothelial cells (Figure 4A). These genes have been verifiably reported. Furthermore, the use of multiple genes to characterize cellular subpopulations improves accuracy. For instance, Meg3, Snhg11, and Malat1 ensure the identification of neuron subpopulations; Atp1a2, Cst3, and Dbi are highly expressed in astrocytes; and Cst3, C1qb, and Ctss exhibit high expression levels in microglia (Figure 4B).

We analyzed single-cell expression profiles containing all genes and separately, the 210 genes, as the basis for constructing a partition-based graph abstraction (PAGA) to describe the spinal cord neuronal cell bioscape. Both displayed the same topological structure, such as a tight association between microglia and astrocytes, indicating that ScnML screened for key molecular markers and removed redundant information (Figures 5A, B). We utilized Scanpy to compare the expression levels of the top 20 genes in each cell subpopulation with their expression levels in the other five clusters. For example, the expression levels of Cst3, Dbi, and Malat1 in the astrocyte subpopulation were each higher than the combined totals from the remaining five cellular subpopulations. In the neuron cell subpopulation, Meg3, Snhg11, and Malat1 showed high levels of expression, suggesting their potential as marker genes (Figure 5C and Supplementary Table S1). These results indicate that ScnML possesses irreplaceable advantages in processing scRNA-seq data without relying on prior biological knowledge.

The single-cell transcriptome dataset of crush-injured adult mouse spinal cord that support the findings of this study are available in figshare with the identifier (https://doi.org/10.6084/m9.figshare.17702045) (Li et al., 2022). Based on the same processing method used by Liu et al. the raw sequence data were aligned to the mm10 (Ensembl 84) reference genome and cell numbers and unique molecular identifiers (UMIs) were estimated using CellRanger (3.1.0). The 6,000 single-cell transcriptome samples were used to classify six spinal cord injury cell subpopulations, including Endothelial, Astrocyte, Microglia, Neuron, Oligodendrocyte (ODC), and Pericyte cells. These single-cell transcriptome samples were randomly divided into a 4800-sample training set and a 2200-sample testing set with a ratio of 7:3. To construct a stringent and robust benchmark dataset, we applied a filtration criterion, excluding genes with unique feature counts of zero or less. This process yielded a final set of 27,998 genes, each expressed in at least one of the 6,000 cells surveyed.

The Mutual Information Coefficient (MIC) is predicated on the idea that the presence of a relationship between two variables allows for the construction of a grid that effectively partitions their scatter plot, encapsulating the essence of their interaction. To enable equitable comparisons across grids of different sizes, the mutual information values derived from these partitions are normalized. This normalization ensures a consistent framework for evaluating the strength and complexity of relationships between variables, irrespective of their scale or the intricacy of their association (Zhou et al., 2004; Reshef et al., 2011; Albanese et al., 2012). I X ; Y = ∑ x , y p x , y log p x , y p x p y = H X − H X | Y (1) Where I X ; Y representing Mutual Information Entropy, is a measure of the information about variable X   o r   Y contained in variable Y   o r   X .

We performed an extensive analysis to assess the represent capability of 210 marker genes in identifying cell subpopulations. The clustering analyses were performed using the Scanpy software (version 1.9.1), and default parameters were used for all analyses (Wolf et al., 2018). Partition-based graph abstraction (PAGA) was also implemented via Scanpy, while uniform manifold approximation and projection (UMAP) visualizations were generated using the umap-learn Python package (version 0.3.9), with parameters set to default values. Furthermore, functional enrichment analysis was executed employing the enrichGO function from the clusterProfiler package (version 4.6.2).

eXtreme Gradient Boosting (XGBoost) is a highly sophisticated and efficient machine learning algorithm that has gained widespread recognition for its performance in various predictive modeling competitions (Chen and Guestrin, 2016; Wang et al., 2023b). XGBoost has gained prominence for its efficiency and effectiveness in various predictive modeling competitions. It operates by constructing a series of decision trees in a sequential manner, where each subsequent tree aims to correct the errors of its predecessors. This approach enables the model to learn complex patterns in the data, enhancing its predictive accuracy. One of the key strengths of XGBoost is its ability to handle large datasets with speed and precision, making it an ideal choice for our study. In addition, compared to models such as KNN and SVM, XGBoost provides a direct way to evaluate the importance of each input variable.

The four classic metrics were used to quantify the performance of model predictions, including Accuracy, Recall, Precision, and F1_measure, defined as (Fu et al., 2019; Wang et al., 2021b; Joshi et al., 2021; Liang et al., 2021; Wang et al., 2023c; Liu et al., 2023; Qian et al., 2023): Accuracy = T P + T N T P + T N + F P + F N (2) Recall = T P T P + F N (3) Precision = T P T P + F P (4) F 1  measure = 2 * p r e c i s i o n * r e c a l l p r e c i s i o n + r e c a l l (5) Where T P , T N , F P  and  F N represent the numbers of true positives, true negatives, false positives and false negatives, respectively. In addition, ROC was used to evaluate the performance of the ScnML (Zeng et al., 2016; Zulfiqar et al., 2024).