The Predictive Role of Immune Related Subgroup Classification in Immune Checkpoint Blockade Therapy for Lung Adenocarcinoma

Background: In lung adenocarcinoma (LUAD), the predictive role of immune-related subgroup classification in immune checkpoint blockade (ICB) therapy remains largely incomplete. Methods: Transcriptomics analysis was performed to evaluate the association between immune landscape and ICB therapy in lung adenocarcinoma and the associated underlying mechanism. First, the least absolute shrinkage and selection operator (LASSO) algorithm and K-means algorithm were used to identify immune related subgroups for LUAD cohort from the Cancer Genome Atlas (TCGA) database (n = 572). Second, the immune associated signatures of the identified subgroups were characterized by evaluating the status of immune checkpoint associated genes and the immune cell infiltration. Then, potential responses to ICB therapy based on the aforementioned immune related subgroup classification were evaluated via tumor immune dysfunction and exclusion (TIDE) algorithm analysis, and survival analysis and further Cox proportional hazards regression analysis were also performed for LUAD. In the end, gene set enrichment analysis (GSEA) was performed to explore the metabolic mechanism potentially responsible for immune related subgroup clustering. Additionally, two LUAD cohorts from the Gene Expression Omnibus (GEO) database were used as validation cohort. Results: A total of three immune related subgroups with different immune-associated signatures were identified for LUAD. Among them, subgroup 1 with higher infiltration scores for effector immune cells and immune checkpoint associated genes exhibited a potential response to IBC therapy and a better survival, whereas subgroup 3 with lower scores for immune checkpoint associated genes but higher infiltration scores for suppressive immune cells tended to be insensitive to ICB therapy and have an unfavorable prognosis. GSEA revealed that the status of glucometabolic reprogramming in LUAD was potentially responsible for the immune-related subgroup classification. Conclusion: In summary, immune related subgroup clustering based on distinct immune associated signatures will enable us to screen potentially responsive LUAD patients for ICB therapy before treatment, and the discovery of metabolism associated mechanism is beneficial to comprehensive therapeutic strategies making involving ICB therapy in combination with metabolism intervention for LUAD.


INTRODUCTION
Lung cancer is one of the most common type of malignancies worldwide, and is the leading cause of cancer-related death among men and women globally (Siegel et al., 2021). Nonsmall cell lung cancer (NSCLC), which includes squamous cell carcinoma, adenocarcinoma and large cell carcinoma, accounts for more than 80% of all primary lung cancers (Kano et al., 2020). Within NSCLC, adenocarcinoma is the most common histological subtype . Despite great improvements in LUAD treatment in recent decades, particularly molecular-targeted therapeutic strategies, such as tyrosine kinase inhibitors (TKIs) treatment targeting epidermal growth factor receptor (EGFR) and/or anaplastic lymphoma kinase (ALK) (Ge and Shi, 2015), the prognosis for LUAD patients remains poor with a 5-years survival rate of only 15% (Siegel et al., 2021). Fortunately, as an emerging therapeutic approach for tumor, immunotherapy, such as immune checkpoint blockade (ICB) therapy, is increasingly approved to be effective for LUAD (Huang et al., 2020a). Cytotoxic T-lymphocyte antigen 4 (CTLA-4) and programmed cell death protein 1/programmed cell death ligand 1 (PD-1/PD-L1) are crucial immune checkpoints to maintain homeostasis for immune response (Meyers and Banerji, 2020). Actually, attenuated anti-tumor immune response or induced immunosuppression in local tumor microenvironment (TME) partially result from excessive negative immune response mediated by immune checkpoints (Anichini et al., 2020). ICB therapy aims to enhance anti-tumor immune response by inhibiting detrimental immunosuppression induced by immune checkpoint in TME.
Owing to heterogeneity existing in LUAD and development of acquired resistance to ICB therapy, the overall performance of ICB therapy in clinical practice for LUAD is far from satisfactory (Pathak et al., 2020). As one of the most immunological cancer type, immunological surveillance, immunoediting and immune escape play a critical role in LUAD development and progression . Screening for potentially responsive LUAD patients to ICB therapy before treatment by using an effective immunoligical biomarker is beneficial to remarkably improve the outcome of LUAD patients with ICB therapy . Tumor-infiltrating lymphocyte (TIL) score and PD-L1 expression in TME are previously suggested as potential biomarkers to select potentially sensitive subpopulation to ICB therapy prior to treatment and to predict survival for LUAD patients (Gascón et al., 2020;Jin et al., 2020;Hashemi et al., 2021). However, evaluations for the status of TIL and PD-L1 are currently non-standardized and limited by tissue samples availability. A comprehensive analysis of the immune associated signature in TME enable a further understanding of the interplay between local immune status and tumor immunotherapy responsiveness (Park et al., 2020;Wang et al., 2020).
"Omics" techniques which are characterized by high-throughput interfaces are able to investigate complex biological systems in order to identify molecular signatures responsible for the complicated biological phenotype (Gillette et al., 2020;Lazarou et al., 2020). In the present investigation, bioinformatics analyses based on ribonucleic acid (RNA) sequencing (RNA-seq) data and clinical information from Cancer Genome Atlas (TCGA) database were performed to comprehensively explore the predictive role of immune associated signature in therapeutic responsiveness to ICB therapy for LUAD.
FIGURE 1 | The workflow of this study. Briefly, immune related subgroup clustering was performed by using LASSO algorithm and K-means algorithm. After characterization of the immune associated signatures of the identified subgroups, TIDE algorithm analysis was performed to predict the potential sensitivities to ICB therapy. Meanwhile, survival analysis and further Cox proportional hazards regression analysis were also performed for LUAD. In the end, GSEA was performed to explore the metabolic mechanism potentially responsible for immune related subgroup clustering. First, immune related subgroup clustering was performed by using the least absolute shrinkage and selection operator (LASSO) algorithm and K-means algorithm. Second, the immune associated signatures of the identified subgroups were characterized by evaluating the status of immune checkpoint associated genes and the immune cells infiltration. Then, potential responses to ICB therapy were predicted via tumor immune dysfunction and exclusion (TIDE) algorithm analysis, and the relationship between the immune associated signature based on the aforementioned immune related subgroup classification and potential sensitivities to ICB therapy were determined. Additionally, survival analysis and further Cox proportional hazards regression analysis were also performed for LUAD, and gene set enrichment analysis (GSEA) was performed to explore the metabolic mechanism potentially responsible for immune related subgroup clustering. In the end, two microarray data sets from the Gene Expression Omnibus (GEO) database were used as validation cohorts in the study. The work flow of this study was shown in Figure 1. Transcriptomics analysis of the association between immune associated signature and ICB therapy in LUAD not only explains for the heterogeneity in the reactivity to ICB therapy partially from an immunological perspective, but also provide potentially promising biomarker or target to direct sensitive LUAD patients screening prior to ICB therapy and combination therapy strategy making involving ICB therapy in combination with metabolism intervention.

Data Acquisition
The RNA-seq data sequenced on the Illumina RNA sequencing platform for LUAD samples from TCGA samples were download from the Cancer Genomics Browser of the University of California Santa Cruz (UCSC) Xena (https://xena.ucsc.edu/public) (Cline et al., 2013). Then, log2 (x+1) transformed HT-seq counts data and Fragments Per Kilobase Million (FPKM) data were selected for further analysis. The corresponding phenotype and survival information were also downloaded from the UCSC Xena. The latest gene ID annotation file (gencode.v32. annotation.gtf) was downloaded from the GENCODE database (http://www. gencodegenes.org) (Frankish et al., 2019) for Entrez gene ID and Ensembl gene ID transformation. Finally, after matching the TCGA sample ID in RNA-seq with the corresponding phenotype and survival information, a total of 572 LUAD samples in TCGA database were included in the study. Meanwhile, a total of 824 genes directly involved in immunological processes were collected using the Immunome database (Breuer et al., 2013). In addition, Microarray data for 398 LUAD samples in GSE72094 and 442 LUAD samples in GSE68465 were also acquired from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The corresponding gene chip annotation messages and clinical messages of this two data sets were downloaded using the R package GEOquery (Davis and Meltzer, 2007). The clinicopathological characteristics of LUAD patients from the training and validation sets were summarized in Table 1.

Data Preprocessing and Immune Related Subgroup Clustering
RNA-seq data and microarray data for LUAD from public database were first standardized for further analysis. "Combat" algorithm (Johnson et al., 2007) of R package sva (Leek et al., 2012) was employed to reduce the batch effect which may lead to deviations and bias to unrelated biological or scientific differences between subgroups (Leek et al., 2010). To filter out the missing values, intersective genes were selected from the TCGA cohort, GSE72094 cohort, GSE68465 cohort and Immunome database in the current study. Based on the expression of intersective genes for LUAD cohort from the TCGA database, LASSO algorithm and 10-fold cross validation method in R package glmnet (Friedman et al., 2010) were used to select the optimal gene set of the immune associated genes for immune related subgroup clustering. The total within sum of square and average silhouette width were calculated using R package factoextra to identify the best number of clustering. K-means algorithm, a classical unsupervised learning algorithm of artificial intelligence, was used for sample clustering in R software version 3.6.0 (https://www.r-project.org/) by 10 iterations with at least 30 samples for each subgroup. Moreover, consensus matrix analysis was performed in each data set to validate the clustering number, and consensus matrices were generated using the R package ConsensusClusterPlus (Wilkerson and Hayes, 2010). The principal component analysis (PCA) plot of the clustered samples were also drawn in the present study.

Evaluation of Immune Cell Infiltration
Scores and Immune Checkpoint Associated Genes Scores in Tumor Microenvironment as the Immune Associated Signature Immune cell Abundance Identifier (ImmuCellAI) (Miao et al., 2020), a gene set signature-based method, was used to evaluate the infiltration scores of immune cells in the TME of LUAD. ImmuCellAI is capable of precisely estimating the abundance of 24 types of immune cell, including 18 T-cell subsets (CD4 + , CD8 + , CD4 + naïve, CD8 + naïve, central memory T (Tcm), effector memory T (Tem), Tr1, induced regulatory T cells (iTreg), natural regulatory T cells (nTreg), Th1, Th2, Th17, Follicular helper T cells (Tfh), cytotoxic T cells (Tc), mucosal-associated invariant T cells (MAIT), exhausted T cells (Tex), gamma delta T (γδ T), and natural killer T (NKT) cells) and six other important immune cells (B cells, macrophages, monocytes, neutrophils, dendritic cell (DC), and natural killer (NK) cells). In addition, it was reported that ImmuCellAI can estimate the abundance of immune cells with superior accuracy to other methods, especially on many T-cell subsets. Immune checkpoint associated genes, such as CTLA4, CD28, CD80, CD86, CD274 (PD-L1) and PD-1 (PDCD1), were selected from previous relevant studies focusing on the correlation between these genes and LUAD development, progression and prognosis.

Prediction of Potential Sensitivity to Immune Checkpoint Blockades Therapy for Lung Adenocarcinoma Patients Based on Immune Related Subgroup Classification
Tumor immune dysfunction and exclusion (TIDE) algorithm (Fu et al., 2020) was used to calculate the potential possibility to respond to ICB therapy for LUAD patients based on immune related subgroup classification. Generally, TIDE analysis mainly consists of scores for TIDE, immune dysfunction, immune exclusion and several immune associated cells and effector molecules. Among which, negative score for TIDE suggests a lack of immune evasion phenotype. Meanwhile, T dysfunction score shows how a gene interacts with cytotoxic T cells to influence patient survival outcome, and the T cell exclusion score assesses the gene expression levels in immunosuppressive cell types that drive T cell exclusion. Scores for suppressive immune cells, such as cancer associated fibroblasts (CAF), myeloid-derived suppressor cell (MDSC), M2 macrophage indicate immune evasion or immunosuppression, suggesting a low possibility to respond to ICB therapy. Whereas, scores for effector immune cells, associated effector molecular and immune checkpoint associated genes, such as CD8 + T cells, interferon-γ (IFN-γ) and PD-L1 (CD274) represent a potential sensitivity to ICB therapy. Additionally, immune related subgroup clustering, immune associated cells infiltration, immune checkpoint associated genes and clinicopathologic parameters, such as age, gender, pathological TNM stages, tumor stages in LUAD were also evaluated and analyzed between different immune related subgroups to perform a Cox proportional hazards regression analysis.
Gene Set Enrichment Analysis (GSEA) to Explore the Underlying Mechanism Responsible for the Immune Related Subgroup Clustering of Lung Adenocarcinoma GSEA is a bioinformatics analysis to determine whether a prior defined set of genes shows statistically significant and concordant differences between two groups (Sun et al., 2020). GSEA version 4.1.0, was used, the number of permutations was set to 1,000, and FDR <0.05 was the screening threshold. Given a close relationship between glucose metabolism reprogramming in tumor and anti-tumor immunomodulation, glucose metabolism process associated gene signatures, including the process of glycolysis, gluconeogenesis, tricarboxylic acid (TCA) cycle and oxidative phosphorylation (OXPHOS) in mitochondria were compared between the identified immune related subgroups (subgroup 1 vs subgroup 3) to explore the underlying mechanism responsible for the immune related subgroup clustering of LUAD.

Statistical Analysis
The differences of immune associated signatures existed between immune related subgroups, such as the expression of immune check point genes and the infiltration scores of immune associated cells, were evaluated by using Kruskal-Wallis test. Before that, Shapiro-Wilk test and Tukey's test were used to evaluate the status of normal distribution, and F test was used to perform homogeneity tests of variances. In addition, a survival analysis (overall survival) using Kaplan-Meier method was performed for LUAD patients, and the log-rank test was used to compare the differences of survival existed between the immune related subgroups aforementioned. Furthermore, univariate Cox proportional hazards regression analysis was performed to determine the correlation between survival and a variety of factors, including clinicopathologic parameters and immune associated signature factors. Afterwards, significantly associated factors were selected for further multivariate Cox proportional hazards regression analysis to determine independent risk factors. A p-value under 0.05 was considered to indicate a statistically significant difference. Data was analyzed using R software version 3.6.0. Multiple testing was corrected using the Benjamini-Hochberg's false rediscovery rate (FDR).

Immune-Associated Subgroup Clustering for Lung Adenocarcinoma From the Cancer Genome Atlas Database
The LASSO algorithm and 10-fold cross-validation were used to extract the optimal subsets of immune associated genes based on Immunome database for immune related subgroup clustering of LUAD cohort from TCGA database. As shown in Figure 2A, the optimal λ which have the minimum mean square error was selected by 10-fold cross validation. LASSO coefficient profile of the selected subsets of immune associated genes (n 11) at the optimal λ for immune related subgroup clustering of LUAD was depicted in Figure 2B. To optimize the average silhouette width and the total within sum of square, the optimal number of clustering was set with k 3 (Figures 2C,D). Based on this clustering, LUAD cohort (n 572) from TCGA was divided into subgroup 1 (n 252), subgroup2 (n 188) and subgroup 3 (n 132). The consensus matrix and the principal component analysis (PCA) plots of this immune related subgroup classification when Frontiers in Genetics | www.frontiersin.org October 2021 | Volume 12 | Article 771830 k 3 was shown in Figure 2E and Figure 2F, respectively. The results about this clustering were further validated in GSE72094 and GSE68465 data sets (Supplementary Figure S1).

Characterization of Immune Associated Signature Based on Immune Related Subgroup Clustering of Lung Adenocarcinoma
In the current investigation, immune checkpoint associated genes and immune cell infiltration scores were used to represent the immune associated signature of each of the immune related subgroups of LUAD. The status of immune checkpoint associated genes, such as CTLA4, CD28, CD80, CD86, PD-L1 (CD274) and PD-1 (PDCD1), were first evaluated for LUAD based on the aforementioned immune related subgroup clustering. As demonstrated in the heatmap ( Figure 3A), the levels of these immune checkpoint associated genes were significantly different between the three subgroups (p < 0.05). Box plots were also used to show the differences in each of these immune checkpoint associated genes between the three subgroups ( Figure 3B). Generally, subgroup 1 tended to have significantly higher expression levels of immune checkpoint associated genes in comparison with other subgroups, particularly with subgroup 3. Next, immune cell infiltration estimation was performed by using ImmuCellAI. As shown in Figure 3C, the general infiltration score was higher in subgroup1 in contrast with other subgroups, and a total of 16 immune cell infiltration scores were found to be statistically different between the three immune-related subgroups. In detail, the infiltration scores for effector immune cells, such as CD8 + cells and cytotoxic cells, were found to be statistically higher in subgroup 1 than that in other subgroups, whereas CD8 naive cell infiltration score was relatively lower in subgroup 1 compared to other subgroups. Meanwhile, the cell infiltration scores for suppressive immune cells, such as natural regulatory T cells (nTreg) and induced regulatory T cells (iTreg) were significantly higher in subgroup 3 than that in the other subgroups. (Figure 3D). Similar results with regard to the characterization of immune associated signature based on immune related subgroup clustering of LUAD were also validated in GSE72094 and GSE68465 data sets (Supplementary Figure S2 and Supplementary Figure S3).

Estimation of Potential Sensitivity to Immune Checkpoint Blockades Therapy for Lung Adenocarcinoma Based on Immune Related Subgroup Clustering
Tumor immune dysfunction and exclusion (TIDE) algorithm was used to evaluate the potential sensitivity to ICB therapy for LUAD patients included in different immune related subgroups. As shown in the heatmap ( Figure 4A), the TIDE analysis associated scores were significantly different between the three subgroups (p < 0.05). Based on the TIDE analysis, a higher potential sensitivity to ICB therapy was suggested for subgroup 1 which was with higher scores for TIDE, dysfunction, CD8 + cells, and interferon-γ (IFN-γ), but with lower scores for exclusion, M2 macrophage and MDSC in comparison with that in subgroup 3 ( Figure 4B). Similarly, this TIDE analysis results were also validated in GSE72094 data sets (Supplementary Figure S3). GSE68465 data set was not used as validation cohort to perform TIDE analysis and survival analysis because of the lack of information for CD274.

Survival Analysis and Cox Proportional Hazards Regression Analysis for Lung Adenocarcinoma
With regard to survival analysis for LUAD, the Kaplan Meier curves were drawn and the log-rank test was performed in Box plots were also shown to demonstrate the differences in each of these included immune checkpoint associated genes between the three subgroups. Generally, subgroup 1 tended to have higher expression levels of immune checkpoint associated genes in comparison with other subgroups (p < 0.05). (C) ImmuCellAI was used to evaluate the immune cell infiltration scores in the TME of LUAD. As shown in the heatmap, immune cell infiltration scores were found to be statistically different between the three immune related subgroups. The general infiltration score was remarkably higher in subgroup 1 in comparison with other subgroups, particularly with subgroup 3. (D) Box plots were also shown to indicate the differences between the three subgroups with regard to the infiltration scores of several representative immune cells. The infiltration scores of positive immune response, such as CD8 + cells and cytotoxic cells, were significantly higher in subgroup 1 than that in the other subgroups. Whereas, the infiltration scores of negative immune response, such as natural regulatory T cells (nTreg) and induced regulatory T cells (iTreg), were found to be statistically higher in subgroup 3 than that in other subgroups.
Frontiers in Genetics | www.frontiersin.org October 2021 | Volume 12 | Article 771830 this study. As demonstrated in Figure 4C, subgroup 3 tended to have an unfavorable prognosis in comparison with that of subgroup 1. Then, univariate Cox proportional hazards regression analysis was performed to identify the significant factors influencing the overall survival (OS) of LUAD. Among all the included factors, including the clinicopathologic parameters, immune checkpoint associated genes, immune cell infiltration scores, TIDE algorithm scores and immune related subgroup classification, a total of eight factors were proved to be significant risk factors influencing survival of LUAD (Table 2). Afterwards, all the eight factors were included in further multivariate Cox proportional hazards regression analysis to identify independent risk factors for LUAD. As shown in the forest plots ( Figure 4D), immune related subgroup clustering, tumor stage and B cell infiltration were suggested as potential independent factors influencing OS of LUAD (p < 0.05). The results of survival analysis and Cox proportional hazards regression analysis in validation data set (GSE72094) were also shown in Supplementary Figure S4.

Potential Metabolism Associated Mechanism Responsible for Immune Related Subgroup Clustering of Lung Adenocarcinoma
Based on the immune related subgroup clustering (subgroup 1 vs subgroup 3), gene set enrichment analyses (GSEA) was performed on LUAD data set from the TCGA database using the gene sets significantly associated with glucose metabolism, including the process of glycolysis ( Figure 5A), tricarboxylic acid (TCA) cycle ( Figure 5B), gluconeogenesis ( Figure 5C), oxidative phosphorylation (OXPHOS) in mitochondria ( Figure 5D). FDR (Q value) < 0.05 was set as the screening threshold. As shown, the upward parabolas indicated that all the included processes of glucose metabolism was enhanced in subgroup 1 in contrast with that in subgroup 3. Glucose metabolic reprogramming was FIGURE 4 | TIDE analysis and survival analysis for LUAD based on immune related subgroup clustering. (A) TIDE analysis was used to evaluate the potential sensitivity to ICB therapy for LUAD patients. As shown in the heatmap, the TIDE analysis associated scores were significantly different between the three subgroups (p < 0.05). (B) Based on the TIDE analysis, a higher potential sensitivity to ICB therapy was suggested for subgroup 1 which was with higher scores for TIDE, dysfunction, CD8 + cells and interferon-γ (IFN-γ), but with lower scores for exclusion, M2 macrophage and MDSC in comparison with that in subgroup 3. (C) Kaplan Meier analysis was performed to estimate the survival of LUAD. As shown, subgroup 3 tended to have an unfavorable prognosis in comparison with subgroup 1 and subgroup 2. (D) A total of 8 factors were included in further multivariate Cox proportional hazards regression analysis to identify independent risk factors for LUAD after an univariate Cox proportional hazards regression analysis. As shown in the forest plots, immune associated subgroup clustering, tumor stage and B cells infiltration were suggested as potential independent factors influencing overall survival (OS) of LUAD (p < 0.05).
Frontiers in Genetics | www.frontiersin.org October 2021 | Volume 12 | Article 771830 suggested as one of the underlying mechanisms for immune related subgroup clustering of LUAD. The results of GSEA analysis in validation data set (GSE72094 and GSE68465) were also shown in Supplementary Figure S5.

DISCUSSION
As an emerging therapeutic approach for malignancies, tumor immunotherapy, particularly for ICB therapy, is increasingly proved to be effective for LUAD patients (Huang et al., 2020a;Kano et al., 2020). However, a remarkable improvement in overall response rate and prognosis for LUAD patients is still not achieved due to the inherent intertumoral and intratumoral heterogeneity and the development of acquired resistance to ICB therapy (Jin et al., 2020;Song et al., 2020). To address this issue, a promising biomarker which is capable of predicting therapeutic efficiency before treatment is needed to screen potential responsive subpopulation prior to treatment and monitor the therapeutic efficiency during the process of treatment (Meyers and Banerji, 2020). Tumor-immune relationship plays an important role in tumor development and tumor progression, and tumor immune microenvironment (TIM) is widely accepted as a significant factor influencing therapeutic efficiency of ICB therapy Wang et al., 2020). Specifically, tumorinfiltrating lymphocytes (TILs) score (Gascón et al., 2020;Hashemi et al., 2021) and PD-L1 status  were previously suggested as potential biomarkers to be applied in clinical practice for LUAD. However, its translation from bench to bedside is largely limited by the dependence on tissue sample availability and the nonstandardization for evaluation of TIL score and PD-L1 expression.
Though previous studies tried to use immunophenotypic subtype classification based on immune signature to address this issue Wang et al., 2020;Xu et al., 2020), a systematic and comprehensive analysis (Seo et al., 2018;Zhang et al., 2020) is still required to determine the correlation between immune landscape based on immune related subgroup clustering and therapeutic reactivity to ICB therapy, and the underlying mechanism is of necessity to be explored (Huang et al., 2020b;Giannone et al., 2020). Previous studies from Chen YS. et al. (Xu et al., 2020) and Chen KX. et al. (Seo et al., 2018) performed immune related subgroup classification by using computational algorithms, however, an elaborated immune landscape characterization for distinct immune related subgroups were inadequate. Even though results from Xing Y. et al.  and Kim Y. et al. (Xu et al., 2020) suggested a potential implication of immune subtype classification for ICB immunotherapy in lung cancer, a comprehensive analysis of the potential response to ICB immunotherapy for lung cancer, such as TIDE algorithm, was actually lacked. In the present investigation, we focused on both the elucidation of different immune signatures and prediction of potential response to ICB therapy for lung adenocarcinoma based on immune related subgroup clustering by using K-means algorithm, a classical unsupervised learning algorithm of artificial intelligence. More importantly, we conducted GSEA analysis to explore metabolism associated mechanism potentially responsible for immune related subgroup clustering of LUAD, particularly emphasized on the glucometabolic mechanism to shed light on comprehensive treatment strategy involving ICB immunotherapy in combination with glucose metabolism intervention. Three distinct immune related subgroups were classified for LUAD in the current study based on RNA-seq data set from TCGA database (n 572) by using a K-means algorithm. Among the classification, subgroup 1 was characterized by higher levels of immune checkpoint associated genes and higher cell infiltration scores for immune associated effector cells, and tended to be more sensitive to ICB therapy and have a favorable prognosis. Whereas, subgroup 3 with lower levels of immune checkpoint associated genes but higher cell infiltration scores for immune associated suppressive cells was found to be less responsive to ICB therapy and have a poor prognosis. Presumedly, subgroup 1 represented an immune-hot or with an immunocompetent TME with a higher infiltration score and an immunocompetent subtype which was possibly associated with a potential response to ICB therapy and a favorable prognosis. Whereas, subgroup 3 was considered as an immunodeficient or immunosuppressive landscape with a lower infiltration score or with an immunosuppressive subtype, suggesting a potential resistance to ICB therapy and an unfavorable prognosis. With respect to subgroup 2, a median subtype with a mixture of characteristics of subgroup 1 and subgroup 3 was considered. After Kaplan Meier analysis and Cox proportional hazards regression analysis, the immune related subgroup clustering was found to be an independent risk factor influencing the OS of LUAD patients. In the end, the GSEA analysis revealed that the metabolic reprogramming status in LUAD is potentially one of the underlying mechanisms for the distinct immune associated signatures based on the immune related subgroup clustering (Hensley et al., 2016;Faubert et al., 2017;Smolle et al., 2020). The enhanced glucose metabolism in subgroup 1 was consistent with the immune-hot landscape and a relatively immunocompetent subtype, whereas the decreased glucose metabolism in subgroup 3 suggested an immunodeficient landscape and/or an immunosuppressive subtype. Validation LUAD cohorts from external GEO database were also used to confirm the aforementioned results. To sum up, the present investigation provided a deep understanding of the interaction between tumor cells and surrounding immune cells (Kareva and Hahnfeldt, 2013;Speiser et al., 2016) and shed light on an improvement in ICB therapy or derived combination treatment for LUAD involving ICB therapy and metabolism intervention treatment.
As we know that, metabolic reprogramming and immunomodulation are two hallmarks of tumor (Hanahan and Weinberg, 2011). From a metabolic perspective, both tumorigenesis and immunoregulation are intricately associated with metabolic reprogramming. Specifically, the metabolic interplay between tumor cells and infiltrating immune cells significantly contributes to tumor FIGURE 5 | Glucose metabolic reprogramming was suggested as one of the underlying mechanisms for immune related subgroup clustering of LUAD. Based on the immune related subgroup clustering (subgroup 1 vs subgroup 3), GSEA was performed by using the gene sets significantly associated with glucose metabolism, including processes of (A) glycolysis (Normalized Enrichment Score (NES) 2.29, p < 0.01, Q < 0.05), (B) tricarboxylic acid (TCA) cycle (NES 1.86, p < 0.01, Q < 0.05), (C) gluconeogenesis (NES 1.74, p < 0.01, Q < 0.05) and (D) oxidative phosphorylation (OXPHOS) in mitochondria (NES 1.98, p < 0.01, Q < 0.05) (NES 1.98, p < 0.01, Q < 0.05). FDR <0.05 was set as the screening threshold. An upward parabola indicated that the indicated process was enhanced in subgroup 1 in contrast with subgroup 3. The barcode plot indicates the position of the genes in each gene set; red and blue colors represent positive and negative Pearson's correlation with subgroup classification.
Frontiers in Genetics | www.frontiersin.org October 2021 | Volume 12 | Article 771830 progression and tumor immunosuppression. As reported previously, metabolic competition between tumor cells and surrounding immune cells (Chang et al., 2015) and an accumulation of a variety of metabolite caused by metabolic reprogramming (Feng et al., 2017) in TME are partially responsible for immune landscape remodeling. Even though improvement in ICB therapy for LUAD in recent decades, a potential marker for effective stratification of LUAD patients before treatment and a promising target for associated molecular targeted therapy in combination with ICB therapy are expected to bring out breakthrough to clinical management for LUAD. The heterogeneity in metabolism status of LUAD was previously described (Hensley et al., 2016) and further confirmed by metabonomics analysis by investigation from others (Lazarou et al., 2020;Zhao et al., 2020). Additionally, multiomics analysis based on single cell sequencing data also recovered a close correlation between immune status and metabolic reprogramming Xiao et al., 2020;Zhong et al., 2021). Therefore, ICB therapy combined with metabolism intervention is expected to improve the prospect of LUAD treatment.
In spite of the innovation and valuable results mentioned above with respect to this study, a few limitations existing in the current investigation is noteworthy. First, the TCGA database mainly comprises Caucasian population, while validation cohort from GEO database mostly consists of Asian patients, thus racial bias was not inevitable in this study. To attenuate this bias, two external validation cohorts from GEO database were used to validate the results. Then, as actual sensitivity to ICB therapy for LUAD was not available in this study because the clinical information regarding to ICB therapy was mostly not provided in TCGA and GEO databases, only potential reactivity to ICB therapy for LUAD was evaluated based on TIDE analysis. In the end, the correlation between immune associated signature and sensitivity to ICB therapy and underlying metabolic reprogramming-associated mechanism were not further validated by basic research in vitro and clinical investigation in vivo, which is what we aim to do in future.

CONCLUSION
In the current investigation, a novel immune related subgroup clustering by an unsupervised learning model was identified for LUAD. Distinct immune associated landscape based on this clustering was significantly correlated with potential sensitivity to ICB therapy and prognosis for LUAD. GSEA analysis revealed that the heterogeneity in metabolic reprogramming is potentially one of the underlying mechanisms responsible for the correlation between immune landscape and potential reactivity to ICB therapy for LUAD. The immune related subgroup clustering based on the transcriptomics analysis will enable us to screen potentially responsive LUAD patients to ICB therapy. Additionally, metabolism intervention is a promising approach to improve the therapeutic efficiency of ICB therapy for LUAD.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
XY and ZW contributed equally to this work. XY, ZW, and XL conceived and designed the study, conducted statistical analysis, and wrote the original draft, WC, LZ, GY, YC, and JL performed the investigation and data interpretation; WX and XL reviewed and revised the manuscript. All authors read and approved the final version of the manuscript for publication.

FUNDING
This study was supported by the Science and Technology Development Fund of Tianjin Education Commission for Higher Education (2018KJ061).