We retrieved RNA-seq transcriptome data of GC cohort from TCGA database¹ on February 28, 2020. These data comprised 375 and 32 samples of GC tissues and marching non-cancer tissues respectively. Also retrieved from the same database were the clinical data of 443 cases, which included gender, age, grade, TNM stages, pathological stage, survival time, and vital status. Subsequently, we merged the RNA-seq data of each sample into a matrix file with a merge script in the Perl language². To change the names of the genes from Ensembl IDs to gene symbols, we utilized the Ensembl database³ in a matrix profile. The mRNAsi and EREG-mRNAsi indices of GC cases in TCGA were acquired from previous studies (Malta et al., 2018). The microarray (GSE29272, GSE27342, GSE26899) results for validation were downloaded from the Gene Expression Omnibus (GEO) database. As the data utilized herein were freely sourced from an open database, approval from the Ethics Committee was not required. As shown in Figure 1, our study design was briefly described in the flow chart.

The beeswarm package in R was employed to compare mRNAsi index in GC tissues versus non-cancer tissues. Similarly, DEGs in the two types of tissues were identified using the limma package (Ritchie et al., 2015). The selection criteria: | log₂ fold change| > 1, P < 0.05 and false discovery rate (FDR) < 0.05. DEGs meeting the criterion were selected for further analysis. Heatmap and volcano plot were drawn using the pheatmap and limma packages, respectively.

The WGCNA package was utilized to build a co-expression network targeting DEGs (Langfelder and Horvath, 2008). The dynamicTreeCut package, doParallel package, fastcluster package, foreach package, GO.db package, Hmisc packages, impute package, matrixStats package, preprocessCore package, and survival package were also used in WGCNA analysis.

We chose mRNAsi and epigenetically regulated mRNAsi (EREG-mRNAsi) as the representative traits to identify the CSC-associated modules and genes. Being a stemness index, mRNAsi, was generated from a group of stemness-associated epigenetically regulated genes. Modules related to the mRNAsi were selected, whose genes were considered to be co-expressed CSC-related genes. Initially, the normal data set and the cases with incomplete data were excluded (Supplementary Figure S1A). Subsequently, based on the gene expression levels of remaining samples, we clustered the data and reduced the outlier. A heatmap was generated to show the global outline of the mRNAsi and the EREG-mRNAsi expression in screened cases (Supplementary Figure S1B). Next, the power-value was chosen to construct a scale-free network based on the Pearson correlation coefficient among genes. The appropriate power-value = 4 was selected based on mean connectivity and scale-free correlation coefficient (Supplementary Figure S1C). Accordingly, we constructed a GeneTree and identified dynamic modules with a minimum size of 60 genes. On the GeneTree, branches of the cluster dendrogram corresponded to distinct gene modules and each piece of the leaves on the cluster dendrogram corresponded to a gene (Supplementary Figure S1D). Through further analysis of modules, the module eigengene (ME) dissimilarity was computed and visualized, then a cut-off (<0.25) was chosen for the module dendrogram and some similar modules were merged (Supplementary Figure S1E). Regarding principal component analysis (PCA), we considered MEs as the principal component of module for every gene. Particularly in some modules, each gene expression model was summarized to have distinct features.

To assess the significance of each module, we calculated the gene significance (GS) and analyzed the interaction between the levels of gene expression and sample characteristics. The calculation of GS was the log₁₀ conversion of the p-value in the linear regression between gene expression and mRNAsi or EREG-mRNAsi (GS = lgp). In addition, the mean GS within the module was defined as Module significance (MS), which was determined to analyze the link between each module and sample characteristics. Among all selected modules, the module with the largest MS was taken as the module that is strongly related to sample characteristic.

Subsequently, GS and module membership (MM, relationship between genes in a given module and their expression profiles) for each gene and set their thresholds for screening key genes in the module as cor. gene GS > 0.5 and cor. gene MM > 0.8.

The interactions between key genes at the level of transcription was analyzed by the Pearson correlation analysis with the corrplot package in R software. The PPI network of key genes was constructed using the online Search Tool for the Retrieval of Interacting Genes (STRING)⁴ (Szklarczyk et al., 2019). The bar-plot showed the nodes in the network with top connectivity. Based on this, we computed the sum of adjacent nodes of every gene in the PPI network. Next, using a bar plot, the genes were classified on the basis of adjacent node number. In addition, the pheatmap package was used to draw a heatmap showing the levels of expressions of key genes, and for plotting box-plots, the ggpubr package was utilized.

The org.Hs.eg.db package was chosen to map the key genes with the Ensemble ID. The clusterProfiler package was performed to carry out GO functional annotation and KEGG analyses so as to explore and determine the potential biological functions of each key gene (Yu et al., 2012). The enrichplot package, colorspace package, stringi package, DOSE package and ggplot2 package were also used and the enriched biological processes (BP), cellular component (CC) and molecular function (MF) were obtained. Statistical significance was set at P < 0.05 and an FDR < 0.05. The bar-plot and the bubble-plot were drawn using R software in order to visualize the top results.

The LASSO Cox regression analysis was conducted by the glmnet package and survival package to choose the most suitable genes for modeling. The LASSO regression is an approach for variable selection in fitting high-dimension generalized linear model. By constructing a penalty function with the LASSO regression, we could get a more refined model to decrease the variable numbers and successfully prevent overfitting. Herein, glmnet package was applied to determine the penalty parameter lambda via the cross-validation and identified the optimal lambda value which corresponded to the minimum value of the cross-validation error mean. Then, we chose the best gene group to construct a risk model and categorized the results into the high-risk or low-risk groups. We calculated the risk score based on a linear combination of the coefficients obtained from the LASSO Cox regression model multiplied with the expression value of each selected gene. The independent prognostic role of the risk model was analyzed using multivariate Cox regression. Finally, in order to offer a quantitative tool for predicting the individual probability of patient prognosis, we used the rms package to establish a prognosis nomogram and draw calibration curve to compare the expected and observed survival probabilities.

We completed the analyses using R software (R Core Team, 2013). Herein, all the cut-offs, comprising mRNAsi, expression levels of key genes, and risk score were the median level of each item. We applied the Wilcox test to assess the difference in mRNAsi scores between GC samples and normal samples. The difference in overall survival in low-versus-high risk score patients was analyzed using Kaplan-Meier analysis and log-rank test. The Kruskal-Wallis test was selected to examine correlation in risk scores versus clinicopathological characteristics. Univariate and multivariate analyses were conducted based on a Cox proportional hazard regression model. All data sets merging was performed with a merge script in the Perl language. Statistical significance was set at P < 0.05.

mRNAsi has been applied effectively to evaluate the tendency of tumor cells to dedifferentiate. As such it is used as a marker for identifying CSCs. A remarkably higher mRNAsi was recorded in GC tissues relative to non-cancer tissues (Figure 2A). The DEGs modulating tumor cell stemness were recognized after examining the RNA-seq data retrieved from the TCGA database. Out of the 6,739 DEGs screened, 5,593 were overexpressed, whereas 1,146 were under-expressed (Figure 2B, Supplementary Figure S2, and Supplementary Table S1).

After the screen of DEGs between GC tissues and non-tumor tissues, we constructed a gene co-expression network for the purpose of identifying the biologically significant gene modules by WGCNA, and to further identify genes strongly linked to GC stemness. In this study, a total of 12 modules were obtained for subsequent analysis (Figure 3A). Module significance (MS) was calculated to analyze the link between mRNAsi scores and gene. Due to that R²-value indicates the stronger the link between GC stemness and gene expression, the nearer the value was to 1. As shown in Figure 3B, three modules were considered the strong correlation to GC stemness, namely brown module, blue module and pink module. The brown module exhibited a positive correlation with mRNAsi (R² = 0.76, P < 0.001) (Figure 3C), while the blue and pink modules reflected a negative correlation with mRNAsi (R² = −0.78, P < 0.001; R² = −0.55, P < 0.001, respectively) (Figures 3D,E). Thus, we chose the brow module for further analyses. The thresholds for selecting key genes in the module were defined as cor. gene MM > 0.8 and cor. gene GS > 0.5. Finally, we checked 17 key genes containing ARHGAP11A, BUB1, BUB1B, C1orf112, CENPF, KIF14, KIF15, KIF18B, KIF4A, NCAPH, PLK4, RACGAP1, RAD54L, SGO2, TPX2, TTK, and XRCC2. The concrete expression values of each key gene were extracted and the corresponding heatmap and boxplot were drawn, showing the upregulation of the key genes in GC tissues (Figures 4A,B). To verify the expression difference of key genes in the TCGA database, we chose the GEO dataset and selected three data sets (GSE26899, GSE27342, GSE29272) with large sample sizes to assess the key genes’ levels of expression in GC. Similar to the results obtained from the TCGA database, upregulation of all the key genes in the GC tissues was observed (Figures 5A–C).

The correlation analyses of the key gene expression were carried out to confirm the relevance of the genes within the brown module. Within the figure, the strength of correlation is displayed on the upper part based on color, whereas the lower part represents the equivalent correlation value (Figure 6A). We found that strong positive correlation among the key genes at the transcription level (Pearson correlation ≥ 0.60). The relationship with the highest correlation was between KIF14 and CENPF (0.88), followed by BUB1 and NCAPH (0.87). BUB1 was highly correlated with NCAPH, BUB1B, and SGO2 (Pearson correlation > 0.80). At the protein level, the interactions between key gene proteins were analyzed using STRING and mapped the PPI network (Figure 6B). The PPI network consisted of 14 nodes and 214 edges, and 6 genes (TTK, TPX2, NCAPH, KIF15, CENPF, BUB1) have the highest node numbers (node number = 15). Except for XRCC2, other key genes showed a closer protein interaction (node number ≥ 7) (Figure 6C).

Regarding gene enrichment, the clusterProfiler package was employed to examine the functional link between the key genes. The GO enrichment analysis included the following three portions: biological process (BP), cell component (CC), and molecular function (MF) (Supplementary Table S2). The principal biological roles of the key genes were nuclear division (GO:0000280), organelle fission (GO:0048285), spindle (GO:0005819), and ATPase activity (GO:0016887) (Figures 7A,B). Regarding signaling pathway enrichment, the enriched pathways included cell cycle (hsa04110) and homologous recombination (hsa03440) (Figures 7C,D). The KEGG results are summarized in Supplementary Table S3.

To investigate the effect of key genes on GC prognosis, 17 key genes were input in LASSO regression to identify robust markers. Then, a prognostic model containing nine genes (BUB1B, NCAPH, KIF15, RAD54L, KIF18B, KIF4A, TTK, SGO2, C1orf112) constructed to evaluate the disease outcome of each patient (Figures 8A,B). We calculated the risk scores of GC patients with the LASSO Cox regression model according to the coefficients of nine genes.

Risk score = (BUB1B × 0.021265698)+(NCAPH × 0.00722087 6)+(KIF15 × −0.100943988)+(RAD54L × −0.133097331) + (KIF 18B × –0.048222701)+(KIF4A × −0.002827726)+(TTK × 0.03503 4586)+(SGO2 × 0.074640650)+(C1orf112 × 0.157469647).

We grouped the patients into two classes (low- and high-risk groups) based on the median risk score. The two groups exhibited substantial difference regarding the survival rate. Specifically, high-risk group had markedly lower 5-year survival rate (23.3%) in comparison to the low-risk group (43.1%) (P < 0.001) (Figure 8C). Notably, the low and high curves exhibited a remarkable intersection at the sixth year. In the first 6 years, the survival probability of high-risk patients was lower in comparison to low-risk patients and in the following 4 years, the survival probability curve was basically flat. The patients who exhibited high risk scores lived beyond 6 years, and, in exceptional cases, they even showed a higher probability of survival. The area under curve (AUC) of risk score was higher than that of clinicopathological characteristics (Figure 8D). The heatmap revealed substantial differences in high-risk versus the low-risk groups in fustat (P < 0.05) (Figure 8E). Also, the risk model performance in GC patients were further assessed by univariate and multivariate Cox regression analysis, which indicated that age, TNM stages, and risk score had a significant association with OS (P < 0.005) (Table 1). Multivariate analyses indicated that the nine-gene risk model is an independent predictor for predicting the disease outcomes GC patients (Table 1 and Figure 8F). A prognostic nomogram was established by incorporating data on risk score and clinical records to offer a quantitative tool that can be applied to predict the probability of relapse for each patient (Figure 9A). Age, gender, grade, TNM stages, and risk score were parameters included in the nomogram. The calibration curve of the prognostic nomogram showed good agreement between prediction and observation (Figure 9B).