A transcriptomic dataset was obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/, accessed on 15 January 2025), which includes RNA sequencing data along with clinical and pathological information for 530 LUAD samples and 59 normal samples (TCGA-LUAD) (Liu et al., 2019). Among these, a total of 517 LUAD samples contained complete survival data, while 222 samples had comprehensive clinical characteristics documented. An independent validation cohort, GSE26939, was retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), comprising 116 LUAD samples, 115 of which included survival information, profiled using the GPL9053 expression profiling by array platform (Song et al., 2022). Furthermore, a curated list of 86 selenium metabolism-related genes (SMRGs) was compiled based on published literature (Fu et al., 2023).

Differentially expressed genes (DEGs) that distinguish LUAD from non-tumor specimens within the TCGA-LUAD dataset were identified utilizing the R package “DESeq2” (v 1.4.2) (Love et al., 2014). For the downloaded raw data, it was recalibrated into a count matrix and filtered to retain genes exhibiting a total expression level greater than ten across all samples, applying criteria of adjusted p < 0.05 and absolute log2-fold change (FC) > 1. Visualization was conducted using the R package “ggVolcano” (v 0.0.2) (Wodrich et al., 2021), resulting in the generation of a volcano plot. Additionally, heatmap visualization was performed employing the R package “ComplexHeatmap” (v 2.14.0) (Gu et al., 2016).

To identify candidate genes related to selenium metabolism, we utilized the R package “VennDiagram” (version 1.7.3) (Chen and Boutros, 2011) to visualize the intersection analysis between DEGs and SMRGs, generating a Venn diagram. Subsequently, a functional enrichment assay was performed to investigate the biological roles of these candidate genes. Gene Ontology (GO) analysis was conducted using the enrichGO function via the “clusterProfiler” package (version 4.2.2) (Wu et al., 2021), employing the “SYMBOL” gene identifier sourced from the “org.Hs.e.g.,.db” annotation database (version 3.18.0) (Qing et al., 2022) (adjusted p < 0.05). Furthermore, for Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, we applied the enrichKEGG function within the same R package (“clusterProfiler”), focusing on human genome data with organism set as “hsa,” while maintaining identical significance thresholds of adjusted p < 0.05.

The STRING database was used to build PPI networks to investigate the functional relationships among candidate genes (confidence >0.4). Following this, Cytoscape software version 3.8.2 (Smoot et al., 2011) was utilized for visualizing these interaction networks.

Utilizing the “survival” package (version 3.7.0) (Lei et al., 2023), a univariate Cox regression analysis was conducted on LUAD samples with comprehensive survival data, using candidate genes as the basis for evaluation. Genes that met significance thresholds (HR ≠ 1, p < 0.05) and satisfied the proportional hazards (PH) assumption, assessed via the cox. zph function from the “survival” package, were retained (p > 0.05). Subsequently, we employed the “survival” package to construct a multivariable Cox regression model utilizing its built-in functions for Cox regression analysis. A bidirectional stepwise variable selection process was implemented based on the Akaike Information Criterion (AIC), employing the step function; additionally, an overall PH assumption test of the model was performed (p > 0.05). Through this methodology, prognostic genes were identified. Risk scores were then calculated by applying coefficients derived from this final model to determine risk associated with each prognostic gene in constructing our risk model. In this context, “coef” refers to the risk coefficient corresponding to each individual prognostic gene, while “expr” indicates the expression intensity of each respective prognostic gene. risk score = ∑ i = 1 n coef gene i × expr gene i

The optimal cutoff value derived from the surv_cutpoint function in the R package “survminer” (version 0.4.9) (Liu et al., 2021) was utilized to classify patients into high-risk groups (HRG) and low-risk groups (LRG). Subsequently, principal component analysis (PCA) was performed using the prcomp function to visualize the separation of risk groups. Kaplan-Meier (KM) survival plots and log-rank tests, facilitated by the R package “survminer,” were employed to compare survival outcomes between these groups. The software package “timeROC” (version 1.18.0) (Blanche et al., 2013) was leveraged to generate receiver operating characteristic (ROC) curves at three time points for evaluating model precision. Following this, expression levels of prognostic genes across both cohorts were examined utilizing the “pheatmap” package (version 1.0.12) (Gu and Hübschmann, 2022). For external validation, data from GSE26939 were analyzed with identical risk stratification methods to assess the robustness of the model.

A univariate Cox regression analysis was conducted on risk scores alongside several clinical variables, including pathological stages, gender, age, and T/N/M stages (p < 0.05, HR ≠ 1); additionally, a proportional hazards (PH) assumption test was performed with p > 0.05 as an acceptance criterion. Variables that met criteria in multivariate Cox regression analysis, specifically those with p < 0.05, and passed PH assumption testing were considered independent predictive factors for prognosis. Based on TCGA-LUAD data, we developed a nomogram to forecast 1-, 3-, and 5-year survival probabilities with the “regplot” package (version 1.1) (Sui et al., 2022). Thereafter, we generated a calibration curve employing bootstrap resampling techniques through two hundred iterations via the R package “rms” (version 6.8.1) (Sachs, 2017), aimed at assessing nomogram accuracy; proximity of slope values to one indicates enhanced predictive accuracy of our model’s estimates for patient outcomes over time.

Furthermore, we utilized the R package “ggDCA” version v (1.1) (Duo et al., 2023) to create a comprehensive nomogram forecasting 1-,3-,and 5-year survival probabilities for patients diagnosed with LUAD.

To investigate the diagnostic efficiency of risk scores in relation to clinical indicators, we analyzed correlations between risk scores and the aforementioned clinical variables. The Wilcoxon test (W-tests) was employed to compare varying risk scores and prognostic genes across different clinical subgroups (p < 0.05). Additionally, the survdiff function within the “survival” package was employed to assess survival differences among risk groups for each clinical indicator (p < 0.05).

Within the TCGA-LUAD data collection, Gene Set Enrichment Analysis (GSEA) was conducted to pinpoint biological pathways linked to HRG and LRG. Initially, differential gene expression analysis was performed between the two risk groups via the “DESeq2″ package (v 1.4.2), resulting in a ranked list of genes ordered by their log₂FC values from largest to smallest. Subsequently, we retrieved a background gene set named “c2. cp.kegg.v7.4. symbols.gmt” from the MSigDB for pathway enrichment analysis utilizing the “clusterProfiler” package (v 4.2.2). Significant pathways were defined as those exhibiting a normalized enrichment score (|NES|) > 1, p < 0.05, and false discovery rate (FDR) < 0.25.

For Gene Set Variation Analysis (GSVA), we employed the “GSVA” package (v 1.42.0) (Hänzelmann et al., 2013) to calculate pathway activity scores across samples denoted as gsva_mat. A design matrix was constructed using the model. matrix function from the “limma” package (v 3.54.0) (Ritchie et al., 2015), followed by an assessment of differential pathway activity between risk groups through linear model fitting via lmFit and empirical Bayes moderation using eBayes methods. Contrast matrices were generated with makeContrasts function, allowing us to extract significantly altered pathways with adjusted p-values <0.05 through topTable function.

The CIBERSORT method was employed to estimate the proportions of 22 distinct immune cell subtypes (Newman et al., 2015) within the tumor microenvironment of LUAD samples. Gene expression matrices underwent preprocessing and were subsequently submitted to the CIBERSORT platform, with 10 permutations and quantile normalization (QN) activated. To compare the levels of immune cell infiltration among different risk groups, W-tests were conducted (p < 0.05). The “psych” package (v 2.1.6) (Correction to Lancet Psych, 2022, 2023) was utilized to calculate Spearman correlation coefficients between differentially abundant immune cells and risk scores, as well as between these immune cells and prognostic genes (|correlation coefficient (cor)| > 0.3, p < 0.05).

The Tumor Immune Dysfunction and Exclusion (TIDE) platform was used to assess mechanisms of immune evasion and predict responses to immunotherapy. Pre-treatment gene expression profiles from LUAD samples were standardized before being submitted to TIDE for generating Dysfunction, Exclusion, and TIDE scores. Dysfunction scores reflect the functional state of effector T lymphocytes within the tumor microenvironment, while Exclusion scores indicate barriers to immune cell infiltration.

Furthermore, a total of 47 immune checkpoint genes (e.g., LAG-3, CTLA-4, PD-1) were extracted from relevant literature (Xue et al., 2022). A differential expression analysis on these genes across risk groups was performed using W-tests (p < 0.05).

The TCGA database was queried to obtain somatic mutation information for LUAD samples. Subsequently, the “Maftools” package (v 2.18.0) (Mayakonda et al., 2018) was employed to analyze the mutation annotation format (MAF) files and generate oncoplots (waterfall plots) that illustrate the mutation landscape of HRG and LRG. Entries with missing values were excluded following the annotation and filtering of variant types and classifications. The distribution of mutation types and allele-specific changes was assessed. Thereafter, using the R package “ggplot2” (v 3.4.1) (Gustavsson et al., 2022), we visualized the distribution of these mutation types across HRG and LRG populations. Tumor mutational burden (TMB), defined as the cumulative count of somatic mutations per megabase, was calculated with statistical significance set at p < 0.05. Meanwhile, the chi-square test was used to analyze the differences in the top 10 mutation classifications, different mutation types, and base changes between HRG and LRG (p < 0.05). In addition, the somaticInteractions function in the “Maftools” package (v 2.18.0) was utilized to perform co - occurrence and mutual exclusivity analysis on the mutation data of the high - risk and low - risk groups (p < 0.05).

Drug sensitivity analyses were conducted utilizing the R package “pRRophetic” (v 0.5) (Geeleher et al., 2014). This method used the Genomics of Drug Sensitivity in Cancer (GDSC, also known as CGP 2016) as the reference training set to construct a drug sensitivity prediction model. Specifically, we input the gene expression matrix (FPKM) of the TCGA-LUAD cohort into the model, with parameters set as tissueType = “lung” and dataset = “cgp 2016”. Subsequently, the pRRopheticPredict () function was applied to calculate the predicted half-maximal inhibitory concentration (IC50) of each drug for each patient. Initially, gene expression data from the TCGA database were integrated with a well-established drug sensitivity dataset known as cgp 2016; drugs classified under LUAD by TCGA with an area under curve (AUC) greater than 0.98 were selected as candidates for further analysis. The half-maximal inhibitory concentration (IC₅₀) for each drug were calculated using the pRRopheticPredict function, followed by a W-test to identify disparities in IC₅₀ values among various drugs between HRG and LRG cohorts, establishing statistical significance at p < 0.01. Additionally, we utilized the R software package “pheatmap” to visually represent relationships between first-line chemotherapeutic agents and prognostic genes in LUAD.

Drawing on TCGA-LUAD data, a W-test was first conducted to evaluate the expression of prognostic genes in LUAD samples relative to normal samples (p < 0.05). Boxplots were created with the “ggplot2″ package.

Ethical approval was obtained from the ethics committee of the First People’s Hospital of Yunnan Province (approval number: KHLL2022-KY159) before this experiment was conducted. Five LUAD tissue samples and five adjacent non-tumor tissue samples were obtained from participants at the First People’s Hospital of Yunnan Province, all of whom provided informed consent. Initially, 50 mg of tissue was extracted from each sample and homogenized with 1 mL of TRIzol (Vazyme, R401-01, China) to ensure thorough mixing and grinding. After standing on ice for 10 min, 200 µL of chloroform was added to facilitate RNA extraction from the aqueous phase. Subsequently, an equal volume of chilled isopropanol was incorporated for RNA precipitation. Following quantification, reverse transcription reactions commenced immediately thereafter. The cDNA synthesis reaction system was established strictly according to the manufacturer’s instructions for the SweScript First Strand cDNA Synthesis Kit (YEASEN, 11141ES60, China). Thereafter, qPCR amplification consisting of 40 cycles was conducted using a CFX96 real-time fluorescence quantitative PCR device (BIO-RAD, XLFZ006, United States). The primer sequences are detailed in Supplementary Table S1. In this study, glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was selected as the reference gene. In terms of experimental design, five biological replicates were set up to reduce the interference of individual sample differences, and three technical replicates were performed for each sample to ensure the repeatability of detection; after the completion of qPCR amplification, the cycle threshold (Ct values) of target genes and the reference gene were obtained simultaneously, and the melting curve was used to verify amplification specificity, while the amplification curve was used to evaluate amplification efficiency.

The expression levels of prognostic genes were evaluated utilizing the 2^−ΔΔCt method (Cheng et al., 2020), with specific steps as follows: first, the ΔCt value was calculated as the difference between the Ct value of the target gene and that of GAPDH (ΔCt = Ct_target gene - Ct_GAPDH); second, the ΔΔCt value was calculated by taking the ΔCt value of the control group as a reference and finding the difference between the ΔCt value of the experimental group and that of the control group (ΔΔCt = ΔCt_experimental group - ΔCt_control group);.finally, the relative expression level of the target gene was obtained through conversion using the 2^−ΔΔCt formula. All experimental data were subjected to statistical analysis and visualization using GraphPad Prism (v8.0) software (Chang et al., 2023): the t-test was used as the statistical method, and a p < 0.05 was considered to indicate a statistically significant difference.

All statistical analyses of data from public databases were performed using R software (version 4.2.2). The Wilcoxon test and chi-square test were used as the significance test method to compare differences between different groups, and a p < 0.05 was considered to indicate statistical significance. All network diagrams were constructed using Cytoscape software (v 3.8.2). For RT-qPCR experiments, the relative mRNA expression levels of prognostic genes were calculated using the 2^−ΔΔCt method, and the t-test was applied to compare differences in expression levels. A p < 0.05 was regarded as statistically significant.

Comparative analysis of LUAD versus normal tissues revealed a total of 14,694 DEGs (adjusted p < 0.05 and |log₂FC| > 1), comprising 11,390 upregulated genes and 3,304 downregulated genes (Figure 1A). The top 50 DEGs exhibiting the largest |log₂FC| values were visualized in a hierarchical clustering heatmap (Figure 1B), which highlighted distinct expression patterns between LUAD and normal tissue groups. Intersection analysis between DEGs and survival-related gene sets identified 40 candidate genes (Figure 1C) (Supplementary Table S2).

GO analysis of the candidate genes uncovered a total of 504 significant entries (adjusted p < 0.05), with antioxidant activity (GO:0016209, p = 1.35e-15) and cellular oxidant detoxification (GO:0098869, p = 8.44e-15) being the most enriched categories (Figure 1D) (Supplementary Table S3). Additionally, KEGG pathway analysis indicated significant enrichment across 29 pathways (adjusted p < 0.05). Notably among these pathways were the complement and coagulation cascades pathway (hsa04610, p = 4.83e-07) as well as the IL-17 signaling pathway (hsa04657, p = 0.000311) (Figure 1E) (Supplementary Table S4). Furthermore, PPI network analysis demonstrated interactions among proteins encoded by the identified candidate genes; specifically, ALB, GPX3, GPX2, IL6, APOB, APOA1, and CRP exhibited interactions with several other candidate gene products (Figure 1F). Collectively, this PPI network provides valuable insights into the dynamic interactions among proteins corresponding to these relevant genes at the protein level.

Cox regression analysis was performed on 517 LUAD samples from TCGA-LUAD to pinpoint genes significantly associated with overall survival. ALOX5AP, CRP, F2, GPX3, KMO, KYNU, SOD3, and TXNRD1 were ultimately retained in the analysis; among these, ALOX5AP, CRP, and F2 were identified as risk factors for LUAD (HR > 1) (Figure 2A). The results of the proportional hazards assumption test were illustrated using Schoenfeld residual plots; the trend lines of the residuals did not show significant differences from one another (p > 0.05) (Supplementary Figure S1A–H). Subsequently, a multivariate Cox regression model was constructed that ultimately screened four prognostic genes: F2, GPX3, KMO, and KYNU (Figure 2B). The overall proportional hazards assumption test for this model yielded p = 0.228, indicating compliance with the proportional risk assumption (Supplementary Figure S1I) (Supplementary Table S5). Risk scores were determined using the coefficients associated with these prognostic genes. As a result, the 517 LUAD patients were divided into a HRG (n = 57) and a LRG (n = 460), based on an optimal cut-off value of 1.705 (Figure 2C). Principal component analysis demonstrated distinct separation between PC1 and PC2 for HRG and LRG groups (Figure 2D).

Conversely, mortality rates were found to be significantly higher within the HRG compared to LRG individuals (Figure 2E). The Kaplan-Meier curve revealed substantial disparities in survival rates between groups; those categorized as high-risk exhibited markedly reduced survival probabilities (p < 0.001) (Figure 2F). Furthermore, the ROC assessment of the risk model validated its effectiveness in forecasting survival probabilities at 1-, 3-, and 5-year for LUAD patients (AUC >0.6) (Figure 2G). The heat map illustrating expression levels of prognostic genes indicated that F2 and KYNU had elevated expression levels in HRG compared to other groups (Figure 2H).

To assess the stability of the risk model, the same set of analyses were carried out using the dataset of the validation set. The obtained results were in line with those from the TCGA-LUAD. This consistency effectively demonstrated that the prognostic model developed in this research is capable of being utilized to predict the prognosis of patients afflicted with LUAD (Figures 2I–N).

A comprehensive univariate Cox proportional hazards regression analysis was systematically conducted to examine the impact of individual factors on prognosis. The analyses revealed a significant association between specific factors and adverse clinical outcomes (HR < 1, P < 0.05). Concurrently, a proportional hazards assumption test was performed, with the resulting p-value exceeding 0.05 indicating that the data met the proportional hazards assumption, thereby ensuring the reliability and validity of the univariate analysis results. Following this thorough evaluation, three variables, namely, risk assessment score, pathologic T stage, and pathologic N stage, were identified as potential prognostic indicators (Figure 3A).

Further multivariate Cox regression analyses were conducted to assess their independent predictive value. The multivariate analysis (p < 0.05) provided statistical evidence underscoring the importance of these variables. Additionally, another proportional hazards assumption test further validated the appropriateness of the model (p > 0.05). Collectively, these analyses robustly confirmed that risk score, pathologic T stage, and pathologic N stage are independent prognostic indicators (Figures 3B,C).

The nomogram illustrated that higher total points corresponded to poorer survival likelihood in LUAD patients (Figure 3D). A calibration plot was employed to assess the predictive accuracy of the nomogram; notably, at time points of 1-, 3-, and 5 years post-diagnosis, the slopes of these calibration plots closely aligned with reference lines, indicating an excellent fit for our model (Figure 3E). Furthermore, decision curve analysis curves demonstrated that the net benefit derived from using this nomogram exceeded zero and surpassed those associated with “all” or “none” options while generally outperforming both pathologic T/N staging and risk scores alone (Figure 3F).

Differences in risk scores were observed across various clinical features. Specifically, male patients exhibited a significantly higher risk score than female patients (p = 0.0058). Notably, box plots demonstrated significant differences in risk scores among different subgroups of clinical features (T1 and T3, stage I and stage III, N1 and N2, stage I and stage II, M0 and M1, N0 and N2, as well as between stages I and IV) (p < 0.05) (Figure 4A). Further analysis of prognostic genes within distinct clinical subgroups revealed that GPX3 and KYNU exhibited higher expression levels in male patients than in female patients; conversely, KMO and F2 showed the opposite trend. Additionally, the expression of KYNU displayed an increasing trend in cases of stage IV LUAD (Figures 4B–E). Remarkably, among LUAD patients at stages III-IV, significant survival disparities persisted between two groups. This was particularly evident among patients with pathologic N0 status as well as those with pathologic T1-2 status across all age groups (>60 years old vs. ≤ 60 years old) (p < 0.001) (Figure 4F).

To elucidate the signaling mechanisms underlying the DEGs between the two groups, GSEA was performed. A total of 50 pathway entries were identified through this analysis. Notably, in the LRG, several pathways associated with immune modulation and disorders were activated, including autoimmune thyroid disease, allograft rejection, asthma, and the intestinal immune network for IgA production. In contrast, the high-risk Group (HRG) exhibited significant enrichment in multiple pathways related to cell proliferation, metabolism, and gene regulation, specifically DNA replication, cell cycle progression, spliceosome activity, and ribosomal function (Figures 5A,B). Subsequently, GSVA was employed to further evaluate the activation or inhibition status of these pathways within both groups (Figure 5C). In particular, certain pathways in the HRG showed notable enrichment that suggests their critical role in intestinal immunomodulation. Core genes within these enriched pathways in the HRG, including MCM7, PCNA, and POLE, are recognized as essential players in DNA replication and repair processes. This observation implies a likely enhancement of cell proliferation activities within the HRG. Furthermore, genes such as CDK6, CCNB1, and CDK1 are known regulators of cell cycle progression; this finding further substantiates that the HRG promotes cellular proliferation.

A heatmap was generated to illustrate the enrichment ratios of 22 distinct immune cell subtypes across various risk levels (Figure 6A). Immune cells exhibiting diverse infiltration levels between the two groups were categorized as differentially infiltrating immune cells. Subsequently, a box-and-whisker plot was utilized to depict the percentage of these differentially infiltrating immune cells in both groups. Significant disparities were observed in the infiltration percentages of ten types of immune cells between the two groups, including CD8⁺ T cells and resting memory CD4⁺ T cells (p < 0.05) (Figure 6B).

Further analysis using Spearman’s rank correlation revealed that plasma cells and naive B cells exhibited a robust positive correlation (cor = 0.5825, p < 0.001), while resting NK cells demonstrated a strong negative correlation with activated NK cells (cor = −0.4887, p < 0.001) (Figure 6C) (Supplementary Figure S1). Notably, a heatmap illustrating the correlation between risk scores and differential immune cell populations indicated that both resting dendritic cells and resting mast cells had positive correlations with monocytes (cor = 0.33; cor = 0.30; p < 0.05). In contrast, activated mast cells showed a negative correlation with resting mast cells (cor = −0.44; p < 0.05) (Figure 6D). This finding suggests that expression levels of prognostic genes are closely associated with the presence of resting mast cells, indicating that these genes may play a significant role in this relationship. These results could enhance our understanding of disease prognosis and underlying immunological mechanisms.

TIDE analysis showed that the dysfunction indicator in the HRG was notably reduced compared with the LRG, suggesting a greater potential for immune escape. Conversely, the value of Exclusion in the HRG was markedly greater than that in the LRG, suggesting a reduced probability of immune escape (Figure 6E). The immune checkpoint mechanism plays a crucial role in evading detection and attack by T cells. We examined and compared the expression patterns of genes associated with immune checkpoints across different risk categories. Among the 47 immune checkpoint genes analyzed, 32 exhibited differential expression among distinct risk groups (p < 0.05) (Figure 6F). Notably, CD276 and TNFSF4 showed elevated expression levels in the HRG compared to those in the LRG. In contrast, the remaining 30 genes demonstrated higher expression levels within the LRG when compared to those in the HRG. The varying expressions of these genes provide a foundation for categorizing LUAD patients who may be suitable candidates for immunosuppressive treatment.

The results from the somatic mutation analysis revealed a higher mutation percentage in the HRG compared to the LRG) with rates of 98.21% versus 92.94%, respectively. Notably, KEAP1, TTN, and USH2A exhibited the highest mutation rates within the HRG at 48%, while TTN displayed the highest mutation rate in the LRG at 52%. This suggests that TTN may be one of the genes most frequently mutated in patients with LUAD (Figures 7A,B). A detailed analysis of the mutation data indicated that synonymous variants and missense variants occurred more frequently in the LRG. Furthermore, single nucleotide polymorphisms (SNPs) were also more prevalent in this group. Interestingly, G > T and C > A base changes were observed to be more common in the LRG as well. These mutations showed statistical significance between HRG and LRG. (Figures 7C–E). In addition, most of the top 20 mutated genes showed significant co-occurrence, and in particular, KRAS and TP53 exhibited mutual exclusivity (p < 0.05) (Figure 7F).

Moreover, drug sensitivity analysis identified a total of 124 compounds exhibiting differential IC50 values between these two groups (Supplementary Table S7). From this screening process, five first-line chemotherapeutic agents for LUAD treatment were highlighted: cisplatin, gemcitabine, docetaxel, etoposide, and vinorelbine. The findings demonstrated a significant correlation between KYNU and Docetaxel (p < 0.01) (Figure 7G). The respective IC50 values for these drugs were lower in the HRG cohort, indicating that individuals within this group exhibited greater sensitivity to these therapeutic agents. Additionally, further confirmation was provided by observing that IC50 values for these drugs remained consistently lower among those classified within HRG; thus reinforcing their heightened responsiveness to such treatments.

Expression validation confirmed the differential patterns of prognostic genes. Compared to normal samples, GPX3 (p < 0.0001), KYNU (p < 0.05), and KMO (p < 0.0001) were found to be downregulated in LUAD, while F2 exhibited significant upregulation (p < 0.0001) (Figure 8A). The experimental results obtained from RT-PCR demonstrated a high degree of consistency with the predictions derived from bioinformatics analyses. In comparison to normal samples, GPX3, KMO, and KYNU were significantly downregulated in LUAD, whereas F2 was upregulated in this context (Figures 8B–E). The amplification curves (obtained after 20 cycles) and melting curves (showing smooth, single peaks) demonstrated good primer specificity and favorable amplification efficiency (Supplementary Figure S2). These findings are consistent with expression profiles derived from TCGA, thereby reinforcing the biological relevance of the prognostic signature.