Data on mRNA sequencing (Fragments Per kilobase of exon model per Million mapped fragments, FPKM) and clinical data of LUAD were downloaded from TCGA as the training set for the next-step analysis. The mRNA sequencing (FPKM) and clinical data of GSE101929, GSE50081, GSE41271, and GSE42127 were downloaded from the Gene Expression Omnibus (GEO) platform. The batch effects between GEO datasets were corrected with the R package SVAR (Irizarry et al., 2003). The processed data were used as the test set for the subsequent analysis. The mRNA sequencing and clinical data of GSE13522 and GSE126044 were also downloaded to evaluate the predictive power of the PRS model for an immunotherapeutic response. The somatic mutation data for the training set were downloaded and analyzed using the R package maftools (Mayakonda et al., 2018). The TMBs and mutation rates of LUAD-related driver genes were calculated. The list of driver genes was derived from Integrative Onco Genomics (https://www.intogen.org/search).

A total of 42 patients (referred as NJDT patients) with stage I or II LUAD, who underwent surgeries in Nanjing Drum Tower Hospital from January 2017 to January 2018 were randomly selected. Paraffin-embedded samples of tumor and normal tissues were collected. Sections of the paraffin-embedded tissues were stained using hematoxylin–eosin and examined by two pathologists. The samples were graded and classified according to the Eighth Edition of TNM Classification for Lung Cancer proposed by IASLC (Goldstraw et al., 2016). mRNA high-throughput sequencing was performed on tumor and matching normal samples, and the FPKM data was used for follow-up analysis.

The Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression Data (ESTIMATE) algorithm was used to evaluate the stromal and immune components of samples in the training set and the stromal score, tumor purity, and immune score were calculated (Yoshihara et al., 2013). Based on signal sample Gene Set Enrichment Analysis (ssGSEA) using the R packages of gsva (Hanzelmann et al., 2013) and GSEABase (Reimand et al., 2019), 24 types of TIICs were classified (Bindea et al., 2013): innate immunity (natural killer cells [NKs], NK CD56^dim cells, NK CD56^bright cells, dendritic cells [DCs], activated DCs [aDCs], immature DCs [iDCs], plasmacytoid dendritic cells [pDCs], neutrophils, eosinophils, mast cells, and macrophages) and adaptive immunity (B, T, T helper 1 [Th1], Th2, T gamma delta [Tgd], CD4⁺ T, CD8⁺ T, T central memory [Tcm], T effector memory [Tem], T follicular helper [Tfh], Th17, regulatory T [Treg], and cytotoxic cells). The training set was clustered hierarchically into high (Immunity_H), medium (Immunity_M), and low (Immunity_L) immunophenotype groups. Then, the CIBERSORT algorithm was used to calculate the relative content of each immune cell subset among 22 types of leukocyte subsets (LM22 signature) with 1,000 permutations (Newman et al., 2015). When the p value of the output for each subset was <0.05, the relative contents were considered accurate and suitable for further analysis.

For genes with multiple probes, the average of the probes was used as the gene expression. The R package limma (Ritchie et al., 2015) was used to identify DEGs between normal and tumor samples (DEGs_NT) in the training set. DEGs between Immunity_H and Immunity_M (DEGs_HM) and DEGs between Immunity_H and Immunity_L (DEGs_HL) groups were also screened in the same method. DEGs were defined by the false discovery rate (FDR) < 0.05 and Log2|FoldChange| > 1. The intersection of DEGs_NT, DEGs_HM, and DEGs_HL was used to determine IDEGs. Differential pathways were enriched using Gene Set Enrichment Analysis (GSEA). With the |normalized enrichment score (NES)| >1, nominal p value < 0.05, and FDR <25%, the enrichment was considered significant.

Univariate Cox regression was used to analyze the correlation between IDEGs and overall survival (OS); genes with p < 0.05 were screened. Then, the above genes were analyzed by LASSO regression (Gao et al., 2010) and lambda (λ) values were calculated. Based on the λ value, which corresponded to the minimum mean standard error in the cross-validation, variables were obtained and regression coefficients were calculated. The regression coefficients multiplied by the mRNA levels of 14 genes were used to construct the formula. The median risk score in the training set was used as the grouping cut-off value. Patients with a risk score greater than the cut-off value were classified into the high-risk group; the rest were classified into the low-risk group. Meanwhile, the test set was divided into high- and low-risk groups by using the same cut-off value. The OS curves of the patients in the two sets were plotted, and Log-rank test was used to analyze the differences. The receiver operating characteristic curves (ROCs) of OS in the two sets were plotted, and the areas under curves were calculated to evaluate the predictive performance of the model. Multivariate Cox regression analysis was performed to construct nomograms in both sets.

The R package SomaticSignatures (Gehring et al., 2015) was used to identify and cluster de novo mutation signatures. The number of these signatures was determined by explained variance and residual sum of squares (RSS). The best number of de novo signatures was chosen for clustering. De novo signatures were then compared to 30 curated signatures in the Cancer Gene Census (COSMIC) by using cosine similarity (Cui et al., 2020), Cochran-Armitage trend test was used to examine the mutation signature contribution among groups.

Immunophenoscores (IPSs) were calculated according to the recently published reports (Charoentong et al., 2017; Hakimi et al., 2016). In brief, consensus determinants including 20 single factors and 6 cell types were divided into four categories: effector cells, suppressive cells, MHC-related molecules, and checkpoints or immunomodulators. The Z scores of the determinants included in the particular category were positively weighted with one and negatively weighted with one. The weighted averaged Z score was then calculated by averaging the Z scores within the respective category leading to four values. The IPSs were calculated on an arbitrary scale of 0–10 based on the sum of the weighted average Z scores of the four categories.

The workflow of this study is shown in Figure 1.

All statistical analyses were conducted with the R software (version 4.0.2). The Wilcoxon test was used to compare continuous variables in two groups. The Kaplan-Meier plotter was employed to generate survival curves for the subgroups in each dataset. The Log-rank test was used to evaluate significant differences in survival. The Chi-square test or Fisher’s exact test were used to analyze the clinicopathological categorical variables between the different PRS subgroups. Spearman correlation analysis was used to compute the correlation coefficient between indicators. The multiple hypothesis test with the Benjamini–Hochberg method was used to control FDR. All statistical tests were two-sided, and p values less than 0.05 were considered statistically significant (∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001).

We analyzed the contents of 24 types of TIICs in both sets and evaluated the results by principal component analysis (PCA). There were significant differences between tumor and normal samples. The differences could be used to distinguish normal and tumor tissues (Supplementary Figures S1A–C). The contents of adaptive immune cells increased in tumor tissues, while those of innate immune cells decreased (Wilcoxon test, p < 0.05) (Supplementary Figures S1B–D).

Furthermore, the ESTIMATE algorithm was used to evaluate mRNA profiles of tumor samples in the training set. The OS of the patients in the high score (greater than the median value) group based on the immune scores was higher than those of the patients in the low score group, and the intergroup difference was significant (Log-rank test, p < 0.05) (Supplementary Figure S2A). It indicates that the prognoses of patients with high immune scores are better than those of the patients with low immune scores. Therefore, hierarchical cluster analysis was performed on the TIICs in tumor samples (Supplementary Figure S2B). According to the immune scores, three clusters were defined as high (Immunity_H), medium (Immunity_M), and low immunophenotypes (Immunity_L) (Figure 2A). The OS of the three immunophenotype groups was statistically different (Log-rank test, p < 0.05) (Figure 2B). The patients in the Immunity_H group had the better OS than others (Log-rank test, p < 0.05) (Figure 2C). The TIICs in each immunophenotype were further compared. The levels of mature immune cells in the Immunity_L group were the lowest. Almost all innate immune cells in the Immunity_H group were more than those in the Immunity_M group, except Tfh, CD8⁺ T, and Treg cells. These three kinds of cell increased in the Immunity_M group (Wilcoxon test, p < 0.05) (Supplementary Figures S3A,B). We also used the CIBERSORT method to quantitate TIICs in each immunophenotype group. Twenty-two types of immune cells were quantified; however, the number of CD4⁺ T naive cells was 0 in all samples. Hence, only 21 types of immune cells were finally analyzed. The contents of most innate immune cells in the Immunity_H group were higher than those in the Immunity_M group (Wilcoxon test, p < 0.05) (Supplementary Figure S4A); however, the numbers of CD8⁺ T and Tfh cells in the Immunity_H group were lower than those in the Immunity_M group (Wilcoxon test, p < 0.05) (Supplementary Figure S4B).

We analyzed the clinical features of patients in the three immunophenotype groups (Figure 3A, Supplementary Table S1). The proportion of female patients in the Immunity_H group was the highest. The patients in the Immunity_H group had the lowest TMB and the highest IPSs (Wilcoxon test, p < 0.05) (Supplementary Figure S5A). We also compared HLA expressions and checkpoints in the three immunophenotype groups. The levels of PD-L1, PD-1, FASL, CTLA4, and CD244 in the Immunity_M group were higher than those in the Immunity_H group (Supplementary Figure S4C). In the Immunity_L group, the expression levels of HLA and checkpoints were lower than those in the Immunity_H group (Wilcoxon test, p < 0.05) (Supplementary Figure S4D).

GSEA analysis of the differential enrichment pathways of Kyoto Encyclopedia of Genes and Genomes (KEGG) showed that, the pathways of nucleotide sugar metabolism, DNA stability and autophagy regulation were significantly upregulated (Supplementary Figure S5B), but immune-related pathways and cell adhesion were significantly downregulated of tumor samples in all three immunophenotypes (Supplementary Figure S5C). The Immunity_H group had the most obvious upregulation of metabolism pathways, such as glucose, lipid and water salt metabolism and lysosome pathways. The Immunity_M group was more strongly related to homologous recombination, DNA replication and repair and gene transcription. The Immunity_L group was specifically associated with lowered immune-related pathways, including B and T cell receptor signaling pathways, NK cell mediated cytotoxicity, cytokine receptor interaction and complement-related pathways (Supplementary Figures S5D,E). Thereafter, we intersected DEG_HL, DEG_HM and DEG_NT and obtained 421 IDEGs for subsequent screening (Supplementary Figure S5F). The molecular function and biological processes of these genes covered immune response, regulation of gene silencing, glucolipid metabolism, cell adhesion and blood coagulation (Supplementary Figure S5G).

Cox regression analysis was performed for the candidate genes among the IDEGs that were specifically associated with OS (Log-rank, p < 0.05), followed by LASSO logistic analysis. The most suitable tuning parameters (λ) and coefficients were calculated by cross-validation (Supplementary Figure S6A). Finally, 14 IDEGs were selected to construct the PRS model. The formula was as follows: Prognostic Risk Score = (−0.0461 × TLR8 mRNA level) + (0.0992 × FGF2 mRNA level) + (0.0467 × F12 mRNA level) + (0.3515 × ST6GALNAC3 mRNA level) + (0.0198 × PTPRH mRNA level) + (0.0368 × EXO1 mRNA level) + (0.0182 × FRMD3 mRNA level) + (0.1891 × E2F7 mRNA level) + (−0.1644 × ABHD6 mRNA level) + (-0.0423 × STK32A mRNA level) + (−0.0203 × COL4A3 mRNA level) + (0.0178 × PLEK2 mRNA level) + (−0.0222 * LIFR mRNA level) + (0.04453 × CYS1 mRNA level). The median score in the training set was considered as the cut-off value, and the patients were divided into high-risk (228 cases) and low-risk (228 cases) groups (Supplementary Table S2). Survival analysis showed that OS (Figures 3A–C), disease-free survival (DFS), progression-free survival (PFS), and disease-specific survival (DSS) (Supplementary Figure S6B) in the high-risk group were significantly worse than those of the low-risk group (Log-rank test, p < 0.001). The expression profiles of 14 genes were visualized as a heatmap (Figure 3D). The area under the ROC was 0.773 (Figure 3E). Then, the variables of age, stage, T, M, N, and PRS were analyzed to establish a nomogram for predicting the 3-years survival rate (Figure 3F).

Further, the PRS model was invalidated in the test set. The above PRS formula, cut-off value, and grouping method were used to divide patients into high-risk (248 cases) and low-risk (221 cases) groups (Supplementary Table S3). The OS of the high-risk group was significantly worse than that of the low-risk group (Log-rank, p < 0.001) (Figures 4A–C). The expression profiles of the PRS model genes were visualized as heatmaps (Figure 4D). The area under the ROC was 0.707 (Figure 4E). Then, the variables of age, stage, and risk score were analyzed to establish a nomogram for predicting the 3-years survival rate (Figure 4F).

After obtaining the reliable PRS model, we analyzed clinical and molecular features of PRS subgroups in the training set. The PRS subgroups showed significant differences in sex, T and N classifications, and stage. The proportion of female, T1, N0, and Stage I patients in the low-risk group was significantly higher than those in the high-risk group (Chi-square test, p < 0.05) (Figure 5A, Supplementary Table S4). GSEA analysis on the enrichment pathways of KEGG between the two PRS subgroups showed that pathways of cell cycle, DNA replication, homologous recombination, mismatch repair, p53 signal pathway, which were associated with gene mutation and chromosome instability, were significantly upregulated in the high-risk patients. In contrast, ABC transporters, B cell receptor signaling pathway, cell adhesion molecules (CAMs), histidine metabolism, and mTOR signaling pathway and other immune-related pathways were significantly upregulated in the low-risk patients (Figure 5B). In addition, PRSs of all patients were positively correlated with TMBs (Spearman correlation, p < 0.001), and negatively correlated with IPSs (Spearman correlation, p < 0.05) (Figures 5C,D) Meanwhile, in combination with literature data, we compared two prediction indicators of neoantigens (Figure 5E): the counts of mutations predicted to yield HLA-binding neopeptides (Predicted NeoAgs) and the ratios of observed versus expected binders per non-silent mutation (Observed/Expected NeoAgs) (Charoentong et al., 2017; Hakimi et al., 2016). The Observed/Expected NeoAgs of the low-risk patients was higher than that of the high-risk ones (Wilcoxon test, p < 0.05). Patients in the low-risk group may have more effective neoantigens to promote immunity against tumor and obtain more benefits from immunotherapy.

Then, we also analyzed mutation rates of genes between PRS subgroups in the training set. The integral mutation rate of the high- and low risk group was 92.61 and 78.51%, respectively. Then, mutation frequencies of 42 driver genes associated with LUAD were calculated. Mutation ratios of the driver genes TP53, LRP1B, CLIP1, EZH2, LRIG3, PIK3CA, RBM10, and KRAS in the high-risk group were higher than those in the low-risk group (Chi-square test, p < 0.05) (Supplementary Table S5). These results provide new insights into mutation biomarkers and immunotherapeutic targets. To understand the effects of these mutations on LUAD development, we conducted a cluster analysis and identified 11 signatures of de novo mutation (S1-S11) (Supplementary Figures S7A,B). Among them, S2, S3, S4, and S5 were similar to the curated signatures in COSMIC (Supplementary Figure S7C; Supplementary Table S6). S2 was related to DNA mismatch repair deficiency (dMMR), S3 was related to APOBEC (apolipoprotein B mRNA editing enzyme catalytic polypeptide-like), S4 was related to age, and S5 was related to tobacco mutagens. The contributions of S2, S3, S4, S5, S10, and S11 in the high-risk group were higher than those in the low-risk group (Cochran-Armitage trend test, p < 0.05) (Supplementary Figure S7D).

In the subsequent analysis, we examined the ability of the PRS model to predict the response to immunotherapy in Asian patients. A total of 27 NSCLC patients (GSE135222) who received anti-PD-1/PD-L1 immunotherapy were selected. On the basis of the cut-off values in this study, the patients were divided into high-risk and low-risk groups. Survival analysis showed that low-risk patients had better OS than high-risk ones (Log-rank test, p < 0.05, Supplementary Figure S8A). Another dataset of NSCLC patients (GSE126044) who were also treated by anti-PD-1/PD-L1 immunotherapy showed that the PRSs of those who responded to immunotherapy were significantly lower than those who did not respond (Wilcoxon test, p < 0.05, Supplementary Figure S8B). The results of these two datasets indicated that the PRS model was also a good predictor for the efficacy of immunotherapy.

Finally, the clinical (Supplementary Table S7) and mRNA sequencing data (Supplementary Table S8) of NJDT patients were analyzed by the PRS model and IPS algorithm (Supplementary Material S1). The PRSs of tumor samples were significantly higher than those of normal samples (Wilcoxon test, p < 0.01) (Supplementary Figure S9A); Among the tumor samples, 16 and 26 cases were categorized into the high- and low-risk groups, respectively. When the clinicopathological features were compared, the immune scores (calculated by ESTIMATE) and IPSs of the high-risk group were lower (Figure 6A), and EGFR mutations of high-risk patients were more frequent. In addition, patients with EGFR mutations had higher PRSs and lower IPSs than wild-type (WT) patients (Wilcoxon test, p < 0.05) (Figure 6B). But KRAS mutations did not showed the similar phenomenon. These results may suggest that the immunity state against tumor of WT patients was superior to that of mutant patients. Considering the mutation results in the training set, we compared the expression of MMR and APOBEC proteins. PMS1 was upregulated in the high-risk group, while APOBEC3A, C, D and G were upregulated in the low-risk group (Wilcoxon test, p < 0.05) (Figure 6C, Supplementary Table S9). Most checkpoints of patients in the low-risk group were higher than those in the high-risk group (Supplementary Figure S9B), and the immune-related pathways in the low-risk group were also upregulated (Supplementary Figure S9C). Furthermore, although the difference between groups was not statistically significant, the low-risk patients had better PFS than high-risk patients (Figure 6D).