Machine learning for identifying tumor stemness genes and developing prognostic model in gastric cancer

Cell culture

The SGC7901 gastric cancer cells utilized in this investigation were procured from the Chinese Academy of Sciences Cell Bank (Shanghai, China). These cells were maintained in RPMI-1640 medium supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin.

Cellular immunofluorescence

Please refer to our previously published literature for specific experimental procedures [13]. The antibodies CD44 and Ki67 in the experiments were from ABclonal (USA).

Flow cytometry for cell cycle detection

The quantification of various cell cycle phases was conducted using a cell cycle detection kit (HY-K1071, MCE, USA) in accordance with the manufacturer’s instructions. Following treatment, cell samples were subjected to analysis via flow cytometry using a BD C6 plus instrument.

Immune infiltration analysis

Tumor purity was assessed using transcriptomic profiles of gastric cancer cohorts from TCGA with the ESTIMATE algorithm. ImmuneScore, StromalScore, and ESTIMATEScore were also compared between the low- and high-TSRS groups. ssGSEA was performed to analyze 28 immune cells using the GSVA package. We then validated these results with pan-cancer immune infiltration data and XCELL, CIBERSORT-ABS, EPIC, TIMER and etc., algorithms [14].

Identification of TSKGs

In this analytical segment, two distinct algorithms were employed to identify tumor stemness genes within gastric cancer, followed by an intersection analysis of their results. Initially, the WGCNA algorithm, coupled with mRNAsi, was utilized to compute tumor stemness gene modules and subsequently extract genes of interest from these modules. Subsequently, the second algorithm leveraged the MSigDB database to download embryonic stem cell genomes, utilizing them as scoring criteria via the GSVAR software package to score gastric cancer samples sourced from TCGA. Following this, the limma package facilitated the classification of samples into high and low groups, enabling the identification of differentially expressed genes. Finally, the outcomes from both algorithmic approaches were subjected to intersection analysis using the R package venn.

Machine learning

In this segment of the analysis, we applied five distinct machine learning algorithms to identify TSKGs. The utilized algorithms encompass Support Vector Machine (SVM), K-Nearest Neighbors (KNN), Decision Tree (DT), Logistic Regression (Logit), and Random Forest (RF). The top 30 genes, ranked by their importance, were visually represented through bar charts for each algorithm. Subsequently, the efficacy of these results was assessed through accuracy testing, employing methodologies such as confusion matrices and ROC curves. Noteworthy R language packages utilized in this process include caret, DALEX, ggplot2, random Forest, kernlab, and pROC.

Analyses of mutated landscapes

Mutational landscapes of the CDC25A high and low expression groups were mapped by applying the website assistant for clinical bioinformatics (https://www.aclbi.com/static/index.html#/).

Statistical analysis

The unpaired T-test was used to calculate the p-value, and p < 0.05 was considered statistically significant.

Availability of supporting data

The data generated during this study are included in this article and its supplementary information files are available from the corresponding author on reasonable request.

ssGSEA and WGCNA were used to screen the tumor stemness genes of gastric cancer

mRNAsi serves as a reference standard for assessing stemness in tumor samples [15]. Initially, we compared the distribution of values between cancer and paracancer tissues, as illustrated in Supplementary Figure 1A, and observed a significant elevation in mRNAsi levels within cancerous tissue. Subsequently, we conducted an analysis of gastric cancer samples from TCGA to identify DEGs, as depicted in Supplementary Figure 1B, 1C. WGCNA was then employed on the DEGs to delineate the module exhibiting the highest positive correlation with mRNAsi. During the WGCNA analysis, we set the threshold for excluding outlier samples at 12,000 (Supplementary Figure 2A), established the R2 value for achieving a scale-free network at 0.9, and set β to 3 (Supplementary Figure 2B). The height for merging similarity modules was set at 0.35 (Supplementary Figure 2C), and the resultant merged modules were visualized on a Cluster Dendrogram (Supplementary Figure 2D). Subsequently, we identified a positive correlation involving nine modules. This correlation was then depicted on a heatmap with respect to mRNAsi and EREG-mRNAsi (Supplementary Figure 2E). The blue module was selected for gene extraction, yielding 694 genes (Gene Significance = 0.2, Module Membership = 0.2). Tumor stem cells represent a distinct class with stem cell-like properties. To assess and score TCGA samples for stem cell characteristics, we utilized the embryonic stem cell dataset from the Gene Set Enrichment Analysis (GSEA) database and applied the ssGSEA algorithm. Following scoring of each sample, the sample population was stratified based on the median score, as depicted in Figure 1A. The t-Distributed Stochastic Neighbor Embedding (tSNE) algorithm was employed to visualize the classification results, revealing distinct clustering between high and low stemness score groups (Figure 1B). Additionally, we analyzed the differentially expressed genes within the clustered sample groups, as presented in Figure 1C.

Machine learning was used to screen the key genes of tumor stemness

We conducted a Venn analysis of the tumour stemness genes acquired via the algorithms WGCNA and ssGSEA, yielding a total of 182 genes (Supplementary Figure 3). Machine learning has extensive use in identifying significant genes. Therefore, we utilized diverse algorithms to identify key genes associated with tumour stemness. We categorised 182 genes into SVM, KNN, DTS, RF and Logit algorithms. The top 30 genes of importance were plotted in histograms for each algorithm (Figure 2A–2E). Subsequently, we evaluated the results from each algorithm to achieve better accuracy in identifying tumour stemness genes. As demonstrated in Figure 3A–3E, the confusion matrix was utilized to calculate the results. The SVM confusion matrix displayed values of 2 in the tumour diagonal module, indicating the SVM algorithm had the lowest error rate among the five machine learning algorithms evaluated. To ensure the reliability of the findings, we used the ROC curve to authenticate our results. The ROC analysis indicated that the SVM had an optimal AUC score of 0.990 (Figure 3F). Based on these results, we chose to investigate the role of tumour stem cells in gastric cancer using the SVM-calculated results for further in-depth study.

Construction and evaluation of tumor stemness gene prognostic model

Cancer stem cells represent the degree of malignancy of tumors and are associated with poor prognosis of patients. Therefore, it is of great significance to construct a prognostic model related to tumor stemness for accurate treatment of patients’ conditions. The outcomes obtained by SVM were classified into the LASSO algorithm (Figure 4A, 4B). We evaluated the model by using the GEO gastric cancer dataset GSE84437 as a test dataset to validate the results.

Firstly, gastric cancer samples were classified into high and low TSRS groups based on their TSRS. The TSRS and survival status of the samples were then plotted as scatter plots to assess whether survival changed with increasing TSRS. Figure 4C, 4D indicated that the number of deaths increased in the high TSRS group, suggesting that patients in this group may have an unfavorable prognosis. To confirm this hypothesis, we conducted a Kaplan-Meier analysis on patients with high and low TSRS. The outcomes demonstrated that those with higher TSRS suffered significantly shorter survival (refer to Figure 4E, 4F). The tSNE and PCA findings in the testing and training sets indicated that TSRS could effectively differentiate the risk of the sample population (refer to Figure 5A–5D). Consequently, we suggested that TSRS could function as an independent prognostic factor. Next, following COX regression analysis, it was discovered that TSRS displayed a significant correlation with OS in both the training (Figure 6A, 6B) and testing sets (Figure 6C, 6D).

To explore the prognostic significance of TSRS in gastric cancer and its clinical implications, an analysis of TSRS distribution across various clinicopathological parameters was conducted. Our findings revealed notable disparities in the N pathological parameter concerning TSRS (Figure 7A). Subsequent assessment delineated the specific distribution of TSRS within the N pathological parameter (Figure 7B). Remarkably, the low TSRS cohort exhibited a predominant proportion of No patients (73 patients, 43%), while the high TSRS group displayed a higher incidence of N1 patients (55 patients, 32%). Furthermore, the prognostic efficacy of TSRS was evaluated across both the training (Figure 8A) and test sets (Figure 8B). Leveraging results from ROC analysis, a nomogram was constructed to harness the prognostic utility of TSRS, facilitating tailored treatment recommendations for patients (Figure 8C).

Analysis of tumour immune cell infiltration in different TSRS subgroups

The tumor microenvironment (TME) plays a pivotal role not only in the initiation of early tumorigenesis and distant metastasis but also in its dynamic alterations as the tumor progresses [16, 17]. Within this complex milieu, cytokines secreted by immune cells associated with the tumor promote the survival and differentiation of tumor stem cells, among other processes. Conversely, tumor stem cells themselves secrete chemokines that recruit immune cells to the tumor site [18]. In our study, we conducted a comprehensive analysis of immune infiltrating cells in TSRS subgroups.

Employing the ESTIMATE algorithm, we first assessed tumor purity across different TSRS groups (Figure 9A). Subsequently, utilizing ssGSEA, we scrutinized and visualized the distribution of 28 immune infiltrating cells within both high and low TSRS groups through heat map. While no significant variations were observed in the distribution of these cells between the two groups, there was a discernible preference of immune-infiltrating cells within the high TSRS group (Figure 9B). These findings further reinforce the outcomes derived from the ESTIMATE ImmuneScore analysis.

We employed CIBERSORT to quantify immune-infiltrating cell content, uncovering significant associations within the Tumor Stemness-Related Signature (TSRS) cohorts. Specifically, T cell CD4+ memory resting and monocytes exhibited elevated levels in the high TSRS group, while macrophage M0 levels were heightened in the low TSRS group (Figure 9C, 9D).

To corroborate these findings, we utilized data from the pan-cancer database and various algorithms. XCELL analysis revealed significant correlations between monocytes and T cell CD4+ memory resting with TSRS, while macrophage M0 showed no such association (Supplementary Figure 4A). Monocytes consistently exhibited an association with TSRS across multiple algorithms including XCELL, QUANTISEQ, CIBERSORT-ABS, and MCPCOUNTER (Supplementary Figure 4A, 4C–4E).

Notably, macrophage M2 displayed a significant association with TSRS according to outcomes from the XCELL, QUANTISEQ, and CIBERSORT-ABS algorithms. Conversely, in TIMER and EPIC analyses, only macrophages and T-cell CD4+ were linked with TSRS, without identifying specific cellular subpopulations (Supplementary Figure 4B, 4F). These comprehensive analyses shed light on the intricate interplay between immune cell compositions and TSRS, offering valuable insights into tumor microenvironment dynamics.

Immune checkpoints play a critical role in regulating immune responses, ensuring autoimmune tolerance, and moderating the duration and intensity of immune reactions within tumor tissues. Consequently, we delved into the expression patterns of immune checkpoints across different subgroups of TSRS. As depicted in Supplementary Figure 5, all immune checkpoint genes exhibited heightened expression levels in the high TSRS group. These included genes associated with immune evasion, exhaustion, and suppression, alongside those activating immune cell cytotoxicity. These findings suggest robust activation of immune cells within the high TSRS subgroup, potentially predisposing them to immune exhaustion.

Expression analysis of CDC25A in gastric cancer and its relationship with malignant phenotype

To probe the functional role of TSRGs, we conducted an analysis examining their expression and correlation with mRNAsi. The scatter plot unveiled a notably stronger correlation between CDC25A and the tumor stemness index compared to ASCL2 (Figure 10A, 10B). Moreover, the distribution patterns of stemness scores and gene expression underscored CDC25A’s closer association with escalating levels of mRNAsi compared to ASCL2 (Figure 10C, 10D). Consequently, CDC25A emerged as the prime candidate for further investigation.

CDC25A expression exhibited a noteworthy upregulation in tumor tissue samples, as depicted in Figure 11A, 11B. Additionally, clinical heatmap analysis indicated a significant gender-based disparity in CDC25A expression, as illustrated in Figure 11C. Subsequent analysis revealed a marked elevation in CDC25A expression among female patients, as depicted in Figure 11D. Given the presence of mutations in multiple genes within cancer cells, somatic mutations in CDC25A were investigated, revealing the presence of five mutations - one nonsense and four missense - in gastric cancer cells, as shown in Figure 12A. Further evaluation of CDC25A’s mutation status within the high and low expression groups uncovered a high prevalence of mutations in TP53 and other oncogenes among the high expression group of CDC25A, as illustrated in Figure 12B. This finding suggests a potential association between the malignant phenotype of the tumor and the heightened expression of CDC25A.

To elucidate the specific biological role of CDC25A within gastric cancer pathways, we conducted an enrichment analysis. Initially, we utilized the median expression value of CDC25A to identify differential genes in gastric cancer samples, as depicted in Figure 13A, 13B. Subsequent pathway enrichment analysis unveiled CDC25A’s involvement in the regulation of cell cycle and cell division in gastric cancer, as shown in Figure 13C. Functional enrichment analysis further underscored CDC25A’s regulatory role in critical cellular processes such as mitosis, sister chromatid segregation, spindle assembly, and DNA replication, indicative of the self-renewal and replication abilities of tumor stem cells, as depicted in Figure 13D.

To validate these findings, we assessed the co-expression of genes with CDC25A CDC25A (Figure 14A). Notably, cell proliferation markers PCNA and MKI67 exhibited significant correlations with CDC25A, as depicted in Figure 14B, 14C. Furthermore, proteins OCR1 and ORC6, pivotal for DNA replication, demonstrated a high correlation with CDC25A expression, as shown in Figure 14D, 14E. Lastly, we conducted a correlation analysis to explore the relationship between CDC25A expression and signaling pathways associated with the cell cycle and proliferation. The results indicated a strong and positive association between CDC25A expression and DNA replication, tumor proliferation, G2M checkpoint, and MYC signaling pathways, as depicted in Supplementary Figure 6A–6D.

CDC25A possesses the capability to regulate the proliferation of cells and the stemness of tumours

In light of the bioinformatic analyses conducted, in vitro studies were carried out to verify the effects of CDC25A on gastric cancer cell proliferation, cell cycle, and tumor stemness. The CCK8 experimental findings indicated that suppressing CDC25A expression led to a noteworthy decline in the viability of gastric cancer cells by 0.4467 ± 0.1209 at 48 h (Figure 15A). The cell cycle results indicated that inhibiting CDC25A expression led to an increase in the proportion of cells in the G0/1 phase and a decrease in the number of cells in the S phase, while the proportion of cells in the G2 phase remained unchanged (Figure 15B). Co-expression analysis unveiled a significant positive correlation between CDC25A and MKI67. This finding was substantiated by our immunofluorescence results, demonstrating a marked decrease in MKI67 expression subsequent to CDC25A inhibition. (Figure 15C). To explore the impact of CDC25A on tumor stemness within gastric cancer, we selected CD44 as a well-established marker for cancer stemness investigation. As illustrated in Figure 15D, inhibition of CDC25A expression led to a notable decrease in CD44 expression. These experimental findings provide compelling evidence that CDC25A possesses the capacity to modulate both cell proliferation and tumor stemness in the gastric cancer.