Introduction

Lung cancer is the leading cause of cancer-related death worldwide.1 Lung adenocarcinoma (LUAD) is the most common histologic subtype of non-small cell lung cancer, accounting for approximately 40% of lung cancers.2,3 There is a need for targeted evaluation and treatment of patients based on different genotypes, particularly those that may affect the immune microenvironment of LUAD. Previous studies of receptor tyrosine kinase ROS proto-oncogene 1 (ROS1) have shown that ROS1 is aberrantly expressed in various cancers and can broadly regulate cell growth, survival, and proliferation through multiple signaling pathways, such as the RAS-RAF-MEK-ERK pathway. The ROS1 fusion gene can produce an oncogenic cancer receptor tyrosine kinase, and mutations in some ROS1 sites can confer resistance to immune checkpoint inhibitors. The study of the effect of ROS1 on LUAD has been focused on its fusion mutations, but the association between the expression level of ROS1 and immunity or prognostic value has yet to be well-studied in LUAD patients. Our study found that the expression level of ROS1 was significantly lower in LUAD compared with normal lung tissue, and higher ROS1 expression was associated with better overall survival (OS) and higher immune scores in LUAD patients. We then established and validated a ROS1-related prognostic risk model. The model had an excellent predictive effect on the ROS1 high-expression group.

In addition, we examined the association between risk scores and immune cell infiltration levels using algorithms and found that risk scores were largely negatively correlated with immune cell infiltration, which may be related to the expression of specific cytokines. The tumor microenvironment (TME) of high-risk score LUAD patients promoted the occurrence of immune escape. These genes downregulated in the high-risk group were significantly associated with immune system function. These findings suggest the potential functional role of the risk score in LUAD and highlight the mechanistic basis for a risk score to influence immune cell infiltration in the TME. Finally, we investigated the regulation role of the model key intercellular adhesion molecule 5 (ICAM5) in LUAD cells using in vitro experiments. These results showed that ICAM5 significantly inhibited the proliferation, invasion and migration of LUAD cells, which suggested that the downregulation of ICAM5 may be one of the mechanisms of the higher malignant degree of ROS1 high-expression patients.

Materials and Methods Data Download and Preprocessing

We downloaded LUAD gene expression profiles from The Cancer Genome Atlas (TCGA) public gene expression database (https://cancergenome.nih.gov/), screening 501 tumors and 59 normal samples for our study. Clinical phenotypes and prognostic information for these samples were extracted from the University of California Santa Cruz (UCSC) Xena database (https://xena.ucsc.edu/). Additionally, we acquired LUAD microarray data from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), including GSE26939, GSE72094, and GSE42127. Aggregated pan-cancer data, encompassing transcripts per kilobase million and clinical information for 33 cancer types and normal tissues, were downloaded from the UCSC Xena database, obtained from TCGA database and the Genotype-Tissue Expression project (https://gtexportal.org/home/). Immunohistochemistry (IHC) staining images of ICAM5 from LUAD and healthy lung tissue were obtained from the Human Protein Atlas database (https://www.proteinatlas.org/).

Prognostic Risk Model Establishment

We used the R package “estimate” (version 1.0.13) to calculate the immune score of tumor samples divided into different groups according to different ROS1 expressions and found that the P value of the immune score between the two groups was the smallest when there were 200 samples in the low-expression group, so based on the difference in immune scores, the samples were divided into control groups according to ROS1 expression: ROS1 high-expression group (H group) and ROS1 low-expression group (L group). 301 patients in the H group provided available prognostic information to develop prognostic characteristics. Differential expression analysis was performed using the R package “limma” (version 3.52.4). To screen for differentially expressed genes, and adjusted were defined, and the same criteria were used to obtain differential genes between high-risk and low-risk groups. Univariate Cox regression analysis was used to determine the genes associated with survival, the LASSO algorithm was used for variable selection, and tenfold cross-validation determined the penalty parameter (λ) for the model according to the minimum criteria. The multivariate Cox regression was then performed to obtain candidate genes. The candidate genes that satisfied the proportional hazards assumption and collinearity test were used to form the final prognostic risk model. The risk score based on the candidate genes was formulated as follows: , where represents the synergistic coefficient, and represents the relative expression of RNA. The median of the risk score is used as a cut-off value to distinguish between high- and low-risk groups. We also group data from GEO database sources similar to the criteria for grouping data from TCGA database sources.

Survival and Prognostic Analysis

Survival differences between high-risk and low-risk groups were estimated by Kaplan-Meier (KM) analysis. TCGA data were divided into different clinical subtypes, and KM survival curves were analyzed separately to verify the model’s reliability. Curves of subject operating characteristics over time were plotted, the area under the curve (AUC) was calculated, and the predictive accuracy of survival time was estimated and compared by risk score and different clinicopathologic factors. Univariate and multivariate Cox regression analyses were performed to search for independent prognostic factors. The relationship between clinicopathologic features and risk score was determined by correlation analysis. Nomograms are used to clinically predict patients’ 1-, 3-, and 5-year survival rates using clinical information and risk score, as both factors may impact survival. Survival probability calibration curves, KM survival curves, and Receiver Operating Characteristic (ROC) curves are used to determine the predictive performance of the nomograms. The dataset in the GEO database is used as the external validation dataset, and the H group obtained by the grouping standard is used as the verification group.

Pathway Dysregulation Analysis

To investigate which biological processes are associated with poor outcomes in patients in the high-risk score group, we performed transcriptome-based pathway dysregulation analysis using the R package “pathifier” (version 1.34.0). We used Spearman correlation analysis to determine the correlation coefficient and significance between pathway dysregulation score (PDS) and risk score. Normal samples were used as a control group. The gene set used (H collection: hallmark gene sets) was obtained from the MSigDB database (https://www.gsea-MSigDB.org/gsea/index.jsp).

Feature Enrichment Analysis

To explore the relationship between immunity and risk score, the relevant immune-related gene set was obtained based on the R package “msigdbr” (version 7.5.1), and the parameter category was set to “C5”. In addition, the functions “gseGO” and “gseKEGG” were used to analyze the pathway enrichment analysis of the differential genes between the high and low ICAM5 expression groups in LUAD. The R package “fgsea” (version 1.22.0) was used to analyze the enrichment of immune-related pathways in different risk score groups.

Immunological Analysis of the Prognostic Risk Model

To investigate the relationship between immune status and risk score, gene signatures of tumor-infiltrating lymphocytes were obtained from the TISIDB database. The ssGSEA analysis based on the R package “GSVA” (version 1.44.5) was used to calculate tumor-infiltrating lymphocytes, the MCP-counter algorithm analyzed the gene expression data, and the infiltration of immune cells was compared. Human leukocyte antigen (HLA) was also used to analyze the differences between the two groups. In addition, a list of immune-relevant genes was downloaded from the ImmPort database (https://www.immport.org/home) to check the expression of immune checkpoints in different risk groups. We aim to explore the model’s ability to differentiate LUAD patients.

Chemotherapy, Radiotherapy, and Immune Response Susceptibility Analysis for Prognostic Risk Model

To assess the relationship between chemotherapy response and risk score, we analyzed drug susceptibility using a “pRRophetic” R package (version 0.5) based on gene expression data.4,5 The parameter “tissueType” was set to the “lung”. The pRRophetic algorithm has been widely used in medical research.6–9 We calculated fractions for X-ray response, ultraviolet radiation (UV) response, and DNA damage repair using ssGSEA. The corresponding gene set was obtained from the MSigDB database (http://www.gsea-MSigDB.org/gsea/index.jsp). To explore the relationship between risk score and immunotherapy efficacy, we used the Tumor Immune Dysfunction and Exclusion (TIDE) score to infer immunotherapy responses in patients with low- and high-risk scores. We downloaded official scoring data from the TIDE database (http://tide.dfci.harvard.edu) to assess their response to anti-PD1 and anti-CTLA4 immunotherapy based on differences in TIDE scores between the two groups.

Model-Specific Analysis of Prognostic Risk Model

We established a prognostic model using defined ROS1 high-expression groups. Survival-related analysis was performed on patients in the ROS1 low-expression group in the TCGA and GEO cohorts, respectively, to determine whether the prognostic risk model had the specificity of the ROS1 high-expression group.

Pan-Cancer Analysis

Correlation plot between pan-cancer levels of ROS1 expression and patient survival obtained from the TIMER database (https://timer.cistrome.org/). Using the R package “easyTCGA”, we integrated pan-cancer data from TCGA database and the Genotype-Tissue Expression project to reveal differences in gene expression between pan-cancer and normal tissues. This data was then combined with pan-cancer survival information for analysis. The R package “TCGAplot” was employed to visualize disparities in HLA genes, immunostimulatory factors, chemokines, and chemokine receptors, ESTIMATE scores, and MSI in pan-cancer, comparing the ICAM5 high- and low-expression groups.

RNA Extracting and Quantitative Real-Time PCR

Total RNA was extracted from cells using AG RNAex Pro Reagent (Accurate Biotechnology (Hunan) Co., Ltd China) following the manufacturer’s protocol. Then, cDNA was synthesized with a reverse transcription kit (Accurate Biotechnology (Hunan) Co., Ltd China) and Quantitative Real-time PCR (qRT-PCR) was performed using the SYBR Premix Ex Tap Kit (Accurate Biotechnology (Hunan) Co., Ltd China). The primer sequences used are listed in Supplementary Table S1.

Cell Culture

A549 and PC9 cell lines were acquired from Procell Life Science & Technology Co., Ltd. A549 cells were grown in F12K medium (Gibco) enriched with 10% FBS, while PC9 cells were maintained in RPMI 1640 medium (Gibco). Both cell lines were incubated at 37 °C in a humidified atmosphere with 5% CO2. 1% penicillin/streptomycin was added to the medium.

Cell Transfection and Treatment

For ICAM5 overexpression, full-length human ICAM5 was cloned into pcDNA 3.1 vector. The empty pcDNA 3.1 vector was used for the negative control of ICAM5. These vectors were purchased from Keyybio (Shandong, China). A549 and PC9 cell lines were transfected with the expression vectors according to the grouping design. Transfection was conducted using jetPRIME® transfection reagent (PolyPlus transfection, Illkirch, France).

Proliferation Assay

Cells were fixed in 10% cold trichloroacetic acid for at least 24 hours, then stained with Sulforhodamine B sodium salt (SRB) (Sigma, USA) for 20 minutes, followed by washing with 1% acetic acid. After drying, 150µL of 10 mmol/L Tris was added, and the absorbance was measured at 562 nm using a microplate reader (Thermo Fisher, USA). For the proliferation assay, cells were seeded in 96-well plates at 2000 cells/well and fixed every 12 hours. Data were normalized to day 1 and presented as mean ± SD, with six replicates for each experiment.

EdU Incorporation Assay

Cells were seeded in 12-well plates and, upon achieving 50% confluence, EdU incorporation rates were assayed using the BeyoClick™ EdU Cell Proliferation Kit labeled with Alexa Fluor 488 (Beyotime, Shanghai, China), following the manufacturer’s guidelines.

Wound Healing Assay

Cells were seeded in 12-well plates and, upon achieving 95% confluence, a scratch was made in the monolayer using a pipette tip. Photographs were taken after 48 hours, and changes to the scratched area were quantified using the ImageJ software.

Transwell Assay

The Transwell assay was carried out using a Transwell plate, where the upper chamber was loaded with 200 µL of cell suspension containing 4×104 cells without the addition of fetal bovine serum (FBS). Concurrently, the lower chamber was filled with 600 µL of medium containing 20% FBS. Following a 36-hour incubation period, the cells present in the lower chamber were fixed and stained using crystal violet. Post-staining, the cells were photographed and subsequently counted using the ImageJ software.

Statistical Analysis

All statistical analyses were performed in R version 4.2.2 and GraphPad Prism 8.0 software. The Wilcoxon test was used to compare the differences between continuous data sets. The Kruskal–Wallis test was used to compare differences between multiple groups. The Shapiro–Wilk test was used to determine the normality of the data; if the data were normal, Pearson correlation analysis was used; otherwise, Spearman correlation analysis was used to determine the correlation between the two data sets. p < 0.05 was considered statistically significant.

Results

This study’s overall design and flow chart are shown in Figure 1. This study encompassed 560 samples of patients with LUAD from TCGA database and 618 LUAD patients from the GEO database. The dataset’s primary demographic and patient characteristics are summarized below: Figure 1.

Figure 1 The overall study design and workflow.

Transcription Levels and Immune Correlation of ROS1 in LUAD Patients

ROS1 expression differs significantly between various cancers and their corresponding normal tissues, and ROS1 expression is significantly downregulated in LUAD compared to normal tissues (Figure 2A). Its expression was significantly correlated with patient outcome, and LUAD patients with high ROS1 expression have more prolonged survival (Figure 2B). We found that there was a significant positive correlation between ROS1 expression and the immune score calculated by the ESTIMATE algorithm (Figure 2C; Supplementary Figure S1A), so we continuously grouped according to ROS1 expression and compared the immune scores between the two groups to obtain the smallest p value between the two groups and the corresponding ROS1 expression (Supplementary Figure S1B). In this way, we distinguished between the high ROS1 expression group (H group) and the low ROS1 expression group (L group). The results of the ssGSEA analysis indicated that multiple anti-tumor immune cells were significantly different between the two groups (Figure 2D). The ESTIMATE analysis showed that the H group had better immunity (Figure 2E).

Figure 2 From the TIMER database, (A) differential expression of ROS1 gene, (B) relationship between ROS1 expression in different tumors and patient prognosis (blue indicates decreased risk and red indicates increased risk). In the TCGA cohort, (C) correlation between ROS1 expression and immune score, (D) ssGSEA analysis result, and (E) ESTIMATE analysis results of H and L groups. *0.01<p<0.05; **0.001<p<0.01; ***0.0001<p<0.001; ****p<0.0001; ns, not significant.

Development and Accuracy Verification of 17 Gene Models

Differential analysis in the H and L groups revealed 316 differentially expressed genes (adjusted P value< 0.05, ) (Figure 3A), of which 276 genes were upregulated, and 40 genes were downregulated. A total of 96 genes were further identified by univariate Cox regression analysis, and 18 candidate genes were obtained by lasso regression analysis to reduce dimensionality further (Supplementary Figure S2A-B). Seventeen genes passed collinearity and time-dependent tests for multivariate Cox regression analysis and forest map (Supplementary Figure S2D). A concordance index of 0.73 indicated good predictive performance for the model. We also obtained the risk score correlation coefficient (Supplementary Table S2). We found a significant negative correlation between ROS1 gene expression and risk score (Supplementary Figure S2C) and a significant correlation between prognostic gene expression and ROS1 expression (Supplementary Figure S2E).

Figure 3 (A) Differentially expressed mRNA between the H and L groups. In the H group: (B) The risk score distribution, OS status, and heat map of the 17 gene-based prognostic risk model. (C) KM survival curves grouped by the median risk score. (D) The ROC curve analysis of the prognostic risk model for predicting OS. (E) ROC curves of the risk score and clinical-pathological factors with OS. (F) Uni-Cox and multi-Cox regression analyses of clinical characteristics and risk score with OS. (G) Nomogram of patients with LUAD.

The risk score-based patient distribution of the TCGA cohort and the heat map of the associated genes for the corresponding patients are plotted (Figure 3B). Analysis of the KM curves of both cohorts clearly showed that patients in the high-risk group had significantly shorter OS than those in the low-risk group (Figure 3C; Supplementary Figure S3A). By plotting ROC curves for 1-, 3-, and 5-year, the AUC maximums for the high ROS1 expression group in the TCGA cohort were 3-year, and the AUC value for 3-year survival was 0.75 (Figure 3D). In addition, the predicted performance of the risk score (AUC = 0.75, 3-year) outperformed the version of the stage, T stage, N stage, M stage, gender, and age (AUC = 0.68, 0.60, 0.62, 0.56, 0.45, and 0.55) (Figure 3E). We performed uni-Cox and multi-Cox regression analyses using risk scores and clinical features to make the model more clinically significant. The results showed that risk scores, T stage, N stage, M stage, and the stage all harmed patient outcomes. In contrast, risk scores and the stage were independent prognostic factors for OS in LUAD patients (Figure 3F). The Nomogram is constructed using risk scores and the stage to facilitate clinical work (Figure 3G). The calibration curves fitting for 1-year, 3-year, and 5-year OS show good agreement (Supplementary Figure S3B). ROC curves predicted with a risk score and stage showed improved predictive performance (Supplementary Figure S3C). These results indicate that risk scores and Nomograms have good predictive power.

To further evaluate the prognostic ability of the model, we grouped LUAD patients in the TCGA cohort by age, gender, stage, T stage, N stage, and M stage to study the relationship between the risk score and OS of LUAD patients in general clinicopathological variables. For different classifications, patients in the low-risk group had significantly longer OS (Supplementary Figure S4A–S4E, Figure S4G, Figure S4I, Figure S4K-S4L) than those in the high-risk group except for T 3–4 stages, N 3–4 stages, and M1 stage (Supplementary Figure S4F, S4H, S4J). This result may be explained by the limited number of patients with advanced non-small cell lung cancer, the considerable heterogeneity within the tumor, the patient’s OS being greatly affected by factors such as tumor metastasis, and the patient’s short survival. Therefore, these results suggest that prognostic models can also predict the prognosis of LUAD patients with different age groups, genders, all stages, T 1–2 stages, N 0–1 stages, and M0 stage.

Risk Score-Related Pathway Dysregulation Analysis and Biological Processes Identification

We obtained a PDS that showed the extent of dysregulation of specific biological processes by transcriptome-based pathifier analysis and sorted patients in the H group by the size of their risk score, resulting in a pathway dysregulation-risk score heatmap (Figure 4A). According to the correlation analysis between PDS and risk score, biological processes with positive correlation coefficients are thought to correlate with poor prognosis in patients in the high-risk score group. The results showed that dysregulated unfolded protein responses were most correlated with risk scores, which could be used by cancer cells to relieve endoplasmic reticulum stress and promote proliferation and metastasis.10 At the same time, the mTORC1 signaling pathway, which promotes multiple anabolic pathways, promotes cell growth.11 The ssGSEA analysis of the two dysregulated pathways showed significant differences between the high and low risk groups in unfolded protein responses and mTORC1 signaling pathway, and this difference might be one of the reasons why the tumors in the high-risk group had a higher malignancy (Supplementary Figure S5A). Other dysregulated pathways such as proliferative-related processes (eg, E2F targets, G2M checkpoint, Myc targets v1, Myc targets v2), metabolism-related processes (eg, glycolysis, fatty acid metabolism, and bile acid metabolism), and homeostasis-related pathways (eg, hypoxia, interferon alpha response, inflammatory response, UV response) were highly correlated with the risk scores, which may also be a potential mechanism for the poor prognosis of the patients in the high-risk group. Through GO analysis, we found that GO analysis of biological processes showed that up-regulated differential genes between high and low risk groups were concentrated in proliferation-related pathways (eg, mitotic nuclear division and chromosome segregation); GO analysis of cellular components showed that spindle, chromosome, centromeric region, and midbody were predominantly enriched; and GO analysis of molecular function showed upregulated gene enrichment in terms of microtubule binding as well as extracellular matrix structural constituents (Figure 4B); In contrast, GO analysis of biological processes showed that downregulated genes were significantly concentrated in immune-related pathways (eg, humoral immune response and leukotriene biosynthetic process), GO analysis of cellular components showed that clathrin-coated endocytic vesicle, multivesicular body, etc. were significantly enriched, and GO analysis of molecular function showed that cargo receptor activity significant enrichment (Figure 4C). The results of the GSEA analysis clearly showed differences in proliferative and metabolic-related processes in the two risk groups (Figure 4D). We conclude that dysproliferation, metabolic dysregulation, and immunosuppression are why the high-risk group has a lower survival rate than the low-risk group.

Figure 4 (A) An overview of the association between risk score and PDS (the correlation coefficient and significance are determined using Spearman correlation.) and the distribution of clinical phenotypes in different risk groups. (B and C) Go enrichment analysis showed the BP (biological processes) and CC (cellular components) of upregulated genes and downregulated genes between the high-risk score group and the low-risk score group. (D) Differences in pathway activities scored by GSEA between high- and low-risk patients. *0.01<p<0.05; ***0.0001<p<0.001; ****p<0.0001.

Immunological Analysis of Prognostic Risk Models

According to ssGSEA (Figure 5A), the enrichment scores of multiple immune cells varied significantly across different risk score groups. ESTIMATE analyses of both groups showed lower immune and stromal scores in the high-risk score group (Supplementary Figure S5B), indicating higher tumor purity. The MCP-Counter results showed that the absolute abundance of multiple immune cells varied significantly in different risk score groups (Figure 5B). HLA is the main histocompatibility complex in humans, which is closely related to the function of the human immune system and is also an essential genetic genome of the human immune system. We further explored the differences in the expression of HLA-related genes between high-risk and low-risk patients. We found that except for HLA-A, HLA-G, and HLA-F-AS1, there were significant differences in the expression levels of other HLA-related genes between high-risk and low-risk groups and that these genes were expressed more in low-risk patients (Figure 5C). These findings suggest that differences in HLA genes may be one of the reasons for differences in antitumor immunity between high and low risk groups. At the same time, we found significant differences between some immune checkpoints and cytokines or their ligands (Figure 5D). Previous articles have reported on their function.12–19 CXCR6 is critical for CD8+ cytotoxic T cell-mediated tumor control and is the strongest indicator of all chemokine receptor genes for a good quality of the immune infiltrate that prolongs patient survival in several immunogenic human cancer types.13 Spearman correlation analysis showed that risk score was significantly correlated with CD8+ T cell infiltration and CXCR6 expression, and CXCR6 was significantly correlated with CD8+ T cells (Figure 5E).

Figure 5 (A) Boxplots of ssGSEA. (B) Boxplots of MCP-counter. (C) HLA-related gene expression level in high- and low-risk score group patients. (D) Boxplots of immune checkpoints and cytokines or their ligands. (E) Spearman correlation analysis between risk score and CD8+ T cells; risk score and CXCR6 expression; and CXCR6 and CD8+ T cells. *0.01<p<0.05; **0.001<p<0.01; ***0.0001<p<0.001; ****p<0.0001; ns, not significant.

Comparison of Radiotherapy, Chemotherapy, and Immune Response Between the Two Risk Groups

We validated seven DNA damage response (DDR) pathways and found they were upregulated in the high-risk group (Figure 6A). The high-risk group was more responsive to ultraviolet, gamma, and X-rays (Figure 6B), indicating that patients in the high-risk group had more excellent radiation resistance. This result was consistent with better DNA damage repair ability in the high-risk score group.

Figure 6 Boxplot of ssGSEA scores for (A) injury repair responses and (B) radiation responses between the high- and low-risk score groups. (C) A box plot of susceptibility scores for ten chemotherapy drugs in two groups. (D) Boxplot of TIDE score between high- and low-risk score groups. *0.01<p<0.05; ****p<0.0001.

A total of 138 chemotherapy drugs were analyzed using the “pRRophetic” R package, and 87 drugs showed a statistically significant difference (P<0.05) between the two groups (Supplementary Table S3). The sensitivity of four commonly used chemotherapy drugs for LUAD (cisplatin, paclitaxel, gefitinib, and docetaxel) and six other oncology chemotherapy drugs (PLX4720, EHT.1864, doxorubicin, vinblastine, cytarabine, and sorafenib) was statistically significant in both risk groups (Figure 6C). The abnormalities in the proliferation of tumor cells in high-risk patients provide a basis for their sensitivity to cisplatin, paclitaxel, docetaxel, and other drugs that inhibit tumor cell proliferation. PLX4720 and EHT.1864 were the two drugs with the greatest difference in sensitivity between the two risk groups.

To explore the role of prognostic risk models in predicting immunotherapy response in the future, we compared the TIDE score (Figure 6D) between the two risk groups. The results showed that the low-risk score group had a lower risk of immune evasion and appeared more sensitive to immunotherapy. We also correlated ROS1 expression with TIDE scores in the TCGA cohort. The results (Supplementary Figure S5C) showed that ROS1 expression significantly correlated negatively with the TIDE score.

Model-Specific Analysis of Prognostic Risk Models

The KM analysis plot (Figure 7A) of patients in the L group in the TCGA cohort shows that the risk score does not significantly distinguish OS from all LUAD patients. The same conclusion was validated in the GEO cohort (Figure 7C). The AUCs for 1-year, 3-year, and 5-year in the TCGA cohort were 0.69, 0.62, and 0.51, respectively (Figure 7B). Univariate cox analysis (Figure 7D) and multivariate cox analysis (Figure 7E) of patients in the L group in the TCGA cohort showed that risk scores had no significant effect on survival and that the N stage was an independent risk factor for prognosis. In addition, we found that the risk score of the L group of the TCGA cohort did not outperform the prediction of the stage, N stage (3-year, AUC = 0.61, 0.69, 0.67) (Figure 7F). Therefore, risk scores are not applied to patients in the L group.

Figure 7 (A) KM survival curves in the L group in the TCGA cohort and (C) in the L group in the GEO cohort. (B) The ROC curve analysis of the prognostic risk model for predicting OS in the L group in the TCGA cohort. (D) Uni-Cox, and (E) multi-Cox regression analyses of clinical characteristics and risk score with OS in the L group in the TCGA cohort. (F) ROC curves of the risk score and clinical-pathological factors with OS.

ICAM5 Affects the Biological Behaviors of LUAD Cells in vitro

By pan-cancer analysis, we found significant differences in ICAM5 gene expression in a variety of tumors (Figure 8A). IHC-stained images of ICAM5 in normal lung tissue and LUAD showed that protein expression of ICAM5 was lower in LUAD (Figure 8B). The ESTIMATE algorithm showed that in LUAD, the ESTIMATE scores were significantly higher in the ICAM5 high-expression group (Figure 8E). By calculating the correlation of gene expression with immune-related genes and the differences in HLA gene expression in different ICAM5 expression groups, we found that ICAM5 gene expression was significantly correlated with the expression of HLA genes, immunostimulatory factors, chemokines, and chemokine receptors in a variety of tumors (Figure 8C and D and Supplementary Figure S6A-B). In LUAD, the expression of ICAM5 gene was significantly reduced compared with normal tissues, meanwhile, the expression level of this gene had a significant positive effect on OS, disease-specific survival (DSS) and progression-free interval (PFI) (Figure 8F). KM survival analysis showed a significant difference in OS, DSS and PFI between the two groups (Figure 8G and Supplementary Figure S6C-D). In addition, the radargram showed that ICAM5 gene expression was significantly positively correlated with MSI in LUAD (Supplementary Figure S6G). We investigated the KEGG, GO pathway enriched between the two groups using GSEA (Figure 8H and I and Supplementary Figure S6E-F). We found that the low-expression ICAM5 group was associated with tumor malignant phenotypes such as cell cycle, DNA replication, and sister chromatid segregation, while the high-expression ICAM5 group was associated with MHC class II protein complex assembly, antigen processing and presentation of exogenous antigen. Thus, we recognized that ICAM5 gene has an important impact on tumor survival and immunity, especially in LUAD, which is of high research value. Next, we performed a series of in vitro experiments to elucidate the effect of ICAM5 gene on the biological behavior of LUAD.

Figure 8 (A) Expression of ICAM5 in normal and tumor tissues of pan-cancer cohorts. (B) IHC staining images of ICAM5 in normal lung tissue and LUAD. (C) Differences in HLA genes expression in LUAD cohort. (D) Differences in immunostimulatory factors expression in LUAD cohort. (E) Differences in the results of the ESTIMATE algorithm in pan-cancer cohort. (F) Forest plot demonstrating the effect of ICAM5 expression in LUAD on OS, DSS, PFI and DFI. (G) KM survival analysis of OS in LUAD. (H, I) Enrichment plots from GSEA in the high-ICAM5 expression group and low-ICAM5 expression group. *0.01<p<0.05; **0.001<p<0.01; ***0.0001<p<0.001; ****p<0.0001; ns, not significant.

We transfected pcDNA3.1 into A549 and PC9 cells. The overexpression effect was detected using qRT-PCR (Supplementary Figure S7A-B). We performed the SRB assay to test cell proliferation, and the overexpression of ICAM5 significantly inhibited the proliferation of A549 and PC9 cells (Figure 9A and B). In addition, the overexpression of ICAM5 reduced the rate of EdU-positive cells in LUAD cells compared to the control group (Figure 9C–F). To study cell invasion, we performed the Transwell assay (Figure 10A–D) and found that the overexpression of ICAM5 significantly reduced the invasion of A549 and PC9 cells. The overexpression of ICAM5 also impaired the migration of A549 and PC9 cells by wound healing assays (Figure 10E and F). In conclusion, the overexpression of ICAM5 significantly inhibited the proliferation, migration and invasion of LUAD cells.

Figure 9 (A, B) SRB and (C-F) EDU assay detected the proliferative capacity of A549 and PC9 LUAD cells. *0.01<p<0.05; **0.001<p<0.01; ***p<0.001.

Figure 10 Migration and invasion of cells were detected employing the (A-D) Transwell assay and (E, F) wound healing assay in A549 and PC9 LUAD cells. *0.01<p<0.05; **0.001<p<0.01; ***p<0.001.

Discussion

Developing individualized treatment strategies and accurately predicting patient outcomes is exceptionally challenging due to the tumor heterogeneity and complex carcinogenic mechanisms in LUAD.20 The prognosis of cancer patients is closely related to the TME. The immune response in the TME is also an important determinant of tumor aggressiveness and progression. Notably, we found that ROS1 is a gene closely associated with immune cell infiltration in TME. In addition, abnormal expression of ROS1 in various cancers has been reported in the literature,21 and ROS1 expression identified by IHC has been found to be a stage-independent predictor of OS improvement.22 Therefore, we conducted this study to explore the immune landscape based on ROS1 expression, construct an immune-related prognostic risk model, and predict individualized treatment methods for patients with high ROS1 expression. Tumor purity and immune score are considered essential factors influencing the prognosis of cancer patients.23,24 In this study, we investigated the relationship between ROS1 expression and immune score and used this as a grouping criterion. A novel 17-gene prognostic risk model with good immune-related robustness was identified and constructed using the TCGA dataset. Its reliability was validated using the clinical subset and the GEO dataset to accurately discriminate between high- and low-risk patients. The AUC value of the ROC curve indicates that our risk model has a better prognostic value than other clinicopathological features, including age, gender, T stage, N stage, M stage, and stage. Based on univariate and multivariate Cox regression analysis, our prognostic risk model is an independent prognostic factor for patients with LUAD in the H group.

Based on this, this study examined the tumor immune microenvironment of high-risk and low-risk patients, revealing that high-risk patients have higher tumor purity and lower immune score. In our study, high-risk patients with high tumor purity had a poor prognosis. The difference in survival between high- and low-risk patients may be due to the higher frequency of mutations in critical pathways and changes in the TME associated with the risk score. Therefore, we calculated the pathway abnormalities of different risk groups and found that the high-risk group was significantly related to abnormalities in proliferation, metabolism, immunity, and other pathways. The immune landscape of the tumors also illustrates the poor degree of immune cell infiltration in patients in the high-risk group, especially CD8+ T cells, which play an essential role in tumor immunity. We have further explored the molecular mechanisms underlying immune microenvironmental differences.

Various immune cells interacting with tumor cells in the TME play an essential role in fighting tumors. The immune score can help us quantify the immune environment in tumors. Several studies25–27 have shown that high immune scores are associated with better prognosis. Similarly, in our study, high-risk patients with poor prognoses had lower immune scores. At the same time, we found significant differences in cell types between high-risk and low-risk patients, such as CD8+ T cells, B lineage, and so on. This further demonstrates the effectiveness and accuracy of the features constructed in this study in identifying high-risk patients.

It has been shown that there is a class of immune checkpoints that help the immune system determine which antigens it needs to respond to14 and that co-stimulatory pathways play a crucial role in regulating effector T cell function and responses to anti-PD-L1/PD-1 therapy.28 For example, CD28 is a co-stimulatory molecule that binds to ligands to induce T cell activation and differentiation, and ICOS promotes follicular helper T cell differentiation.14 Their high expression in the low-risk group may be one of the reasons why immunity to tumors is enhanced in the low-risk score group. The TIDE score results suggest that the low-risk group is more likely to benefit from immunotherapy. However, we found that patients in the low-risk group had higher TIGIT and no significant difference in CD274 between the two groups, which also underscores that the regulation of tumor immunity results from a multifaceted combination.

Chemokines regulate immune cell trafficking in tumors and are involved in tumor initiation, progression, and angiogenesis.29 IL-7 has various immune effects, including inducing the proliferation of naive and memory T cells without causing Treg cells,18 and CCL19 exerts tumor suppressive function by inducing an anti-tumor immune response and inhibiting angiogenesis.16 Studies30 have shown that the simultaneous production of IL-7 and CCL19 by genetically engineered tumor-responsive T cells can work synergistically with PD-1 blockade therapy to produce effective and long-lasting anti-tumor immunity. CCL5 exerts a chemotactic effect on T cells via CCR5,15,31 and the high expression level of CCL5 may be one of the reasons why CD8+ T cells are enriched in high-risk groups.

The radiation tolerance of cancer cells remains a significant limitation of radiotherapy. Ionizing radiation can cause DNA double-strand breaks, a lethal injury. When damage occurs, it triggers the DDR, which helps cells recover from radiation damage.32 These DDRs confer radiation resistance on tumors, leading to a poor prognosis.33 Our results indicate that both the DDR pathway and radiation response are upregulated in the high-risk group, suggesting that high radiation resistance may contribute to insufficient benefit from radiotherapy in patients in the high-risk group. It also shows that our prognostic risk model can predict radiotherapy response in LUAD patients.

We used drug susceptibility analysis to compare predicted IC50 values in high- and low-risk groups to confirm the ability of risk scores as treatment-guided biomarkers. For the four standard LUAD chemotherapy drugs cisplatin, paclitaxel, gefitinib, and docetaxel, there were statistically significant differences in the high-risk group compared to the low-risk group, except for gefitinib, which was more effective in the low-risk group. Our pathway analysis showed a significant dysregulation in proliferation in the high-risk group. This may be why high-risk groups are more sensitive to cisplatin, paclitaxel, and docetaxel. Based on the value of IC50, the correlation between drug target and risk score, the drug EHT.1864, which targets RAC, is a suitable treatment for patients in the low-risk group compared with the high-risk group. At the same time, we found that the IC50 of the high-risk group was often lower than that of the low-risk group, meaning that patients in the high-risk group were more sensitive to certain drugs despite having more malignant features. This finding highlights the importance of identifying high-risk patients using risk scores to identify practical pharmacotherapeutic approaches highly targeted to different patients.

Currently, the best-known biomarkers (eg, PD-1, CTLA4, etc). do not reliably guide the use of immune checkpoint blockades, resulting in minimal therapeutic benefit for cancer patients.34,35 The TIDE algorithm is a computational method developed by Jiang et al36 to predict the effect of immunotherapy by integrating immune dysfunction and tumor rejection, which is considered an alternative to a single biomarker. TIDE is superior to other biomarkers or indicators in predicting immune checkpoint inhibitor response.35,36 In this study, we evaluated the response of high- and low-risk patients to immunotherapy in terms of the TIDE score, and we found that the TIDE score was significantly higher in the high-risk group than in the low-risk group. Thus, low-risk patients may benefit more from immunotherapy. It is worth noting that although immunotherapy can benefit some lung cancer patients, there are still some patients who do not achieve the expected results after the use of immune checkpoint inhibitors.37–39 We also found that the expression of ICAM5 gene in the model was significantly different between normal and tumor tissues. Higher ICAM5 expression was associated with better prognosis, and ICAM5 was significantly associated with the immune microenvironment. As a member of the ICAM adhesion protein family, ICAM5 plays a role in both the immune system and the nervous system, and it is mainly expressed in the mammalian telencephalon.40 In previous understanding, DNMT1 and DNMT3a promote thyroid cancer development by promoting hypermethylation and transcriptional activation of ICAM5;41 and the overexpression of UCHL5 enhances bladder cancer cell proliferation and migration through the AKT/mTOR pathway via c-Myc, SLC25A19, and ICAM5 transformation.42 However, studies on the role of ICAM5 in LUAD are still lacking, although published Mendelian randomization studies43,44 have confirmed that genetically predicted levels of the ICAM5 gene, as well as plasma ICAM5 protein levels, are significantly associated with a reduced risk of LUAD. Our study showed that LUAD patients with high-expression ICAM5 had higher immune scores, more expression of immune-stimulation-related molecules, and significantly higher expression of most chemokines and their receptors. Meanwhile, in this study, we demonstrated by in vitro experiments that overexpression of ICAM5 was found to inhibit the malignant progression of LUAD, including proliferation, invasion, and migration. GSEA showed that the ICAM5 high-expression group was not significantly enriched in pathways related to tumor malignant phenotypes as compared to the ICAM5 low-expression group, which may explain the above phenomenon.

Our study has several limitations and drawbacks. Our study of the role of ICAM5 in tumor development was limited to in vitro experiments and lacked animal experiments for further exploration. Also, studies addressing the differences in sensitivity to chemotherapeutic agents between high and low risk groups were not supported by conducting relevant experiments. Finally, our study has inherent limitations in terms of retrospective study design, and we need more clinical databases for external validation. Therefore, large-scale prospective studies as well as additional in vivo and in vitro experimental studies are needed to confirm our findings.

Conclusion

The prognostic risk model we constructed provides an accurate clinical application for the prognosis of LUAD patients with high ROS1 expression, predicts the survival, immunotherapy, and chemoradiotherapy response of LUAD patients, and provides some guidance for the treatment and prognosis evaluation of LUAD. Finally, ICAM5, the key model gene, was shown to inhibit the malignant development of LUAD cells.

Data Sharing Statement

The datasets generated or analyzed during this study are freely available in the Cancer Genome Atlas (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga) and the UCSC Xena database (https://xenabrowser.net/datapages/?cohort=TCGA%20Lung%20Adenocarcinoma%20LUAD&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443). The GSE26939 database was downed from the website as follows, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE26939. The GSE72094 dataset was downloaded from the website as follows, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE72094. The GSE42127 dataset was downloaded from the website as follows, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE42127.

Ethics Approval and Consent to Participate

This study was approved by the Biomedical Research Ethic Committee of Shandong Provincial Hospital and conducted under the Declaration of Helsinki. Ethics approval number: NO.SWYX2024-460.

Acknowledgments

The authors are especially grateful to the TCGA platform and GEO platform for providing a large number of valuable datasets as well as funding from the Shandong Provincial Natural Science Foundation [Grant number ZR2021MH192].

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Funding

Shandong Provincial Natural Science Foundation [Grant number ZR2021MH192].

Disclosure

Gemu Huang and Qingtao Song are affiliated with Amoy Diagnostics Co., LTD. The authors declare that they have no other competing interests in this work.

References

1. Thai AA, Solomon BJ, Sequist LV, Gainor JF, Heist RS. Lung cancer. Lancet. 2021;398:535–554. doi:10.1016/s0140-6736(21)00312-3

2. Chen Z, Fillmore CM, Hammerman PS, Kim CF, Wong KK. Non-small-cell lung cancers: a heterogeneous set of diseases. Nat Rev Cancer. 2014;14:535–546. doi:10.1038/nrc3775

3. Xu JY, Zhang C, Wang X, et al. Integrative proteomic characterization of human lung adenocarcinoma. Cell. 2020;182:245–261.e217. doi:10.1016/j.cell.2020.05.043

4. Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One. 2014;9:e107468. doi:10.1371/journal.pone.0107468

5. Geeleher P, Cox NJ, Huang RS. Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines. Genome Biol. 2014;15:R47. doi:10.1186/gb-2014-15-3-r47

6. Liu Z, Guo C, Li J, et al. Somatic mutations in homologous recombination pathway predict favourable prognosis after immunotherapy across multiple cancer types. Clin Transl Med. 2021;11:e619. doi:10.1002/ctm2.619

7. Liu Z, Guo Y, Yang X, et al. Immune landscape refines the classification of colorectal cancer with heterogeneous prognosis, tumor microenvironment and distinct sensitivity to frontline therapies. Front Cell Develop Biol. 2021;9:784199. doi:10.3389/fcell.2021.784199

8. Liu Z, Liu L, Weng S, et al. Machine learning-based integration develops an immune-derived lncRNA signature for improving outcomes in colorectal cancer. Nat Commun. 2022;13:816. doi:10.1038/s41467-022-28421-6

9. Liu Z, Xu H, Weng S, Ren Y, Han X. Stemness refines the classification of colorectal cancer with stratified prognosis, multi-omics landscape, potential mechanisms, and treatment options. Front Immunol. 2022;13:828330. doi:10.3389/fimmu.2022.828330

10. Madden E, Logue SE, Healy SJ, Manie S, Samali A. The role of the unfolded protein response in cancer progression: from oncogenesis to chemoresistance. Biol Cell. 2019;111:1–17. doi:10.1111/boc.201800050

11. Ben-Sahra I, Manning BD. mTORC1 signaling and the metabolic control of cell growth. Curr Opin Cell Biol. 2017;45:72–82. doi:10.1016/j.ceb.2017.02.012

12. Chauvin JM, Zarour HM. TIGIT in cancer immunotherapy. J ImmunoTher Cancer. 2020;8:e000957. doi:10.1136/jitc-2020-000957

13. Di Pilato M, Kfuri-Rubens R, Pruessmann JN, et al. CXCR6 positions cytotoxic T cells to receive critical survival signals in the tumor microenvironment. Cell. 2021;184:4512–4530.e4522. doi:10.1016/j.cell.2021.07.015

14. Edner NM, Carlesso G, Rush JS, Walker LSK. Targeting co-stimulatory molecules in autoimmune disease. Nat Rev Drug Discov. 2020;19:860–883. doi:10.1038/s41573-020-0081-9

15. González-Martín A, Gómez L, Lustgarten J, Mira E, Mañes S. Maximal T cell-mediated antitumor responses rely upon CCR5 expression in both CD4(+) and CD8(+) T cells. Cancer Res. 2011;71:5455–5466. doi:10.1158/00