Identification of differentially expressed MRGs in LUAD
According to the Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathway-related gene sets, a total of nine hundred and forty-four MRGs were obtained from the gene sets of “c2.cp.kegg.v7.0.symbols”. We matched these genes with the sequence data of LUAD related mRNA in the TCGA database and GSE31210 dataset, and only common mRNAs were used. Considering the cutoff criteria (adjusted P value <0.05 and |log FC| > 1.0), 116 differentially expressed MRGs, which consist of 31 downregulated and 85 upregulated MRGs(Fig.1), were selected for subsequent analysis.
Functional enrichment of the differentially expressed MRGs
To investigate the potential functional implication of these MRGs, 116 differentially expressed MRGs were further analyzed by GO functional enrichment analysis and KEGG pathway enrichment analysis. A total of 431 GO terms and 42 pathways were identified (adjusted P<0.05).The top 30 enrichment GO analysis and top 30 enrichment KEGG analysis were displayed in Figure 2. The top enriched GO terms in biological processes were carboxylic acid biosynthetic process and organic acid biosynthetic process, and those in cellular components were mitochondrial matrix, ficolin-1-rich granule lumen, and ficolin-1-rich granule, in terms of molecular function, genes were mostly enriched in terms of co-factor binding. In the KEGG pathway enrichment analysis, these genes were shown to be significantly associated with signaling pathway related to material synthesis and material metabolism, such as “biosynthesis of amino acids”, “arginine and proline metabolism”, “glycolysis/gluconeogenesis”, “carbon metabolism” and et al.
Establishment of metabolism-related prognostic signature for LUAD
To identify MRGs associated with OS, a univariate Cox proportional hazard regression analysis was initially performed on 116 differentially expressed MRGs in the TCGA database. The result showed that 12 MRGs were significantly associated with the OS (Fig.3a; P <0.001). Of the survival-related MRGs, 10 genes (ALDOA, TPI1, PKM, LDHA, GPI, PFKP, RRM2, TYMS, GNPNAT1, and ENTPD2) were considered risk factors (all P<0.001; HRs, 1.0026-1.1103) and that their overexpression may reduce survival; overexpression of the remaining two genes (CAT and FBP1) (all P<0.001; HRs,0.9747 and 0.9907, respectively) may improve the survival of patients. The Lasso regression analysis was then used to remove MRGs that may be highly related to other MRGs(Fig.3b-c). Furthermore, a prognostic signature model was established based on multivariate Cox regression analysis. Finally, six MRGs were confirmed and applied to establish a metabolism-related signature(Fig.3d). A prognostic model was constructed to evaluate the prognosis of each patient as follows:Risk score = (0.001709×expression value of ALDOA) + (-0.01187×expression value of CAT) + (0.082279×expression value of ENTPD2)+(0.030344×expression value of GNPNAT1)+( 0.003499×expression value of LDHA)+( 0.018476×expression value of TYMS).
Then, the risk score of each patient was calculated according to this prognostic modle. Based on the median risk score, 426 LUAD patients were classified into a high risk group (n = 213) and low risk group (n = 213).The risk score, survival status and gene expression heatmap of these prognostic MRGs are presented in Fig 4a-c. Kaplan-meier log-rank test indicated that patients in the high risk group showed markedly poorer overall survival than those in the low risk group (Fig. 4d). Areas under the curve value of the signature predicting the 1-, 3- and 5-year OS rates were 0.73, 0.703 and 0.854, indicating that this prognostic model exhibited a good sensitivity and specificity (Fig.4e).
Validation of the metabolism-related prognostic signature for LUAD
The GSE31210 dataset including 174 LUAD samples were used for the validation of the metabolism-related signature. According to the median risk score, we divided patients into high risk (n = 78) and low risk groups (n = 96). Consistent with the results derived from the TCGA database, the Kaplan-Meier curve demonstrated that patients in the high risk group exhibited markedly poorer overall survival than those in the low risk group (Fig. 5d). The risk score, survival status and gene expression heatmap of these prognostic MRGs are shown in Fig 5a-c.The AUCs for 1-, 3- and 5-year OS rates were 0.654, 0.705 and 0.725 (Fig.5e). A nomogram for predicting 1-, 3- and 5-year OS of patients with LUAD was constructed with the six prognostic genes that had most significant values in multivariate analysis(Fig.6a). In addition, the calibration curve of the prognostic nomogram demonstrated good agreement between the predicted and observed survival rates for each of OS(Fig.6b-d).
Analysis of these crucial MRGs expression level
To explore the six hub genes at the transcription level, the mRNA expression levels were analyzed using the TCGA database. The results demonstrated that the expression of ALDOA, ENTPD2, GNPNAT1, LDHA, and TYMS in LUAD tissues were all higher than those of adjacent tissues, while the expression of CAT was lower than those of adjacent tissues(Fig.7).To assess the six hub genes at the translational level, the protein expression levels were analyzed using the HPA database. The results showed that the protein level of ALDOA, ENTPD2, GNPNAT1, LDHA, and TYMS were higher in LUAD tissues than in normal tissues, matched their mRNA expression levels(Fig.8).There is no difference between LUAD tissues and normal tissues for the protein level of CAT(Fig.8b).
Clinical value of prognostic signature
Univariate and multivariate Cox regression analysis was conducted to evaluate the independent prediction ability of metabolism-related prognostic signature between the signature and other common prognostic factors, including age, gender, histological grade, pathological stage and TNM stage. Although univariate Cox analysis indicated that pathologic stage, T stage, N stage and our model were markedly associated with OS (Fig.9a; p<0.001), after the multivariate analysis, only metabolism-related prognostic signature(p<0.001) and pathological stage (p<0.007)can be used as an independent prognostic factor (Fig.9b). To further evaluate the clinical value of MRGs, the relationship between MRGs prognostic indicators and clinical features were investigated, and the results indicated that ALDOA, ENTPD2, GNPNAT1, LDHA, and CAT were differentially expressed in patients with various clinical features (Fig.10). To validate the clinical value of the metabolism‑related prognostic signature, the association between the risk score and clinical characteristics were subsequently assessed, and the results demonstrated that high risk scores were positively associated with survival status, gender, N stage, and pathologic stage in patients with LUAD (Fig.10).
Experimental validation
The protein expression levels of ALDOA, ENTPD2, LDHA, TYMS and CAT were investigated in 5 lung cancer cell lines (A549, H460, H1299, H1975, PC9), normal airway epithelial cells (16HBE) as control. The same with the results we have got from bioinformatics, ALDOA, ENTPD2, LDHA, TYMS were significantly increased in 5 lung cancer cell line, comparing with in 16HBE. CAT was significantly decreased in 5 lung cancer cell line, comparing with in 16HBE. The gray value of protein bands were quantified (Figure 11).
The mRNA expression levels of CAT , ENTPD2 and LDHA were also investigated in 6 lung cancer cell lines (A549, H460, H1299, H1975, PC9, Lewis), normal airway epithelial cells (16HBE) as control. The results were also same with protein, ENTPD2 (Figure 12a) and LDHA (Figure 12b) were significantly increased, while CAT (Figure 12c) was significantly decreased, comparing with in 16HBE .
Since ENTPD2 may be a good prognostic marker and therapeutic target for cancer patients, especially those receiving immune therapy[25]. We used ENTPD2 inhibitor POM-1 in 5 lung cancer cell line, found that it could inhibit the formation of colonies in A549 and PC9, decreased colony-forming were also observed in H460, H1299 and H1975, Which were lung adenocarcinoma cell lines (Figure 13b). The protein expression of ENTPD2 in 4 cell lines were confirmed by western blot after use POM-1 (Figure 13a). Most important, we found POM-1 can inhibit the migration of 5 lung cancer cell lines (Figure 13c).