The mRNA sequencing and clinical data of 535 LUAD samples and 59 normal samples were downloaded from the TCGA data portal. The metabolism-related genes were downloaded from published work (19). The immune-related genes were downloaded from an online website (https://www.immport.org/). Low expressed genes were excluded from the study and the data was normalized to log2 (tpm+1) (average expression after normalization <0.5). Finally, 346 immune-related metabolic genes were selected by cor test using the Pearson correlation method (P<0.05, |R|>0.2).

We divided 346 differentially expressed genes (DEGs) into up-regulated and down-regulated genes. R was used to perform GO and KEGG enrichment analysis. The “clusterProfiler”, “richplot”, and “ggplot2” packages were used for analysis (20, 21). The GO analysis was performed to annotate genes and classify up-regulated and down-regulated DEGs. The GO terms consisted of 3 parts: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). The KEGG database included the systematic analysis, annotation, and visualization of gene functions (22). STRING online website was used to construct a protein-protein interaction (PPI) network for the selected DEGs (23). For PPI analysis, the confidence score was set to > 0.9, and only terms with both p- and q-value of <0.05 were considered significantly enriched. Cytoscape software further analyzed the most closely connected modules and identified the top 10 central genes (24).

The “limma” package (25) in R was used to identify DEGs between Cancer and adjacent tissue samples. Merely genes with | log2fold change | > 1 and P < 0.05 were considered as DEGs. The “pheatmap” package was used to draw heat maps, and “ggplot2” was used to draw volcano maps.

The immune-related metabolic genes were divided into different clusters by the cell consistency clustering method. We used the “ConsensusClusterPlus” package (100 iterations and 80% resampling rate, http://www.bioconductor.org/) to classify patients with LUAD into different subtypes. The heat map and dela diagram established the optimal number of clusters. The cumulative distribution function (CDF) was plotted to identify the number of best clusters. The Progress Free Survival (PFS) between various clusters was compared. The survival analysis was analyzed by the R package “survival”, and the “ggplot2” package was used for plotting.

The expression of immune-related molecules among these clusters with the ESTIMATE algorithm was analyzed by the R “ESTIMATE” package. These immune-related genes regulate four immune functions these included, antigen presentation, chemokine-related genes, cytokines, and immune checkpoints. “ggplot2” package was used to draw box plots.

Four methods were used to assess the infiltration of immune cells in three clusters. These methods were single-sample gene set enrichment analysis (ssGSEA), Microenvironment Cell Populations (MCP)-counter, CIBERSORT, and Xcell (26–29). The three different immune cell infiltrating clusters were also compared and found their immune scores. We compared the pathological classification proportions between different clusters to further distinguish the differences between different clusters, including T, N, M, clinical drug treatment response, and the pathological stage.

The expression matrix of LUAD was randomly divided into training and test sets. The training set was 70% and the test set was 30%. We used single-factor analysis on the two groups of genes. The genes with p<0.05 were selected. The Least Absolute Shrinkage and Selection Operator (LASSO) cox regression method was further optimized through multi-factor COX regression analysis to help us to determine the best number of genes to build a model (30, 31). Moreover, we collected 80 lung cancer samples along with complete survival information and constructed a prognostic model as a control. Finally, the gene’s risk score was screened to get a good predictive ability on the patient’s survival. The area under the ROC curve (AUC) was used to judge the prognostic model’s predictive power. The ten-fold cross-validation based on the “glmnet” package in R was used for lasso penalty Cox regression analysis. The survival analysis was analyzed by R package “survival”, while AUC was analyzed by R package “survivalROC”.

First we selected the immune metabolism genes that were highly expressed in the tumor site (logFC>1.5, FDR<0.05), and the immune metabolism genes in the training set that were negatively related to survival (p ≤ 0.01). Through pan-cancer analysis (https://cistrome.shinyapps.io/timer/) and the expression level analysis of collected clinical samples, we screened out genes with higher expression abundance in tumors. The signaling pathway was determined with criterion (spearman r>2 or<-2, q value< 0.05) using online website cBioPortal (https://www.cbioportal.org/). The selected genes were used to perform GO and KEGG enrichment analyses.

The different cell lines used in this study were, Human normal lung epithelial cell line (BEAS-2B) A549, NCI-H292, and Calu-3. BEAS-2B cells were purchased from Procell Life Science and Technology Co., Ltd. (Wuhan, China). A549, NCI-H292, and Calu-3 were purchased from the cell bank of the Chinese Academy of Sciences (Shanghai, China). The BEAS-2B cell line was cultured in Dulbecco’s Modified Eagle’s medium (DMEM; Gibco, Grand Island, NY, USA). The A549 cell line was cultured in Ham’s F-12K medium (Gibco) The NCI-H292 cell line was cultured in Roswell Park Memorial Institute-1640 medium (RPMI-1640; Gibco), and the Calu-3 cell line was maintained in modified eagle medium (MEM; Gibco). All media were supplemented with 10% fetal bovine serum (FBS; Gibco) and antibiotics (100 units/ml penicillin and 100 ug/ml streptomycin; Gibco). All cells were incubated in a humidified atmosphere of 5% CO2 at 37°C.

Total RNA was isolated from tissues and cells using TrIzol reagent (Gibco) according to the manufacturer’s instructions. The extracted RNA was reverse transcribed into complementary DNA using a reverse transcription kit (Takara, Dalian, China). Quantitative real-time PCR (qRT-PCR) was performed using the SYBR-Green PCR kit (Roche Diagnostics, Indianapolis, IN) on a Step One Plus Real-Time PCR system (Applied Biosystems, Foster City, CA). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as an internal control. The results were analyzed using the 2^-ΔΔCt method. Primers were synthesized by Sangon Biotech (Shanghai, China). All the primer sequences were listed ( Supplementary Table 1 ).

To identify the immune-related metabolic genes, we obtained 1041 immune genes and 1613 metabolic genes. The general research design and flow of the study was shown ( Figure 1 ). The correlation analysis identified 346 immune-related metabolic genes ( Figure S1 ). The GO analysis consisted of three parts: BP, CC, and MF. Our results indicated that the immune-related metabolic genes were significantly enriched in the BP-associated organic acid biosynthetic process, carboxylic acid biosynthetic process, and monocarboxylic acid biosynthetic process. For the CC, the immune-related metabolic genes were mainly enriched in the Golgi, lysosomal, and vacuolar lumens. Furthermore, the MF analysis showed that the immune-related metabolic genes were significantly enriched in cofactor binding, oxidoreductase activity, acting on the CH−OH group of donors, and carboxylic acid-binding. The immune-related metabolic genes were found to be involved in Arachidonic acid metabolism, PPAR signaling pathway, and Biosynthesis of amino acids ( Figures 2A, B ). PPI network was established to further dissect the potential mechanism of these genes ( Figure 2C ). Top10 core genes were identified by Cytoscape plug-in cytoHubba: these genes included, SDC2, GPC3, GPC1, HSPG2, AGRN, GPC2, GPC5, GPC4, GPC6, and VCAN ( Figure 2C ).

The heat map clearly distinguished the immune related metabolic genes in tumor and normal tissues ( Figure S2A ). A total of 141 DEGs were identified (|log2fold change |> 1, P < 0.05). Among these, 72 genes were up-regulated and 69 genes were down-regulated ( Figure S2B ). Then, we performed GO and KEGG enrichment analysis on 141 differential genes ( Figures S2C–F ). It was observed that DEGs were mainly involved in the fatty acid metabolic process; organic hydroxy compound metabolic process and small molecular metabolic process.

Consistent clustering of immune-related genes was performed to unwind metabolic patterns of tumor cells. Tumor samples were divided into different clusters according to the expression patterns of immune-related gens. To determine the optimal cluster number, the cumulative distribution function (CDF) was plotted and three different clusters were identified. Moreover, heat maps were drawn to compare the expression of immune-related metabolic genes among the various clusters ( Figures 3A–C ). Furthermore, the survival status of the three clusters was evaluated by comparing progression-free survival (PFS) and clinicopathological parameters. Our results showed that cluster 2 have prolonged survival in early times ( Figure 3D ). Consistent with these findings, patients in the cluster 2 have lower T, N, M, and stage as well as the status of lymph node metastasis ( Figures S3A–D ). Similarly, a comparatively more proportion of patients had a complete response to treatments in this cluster ( Figure S3E ).

Furthermore, we were interested to determine immune-related genes in these clusters. We determined those immune-related genes that were involved in antigen‐presentation (B2M, HLA-A, HLA-B, HLA-C, HLA-DPA1, HLA-DQA1, TAP1, TAP2), chemokine-related genes (CCL4, CCL5, CXCL10, CXCL13, CXCL9), immune checkpoint genes (CD226, CD274, CD276, CD40, CTLA4, HAVCR2, LAG3, PDCD1) and genes responsible for the production of cytokines (GZMB, GZMH, IFNG, IL2, PRF1, TNF) expressions. We used a box plot for comparison and found that HLA-DPA1, HLA-DQA1, HLA-B, CXCL13 and CD226 had a high expression within C2 cluster ( Figures 4A–D ).

Moreover, we determined infiltration levels of immune cells in these clusters. We used four reported methods (CIBERSORT, MCP-counter, ssGSEA, and Xcell) for this purpose. Two aspects were explored for these clusters i.e. immune effector cells ( Figures 5A–D ) and immunosuppressive cells ( Figures 5E–G ). Our analysis delineated that Cluster 1 had the least infiltration of immune effector cells and immunosuppressive cells ( Figures 5A–G ). This suggests that cluster1 might be the immunologically-cold tumors. Cluster 2 and cluster 3 were found to be enriched in the immunologically-hot tumor immune microenvironment. Both of these clusters were enriched in both immune effector and immunosuppressive cells. Activated B cells, dendritic cells (DC), and monocytes were significantly enriched in cluster 2 ( Figures 5A–G ). We compared the three different clusters and reached their immune score. The results showed that cluster 2 had higher immune and stromal scores ( Figure 5H ).

Next, we were interested in whether these immune-related metabolic genes could be used to predict survival. The LUAD matrix was divided into training (70%) and test (30%) sets. We selected 80 genes having p<0.05 and performed Unicox analysis to compute the regression coefficient ( Figures 6A, B ). Moreover, multivariate regression was performed to calculate the formula. We identified nine genes that were used to construct a prediction model. These nine genes were, TK1, TCN1, CAV1, ACMSD, HS3ST2, HS3ST5, AMN, ADRA2C, and ACOXL ( Figure 6C ). Patients were categorized into high and low-risk groups in training and test sets. A survival curve was plotted according to the clinical information of two groups of patients ( Figures 6D–F ). The results showed that training and test sets with high-risk score patients had a worse Overall Survival (OS) rate than those of low score patients (p <0.0001) ( Figures 6D–F ). The area under the ROC curves of the predictive model for LUAD has the same performance in the first year, third year, and fifth-year (Training set: AUC at one year: 0.83, AUC at three years: 0.72, AUC at five years: 0.71; Test set: AUC at one year: 0.68, AUC at three years: 0.76, AUC at five years: 0.61) ( Figures 6E–G ). Moreover, the prediction model was validated using our clinical specimens. The validation results confirmed that the high-risk score group had a worse survival (p <0.0001) ( Figures 6H, I ). Meanwhile, our prediction model had high accuracy (AUC at one year: 0.74, AUC at three years: 0.83, AUC at four years: 0.78).

As the highly expressed genes in the tumor could be a potential factor to promote tumor growth, therefore, we selected those immune-related metabolic genes whose expression was high in the tumor site (logFC > 1.5, FDR < 0.05). Five potential targets i.e., HMMR, PFKP, RRM2, TCN1, and TK1 were obtained that were negatively correlated with survival rate ( Figure S4 ). The expression of five hub genes in pan-carcinoma was shown ( Figure S5 ). The correlation of the four genes HMMR, PFKP, TCN1, and TK1 with tumor-infiltrating immune cells and the survival curve in lung cancer was shown ( Figure S6 ). Furthermore, we determined the expression levels of five hub genes in our clinical specimens. We found that the expression of RRM2 was higher in tumor tissues ( Figures 7A–E ). Similarly, the expression level of RRM2 was significantly higher in lung cancer cell lines (A549, NCL-H292, and Calu-3) compared to normal lung epithelial cell line (BEAS-2B) ( Figure 7F ). The survival curves of RRM2 and immune cell infiltration in lung cancer patients were determined. The overall survival of lung cancer patients showed that low RRM2 expression had a better prognosis (p=0.000015) ( Figure 7G ), and disease-free survival also suggested that patients with low RRM2 expression had a better prognosis (p=0.019) ( Figure 7H ). The relationship with immune cells showed that RRM2 was associated with tumor infiltration by B cells and Neutrophils ( Figure 7I ).

To determine the function of RRM2, GO and KEGG analyses were performed to find correlated genes with RRM2 ( Figure S7 ). The RRM2-related genes were mainly enriched in catalytic activities, acting on DNA as determined by GO analysis. While KEGG pathway enrichment analysis showed that RRM2-related genes were mainly enriched in cell cycle regulation. Moreover, the correlation between the RRM2 gene and the CDK family of genes was analyzed. This result showed that the RRM2 gene was highly related to the CDK family of proteins ( Figure 7J ). We also found the same results in tumor samples analyzed by qPCR, showing that this gene was associated with CDK family of proteins. Our results delineated that the expression levels of RRM2, CDK2 (r=0.492, p<0.001), CDK4 (r=0.365, p<0.01), CDK6 (r=0.406, p<0.01) and CDK8 (r=0.440, p<0.01) were positively correlated, which means that RRM2 was significantly correlated with cell cycle signaling ( Figures 7K–N ).