The list of redox-associated mRNAs (ramRNAs) was extracted from 20 gene sets of the Molecular Signatures Database (MSigDB, http://www.gsea-msigdb.org/gsea/index.jsp), an integrated web server providing annotated gene sets of distinct species. The transcriptional profile of LUAD patients was downloaded in The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) including 59 normal tissues and 535 tumor tissues along with the corresponding clinicopathological phenotypes of patients with cancer single-nucleotide variant (SNV) information. The RNA sequencing data were formatted in transcripts per kilobase million and processed by log2 transformation. Since the prognostic information of CPTAC-3 is not available, we excluded the cohort from our study. Another independent cohort (GSE36471) of Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) was chosen to validate the results externally. The rationale of selecting this cohort included the two reasons. First, the follow-up information can be obtained, which can be used for survival analysis. Second, the overlapped genes sequenced in this dataset with TCGA contained the candidates responsible for the redox signature. Statistical and bioinformatics analyses in this study were conducted using the R program (version 3.6.3 and 4.0.5).

A univariate Cox proportional hazard regression was performed to screen for survival-related ramRNAs in LUAD with the threshold of p < 0.05. The DEGs between normal and tumor tissues were identified with |log(Fold Change) (log FC)|>1.5 and p < 0.01. Then, the overlapped ramRNAs were included in K-means consensus clustering to classify LUAD patients into different subtypes based on their expression levels by using the ConsensusClusterPlus package.

Kaplan–Meier (KM) analysis was conducted to estimate the statistical differences of survival outcomes among the clustering subtypes. To assess relationships between clustering and tumor purity, the “Estimation of Stromal and Immune cells in Malignant Tumors using Expression data” (ESTIMATE) algorithm (Yoshihara et al., 2013) and Wilcoxon rank-sum test were carried out to predict the differences among clustering subtypes and immune and stromal components in LUAD. This study also investigated the distribution of clinical phenotypes including the follow-up status, age, gender, pathologic stage, T stage, N stage, and tumor status in different subtypes by the chi-squared test. The subtle biological pathways were identified between each of the two clustering subtypes by gene set variation analysis (GSVA) (Hanzelmann et al., 2013) based on “c2. cp.kegg.v7.4. symbols.gmt” and “c2. cp.reactome.v7.4. symbols” gene sets retrieved from MSigDB for Kyoto Encyclopedia of Genes and Genomes (KEGG) and RACTOME enrichment, respectively. In addition, tumor immune-infiltrating levels of 28 cell types among the clustering subtypes were evaluated by single-sample Gene Set Enrichment Analysis (ssGSEA) (Barbie et al., 2009) and the Kruskal–Wallis (KW) test.

The overlapped DEGs between each of the two subtypes with |log FC|>0.5 and p < 0.001 were included in the univariate Cox model. Significant ramRNAs were used to perform further dimension reduction through a weighted random forest and sliding window sequential forward feature selection by curtailing out-of-bag (OOB) error rate (Zhang et al., 2019) in the ranger and randomForest packages. Then, the TCGA cohort was divided into a training set and a validation set randomly at a ratio of 6:4. To validate the randomness, Pearson’s chi-squared test was performed with clinicopathological information (Bao et al., 2020). LASSO Cox regression was used to screen for optimal redox-associated signatures based on the results of feature selection with the training set. The risk score was formulated as follows: R i s k   s c o r e = α 1 α 2 … α n 1 × n × β 11 β 12 … β 1 m β 21 β 22 … β 2 m ⋮ ⋮ ⋱ ⋮ β n 1 β n 2 … β n m n × m , (1) where the row vector represents the coefficient of the n signatures, and the matrix represents the corresponding expression levels of m samples. Patients from both TCGA cohort and GEO cohort were distinguished into a high-risk and a low-risk group at the median cutoff of risk score in the training set.

A log-rank test was taken to analyze the survival differences between the high- and low-risk groups in the training set, validation set, whole TCGA cohort, and GEO cohort. Moreover, the survival outcomes were compared between the high- and low-expression groups of patients classified at the optimal cutoff of transcriptional expression in each ramRNA signature. Receiver operating characteristic (ROC) curves and AUC values of 1-, 2-, and 3-year survival probability were adopted to reveal the accuracy and sensitivity of the prediction. The relationships between the survival status and risk score were investigated with the strategy of statistical difference measured by the Wilcoxon rank-sum test. Principal component analysis (PCA) was performed to visualize the discrimination of the model. Independent prognostic factors were determined by a univariate and a multivariate Cox proportional hazard regression with the clinicopathological phenotype information in the training set, validation set, and TCGA cohort. In addition, a nomogram model was developed based on the TCGA phenotype information and validated by calibration and concordance index (C-index). Decision curve analysis (DCA) was conducted to compare the net benefit between the risk score and traditional staging system.

The top 20 mutated genes in the low- and high-risk groups were analyzed, respectively. Furthermore, the correlation between the risk and tumor mutational burden (TMB) of each sample was estimated by Spearman’s correlation coefficients. Patients with TMB information were divided into a low-TMB and a high-TMB group at the optimal cutoff for survival. Afterward, risk and TMB groups were integrated as four types of combinations: L-Risk + L-TMB, L-Risk + H-TMB, H-Risk + L-TMB, and H-Risk + H-TMB. The survival differences among the re-constructed subtypes were tested by the KM method. Stemness indices calculated by machine learning were recognized as a novel potential biomarker in anti-cancer therapies recently (Malta et al., 2018). In the present study, the relationship between the risk and RNA stemness score (RNAss) or DNA stemness score (DNAss) obtained from the UCSC Xena website (http://xena.ucsc.edu/) was examined. Moreover, GSEA was employed to screen for a significantly enriched biological process in the low- and high-risk groups.

To explore the immune differences between the low- and high-risk groups, the KW test was used to investigate the varieties of risk scores among multiple immunophenotypes, including wound healing (C1), IFN-gamma dominant (C2), inflammatory (C3), lymphocyte depleted (C4), and TGF-beta dominant (C6). The information of tumor immune-infiltrating levels estimated by TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC was acquired from the TIMER database (http://timer.cistrome.org/). The immune-infiltrating differences and PD-L1 expression levels between risk subtypes were evaluated by the Wilcoxon rank-sum test. Furthermore, immune dysfunction and exclusion of immune escape mechanisms in LUAD according to Tumor Immune Dysfunction and Exclusion (TIDE, http://tide.dfci.harvard.edu/) were explored along with risk groups.

The pRRophetic package (Geeleher et al., 2014) was adopted to predict the chemosensitivity measured by half maximal inhibitory concentration (IC₅₀) based on transcriptional expression profiles between the low- and high-risk groups.

The workflow of the present study is shown in Figure 1. A total of 153 ramRNAs were examined by univariate Cox proportional hazard regression to identify the prognostic targets (Supplementary Table S1). A total of 89 DEGs including 39 downregulated genes and 50 upregulated genes compared within normal and tumor tissues were determined (Figure 2A, Supplementary Table S2). Then, 31 aggregated prognostic differentially expressed ramRNAs were obtained (Figure 2B). Unsupervised consensus clustering with the K-means method measured by the Euclidean metric demonstrated k = 3 was an optional clustering number: 125 samples in Cluster A (low immunity), 203 samples in Cluster B (median immunity), and 185 samples in Cluster C (high immunity) (Figure 2C, Supplementary Table S3). Survival analysis indicated that patients in Cluster A had the poorest prognosis while those in Cluster B had relatively favorable outcomes (p < 0.001, Figure 2D). The chi-squared test indicated that survival status (p < 0.001), gender (p < 0.001), pathologic stage (p < 0.05), and tumor status (p < 0.01) had significant differences among the clustering subtypes (Figure 2E). As shown in Figures 2F–H, immune score, stromal score, and ESTIMATE score increased gradually and significantly from Cluster A to Cluster C (p < 0.001).

GSVA functional enrichment revealed that Cluster A was markedly enriched in pathways such as xenobiotic metabolism mediated by cytochrome P450, phenylalanine metabolism, steroid hormone biosynthesis, and drug metabolism induced by cytochrome P450 compared with Cluster B (Figure 3A); pathways including pentose and glucoronate interconversions, porphyrin metabolism, and proteasome were active in Cluster A, while several pathways of asthma, intestinal immune networks for immunoglobulin A production, and cell adhesion molecule cams were active in Cluster C (Figure 3B). DNA replication, cell cycle, and homologous recombination pathways were active in Cluster B, while Cluster C presented enrichment pathways of arachidonic acid metabolism, linoleic acid metabolism, and primary bile acid biosynthesis (Figure 3C). The results of REACTOME pathway enrichment are provided in Supplementary Figures S1A–C. Variations in tumor immune-infiltrating levels of 28 cell types were estimated by the KW test. It was demonstrated that Cluster C had significantly higher infiltrating levels of the most cell types, such as CD8+cells, regulatory T cell, activated B cell, and natural killer cells (Figure 3D). The differentially expressed strategy determined 112 DEGs between each of the two clusters was also acquired (Figure 3E, Supplementary Figure S2A; Supplementary Table S4).

A univariate Cox proportional hazard regression screened for 79 signatures with p < 0.05 out of 112 DEGs (Figure 4A). Dimension reduction was further carried out by weighted random forests and sliding window sequential forward feature selection. It was demonstrated that the OOB error reached its minimum with 49 genes (Figure 4B). The relative importance of variables is illustrated in Supplementary Figure S2B. The entire TCGA cohort was randomly divided into a training set and a validation set. Pearson’s chi-squared test was adopted to examine the clinicopathological difference between the two inner separate sets (Table 1). Based on the training set, LASSO Cox regression determined five optimal signatures (ANLN, HLA-DQA1, RHOV, TLR2, and TYMS) at the minimum partial likelihood deviance for LUAD prognosis (Figure 4C). Then, patients from both TCGA cohort and GEO cohort were classified into low- and high-risk groups at the cutoff of 2.198, which was the median value of the risk score in the training set.

Survival differences between the low- and high-risk groups in the training set (p < 0.001, Figure 4D), validation set (p = 0.005, Figure 4E), entire TCGA cohort (p < 0.001, Figure 4F), and GSE36471 (p = 0.047, Figure 4G) revealed that patients in the high-risk group experienced poor prognosis compared with the low-risk group. Simultaneously, the log-rank test of low- and high-expression groups with each optimal signature suggested that overexpressed ANLN, RHOV, and TYMS led to worse survival, while upregulation of HLA-DQA1 and TLR2 led to favorable prognosis (Figures 5A,B). ROC curves and AUC values of 1-, 2-, and 3-year survival in the training set (Figure 5C), validation set (Figure 5D), entire TCGA cohort (Figure 5E), and GSE36471 (Figure 5F) are shown. PCA demonstrated that the cumulative value of principal component 1 (PC1) and principal component 2 (PC2), which was over 70%, sufficiently explained the variance of the signature expression levels (Figures 6A–D). Hence, it was notable that patients in the different risk groups were stratified effectively. Patients in the high-risk group had higher mortality rates (training set: p < 0.001; validation set: p < 0.01; TCGA cohort: p < 0.001; GSE36471: p < 0.01, Figures 6E–H), and ANLN, RHOV, and TYMS overexpressed in the high-risk group, while the expression of HLA-DQA1 and TLR2 were reduced in the high-risk group (Supplementary Figures 3A–D). Independent prognostic analysis demonstrated that the tumor status and risk score could serve as an independent factor for predicting LUAD patients’ prognosis in the training set, validation set, and entire TCGA cohort (Table 2). Also, KM survival analysis of different clinicopathological variables showed that patients of the high-risk group had significantly poor prognosis in all subgroups including age, gender, pathological stage, and tumor status (Supplementary Figure S4). It should be noted that the survival of patients with advanced stages showed no statistical significance (Supplementary Figures 4C–E). We also found that male or advanced patients exhibited significantly higher risk scores than the other populations (Supplementary Figures 5A–D).

A nomogram model with age, gender, pathological stage, T stage, N stage, tumor status, and risk score were established (Figure 7A). The predictive accuracy of the nomogram was determined by 1-, 2-, and 3-year survival probability calibration curves (Figures 7B–D) and C-index of 0.786 (95% confidence interval (CI) = 0.765–0.808), showing favorable fitness and discrimination ability. In addition, compared with the conventional TNM stage, the risk score exhibited higher net benefits of 1-, 2-, and 3- year survival probability (Figures 7E–G), which indicated that risk scores might have optimal prognostic values.

The KW test demonstrated that Cluster C exhibited the lowest risk score (p < 0.001, Figure 8A). The distribution of different clusters among risk subtypes and the follow-up status are shown by an alluvial diagram in Figure 8B. By summarizing the mutation information including frequency and types, the mutation landscape in the low- and high-risk groups was manifested, respectively, by an oncoplot (Figures 8C,D). It was notable that MUC16, TP53, and TNN exhibited the highest frequency in the low-risk group. MUC16 and TP53 showed mostly non-sense mutation, while TP53 had the most missense mutation. Then, TMB of each sample was calculated, and Figure 8E evinced that patients in the high-risk group had higher TMB compared with the low-risk group (p < 0.001). TMB was also positively correlated with the risk score with the coefficient of 0.38 and p < 0.001 (Figure 8F). After separating samples into four novel TMB-risk subtypes, their survival differences were analyzed. It was noticeable that patients with high-TMB and low-risk showed better prognosis while patients with low-TMB and high-risk experienced the poorest survival outcomes (Figure 8G). Furthermore, stemness indices measured by DNAss (Figure 8H) and RNAss (Figure 8I) also had a significantly positive connection with the risk score (p < 0.001).

GSEA functional enrichment for biological processes indicated that the term “adenylate cyclase activating adrenergic receptor signaling pathway” was active in the low-risk group, while high risk manifested significant terms such as chromosome segregation, DNA dependent DNA replication, organelle fission, and innate immune response activating signal transduction (Figure 9A). With the information regarding immunophenotypes of TCGA cohort retrieved from UCSC Xena, it could be observed that C2 held the highest risk score, and the five immune subtypes had significantly different score levels (p < 0.001, Figure 9B). The evidence from the TIMER database demonstrated that the immune-infiltration abundance of cell types such as macrophage M0/1 (by CIBERSORT) and T cell CD4⁺ Th1/2 (by XCELL) was positively associated with the risk score while that of B cells (by TIMER), activated mast cells (by CIBERSORT), and M2 macrophages (by QUANTISEQ) were negatively associated with the risk score (Figure 9C). To explore anti-PD-L1 therapeutic effectiveness, PD-L1 expression was analyzed between the low- and high-risk groups. The Wilcoxon rank-sum test implied that PD-L1 was overexpressed in the high-risk group than the low-risk group from both TCGA cohort (p < 0.001, Figure 9D) and GEO cohort (p = 0.001, Figure 9E). In addition, immunotherapy examination according to the TIDE prediction score represented the possibility that patients with high-risk enjoying significantly lower TIDE scores (p < 0.001, Figure 9F) were insusceptible to immune evasion, which suggested that the patients were more likely to benefit from therapy with ICIs. Furthermore, the IFNG score (p < 0.01, Figure 9G) and T-cell exclusion score showed higher level in the high-risk group (p < 0.001, Figure 9H), while T-cell dysfunction score exhibited lower level in the high-risk group (p < 0.001, Figure 9I).

Chemosensitivity of drugs for anti-LUAD therapies such as cisplatin (p < 0.001), docetaxel (p < 0.001), and gemcitabine (p = 0.006) had significantly lower IC₅₀ in the high-risk group (Figures 10A–C), while erlotinib (p < 0.001) held higher IC₅₀ in the high-risk group (Figure 10D).