Bioinformatics analysis of an immunotherapy responsiveness-related gene signature in predicting lung adenocarcinoma prognosis

Introduction

Non-small cell lung cancer (NSCLC) is the major type of lung cancer and can be further subdivided into four pathological types: lung squamous cell carcinoma (LUSC), lung adenocarcinoma (LUAD), large cell carcinoma (LCC), and pulmonary sarcomatoid carcinoma (PSC). LUAD accounts for 55–60% of all NSCLC cases (1). At present, the survival rate of NSCLC remains unsatisfactory, making it the leading cause of cancer-related death in both men and women (2,3). Recently, the development of immunotherapy has provided substantial improvements to the overall survival (OS) and progression-free survival (PFS) rates of patients with lung cancer (4).

The basic therapeutic options for NSCLC currently include surgery, radiotherapy, chemotherapy, targeted therapy, and immunotherapy (3-5). Immunotherapy is based on the notion that cancer cells can escape tumor immunosurveillance via antigen modulation or mutation; thus, reactivation of the immune system defense with immune checkpoint inhibitors (ICIs) represents a novel pathway to treating cancer. Immunotherapy has served as a major front-line treatment for cancer therapy over the last decade. However, some patients respond poorly to immunotherapy, especially among patients with LUAD, possibly due to the differential composition of the tumor immune microenvironment (TIME) (6,7). Therefore, the response to immunotherapy may provide novel insights into the prognosis of LUAD.

Many efforts have been made to predict the responsiveness to immunotherapy in LUAD patients. Currently, it is known that lack of immune cell infiltration in the tumor (“immune-cold” tumor) or a suppressive TIME are correlated to poor responses to immunotherapy (8,9). Multiple models have been reported to be able to predict the prognosis of LUAD, but most have been based on certain pathogenesis pathways rather than on responsiveness to therapies, thereby limiting their potential in predicting the efficacy of treatments (10-12). With the advancement in next-generation sequencing (NGS) and bioinformatics, we are now able to obtain clinical prognosis cohorts and their tumor sequencing information from open databases. The Cancer Genome Atlas (TCGA)-LUAD database and GSE72094 database are two most widely used sources of LUAD, they are highly rated for their completeness of clinical data, standardization in formats, and massive sample sizes (13,14). In this study, we aimed to establish a prognostic prediction model for LUAD based on the differentially expressed genes (DEGs) between responders and nonresponders to immunotherapy [anti-programmed cell death protein 1/anti-programmed death-ligand 1 (anti-PD-1/anti-PD-L1)] with data obtained from TCGA and Gene Expression Omnibus (GEO) databases. Our model could not only precisely predict the overall prognoses in patients with LUAD but also provide satisfactory predictions in subgroup analyses. Mechanistically, we found that the most significantly varied pathway between each group was ribonucleoprotein complex biogenesis, which correlated closely with TP53 and TTN mutation probabilities, suggesting that this pathway is the kernel mechanism determining the prognosis of LUAD. Moreover, the immune microenvironments and drug sensitivities varied significantly between each group. Additionally, a nomogram was also established to facilitate the application of the model in clinical practice. We present this article in accordance with the TRIPOD reporting checklist (available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-24-309/rc).

Methods

Data acquisition

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The TCGA-LUAD dataset was downloaded from TCGA database (https://portal.gdc.cancer.gov), which contained the transcriptome data of 598 patients with LUAD, including 539 tumor samples and 59 normal samples. Besides the transcriptome data, clinical information such as age, sex, clinical stages, T-stage, N-stage, M-stage, and follow-up status, were included in this cohort. After the removal of normal sample data and data without survival status, data from 508 tumor samples were included. Other data sources included the GSE126044, GSE72094, and GSE135222 datasets obtained from the GEO database, and data without survival status were also removed. The GSE72094 dataset contained gene expression array information of 442 LUAD patients along with the clinical information including age, sex, and follow-up status. The detailed clinical features of the patients in the GSE126044 and GSE72094 cohorts are displayed in Table 1 (clinical feature information was not available for the GSE135222 cohort). The study design is illustrated in Figure 1.

Figure 1 Design and flowchart of the study. GEO, Gene Expression Omnibus; LUAD, lung adenocarcinoma; TCGA, The Cancer Genome Atlas; DEG, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes.

Screening for prognosis-related genes

In the GSE126044 cohort, the R package “limma” (The R Foundation for Statistical Computing) was used to identify 1,636 DEGs between the responders and nonresponders to immunotherapy, with the filter criteria being P<0.05 and |log₂ fold change (FC)|>1. TCGA cohort and GSE72094 cohort were merged using the R package “sva”, and then the batch effects between the two cohorts were eliminated with “combat” algorithm. A total of 1,514 intersecting genes were obtained after the DEGs and merged cohort sets were mixed. Based on the clinical information, survival analysis (Kaplan-Meier analysis, P<0.05) of the 1,514 intersection genes was conducted in TCGA cohort, and 1,163 prognostic-related genes were identified.

Establishment and verification of the prognostic model

Univariate Cox regression analysis was applied using TCGA cohort as the training set and the GSE72094 cohort as the validation set. The 1,163 prognosis-related genes were first screened with P<0.05 as the filter criterion, which was followed by an analysis with least absolute shrinkage and selection operator (LASSO) Cox regression to subsequently build the risk prediction model. The R package “survminer” was used to calculate the optimal cutoff value of the risk score, and samples were divided into two groups according to the optimal cutoff value: a high immunotherapy responsiveness-related risk score (IRRS) group and a low IRRS group. The survival curves of the high- and low-IRRS groups were generated using the “survminer” and “survival” R packages. Risk score stability was analyzed in the validation set, and the R package “survivalROC” was used to evaluate the performance of the prognostic formula using time-dependent receiver operating characteristic (ROC) analysis. Univariate and multivariate Cox regression analyses were used to evaluate whether the risk score was an independent prognostic predictor for OS.

Enriched pathway analysis

The R package “clusterProfiler” was used to identify enriched pathways in Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses.

TIME and treatment effect analysis

Immune cell infiltration was calculated using the CIBERSORT algorithm. The relationship between the model and the infiltration of immune cells was analyzed using risk score and proportional immune cell infiltrations via the R package “psych”. Moreover, immune checkpoints were acquired as previously described (11,12). The immunophenoscore (IPS) were obtained from the Cancer Immunome Atlas (https://tcia.at/). Finally, whether the formula could predict survival in patients receiving immunotherapy was tested in the GSE135222 cohort (15,16).

Somatic mutation analysis

Genome mutation information of TCGA-LUAD was obtained from TCGA database. The R package “maftools” was used to present the mutational differences between the low- and high-IRRS groups.

Drug sensitivity tests analysis

For anticancer drugs, the half maximal inhibitory concentration (IC₅₀), which can indicate the efficacy of a substance in prohibiting an organism or a biochemical process, was predicted in the low-IRRS group and high-IRRS group with the R package “pRRophetic”.

Establishment of the nomogram

Clinical information, including age, sex, and clinical stages of TCGA-LUAD patients, was extracted. To evaluate whether the risk score was an independent prognostic factor for OS, the clinical information and risk scores were analyzed using univariate and multivariate Cox regression. A Cox regression model along with the R package “rms” were used to establish an OS prediction nomogram incorporating prognostic factors and risk score. Survival probabilities at 1, 3, and 5 years were set as endpoints. To assess the performance of the nomogram, calibration plots were constructed to visualize the consistencies between the actual and predicted survival probabilities at 1, 3, and 5 years.

Statistical analysis

R software (version 4.1.0) was used to perform all statistical analyses. The Wilcoxon signed-rank test was used for TIME analysis, while the log-rank test was used for intergroup Kaplan-Meier survival analysis. Differences were considered statistically significant at P<0.05.

Results

Screening and functional analysis of immunotherapy responsiveness-related DEGs

First, we aimed to identify DEGs between immunotherapy responders and immunotherapy nonresponders. A total of 1,635 DEGs were obtained, 1,161 of which were upregulated (log₂FC>1 and P<0.05) in immunotherapy responders compared to nonresponders and 475 of which were downregulated (log₂FC>1, P<0.05) (Figure 2A,2B). In exploring the functions of DEGs using GO enrichment analysis, we first identified that the immunotherapy responsiveness-related DEGs with the greatest magnitude were involved in neutrophil activation (Figure 2C). In KEGG enrichment analysis, we found that immunotherapy responsiveness-related DEGs were substantially associated with cell adhesion molecules (Figure 2D).

Figure 2 Screening and functional analysis of immunotherapy responsiveness-related DEGs. (A) Heatmap of differentially expressed immunotherapy responsiveness-related genes. (B) The volcano plot of DEGs with a cutoff at P<0.05 and |log2FC| >1. (C) GO enrichment of immunotherapy responsiveness-related DEGs. (D) KEGG pathway of immunotherapy responsiveness-related DEGs. DEG, differentially expressed gene; FC, fold change; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MHC, major histocompatibility complex; BP, biological process; CC, cellular component; MF, molecular function.

Establishment of the risk score formula incorporating immunotherapy responsiveness-related DEGs

A Venn diagram was drawn showing the 1,514 coexpressed immunotherapy responsiveness-related DEGs between GSE126044, TCGA-LUAD, and GSE72094 (Figure 3A). To assess the prognostic value of these immunotherapy responsiveness-related DEGs, we first performed a univariate Cox regression analysis plus LASSO regression analysis (Figure 3B,3C). A total of 34 genes (Figure S1) were selected and incorporated into the following risk score formula: IRRS = LRRC37A3 × 0.376 + PABPC3 × 0.327 + TRIM28 × 0.318 + CEBPB × 0.309 + NGRN × 0.269 + SQLE × 0.215 + EFNB2 × 0.188 + HPSE × 0.134 + PTPRH × 0.137 + KRT18 × 0.121 + INTS7 × 0.112 + GPR37 × 0.108 + LGR4 × 0.106 + RAB27B × 0.086 + KCNK1 × 0.080 + DDX56 × 0.079 + ANKRD65 × 0.062 + HS3ST2 × 0.018 + TFAP2A × 0.017 + KLHDC9 × 0.009 − CASZ1 × 0.023 − TP53I3 × 0.025 − ZNF750 × 0.033 − STAP1 × 0.067 − IVD × 0.078 − HTRA − 4 × 0.104 − IQCC × 0.108 − CCT6A × 0.168 − FCRL6 × 0.171 − MTUS1 × 0.254 − ZNF77 × 0.348 − MYO6 × 0.354 − CAMSAP3 × 0.393 − CD68 × 0.415. According to the optimal cutoff value (Figure 3D-3F), 323 patients were categorized into the low-IRRS group while 185 patients were categorized into the high-IRRS group. Kaplan-Meier survival analysis indicated that patients in the low-IRRS group had significantly higher OS than did those in the high-IRRS group [hazard ratio (HR) =4.86, 95% confidence interval (CI): 3.59–6.59; P<0.0001] (Figure 3G). The area under the curve (AUC) for predicting 1-year, 3-year, and 5-year OS was 0.798, 0.796, and 0.773, respectively (Figure 3H). These results showed that the risk score formula had excellent performance in predicting the OS of patients with LUAD.

Figure 3 Construction of the IRRS formula. (A) The Venn diagram displays the intersection of common genes among the three cohorts. (B) LASSO coefficient profile plots of prognostic-related genes showing that the variations in the size of the coefficients of parameters decrease with an increasing value of the k penalty. (C) Penalty plot for the LASSO regression analysis. (D) Distribution patterns of survival time and status in the high-IRRS group and low-IRRS group in the training set. (E) Distribution of IRRS in the training set. (F) Heatmap of the 34 prognostic-related genes for each patient in the training set. (G) Kaplan-Meier survival curve of patients in the high-IRRS group and low-IRRS group in the training set. (H) Time-related ROC analysis to evaluate the prognostic value of the IRRS in the training set. DEG, differentially expressed gene; IRRS, immunotherapy responsiveness-related risk score; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; AUC, area under the curve; TPR, true positive rate; FPR, false positive rate.

Validation for the stability of the IRRS formula

Data from the validation set (GSE72094 cohort) were used to verify the stability of the IRRS formula established from TCGA cohort (training set). Patients in the GSE72094 cohort were also categorized into a high-IRRS group (149 patients) and low-IRRS group (249 patients) in accordance with the cutoff value developed in TCGA cohort (Figure 4A-4C). The OS of patients in the high-IRRS group was significantly lower than that in the low-IRRS group (HR =1.88, 95% CI: 1.30–2.72; P<0.001) (Figure 4D). The AUC for predicting 1-year, 3-year, and 5-year OS were 0.64, 0.627, and 0.66, respectively (Figure 4E). These results, in combination with the results presented above, indicate that the IRRS formula indeed had good performance in predicting OS for LUAD.

Figure 4 Validation of the stability of the IRRS formula. (A) Distribution patterns of survival time and status in the high-IRRS group and low-IRRS group in the validation set. (B) Distribution of IRRS in the validation set. (C) Heatmap of the 34 prognostic-related genes for each patient in the validation set. (D) Kaplan-Meier survival curve of patients in the high-IRRS group and low-IRRS group in the validation set. (E) Time-related ROC analysis for evaluating the prognostic value of the IRRS in the validation set. IRRS, immunotherapy responsiveness-related risk score; ROC, receiver operating characteristic; AUC, area under the curve; TPR, true positive rate; FPR, false positive rate.

Subgroup analyses to evaluate the performance of the IRRS formula

We further carried out Kaplan-Meier survival analyses in patients from TCGA cohort with different LUAD clinical stages, ages (>65 or ≤65 years), and sexes to test the stability of the IRRS in predicting prognoses in subgroups. The IRRS excellently predicted OS in patients with LUAD from stage IA to stage IV (Figure 5A-5G) (patients who lack the information for subdivision into IA/IB or IIA/IIB were excluded). Moreover, the IRRS accurately identified patients with discrepant prognoses in both age subgroups and sexes (Figure 5H-5K). These results suggest that the IRRS formula could not only predict prognoses in patients with LUAD but could also predict OS in the different LUAD subgroups (categorized by clinical stage, age, and sex).

Figure 5 Subgroup analyses with the IRRS formula. (A) Kaplan-Meier survival curve of patients in the high-IRRS group and low-IRRS group in patients with stage IA with LUAD. (B) Kaplan-Meier survival curve of patients in the high-IRRS group and low-IRRS group in patients with stage IB LUAD. (C) Kaplan-Meier survival curve of patients in the high-IRRS group and low-IRRS group in patients with stage IIA LUAD. (D) Kaplan-Meier survival curve of patients in the high-IRRS group and low-IRRS group in patients with stage IIB LUAD. (E) Kaplan-Meier survival curve of patients in the high-IRRS group and low-IRRS group in patients with stage IIIA LUAD. (F) Kaplan-Meier survival curve of patients in the high-IRRS group and low-IRRS group in patients with stage IIIB LUAD. (G) Kaplan-Meier survival curve of patients in the high-IRRS group and low-IRRS group in patients with stage IV LUAD. (H) Kaplan-Meier survival curve of patients with age >65 years in the high-IRRS group and low-IRRS group. (I) Kaplan-Meier survival curve of patients with stage IA LUAD and age ≤65 years in the high-IRRS group and low-IRRS group. (J) Kaplan-Meier survival curve of female patients in the high-IRRS group and low-IRRS group. (K) Kaplan-Meier survival curve of male patients in the high-IRRS group and low-IRRS group. IRRS, immunotherapy responsiveness-related risk score; LUAD, lung adenocarcinoma.

Identification of IRRS-associated mutation landscape and biological functions

Next, we examined whether the IRRS correlated with the tumor mutation landscape in LUAD. We compared the mutational landscapes of the two IRRS groups. In comparison to that of low-IRRS group, the overall tumor mutational burden was significantly higher in the high-IRRS group, with the TP53 and TTN being the most predominant mutations. The other 13 LUAD-related genes also had significantly higher mutation burdens in the high-IRRS group, whereas the TP53 and TTN genes had higher mutation burdens in the low-IRRS group (Figure 6A). These results indicated that IRRS was correlated with the LUAD mutation landscape. We further analyzed the differences in pathways between the low- and high-IRRS groups using GO enrichment analysis. The biological processes with the most significant differences were ribonucleoprotein complex biogenesis, antigen processing and presentation, and ribosome biogenesis (Figure 6B). KEGG enrichment analysis was also performed, and the most significantly different pathways were the cell cycle, focal adhesion, and nucleocytoplasmic transport (Figure 6C). These differences may explain the lower OS in the high-IRRS group.

Figure 6 Mutation landscape and functional analysis of IRRS-related DEGs. (A) Mutation burdens of various oncogenes in the high-IRRS group and low-IRRS group. The genes with significant differences between the two groups were highlighted in red font. (B) GO enrichment analysis of DEGs between the high-IRRS and low-IRRS groups. (C) KEGG pathway enrichment of DEGs between the high-IRRS and low-IRRS groups. TMB, tumor mutation burden; IRRS, immunotherapy responsiveness-related risk score; DEG, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes.

Comparison of immune features and immunotherapy responsiveness of the two IRRS groups

We subsequently explored the correlation between IRRS and immune features in TCGA cohort. The association between IRRS and the TIME was determined using the stromal score (SS), immune score (IS), and ESTIMATE score. The results revealed a negative correlation between the IRRS and immune status (Figure 7A). There were also significant differences in the proportional infiltration of immune cells between the high- and low-IRRS groups. Naïve B cells, activated natural killer (NK) cells, macrophages (M0, M1, and M2), resting dendritic cells, eosinophils, and neutrophils were enriched in the high-IRRS group, whereas CD4 memory T cells (activated), plasma cells and mast cells (activated) were enriched in the low-IRRS group (Figure 7B). The correlations between IRRS and immune cell types and between different types of immune cells were also characterized (Figure 7C). We found that the expression of immune checkpoints differed significantly between the two IRRS groups (Figure 7D). To explore the usefulness and application value of IRRS in predicting the prognoses of patients with LUAD receiving immunotherapy, Kaplan-Meier survival analysis was performed in the GSE135222 cohort. The results showed that the IRRS could predict the prognosis of patients with LUAD receiving immunotherapy (Figure 7E). Moreover, we found that the IRRS negatively correlated with the IPS, the IPS of cytotoxic T-lymphocyte-associated protein 4 (CTLA4) was significantly higher in the low-IRRS group, but the IPS score of PD-1 did not quite reach statistical significance between these groups (P=0.06) (Figure 7F,7G). Collectively, these results indicated that IRRS correlated with the immune status of LUAD and had prognostic value in patients with LUAD receiving immunotherapy. In addition, we also used GSE135222 to test whether the IRRS could evaluate the PFS of patients with LUAD treated with PD-1/PD-L1. However, we found that the PFS in the high-IRRS group was longer than that in the low-IRRS group, which may indicate that low IRRS may indicate malignancy in patients with LUAD.

Figure 7 Immune status analysis and anticancer drug sensitivity analyses. (A) The relationships between IRRS, immune score, stromal score, and ESTIMATE score. (B) Infiltration of various types of immune cells in the high-IRRS group and low-IRRS group. (C) The correlations between IRRS and immune cells. (D) The expressions of immune checkpoints in the high-IRRS group and low-IRRS group. (E) Kaplan-Meier survival curve of patients in the high-IRRS group and low-IRRS group in the GSE135222 cohort. (F,G) IPS of the high-IRRS group and low-IRRS group. (H) The IC₅₀ of cisplatin in the low-IRRS group and high-IRRS group. (I) The IC₅₀ of erlotinib in the low-IRRS group and high-IRRS group. (J) The IC₅₀ of gemcitabine in the low-IRRS group and high-IRRS group. (K) The IC₅₀ of vinorelbine in the low-IRRS group and high-IRRS group. (L) The IC₅₀ of paclitaxel in the low-IRRS group and high-IRRS group. *, P<0.05; **, P<0.01; ***, P<0.001. IRRS, immunotherapy responsiveness related risk score; IPS, immune cell proportion score; TME, tumor microenvironment; NK, natural killer cell; PD-1, programmed cell death protein 1; CTLA4, cytotoxic T-lymphocyte associated protein 4; IC₅₀, half maximal inhibitory concentration.

Anticancer drug sensitivity analysis in high-IRRS and low-IRRS groups

To further explore whether IRRS could predict anticancer drug sensitivity, we compared the IC₅₀ of different anticancer drugs in the high-IRRS and low-IRRS groups. The anticancer drugs tested covered almost all the common drugs recommended for LUAD treatment according to the National Comprehensive Cancer Network (NCCN; https://www.nccn.org). The IC₅₀ of cisplatin, erlotinib, gemcitabine, vinorelbine, and paclitaxel was higher in the low-IRRS group (Figure 7H-7L). These results suggest that the IRRS is also valuable in predicting anti-LUAD drug sensitivities.

Construction and validation of a nomogram incorporating the IRRS

We additionally performed univariate Cox regression analysis to verify whether the IRRS could independently predict the OS of patients with LUAD. The results showed that the IRRS, together with clinical stage, T stage, N stage, and M stage, were independent predictors of OS in patients with LUAD, while age and sex were not (Figure S2A). After other confounding factors were controlled for via multivariate regression analysis, the IRRS remained an independent predictor for OS of LUAD, while the other parameters, including clinical stage, T stage, N stage, and M stage, were not (Figure S2B). These results revealed that IRRS was an independent predictor for OS in patients with LUAD.

We further established a nomogram incorporating certain parameters, including age, sex, T stage, N stage, M stage, clinical stage, and IRRS. Each parameter was assigned a specific score (Figure 8A), and calibration plots verified that the calculated scores could accurately predict the prognoses of 1-year survival probability, 3-year survival probability, and 5-year survival probability in patients with LUAD, with an area under curve (AUC) of 0.718, 0.702, and 0.68, respectively (Figure 8B-8E). The establishment of a nomogram should facilitate the application of the IRRS in clinical practice.

Figure 8 Construction and validation of IRRS-related nomogram. (A) A nomogram incorporating parameters including age, sex, T stage, N stage, M stage, clinical stage, and the IRRS. (B) The calibration curves for 1-year OS prediction. (C) The calibration curves for 3-year OS prediction. (D) The calibration curves for 5-year OS prediction. (E) The ROC curve for the nomogram. IRRS, immunotherapy responsiveness related risk score; OS, overall survival; ROC, receiver operating characteristic curve; AUC, area under curve.

Discussion

LUAD, the most common form of lung cancer, is one of the leading causes of death worldwide, and its prognosis remains poor (17,18). Immunotherapy is a novel treatment option and has been approved for treating LUAD, exhibiting promising outcomes with both monotherapy and combination therapy regimens (19). However, a large proportion of patients respond poorly to immunotherapy, and currently the cost of this therapy hinders its accessibility; therefore, it is critical to identify novel means to predict prognoses and treatment responsiveness in patients with LUAD (20). Based on previous data on immunotherapy response, in the present study, we constructed and validated a prognostic signature incorporating immunotherapy responsiveness-related DEGs and established an IRRS formula. We found that the IRRS reliably predicted the 1-, 3-, and 5-year OS in patients with LUAD and was an independent predictor of OS for LUAD. Subgroup analyses further revealed that the IRRS had a good performance in predicting OS in different LUAD subgroups (categorized by clinical stage, age, and sex). We also validated the usefulness of the IRRS in predicting immunotherapy responsiveness in another LUAD cohort (GSE135222). Mechanistically, we found that the IRRS correlated well with the immune status and mutation landscape in LUAD, and the DEGs between the high-IRRS and low-IRRS groups were enriched mainly in biological processes including ribonucleoprotein complex biogenesis, antigen processing and presentation, and ribosome biogenesis and were enriched in pathways including cell cycle, focal adhesion, and nucleocytoplasmic transport. The tumor mutational burden was significantly higher in the high-IRRS group, among which the TP53 and TTN were the mutations with the most significant differences. Moreover, we found that IRRS could predict anticancer drug sensitivity in LUAD, and we further constructed a nomogram incorporating parameters including sex, age, T stage, N stage, M stage, clinical stage, and IRRS. These provide examples of the applications of IRRS in clinical practice.

In the prognosis prediction model, we incorporated 34 prognostic-related DEGs into the risk score formula, among which Leucine-rich repeat containing 37 member A3 (LRRC37A3) and CD68 had the greatest positive and negative coefficient values, respectively. LRRC37A3 is a newly identified paralog of the core duplicon LRRC37A located on chromosome 17, and its functions remain largely unknown; however, one study did report that single-nucleotide polymorphisms (SNPs) of LRRC37A3 are involved in growth retardation, hepatopathy, and intellectual disability (21). The functions of other paralogs of LRRC37A are also ambiguous although a few studies have suggested their participation in the pathogenesis of Parkinson disease (22), coronary heart disease (23), and antibody reaction (24). LRRC37A was found to be a membrane-associated protein implicated in inflammation, cellular migration, and chemotaxis (22). Further investigations are needed to clarify the role and mechanisms of LRRC37A3 in LUAD or in immunotherapy responsiveness. CD68 is a macrophage-myeloid cell-associated antigen, and macrophages are the predominant cells expressing CD68 in the TIME (25). CD68-positive tumor-associated macrophages (CD68⁺ TAMs) are one of the hallmarks of TIME in LUAD (26), indicating that our model reflects differences in TIME between immunotherapy responders and nonresponders. Our model also revealed that the biological process with the most notable difference between the high- and low-IRRS groups was ribonucleoprotein complex biogenesis. Ribonucleoprotein complexes are RNA-binding proteins that are conjugated with RNA and exert broad spectral functions inside the cells, including tumorigenesis (27) and anticancer drug sensitivity (28,29). The IRRS established in the present study could also predict the anticancer drug sensitivities of LUAD. We found that patients in the high-IRRS group had a higher IC₅₀ for cisplatin, erlotinib, gemcitabine, vinorelbine, and paclitaxel, which have been proven to be effective and approved for treating LUAD (3). This indicates that the IRRS could be used to guide clinicians in choosing anticancer drugs for patients with LUAD.

The exploration of predictive models and the underlying mechanisms of immunotherapy responsiveness or resistance in LUAD is currently an area of intense research activity. Our recent study suggested that a Golgi apparatus-related signature could satisfactorily predict prognosis and immunotherapy responsiveness in LUAD (11), and recent studies have revealed that other signatures, including exosomes, lipid metabolism and immunity, basement membranes, one-carbon metabolism, tumor microenvironment, and epidermal growth factor receptor inhibitor resistance also correlate with immunotherapy efficacy (30-32). To our knowledge, this study is the first to predict the OS and immunotherapy efficacy in patients with LUAD based on immunotherapy responsiveness-related DEGs. We found that IRRS, the signature we established in this study, correlated with the TIME status in LUAD. The stromal, immune, and ESTIMATE scores all revealed significantly lower tumor micro environment (TME) scores in the high-IRRS group, which indicated lowered immune response in this group, and thus it is currently clear that a lower immune response correlates with poorer immunotherapy response in LUAD (33-35). We also observed remarkable differences in the infiltration proportions of various types of immune cells. For example, our model revealed that there was significantly higher infiltration of neutrophils in the high-IRRS group, and previous studies have demonstrated that the accumulation of neutrophils could lead to immunotherapy resistance (36,37). In addition, we found substantial differences in macrophage infiltration between the high- and low-IRRS groups. M0, M1, and M2 macrophage infiltration was significantly higher in the high-IRRS group. TAMs are one of the most abundant immune cells in the TIME, and M1 macrophages are generally believed to have antitumor properties and may increase after immunotherapy, while M2 macrophages are believed to exert protumor effects and may decrease in immunotherapy responders (38). The explanation for the results we obtained could be that the increased M2 macrophage infiltration overwhelms the antitumor effects of increased M1 macrophage infiltration, which addresses the speculation that suppressing M2 macrophage function should be more effective than promoting M1 macrophage function in enhancing immunotherapy responses (36). Moreover, we found that the expression levels of multiple immune checkpoints and the levels of tumor-related mutation burden varied significantly between the high- and low-IRRS groups. Numerous studies have confirmed that expression levels of immune checkpoints and tumor-related mutation burdens correlate with OS and immunotherapy responses in LUAD (39-42). Collectively, these mechanisms may explain the lower OS in the high-IRRS group.

Our risk formula could effectively distinguish the differences in OS and responses to various anticancer drugs between the high- and low-IRRS groups. We also found that patients with high-IRRS in the GSE135222 cohort had a longer PFS than did those with low IRRS. Further, we have found that the IRRS correlated with drug sensitivities to cisplatin, erlotinib, gemcitabine, vinorelbine, and paclitaxel. In addition, we have found that the IRRS correlated with TP53 mutation burdens. It is well-known that TP53 mutation mediates multiple anticancer drug resistances (43), it is reasonable to deduce that TP53 mutation burden difference was one of the main mechanisms for the correlation between the IRRS and anticancer drug sensitivities. One possible explanation for this result is the patients with low IRRS in the cohort we examined might have had an overall higher severity of malignancy than did those in the high-IRRS group; therefore, although they might respond better to immunotherapy, the OS benefits were overwhelmed by the higher extent of malignancy. Incorporating more cohorts would clarify this discrepancy. We also found that the IC₅₀ of multiple anticancer drugs were higher in the low-IRRS group than in the high-IRRS group, and also the IPS score of CTLA4 is higher in low-IRRS group. These indicating that despite patients in the low-IRRS group potentially responding to better to immunotherapy, they are more likely to be resistant to other anticancer drugs. This result was not unexpected because immunotherapy is based on different mechanisms from chemotherapy or targeted therapies, and patients who might benefit less from immunotherapy might respond better to other anticancer modalities. This result highlights the importance of a comprehensive evaluation and assessment of each patient in clinical practice.

This study had some limitations which should be addressed. First, the data sources were publicly available databases, and the lack of local data might influence the application of this model in actual clinical practice. Second, the data sources used were open-source datasets. Given that the LUAD guidelines are updated rapidly nowadays, it is possible that the therapy regimens used in the datasets vary from current recommendations, and the results might differ if real-time data are used to construct the same prognostic model. Third, as mentioned above, the functions of the prognostic-related genes incorporated in the IRRS formula were beyond the scope of the present study, and further investigations are needed to clarify their roles and mechanisms in LUAD.

Conclusions

In conclusion, we successfully established and validated a prognostic prediction model for LUAD based on immunotherapy-responsive DEGs. The IRRS model could predict the overall and subgroup prognoses and was an independent predictor of OS in patients with LUAD. The most significantly varied biological process between the high- and low-IRRS groups was ribonucleoprotein complex biogenesis, and the IRRS correlated with TP53 and TTN mutation burdens. Moreover, the immune status and anticancer drug sensitivities between the high- and low-IRRS groups were significantly different. A nomogram based on the IRRS was also established. Our novel model can be practically used for prognosis prediction and clinical decision-making.

Acknowledgments

Funding: This study was supported by the Natural Science Foundation of Hunan Province (Nos. 2023JJ40828, 2023JJ60081, 2021JJ40879, 2023JJ40876); The Science Project of the Hunan Health Commission (No. 202203022422); The National Natural Science Foundation of China (No. 82150006); Chinese Society of Clinical Oncology Research Foundation (Nos. Y-HR2017-117, Y-HH202102-0060); Beijing Medical and Health Foundation (No. YWJKJJHKYJJ-F3046D); Qujiang District Quzhou City Life Oasis Public Service Center (No. BJHA-CRP-040); Wu Jieping Medical Foundation (No. 320.6750.2023-05-10); The National Natural Science Foundation of China (No. 82360377); Guangxi Natural Science Foundation (No. 2024GXNSFAA010047).

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-24-309/rc

Peer Review File: Available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-24-309/prf

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-24-309/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Chen P, Liu Y, Wen Y, et al. Non-small cell lung cancer in China. Cancer Commun (Lond) 2022;42:937-70. [Crossref] [PubMed]
2. Maomao C, He L, Dianqin S, et al. Current cancer burden in China: epidemiology, etiology, and prevention. Cancer Biol Med 2022;19:1121-38. [Crossref] [PubMed]
3. Ettinger DS, Wood DE, Aisner DL, et al. Non-Small Cell Lung Cancer, Version 3.2022, NCCN Clinical Practice Guidelines in Oncology. J Natl Compr Canc Netw 2022;20:497-530. [Crossref] [PubMed]
4. Wang MM, Zhang Y, Wu S, et al. Clinical outcomes of KRAS-mutant non-small cell lung cancer under untargeted therapeutic regimes in the real world: a retrospective observational study. Transl Lung Cancer Res 2023;12:2030-9. [Crossref] [PubMed]
5. Smyth MJ, Ngiow SF, Ribas A, et al. Combination cancer immunotherapies tailored to the tumour microenvironment. Nat Rev Clin Oncol 2016;13:143-58. [Crossref] [PubMed]
6. Ahmed KM, Veeramachaneni R, Deng D, et al. Glutathione peroxidase 2 is a metabolic driver of the tumor immune microenvironment and immune checkpoint inhibitor response. J Immunother Cancer 2022;10:e004752. [Crossref] [PubMed]
7. Lim JU, Lee E, Lee SY, et al. Current literature review on the tumor immune micro-environment, its heterogeneity and future perspectives in treatment of advanced non-small cell lung cancer. Transl Lung Cancer Res 2023;12:857-76. [Crossref] [PubMed]
8. Chen K, Liu S, Lu C, et al. A prognostic and therapeutic hallmark developed by the integrated profile of basement membrane and immune infiltrative landscape in lung adenocarcinoma. Front Immunol 2022;13:1058493. [Crossref] [PubMed]
9. Lin T, Wang H, He X. Identification of an Exosome-Related Signature for Predicting Prognosis, Immunotherapy Efficacy, and Tumor Microenvironment in Lung Adenocarcinoma. J Oncol 2022;2022:1827987. [Crossref] [PubMed]
10. Wang M, Zhang Y, Liu M, et al. Exploration of a novel prognostic model based on nomogram in non-small cell lung cancer patients with distant organ metastasis: implications for immunotherapy. Transl Lung Cancer Res 2023;12:2040-54. [Crossref] [PubMed]
11. Jiang Y, Ouyang W, Zhang C, et al. Prognosis and Immunotherapy Response With a Novel Golgi Apparatus Signature-Based Formula in Lung Adenocarcinoma. Front Cell Dev Biol 2021;9:817085. [Crossref] [PubMed]
12. Ouyang W, Jiang Y, Bu S, et al. A Prognostic Risk Score Based on Hypoxia-, Immunity-, and Epithelialto-Mesenchymal Transition-Related Genes for the Prognosis and Immunotherapy Response of Lung Adenocarcinoma. Front Cell Dev Biol 2021;9:758777. [Crossref] [PubMed]
13. Wang L, Wu P, Shen Z, et al. An immune checkpoint-based signature predicts prognosis and chemotherapy response for patients with small cell lung cancer. Int Immunopharmacol 2023;117:109827. [Crossref] [PubMed]
14. Wang Y, Xu Y, Deng Y, et al. Computational identification and experimental verification of a novel signature based on SARS-CoV-2-related genes for predicting prognosis, immune microenvironment and therapeutic strategies in lung adenocarcinoma patients. Front Immunol 2024;15:1366928. [Crossref] [PubMed]
15. Wu Z, Wu H, Dai Y, et al. A pan-cancer multi-omics analysis of lactylation genes associated with tumor microenvironment and cancer development. Heliyon 2024;10:e27465. [Crossref] [PubMed]
16. Han Y, Liu SM, Jin R, et al. A risk score combining co-expression modules related to myeloid cells and alternative splicing associates with response to PD-1/PD-L1 blockade in non-small cell lung cancer. Front Immunol 2023;14:1178193. [Crossref] [PubMed]
17. Wu Y, Verma V, Gay CM, et al. Neoadjuvant immunotherapy for advanced, resectable non-small cell lung cancer: A systematic review and meta-analysis. Cancer 2023;129:1969-85. [Crossref] [PubMed]
18. Provencio M, Carcereny E, López Castro R, et al. Real-world treatment patterns and survival outcomes for patients with stage III non-small cell lung cancer in Spain: a nationwide cohort study. Transl Lung Cancer Res 2023;12:2113-28. [Crossref] [PubMed]
19. Wan C, Zhang Y, Liu P, et al. Efficacy and safety of anti-programmed cell death protein 1 antibody combination therapy in patients with advanced experienced epidermal growth factor receptor-tyrosine kinase inhibitor-resistant lung adenocarcinoma: a retrospective cohort study. J Thorac Dis 2023;15:5648-57. [Crossref] [PubMed]
20. Gao Y, Sun X, Hou Y, et al. Efficacy and safety of immunoradiotherapy for advanced non-small cell lung cancer: a retrospective comparative cohort study. J Thorac Dis 2023;15:3182-96. [Crossref] [PubMed]
21. Diao H, Zhu P, Dai Y, et al. Identification of 11 potentially relevant gene mutations involved in growth retardation, intellectual disability, joint contracture, and hepatopathy. Medicine (Baltimore) 2018;97:e13117. [Crossref] [PubMed]
22. Bowles KR, Pugh DA, Liu Y, et al. 17q21.31 sub-haplotypes underlying H1-associated risk for Parkinson's disease are associated with LRRC37A/2 expression in astrocytes. Mol Neurodegener 2022;17:48. [Crossref] [PubMed]
23. Zheng Q, Ma Y, Chen S, et al. Identification of genetic loci jointly influencing coronary artery disease risk and sleep traits of insomnia, sleep duration, and chronotype. Sleep Med 2020;74:116-23. [Crossref] [PubMed]
24. Falola MI, Wiener HW, Wineinger NE, et al. Genomic copy number variants: evidence for association with antibody response to anthrax vaccine adsorbed. PLoS One 2013;8:e64813. [Crossref] [PubMed]
25. Xue Q, Wang Y, Zheng Q, et al. Prognostic value of tumor immune microenvironment factors in patients with stage I lung adenocarcinoma. Am J Cancer Res 2023;13:950-63. [PubMed]
26. Shakfa N, Li D, Conseil G, et al. Cancer cell genotype associated tumor immune microenvironment exhibits differential response to therapeutic STING pathway activation in high-grade serous ovarian cancer. J Immunother Cancer 2023;11:e006170. [Crossref] [PubMed]
27. Qu R, Ye F, Hu S, et al. Distinct cellular immune profiles in lung adenocarcinoma manifesting as pure ground glass opacity versus solid nodules. J Cancer Res Clin Oncol 2023;149:3775-88. [Crossref] [PubMed]
28. Kobayashi R, Kawabata-Iwakawa R, Sugiyama M, et al. Multiplexed genome editing by in vivo electroporation of Cas9 ribonucleoproteins effectively induces endometrial carcinoma in mice. Int J Cancer 2023;152:2331-7. [Crossref] [PubMed]
29. Shi ZD, Hao L, Han XX, et al. Targeting HNRNPU to overcome cisplatin resistance in bladder cancer. Mol Cancer 2022;21:37. [Crossref] [PubMed]
30. Mir C, Garcia-Mayea Y. LLeonart ME. Targeting the "undruggable": RNA-binding proteins in the spotlight in cancer therapy. Semin Cancer Biol 2022;86:69-83. [Crossref] [PubMed]
31. Wang Y, Xu J, Fang Y, et al. Comprehensive analysis of a novel signature incorporating lipid metabolism and immune-related genes for assessing prognosis and immune landscape in lung adenocarcinoma. Front Immunol 2022;13:950001. [Crossref] [PubMed]
32. Zhou N, Tang Q, Yu H, et al. Comprehensive analyses of one-carbon metabolism related genes and their association with prognosis, tumor microenvironment, chemotherapy resistance and immunotherapy in lung adenocarcinoma. Front Mol Biosci 2022;9:1034208. [Crossref] [PubMed]
33. Jie Y, Wu J, An D, et al. Molecular characterization based on tumor microenvironment-related signatures for guiding immunotherapy and therapeutic resistance in lung adenocarcinoma. Front Pharmacol 2023;14:1099927. [Crossref] [PubMed]
34. Zhou E, Wu F, Guo M, et al. Identification of a novel gene signature of lung adenocarcinoma based on epidermal growth factor receptor-tyrosine kinase inhibitor resistance. Front Oncol 2022;12:1008283. [Crossref] [PubMed]
35. Jia H, Tang WJ, Sun L, et al. Pan-cancer analysis identifies proteasome 26S subunit, ATPase (PSMC) family genes, and related signatures associated with prognosis, immune profile, and therapeutic response in lung adenocarcinoma. Front Genet 2022;13:1017866. [Crossref] [PubMed]
36. Ma S, Zhu J, Wang M, et al. A cuproptosis-related long non-coding RNA signature to predict the prognosis and immune microenvironment characterization for lung adenocarcinoma. Transl Lung Cancer Res 2022;11:2079-93. [Crossref] [PubMed]
37. Mousset A, Bellone L, Gaggioli C, et al. NETscape or NEThance: tailoring anti-cancer therapy. Trends Cancer 2024;S2405-8033(24)00056-6.
38. Hu J, Zhang L, Xia H, et al. Tumor microenvironment remodeling after neoadjuvant immunotherapy in non-small cell lung cancer revealed by single-cell RNA sequencing. Genome Med 2023;15:14. [Crossref] [PubMed]
39. Siolas D, Vucic E, Kurz E, et al. Gain-of-function p53(R172H) mutation drives accumulation of neutrophils in pancreatic tumors, promoting resistance to immunotherapy. Cell Rep 2021;36:109578. [Crossref] [PubMed]
40. Wang C, Lu T, Xu R, et al. A bioinformatics-based immune-related prognostic index for lung adenocarcinoma that predicts patient response to immunotherapy and common treatments. J Thorac Dis 2022;14:2131-46. [Crossref] [PubMed]
41. Sadhukhan P, Seiwert TY. The role of macrophages in the tumor microenvironment and tumor metabolism. Semin Immunopathol 2023;45:187-201. [Crossref] [PubMed]
42. Zhang X, Wang Y, A G, et al. Pan-Cancer Analysis of PARP1 Alterations as Biomarkers in the Prediction of Immunotherapeutic Effects and the Association of Its Expression Levels and Immunotherapy Signatures. Front Immunol 2021;12:721030. [Crossref] [PubMed]
43. Zong Z, Xie F, Wang S, et al. Alanyl-tRNA synthetase, AARS1, is a lactate sensor and lactyltransferase that lactylates p53 and contributes to tumorigenesis. Cell 2024;187:2375-2392.e33. [Crossref] [PubMed]