Joint multi-omics screening of neutrophil-related biomarkers for lung adenocarcinoma

Data acquisition

Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) [18] and TCGA (https://portal.gdc.cancer.gov/) databases were accessed to source related data. The scRNA-seq dataset GSE253013 with four tumor samples and three datasets (GSE50081, GSE37745, and GSE31210) containing the expression profiles and corresponding survival information were downloaded from GEO database. All the samples included in analyses had valid survival times between 0 and 10 year(s) and a clear survival status. Corresponding platform files were downloaded to convert probes to gene names after eliminating data with ambiguous or one-to-many correspondences. When multiple probes were mapped to the same gene, their average expression value was computed. Subsequently, GSE50081, GSE37745, and GSE31210 with a sample number of 126, 90, and 225, respectively, were integrated into one dataset. The TCGA database, a source of clinical phenotype data of LUAD, was similarly processed as above, and a total of 491 tumor samples were obtained for subsequent analyses.

Single-cell clustering, dimensionality reduction, and cell annotation

Data filtering and quality control were used to extract high-quality single-cell data. Under stringent filtering criteria, each gene was required to be expressed in at least 3 cells and each cell expressed a minimum of 250 genes. Subsequently, the correlation between unique molecular identifiers (UMIs)/mRNA ratio and mitochondrial gene content was evaluated to prevent significant impacts of mitochondrial gene expression on the overall analysis of cellular transcriptomes. The proportion of mitochondrial and rRNA genes in each cell was quantified by the PercentageFeatureSet function [19] in the ‘Seurat’ R package (version 5.0.3), ensuring that cells expressed gene numbers between 500 and 5,000 genes with less than 10 % mitochondrial gene content and non-normal cells were excluded. To enhance the purity and accuracy of the dataset, the ‘DoubletFinder’ R package (version 2.0.4) [20] was applied to identify and remove double cells. Next, data from the four samples in GSE253013 were standardized using the SCTransform function in the ‘Seurat’ package (version 5.0.3) [19], and batch effects were removed using the Harmony algorithm [21].

Principal component analysis (PCA) with the RunPCA function [19] and the uniform manifold approximation and projection (UMAP) algorithm were employed to reduce dimensionality at dim=40. Based on the low-dimensional data, cell similarities were calculated by the FindNeighbors function [19], after which the cells were clustered using the FindClusters function [19] at resolution=0.3. Applying the ‘singleR’ package in R (version 1.10.0), a total of 24 cellular subpopulations were annotated using a combination of classical marker genes reported by a previous study [22].

Then, marker genes were identified using the FindAllMarkers function [19] under logic=0.35, Minpct=0.35, and a corrected p<0.05. Differences in cell composition across samples were compared by quantifying distinct cell types in each sample. Myeloid cells were subjected to the same dimensionality reduction and clustering analysis for a more detailed classification of subpopulations. Dimensionality reduction and clustering were also performed to further subdivide neutrophils for analyzing their heterogeneity.

Subpopulation analysis of neutrophils

Malignant cells were distinguished from normal cells by the ‘inferCNV’ package in R (version 1.14.2) [23], with normal immune cells and mesenchymal cells as references. Cell populations that were significantly different from the reference cells were defined as malignant. Next, the ‘monocle2’ package (version 2.26.0) was employed for pseudotime analysis to reveal the developmental trajectory and subpopulation characteristics of neutrophils, followed by developing a cell communication network using the ‘CellChat’ package (version 1.1.3) [24]. Finally, the regulatory role of transcription factors (TFs) in neutrophil subpopulation was predicted by pySCENIC (version 0.12.1) [25].

Classification of molecular subtypes

The potential influence of diverse neutrophil subsets on the prognosis of LUAD was analyzed based on the ssGSEA scores assigned to the marker genes specific to various neutrophil types. Kaplan-Meier (KM) survival curves were plotted to visualize the survival results associated with these neutrophil subsets. To further explore the distribution of these subsets in tumor tissues and their potential clinical significance, signature genes of Neu1 and Neu2 were selected to perform consensus clustering analysis using the ‘ConsensusClusterPlus’ package in R (version 1.62.0) [26] in R based on gene expression in tumor tissue samples from the TCGA dataset. The partitioning around medoids (PAM) algorithm was employed with “pearson” as the distance metric and 500 bootstraps were used to enhance the robustness of the clustering results. The optimal number of clusters was determined based on the consensus matrix and consensus cumulative distribution function (CDF). To uncover the underlying biological mechanisms of the neutrophil subsets, gene sets related to pathways were downloaded from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database via the GSEA official website (https://www.gsea-msigdb.org/gsea/index.jsp), and pathway scores for each patient were calculated using the ‘GSVA’ package in R (version 1.34.8) [27]. According to the median pathway scores, patients were divided into high-risk and low-risk groups and significantly enriched pathways between the groups were identified using the KEGG database.

Weighted correlation network analysis (WGCNA) analysis

Based on the transcriptome expression profile data of all the samples, similarity and potential outlier samples were identified by weighted correlation network analysis (WGCNA) [28]. Subsequently, a network topology analysis was conducted to develop a gene co-expression network. According to the results of the topology analysis, gene networks were constructed and modular structures within the identified network.

Screening of key genes and biomarkers

Univariate Cox regression analysis was conducted to screen genes with p<0.05. Subsequently, these genes were intersected with the hub genes identified by WGCNA to obtain key genes in LUAD. Further, the LUAD cohort from the Shengxindouyacai Database (http://www.sxdyc.com/) was accessed to develop a prognostic model using integrated machine learning methods. Leave-One-Out Cross-Validation (LOOCV) framework was employed to fit 101 predictive models, and the c-index for each model was calculated on all validation datasets to assess the prediction performance [29]. Key genes were selected based on the optimal model.

To refine the gene selection, an in-depth analysis of their functional contributions across the three GO terms was conducted using the ‘GOSemSim’ package in R (version 2.24.0) [30]. The top five genes with the greatest functional importance were determined as potential biomarkers.

Cell culture and transfection

Human lung cancer cell A549 (BNCC337696), human bronchial epithelial cells BEAS-2B (BNCC359274) and human NSCLC cells H2228 (BNCC340351) purchased from BeiNa Culture Collection (Henan, China) were respectively cultured in 90 % high-glucose Dulbecco’s modified Eagle’s medium (C3260, Solarbio Lifesciences, Beijing, China), Ham’s F-12K medium (C3340, Solarbio Lifesciences), RPMI-1640 medium (C3250, Solarbio Lifesciences) added with 10 % fetal bovine serum (S9020, Solarbio Lifesciences) at a constant temperature of 37 °C in 5 % CO₂. The cells were authenticated by short tandem repeat (STR) profiling and the mycoplasma detection results were negative.

The small intern RNA (siRNA) synthesized by Merck KGaA (Shanghai, China) were transfected into HMOX1 cells using Lipofectamine™ 2000 (11,668,027, Invitrogen, Carlsbad, CA) and the negative control with scramble target sequence (target sequence: 5′-GCTGAGTTCATGAGGAACTTT-3′) was also purchased. Immunohistochemical staining information of HMOX1 was provided by Human Protein Atlas (HPA, https://www.proteinatlas.org/).

qRT-PCR analysis

Total RNAs were extracted from BEAS-2B, A549, and H2228 cells utilizing the TRIzol™ kit (15596026CN, Invitrogen, Carlsbad, CA), and the concentration of the extracted RNA was detected using a spectrophotometer. A cDNA synthesis kit (B300902, Sangon Biotech, Shanghai, China) was employed to synthesize the corresponding cDNA. Subsequently, PCR analysis was performed using SYBR Green Abstart PCR Mix (B110031, Sangon Biotech, Shanghai, China) and the CFX396 Touch™ PCR system (1845097, Bio-Rad, CA, USA) with an initial denaturation at 95 °C for 3 min (min), followed by 40 cycles at 94 °C for 15 s (s), at 60 °C for 30 s, and at 72 °C for 1 min. Finally, the relative mRNA expressions were calculated using the 2^−ΔΔCT method, with GAPDH as an internal control reference [31]. The primers were designed using Primer Premier 6, and that of GAPDH was obtained from OriGene Technologies, Inc (Rockville, MD, USA). See Supplementary Table S1 for the primer sequences used in this study.

Wound healing assay

A549 and H2228 cells were inoculated into a 6-well plate at 1 × 10⁵ cells/well and cultivated until complete confluence. Subsequently, a controlled artificial wound on the cell monolayers was scratched using a 10 μL pipette tip, after which the cells were maintained in a serum-free medium at 37 °C. Cell images were captured at 0 and 48 h(s) under an Eclipse Ts2R microscope (Nikon, Tokyo, Japan).

Transwell assay

The upper Transwell chamber (pore size: 8.0 μm; 3422, Corning, Inc, NY, USA) coated with Matrigel matrix (M8370, Solarbio Lifesciences, Beijing, China) contained 200 μL of serum-free medium and added with A549 and H2228 cells (2 × 10⁵ cells each). Meanwhile, the corresponding lower chamber was added with 700 μL of complete medium. After culturing for 48 h, 4 % paraformaldehyde (P1110, Solarbio Lifesciences) was used to fix the cells invading the lower chamber, followed by cell staining with 0.1 % crystal violet staining solution (G1063, Solarbio Lifesciences) for 15 min. The invaded cells were observed and quantified under the same optical microscope used in the previous section.

Verification of the prognostic model and development of a nomogram

Univariate and multivariate Cox regression analyses were performed on the RiskScore and clinical characteristics, based on which a nomogram was developed. The prediction accuracy of the model was tested according to calibration curve analysis and decision curve analysis (DCA).

Mutation characteristics of the risk model

The ‘maftools’ package (version 2.14.0) in R [32] was utilized to identify the top 20 frequently mutated genes in LUAD samples. Genomic differences between high-risk and low-risk groups were compared in terms of distributions of tumor mutation burden (TMB), homologous recombination defects (HRD), number of segments, and fraction altered [33]. Additionally, copy number variation (CNV) profiles of specific genes in the risk model were analyzed.

Prediction of potential therapeutic target drugs and drug sensitivity analysis

The Cancer Cell Line Encyclopedia (CCLE) database (https://sites.broadinstitute.org/ccle) was utilized as a training set to identify key target drugs for patients in the two risk groups. The Cancer Therapeutics Response Portal (CTRP) and Profiling Relative Inhibition Simultaneously in Mixtures (PRISM) methodologies were employed to predict the sensitivity of patients in the TCGA dataset to various drugs. Effective regulatory drugs were screened under the criterion of |log₂(Fold Change)|>0.2, ensuring that the selected drugs exhibited significant differences in their therapeutic effects.

Statistical analysis

All the statistical analyses were conducted using R language (version 3.6.0) or GraphPad Prism (version 8.0). Statistical analysis was performed with unpaired Student’s t-test or one-way ANOVA followed by post hoc Tukey’s multiple comparison tests. Differences between two groups of continuous variables were compared by the Wilcoxon rank-sum test. Spearman correlation analysis was used to compute correlations. p<0.05 was considered to be statistically significant. SangerBox (http://sangerbox.com/) provided analytical assistance for this study [34].

Single-cell clustering results

After filtering and quality control of the obtained data, we successfully merged a dataset containing 54,301 high-quality cells. Correlation analysis revealed a significant correlation between UMI and mRNAs, while the UMI/mRNA ratio was not significantly correlated with the content of mitochondrial genes, indicating that mitochondrial genes did not significantly affect the results of the current study (Supplementary Figure S1A). Additionally, the violin plots of the cell quality before and after filtering and quality control (Supplementary Figure S1B) showed significant improvement in the data quality. The results of PCA dimensionality reduction are displayed in Supplementary Figure S1C. Data from the single-cell samples were normalized and standardized (Supplementary Figure S2A, B). After dimensionality reduction and clustering, a total of 24 cellular subsets were identified (Supplementary Figure S2C). In particular, neutrophils were annotated using known marker genes S100A8, S100A9, and CSF3R (Supplementary Figure S2D).

A total of 11 major clusters were classified after data normalization, batch effect removal, dimensionality reduction, clustering, and cell type annotation of these 54,301 single cells (Figure 1A and B). Specifically, these 11 clusters were B_cells, T_cells, CMP, DC, Endothelial, Epithelial, Macrophage, Mesenchymal, Monocyte, NK, and Neutrophils (the focus of this study). The expressions of the top 5 marker genes that contributed most prominently to these 11 cell clusters were analyzed (Figure 1C), and it was observed that different types of cells were specifically enriched in different samples (Figure 1D). Further dimensionality reduction and clustering of neutrophils divided 4 neutrophil subsets (Figure 1E), namely, Neu0, Neu1, Neu2, and Neu3. Functional enrichment analysis showed unique functional characteristics of these four neutrophil subsets (Figure 1F). Analysis of the enriched pathways of these four subsets revealed that the Neu0 subset was closely related to cellular energy metabolism, specifically affecting ATP-related energy pathways. The Neu1 and Neu3 subsets had similar pathways that mainly influenced immune system function, antigen processing, and presentation as well as cellular interactions. The Neu2 subset was primarily enriched in pathways closely related to intracellular protein synthesis and metabolism.

Figure 1:

The classification and functional enrichment analysis of neutrophil subsets. (A, B) Distribution of samples and cell types, resulting in the identification of 11 major clusters through annotation. (C) Exhibition of the expression patterns of the top 5 marker genes across the 11 cell types. (D) The percentage of each cell type present in each sample. (E) UMAP dimensionality reduction plot reveals four neutrophil subsets, namely Neu0, Neu1, Neu2, and Neu3. (F) Functional enrichment analysis bubble plot of the four neutrophil subsets. Note: 100646: MRC003; 101806: MRC007; 108745: MRC009; 116023: MCR001.

Neutrophil subpopulation and TF analysis

Cluster analysis successfully distinguished malignant cells (Figure 2A and B). By selecting the top 500 genes with the most significant expression alterations (Supplementary Figure S3A), five distinct functional states of neutrophils were characterized (Figure 2C). Notably, these neutrophil subpopulations played distinct roles in the temporal progression framework. Meanwhile, we developed a cell communication network to analyze the interactions between cells and their signal transmission patterns (Figure 2D and E). Finally, TF prediction analysis revealed specific TF activation and inhibition in different neutrophil subpopulations (Figure 2F), and that the expressions of these TFs in cells varied in different subpopulations (Figure 2G). Furthermore, we analyzed the patterns of cellular signals to reveal the signaling pathways and mechanisms of the signals in cellular networks (Supplementary Figure S3B). Complex interactions of the receptor-ligand pairs between neutrophil subsets and other cellular subsets were visualized in bubble diagrams (Supplementary Figure S3C). Prediction results of the amount and intensity of signals exchanged between various cell types in the cell communication network are presented in Supplementary Figures S4, 5. It was found that the Neu0 subpopulation was most inhibited by specific TFs, while the Neu3 subpopulation was most activated by specific TFs, and that the expression levels of these specific TFs were the highest in the Neu3.

Figure 2:

Analysis of neutrophil subpopulations and TFs. (A) CNV clustering heatmap. (B) UMAP dimensionality reduction plot annotated with malignant cells. (C) Pseudotime analysis of neutrophils, which were divided into five distinct states. (D) Chord diagram presenting the number of cell–cell communications. (E) Chord diagram illustrating the intensity of cell–cell communications. (F) Predicted activation and inhibition scores of TFs in the four neutrophil subpopulations. (G) Expression levels of TFs across the four neutrophil subpopulations.

Molecular subtypes and high- and low-risk groups

Based on the KM curve analysis, it was observed that the Neu1 and Neu2 subpopulations had a significant impact on the prognosis of LUAD patients in the TCGA dataset (p<0.05, Figure 3A–D), suggesting that these two neutrophil subpopulations played crucial roles in LUAD progression and patient prognosis. Consensus clustering analysis determined the optimal number of clusters to be 2, which was further supported by the CDF Delta area curve. Hence, the patients were divided into two molecular subtypes (LUAD1 and LUAD2) (Figure 3E). Further analysis showed significant prognostic differences between LUAD1 and LUAD2 (p<0.05) (Figure 3F). By scoring KEGG pathways in the two subtypes, we also discovered significant variations in multiple key pathways. Specifically, LUAD1 was primarily enriched in metabolism, cellular processes, and genetic information processing, while LUAD2 was primarily enriched in environmental information processing and organismal systems (Figure 3G). Subsequently, we conducted a differential analysis between patients in the high-risk and low-risk groups (Supplementary Figure S6) and performed KEGG enrichment analysis on the differentially expressed genes (DEGs). The results of KEGG enrichment analysis revealed that metabolism, cellular processes, and genetic information processing were more inhibited in the high-risk group but more activated in the low-risk group. Conversely, pathways associated with environmental information processing were more suppressed in the low-risk group demonstrated but more activated in the high-risk group (Figure 3H).

Figure 3:

Prognosis and enrichment analysis of the neutrophil subpopulations. (A–D) KM curves depicting patient prognosis in (A) Neu0 neutrophil subset, (B) Neu1 neutrophil subset, (C) Neu2 neutrophil subset, and (D) Neu1 neutrophil subset. (E) CDF delta area indicates relatively stable clustering results when the number of clusters is set to 2. (F) KM curves comparing the two subtypes of LUAD. (G) KEGG enrichment analysis between the two subtypes of LUAD. (H) KEGG enrichment analysis of DEGs between high-risk and low-risk groups in LUAD. *p<0.05, **p<0.01, ***p<0.001.

The results of WGCNA

Consistent clustering patterns across all the samples were confirmed, without significant detectable outliers (Figure 4A). Therefore, all the samples were included in subsequent analyses. By comparing the network topology characteristics under different soft-thresholding powers, the optimal soft-thresholding was determined to be 7 (Figure 4B and C). Next, we developed a gene network and ultimately sectioned 11 network modules, including 10 valid modules and 1 grey module (Figure 4D). Calculation of the correlations between each module and sample phenotypes revealed that the turquoise module had the strongest correlation with LUAD (Figure 4E). Finally, 381 hub genes, which were shown by the red-line box, were selected from the turquoise module under the criteria of gene significance (GS)>0.5 and module membership (MM)>0.5 (Figure 4F).

Figure 4:

The results of WGCNA. (A) Sample hierarchical clustering plot. (B) The relationship between soft threshold and scale-free topology model fit signed R². (C) The relationship between soft threshold and mean connectivity. (D) Module clustering plot and heatmap. (E) Heatmap of correlations between modules and traits. Numbers in parentheses represent p-values. (F) Correlation scatter plot between MM and GS. In the red-line box, are 381 hub genes.

Screening results for key genes

Univariate Cox analysis identified 3,487 genes closely related to the prognostic outcomes of LUAD, and 47 key genes were obtained by intersecting the 3,487 genes with the 381 hub genes screened by WGCNA analysis (Figure 5A). These 47 key genes were then subjected to a LOOCV framework to fit 101 predictive models. The optimal model was established by calculating the c-index for each model on the validation dataset. Ultimately, the Lasso+plsRcox model was determined as the optimal model as it showed the highest average c-index (0.668) in the test set (Figure 5B). Based on the Lasso+plsRcox model, we further refined the gene list to 15 key genes (CD40LG, CLEC10A, CPA3, CXorf21, CYTIP, DPEP2, FGD3, GBP1, GPR65, HMOX1, KLRB1, P2RY13, PPM1M, TM6SF1, and ZYX). The results of the univariate Cox regression analysis of these 15 genes are shown in a forest plot (Figure 5C). All 15 genes significantly influenced the prognosis of LUAD (p<0.05), with ZYX, HMOX1, and GBP1 having negative effects on the prognosis of LUAD with a Hazard Ratios (HR)>1, while the remaining 12 genes exhibited positive effects (HR<1). Using these 15 key genes, we calculated the prognostic RiskScores for each sample in the training and validation sets. Consistently, the results demonstrated that in the four cohorts of TCGA, GSE37745, GSE50081, and GSE31210, the high-risk group had noticeably poorer prognoses, while the low-risk group had relatively better prognoses (p<0.01) (Figure 5D–G).

Figure 5:

Identification and validation of key genes and biomarkers in LUAD. (A) Venn diagram intersection of Cox genes and hub genes. (B) Performance evaluation of 101 combined machine learning models on training and testing sets. (C) Forest plot of univariate Cox analysis for the 15 selected genes. (D–G) KM curve analysis comparing high- and low-risk groups across cohorts in (D) TCGA Cohort, (E) GSE37745 Cohort, (F) GSE50081 Cohort, (G) GSE31210 Cohort. (H) A bar chart displaying the percentage of functional representation for the 15 key genes across BP, CC, and MF terms in GO analysis.

To further optimize the selection of key genes, we analyzed the proportion of these 15 genes in the three GO terms and the top 5 genes (P2RY13, CYTIP, ZYX, KLRB1, and HMOX1) with the most significant functions were considered as potential biomarkers for LUAD (Figure 5H). Additionally, we examined the association between the five key neutrophil-related genes and immune cells in the TME of LUAD by applying the ssGSEA algorithm. As shown in Supplementary Figure S7, significant associations were observed between all 5 genes and most immune cells (e.g., mast cells, regulatory T cells, macrophages, etc.).

Validation of the potential biomarkers by in vitro experiments

To verify the expression profiles of P2RY13, CYTIP, ZYX, KLRB1, and HMOX1 across different cell types, we first performed quantitative experiments in BEAS-2B cells (a control cell line), A549 and H2228 cells (two lung cancer cells). It was found that the relative mRNA expression levels of P2RY13, ZYX, and KLRB1 were significantly higher in BEAS-2B cells compared to A549 and H2228 cells (Figure 6A, p<0.01). In contrast, the expressions of CYTIP and HMOX1 were markedly upregulated in A549 and H2228 cancer cells (Figure 6A, p<0.01), with HMOX1 showing a particularly high expression in A549 and H2228 cells. Hence, HMOX1 was further analyzed in subsequent experiments. We obtained immunohistochemical staining results of HMOX1 based on HPA (Supplementary Figure S8). Wound healing assay and transwell assay were conducted to investigate the role of HMOX1 in cancer cell behaviors. The results demonstrated that silencing HMOX1 in A549 and H2228 cells significantly impaired the wound healing capacity (Figure 6B, p<0.01), which was accompanied by a notable decrease in the number of migrating and invading cells (Figure 6C, p<0.01). The finding suggested that HMOX1 promoted the migration and invasion of LUAD cells.

Figure 6:

The results of in vitro validation experiments for the five key genes. (A) Relative mRNA expressions of P2RY13, CYTIP, ZYX, KLRB1, and HMOX1 in BEAS-2B, A549, and H2228 cells. (B) The impact of HMOX1 silencing on wound healing of A549 and H2228 cells. (C) Effects of HMOX1 silencing on the invasion of A549 and H2228 cells tested via transwell assay. All the experimental data of independent triplicates were expressed as mean±standard deviation. **p<0.01, ***p<0.001, ****p<0.0001.

Further improvement of the prognostic models and survival prediction by combining RiskScore with clinicopathological features

As shown by the results of univariate and multivariate Cox regression analyses, the RiskScore (HR=2.718 (95 % CI: 2.042, 3.619)) was an independent prognostic factor that was significantly more effective than other clinical features (HR=2.245 (95 % CI: 1.581, 3.187)) (Figure 7A and B). This indicated that the RiskScore had a significant influence on predicting the prognosis for LUAD patients. Further development of a nomogram also confirmed the influence of the RiskScore on the survival prediction in LUAD (Figure 7C). The calibration curve analysis showed that the predictive calibration curves of the model were close to the standard ones at the 1-, 3-, and 5-year calibration points, validating a strong predictive performance of the nomogram (Figure 7D). Furthermore, according to the DCA, both the RiskScore and the nomogram manifested significantly greater benefits than the extreme curves (Figure 7E). The analysis of the area under the curve (AUC) values also showed that both the nomogram and the RiskScore exhibited the strongest survival prediction capabilities than other clinicopathological features (Figure 7F).

Figure 7:

RiskScore and nomogram validation. (A) Univariate Cox analysis of RiskScore and clinical features. (B) Multivariate Cox analysis of RiskScore and clinical features. (C) Construction of the nomogram. (D) Calibration curves for 1-, 3-, and 5-year survival. (E) DCA of the nomogram. (F) The AUC predicted that the nomogram exhibited the most powerful capacity for survival prediction in comparison to other clinicopathological features. *p<0.05, **p<0.01, ***p<0.001.

Genetic mutation and genomic variation analysis on high- and low-risk groups

A comparison of the distribution of the top 20 frequently mutated genes in the two risk groups (Figure 8A and B) showed that TTN and TP53 were the most frequently mutated in the high-risk group, while TTN was also the most frequently mutated in the low-risk group, followed by MUC16 and RYR2. We observed that multiple gene pairs had significant mutation co-occurrence (e.g., TTN) in the high-risk group, whereas the mutations of some genes were exclusive (e.g., EGFR) in the low-risk group (Supplementary Figure S9). As for the genomic characteristics between the two risk groups, it was found that the high-risk group had significantly higher levels of HRD, number of segments, TMB, and fraction altered than the low-risk group (Figure 8C–F). These differences suggested a greater genetic instability of the tumors in the high-risk, which also explained worse treatment response and prognosis in the group. Further analysis of the CNV dataset revealed that the 15 genes in the risk model exhibited either deletions or amplifications in the samples (Figure 8G), implying that CNVs played an important role in the initiation, progression, and prognosis of LUAD by influencing gene expressions and functions.

Figure 8:

Mutation profiles and genomic characteristics in high- and low-risk groups. (A, B) Waterfall plots presenting the top 20 mutated genes in (A) high-risk and (B) low-risk groups. (C–F) Comparison of differences in (C) HRD score, (D) fraction altered, (E) number of segments, and (F) TMB between high- and low-risk groups in the TCGA cohort. (G) Distribution of patient proportions with CNV mutations in module genes within the TCGA dataset. ****p<0.0001.

Prediction of potential therapeutic targets and drug sensitivity analysis

The CTRP and PRISM methods were employed for drug sensitivity analysis of the TCGA dataset. The CTRP method identified SB-743921, paclitaxel, BI-2536, vincristine, and GSK461364 as potential therapeutic drugs for the treatment of lung cancer (Figure 9A). Meanwhile, the PRISM method discovered cabazitaxel, epothilone-b, vincristine, litronesib, rubitecan, and volasertib as the six potential drugs (Figure 9B). Vincristine emerged as a common potential drug predicted by both methods. Further analysis with the CTRP method on the differences in drug sensitivity between the two risk groups revealed that ML006 exhibited stronger therapeutic effects on the low-risk group, while the high-risk group was more sensitive to ML162 (Figure 9C). In contrast, using the PRISM method, it was found that mepivacaine had stronger therapeutic effects on the low-risk group, while XL888, 8-bromo-cGMP, and prednisolone-acetate showed more potent effects on the high-risk group (Figure 9D).

Figure 9:

Drug sensitivity prediction by CTRP and PRISM. (A) CTRP predicts drugs with significant differences in therapeutic efficacy based on their sensitivity. (B) PRISM predicts drugs exhibiting notable disparities in treatment outcomes according to their sensitivity. (C) Heatmap illustrating the fold changes in drug sensitivity between high- and low-risk groups, as well as sensitivity predictions, obtained from CTRP. (D) Heatmap depicting the fold changes in drug sensitivity between high- and low-risk groups and sensitivity predictions by PRISM.