Cuproptosis and Immune-Related Gene Signature Predicts Immunotherapy Response and Prognosis in Lung Adenocarcinoma

Cuproptosis and associated immune-related genes (IRG) have been implicated in tumorigenesis and tumor progression. However, their effects on lung adenocarcinoma (LUAD) have not been elucidated. Therefore, we investigated the impact of cuproptosis-associated IRGs on the immunotherapy response and prognosis of LUAD using a bioinformatical approach and in vitro experiments analyzing clinical samples. Using the cuproptosis-associated IRG signature, we classified LUAD into two subtypes, cluster 1 and cluster 2, and identified three key cuproptosis-associated IRGs, NRAS, TRAV38-2DV8, and SORT1. These three genes were employed to establish a risk model and nomogram, and to classify the LUAD cohort into low- and high-risk subgroups. Biofunctional annotation revealed that cluster 2, remarkably downregulating epigenetic, stemness, and proliferation pathways activity, had a higher overall survival (OS) and immunoinfiltration abundance compared to cluster 1. Real-time quantitative PCR (RT-qPCR) validated the differential expression of these three genes in both subgroups. scRNA-seq demonstrated elevated expression of NRAS and SORT1 in macrophages. Immunity and oncogenic and stromal activation pathways were dramatically enriched in the low-risk subgroup, and patients in this subgroup responded better to immunotherapy. Our data suggest that the cuproptosis-associated IRG signature can be used to effectively predict the immunotherapy response and prognosis in LUAD. Our work provides enlightenment for immunotherapy response assessment, prognosis prediction, and the development of potential prognostic biomarkers for LUAD patients.


Introduction
Lung adenocarcinoma (LUAD), the primary histological type of lung cancer, is one of the most common malignancies and the leading cause of cancer morbidity and mortality in humans [1]. Nearly 70% of lung cancer patients have locally advanced or metastatic disease at diagnosis [2]. With advances in cancer genomics, a group of genes has been identified as drivers of LUAD, including mutations in the epidermal growth factor receptor (EGFR), c-MET, KRAS, and anaplastic lymphoma kinase (ALK) [3]. In patients with stage IV LUAD who do not have driver gene status changes, if there are no obvious contraindications, the use of immune checkpoint inhibitors (ICIs) is strongly recommended in clinical practice [4]. At the same time, LUAD is generally suitable for ICIs therapy due

Unsupervised Clustering for Cuproptosis-Related IRGs
First, IRGs and cuproptosis-related gene expression matrices were extracted from the TCGA-LUAD gene expression matrix of 516 tumor samples, and 64 cuproptosis-related IRGs were screened by coexpression analysis (Supplementary Materials Table S1). Then, univariate Cox regression analysis was applied to filter 17 cuproptosis-related IRGs that were significantly associated with LUAD prognosis (p < 0.05). Moreover, consensus classification of the 17 cuproptosis-related IRGs was conducted by using the R package "Consen-susClusterPlus" [25], and the parameters clusterAlg, distance, reps, pltem, and pFeature were set to pam, euclidean, 100, 0.8, and 1, respectively. Next, according to the optimal cluster number, survival analysis was employed to clarify the survival differences between the two clusters. Finally, based on the clustering results, data downscaling was performed using principal component analysis (PCA), principal coordinate analysis (PCoA), and t-distributed stochastic neighbor embedding (tSNE).

Derivation of the Cuproptosis Prognostic Signature
First, the clinical data of 476 cases in the TCGA-LUAD cohort and the expression matrix of cuproptosis-related IRGs were merged by sample name and split into a training cohort (70%) and testing cohort (30%). The cuproptosis-associated IRGs were then filtered by univariate Cox analysis with p-values less than 0.05 in the training set. The genes required for modeling were further selected by least absolute shrinkage and selection operator (LASSO) regression analysis using the R package "glmnet" [31]. Stepwise regression was employed to construct an optimal multivariate Cox regression model for the cuproptosisrelated IRGs screened by LASSO, and then the risk scores of each sample in the training and testing sets were calculated as follows: β and exp represent the regression coefficient and expression of each cuproptosisassociated IRG, respectively. i represents the number of each gene. To prevent bias and guarantee the robustness of the model, the testing set and whole cohort were divided into low-and high-risk groups based on the median risk scores of the training set.

Real-Time Quantitative PCR (RT-qPCR)
The tumor specimens used in this study were approved (Approval No. 69, 2022) by the Human Ethics Committee of the First Affiliated Hospital of Guangdong Pharmaceutical University. Lung cancer specimens and associated clinical data including survival data from 12 patients with LUAD were collected for the study. RT-qPCR was performed to examine the expression levels of three key genes in these 12 samples. Primer sequences of NRAS (F GTGGAGCTTGAGGTTCTTGC; R CTGGATTGTCAGTGCGCTTT), TRAV38-2DV8 (F CCTGTCTTGAATTTAGCATGGCTC; R GCGAATAACGAGAATCATCTGCC), and SORT1 (F GACCTTGGGGCTCTGGAATTATG; R CCCTTGATCTGTTGAAACGTGGA) were designed based on the gene sequences on NCBI and detected by RT-qPCR using cDNA as a template. NAnodrop2000c was used to detect the RNA concentration in the samples. The samples were heated at 70 • C for 5 min to disrupt the RNA secondary structure. The reagents were then sequentially added using the First Strand cDNA Synthesis Kit, and the samples were kept at 50 • C for 30 min, followed by heating at 85 • C for 5 min and cooling on ice to prepare the cDNA. The qPCR was performed using the Ipure SYBR Green qPCR Master Mix kit, and the gene expression levels were normalized to GAPDH levels. The prognostic risk score of patients was then calculated based on the results of RT-qPCR using Equation (1), and the median value of the risk score of the training set was applied to split the 12 patients into low-and high-risk groups.

Establishment and Assessment of the Nomogram
We first executed a univariate Cox regression analysis on risk scores and clinical signatures in the TCGA-LUAD cohort, followed by a multivariate Cox regression analysis. The nomogram was then plotted by the "regplot" R package [32]. The area under the curve (AUC) and C-index were computed by the R package "timeROC" to assess whether the predicted values of the model were aligned with the actual values [33]. The calibration and decision curves were plotted with the R packages "rms" and "ggDCA", respectively.

Biological Features Analyses
To further identify the prognostic differences between the two subgroups, we conducted gene set enrichment analysis (GSEA) [40] on the samples using "c2.cp.kegg.v7.5.1.symbols.gmt" and gene sets compiled from references (Supplementary Materials Table S3) as the annotated gene sets. The statistically significant p-value was less than 0.05. Gene Ontology (GO) analysis was conducted using the "clusterProfile" software (Version 4.2.2).

Gene Mutation and Immunotherapy Response Analysis
Mutation data of LUAD were extracted from the official TCGA website, and TMB was calculated separately for the two subgroups. The "maftools" [41] software (version 2.10.05) was employed to display the top 10 gene mutation maps for tow subgroups. Furthermore, we assessed statistical significance in immune checkpoint-related genes' expression and immunotherapy sensitivity in both subgroups.

Statistical Analysis
In this study, the Wilcoxon test was performed to analyze statistical differences between the two groups that were not normally distributed. Correlation between two variables with non-normal distribution was assessed by Spearman's rank correlation test. The results of the RT-qPCR were compared between two groups using an unpaired t-test. Difference in survival curves was evaluated by a log-rank test. The value 0.05 was determined to be the significance threshold for p-value. Data analysis and visualization of the results were implemented with R 4.1.3 software.

Confirmation of Novel Subtypes of Cuproptosis-Associated IRGs in LUAD
The flow chart of our study is shown in Figure 1. Through coexpression analysis of cuproptosis-associated genes and IRGs, 64 cuproptosis-associated IRGs were filtered. Then, 17 cuproptosis-associated IRGs with prognostic value were screened by one-way Cox regression analysis in 476 patients with LUAD ( Figure 2A). To explore the classification of cuproptosis subtypes in TCGA-LUAD, an expression matrix of 17 cuproptosis-related IRGs was analyzed using unsupervised clustering. A total of 9 clusters were applied to consensus clustering analysis (Figures 2B and S1A-G), and the optimal number of clusters identified by the delta area plot and cumulative distribution plot of consensus scores was 2 ( Figures 2D and S1H). We further analyzed the distribution of clinical features of LUAD patients in the new LUAD classification (2C). To better elucidate the clinical significance of the LUAD classification, we split the TCGA-LUAD cohort into two groups according to the LUAD classification for survival analysis, and the results indicated that cluster 1 had a worse prognosis compared with cluster 2. ( Figure 2E). Meanwhile, PCA was performed on the TCGA-LUAD expression matrix, and the results further confirmed the significant differences between the two subtypes ( Figure 2F), which was in accordance with the results of PCoA ( Figure S1I) and tSNE ( Figure S1J) analyses. Additionally, we found that 6 out of 17 cuproptosis-related IRGs were differentially expressed in the subtypes ( Figure S1K).

Differences in TME of Cuproptosis-Associated Subtypes
To better understand the survival differences in both subtypes, we next investigated the differences in the biological signatures and tumor microenvironment (TME) of the two clusters. The GSVA algorithm performed on the two subtypes demonstrated that immune activation-associated pathways activity was upregulated in cluster 2, such as antigen processing and presentation, T-cell and B-cell receptor-signaling pathways, and natural killer cell-mediated cytotoxicity ( Figure 3A). At the same time, pathways activity of tumor-associated biological processes, such as DNA damage repair, proliferation (cell cycle progression and tumor proliferation rate), and stemness (RAMALHO stemness UP and 2019 PNAS stemness), were upregulated in cluster 1, while cluster 2 had a significantly higher enrichment score of interstitial activation pathways (endothelium, cancer-associated fibroblasts, and pan-fibroblast TGF-β response signature (pan-F-TBRS)) than cluster 1 ( Figure 3B,C). ssGSEA was performed on two clusters to elucidate the characteristics of TME. The two subtypes demonstrated significant differences in the abundance of immunoinfiltrating cells, with cluster 2 having a dramatically higher infiltration abundance than cluster 1 ( Figures 3D and S2A). The 17 cuproptosis-related IRGs were also significantly associated with the abundance of these 28 immunoinfiltrating cells ( Figure S2B). To ensure that the results of the ssGSEA were not biased and to demonstrate that the above results are robust and accurate, five other algorithms, namely, ESTIMATE, TIMER2.0, MCP-counter, CIBERSORT, and xCell, were used (Figures 3E and S2C-G). two subtypes demonstrated significant differences in the abundance of immunoinfiltrating cells, with cluster 2 having a dramatically higher infiltration abundance than cluster 1 ( Figures 3D and S2A). The 17 cuproptosis-related IRGs were also significantly associated with the abundance of these 28 immunoinfiltrating cells ( Figure S2B). To ensure that the results of the ssGSEA were not biased and to demonstrate that the above results are robust and accurate, five other algorithms, namely, ESTIMATE, TIMER2.0, MCP-counter, CIBER-SORT, and xCell, were used (Figures 3E and S2C-G). (E) The ESTIMATE calculated the immune scores between the two subtypes. "**": p < 0.01; "***": p < 0.001; "****": p < 0.0001. The ESTIMATE calculated the immune scores between the two subtypes. "**": p < 0.01; "***": p < 0.001; "****": p < 0.0001.

Creation of the Cuproptosis-Associated Prognostic Signature
First, to obtain cuproptosis-associated signatures that could be employed to predict LUAD prognosis, 70% of TCGA-LUAD patients served as a training cohort and 30% as a test cohort (Table 1). Univariate Cox regression analysis and LASSO analysis were conducted to identify seven key genes in the training cohort ( Figure 4A,B). Then, the best prognostic risk model was screened using stepwise regression, and three prognostic genes were selected, including NRAS, T-cell receptor alpha variable 38-2/delta variable 8 (TRAV38-2DV8), and SORT1. A prognostic risk model for LUAD was developed using these three critical genes, and cuproptosis-related risk scores were calculated using the following equation:  In the TCGA-LUAD training cohort, the probability of OS was remarkably higher in the low-risk subgroup, and the prognosis was clearly worse for the high-risk subgroup ( Figure 4C,F). Likewise, the same results were found in the testing set ( Figure 4D,G). In the entire TCGA-LUAD cohort, progression-free survival was compared in both subgroups, demonstrating superior performance in the low-risk subgroup ( Figure 4E). The risk scores also accurately predicted the probability of OS for the entire cohort ( Figure  S3A,B). Furthermore, we performed the Spearman's test to examine whether the three prognostic genes' expression levels had a strong relevance to the risk score. The findings revealed that TRAV38-2DV8 and SORT1 had a negative correlation, while NRAS expression level was increased with rising risk scores ( Figure S3C-E). The median risk score (cutoff value: 1.00009) was applied to classify the training cohort, the test cohort, and the whole cohort into two groups.
In the TCGA-LUAD training cohort, the probability of OS was remarkably higher in the low-risk subgroup, and the prognosis was clearly worse for the high-risk subgroup ( Figure 4C,F). Likewise, the same results were found in the testing set ( Figure 4D,G). In the entire TCGA-LUAD cohort, progression-free survival was compared in both subgroups, demonstrating superior performance in the low-risk subgroup ( Figure 4E). The risk scores also accurately predicted the probability of OS for the entire cohort ( Figure S3A,B). Furthermore, we performed the Spearman's test to examine whether the three prognostic genes' expression levels had a strong relevance to the risk score. The findings revealed that TRAV38-2DV8 and SORT1 had a negative correlation, while NRAS expression level was increased with rising risk scores ( Figure S3C-E). Furthermore, we discovered that patients with different clinical stages and genders had considerably different risk scores ( Figure S3F-G). The stratified analyses illustrated that three cuproptosis-associated IRGs' signatures, regardless of clinical stage and gender, were utilized to reliably detect prognostic differences for patients in the high-risk subgroup ( Figure S3H-K). This demonstrated that the cuproptosis-related signature distinguished prognostic differences based on clinical characteristics. The results of the multivariate Cox regression analyses supported those of the univariate Cox regression analyses in that the prognostic influence of age, gender, stage, and TMB in patients declined risk stratification. However, risk stratification resulted in fluctuating prognostic implications for risk scores, indicating that risk stratification was linked to the prognostic value of this feature (Table S4). We performed RT-qPCR on tumor samples from 12 LUAD patients, and the results suggested significant differences in the expression of TRAV38-2DV8 and SORT1 in the low-and high-risk groups and were consistent with the trend in the TCGA cohort, while NRAS expression was not significantly different in the two risk subgroups ( Figure 4H). The possible reason for this biased result was the small sample size.

Gene Expression Pattern Analysis
The result of the single-cell data integration is shown in Figure 5A. We performed dimension reduction and clustering in lung adenocarcinoma tissue samples ( Figure 5B) and malignant pleural effusion samples ( Figure 5C), respectively, and then manually annotated the cell clusters of these two groups of samples according to marker genes. Cell types of lung adenocarcinoma tissue samples include immune cells (DCs, macrophages, mast cells, and T cells), endothelial cells, epithelial cells, and fibroblasts; B, DCs, epithelial cells, fibroblasts, Mono/Mac, neutrophils, plasma cells, and T cells are present in malignant pleural effusion samples ( Figure S4). The expression of SORT1 and NRAS was also observed to be higher in macrophages; however, only pDCs from samples of malignant pleural effusion showed considerable expression of TRAV38-2DV8 ( Figure 5D-F). Immunohistochemistry from The Human Protein Atlas (THPA) also confirmed that NRAS and SORT1 expression levels were upregulated in LUAD tissues ( Figure 5G,H).

Establishment of a Cuproptosis-Associated IRG Prognostic Risk Model
Univariate and multivariate Cox regression analyses of clinical characteristics and the risk score of LUAD patients ( Figure 6A,B) indicated that risk score had an independent influence on patient prognosis. A nomogram was set up using all the above features ( Figure 6C). The nomogram shows the predicted probability of survival for patient number 10. The total score was determined based on the score for each item calculated using the nomogram. The AUC, which were, respectively, 0.729, 0.687, and 0.632 on the ROC curve, measured how well the model predicted the OS probability of LUAD patients at 1, 3, and 5 years ( Figure 6D). The C-index of all features in the model for 5-year prognosis prediction of LUAD patients was calculated separately, indicating that the risk score had sufficient predictive power for the patients' prognosis ( Figure 6E). Decision curve analysis (DCA) demonstrated that the nomogram exhibited excellent predictive power with high clinical benefit ( Figure 6F). The calibration curves confirmed that the model was highly accurate in projecting the odds for LUAD patients' 1-, 3-, and 5-year OS ( Figure 6G). Univariate and multivariate Cox regression analyses of clinical characteristics and the risk score of LUAD patients ( Figure 6A,B) indicated that risk score had an independent influence on patient prognosis. A nomogram was set up using all the above features (Figure 6C). The nomogram shows the predicted probability of survival for patient number 10. The total score was determined based on the score for each item calculated using the nomogram. The AUC, which were, respectively, 0.729, 0.687, and 0.632 on the ROC curve, measured how well the model predicted the OS probability of LUAD patients at 1, 3, and 5 years ( Figure 6D). The C-index of all features in the model for 5-year prognosis prediction of LUAD patients was calculated separately, indicating that the risk score had sufficient predictive power for the patients' prognosis ( Figure 6E). Decision curve analysis (DCA) demonstrated that the nomogram exhibited excellent predictive power with high clinical benefit ( Figure 6F). The calibration curves confirmed that the model was highly accurate in projecting the odds for LUAD patients' 1-, 3-, and 5-year OS ( Figure 6G).

Biological Characteristics Analysis of the Prognostic Model
Genes that had differential expression between the two subgroups were capitalized on for GO enrichment analysis, and the top 10 terms were retrieved for visualization ( Figure 7A). A powerful connection between immune activation and prognostic risk models are suggested by the fact that GO terms were primarily enriched in molecular mechanisms of immune activation, such as pathways related to leukocyte, B cell, and lymphocyte-mediated immunity; MHC class II protein complexes; and pathways related to receptor activity. DNA replication, folate biosynthesis, and systemic lupus erythematosus (SLE) were among the pathways that were upregulated in the high-risk subgroup according to GSEA results for the KEGG pathway ( Figure 7B), whereas immune enterocolitis, asthma, allograft rejection, and immunodeficiency-related pathways were dramatically downregulated. (Figure 7C). Furthermore, we performed GSEA using 14 pathways selected from pathway databases and publications and indicated that the high-risk subgroup presented elevated activities of DNA damage repair and cell cycle pathways, while the low-risk subgroup highly expressed immune, stromal, and partial carcinogenesis-related pathways ( Figure 7D,E). We further investigated the relevance between enrichment levels of these 14 pathways and risk scores, and the results confirmed this finding, which largely illuminated the reason for the worse prognosis of patients in the high-risk subgroup ( Figure S5A,B).

PEER REVIEW 13 of 19
Genes that had differential expression between the two subgroups were capitalized on for GO enrichment analysis, and the top 10 terms were retrieved for visualization (Figure 7A). A powerful connection between immune activation and prognostic risk models are suggested by the fact that GO terms were primarily enriched in molecular mechanisms of immune activation, such as pathways related to leukocyte, B cell, and lymphocyte-mediated immunity; MHC class II protein complexes; and pathways related to receptor activity. DNA replication, folate biosynthesis, and systemic lupus erythematosus (SLE) were among the pathways that were upregulated in the high-risk subgroup according to GSEA results for the KEGG pathway ( Figure 7B), whereas immune enterocolitis, asthma, allograft rejection, and immunodeficiency-related pathways were dramatically downregulated. (Figure 7C). Furthermore, we performed GSEA using 14 pathways selected from pathway databases and publications and indicated that the high-risk subgroup presented elevated activities of DNA damage repair and cell cycle pathways, while the low-risk subgroup highly expressed immune, stromal, and partial carcinogenesis-related pathways ( Figure 7D,E). We further investigated the relevance between enrichment levels of these 14 pathways and risk scores, and the results confirmed this finding, which largely illuminated the reason for the worse prognosis of patients in the high-risk subgroup ( Figure  S5A,B).

Gene Mutation Landscape and Immunotherapy Susceptibility
We analyzed the difference in gene mutation distribution between the two subgroups in the TCGA-LUAD cohort, which demonstrated that tumor mutation burden (TMB) was higher in the high-risk population ( Figure 8A-C). In addition, our study found that im-

Gene Mutation Landscape and Immunotherapy Susceptibility
We analyzed the difference in gene mutation distribution between the two subgroups in the TCGA-LUAD cohort, which demonstrated that tumor mutation burden (TMB) was higher in the high-risk population ( Figure 8A-C). In addition, our study found that immune checkpoint-related genes had higher expression in the low-risk population ( Figure 8D). Lower clinical efficacy scores for ICIs indicated that patients are less sensitive to immunotherapy. Comparison of efficacy scores in two subgroups of the four treatment groups, either with anti-PD-1 and anti-CTLA4 alone or in combination, revealed higher immunotherapy scores in the low-risk subgroup, indicating that this subgroup was better suited for immunotherapy ( Figure 8E-H). immunotherapy. Comparison of efficacy scores in two subgroups of the four treatment groups, either with anti-PD-1 and anti-CTLA4 alone or in combination, revealed higher immunotherapy scores in the low-risk subgroup, indicating that this subgroup was better suited for immunotherapy ( Figure 8E-H).

Discussion
The majority of patients who die from malignancies are patients with lung cancer, and LUAD has the highest proportion of lung cancer occurrence [42]. Analysis of exon copy number profiling, mutation, rearrangements, and DNA methylation in LUAD samples, as well as further assessment of mRNA, miRNA, and protein expression, have demonstrated that LUAD lesions are dramatically heterogeneous [43,44]. Although deep learning and machine learning have been extensively employed for tumor classification and prognostic assessment in recent years [45], markers with accuracy of predicting patient prognosis are still limited. Therefore, an accurate identification of the molecular subtypes of LUAD is crucial for guiding antitumor therapy.
The process of tumor proliferation, invasiveness, metastasis, and migration is strongly linked to TME, and immunoinfiltrating cells play an obviously crucial role in tumor escape and antitumor therapy. It has been ascertained that both innate and adaptive immune cells are present in TME [46]. TME is closely associated with lung cancer heterogeneity and has a major impact on lung carcinogenesis and development. Its dy-

Discussion
The majority of patients who die from malignancies are patients with lung cancer, and LUAD has the highest proportion of lung cancer occurrence [42]. Analysis of exon copy number profiling, mutation, rearrangements, and DNA methylation in LUAD samples, as well as further assessment of mRNA, miRNA, and protein expression, have demonstrated that LUAD lesions are dramatically heterogeneous [43,44]. Although deep learning and machine learning have been extensively employed for tumor classification and prognostic assessment in recent years [45], markers with accuracy of predicting patient prognosis are still limited. Therefore, an accurate identification of the molecular subtypes of LUAD is crucial for guiding antitumor therapy. The process of tumor proliferation, invasiveness, metastasis, and migration is strongly linked to TME, and immunoinfiltrating cells play an obviously crucial role in tumor escape and antitumor therapy. It has been ascertained that both innate and adaptive immune cells are present in TME [46]. TME is closely associated with lung cancer heterogeneity and has a major impact on lung carcinogenesis and development. Its dynamic changes depend on infiltrating lymphocytes, cellular regulatory factors, immune-related genes, and protein expression profiles [47]. Todd Golub's team first proposed the concept and mechanism of cuproptosis, and many studies have found that cuproptosis-related genes are inextricably interconnected with different types of immune cell infiltration, such as LIAS and FDX1 [48,49]. LIPT1 is one of the genes associated with cuproptosis, and Lv H et al. found that in the TME of melanoma patients, the abundance of resting CD4 memory T cells had a positive correlation with the expression of this gene, whereas the abundance of effector T cells and natural killer (NK) cells was inversely correlated. [50]. Furthermore, cuproptosis also has a crucial influence on tumors in other ways, such as increasing glycolysis by downregulating PDHA1 expression and thus promoting gastric cancer development [51].
IRGs are closely bound up with the prognostic evaluation and treatment of tumors. Currently, multiple immune-associated characteristics are applied to identify different subgroups of prognostic patients with LUAD and to predict their prognosis. [3,[52][53][54]. Ma KY et al. explored the heterogeneity of IRGs expression in tumors using scRNA-seq data and demonstrated that it has paramount implications for immunotherapy efficacy [55]. However, the impact of cuproptosis-associated IRGs on oncogenesis, invasion and metastasis, prognosis, and treatment of LUAD is currently unclear. Therefore, our study attempted to investigate the comprehensive effects of cuproptosis-associated IRGs in LUAD. We also designated specific cuproptosis-associated IRGs that could be employed to identify subtypes and estimate the prognosis of LUAD.
In this study, we first screened for cuproptosis-related IRGs and separated the TCGA-LUAD cohort into two clusters by consensus clustering. Survival analysis elucidated a considerable difference in prognosis in both subtypes, with the cluster 2 cohort enjoying a lopsided survival advantage. To explore the reasons for the difference, GSVA enrichment analysis found that immune activation-related pathways were upregulated in cluster 2, such as NK cell-mediated cytotoxicity, antigen processing and presentation, T-and B-cell receptor signaling, and cell adhesion molecules. At the same time, various algorithms were conducted to estimate the immune composition of LUAD. The results illuminated that NK cells, CD8+ T cells, natural killer T (NKT) cells, M1 macrophages, and γδ T cells were more abundantly infiltrating in cluster 2. Studies have shown that M1 macrophages, NKT cells, NK cells, and γδ T cells contribute to the inhibition of tumor development [56][57][58]. We also found that compared with cluster 1, 6 of 17 cuproptosis-related IRGs used to construct cuproptosis subtypes, namely, CCL13, TLR7, HLA-DRA, TRAV38-2DV8, HLA-DMB, and TRBV25-1, were upregulated in cluster 2. Zhao W et al. found that CCL13 can be utilized to indicate the intratumoral heterogeneity of immunoinfiltration in lung carcinoma and its association with OS [59]. TLR7 stimulation leads to decreased expression of CD200R in immunoinfiltrating cells (CD45+), resulting in anti-tumor effects in multiple tumors [60,61]. Jie Mei et al. found that human leukocyte antigen-DR alpha (HLA-DRA) expression levels were reduced in NSCLC tissues, related to TME inflammation, and predictive of the response of NSCLC to ICIs. This study also found that HLA-DRA was expressed in tumor and immune cells [62]. Fling SP et al. showed that HLA-DMB, encoding a component required for the assembly of MHC class II intramolecular peptides, is a gene that functionally maps between HLA-DP and HLA-DQ, and is involved in class II antigen presentation [63].
Next, Cox regression analysis and LASSO analysis were executed to screen for three cuproptosis-associated IRGs (NRAS, TRAV38-2DV8, and SORT1) and to construct a prognostic risk model. Validation with internal and external data ascertained that the model can precisely distinguish risk subgroups of LUAD patients and provide a clinical reference to facilitate quantitative risk management of LUAD patients. In addition, we performed RT-qPCR using tumor samples from 12 patients and initially verified that TRAV38-2DV8 and SORT1 expression was remarkably elevated in the low-risk group. Due to the small sample size, NRAS expression had no significant differences in the two subgroups. According to a recent study, NRAS has a responsibility in immune cell infiltration in the TME and has an independent effect on the prognosis of LUAD patients. Its expression is strongly inversely correlated to the prognosis of LUAD [64]. Mutated or overexpressed NRAS promotes tumor lung colonization by regulating the expression of IL-8-related chemokines and initiating interactions between tumor cells, pulmonary blood vessels, and myeloid cells [65]. Several studies have demonstrated that SORT1 expression is elevated in several types of tumors, which is correlated with a poor prognosis, such as liver [66], gastric [67], prostate [68], and colorectal cancers [69]. We demonstrated that the expression levels of SORT1 and NRAS were obviously elevated in macrophages in LUAD samples by scRNA-seq data analysis. To date, the prominence of TRAV38-2DV8 and SORT1 for prognostic prediction and clinical therapeutic efficacy assessment of LUAD have not yet been thoroughly investigated in relevant studies, and more research is required to evaluate the effect of SORT1 and NRAS expression levels in macrophages on the occurrence and progression of LUAD. Nomograms can calculate event probabilities based on features in prognostic models to assess patient prognostic risk levels and are widely applied in prognostic prediction in cancer [70]. In this study, a nomogram was established to predict OS of patients at 1, 3, and 5 years. The calibration curve and the DCA curve revealed that the nomogram has a great clinical predictive value.
Moreover, GO analysis and GSEA were carried out in two subgroups in order to clarify the molecular mechanisms of the prognostic model. Cell cycle and DNA damage-repairrelated pathways were strongly expressed in the high-risk subgroup with a poor prognosis, whereas immunological-, stromal-, and certain oncogenic-associated pathways were highly expressed in the low-risk subgroup with improved prognosis. We also found that the low-risk subgroup was more responsive to immunotherapy and had a lower frequency of mutations. This largely explained why the low-risk subgroup had a better prognosis.
This research has several restrictions. First, the statistical study used data from open databases, and more research is required to determine the function of cuproptosisassociated IRGs, including NRAS, TRAV38-2DV8, and SORT1, in the TME of LUAD. Furthermore, we have not yet assessed the risk model using a clinical cohort due to the tiny clinical sample that was gathered.

Conclusions
In conclusion, using the cuproptosis-associated IRG signature, we identified three key cuproptosis-associated IRGs, NRAS, TRAV38-2DV8, and SORT1, and meaningfully categorized LUAD into two subgroups with prominent differences in prognosis, molecular features, and immunoinfiltration abundance. More importantly, we successfully established a prognostic prediction model based on the three IRGs connected to cuproptosis, and the model has superior performance in predicting the immunotherapy response and prognosis of LUAD. These three key genes expression profile in clinical samples from patients with LUAD is validated. Furthermore, we elucidated how cuproptosis-associated IRGs functioned through multiple biological pathways. Our findings provide meaningful insights into issues related to immunotherapy response assessment, prognosis prediction, and development of potential prognostic biomarkers in patients with LUAD.