To identify key genes associated with subpopulations of breast cancer cells, we employed three feature selection methods F-score, coefficient of variation squared (CV²), and principal component analysis (PCA) to evaluate the significance of 29,733 genes and rank them according to their contribution (Figures 2A–C). Genes with importance scores less than or equal to zero were excluded from further analysis. The CV², PCA, and F-score extracted 22,000 important genes. Subsequently, machine learning models combined with incremental feature selection (IFS) were utilized to identify the optimal subset of genes. Using five-fold cross-validation, the machine learning models (SVM, RFC, XGBoost, and KNN) were trained with single-cell gene expression matrices as input features.

The analysis of the training dataset showed that the combination of F-score and XGBoost model (BreCML) using the top 360 genes achieved the best predictive performance, successfully classifying breast cancer cell subpopulations with 98.92% accuracy (Supplementary Table S3). Notably, significant predictive performance was also achieved when the four machine learning models were combined with PCA. However, BreCML uses only 360 feature genes, while XGBoost combined with PCA uses 20,000 feature genes, meaning that the complexity of the model was greatly reduced. Therefore, BreCML was selected as the classifier by us. To prevent the feature selection methods from exhibiting similar scoring preferences, we compared the top 100 genes ranked by each feature selection method. As demonstrated in Figure 2D, the top 100 genes selected by PCA, CV², and F-score exhibit minimal overlap, thereby validating the distinct effectiveness of each feature selection method. These findings underscore the utility of combining diverse feature selection techniques to enhance the robustness and accuracy of predictive models in breast cancer research.

BreCML demonstrated exceptional performance on the test dataset, achieving outstanding results across several key metrics: accuracy of 98.67%, precision of 99.15%, recall of 99.49%, and F1-score of 99.79% (Table 1). To further evaluate the model’s effectiveness, we assessed its predictive capabilities using Receiver Operating Characteristic (ROC) curves and confusion matrices. The ROC analysis revealed an impressive area under the curve (AUC) of 0.97 for the BreCML model, as shown in Figure 3A. Additionally, the confusion matrix provided a detailed breakdown of the model’s performance across different breast cancer subgroups, highlighting a notably low misclassification rate (Figure 3B). This strong performance underscores the robustness and reliability of the BreCML model in clinical diagnostics.

To evaluate the robustness of the proposed BreCML for breast cancer cell subpopulation prediction, we assessed the performance of BreCML in an independent dataset and compared it with two state-of-the-art methods: eHSCPr, and HelPredictor. To ensure a fair comparison, these models were executed and evaluated using the same independent test set containing 360 genes. As shown in Table 2, BreCML achieved the best performance among all of the tested methods, with an accuracy of 94.78%, precision of 94.98%, recall of 95.54%, and F1-measure of 95.79%. Specifically, compared to other existing methods, the accuracy of our method is higher by 4.66 to 6.53%. The AUC of the three methods is shown in Figure 3C, and the AUC of BreCML is 0.97, which is outperformed by the other prediction models. Furthermore, the confusion matrix further validated the predictive performance of the model for each cell subpopulation, and the low misclassification rate demonstrated the power of the BreCML model (Figure 3D). Therefore, we conclude that our method is more effective than eHSCPr, and HelPredictor in predicting breast cancer cell subpopulations.

Further analysis was performed using Uniform Manifold Approximation and Projection (UMAP) on 4,874 single cells to evaluate the comparative performance of the 360 selected genes against the entire gene set. The results demonstrated that the 360 selected genes significantly outperformed the full gene set in terms of clustering efficiency and specificity. When clustering was conducted using all genes, samples from different categories were almost entirely intermingled, resulting in poor classification outcomes (Figure 4A). In contrast, the application of the top 360 genes produced a clear and distinct distribution of cell subpopulations (Figure 4B). This enhanced clustering not only improved the visual differentiation of categories but also underscored the effectiveness of selecting key marker genes for precise subpopulation identification.

Additionally, we explored the representation of the 360 marker genes across the biological landscape, identifying several key genes that serve as robust markers for specific cell types within the immune system. Notably, genes such as CD4, IL7R, and CD3D were found to be highly expressed in T-cells, underscoring their importance in cellular immunity functions. Similarly, CD68 was predominantly expressed in myeloid cells, while MS4A1 was identified as a characteristic gene of a B-cell subpopulation (Figure 5). These genes have undergone rigorous validation, and their expression patterns have been corroborated by extensive literature, highlighting their biological relevance and utility in cellular characterization.

Using multiple genes to characterize cellular subpopulations significantly enhances accuracy. For instance, GAPDH, RPL22, RPS12, RPS6, and RPS18 were crucial in identifying epithelial subpopulations. Similarly, CD74, HLA-DRA, and HLA-DPA1 exhibited high expression levels in B-cells, while SSR4, B2M, MZB1, HERPUD1, and XBP1 were prominently expressed in plasmablasts. Additionally, TMSB4X, ZFP36L2, HLA-A, and CCL5 were highly expressed in T-cells (Figure 6). This multi-gene approach not only improves the precision of cell type identification but also provides a more comprehensive understanding of the molecular signatures associated with different cellular subpopulations.

Single-cell transcriptome data for human breast cancer were obtained from the National Center for Biotechnology Information (GSE176078) and include 4,874 cells (13). GSE176078 is one of the most comprehensive scRNA-seq dataset specifically focused on breast cancer, including a wide variety of breast cancer subpopulations. This dataset offers detailed transcriptional profiles of thousands of individual cells from multiple patients, making it highly suitable for studying intratumoral heterogeneity, identifying distinct cell subpopulations, and exploring cell-specific biomarkers. The BCL files were demultiplexed and aligned to the GRCh38 reference genome using Cell Ranger Single Cell software v2.0 (10× Genomics). Cell filtering was performed with the EmptyDrops method from the DropletUtils package v1.2.2, applying additional criteria: cells with more than 200 genes and 250 unique molecular identifiers, and a mitochondrial gene percentage below 20%. The dataset included five different cell subpopulations: B-cells (773), Cancer Epithelisl (1,184), Myeloid (897), Plasmablasts (1,020), and T-cells (1,000). The dataset was split into a training dataset and a test dataset at a 8:2 ratio. More dataset details are provided in the Supplementary Table S1. The Python packages Numpy (version 1.21.6), Pandas (version 1.3.5) and Scanpy (version 1.9.1) were used to read and process the data.

To further validate the robustness of BreCML, we collected single-cell transcriptome data of breast cancer from Wu et al. (14) in the NCBI database (GSE158677). The dataset also included five different cell subpopulations: B cells (598), Cancer epithelial cells (600), myeloid cells (601), Plasmablasts cells (600), and T cells (599). This dataset was used as an independent test set to evaluate the BreCML performance (Supplementary Table S2).

In this study, we performed an extensive analysis to evaluate the predictive capability of 360 marker genes in identifying cell subpopulations. For the clustering analysis in Figure 4, UMAP visualization was executed using the python package umap-learn (version 0.3.9), with all settings maintained at default values. For the clustering analysis in Figure 5, we used the Preprocessing and clustering module in Scanpy (version 1.9.1), which facilitated to identify specific subpopulations of cells associated with these marker genes; default parameters were used throughout. Pearson correlation analysis was conducted on five distinct human breast cancer cell populations, based on the expression profiles of the 360 marker genes, using Pandas (version 1.4.4).

Feature-scML is a scalable and friendly toolkit that allows users to comprehensively score and rank each feature in scRNA-seq data. The PCA module of Feature-scML was used to assess the feature importance of each gene. The source code is available at https://github.com/liameihao/Feature-scML.

The F-score can be used to measure the degree of differentiation of features in different categories and has been shown to be a simple and effective method for feature selection. This method significantly improves the interpretability and classification performance of the model while reducing the bias (15). The F-score of the ith feature is defined as (16, 17):

where X i ¯ represents the average of the ith feature of the whole X i ¯ + is the number of positive samples, X i ¯ − is the number of negative samples. X ¯ k , i + , X ¯ i , i − are the ith feature of the kth positive and negative instances, respectively. The larger the F-score value, the stronger the distinguishing degree of the feature among different categories.

The squared coefficient of variation (CV²) is a quantitative statistical method for quantifying technical variation at the gene level and assessing variability in cell biology, and is widely used in the field of single-cell experiments (18). The CV² method operates by calculating the squares of the coefficients of variation and curve-fitting the observations using the generalized linear model (GLM) in the R package statmod,

Extreme gradient boosting (XGBoost) is recognized as an exceedingly complex and efficient machine learning algorithm, widely acknowledged for its outstanding performance in predictive modeling competitions (19). XGBoost attracts significant attention primarily due to its efficiency and effectiveness demonstrated in various competitive settings. The algorithm operates by sequentially constructing a series of decision trees, each designed to correct the errors of its predecessor. This approach allows the model to capture complex patterns in the data, thereby enhancing predictive accuracy. One major advantage of XGBoost is its ability to quickly and accurately process large datasets, making it an ideal tool for our research. Furthermore, compared to other models such as KNN and SVM, XGBoost also offers a direct method for assessing the importance of each input variable.

During the exploratory data analysis, it was crucial to identify key relationships and assign appropriate weights to features to filter out less relevant or weaker information. We employed three feature selection techniques PCA, CV², and F-score—to evaluate and rank the importance of genes in descending order. Genes with weights equal to or below zero were excluded from further analysis. The sorted gene expression profiles of breast cancer cell subpopulations served as input features for training machine learning models. Utilizing the incremental feature selection (IFS) strategy, we formed 12 combinations by integrating the three feature selection methods with four machine learning models: KNN, XGBoost, SVM, and RFC. Grid search was used to determine the optimal parameters for each combination. The optimal gene set for each combination was identified when the accuracy no longer showed improvement with the addition of more genes. Ultimately, the combination of F-score and XGBoost proved to be the most effective and was employed to develop the BreCML model.

The four classic metrics were used to quantify the performance of the model predictions, namely, the accuracy (Acc), recall (Re), precision (Pre), and F1 measure (F1), defined as (20–28):

where TP, TN, FP, and FN represent the numbers of true positives, true negatives, false positives and false negatives, respectively. In addition, the ROC curve was used to evaluate the performance of the BreCML (29–32).