A New Prognostic Indicator of Immune Microenvironment and Therapeutic Response in Lung Adenocarcinoma Based on Peroxisome-Related Genes

Lung adenocarcinoma (LUAD) has been the major cause of tumor-associated mortality in recent years and exhibits a poor outcome. New data revealed that peroxisomes have a function in the regulation of the development and progression of several tumors. However, the prognostic values of peroxisome-related genes (PRGs) were rarely reported. Genomic sequence, mutation, and clinical data of 535 LUAD tissues were obtained from TCGA data sets. Within the TCGA cohort, a multigene signature was constructed with the assistance of the LASSO Cox regression model. Three GEO data sets, including GSE3141, GSE31210, and GSE72094, were obtained as validation cohorts. ROC assays, Kaplan-Meier methods, and multivariate assays were applied to examine the prognostic capacities of the novel signature. Gene Set Enrichment Analysis (GSEA) was performed to further understand the underlying molecular mechanisms. In this study, we identified 47 differentially expressed peroxisome-related genes (PRGs), including 25 increased and 22 decreased PRGs. A prognostic model of six PRGs was established. The univariate and multivariate Cox analyses both showed that the p value of risk score was less than 0.05. In LUAD patients, the strong connection between the risk score and overall survival was further verified in three other GEO data sets. TMB and cancer stem cell infiltration were shown to be significantly higher in the high-risk group in comparison to the low-risk group. The TIDE score of the group with the low risk was considerably greater than that of the group with the high risk. Several drugs, targeting PRG-related genes, were available for the treatments of LUAD. Overall, we developed a novel peroxisome-related prognostic signature for LUAD patients. This signature could successfully indicate LUAD patients' chances of survival as well as their immune system's responsiveness to treatments. In addition, it has the potential to serve as immunotherapeutic targets for LUAD patients.


Introduction
Lung carcinoma is one of malignant cancers with the highest incidence and the worst prognosis worldwide [1]. Among all lung carcinoma subtypes, lung adenocarcinoma (LUAD) is the most common and aggressive subgroup [2]. Despite recent breakthroughs in the diagnosis and treatment of malignancies, the clinical outcome for LUAD remains unsatisfactory [3,4]. Tumor excision, for example, can result in systemic inflammation and the release of cancer cells into the bloodstream, both of which aid cancer dissemination [5,6]. In this regard, a thor-ough knowledge of the potential mechanism underlying LUAD formation and progression is essential to enhance diagnosis and prognosis at the beginning.
Peroxisomes are metabolic organelles which play an important role in cellular redox balance and lipid metabolism [7]. Peroxisomal functions are very critical to reactive oxygen species homeostasis, bile acid synthesis, ether phospholipid synthesis, and fatty acid oxidation [8,9]. The regular functioning of cells depends on the peroxisome's free radical production and scavenging process being in equilibrium [10]. Peroxisomal disorders can be caused by mutations in genes coding for peroxisomal proteins, many of which demonstrate metabolic abnormalities [11,12]. Peroxisomes have recently been implicated in the formation and progression of a variety of cancers, including LUAD, according to recent researches [13,14]. In LUAD and other malignancies, a number of peroxisomal enzymes and metabolic processes are changed. Therefore, an in-depth understanding of the peroxisome process in LUAD could provide an important solution for the development of a new treatment method. Gene chips and high-throughput sequencing technology have advanced significantly in the last few years, indicating that the peroxisome's genetic signature can be utilized to predict LUAD's overall survival [15,16].
In this study, we identified the differentially expressed peroxisome-related genes (PRGs) between LUAD specimens and nontumor specimens. Then, using univariate and least absolute shrinkage and selection operator (LASSO) Cox regression analysis, we established a unique predictive signature from PRGs in LUAD cohorts. Moreover, we confirmed the prognostic value of the novel model and its association with immune cell infiltration. On the basis of our findings, we suggested that this novel signature could aid in developing personalized cancer treatment strategies and that the relevant PRGs could serve as potential LUAD therapeutic targets.

LUAD Sequencing Data and Peroxisome-Related Gene
Collection. TCGA database (https://cancergenome.nih.gov/) provides genomic sequence, mutation, and clinical data of more than 30 different types of cancer. For this study, genomic sequence, mutation, and clinical data of 535 LUAD tissues were obtained from it. In addition, the sequencing data of 59 normal tissues were downloaded for differential expression analysis. Of the 535 LUAD patients, 49 patients, which lacked follow-up time or were followed up for less than 1 month, were removed.

Differential Expression Analysis of Peroxisome-Related
Genes. The differentially expressed PRGs were screened according to FDR < 0:05 using "limma" package, and the result was demonstrated with boxplot using "ggpubr" package. Functional enrichment analyses of differentially expressed PRGs, including KEGG and GO analysis, were conducted using org "org.Hs.eg.db," "enrichplot," and "clusterprofiler" packages. The interaction between protein-protein was analyzed in STRING website (http://string.embl.de/).

Construction of a Prognostic Signature with LASSO Logistic
Regression. Based on differentially expressed PRGs, prognosisrelated PRGs were selected with the univariate Cox regression analysis using "survival" package (p value < 0.05). Based on prognosis-related PRGs, LASSO logistic regression was utilized to construct the prognostic signature using "glmnet" R package. The penalty parameter (λ) of LASSO logistic regression was estimated. The formula was displayed: Risk score = sum ð Gene i ' s coefficient * Gene i ' s expressionÞ. According to the median of risk score, the risk score of patients was classified into binary variables, including high risk and low risk. Survival analysis was performed between the two groups to assess whether there were differences in survival. ROC, PCA, and t-SNE analyses were utilized to evaluate the effectiveness of the signature. The Cox regression analysis was conducted to see if the signature was a significant independent predictor.

Validation of the Signature with GEO Data Sets.
To observe whether the performance of the signature was robust, three GEO data sets, including GSE3141, GSE31210, and GSE72094, were obtained as validation cohorts.
2.5. Analysis of Immune Microenvironment. To find out the difference of immune microenvironment between two risk groups, infiltration of immune cells and expressions of immune checkpoints were analyzed. Using the "estimate" package, the ESTIMATE method was used to estimate the stromal and immune cell proportions in each of the LUAD samples. ssGSEA was applied to quantify the enrichment degrees of sixteen immune cells in TCGA data set and GSE72094 using "GSEABase" and "GSVA" packages. The expression of 47 immune checkpoints was analyzed between the two risk groups, and the differentially expressed checkpoints were exhibited. TIDE (tumor immune dysfunction and exclusion) is a computational method that simulates tumor immune evasion. A high TIDE score indicates the greater tendency of tumor immune escape and the worse effect of immune checkpoint blockade treatment (ICB) [18]. TIDE score of patients in TCGA data set was downloaded from TIDE website (http://tide.dfci.harvard.edu/).

2.6.
Relationship between Risk Score and Ferroptosis. It has been reported that peroxisome induced ferroptosis by synthesizing polyunsaturated lipids [19]. To observe whether there were differences in the activity of ferroptosis between the two risk groups, we obtained four representative ferroptosisrelated genes from the literature, including ACSL4, TFRC, PTGS2, and CHAC1 [20]. Ferroptosis in different risk groups was observed by comparing the expression of these four genes in different risk groups.

Relationship between Risk Score and Cancer Stem Cell
Infiltration. Cancer stem cells may be the source of tumor cells and could cause drug resistance and lead to distant metastasis [21]. The cancer stem cell index, representing the infiltration degrees of cancer stem cells, was calculated at the DNA and RNA levels, respectively, to examine the differences in stem cell infiltration between risk groups.    Journal of Immunology Research

Relationship between Risk Score and Tumor Mutation
Burden. Tumor mutation burden (TMB) has been proved to be a useful biomarker for the selection of ICB in some cancer types, including lung cancer [22]. We compared TMB of two risk groups by calculating the TMB of each sample in the TCGA data set. The information about signa-ture genes mutation was obtained from cBioPortal database (http://www.cbioportal.org/).

Gene Set Enrichment Analysis.
To find out different biological processes and possible pathways between two risk groups, GSEA was conducted, including GO and KEGG assays.

Differential Expression Analysis of PRGs.
There were 47 differentially expressed PRGs, including 25 increased and 22 decreased PRGs (Figure 1(a)). The proteins expressed by these 47 PRGs were interrelated (Figure 1(b)). The results of GO showed that these differentially expressed PRGs were largely related to the biological processes related to peroxisome and lipid metabolism, such as peroxisome organization, peroxisomal   Journal of Immunology Research matrix, and fatty acid ligase activity (Figure 1(c)). KEGG showed that these genes were involved in a variety of pathways, in addition to fatty metabolism-related pathways, but also PPAR signaling pathway, glyoxylate and dicarboxylate metabolism, ferroptosis, and so on (Figure 1(d)). These results indicated that these differentially expressed genes were indeed peroxisomerelated genes and played some role in a variety of pathways. It was worth further studying their role in LUAD.

Construction of a Prognostic PRG Signature.
To screen the prognostic genes in LUAD, we performed the univariate Cox regression analysis. The results confirmed that total nine PRGs were associated with clinical outcomes of LUAD (Figure 2(a)). Six PRGs were screened to construct the signature when λ was at minimum (Figures 2(b) and 2(c)). The formula was as follows: Risk score = ½ð0:180594464 * ABCD1Þ + ð−0:14906947 * ACAT1Þ + ð− 0:080152934 * ACSL1Þ + ð0:26612956 * ACSL3Þ + ð−0:052574859 * CATÞ + ð0:431596295 * LDHAÞ . The patients with high risk score showed a shorter overall survival than those with low risk score (Figure 2(d)). The AUCs at 1, 2, and 3 years were 0.718, 0.721, and 0.707, respectively ( Figure 2(e)). Compared with AUCs of multiple clinical features, the AUC of risk score had obvious advantages (Figure 2(f)). The PCA and t-SEN analyses showed that the signature could well distinguish low-risk patients from highrisk patients (Figures 2(g) and 2(h)). The relationship between risk score and survival time and the distribution of risk score were visualized as scatter plots (Figures 2(i) and 2(j)). The univariate and multivariate Cox analyses both showed that the p value of risk score was less than 0.05, indicating that risk score was an independent indicator for patients with LUAD (Figures 2(k) and 2(l)). Our findings suggested that this new model may be used to predict the clinical outcome of LUAD patients.

Analysis of Immune Microenvironment.
It showed that the infiltration of immune cells was negatively associated with the risk score, while there was no correlation between the infiltration of stromal cell and risk score (Figures 4(a) and 4(b)). The TIDE score of low-risk group was significantly higher than high-risk group, suggesting that the high-risk group had a better response to ICB treatment (Figure 4(c)). Furthermore, we analyzed which cells were associated with the risk score. The ssGSEA results of TCGA data set and GSE72094 showed that the infiltration of DCs, B cells, mast cells, neutrophils, T helper cells, and TIL cells were higher in low-risk than high-risk group (Figures 4(d) and 4(e)). There were 18 immune checkpoints in TCGA data set, and 23 immune checkpoints in GSE72094 exhibited a dysregulated level between two risk groups (Figures 4(f) and 4(g)). And there were 16 common differentially expressed immune checkpoints in these two data sets, such as CD276, BTLA, and CD28.

Relationship between Risk
Score and Ferroptosis. The expression of ACSL4, TFRC, PTGS2, and CHAC1 was higher in high-risk group than in the low-risk group (Figures 5(a)-5(d)), indicating that ferroptosis in high-risk group is more active.

Relationship between Risk Score and Cancer Stem Cell
Infiltration. The R value was 0.31 (p = 8:2e − 11) in the RNA score, and 0.17 (p = 0:00043) in DNA score, indicating that risk score was positively associated with the degree of tumor stem cell infiltration (Figures 5(e) and 5(f)).

3.7.
Relationship between Risk Score and TMB. The TMB of high-risk group was higher than low-risk group (Figure 5(g)).
The gene mutation rate was 90.68% in the high-risk group and 86.15% in the low-risk group (Figures 5(h) and 5(i)). The gene with the highest mutation rate in the high-risk group was TNN. Among six signature genes, the mutation rates of ACSL1 were the highest, and ACSL1 was the lowest (Figure 5(j)). The common mutation type in six genes was deep deletion (Figures 5(k)-5(p)).

Gene Set Enrichment
Analysis. KEGG showed that pathways of high-risk group enriched in proteasome, cell cycle, DNA replication, and so on and pathways of low-risk group enriched in α-linolenic acid metabolism, asthma, drug metabolism cytochrome P450, and so on (Figures 6(a) and  6(b)). GO showed that biological processes of high-risk group were chromosome segregation, cell cycle G2-M phase transition, DNA replication, and so on, and low-risk group were axoneme assembly, cilium movement, ciliary plasm,  Tables S2 and S3. 3.9. Screening of Targeted Drugs. If the correlation coefficient was greater than 0, the gene was sensitive to drugs, and if the correlation coefficient was less than 0, the gene was resistant to drugs. It showed that ABCD1, ACSL1, ACSL3, and LDHA were sensitive to dabrafenib (Figure 7). ACAT1 was sensitive to parthenolide and CAT was sensitive to crizotinib (Table S4).

Validation of Signature Genes with Immunohistochemistry.
The result of immunohistochemistry showed that the expression of ABCD1, ACSL3, and LDHA was higher in LUAD tissue than normal tissue and the expression of ACAT1, ACSL1, and CAT was lower in LUAD tissues, consistent with the results of differential expression analysis (Figures 8(a)-8(f)).

Discussion
LUAD is a public health concern for its high morbidity and mortality [23]. There are numerous difficulties in designing personalized treatment plans for LUAD because of its high tumor heterogeneity and complex tumorigenic mechanism [24,25]. Therefore, robust prognostic signature is necessary. In recent years, novel prognostic biomarkers based on PRGs have been developed in hepatocellular carcinoma. Wu et al. developed a prognostic model of 9 PRGs. They also con-firmed that overall survival was much greater in the lowrisk group compared to the high-risk group [26]. However, in other types of tumors, the diagnostic and prognostic models based on PRGs have not been reported.
In this study, we analyzed TCGA data sets and identified 47 differentially expressed PRGs, including 25 increased and 22 decreased PRGs. KEGG showed that these genes were involved in a variety of pathways, in addition to fatty metabolism-related pathways, but also PPAR signaling pathway, glyoxylate and dicarboxylate metabolism, ferroptosis, and so on. These results indicated that these differentially expressed genes were indeed peroxisome-related genes and played some role in a variety of pathways. After that, a prognostic risk model based on the univariate Cox and LASSO Cox regression was developed, which included nine genes. Finally, six PRGs (ABCD1, ACAT1, ACSL1, ACSL3, CAT, and LDHA) were screened to construct the prognostic signature. Survival assays revealed that the patients with high risk score showed a shorter overall survival than those with low risk score. The multivariate Cox analyses further confirmed that risk score was an independent predictor of LUAD, highlighting the potential of the new model used as a novel prognostic biomarkers for LUAD patients. Moreover, the risk model was effective and can be applied to predict clinical outcome of LUAD patients via external validations, including GSE3141, GSE31210, and GSE72094 data sets.  of lipids and polyunsaturated fats in DCs in tumor-bearing animals, which reduced their ability to deliver antigen and trigger CD8+ T cell responses. DCs and CD8+ T lymphocytes may be involved in ferroptosis via modulating lipid and PUFA levels, according to previous studies. In addition, experimental evidence shows that T cell lipoperoxidation can impair immunity to infection by inducing ferroptosis in the cells. Consistent with the results of previous studies, CD8+ T lymphocytes have been shown to induce ferroptosis-specific lipid peroxidation and ferroptosis, therefore enhancing immunotherapy's efficacy in preclinical models. The PRG signature has emerged as an important prognostic marker and therapy option for LUAD patients, as evidenced by the correlation between risk scores and immunity. A subset of tumor cells known as cancer stem cells (CSCs) has been identified as having the ability to control self-renewal and differentiation, making it difficult to eradicate the tumor [29,30]. Breast cancer, glioma, and LUAD are among the many solid tumors where CSCs have been found [31][32][33]. In addition, we checked whether there was a link between the expression of PRG and RNAss in LUAD. Intriguingly, the outcomes were shown to be negatively related to the expression of stem cells in distinct malignancies. As a result, it was clear that the model we developed was capable of identifying CSC scores and that focusing on those genes might alter CSCs in order to prevent LUAD.
Drug resistance in LUAD treatment is on the rise, according to mounting evidence, and this is resulting in less-than-optimal therapeutic outcomes [34,35]. The increased expression of genes associated with multidrug resistance and the decrease in drug sensitivity have piqued interest. As a result, the LUAD chemotherapy prognostic model was given additional attention. Consistently, the findings exhibited that ABCD1, ACSL1, ACSL3, and LDHA were sensitive to dabrafenib. ACAT1 was sensitive to parthenolide and CAT was sensitive to crizotinib. Our findings suggested that chemotherapy for LUAD patients with lower expression of prognostic genes was beneficial.
However, our study has several limitations. First, using a single characteristic (peroxisome-related genes) to establish the predictive model was an intrinsic weakness. In fact, LUAD progress and development were influenced by a wide range of factors. Second, retrospective data sets from a public database were used to both create and validate the current model. In order to confirm its clinical value, additional data from prospective studies were still required. Third, these prognostic PRGs in LUAD necessitated more investigation to better understand their functions and processes.

Conclusion
A stable prognostic signature was established based on peroxisome-related genes to help the classification of LUAD patients. In addition, we conducted a complete review of immunology, TMB, and tumor stem cell invasion to find new therapy options. A wide range of medications targeting various genes and patient populations were examined. Novel prognostic biomarkers and therapeutic targets can be derived from the genes found in the prognostic signature.

Data Availability
The data sets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Consent
No consent was necessary.

Conflicts of Interest
The authors declare no conflicts of interest.