Prognostic Prediction Using a Stemness Index-Related Signature in a Cohort of Gastric Cancer

Abstract

Background: With characteristic self-renewal and multipotent differentiation, cancer stem cells (CSCs) have a crucial influence on the metastasis, relapse and drug resistance of gastric cancer (GC). However, the genes that participates in the stemness of GC stem cells have not been identified.

Methods: The mRNA expression-based stemness index (mRNAsi) was analyzed with differential expressions in GC. The weighted gene co-expression network analysis (WGCNA) was utilized to build a co-expression network targeting differentially expressed genes (DEG) and discover mRNAsi-related modules and genes. We assessed the association between the key genes at both the transcription and protein level. Gene Expression Omnibus (GEO) database was used to validate the expression levels of the key genes. The risk model was established according to the least absolute shrinkage and selection operator (LASSO) Cox regression analysis. Furthermore, we determined the prognostic value of the model by employing Kaplan-Meier (KM) plus multivariate Cox analysis.

Results: GC tissues exhibited a substantially higher mRNAsi relative to the healthy non-tumor tissues. Based on WGCNA, 17 key genes (ARHGAP11A, BUB1, BUB1B, C1orf112, CENPF, KIF14, KIF15, KIF18B, KIF4A, NCAPH, PLK4, RACGAP1, RAD54L, SGO2, TPX2, TTK, and XRCC2) were identified. These key genes were clearly overexpressed in GC and validated in the GEO database. The protein-protein interaction (PPI) network as assessed by STRING indicated that the key genes were tightly connected. After LASSO analysis, a nine-gene risk model (BUB1B, NCAPH, KIF15, RAD54L, KIF18B, KIF4A, TTK, SGO2, C1orf112) was constructed. The overall survival in the high-risk group was relatively poor. The area under curve (AUC) of risk score was higher compared to that of clinicopathological characteristics. According to the multivariate Cox analysis, the nine-gene risk model was a predictor of disease outcomes in GC patients (HR, 7.606; 95% CI, 3.037-19.051; P < 0.001). We constructed a prognostic nomogram with well-fitted calibration curve based on risk score and clinical data.

Conclusion: The 17 mRNAsi-related key genes identified in this study could be potential treatment targets in GC treatment, considering that they can inhibit the stemness properties. The nine-gene risk model can be employed to predict the disease outcomes of the patients.

Keywords: LASSO regression; TCGA; WGCNA; cancer stem cells; gastric cancer; mRNAsi; prognosis.

FIGURE 1

Flow chart of the study design. TCGA, The Cancer Genome Atlas; mRNAsi, mRNA expression-based stemness index; DEGs, differentially expressed genes; GC, gastric cancer; WGCNA, Weighted gene co-expression network analysis; GEO, Gene Expression Omnibus; PPI, protein-protein interaction; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic.

FIGURE 2

The mRNAsi and DEGs in GC (375 tumor tissues and 32 non-tumor tissues) based on TCGA database. (A) Differences in mRNAsi in GC tissues vs. non-tumor tissues. (B) volcano plot showing differential expression in GC tissues vs. non-tumor tissues. The upregulated gene is displayed in red dot and the downregulated gene is in blue. In total, 6,739 DEGs were identified, of which 5,593 were upregulated, and 1,146 were downregulated. mRNAsi, mRNA expression-based stemness index; DEGs, differentially expressed genes; GC, gastric cancer; TCGA, The Cancer Genome Atlas.

FIGURE 3

Weighted gene co-expression network analysis. (A) Co-expression module identification in GC. The branches of the cluster dendrogram represent the 12 different gene modules. Each module denotes a collection of co-related genes and was given a unique color. Each piece of the leaves on the cluster dendrogram represent a gene. (B) Heatmap displaying the correlations and significant differences between the gene modules and mRNAsi scores or EREG-mRNAsi. The upper row in each cell represent the correlation coefficient ranging from -1 to 1 of the correlation between a certain gene module and mRNAsi or EREG-mRNAsi. P-values are shown in brackets. Scatter plot of module eigengenes in the brown (C), pink (D), blue (E) modules. Each circle denotes a gene, and the circles in the upper right stand for the key genes in the modules. GC, gastric cancer; mRNAsi, mRNA expression-based stemness index.

FIGURE 4

The differential expression of the key genes in GC (375 tumor tissues and 32 non-tumor tissues) based on TCGA database. (A) Heatmap of the key genes in the two groups. Red indicates upregulation, and green indicates downregulation. (B) Box-plot of the key genes in the two groups. ***P < 0.001; GC, gastric cancer.

FIGURE 5

Validation of the key genes in the GEO microarray database. (A) GSE29272, 268 samples. (B) GSE27342, 160 samples. (C) GSE26899, 108 samples. *P < 0.05, **P < 0.01, *** P < 0.001. GEO, Gene Expression Omnibus.

FIGURE 6

Correlation analysis of the key genes. (A) Correlation among the key genes at transcriptional level. Within the figure, the upper part shows the correlation strength based on the color, whereas, the lower part represents the matching correlation value. (B) Protein-protein interactions of the key genes. Circles denotes genes, lines indicate interactions between gene-encoded proteins, and line colors denote proof of associations between proteins. (C) Edge number of each key gene.

FIGURE 7

Functional annotation and pathway enrichment analysis. (A,B) The GO enrichment analysis of the key genes. Bubble-plots and Bar-plots showing the top terms in groups of BP, CC, and MF. (C,D) The top terms of the KEGG pathway analysis for the key genes are also shown by bar-plot and bubble-plot. Count: Number of genes linked to the enriched GO or KEGG pathway. BP, biological process; CC, cell component; MF, molecular function.

FIGURE 8

Construction of a nine-gene risk model and its prognostic value for GC patients based on TCGA database. (A) Selection of the optimal genes used to construct the final prediction model by LASSO regression analysis. Ten-fold cross-validation for tuning parameter selection. The number on top of the plot represents the total number of genes. Partial likelihood deviance is plotted against log lambda. Dotted vertical lines were drawn at the optimal values. The optimal gene group was chosen by 10-fold cross-validation and the minimal value of lambda. (B) LASSO coefficient profiles of the key genes. The number on top of the plot represents the total number of genes. Each curve represents corresponding key gene and the number next to it is the serial number of each gene. (C) Kaplan-Meier curve of the relationship between risk score and OS of GC patients. (D) The receiver operating characteristic (ROC) curve of the risk model for survival prediction. (E) The heatmap displays the expression of the nine genes and the correlation of clinicopathological parameters with different risk groups. Red indicates upregulation, and green indicates downregulation. (F) Forest plot for the multivariate Cox proportional hazard regression model of the risk score and clinicopathological parameters. GC, gastric cancer; OS, overall survival. *P < 0.05, ***P < 0.001.

FIGURE 9

(A) Prognostic nomogram on the basis of risk score and clinical information. (B) The calibration curve of the prognostic nomogram. Dashed line at 45° represents perfect prediction.