Our study design is illustrated in Figure 1. The gene expression data and clinicopathological parameters of lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC) were downloaded from The Cancer Genome Atlas (TCGA) database (https://gdc-portal.nci.nih.gov/), including 994 cancerous and 107 normal samples with prognostic information. Additionally, gene expression patterns of lung cancer with prognostic information were also obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), namely, GSE50081 (16), GSE37745 (17), and GSE248378 (18). Meanwhile, GSE186143 (19), an independent cohort of malignant melanoma receiving ICI treatment, was also downloaded for the subsequent study. All data were assessed using the Illumina HiSeq 2000 RNA Sequencing platform.

Gene expression data for LUAD (n = 585) and LUSC (n = 550) were downloaded from TCGA, generated on the Illumina HiSeq 2000 platform as standardized log₂(FPKM + 1)values. After integrating clinical information, a total of 994 NSCLC tumor samples and 107 normal control samples with complete survival data were included in the final analysis (Supplementary Tables 1, 2). To leverage transcriptomic data reflecting immune dysregulation in irAEs—a strategy supported by prior studies—we prioritized TCGA-LC for analysis. This public dataset, with comprehensive lung cancer transcriptomes and clinical annotations, allowed cost-effective validation of gene expression in large cohorts. Firstly, we curated an initial set of irAE-related genes from the literature (20, 21) and applied objective statistical criteria [FDR < 0.05, p < 1 × 10⁻⁴ and RS (Spearman’s correlation coefficients) ≥0.7] to select the 30 most significant immunologic gene sets for further analysis (Supplementary Table 3), which represent key immune pathways in irAE pathogenesis. Subsequently, using GSVA (R v1.48.0), we assessed enrichment scores (ESs) for these gene sets in tumors vs. normal lung tissues to quantify their tumor microenvironment activity. Then, the optimal ES cutoff was determined via the X-tile software (Yale University, version 3.6.1) by maximizing the separation of high-/low-ES samples in TCGA-LC, as defined by the cutoff that generated the maximum log-rank χ² statistic. This cutoff was further validated by comparing overall survival between high- and low-ES groups. Next, differentially expressed genes (DEGs) were identified using Student’s t-test on the log₂(FPKM + 1) normalized expression values, with adjustment for multiple testing using the Benjamini–Hochberg false discovery rate (FDR) procedure, with significance thresholds set at FDR <0.05 and |log₂FC| >1. This analysis was performed between high- and low-ES specimens, as well as between cancer and non-cancer samples. By intersection analysis of DEGs between the above different groups, the correlations between the overlapping DEGs and the 30 genes related to irAEs were analyzed by applying the “cor” function in R. The genes with threshold values of |Pearson R| >0.7 and p <0.05 were identified as irAE-related DEGs. Then, to explore the biological significance and pathway enrichment, Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) (22).

Based on TCGA-LC samples, we first identified prognosis-related genes from irAE-related DEGs using univariate Cox regression (R package survival) with a significance threshold of p <0.05. Next, we applied the least absolute shrinkage and selection operator (LASSO) regression (R package lars) to derive a parsimonious subset of prognostic genes, selecting the optimal regularization parameter (λ) via 10-fold cross-validation under the minimum-criteria rule. These candidates were then entered into a multivariate Cox regression to construct a multigene predictive model. The irAE-related prognostic risk score was calculated as follows:

RS=Coef_Gene1×Gene1 expression+Coef_Gene2×Gene2 expression+…+Coef_{Gene n} ×Gene n expression.

Where, Coef_{Gene n} is the regression coefficient of Gene n.

Using this expression, the RSs of samples in the TCGA lung cancer dataset and the GEO validation dataset were analyzed. According to the median of RSs, all specimens were divided into high- and low-risk groups in the training and validation datasets of lung cancer, and survival analysis was conducted to assess the prognosis between different risk groups via the “survival” package (p <0.05). Additionally, the differences of RS values between irAE and non-irAE groups in the GSE186143 dataset were analyzed via the Kruskal–Wallis rank sum test.

Based on the expression levels in each sample of the TCGA-LC cohort, the proportion of immune cells in the tumor microenvironment was evaluated via the CIBERSORT algorithm, and the differences in immune cells (DICs) between high- and low-risk subtypes were analyzed via the “preprocessCore,” “e1071,” and “limma” packages. Additionally, the estimate package in R was applied to estimate the immune and stromal components via the ESIMATE scores, and the differences in stromal and immune scores between different risk groups were estimated using the “estimate” and “limma” packages. Furthermore, between the two risk groups, Hallmark signaling pathways were identified using gene set enrichment analysis (GSEA), which were analyzed by the “org.Hs.eg.db,” “clusterProfiler,” and “enrichplot” packages. Moreover, the expression of human leukocyte antigen (HLA) genes was analyzed via the “limma” package, with a threshold of p <0.05. The tumor immune dysfunction and exclusion (TIDE) scores of NSCLC samples were calculated utilizing the Tumor Immune Dysfunction and Exclusion (TIDE, http://tide.dfci.harvard.edu/login/) (23).

To study the independent prognostic factors, RSs and clinicopathological parameters in the TCGA-LC dataset were subjected to univariate and multivariate Cox regression analysis through the “survival” package, with a threshold value of p <0.05. Then, a nomogram model was constructed based on the independent parameters via the “rms” package. Concordance index (C-index) and area under the curve (AUC) for 1-, 3-, and 5-year survival probabilities were analyzed to evaluate the predictive accuracy of the nomogram model. Furthermore, to specifically assess the model’s stability and correct for overoptimism, we performed an internal validation using bootstrap resampling with 1,000 repetitions. The optimism-corrected C-index was calculated to evaluate the calibrated predictive performance.

A total of 41 patients with histologically confirmed NSCLC who were treated with ICIs were included from January 2020 to April 2023 at the 920^th Hospital (Yunnan, China). All data, including patient history, laboratory results, and irAEs, were obtained retrospectively. In our study, ICIs for NSCLC included sintilimab, camrelizumab, and tislelizumab, one of which was administered at a fixed dose every 3 weeks, as there was no evidence which ICI is superior to the others. All the patients were followed up for more than 6 months, and computed tomography (CT) was performed to evaluate the patients’ therapeutic response according to the Response Evaluation Criteria in Solid Tumors (RECIST) V1.1 (24). Furthermore, their irAEs were evaluated at every visit according to the Common Terminology Criteria for Adverse Events (CTCAE) version 5.0 (25). Two independent oncologists jointly determined irAEs based on the patients’ clinical symptoms, laboratory tests (e.g., liver function indicators, thyroid function), and imaging reports. In cases of disagreement in grading, a third expert was consulted for arbitration. irAEs were classified into grades 1–5 based on severity, with events of grade ≥3 defined as severe toxic reactions. The Institutional Ethics Committee at the 920^th Hospital approved this research, and informed consent based on the Declaration of Helsinki was obtained from all patients.

The specimens from the patients were fixed in 10% formalin, dehydrated with a series of acetone, and embedded in paraffin. The specimens were sliced into 4-μm slices and analyzed by hematoxylin and eosin staining. Immunohistochemistry (IHC) was performed to assess the protein expression of CRTAM, ANKRD44, and POU2AF1 on formalin-fixed, paraffin-embedded (FFPE) tissue sections. The following primary antibodies were used: rabbit monoclonal anti-CRTAM (Cat. ab820487, Zen bioscience, Chengdu, China; dilution 1:100), rabbit polyclonal anti-ANKRD44 (Cat. AP50195PU-N, OriGene, Beijing, China; dilution 1:200), and rabbit recombinant monoclonal anti-POU2AF1 (Cat. ab76010, Hua’an Biotechnology, Hangzhou, China; dilution 1:200). Antigen retrieval was performed in citrate buffer (pH 6.0) at 95°C for 20 min. Sections were then incubated with primary antibodies overnight at 4°C. After washing, immunostaining was visualized using a DAB detection system and counterstained with hematoxylin. Two senior pathologists reviewed the slides independently. The intensity scoring criteria were as follows: 0, no staining; 1, <10% staining; 2, 10% to 50% staining; and 3, 50% to 100% staining. The final IHC score was the sum of both scores: 0–2 for negative and 3–6 for positive.

All data processing was performed using R version 3.6.1. The normality of continuous variables was assessed using the Shapiro–Wilk test. Based on the results of the normality test and the exploratory nature of the data, non-parametric methods were primarily employed for group comparisons. Continuous variables were analyzed using the Kruskal–Wallis test and the Wilcoxon rank-sum test, while categorical variables were compared using the chi-square (χ²) test. Spearman’s correlation test was applied to determine correlation coefficients. Independent predictive factors were identified through Cox proportional hazards regression analysis. Statistical significance was defined as p <0.05.

The screening strategy of irAE-related DEGs in our study is shown in Figure 2A. GSVA was conducted to assess the ESs of 30 immunologic gene sets reported in the literature based on the TCGA-LC dataset. Normal lung tissues exhibited significantly higher ES values than NSCLC samples (Supplementary Figure S1A, p < 0.05), confirming immune pathway suppression in tumors. According to the optimal cutoff value of 0.23, the lung cancer samples were divided into high-ES groups (n = 761) and low-ES groups (n = 233). Notably, the OS of the high-ES group was significantly longer than that of the low-ES group (p < 0.05) (Supplementary Figure S1B). Then, 1,057 DEGs between the high- and low-ES groups were identified. After intersection with 3,263 DEGs between cancer and normal groups, 565 overlapping DEGs were obtained (Figures 2B, C; Supplementary Table 4). By conducting further Pearson correlation analysis, 132 irAE-related DEGs were identified (Figure 2D), with a total of 607 pairs of connections (|Pearson R| > 0.7 and p < 0.05, Supplementary Table 5).

To elucidate the biological significance and pathway enrichment underlying the irAE-related DEG set, GO and KEGG pathway analyses were conducted, respectively. As shown by GO analysis, 50 biological processes, including cytokine activation, cell adhesion, and receptor recognition, were significantly enriched. Meanwhile, some immunologic regulatory pathways, such as immune response, inflammatory response, and positive regulation of T-cell proliferation, were identified via KEGG analysis (Figures 2F, G; Supplementary Table 6).

Moreover, according to the prognostic information in the TCGA-LC dataset, 14 prognosis-related genes were identified from 132 irAE-related DEGs. To avoid overfitting, LASSO regression analysis was further performed. A panel of eight feature genes was selected and used to generate the RS model as follows (Figure 2E):

RS = (−0.0443) × ANKRD44 + (−0.0134) × POU2AF1 + (−0.0056) × P2RY13 + (−0.0056) × CXCL13 + (−0.0036) × CD69 + (−0.0025) × HLA-DMB + (−0.0002) × ARHGAP25 + (−0.0160) × CRTAM

Furthermore, the TCGA-LC training dataset was clustered into different subtypes based on the medium expression of each feature gene. Kaplan–Meier (KM) survival analysis revealed that there was a significant difference in prognosis between the high and low subtypes (Figure 2H, p < 0.05).

The RS of each sample was analyzed in the TCGA lung cancer dataset. According to the median level of RS, patients were stratified into high- and low-risk groups. Compared to the TCGA lung cancer cases in the low-risk group, the OS rate in the high-risk group was shorter (Figure 3A). ROC analysis showed that the predictive value (AUC) of this RS for the 5-year survival rates of lung cancer patients in the TCGA dataset was 0.770 (Figure 3D). We also observed that the proportion of deaths in the high-risk group was higher than in the low-risk group (Figures 3B, C). Furthermore, the prognostic value of the RS was also validated in the GSE50081 and GSE37745 datasets. In these two cohorts, the OS rates of patients with high RSs were significantly lower than those with low RSs (Figures 3E, I), and the AUC values for predicting the 5-year survival rate were 0.767 and 0.758, respectively (Figures 3H, L). Meanwhile, a higher proportion of deaths was also identified in the high-risk group (Figures 3F, G, J, K).

To investigate the ability of RS to predict immunotherapy efficacy, RS was applied to the GSE248378 cohort, in which all lung cancer patients had complete clinical information of ICI treatment. Firstly, the expression of eight signature genes of RS (ANKRD44, POU2AF1, P2RY13, CXCL13, CD69, HLA-DMB, ARHGAP25, and CRTAM) was investigated between the two RS groups. In the low-RS group, eight feature genes were upregulated (Figures 4A–H). When compared to the high-risk group, the patients in the low-risk group had significantly better OS (Figure 4I). Time-dependent ROC analysis showed that the AUC values at 1, 3, and 5 years for predicting OS associated with ICI treatment were 0.852, 0.784, and 0.763, respectively (Figures 4J). The proportion of patients with no benefits was lower in the low-risk group than in the high-risk group (Figure 4K). Moreover, by TIDE analysis, we found that the patients with low RSs had lower TIDE scores than those with high RSs (Figure 4L). These results indicated that patients in the high-risk group exhibited resistance to ICIs.

Considering that the RS was conducted based on irAE-related genes, we evaluated the predictive value of irAEs. In the melanoma cohort GSE186143 with clinical information of irAEs, we found that the incidence of irAEs in low-risk group was significantly higher than in the high-risk group, and patients with irAEs had significantly lower RS than patients without irAEs (Figures 5A, B). Further analysis revealed that the AUC value of the RS to predict irAEs was 0.807 (Figure 5C). In addition, we investigated the potential of irAE prediction in our cohort with a total of 41 NSCLC patients with ICI treatment (Supplementary Table 7). Among them, 16 patients (39.02%) experienced at least one adverse event, and 25 patients (60.98%) had no irAEs. The most common irAEs were thyroiditis (n = 5, 31.25%), followed by reactive cutaneous capillary endothelial proliferation (RCCEP) (n = 4, 25%) and immune-related pneumonia (n = 2, 12.5%). Among the patients with irAEs, 4 patients (25%) presented with a complete response (CR), 10 patients (62.5%) developed a partial response (PR), 1 patient (6.25%) had stable disease (SD), and 1 patient (6.25%) had progressive disease (PD). Meanwhile, in the non-irAE group, the number of patients with CR, PR, SD, and PD was 2, 9, 8, and 6, respectively (Supplementary Table 8). The objective response rate (ORR, CR+PR) and the disease control rate (DCR, CR+PR+SD) were higher in the irAE group than in the non-irAE group (87.50% vs. 44.00%, 93.75% vs. 76.00%) (Figure 5D). Furthermore, patients with irAEs had a significantly longer PFS and OS than patients without irAEs (median PFS: 20.00 vs. 9.86 months; median OS: 21.37 vs. 12.08 months) (Figures 5E, F). Next, we evaluated the expressions of three risk genes with the highest weights (POU2AF1, ANKRD44, and CRTAM) in patients with irAEs versus those without irAEs. As compared to patients without irAEs, the protein level and positive rate of CRTAM, POU2AF1, and ANKRD44 were increased in patients with irAEs (Figures 5G, H). We further evaluated the correlation between survival period and protein expression and found that the OS in the ANKRD44-positive group was significantly improved than in the negative group (Figure 5I). The expression of POU2AF1, CRTAM, and ANKRD44 was significantly upregulated in lung tumor tissues with irAEs (10.5-fold, 1.97-fold, and 2.62-fold vs. no irAEs, p < 0.05) and showed a positive correlation with irAE occurrence (Supplementary Table 9). Immunohistochemical staining (10 × 20 magnification) further validated nuclear localization differences in the three risk gene expression patterns (Figure 5J).

Because dysregulation of immune cell infiltration in the TME plays a crucial role in immunotherapy, we investigated the correlation between immune cell infiltration and RS based on the CIBERSORT algorithms. The proportion of 22 types of immune cells was quantified on the expression data in TCGA-LC samples. As compared with the high-risk group, the abundance of CD8⁺ T cells, CD4⁺ memory resting T cells, M1 macrophages, and activated mast cells was obviously elevated in the low-risk group (Figure 6A). Based on the estimate algorithm, we found that the stromal, immune, and ESTIMATE scores were significantly higher in the low-risk group, along with a dramatic decrease in tumor purity (Figure 6B). We further used GSEA to search for Hallmark terms across the TCGA-LC dataset and found that the significantly enriched Hallmark terms in the low-risk group were cytokine–cytokine receptor interaction, cell adhesion molecules, etc. (Figures 6D, E). Meanwhile, the differences of human leukocyte antigen (HLA) genes were assessed between the two risk groups, which were essential in antigen presentation. As compared with the high-risk group, most HLA genes, such as HLA-I and HLA-II, had higher expression levels in the low-risk group (Figure 6C). Our results indicated that risk score was closely related to TME in lung cancer patients, in which the low-risk group tended to have more immune cells with enhanced antitumor immune responses.

To explore the independent prognostic factors for patients with lung cancer, we combined clinical features and the precalculated score value into an integrated analysis. As shown by the univariate Cox regression analysis, some factors, including age, pathological T, clinical stage, tumor recurrence, enrichment score of irAEs, and RS, were significantly related to OS (p < 0.05) (Figure 7A). The multivariate Cox regression analysis further revealed that age, clinical stage, tumor recurrence, and RS were independent prognostic factors (p < 0.05) (Figure 7B). The detailed statistical results are shown in Supplementary Table 10. The effect of age, clinical stage, and tumor recurrence on survival was shown by KM plot analysis (Figure 7C). To improve the clinical unity of RS, a nomogram was established based on the four independent factors (Figure 7D). The C-index values of our model at 1, 3, and 5 years were 0.796, 0.784, and 0.763, respectively (Figure 7E). As shown by the decision curve analysis (Figure 7F), the nomogram combined with four factors had a better clinical benefit net than that combined with a single factor. In addition, Figure 7G reflects the clinical benefit of the predictive model at different decision thresholds, indicating that the model has significant clinical value in identifying high-risk patients and guiding early intervention. Finally, bootstrap internal validation (n = 1,000) confirmed the model’s robustness. Visual comparison of the corrected calibration curves (Supplementary Figure S2) showed that our RS model’s performance surpassed all alternative predictors, highlighting its superior and stable predictive ability.