Assessing the Prognostic Capability of Immune-Related Gene Scoring Systems in Lung Adenocarcinoma

Background Lung adenocarcinoma (LUAD) is the commonest of the subtypes of lung cancer histologically. For this study, we intended to analyze the expression profiling of the immune-related genes (IRGs) from an independently available public database and developed a potent signature predictive of patients' prognosis. Methods Gene expression profiles and the clinical data of lung adenocarcinoma were gathered from the Gene Expression Omnibus database (GEO) and The Cancer Genome Atlas (TCGA), and the obtained data were split into a training set (n = 226), test set (n = 83), and validation set (n = 400). IRGs were then gathered from the ImmPort database. A prognostic model was constructed by analyzing the training set. Then the GO and KEGG analysis was performed, and a gene correlation prognostic nomogram was constructed. Finally, external validation, such as immune infiltration and immunohistochemistry, was performed. Results The 110 genes were significant by univariate Cox regression analysis and randomized survival forest algorithm for the training set and showed a good distinction between the low-risk-score and high-risk-score groups in the training set (P < 0.0001) by screening for four prognosis-related genes (HMOX1, ARRB1, ADM, PDIA3) and validated by the test set GSE30219 (P=0.0025) and TCGA dataset (P=0.00059). Multivariate Cox showed that the four gene signatures were an individual risk factor for LUAD. In addition, the genes in the signatures were externally verified using an online database. In particular, PDIA3 and HMOX1 are essential genes in the prognostic nomogram and play an important role in the model of immune-related genes. Conclusion Four immune-related genetic signatures are reliable prognostic indicators for patients with LUAD, providing a relevant theoretical basis and therapeutic rationale for immunotherapy.


Introduction
Lung cancer is major cancer in the world and adenocarcinoma patients lack tumor-specific clinical symptoms in the early stage, and local infiltration and even distant metastasis occur in the middle and late stages of lung cancer, with poorer efficacy and overall survival rate. Currently, 60% of NSCLC is LUAD, of which NSCLS is the major component of lung cancer [1,2]. In the United States, in 2021, the number of new incidences in a year was 235,760, and the number of deaths reached 131,880 [3]. For treating patients with early LUAD, surgical lobectomy is the most frequently used method. However, 10-44% of them have a less favorable prognosis five years after surgery [4]. erefore, it is necessary to develop more effective biomarkers to obtain more effective prognostic models of LUAD patients to help early diagnosis of LUAD and treat different patients with reasonable treatment plans so that patients can receive more appropriate treatment for their own conditions and get the best treatment results.
It is generally accepted that cancer is an incredibly complex disease that involves the interplay between the tumor and the immune system [5]. e human immunerelated systems have been confirmed to play a critical part in the development and progression of aggressive cancer [6,7]. More recently, immune checkpoint blockade therapy has shown remarkable efficacy in treating solid tumors (melanoma, lung cancer, and so on). is treatment is based on a dynamic process between tumor and immunity, and most current studies predict efficacy based on marker levels before or at a point in time during treatment. Programmed deathligand (PD-1) and tumor mutation burden (TMB) are widely used as prognostic biomarkers for immunotherapy [8]. However, immune-related therapies are only available for some patients. ere are remarkable personal differences in therapeutic effectiveness, illustrating the complexity of cancer-causing mechanisms and the existence of tumor heterogeneity [9,10]. A few studies have reported that immune-related genes can predict prognosis in tumor survivors and provide potential targets for immunotherapeutic treatment [11][12][13]. However, the prognostic model of immune-related genes in early-stage LUAD patients is relatively rare.
is study explores early LUAD-related prognostic models based on immune-related genes based on gene expression datasets from GEO and TCGA. After we matched the immune-related genes list to the three cohorts, the resulting genes were screened by a series of methods. An immune signature based on four genes for prognosis with outstanding predictive capability was established. en enriched for the function of related genes. e infiltration of tumors into immune cells was also analyzed, and finally, the differential expression of the four genes in different cancers was analyzed. It helps physicians in the prognosis of early LUAD patients. Furthermore, it provides some theoretical basis for personalized treatment. Figure 1 shows the flow chart of this work.

Expression Data.
e Gene Expression Omnibus (GEO, GSE31210, GSE30219) and e Cancer Genome Atlas (TCGA) databases, which were widely used and generally recognized, were selected as the data sources for the study, from which clinical information and gene expression profiles of lung adenocarcinoma patients. GSE31210 were selected as training and GSE30219 and TCGA were selected as external validation sets, respectively. e total number of patients of the three cohorts was 709. e distribution of cases in the three data was training set GSE31210 (N � 226) and the validation sets are GSE30219 (N � 83) and TCGA (N � 400), respectively. Clinical information and gene expression profiles related to the dataset were collected according to the following methods. In addition, the GSE50081 (N � 128) dataset was also collected for the validation of the Kaplan-Meier analysis. In addition, the GPL570 microarray platform was annotated by probes to obtain the final expression profiles of the GEO data [13].

Selection of IRGs set. Immunology Database and
Analysis Portal (ImmPort) is an open repository of discipline-level human immunology databases for translation and clinical research.
e IRGs set was obtained from ImmPort, and the related R package "clusterProfiler" was used to match the corresponding ensemble and other related information. Finally, 1455 IRGs were obtained.

Development of the Prognostic Gene Signature.
Prognostic immunity-related genes were obtained by univariate Cox regression screening on survival status for each gene in the training set. Genes that showed significant differences (P < 0.01) in the Cox regression were subsequently analyzed using the stochastic survival forest algorithm (RSFA) for dimensionality reduction [14]. In addition, the equation is shown below: e meaning of each parameter is as follows: in which N is the amount of IRGs, Expr is the IRGs expression value, and coeff is the coefficient value of IRGs. In the following section, patients are classified into low and high-score-risk groups according to the median value.
We combined the genes generated in RSFA with sufficient scale to prevent overfitting during modeling. e obtained combinations of IRGs were tested for their performance using time-correlated receiver operating characteristic (time-ROC) analysis. e gene combinations with the largest area under the curve (AUC) were used for subsequent predictions, and the model is statistically significant.
e accuracy of the model for prognosis was subsequently validated in internal and external validation.
Clinical information about cancer patients in the dataset may also be an influential factor in prognosis. So we analyzed the relationship between the two in the following way. Evaluation of patient prognostic outcomes and their relationship to clinical characteristics Chi-square analysis was used to determine the association between the IRGs prediction model and clinical information. KM survival curves and log-rank test were used to characterize the interaction between the model and survival time. In addition, the clinical information in the training and validation sets was related to OS using multivariate Cox regression analysis. e final univariate Cox analysis was utilized to determine whether clinical features could assist immune-related gene models used together for prognosis.

Construction and Evaluation of an IRG Nomogram.
A predictive nomogram for gene expression correlation was developed. Subsequently, we used the calibration curve to detect its accuracy in the GEO and TCGA cohorts. In addition, the prediction bias of the nomogram was assessed for evaluation.

Detection of the Infiltration of Immune Cells with Prognosis-Related Irgs in Tumors.
Tumor Immune Estimation Resource (TIMER) is a website that uses gene expression profiling data to detect the infiltration of immune cells in tumor tissue. We applied it to explore the infiltration of immune cells with the four genes (HMOX1, ARRB1, ADM, and PDIA3) in LUAD and six infiltrating immune cells (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells) were available for analysis [9,15].

External Validation of the IRGs Signature.
Furthermore, TIMER database was utilized to validate the differential expression of the four prognostic IRGs in different cancers. e Human Protein Atlas (HPA) is based on proteomic, transcriptomic, and systems biology data and covers protein expression not only in tumors but also in normal tissues. It was used to verify the protein levels in the model.

Functional Enrichment Analysis.
Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were applied using the cluster-Profiler package to analyze the underlying biological processes of the prognostic IRGs [16].

Statistical
Analysis. Data analysis was performed by R (version 4.1.2). KM and cox analyses were applied in three separate datasets using the R package "survival." Cox analysis was utilized to identify prognostic IRGs and detect the models. Chi-square or rank-sum tests were applied to stratified variables. "ROC" and "TimeROC" can be used to validate the model's viability. Functional annotation is performed using the "Clusterprofiler" package. e P-values involved in the analysis were < 0.05, which was a significant statistical value.

Patient Population Information.
e number of LUAD patients was obtained from GSE31210 and GSE30219, which were 226 and 83, respectively, and 400 LUAD patients were gathered from the TCGA database. A sum of 1084 IRGs was Kaplan-Meier analysis identified for expression in the GSE31210 dataset. In Table 1, we can easily obtain the median age of the patients as 61 years. e number of male and female cases in the training set was 105 and 121, respectively, with 35 deaths and 191 living cases, and the median OS was 5.33 years. Data per sample from GSE31210 and TCGA data were distributed in pathological stages I-II of LUAD, and data per sample from GSE30219 data were distributed in T stages T1-T2 of LUAD.

Identifying four Prognostic IRGs in the Training Set.
After univariate cox analysis and ROC curve analysis, 110 prognosis-related genes were obtained by screening based on P < 0.01 and AUC > 0.6, and the results are shown in Supplementary Table S1. Nine IRGs were next obtained by importance scoring in a randomized survival forest, and subsequent permutations of these nine genes yielded 2 9 -1 � 511 prediction models (Figure 2(a) and 2)(b). en, these models were evaluated by AUC, and the combination of four genes, HMOX1, ARRB1, ADM, and PDIA3, was found to have the most considerable AUC value of 0.779 ( Figure 2)(c), which had the optimal predictive power. e risk score was shown in the following formula: where the gene name represents the expression level of this gene in a particular sample. From the equation, we can get the coefficients of 1.32004, −1.36259, 1.17434, and 1.44023 for HMOX1, ARRB1, ADM, and PDIA3, respectively.
We computed the risk scores for each patient using the RSF formula and plotted heat maps for the four genes (Figures 3(a)-3)(c). e number of deaths increased with the increase of risk score in the three datasets. In the high-riskscore group, two genes, HMOX1, ADM, were highly expressed and ARRB1 was lowly expressed, with the same pattern in all three data sets. However, there was no significant pattern in the expression of the PDIA3 gene. Among the relevant clinical features, relapses become more frequent as risk scores increase, but no significant trends were found for age and gender following the score during the three cohort. By cox analysis, ARRB1 was highly expressed in the low-risk-score group, which was a protective factor ( Table 2).

Evaluation of a Prognostic IRGs Signature.
A risk score was computed for each individual with the prognostic models of immune-related genes. In the training cohort, the KM analysis was used to confirm the difference in survival between the high-risk-score (N � 113) and low-risk-score (N � 113) populations classified by a median score [17,18]. In Figure 4(a), the 5-year survival incidence was 58.07% (N � 66) in the low -risk-score group vs. 32.74% (N � 37) in the high-risk-score group. In the training cohort, significant OS was observed in the low-risk-score group. To explore this in the validation set, the same methodology was then adopted for the GSE30219 ( Figure 4)(b), and the model showed significant differentiation capability (the 5-year survival rate was 69.05% in the low-risk-score group (N � 42) vs. 36.59% in the high-score-risk group (N � 41), log rank P � 0.0025). In Figure 4(c), TCGA data (N � 400) with a large sample size was used for survival prediction. LUAD in this cohort was also classified into high-or low-risk-score groups (the 5-year survival incidence was 13.50% in the lowrisk-score group (N � 200) vs. 9.00% in the high-risk-score group (N � 200), log rank P � 0.0025). In Figure 4(d), additional selected GSE50081 data were used for the validation of the survival model (the 5-year survival incidence was 50.00% in the low-score-risk group (N � 64) vs. 35.94% in the high-score-risk group (N � 64), log rank P � 0.16). Although not significant, it still has some predictive power. From the results of the four cohorts, it can be obtained that the four prognostic IRGs model had excellent prognostic ability.

Correlation Analysis Between Signature and Clinical
Features. In Table 3, the association between the IRGs signature and clinical features in the three datasets was validated by the chi-square test.
ere was a significant correlation between the pathological stage and the prognostic IRGs signature in the GSE31210 and TCGA cohort (P < 0.05). In the GSE31210 cohort, gender and age both had no obvious correlation with the prognostic IRGs signature. In the TCGA dataset, gender and the prognostic IRGs signature had a significant relationship (P < 0.05), but the age had not with this model. e GSE30219 cohort has a different clinical characteristic T-stage than the GSE31210 and TCGA cohort, and no statistically significant between the three clinical features and the prognostic IRGs signature in the GSE30219.
In Table 4, the IRGs signature was shown to be statistically significant by multivariate COX regression, affirming it as an independent adverse predictor. In all three data sets, signatures proved to be powerful predictors of clinical outcomes in LUAD patients (high-vs. low-risk, GSE31210, HR � 32.11, 95% CI 4.32-238.91, P < 0.001, n � 226;   Univariate Cox also suggests that prognostic IRGs signature is a risk factor. Interestingly, the P-values for the pathological stage in GSE31210 and TCGA were less than 0.5. Unfortunately, this clinical feature was not available in GSE30219, so this information could not be further judged by the three datasets. Gender and age cannot be considered as risk factor for prognosis in the three cohorts.

Exploring the Functions of the IRG Panel.
To begin with, we obtained 110 IRGs by survival analysis and AUC analysis. en we analyzed the pathways and functions of IRGs using GO and KEGG. Top 10 functional annotations of biological processes (BPs), cellular components (CCs), and molecular functions (MFs) were selected, among which the main outcomes of BPs were regulation of chemotaxis, leukocyte migration, and cytokine production. e preliminary results of CCs were linked to the membrane. MFs enrichment   3.6. IRGs Survival Prediction Nomogram. In Figure 6(a), to enable comprehensive prediction, predictive models for the four genes were translated into a nomogram to provide a visual projection of OS at 1, 3, and 5 years. For instance, different scores were obtained from the expression of the four genes, and the total scores were summed to calculate the survival rates at different years. As can be seen from the nomogram, the two genes HMOX1 and PDIA3 in the prognostic IRGs model require focused attention in the model for prediction.
To evaluate how well this nomogram simulates the real situation, calibration curves using 1000 bootstrap tests were plotted. As shown in Figure 6(b), the actual case and the predicted situation show a good agreement in the training set. Furthermore, the calibration curve still showed a good deal in the validation set (Figures 6(c) and 6)(d).
ese results indicate that our nomogram is a good predictor of reality, and both independent validation databases show that the nomogram has excellent utility.
In the Figure 6      10 Journal of Oncology ( Figure 6)(g). e results from these two datasets show that our prognostic IRGs model has excellent prognostic accuracy. immune cells. We utilized TIMER to explore the potential association of infiltration of immune cells with prognosisrelated IRGs in tumors. Interestingly, HMOX1, ARRB1, ADM, and PDIA3 were positively associated with B cells and macrophages (P < 0.05). ese results validate that the four IRGs are closely associated with the level of B cells and macrophages infiltration. Conversely, there is no significant association between these four genes and infiltration of purity, CD8+ T cells, CD4+ T cells, neutrophils, and dendritic cells (Figure 7(a)-7(d)).

IRGs Expression and Immunohistochemistry in Online
Databases. In Figure 8(a), PDIA3 was highly expressed, and the results were consistent with the training set and TCGA set. ARRB1 was lowly expressed, and the results were consistent with the three independent data sets, interestingly, in LUAD, statistically significant differences in the expression of HMOX1, ARRB1, and PDIA3. However, ADM were not differentially expressed between LUADnormal and LUAD-cancer. e differential expression levels of HMOX1, ARRB1, and PDIA3 in the prognostic IRGs model in LUAD also provide positive evidence that the model and LUAD have a strong correlation and can be used effectively for prognosis. According to the immunohistochemical analyses in the HPA database, representative protein expressions of HMOX1, ARRB1, ADM, and PDIA3 in normal and cancerous tissues were compared (Figures 8(b) and 8)(c). e staining intensity of these two genes (HMOX1, ADM) in tumor cells is higher than in normal pneumocytes. e staining intensity of ARRB1 in tumor cells is lower than in normal pneumocytes, while PDIA3 did not show a striking difference.

Discussion
When immunotherapy was not well established, early treatment of LUAD was importantly based on surgical lobectomy, with suboptimal outcomes because of the malignant nature of LUAD and the limited results of surgery. However, with the increasing maturity and clinical use of immune checkpoint inhibitors and targeted therapies, the prognosis of patients has improved significantly [19,20]. In this study, we choose three separate datasets, GSE31210, GSE30219, and TCGA, for the establishment and validation of an IRGs prognostic model. First, we carried out a dimensionality reduction on the obtained IRGs and subsequently avoided overfitting the model. Finally, we obtained 511 combinations by recombination and performed screening to obtain four gene prediction signatures.
Several recent studies have used IRG characteristics to predict the prognosis of LUAD patients. However, not all models showed an excellent prognosis. Wu et al. developed a 21-IRGs prognostic model. In this study, multivariate regression was used to screen for genes of relevance, but there is no practical method to prevent overfitting. Unfortunately, the AUC was only 0.61, 0.65, and 0.62 for 1, 3, and 5-year survival in the training set [21]. While the present study avoided overfitting and predicted 1, 3, and 5-year survival by four immune-related genes models, the AUC obtained was 0.755, 0.752, and 0.803 in the training set. e included IRGs for signature were HMOX1, ARRB1, ADM, and PDIA3. HMOX1 (Heme Oxygenase 1) is a protein-coding gene whose related functions include protein homodimerization activity and oxidoreductase activity [22]. Several published pieces of evidence support the overexpression of HMOX-1 in several human malignant tumors, including kidney, gastrointestinal, lung, and breast cancers [23]. Tsai et al. research discovered that HMOX1 is a negative prognostic NSCLC gene. Its high expression may increase the metastasis ability of cancer cells of NSCLC patients, and HMOX1 has the potential as a therapeutic target for NSCLC in the future [24]. ARRB1 (Arrestin Beta 1) plays an important role in transcription factor binding and ubiquitin protein ligase binding functions, and it is also a protein-coding gene [25]. Li et al. suggest that ARRB1 may be associated with the prognosis of LUAD, and it is considered as a new molecular biomarker for the diagnosis and prognosis of LUAD [26].
ADM (Adrenomedullin) gene encodes a protein that is a prohormone with several functions, including vasodilation, hormone secretion regulation, angiogenesis promotion, and antimicrobial activity [27]. Some studies suggest that NSCLC cells are essential targets for ADM, and it may regulate the activity of these malignant lung cells through differential induction of different early response genes [28].
Using nomogram to show the prediction of survival by IRGs signatures, it visually exhibits excellent prognostic capability. e calibration curves also show good agreement between the predictions value and the actual situation in training. We also present the cellular functions and pathway enrichment of the selected IRGs based on GO analysis and KEGG analysis. According to GO analysis, regulation of chemotaxis, leukocyte migration, and cytokine production were significantly enriched functions and most genes associated with neuroactive ligand-receptor interactions based on KEGG-Gene-Concept Network.
By exploring the potential association of infiltration of immune cells with prognosis-related IRGs in tumors, we found that the four IRGs were closely associated with the level of B-cell and macrophage infiltration. HMOX1, ARRB1, and PDIA3 were significantly differentially expressed in LUAD-normal and LUAD-cancer. e proteins encoded by the three genes HMOX1, ADM, and ARRB1 also differ significantly in tumor cells and normal pneumocytes.
ese results suggest that the four IRGs have prognostic potential.
is study also has some restrictions. Firstly, the cohorts of GEO and TCGA databases we used for the construction of immune-related gene prognostic models are retrospective experiments and lack large clinical samples to prospectively validate the prognostic value of LUAD patients. Additionally, due to the lack of experimental exploration in these immune genes' potential functions and mechanisms, the potential functions and mechanisms of these immune genes need to be further validated in clinical trials.

Conclusion
In summary, our study systematically analyzed the expression, prognostic value, and potential functions of IRGs. Ultimately, a four-immune-related genes signature was constructed that can be independently used as a biomarker for OS prediction in LUAD. Eventually, our study provides an essential theoretical basis for further research on the role of IRGs in LUAD.

Conflicts of Interest
e authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as potential conflicts of interest.