3.1 Characteristics of GC patients
GC cohorts in TCGA database includes totally 375 stomach adenocarcinoma patients. The clinical characteristics of GC patients were presented in Table S1. We plotted Kaplan-Meier curves of tumor (T), metastasis (M), lymph (N), and stage for GC cohorts (Figure. S1)
3.2 Differentially expressed autophagy-associated genes between GC and adjacent non‑tumor tissues
Totally 375 GC and 32 non-tumor tissues with RNA-seq data were analyzed in this study, we eventually screened 28 differentially expressed autophagy-associated genes, including 18 upregulated genes (ATIC, BID, BIRC5, CASP8, CD46, CDKN2A, CTSB, DRAM1, ERBB2, HSP90AB1, IFNG, IL24, IRGM, ITGA6, NLRC4, SPHK1, VEGFA, VMP1) and 10 downregulated genes (CDKN1A, FOS, GRID2, HSPB8, NKX2-3, NRG2, NRG3, PINK1, PRKN, TMEM74) with |log2FC| > 1 (Figure. 1). We used “ggpubr” package in R software to present the 28 differentially expressed autophagy-associated genes in GC and normal controls [18] (Figure. 2).
3.3 GO and KEGG analyses of autophagy-associated genes
To explore the potential mechanisms of the above 28 autophagy-associated genes in GC, GO and KEGG enrichment analyses were performed in R software (Figure. 3, Figure. 4). The results of GO analysis showed that the top enriched terms for biological processes (BP) include cell growth, neuron death, regulation of neuron death, regulation of response to cytokine stimulus, and regulation of cytokine-mediated signaling pathway (Figure. 3a-3b). As for KEGG analysis, the selected differentially expressed autophagy-associated genes were mainly associated with platinum drug resistance, bladder cancer, apoptosis, p53 signaling pathway, pancreatic cancer, hepatitis B, ErbB signaling pathway, apoptosis-multiple species, IL-17 signaling pathway and endocrine resistance (Figure. 4a-4b).
3.4 Prognostic signature for GC cohorts
A total of 10 autophagy-associated genes (MAP1LC3B, IRGM, GRID2, ATG4D, GABARAPL1, and HSPB8, CTSL, GABARAPL2, CXCR4, and DLC1) were found to be markedly associated with OS in GC patients via conducting a univariate Cox regression analysis (Table 1). Among these 10 survival-related autophagy-associated genes, 9 genes (MAP1LC3B, IRGM, GRID2, GABARAPL1, and HSPB8, CTSL, GABARAPL2, CXCR4, and DLC1) were identified as risk factors (HRs, 1.112-5.214; P < 0.05) and their overexpression predict worse outcome. However, overexpression of the ATG4D (HR=0.654 (0.475-0.903), P < 0.05) may improve the OS in GC patients. Then, these 10 genes were entered into a Lasso regression analysis. Figure. 5a illustrated the regression coefficient of these 10 autophagy-associated genes in GC. As shown in Figure. 5b, when 8 genes (MAP1LC3B, IRGM, GRID2, ATG4D, CTSL, GABARAPL2, CXCR4, and DLC1) were included, this model achieved the best performance. Table 2 listed the functions and coefficients of these 8 genes, which mainly involved in the formation of autophagosomal vacuoles, apoptosis regulation, as well as autophagosome maturation.
Table 2. Functions of 8 survival-related autophagy-associated genes.
No
|
Gene symbol
|
Full name
|
Function
|
Risk coefficient
|
1
|
MAP1LC3B
|
Homo sapiens microtubule-associated protein 1 light chain 3 beta
|
Play a role in the formation of autophagosomal vacuoles
|
0.17160618
|
2
|
IRGM
|
Immunity-related GTPase family M protein
|
Autophagy-related protein
|
0.57323373
|
3
|
GRID2
|
Glutamate receptor ionotropic, delta-2
|
Autophagy-related protein
|
1.21739174
|
4
|
ATG4D
|
Autophagy-related protein 4 homolog D
|
Play a role as an autophagy regulator that links mitochondrial dysfunction with apoptosis
|
-0.67753721
|
5
|
CTSL
|
Cathepsin L1
|
Regulate autophagy and apoptosis
|
0.09947603
|
6
|
GABARAPL2
|
Gamma-aminobutyric acid receptor-associated protein-like 2
|
Essential for a later stage in autophagosome maturation
|
1.86336532
|
7
|
CXCR4
|
C-X-C chemokine receptor type 4
|
Play a role in immune signaling and autophagy
|
0.40946159
|
8
|
DLC1
|
Deleted in liver cancer 1 protein
|
Regulate autophagy and apoptosis
|
0.05139402
|
Table 1 Univariate and multivariate Cox regression analyses of OS in GC patients.
Genes
|
HR(95% CI)
|
P
|
HR(95% CI)
|
P
|
Coef
|
MAP1LC3B
|
1.5880(1.0568-2.3863)
|
0.0260
|
|
|
|
IRGM
|
5.2142(1.4057-19.3403)
|
0.0135
|
|
|
|
GRID2
|
4.6927(1.5849-13.8947)
|
0.0052
|
6.4603(2.1195-19.6910)
|
0.0010
|
1.8657
|
ATG4D
|
0.6545(0.4746-0.9026)
|
0.0097
|
0.6850(0.4966-0.9448)
|
0.0211
|
-0.3784
|
GABARAPL1
|
1.2709(1.0078-1.6028)
|
0.0428
|
|
|
|
HSPB8
|
1.1123(1.0053-1.2307)
|
0.0392
|
|
|
|
CTSL
|
1.2377(1.0045-1.5250)
|
0.0453
|
|
|
|
GABARAPL2
|
2.0348(1.2877-3.2154)
|
0.0023
|
1.9454(1.2091-3.1301)
|
0.0061
|
0.6655
|
CXCR4
|
1.2270(1.0657-1.4127)
|
0.0044
|
1.1716(1.0164-1.3504)
|
0.0289
|
0.1583
|
DLC1
|
1.3163(1.0591-1.6359)
|
0.0132
|
|
|
|
To understand the contributions of the 8 genes to gastric carcinogenesis, we then explored the genetic alteration these 8 genes via the cBio Cancer Genomics Portal. The Firehose Legacy dataset (393 samples) and PanCancer Atlas dataset (434 samples) of GC were both analyzed in this study. Genes of interest are altered in 100 (25%) of 393 sequenced patients (TCGA, Firehose Legacy dataset) (Figure. S2), compared with that altered queried genes were detected in 103 (24%) of 434 sequenced patients (TCGA, PanCancer Atlas dataset) (Figure. 5c). These findings indicated that these 8 autophagy-associated genes play an essential role in gastric carcinogenesis.
The 8 survival-related autophagy-associated genes were finally analyzed by a multivariate Cox proportional hazards model, and 4 candidate genes (GRID2, ATG4D, GABARAPL2, and CXCR4) were screened as the important prognostic indicators for GC patients (Table 1). The formula of risk score = (1.8657 * expression value of GRID2) + (-0.3784 * expression value of ATG4D) + (0.6655 * expression value of GABARAPL2) + (0.1583 * expression value of CXCR4).
3.5 Validation of the model performance
375 GC cases were divided into low-risk and high-risk groups based on the median values of risk score. To evaluate the performance of this model, Kaplan-Meier survival curves were plotted to compare the GC survival in two groups. We found that patients in high-risk group present lower survival than those in low-risk group (5-year OS, 27.7% vs 38.3%; P=9.524e-07) (Figure. 6a). The ROC curve of the predictive model was illustrated in Figure. 6b, with AUC of 0.67. Besides, Figure. 6c-6e showed that the mortality rate of GC patients increases along with the increasing of risk score.
3.6 The nomogram for predicting survival rate of GC
The clinical nomogram was used to quantitatively evaluate the individuals’ risk via integrating several risk factors. In the nomogram, the individuals’ 3-year and 5-year survival rates were assessed by the total points of risk factors (Figure. 7a). The C-index reaches 0.70 (95% CI: 0.65-0.72). Moreover, calibration curves illustrated good concordance between nomogram-predicted survival and actual survival (Figure 7b and 7c), especially for 3-year survival.