A total of 1646 patients’ mRNA expressions with detailed survival information were downloaded from TCGA-LUAD from TCGA database (https://portal.gdc.cancer.gov/)and GEO datasets including GSE31210, GSE30219, GSE68465, GSE72094, TCGA-LUAD, GSE30219, and GSE72094.

After merging multiple data sets across platforms, the normalization of data was acquired. Expression values were log-transformed and we used the “ComBat” algorithm to reduce the possibility of batch effects due to non-biotech bias between different data sets (Johnson et al., 2007).

The gene expression matrix data was uploaded to the CIBERSORT by employing the LM22 signature and 100 permutations were set to obtain the immune cell infiltration matrix (Newman et al., 2015). Meanwhile we used the ESTIMATE package to evaluate the immune and stromal score per LUAD sample (Yoshihara et al., 2013). These samples were stratified and clustered based on the ICI pattern. In the analysis, we used the “hc” method of the unsupervised clustering with Pearson, and both innerLinkage and finalLinkage were set to Ward. D2, executed by using the “ConsensusClusterPlus” R package (Yu et al., 2012). In order to make the classification steady, the maxK was set to three and the reps was set to 100 times.

Based on the ConsensusClusterPlus analysis, these patients were grouped into three ICI groups. Meanwhile, we used the limma package to identify the differential expression genes (DEGs) in those ICI clusters with significance cutoff criteria p < 0.05 (adjusted) and absolute fold-change > 1.5.

Basing on the three ICI clusters, we used consensus clustering (parameters: reps:100, pItem:1, pFeature:1; Ward. D2 and euclidean distance, k:4) using expression values of those DEGs to identify the gene clusters; meanwhile, in order to reduce the relevant variables and unnecessary noise, we used the Boruta package to reduce dimension. When the DGEs were positive the cluster was marked as signature A and when it was negative it was marked signature B. Then we identified the principal component 1 as the signature score by using the PCA analysis. Finally, we constructed an approach that is similar to the gene expression grade index to identify the ICI score of various patients: ICI score = PC1A- PC1B:

To confirm the tumor mutational burden of lung adenocarcinoma, we download the mutation data from cBioPortal (https://www.cbioportal.org/; luad_tcga_pan_can_atlas 2018) and evaluated the entire number of non-synonymous mutations in LUAD. In both the high ICI score and low ICI score groups, we used the “maftool” R package to analyze the somatic alterations, especially the top 25 driver genes with the highest alteration frequency in lung adenocarcinoma.

For further immunotherapy research, we download the IMvigor210 dataset to analyze the value of ICI scores in the PD-1 response predicted, which can be acquired from http://research-pub.gene.com/IMvigor210CoreBiologies with completed information about the response to PD-L1 blockade.

We downloaded a data set including 535 lung adenocarcinoma patients from TCGA database (https://portal.gdc.cancer.gov/cart), and then used TIDE website (Jiang et al., 2018) to predict the cancer immunotherapy response for the 535 patients. The 535 lung adenocarcinoma patients were then divided into two groups according to the TIDE prediction: Response group and Non-response group. Finally, we performed the ICI score analysis for these two groups and compared the results with the TIDE results.

We utilized Graphpad prism and R 4.0.0 software to conduct all statistical analyses. The Kruskal–Wallis test and Wilcoxon test was used to compare two groups and more than two groups. In each dataset, survival curves for subgroups were analyzed by using the Kaplan-Meier plotters. Chi-square test was used to analyze the correlation between ICI score subgroup and somatic mutation frequency, and Spearman analysis was used to calculate the correlation coefficient. p < 0.05 was considered as statistically significant.

Firstly, we used the “ComBat” algorithm to reduce the possibility of batch effects due to non-biotech bias between different data sets (Supplementary Figures S1A, B). Then, the CIBERSORT and ESTIMATE algorithms were employed to identify the immuno-cell infiltration in LUAD tumor tissues. Based on 1646 tumor samples with matched immune cell infiltration (ICI) profiles from the meta-cohort (Array express database: GSE31210, GSE30219, GSE68465 and GSE72094; The Cancer Genome Atlas TCGA-LUAD), unsupervised clustering of LUAD patients into different subtypes was performed as described in the MATERIAL and METHOD section.

To our knowledge, the immune cell infiltration condition is highly associated with the patients’ prognosis. To identify the role of specific immune cells in tumor progression, classification of LUAD based on the immune cell infiltration was employed. Here, three ICI clusters were identified with various infiltrated immune cells (Figure 1A). We also found that the ICI cluster 2 has a favorable prognosis (Figure 1B). Considering the immune landscape of LUAD (Figure 1A), the infiltration percentage of Dendritic cells resting, Dendritic cells activated, Macrophages M2, Mast cells resting, NK cells activated, and T cells CD4 memory resting in the ICI cluster 2 were significantly higher than ICI cluster 1 and 3.

The landscape of immune cell interaction in TME are pictured in Figure 1C. We can observe that the correlation between T CD8⁺ cell and T CD4⁺ cell had more memory resting and activated cells than other infiltration immune cells. Moreover, we also tested the PD-1 and CLAT4 expression level in each ICI cluster by using the Kruskal–Wallis test, and the results showed that cluster 3 has higher expression levels with significant difference (Figures 1D,F). However, we cannot observe any difference between cluster 1 and 2, in which the expression level of PD-1 and CTLA4 is similar.

To confirm the underlying biological functions in different immune-phenotypes, we employed the Limma R package to identify the differential expression genes. Because the clinical information of the TCGA-LUAD cohort is not complete, we mainly used LUAD patients with entire clinical information to further perform the analysis in the following study.

Firstly, we analyzed the aforementioned differentially expressed genes (DEGs) by using the ConsensusClusterPlus package with the unsupervised clustering analysis functions, and these genes can be classified into three genomic clusters: gene clusters A, B, and C. To reduce the noises of cluster genes, we used the Boruta package to reduce the dimension in both signature A and B. Here, we identified 144 genes, which were positively correlated with the gene cluster A. The remaining 64 genes were identified as signature B. All the gene expression in LUAD patients were pictured in Figure 2A, mainly based on the gene cluster types. To analyze by the expression level, gene cluster A is significantly different from the other two types, in which ICI signature gene A is lower and ICI signature B is higher. These results are consistent with the prognosis of LUAD patients based on the gene cluster classification (Figure 2B).

Moreover, we also performed the Gene Ontology (GO) analysis to identify the role of differentially expressed genes. As shown in Figures 2C,D, some signal pathways associated with T regulation processes and cytokines secretion were enriched, for example positive regulation of T cell activation and T cell activation. These results indicated that the T cell regulation related pathways played a critical role in promoting the tumorigenesis.

Previous investigations have proven that the immune system may display different roles with both favorable and adverse outcomes, which is determined by the role of immune cells as either pro-tumor or anti-tumor, as we observed in our previous analysis. As shown in Figure 2E, gene cluster A exhibited higher B cells memory and more Dendritic cells activated, Dendritic cells resting, Macrophages M2, Mast cells resting, Monocytes, and T cells CD4 memory resting. These immune cells displayed the highest level of infiltration, and may evolve into improving the therapy. PD-1 is the critical biomarker for immunotherapeutic intervention and is always expressed on activated T, natural killer and B lymphocytes, macrophages, dendritic cells, and monocytes (Ahmadzadeh et al., 2009). In addition, the expression level of PD-1 was significantly different among the three genomic clusters (Kruskal–Wallis, p < 0.01; Figure 2F). ICI gene clusters of A and B were associated with lower PD1/PD-1 expression levels, while ICI gene clusters C were associated with higher PD1/PD-L1 expression levels. Our investigation showed that the consistency of the immune profiles of different gene groups indicated the different prognostic profile.

We used the principal component analysis to count the total score of both the ICI score A and B to obtain appropriate indicators, and this score was regarded as the summed index of each patient. We utilized the X-tile software to divide the TCGA cohort into high and low ICI scores, which was shown in Figure 3A. Here, we chose six biomarkers (CD274, CTLA4, HAVCR2, IDO1, LAG3, and PDCD1) as immune-checkpoint-relevant signatures, while eight other biomarkers (CD8A, CXCL10, CXCL9, GZMA, GZMB, IFNG, PRF1, TBX2, and TNF) were treated as immune-activity-related signatures. As shown in Figure 3B, these biomarkers are overexpressed in the high ICI group by using the Wilcoxon test. Then, we used the gene set enrichment analysis (GSEA) to explore the underlying biology functions. The GSEA results indicated that the enriched pathways in the high ICI score group was associated with various cell proliferation signaling (Figure 3C) and the low ICI score group was attributed to cell cycle pathways (Figure 3D). Meanwhile, the Kaplan-Meier plot showed that the high ICI score group had a much better survival rate than the low group (Figure 3E).

More evidence has demonstrated that high mutation load (non-synonymous variation) in tumors were associated with increasing infiltration of CD8⁺ T cells, which showed that the tumor load mutation (TMB) may play an important role in the immunotherapy response. Based on the important role of TMB in clinical diagnosis, we aimed to study the correlation between TMB and ICI score to clarify the genetic imprinting of each ICI subgroup. Thus, we investigated the TMB between the high and low ICI scores, which showed that the high ICI score group suffered significantly higher TMB than the low ICI score group (Wilcoxon test p < 0.001, Figure 4A). Meanwhile, the correlation results (Figure 4B) showed that the ICI score was positively correlated with the TMB (R = 0.359, p < 0.001). Further study about TMB in survival was employed, and the result showed that the survival predication of high TMB was better than the low TMB group, to some extent (Figure 4C). In order to consider the association between the TMB and ICI scores, we analyzed the synergistic effect of these scores in prognostic stratification and the result showed that the overall survival of high ICI score (including ICI^high -TMB^high and ICI^high-TMB^low) was much better than the low ICI score (including ICI^low -TMB^high and ICI^low-TMB^low), and the overall survival of ICI^high -TMB^high was comparable with that of ICI^high-TMB^low (Figure 4D), which indicated that high ICI score played a leading role in predicting survival. Meanwhile, the various ICI score groups in different TMB groups had significant survival differences. In conclusion, our study identified that the ICI score may be an important diagnoses index in predicted survival and evaluates the immunotherapy response.

Moreover, we assessed the driver genes’ distribution in both high ICI score and low ICI score using the MAFtools. The top 25 genes with highest frequency were pictured in both the high ICI score group and low ICI score group (Figures 4E,F), in which those genes were significantly different in both groups. These results provide one new analysis approach for studying the mechanism of tumor ICI components and gene mutations in immune therapy.

Immune therapy was applied to improve the tumor therapy outcomes, especially by using PD-1 or PD-L1 specific monoclonal antibody (mAb). However, the remission rate remains lower than expected, rarely exceeding 40 percent. In order to analyze the application of ICI score in the assessment of immune therapy, we use the IMvigor210 cohort, in which these patients received the anti-PD-L1 immunotherapy, to be divided into the high ICI score group and low ICI score group. Then, we assessed the survival rates in the different ICI score groups, and the results showed that the high scoring group had much better survival than the low scoring group (Figure 5B). The objective response rate of anti-PD-L1 therapy was higher in the high ICI score group than in the low ICI group (Figure 5A). We also found that higher ICI scores are correlated with objective response to anti-PD-L1 therapy (Figure 5C).

We also performed a validation set, and we used TIDE website (Jiang et al., 2018) to predict the cancer immunotherapy response for the 535 patients from TCGA database, and our ICI score is consistent with the TIDE index from the aspect of predictive results ((Figure 5D), which can be an effective immunotherapy index for identifying suitable candidates to receive the immunotherapy. In order test whether this ICI score signature can be used as an independent prognostic factor, Cox analysis was used, and the result indicated that ICI score can be a better marker for prognosis (S2).