The transcriptome data and corresponding clinicopathological information of LUAD were collected from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) (Tomczak et al., 2015), comprising 535 LUAD samples and 59 normal tissues. To obtain oxidative stress-related genes accurately, 149 oxidative stress-related genes were contained from Gene Cards (https://www.genecards.org) (Safran et al., 2010) with a relevance score ≥ 16. In addition, the RNA-sequencing dataset with its clinical information downloaded from Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) (Barrett et al., 2013) were used for validation.

Comprehensive analysis of oxidative stress-related genes and LUAD transcriptome data from TCGA was conducted with the help of the “limma” R package. The filter criteria were set as |Fold Change|>2, padj <0.05. Finally, 34 genes met the filter conditions. Spearman’s correlation analysis was performed to identify the correlation between screened genes.

Survival R package and a univariate Cox regression analysis were used to analyze the association between 34 screened genes and overall survival, respectively. The least absolute shrinkage and selection operator (LASSO) regression was conducted to construct the optimal prognostic risk model. Four OxS-related oxidative stress genes were selected. The risk score model was calculated as follows: risk score = (−0.168∗ expression value of CYP2D6) - (0.167 ∗ expression value of FM O 3)—(0.059 ∗ expression value of CAT) + (0.306 ∗ expression value of GAPDH).

First, The Kaplan Meier plotter database (http://kmplot.com/analysis/) (Lánczky et al., 2016) was conducted to evaluate the association of four genes and overall survival. The samples were separated into high- and low-risk groups based on the median of risk score. The “survminer” and “survival” R Bioconductor packages were used to analyze the overall survival between the two groups using the Kaplan-Meier method with a log-rank test. The receiver operating characteristic curve (ROC) was calculated to analyze the predictive power of the risk model via the R “timeROC” package. Both univariate and multivariate Cox regression analysis were also performed to analyze whether the risk score and clinical factors could be independent prognostic factors for LUAD. The results were displayed in forest plots. Finally, a prognostic nomogram based on the clinical characteristics (gender, pathological stage, and age) and risk score was developed to predict the one-, two-, and three-year survival of patients with LUAD through the R “rms” package.

All cell lines were obtained from the American Type Culture Collection (ATCC). BEAS-2B, H1299, and H2030 were cultured in RPMI-1640 medium with 10% fetal bovine serum. All cells were cultured in a humid environment containing with 5% CO₂ at 37°C.

Total RNAs were obtained from cells by using TRIzol reagent (Invitrogen). Then, cDNA was synthesized by using Prime Script RT Master Mix (TaKaRa, Dalian, China), followed by quantification by TB Green Premix Ex Taq (TaKaRa) on the 7900HT Fast Real-Time PCR System (Applied Biosystems, Foster City, CA, United States ; Thermo Fisher Scientific). The relative gene expression was calculated by using the 2^−ΔΔCT method. The mRNA levels were normalized by β-actin. The primer sequences were listed in Table 1.

H1299 cells (6 × 10⁵) were seeded on each well of 6-well plates the day before transfection. According to the manufacturer’s instructions, the CAT overexpressing vector pcDNA-CAT (OE-CAT) and corresponding negative control (NC) were transfected into H1299 cells by Lipofectamine 2000. After 48 h, cells were harvested for total RNAs and following functional explorations.

To evaluate the proliferation of LUAD cells, a Cell counting kit 8 (CCK-8) (Dojindo, Japan) was adopted. In short, cancer cells (5,000 cells/200 μl) were seeded into 96-well plates, and CCK-8 reagent (20 μl) was added into each well. After incubating for 2 h, the absorbance values (450 nm) were detected. For 5-ethynyl-2′-deoxyuridine (Edu) assay, cancer cells were seeded in 96-well plates, incubated with Edu (50 µM) (Sigma) for 2 h, then fixed and permeabilized. The cells were observed analyzed through a fluorescence microscopy.

For invasion assay, we seeded cells (8 × 10⁴) in 200 μl serum-free medium into the upper chambers coated with matrigel (BD Biosciences, San Diego, CA, United States ). Then we added 500 μl medium with 10% FBS into the lower chamber. After 24 h of incubation, the cells invading the matrigel were fixed with 4% paraformaldehyde, stained with 1% crystal violet, and imaged under a light microscope. For migration assay, cells (3 × 10⁴) were seeded into the upper chamber in a similar manner without matrigel.

To quantify the proportions of infiltrating immune cells through the transcriptome data profiles from TCGA-LUAD patient cohort, the CIBERSORT algorithm was used to analyze the level of immune cell infiltration. The level of gene expression matrix of tumor-infiltrating immune cells was downloaded from the CIBERSORT database (https://cibersortx.stanford.edu/) (Chen et al., 2018). The CIBERSORT output of infiltrating immune cells proportion were accurate with p < 0.05, and only cases with p < 0.05 were eligible for further study. The R package “corrplot” was used to evaluate the correlation between immune cells and provide a visual representation.

The somatic mutation data of the LUAD patients were downloaded from TCGA. The R “maftools” package was obtained to analyze the mutation information. The Tumor Mutation Burden (TMB) considered as the mutation density of tumor genes was calculated as the transformation of total non-synonymous mutations per megabase (Schumacher et al., 2015). The difference in the TMB between the two subgroups was analyzed using a Wilcox test.

To explore the differences in the chemosensitivity between the low-risk and high-risk groups, The R “pRRophetic” package was used to predict half-maximal inhibitory concentration (IC50) values (Geeleher et al., 2014a; Geeleher et al., 2014b). Various common anticancer drugs and the cell line expression spectrum were downloaded from Genomics of Drug Sensitivity in Cancer (GDSC) (www.cancerrxgene.org/) (Yang et al., 2013). A ridge regression model was constructed to predict the IC50 of chemotherapy drugs.

R software (version 4.1.3) and GraphPad Prism 9 were used for statistical analysis. Subgroup differences were analyzed by a Wilcox test. The threshold of statistical significance was set at p < 0.05.

A total of 148 OxS genes were included to conduct a differential expression analysis in TCGA-LUAD and the adjacent tissues. The heatmap displayed 34 differentially expressed OxS genes between the normal and tumor tissues, including 13 upregulated and 21 downregulated genes (|Fold Change|>2, padj <0.05) (Figure 1A). A correlation analysis was performed to further explore the intrinsic association between these genes. As shown in Figure 1B, the level of CAT expression was most negatively correlated with GAPDH, whereas GSR expression was most positively correlated with NQ O 1. The univariate Cox regression analysis revealed that six out of thirty-four differentially expressed genes were significantly related to OS in LUAD patients. The HBG2 (p = 0.014), MRPL12 (p = 0.038) and GAPDH (p = 0.000) were considered as risky genes (HR > 1), whereas the CYP2D6 (p = 0.038), FM O 3 (p = 0.002), and CAT (p = 0.017) were considered protective genes (HR < 1) (Figure 1C). Chord diagram of Figure 1D showed that among the six genes, the correlation between GAPDH and CAT was the most significant. The level of GAPDH expression was the most likely to be negatively correlated with CAT.

The LASSO regression algorithm was conducted for the OxS-related oxidative stress genes with the optimal prognostic power. Four optimal genes, CYP2D6, FM O 3, CAT, and GAPDH, were screened as factors to establish the prognostic risk model for LUAD (Figures 2A,B). The bar chart displayed the coefficients analyzed from the LASSO regression (Figure 2D). Moreover, the risk score distribution was accurate as checked in Figure 2C. The risk score level was significantly different between the survival status groups (Figures 2E,F).

To further confirm the prognostic value of the four OxS-related genes in LUAD, the correlation of gene expression with prognosis was investigated respectively in the Kaplan-Meier plotter database. The results showed the expression of CYP2D6 (p = 0.012) and GAPDH (P < 1e-16) were negatively correlated with OS, whereas the expression of FM O 3 (p = 9.6e-14) and CAT (p = 9.4e-14) were positively correlated with the OS (Figures 3A–D). Next, the prognostic value of the model was explored in the TCGA and GEO datasets. The risk scores of all LUAD patients in the two cohorts were acquired based on the risk calculation formula. Based on the median value of risk score, the LUAD patients in the TCGA cohort were divided into a low-risk group (n = 249) and high-risk group (n = 248). The survival analysis revealed that the patients in the high-risk group exhibited a significantly worse prognosis than the low-risk group (p < 0.001) (Figure 3E). The prognostic value was further verified in the GEO cohort (GSE68465). Compared with the TCGA cohort, the survival analysis in the GSE68465 cohort showed similar results to that of the OS in the high-risk group (n = 221) were obviously shorter than in the low-risk group (n = 221) (Figure 3G). The analysis indicated that the risk score model had good power for predicting the OS of LUAD patients. The AUC at 1, 3, and 5 years achieved 0.64, 0.67, 0.61, and 0.63, 0.51, 0.63, respectively for the TCGA and GSE68465 cohort (Figures 3F,H). Despite the limited 5-year AUC in the GSE68465 cohort, the model had a good ability for predicting patient survival.

The qRT-PCR assay was employed to validate the expression levels of 4 genes in cell lines. The result showed that CAT and FM O 3 expressions were downregulated, whereas CYP2D6 was upregulated in H1299 and H2030 cell lines compared with BEAS-2B. The expression of GAPDH was upregulated in H2030 cell line but similar in H1299 compared with BEAS-2B (Figure 4A). Then we investigated the effect of the aberrant expression of CAT on LUAD cells, because which possess the strongest differential expression. A gain-of-function analyses was performed in vitro by transfecting a vector harboring CAT in H1299 cells. We found that increased CAT expression significantly inhibited the proliferation of H1299 cells (Figures 4B,C). The similar results was also observed in Edu assays (Figure 4D). Furthermore, transwell assays with or without Matrigel results demonstrated that ectopic expression of CAT remarkably suppressed the invasion and migration of H1299 cells (Figures 4E,F).

A total of 480 cases with complete clinical information were screened for the univariate and multivariate Cox regression analysis. The univariate Cox regression analysis revealed that the OS of LUAD patients was related to stage (p < 0.001, HR > 1) and risk score (p < 0.001, HR > 1) (Figure 5A). The multivariate Cox regression demonstrated that stage (p < 0.001) and risk score (p < 0.001) could be used as independent prognostic factors for predicting the prognosis of LUAD patients (Figure 5B). Nomogram is a quantitative model for predicting patient clinical outcomes. Based on gender, age, stage, and risk score, a novel prognostic nomogram was constructed to predict the survival of LUAD patients. Each variable had its normalized corresponding point. We calculated the total points of each patient by totaling the points of all variables (Figure 5C). The 1-, 2-, and 3-year survival probabilities of the patients was estimated by drafting a vertical line from the total point axis to the survival axis, which may help clinical workers make clinical decisions.

The immune microenvironment of the tumor tissue was composed of fibroblasts, stromal cells, and various immune cells, which influenced the prognosis and treatment in LUAD. To explore the correlation of risk score and level of infiltrating immune cells, CIBERSORT was used as a tool for analyzing the immune cell distribution within the TCGA-LUAD dataset. The landscape of the relative percentage of infiltrating immune cells was displayed in Figure 6A. The correlation between immune cells was analyzed in Figure 6B. All samples were divided into two groups based on the median value of the risk score. The results showed that there was a difference in various infiltrating immune cells between the two groups, including memory B cells (p < 0.0001), resting memory CD4 T cells (p < 0.0001), activated memory CD4 T cells (p < 0.0001), resting NK cells (p < 0.01), activated NK cells (p < 0.05), Monocytes (p < 0.0001), M0 macrophages (p < 0.0001), M1 macrophages (p < 0.05), resting dendritic cells (p < 0.0001), and resting mast cells (p < 0.0001) (Figure 6C). In contrast, there was no statistical significance in the number of naive B cells, plasma cells, CD8 T cells, follicular helper T cells, regulatory T cells, gamma delta T cells, M2 macrophages, activated dendritic cells, activated mast cells, eosinophils and neutrophils between the two groups. The proportion of each immune cell type was provided in Supplementary Table S2.

Moreover, the somatic mutation information of TCGA-LUAD patients was used to explore the correlation between the risk score and mutational landscape. After excluding samples with incomplete mutation information, 451 LUAD patients from the TCGA database were incorporated into the analysis. The samples were divided into a low- and high-risk group based on the median risk score. Detailed mutation information of each gene in all samples was exhibited by a waterfall plot. The top 20 mutated genes were displayed in the plot, whereas the different colors represent different mutation types. Missense mutations were the most common mutation type among the 20 genes shown. Moreover, the mutation frequency significantly differed between the low- and high-risk groups in each gene cohort. We found that more mutation events occurred in the high-risk group than in the low-risk group (Figure 7A). The tumor mutation burden (TMB) was visualized in Figure 7B. The TMB in the high-risk group was appreciably higher compared with the low-risk group, indicating the superior effect of immunotherapy. (Figure 7C)

Immune checkpoint molecules are regulatory molecules that play a role in the immune system and play a critical role in maintaining self-tolerance, preventing an autoimmune response, and controlling the timing and intensity of the immune response. The difference in the levels of immune checkpoint expression in the low- and high-risk groups was compared according to the median value. As shown in Figure 8A, the level of TIGIT (p = 0.05), CTLA4 (p < 0.05), HAVCR2 (p < 0.05), IL-10 (p < 0.0001), and TGFB1 (p < 0.01) differed between the two groups, whereas no statistical difference was observed in the level of LAG3, PDCD1, and CD274 between the two groups. The results suggested that the risk score may be correlated with some immune checkpoints. It may be helpful in the immunotherapy of LUAD. The distribution of immune checkpoint expression levels was also visualized by heatmap in Figure 8B.

We identified some chemotherapeutic drugs and immunosuppressors via pRRophetic according to the OxS-related risk score. Lenalidomide (p < 0.001), nilotinib (p < 0.001), shikonin (p < 0.001), methotrexate (p < 0.001), displayed lower sensitivity in the high-risk group, whereas epothilone (p < 0.001), thapsigargin (p < 0.001), rapamycin (p < 0.001), vinblastine (p < 0.001), elesclomol (p < 0.001), docetaxel (p < 0.001), parthenolide (p < 0.001), and Paclitaxel had higher sensitivity in the high-risk group (Figure 9).