mRNA expression data (normal: 59, tumor: 535) in fragments per kilo-base of exon per million fragments mapped (FPKM) and count formats (normal: 59, tumor:535), clinical data, and single-nucleotide variant (SNV) data (VarScan2 Annotation, sample number: 561) were downloaded from TCGA (https://portal.gdc.cancer.gov/; October 20th, 2020). Dataset GSE72094 was accessed from GEO (https://www.ncbi.nlm.nih.gov/geo/) as the validation set. Raw data were provided by GPL15048 platform.

Differential expression analysis was undertaken on the normal group and tumor group in the training set using “edgeR” package to screen differentially expressed genes (DEGs). The threshold value was set as |logFC| > 1.5 and false discovery rate (FDR) < 0.05 (Robinson et al., 2010). Tumor metabolism-associated gene sets compiled by Possemato et al. (2011) were downloaded from Pubmed (Supplementary Table S1). DEGs were intersected with tumor metabolism-associated genes to obtain DEGs associated with LUAD metabolism. Thereafter, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed on metabolism-associated DEGs using “clusterprofiler” package (q value < 0.05) (Yu et al., 2012).

Samples whose survival time is less than 30 days in TCGA-LUAD were removed. Univariate COX regression analysis was undertaken on metabolism-associated DEGs using “survival” package to obtain survival-related DEGs in LUAD (p < 0.05) (Modeling Survival Data, 2013). To avoid overfitting of the statistical model, “glmnet” package was used to perform LASSO COX regression analysis on the above-screened DEGs (Friedman et al., 2010). Penalty parameter “λ” was selected to remove genes with strong relevance through cross validation to reduce the complexity of the model. Finally, “survival” package was used to undertake multivariate COX regression analysis on the above genes. Prognostic feature genes associated with LUAD metabolism were identified. A risk score model was established, and the risk score was calculated by using the following formula: R i s k   s c o r e = ∑ i = 1 n e x p i * β i (1)

The number of prognostic feature genes associated with metabolism is denoted by n; the expression level of gene i is denoted by exp _i ; the regression coefficient of gene i is denoted by β _i .

The risk scores of patients in TCGA-LUAD were calculated based on the expression levels of prognostic feature genes associated with metabolism. The patients were divided into high- and low-risk groups with median risk score as the threshold value. Survival curves of the two groups were drawn using “survival” package. Receiver’s operative curve (ROC) of patient’s 1-, 3-, and 5 years OS was drawn with “timeROC” package. The area under the curve (AUC) was calculated. The results were validated in the validation set to evaluate the predictive performance of the model (Blanche et al., 2013).

Gene Set Enrichment Analysis (GSEA) enrichment analytics tool was accessed from http://www.gsea-msigdb.org/gsea/index.jsp. The signaling pathway enrichment in high- and low-risk groups was analyzed using GSEA software (p < 0.05) to differentiate biological functions in the two groups. The significance of the enrichment score was analyzed by permutation test (permutation test time: 1,000) (Subramanian et al., 2005).

Tumor mutation bearing (TMB) is defined as the total number of detected somatic cell gene coding errors, base substitutions, errors in gene insertion, or deletions per million bases (Yarchoan et al., 2017). The significance of TMB in the two groups in TCGA-LUAD was analyzed using Wilcoxon test. Mutation genes in the high- and low-risk groups were analyzed, combining SNV mutation data. Waterfall plots of the top 30 gene mutations in the two groups were drawn by R package “GenVisR” (Skidmore et al., 2016).

R package “estimate” was used to assess the stromal score, immune score, and tumor purity in LUAD samples in TCGA. Single simple GSEA (ssGSEA) analysis was performed on 29 immune cells using “GSVA” package to assess the immune infiltration levels of each tumor sample. Differential expression analysis was performed on immune infiltration levels in the two groups using Wilcoxon test (Barbie et al., 2009).

Altogether, 3,591 DEGs were acquired through differential expression analysis on normal and tumor groups in TCGA-LUAD in the training set (|logFC| > 1.5, FDR <0.05), including 2,553 upregulated and 1,038 downregulated genes (Figure 1A). As shown in Figure 1B, 562 LUAD metabolism-associated DEGs were acquired by overlapping DEGs and tumor metabolism-associated gene sets. GO and KEGG enrichment analyses were undertaken on metabolism-associated DEGs in LUAD. GO enrichment analysis showed that these genes were mostly enriched in biological functions including regulation of membrane potential, small molecule catabolic process, organic acid transport, and cellular response to xenobiotic stimulus (Figure 1C). KEGG enrichment analysis showed that these genes were mostly enriched in signaling pathways including the metabolism of xenobiotics by cytochrome P450, retinol metabolism, drug metabolism-other enzymes, arachidonic acid metabolism, and purine metabolism (Figure 1D).

Combining patient’s survival data in TCGA-LUAD in the training set, 562 DEGs associated with metabolism of LUAD were subjected to univariate COX regression analysis. Altogether, 117 genes relevant to survival were acquired (Supplementary Table S2). Optimal penalty parameter “λ” was chosen through cross validation. Eight metabolism-associated prognostic feature genes were acquired (Figures 2A,B). These eight feature genes were subjected to multivariate regression analysis. Lastly, five optimal prognostic feature genes associated with LUAD metabolism were obtained to establish a risk score model (Supplementary Table S3). Protective factors were CYP4B1 and SLC24A4. Hazard ratio (HR) was 0.94 and 0.89. Risk factors were CRIK2 (1.09), ABCC2 (1.05), and glyceraldehyde 3-phosphate dehydrogenase (GAPDH) (1.27) (Figure 2C).

Risk scores of samples in TCGA-LUAD in the training set were calculated. Samples were then divided into high- and low-risk groups according to the median score. Meanwhile, we drew survival status plots, survival curves, and ROC curves of the two groups (Figures 3A–D). Survival analysis suggested poorer survival status in the high-risk group in comparison with the low-risk group. ROC curve showed that AUC values of 1-, 3-, and 5 years survival curves were 0.7, 0.7, and 0.66. The favorable prognosis predictive performance of the model was further proved by survival curve and ROC curve of GSE72094 in the validation set (Figures 3E,F). As shown by heat map of expression levels of five feature genes in the two groups, with the increasing of risk scores, the expression of risk factors (CRIK2, ABCC2, GAPDH) were gradually elevated, while the expression of protective factors (CYP4B1, SLC24A4) was decreased (Figure 3G). Overall, the constructed model could predict the LUAD patient’s prognosis well.

Based on KEGG pathway enrichment analysis, the high- and low-risk groups displayed significant differences in pathways like pyrimidine metabolism, glycolysis gluconeogenesis, P53 signaling pathway, glyoxylate and dicarboxylate metabolism, riboflavin metabolism, and purine metabolism (Figures 4A–F). These pathways were mostly relevant to signaling pathways like cell carbohydrate metabolism pathway, lipid metabolism pathway, and P53 signaling pathway relevant to cell cycle, apoptosis, and aging.

As indicated by Wilcoxon test, high-risk groups exhibited significantly higher TMB (Figure 5A). Further analysis on gene mutation revealed differences in the top30 mutation genes in the two groups (Figures 5B,C). GSEA showed that high- and low-risk groups had differences in the P53 signaling pathway. Combining clinical data and SNV data in TCGA-LUAD and GSE72094 datasets, we acquired mutation of TP53 genes in the two groups in two datasets. Chi-square test indicated that TP53 mutation frequency in the high-risk group was evidently higher than that in the low-risk group in two datasets (p < 0.001, Table 1 and Figures 5D,E.

R package “estimate” was used to evaluate the infiltration levels of stromal cells, immune cells in TCGA-LUAD samples to acquire stromal score, immune score, and ESTIMATE score. Stromal score, immune score, and ESTIMATE score in the high-risk group were evidently lower than those in the low-risk group (Figure 6A). Subsequently, ssGSEA method was used to analyze the immune activity of LUAD samples. Enrichment levels of 29 types of immune cell sets were acquired. Differences in immune infiltration and activity of these 29 cells in the two groups were also compared. Stromal score, immune score, and ESTIMATE score were decreased with the elevation of risk score, whereas tumor purity was increased. the low-risk group showed higher immune infiltration levels (Figure 6B). In detail, immune cells like T helper cells in the low-risk group had higher infiltration levels (p < 0.001, Figure 6C), and most immune function products such as human leukocyte antigen (HLA) had higher expression level (Figure 6D). In summary, the low-risk group showed higher immune activity, which may lead to better prognosis.