LUAD multi-omics data were downloaded from the TCGA database (https://portal.gdc.cancer.gov), encompassing clinical details (n=503), full transcriptome expression (FPKM format, n=576), DNA methylation (Methylation450k format, n=492), somatic mutations (mask format, n=526), and copy number variations (gistic2 format, n=516). Gene categories and names (lncRNA and mRNA) were annotated using official website files. The mature miRNA IDs from TCGA were identified using the “miRBaseVersions.db” package, while somatic mutations were analyzed with the “maftools” package. For DNA methylation data, β-values were logit-transformed, adjusted using ComBat, and then reverse logit-transformed (n=503). In addition, three other LUAD cohorts (GSE72094 (n=442), GSE68465 (n=462), and GSE31210 (n=246)) obtained from the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo) were served as external validation cohorts. Gene symbols were assigned to the microarray data probe IDs based on the GPL15048, GPL96, and GPL570 platforms, and the expression profiles were deduplicated and normalized using the robust multichip average (RMA) algorithm (18). For genes associated with several probes, the average expression value was utilized. To eliminate potential batch effects across datasets, the ComBat function from the “sva” R package was applied (19). LUAD patients with OS information and shared gene expression data across all cohorts were included in this study. The number of patients included in each cohort is as follows: TCGA (n=503), GSE72094 (n=398), GSE68465 (n=442), and GSE31210 (n=226). Supplementary Table 1 offers an overview of survival data and clinical characteristics for the patients. Notably, for multi-omics data analysis within the TCGA cohort, only patients with shared multi-omics profiles were selected as study subjects. This research flow chart is summarized in Figure 1 .

The “MOVICS” R package was created to enable thorough multi-omics clustering and visualization for cancer classification research (20). For integrative clustering analysis, we utilized TCGA multi-omics data to create two distinct data matrices, where columns correspond to shared samples (n=429) and rows denote omics characteristics. At first, a univariate Cox regression analysis was performed to pinpoint elements linked to overall survival (OS), using the pertinent data mentioned earlier. Genes with a mutation frequency exceeding 10% were classified as mutated. To precisely identify the subtypes, the cluster prediction index (CPI) and the gap statistic were utilized to assess the optimal number of clusters (21). The best cluster count for the given data was identified by choosing the number that maximized the gap statistic and CPI. Subsequently, ten different clustering techniques were utilized to categorize patients into unique subgroups. A consensus-based categorization was subsequently applied to ensure the robust identification of each subtype.

To quantify pathway activity like the EGFR network, immune-suppressed oncogenic pathways, and radiotherapy-anticipated pathways, a sample-based gene set variation analysis (GSVA) was conducted on each enriched pathway to determine patient-specific GSVA scores (22). Using the “Reconstruction of Transcriptional Regulatory Networks and Analysis of Regulons (RTN)” R package, we built transcriptional regulatory networks (regulons) that included 23 transcription factors (TFs) linked to activated or suppressed targets, along with 71 potential regulators connected to chromatin remodeling in cancer (23). Following this, the expression levels of immune checkpoint genes were analyzed among various subtypes, and the ESTIMATE algorithm was utilized to calculate stromal and immune scores for each sample (24). The MeTIL score for tumor-infiltrating lymphocytes was determined using standard procedures for DNA methylation. To assess the relative presence of immune cells, we employed single-sample gene set enrichment analysis (ssGSEA) using the “GSVA” R package. To confirm the consistency of subtypes, we initially verified the clustering outcomes with subtype-specific biomarkers in the other group. Subsequently, we evaluated the reliability of consensus clustering by contrasting it with the Nearest Template Prediction (NTP) and Partition Around Medoids (PAM) methods, conducting these analyses on both the training and test groups (20).

The “limma” R package was utilized to examine the genes with varying expression levels (DEGs) between the two MOVICS subtypes (25). Prognosis-related DEGs, identified as SPRGs, were determined in the training sets through univariate Cox regression (P<0.01). Following this, a comprehensive analysis was performed to develop a unified signature by employing 10 machine-learning techniques and 101 algorithm combinations for SPRGs. These included random survival forest (RSF), elastic network (Enet), Lasso, Ridge, stepwise Cox, CoxBoost, partial least squares regression for Cox (plsRcox), supervised principal components (SuperPC), generalized boosted regression modeling (GBM), and survival support vector machine (survival-SVM). The hyperparameters of all algorithms used the developer’s default settings. The procedure for generating the model included the following steps: (1) Prognostic biomarkers were pinpointed using Univariate Cox regression on the TCGA dataset; (2) Next, 101 different algorithm combinations were employed to create prediction models within a leave-one-out cross-validation (LOOCV) setup in the TCGA dataset; (3) These models were subsequently validated with the GEO cohorts (GSE72094, GSE68465, and GSE31210); (4) For each model, Harrell’s concordance index (C-index) was computed across all TCGA and GEO datasets, and the model with the highest average C-index was deemed the best. Comparable machine learning algorithms have been employed in prior studies (26). Parameter tuning details for the R scripts used in this study are available on GitHub (https://github.com/Zaoqu-Liu/IRLS). The detailed procedures for model selection and construction are described in the Supplementary Methods . Using the developed prognostic model, we computed a riskscore for each patient, categorizing them into high- and low-risk groups dependented on the median riskscore from whole datasets. Additionally, the performance of the riskscore was compared with that of 58 published signatures for predicting patient prognosis. To determine prognostic risk factors for LUAD, both univariate and multivariate Cox regression analyses were performed, and a predictive nomogram was created using the “rms” package in R, incorporating riskscore and clinical features.

The “maftools” R package was employed to assess the mutation status of individuals and to analyze patterns of mutually exclusive or coexisting mutations (27). GISTIC2.0 identified and pinpointed recurrent focal somatic copy number alters (CNAs) from genotype data, using a threshold of ±0.3 for amplifications and deletions (28). The TCGA-LUAD cohort’s scores for fractions of genome altered (FGA), genome gained (FGG), and genome lost (FGL) were derived from copy number fragment data. Data for tumor mutational burden (TMB) and tumor neoantigen burden (TNB) were retrieved from the UCSC database.

Charoentong et al. identified 28 signatures associated with immune cells in their research (29). The “GSVA” R package was utilized to conduct single-sample gene set enrichment analysis (ssGSEA) in order to measure enrichment scores for each gene set and sample. Additionally, immune cell infiltration abundances were assessed using three other methodologies: TIMER (30), MCP-counter (31) and ESTIMAT (31). Data on the activation stages of the seven-phase Cancer Immunity Cycle were sourced from the tracking tumor immunophenotype (TIP) database (32). Additionally, the tumor immune microenvironment (TIME) was defined by the presence of 35 immune checkpoint inhibitor (ICI) genes, as outlined in our prior research (26). TIME and metabolic signatures were also gathered from earlier research and calculated using GSVA (33). To investigate variations in 50 hallmark pathways among different risk categories, we performed GSVA enrichment analysis using the “GSVA” R package, with pathways obtained from MSigDB database (34). To confirm the crucial outcomes from the GSVA analyses, gene set enrichment analysis (GSEA) was utilized.

The SubMap technique, which operates without supervision, was utilized to evaluate the expression similarities among patients with different responses to immunotherapy. Greater similarity in expression profiles indicated a higher likelihood of similar clinical outcome (35). The TIDE framework (http://tide.dfci.harvard.edu/) was employed to predict the likelihood of tumor immune evasion by analyzing gene expression data from cancer specimens (36, 37). The IMvigor210, GSE103668, and GSE79671 datasets were analyzed to predict immunotherapy response, with the riskscore calculated for each dataset. Drug sensitivity profiles were created using data from the Cancer Therapeutics Response Portal (CTRP) and Profiling Relative Inhibition Simultaneously in Mixture (PRISM) databases, which include sensitivity information for more than 1,000 compounds (38, 39). Both databases report AUC values as indicators of drug sensitivity, where higher AUC values correspond to reduced sensitivity to specific compounds. Substances with over 20% of data missing were omitted from the inferential study (40).

Human LUAD cell lines A549, H838, and human bronchial epithelial cells BEAS-2B were sourced from Procell from Procell (Wuhan, China). H838 cells were cultured in RPMI-1640 medium supplemented with 10% fetal bovine serum (FBS), whereas A549 and BEAS-2B cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM) with 10% FBS. Every cell culture was maintained in a moist environment containing 37°C and 5% CO₂. Cells were transfected with SLC2A1 siRNAs (Hanheng, Shanghai, China) utilizing Lipofectamine 3000, adhering to the provided guidelines. Supplementary Table 2 contains the sequences for SLC2A1 siRNAs.

Total RNA was isolated from LUAD cells with the AG RNAex Pro Reagent (AG, Changsha, China) and subsequently converted into cDNA using the Reverse Transcription Kit Mix (Promega, Madison, Wisconsin, USA). cDNA amplification followed, employing SYBR Premix Ex Taq II (Promega, Wisconsin, USA). mRNA levels were quantified via qRT-PCR on the Roche LightCycler 480II, using the 2^-ΔΔCt method, with Beta-actin as the endogenous control. Primer sequences are provided in Supplementary Table 3 .

Outdo Biotech (Shanghai, China) provided fifteen sets of LUAD tumor and nearby normal tissue microarrays (HLugA030PG04). Immunohistochemical staining was conducted on theses microarrays. The microarrays were incubated overnight at 4°C with rabbit anti-SLC2A1 antibodies (ABclonal, Wuhan, China, Cat. A6982, 1:100). The evaluation of SLC2A1 expression utilized a scoring method that considered staining intensity (0 for no staining, 1 for weak, 2 for moderate, and 3 for strong) and the proportion of cells showing positive staining (<5% = 0, 5% to <25% = 1, 25% to 50% = 2, >50% to <75% = 3, >75% = 4). The ultimate score was calculated by multiplying the extent rating with the intensity rating.

For the CCK-8 test, 10,000 cells per well were seeded in triplicate into 96-well plates and incubated at 37°C with 5% CO₂ in a humidified atmosphere. The assay was conducted at 24, 48, and 72 hours post-seeding, following the manufacturer’s instructions. During the plate clone formation test, cells in the exponential growth phase were plated at a concentration of 500 cells per milliliter in 6-well dishes. The cells were cultured at 37°C in a humid atmosphere containing 95% air and 5% carbon dioxide for approximately two weeks. Colonies were then imaged and quantified using ImageJ software.

Transwell experiments were conducted to assess both cell migration and invasion. The transwell inserts were placed into a 24-well culture plate, with the insert referred to as the upper chamber and the plate as the lower compartment. Cells were broken down in a medium without serum, set to a concentration of 10⁷ cells per milliliter, and subsequently placed in the upper chamber. The bottom compartment was occupied by a base medium with 10% FBS. For invasion assays, the upper membrane was precoated with 40 µl of 8% matrigel matrix. Following a 24-hour period, the cells were treated with 4% paraformaldehyde, rinsed with PBS, stained using 0.1% crystal violet, and observed under an optical microscope at suitable magnification.

A 6-well plate was inoculated with a cell mixture at a concentration of 2 × 10⁵ cells per well. Once the cells grew more than 80% fusion, a wound was created using a 200 µl pipette tip. The cells were then washed twice with PBS and incubated in a serum-free culture medium for 24 hours. Wound closure was photographed at 24 hours using a light microscope, and cell migration ability was assessed by calculating the wound-healing rate.

R version 4.3.0 was utilized for data handling, statistical computations, and visualization. To assess relationships between continuous variables, Pearson’s correlation coefficients were used, and the Wilcoxon test was applied for differential analysis. A p-value of less than 0.05 was considered indicative of statistical significance.

We independently identified two cancer subtypes (CSs) using 10 different multi-omics ensemble clustering algorithms. The quantity of subtypes was established by combining multiple algorithms and findings from prior studies ( Supplementary Figures 1, 2 ). The clustering results were then consolidated by a consensus ensemble approach using multi-omics data, with the top 15 items for each omic presented in Figures 2A-C . Notably, CS2 exhibited a significantly higher survival probability compared to CS1 ( Figure 2D ).

To explore the potential biological functions of the different subtypes, we characterized their molecular features. Using the ssGSEA algorithm, we assessed the enrichment of various molecular signatures within the samples. Notably, CS2 had a significant enrichment of immune-suppressive cancer pathways, while CS1 showed a higher presence of pathways indicative of radiotherapy effectiveness. This suggests that CS1 may be more responsive to radiotherapy, while CS2 shows a greater sensitivity to immunosuppressive therapy ( Figure 3A ). In order to delve deeper into transcriptomic variations, we examined possible regulators linked to cancer chromatin modification and 23 LUAD-specific transcription factors (TFs) ( Figure 3B ). The strong correlation between regulator activity and CSs underscores the biological significance of these subtypes. Specifically, ERBB2, RARA, FGFR, RXRA, ERBB3, RXRB, ARM, STAT3, GATA6, and PGR were markedly activated in CS2, while FOXA1, RARB, FOXM1, and HIF1A were activated in CS1. The activity profiles of regulons associated with cancer-related chromatin changes highlight possible patterns of varied regulation among subtypes, suggesting that transcriptional networks influenced by epigenetics could be crucial distinguishing elements between these molecular subtypes. Additionally, we quantified microenvironmental cell infiltration levels and observed a significant increase in immune cell infiltration in CS1 ( Figure 3C ). Through differential expression analysis among subtypes, we choosed 50 genes that are distinctly upregulated in each subtyp act as classifiers and confirmed their consistency across several external cohorts. Using the NTP technique, each sample from the external groups was categorized into one of the determined subtypes. In line with these results, CS2 in the meta cohorts showed improved outcomes, a pattern also seen in additional external cohorts ( Figures 3D, E ). We also assessed the concordance of subtype classifications using the NTP and PAM algorithms ( Figures 3F-H ).

More than 35 possible prognostic biomarkers underwent a comprehensive analysis employing ten machine learning techniques, enabling the creation of a precise and reliable prognostic model. A total of 101 predictive models were generated, and their C-index values for the training and testing groups are shown in Figure 4A . Out of all the models, the one built with the RSF technique was deemed the best, attaining the top average C-index of 0. 724 (see Figure 4A ). Within the RSF framework, the best trees were achieved when the partial likelihood deviance hit its lowest point ( Figures 4B, C ). Seventeen genes were ultimately selected and used to construct the model. Kaplan-Meier analysis of the training group revealed that individuals identified as low-risk experienced notably improved outcomes compared to those labeled as high-risk ( Figure 4D ). This trend was also consistently seen in the testing and meta groups ( Figures 4E-H ). Supplementary Figure 3 demonstrates significant survival disparities between high- and low-risk groups across different subcategories such as age, gender, T, N, M, and stage. The outcomes of the chi-squared analysis reveal a notable link between the riskscore and clinical features like status, stage, N, and T ( Supplementary Figure 4 ). To further validate the prognostic significance of the genes included in the riskscore, we conducted a Kaplan-Meier analysis across pan-cancer datasets, which produced results largely consistent with those derived from the Cox algorithm. Additionally, These genes exhibited significant associations with the majority of tumors, highlighting their strong prognostic relevance for cancer patients ( Supplementary Figure 5 ). In the TCGA cohort, the AUC for predictions spanning 1 to 5 years varied between 0.94 and 0.99 ( Figure 5A ), indicating that the prognostic model successfully distinguished between positive and negative outcomes for LUAD patients. The signature’s predictive accuracy stayed consistent and similar across both the testing and meta cohorts ( Figures 5B-E ). Notably, the C-index of the riskscore surpassed that of nearly all clinical characteristics in both the training and testing cohorts ( Figures 5F-J ). To evaluate the predictive power of the riskscore relative to other prognostic signatures, we randomly selected 58 previously established LUAD prognostic signatures ( Supplementary Table 2 ) and calculated their C-index values. As illustrated in Figure 6 , our riskscore’s C-index surpassed that of the majority of other signatures in the training, testing, and meta cohorts. Both univariate and multivariate Cox regression analyses determined that the riskscore independently influences LUAD patient outcomes in the training and testing groups (refer to Tables 1 – 4 ).

Afterward, we utilized the GSCALite public platform (http://bioinfo.life.hust.edu.cn/web/GSCALite/) to comprehensively examine the multi-omics characteristics of the riskscore across 31 different cancer types in the TCGA dataset. This study found that in cancer types with over 12 tumor and para-cancer tissues, the genes RRM2, TK1, CCNB1, DLGAP5, CCNA2, MYBL2, and HJURP were repeatedly overexpressed across various cancer tissues ( Supplementary Figure 6A ). Furthermore, we noticed a positive relationship between mRNA expression levels and CNVs of riskscore genes in the majority of cancer types ( Supplementary Figure 6B ). Further examination of CNV frequency variations revealed notable disparities in the CNVs of riskscore genes across different cancer types ( Supplementary Figures 6C, D ). Additionally, we noticed that the methylation levels of riskscore genes varied considerably between cancerous and normal tissues in the majority of cancer samples ( Supplementary Figure 7A ). Additionally, the methylation levels of these genes were inversely correlated with their mRNA expression levels across most cancers ( Supplementary Figure 7B ). Furthermore, riskscore genes were found to activate the pan-cancer cell cycle while significantly inhibiting the hormone ER pathway ( Supplementary Figures 7C, D ).

Examining somatic mutations and CNVs uncovered notable distinctions between the low- and high-risk categories. We charted the 20 genes with the highest mutation frequencies across both risk categories, pinpointing TP53, TTN, MUC16, CSMD3, and RYR2 as the leading five based on mutation rates ( Figures 7A, B ). Figures 7C, D illustrate that the TNB and TMB levels were markedly elevated in the high-risk group relative to the low-risk group. According to Kaplan-Meier analysis, patients grouped by TNB and riskscore showed that individuals with elevated TNB and low risk had the most favorable prognosis, while those with reduced TNB and high risk had the poorest outcome ( Figure 7E ). A comparable trend was observed when integrating TMB with riskscore; patients with elevated TMB and minimal risk had the most favorable prognosis, whereas those with reduced TMB and high risk experienced the worst outcomes ( Figure 7F ). Analysis of the somatic mutation spectrum showed that the high-risk group had increased numbers of synonymous and non-synonymous mutations, along with a higher overall mutation count, in comparison to the low-risk group ( Supplementary Figures 8A–C ). Additionally, an increase in mutations was positively linked to the riskscore ( Supplementary Figures 8A–C ). Importantly, 27 genes exhibited markedly different mutation rates between the two cohorts, with a notable presence of co-mutations ( Supplementary Figures 8D, E ). Analyzing CNV between the two risk groups revealed a greater occurrence of CNV events in the high-risk group ( Figure 7G ). These findings were further corroborated by the FGA, FGG, and FGL rates ( Figure 7H ).

As previously noted, we assessed the infiltration levels of diverse immune cells in LUAD using various TIME contexture decoding algorithms. Among the high-risk cohort, TIME exhibited notably increased immune cell infiltration relative to the low-risk cohort ( Figure 8A ). In the high-risk group, the cancer immune cycle exhibited greater dynamism, highlighted by elevated activities like the release of antigens from tumor cells and the augmented recruitment of basophils, CD8 T cells, neutrophils, and natural killer (NK) cells ( Figure 8B ). Additionally, the riskscore was positively correlated with the levels of various immune checkpoints like CD274, CTLA4, and PDCD1, along with the enrichment scores of gene signatures linked to immunotherapy effectiveness ( Figure 8C ). Conversely, the riskscore showed a negative correlation with many metabolic pathways ( Figure 8D ). Figure 8E illustrates that the high-risk group showed notable enrichment in pathways associated with tumor progression, including PI3K-AKT-MTOR signaling, MYC targets, MTORC1 signaling, and the G2M checkpoint. These findings were further confirmed by GSEA analysis ( Supplementary Figure 9 ).

Considering the high rate of genomic changes and tumor mutational burden in high-risk LUAD patients, along with their active TIME and elevated levels of immune checkpoint molecules, we proposed that these individuals could show greater responsiveness to immunotherapy. In support of this, the TIDE online tool showed that the low-risk group had notably reduced TIDE scores ( Figure 9A ), and Submap analysis demonstrated that the gene expression patterns of low-risk individuals were more similar to those of melanoma patients who responded positively to immune checkpoint inhibitors (ICIs) ( Figure 9B ). These findings suggest that patients with a low riskscore may derive greater benefit from immunotherapy. Furthermore, we assessed the prognostic significance of the riskscore across three immunotherapy datasets. Importantly, there were no notable differences in survival rates between the high- and low-risk groups in the IMvigor210 study ( Figure 9C ). However, an analysis across the IMvigor210, GSE103668, and GSE79671 cohorts demonstrated that a higher riskscore was associated with increased immunotherapy response rates ( Figures 9D-F ). Additionally, we utilized a formula to pinpoint agents that might be effective for the high-risk group, leading to the identification of two agents from CTRP (paclitaxel and SB-743921) and four from PRISM (epothilone-b, litronesib, cabazitaxel, and daunorubicin). Figures 9G, H illustrate that the predicted AUC values for these agents have a statistically significant inverse relationship with the riskscore.

To evaluate the generalizability of the riskscore across various tumor types, we applied the same model to pan-cancer data. Using the defined formula, we derived the distribution of the riskscore across multiple cancers, with skin cutaneous melanoma displaying the highest riskscore ( Figure 10A ). The riskscore was also identified as a major risk factor in various cancers, such as glioma, lung squamous cell carcinoma, kidney renal clear cell carcinoma, bladder cancer, adrenocortical carcinoma, diffuse large B-cell lymphoma, pheochromocytoma and paraganglioma, kidney chromophobe, prostate adenocarcinoma, and uveal melanoma ( Figure 10A ). Furthermore, we assessed the riskscore in 31 different tumor tissues. Analysis revealed that males had significantly higher riskscores in esophageal, stomach, colon, gallbladder, ovarian, and uterine cancers. Conversely, in females, elevated riskscore were observed in skin, esophagus, stomach, colon, gallbladder, prostate, and testicular cancers ( Figure 10B ). This gender-based variation in riskscore profiles underscores the need for tailored approaches in cancer risk assessment and management.

To further elucidate the expression and functional implications of the riskscore, we initially conducted RT-qPCR analyses on nine genes in cell lines. In LUAD cells, the expression levels of SLC2A1, SFTPD, RRM2, CCNB1, CACNA2D2, SFTPB, DLGAP5, MYBL2, and HJURP were significantly higher compared to normal human lung cells, while PGC, CYP4B1, SFTPC, and CCNA2 showed marked decreases ( Figure 11 ). Due to its pronounced importance and marked upregulation, SLC2A1 was selected for in-depth experimental validation. Figures 12A, B demonstrate that immunohistochemistry (IHC) revealed a marked overexpression of SLC2A1 in LUAD tissues relative to normal tissue. To explore SLC2A1’s specific role in LUAD, we engineered cell lines with stable SLC2A1 knockdown. Post-siRNA treatment targeting SLC2A1, RT-qPCR confirmed a significant reduction in its expression in these LUAD cells relative to controls ( Figure 12C ). Functional assays, including CCK8 and colony formation tests, demonstrated that SLC2A1 knockdown markedly inhibited LUAD cell proliferation ( Figures 12D, E ). Furthermore, Transwell assays showed significant reductions in cell migration and invasion following SLC2A1 suppression ( Figure 12F ), a finding that was supported by scratch wound healing assay ( Figure 12G ).