Integration of the bulk transcriptome and single-cell transcriptome reveals efferocytosis features in lung adenocarcinoma prognosis and immunotherapy by combining deep learning

Abstract

Background

Efferocytosis (ER) refers to the process of phagocytic clearance of programmed dead cells, and studies have shown that it is closely related to tumor immune escape.

Methods

This study was based on a comprehensive analysis of TCGA, GEO and CTRP databases. ER-related genes were collected from previous literature, univariate Cox regression was performed and consistent clustering was performed to categorize lung adenocarcinoma (LUAD) patients into two subgroups. Lasso regression and multivariate Cox regression analyses were used to construct ER-related prognostic features, and multiple immune infiltration algorithms were used to assess the correlation between the extracellular burial-related risk score (ERGRS) and tumor microenvironment (TME). And the key gene HAVCR1 was identified by deep learning, etc. Finally, pan-cancer analysis of the key genes was performed and in vitro experiments were conducted to verify the promotional effect of HAVCR1 on LUAD progression.

Results

A total of 33 ER-related genes associated with the prognosis of LUAD were identified, and the prognostic signature of ERGRS was successfully constructed to predict the overall survival (OS) and treatment response of LUAD patients. The high-risk group was highly enriched in some oncogenic pathways, while the low-ERGRS group was highly enriched in some immune-related pathways. In addition, the high ERGRS group had higher TMB, TNB and TIDE scores and lower immune scores. The low-risk group had better immunotherapeutic response and less likelihood of immune escape. Drug sensitivity analysis revealed that BRD-K92856060, monensin and hexaminolevulinate may be potential therapeutic agents for the high-risk group. And ERGRS was validated in several cohorts. In addition, HAVCR1 is one of the key genes, and knockdown of HAVCR1 in vitro significantly reduced the proliferation, migration and invasion ability of lung adenocarcinoma cells.

Conclusion

Our study developed a novel prognostic signature of efferocytosis-related genes. This prognostic signature accurately predicted survival prognosis as well as treatment outcome in LUAD patients and explored the role of HAVCR1 in lung adenocarcinoma progression.

Introduction

Lung cancer mortality worldwide is at a high level among cancers [1], of which lung adenocarcinoma (LUAD) is the most common histologic subtype [2, 3]. Due to the insidious nature of the disease and lack of specificity, most lung cancer patients are in advanced stages at the time of diagnosis. Therefore, non-surgical radiotherapy and chemotherapy have become the mainstay of treatment for patients with advanced lung cancer [4]. Many new therapies have also been discovered that can be used to treat LUAD, such as immunotherapy [5]. Despite advances in timely diagnosis and drug formulation, clinical outcomes for patients with advanced LUAD are poor, with five-year survival rates of 4–17% for patients with localized LUAD [6]. Due to lung adenocarcinoma has no obvious symptoms in its early stages, the majority of patients are not diagnosed until the disease has progressed to advanced stages or metastasized [7, 8]. Therefore, the identification of novel molecular markers and prognostic indicators for LUAD is essential to improve treatment outcomes and minimize disease burden [9].

Efferocytosis, the process of phagocytosis and elimination of apoptotic cells, is a key biological process [10]. During cancer progression, tumor cells undergo various forms of cell death, such as apoptosis, necrotic apoptosis, ferroptosis, and pyroptosis, as a result of mutations, hypoxia, and radiotherapy and chemotherapy [11]. Apoptotic tumor cells are rapidly and efficiently detected and removed by professional and non-professional phagocytes through efferocytosis, preventing secondary necrosis [12]. Most of the genes that promote efferocytosis are involved in tumor progression and metastasis and are frequently overexpressed in a variety of cancers, including lung, breast, and leukemia [13, 14]. The process of efferocytosis contributes to the formation of an immunosuppressive microenvironment, thereby assisting in tumor immune escape [15]. The critical role of efferocytosis action in cancer development and progression is attributed to its effects on tumor cell growth, metastasis, epithelial mesenchymal transition (EMT), and angiogenesis [16]. Although conventional oncology treatments such as chemotherapy and radiotherapy trigger apoptosis and efferocytosis effects, they also lead to tumor inflammation and limit antitumor immunity [17]. Studies have shown that blocking efferocytosis effects alone does not completely inhibit the production of tumor immunosuppressive cells and mediators [18]. In addition, the efferocytosis action process may act as an immune checkpoint similar to PD-1/PD-L1, which could be a target for therapeutic intervention [19]. However, studies examining the effect of efferocytosis action on LUAD progression are still lacking.

The aim of this study was to characterize the prognostic potential and tumor microenvironment (TME) landscape of LUAD based on efferocytosis (ER) features obtained by integrating scRNA-seq and bulk transcriptomic data, which would provide important theoretical guidance for prognostic stratification and precision treatment of LUAD. A prognostic scoring system for patients was constructed based on the obtained efferocytosis features through deep learning validation. We investigated the potential mechanisms underlying the different prognoses of patients with different ER risks in terms of coherent clustering, immune infiltration, model comparison, immunotherapy prediction, and drug sensitivity analysis from single-cell transcriptomes and bulk transcriptomes. Nomogram models were also developed to accurately predict the survival of LUAD patients, which helps to tailor the treatment plan for patients according to their ER risk class. The workflow is shown in Fig. 1.

Fig. 1

Work flow

Methods

Data collection and processing

Clinical information, transcriptome expression, and data on copy number variations (CNVs) and single nucleotide variants (SNVs) of TCGA-LUAD patients were obtained from the TCGA database (https://portal.gdc.cancer.gov/) [20]. The pan-cancer cohort contains data on more than 10,000 patients with 33 different cancers, which were also obtained from the TCGA database. From the GEO database (https://www.ncbi.nlm.nih.gov/geo/) [21] five GEO datasets of lung cancer patients were obtained from the GEO database: GSE68465, GSE31210, GSE37745, GSE50081, and GSE30219. We also collected transcriptomic and clinical information on cancer patients in the IMvigor210 cohort who received anti-PD-L1 therapy. The source of the data is as follows: http://research-pub.gene.com/IMvigor210CoreBiologies [22]. Genes associated with efferocytosis were based on previously published studies [23, 24] (Table S1). TIDE scores for LUAD patients predicting ICB response were calculated by TIDE (http://tide.dfci.harvard.edu) [25].

Single-cell RNA-seq analysis data collection and processing

In our study, we incorporated 10 samples from the dataset by Philip Bischoff and colleagues [26], which comprised 5 samples from individuals with normal lung tissue and 5 from patients with Lung Adenocarcinoma (LUAD). We leveraged the Seurat R package for the analysis of single-cell RNA sequencing data. To mitigate batch effects, we employed the Harmony software package. Clustering of cells was accomplished using Seurat's "FindClusters" and "FindNeighbors" functions, with the results visualized through t-SNE mapping. Cellular annotation was carried out by referencing the signature genes associated with distinct cell types. The Seurat package's "AddModuleScore" function was utilized to assess the activity levels of a particular gene set within individual cells. For the comparison of differential gene expression between the two cohorts, we relied on Seurat's "FindMarkers" function. The Wilcoxon test was applied to determine the statistical significance of differentially expressed genes (DEGs), with a threshold of p.adj < 0.05; all other parameters were left at their default settings.

Consistency clustering

In this study, we identified 33 efferocytosis-related genes with prognostic relevance. We used the "ConsensusClusterPlus" package for consistency clustering to categorize lung adenocarcinoma patients into two subgroups [27]. In the two independent cohorts, there were significant differences in clinical and prognostic characteristics between the two subgroups of LUAD patients.

Functional enrichment analysis

To elucidate the potential biological pathways associated with efferocytosis-related genes (ERGs) and efferocytosis models (ERGRS), we used the R package ClusterProfiler [28] performed GO, KEGG and GSEA enrichment analysis. In addition, we utilized the GSVA software package was performed to further reveal the relevant biological functions [29].

Quantifying immune infiltration in LUAD patients

The ESTIMATE algorithm was utilized to estimate infiltration abundance between different groups using transcript expression data from the TCGA database. We used six different algorithms to assess immune cell infiltration in LUAD patients in the TCGA dataset. These included TIMER, CIBERSORT, MCP-counter, EPIC, xCell and quanttiseq, all using the R package “IOBR” [30] calculation; in addition, the level of immune cell infiltration in LUAD TME was assessed by the single sample gene set enrichment analysis (ssGSEA) algorithm [31].

Construction of the ER risk scoring system

First, efferocytosis-related genes with prognostic value were obtained by univariate Cox analysis. Subsequently, prognostic hub efferocytosis genes were further screened by the least absolute shrinkage and selection operator (LASSO). Finally, the efferocytosis gene signature consisting of seven genes was obtained by stepwise multivariate Cox regression analysis. Based on the efferocytosis gene signature, the ER risk scoring system ERGRS was established: − 0.403*ENTPD1 + 0.248*HAVCR1 + 0.143*IL1A + − 0.221*MERTK + − 0.346*P2RX1 + 0.305*PANX1 + 0.331*SPHK1. The model performance is compared with that of published papers.

ER risk prediction by a deep learning model

The random forest model determines the optimal number of variables by calculating the average error rate of candidate central genes [32]. We then calculated the error rate for each of the 1 to 500 trees and determined the optimal number of trees based on the lowest error rate. In determining the above parameters, a random stand tree model is established. Finally, the feature importance score of each candidate central gene was determined, and 5 genes with an importance value greater than 0.25 were selected. We used the random forest (RSF) algorithm to screen five gene features and took the intersection with the genes screened by the LASSO algorithm to obtain a total of three genes. A deep neural network (DNN) is also called a Multi-Layer perceptron (MLP), which consists of an input layer, hidden layer, and output layer. In the study, we choose python to build the DNN model and TensorFlow as the back-end tool of Keras. The first layer in front is the input layer, the middle is the hidden layer, and the last layer is the output layer. We randomly divided the TCGA-LUAD samples into a training set (80%) and a test set (20%), and DNN classification model to investigate the obtained three efferocytosis gene features' ability to predict ER risk in LUAD patients. The performance of the model in the training set was evaluated by accuracy, sensitivity, specificity, and area under the ROC curve (AUC) [33].

Assessment of response to neoadjuvant therapy in patients with different ER risks

Tumor mutational load (TMB), a potential biomarker of immunotherapy response, was calculated based on somatic non-synonymous mutations. Tumor immunophenotype scores (IPS) were obtained from The Cancer Immunome Atlas (TCIA) and analyzed. Predictions of immunotherapy efficacy were validated using risk modeling using three immunotherapy cohorts: the GSE78220, an anti-PD-L1 antibody (IMvigor210 cohort) [34], anti-CTLA4 and anti-PD1 therapy (GSE91061) [35]. Responses and gains for the TCGA cohort were based on the TIDE algorithm. Immunotherapy responses were then predicted using the SubMap online tool [36]. Drug sensitivity data were obtained from the CTRP (https://portals.broadinstitute.org/ctrp/) and PRISM (https://www.theprismlab.org/knowledge-base/) databases.

Cell culture and transfection

Human lung adenocarcinoma cell lines A549 was mainly purchased from the cell bank of the Chinese Academy of Sciences (Shanghai, China). We used A549 cells for in vitro culture experiments in DMEM medium (Gibco, ThermoFisher Scientific, United States) supplemented with 10% fetal bovine serum, 1% penicillin and streptomycin (Gibco). Small interfering RNA (siRNA) targeting HAVCR1 and interfering RNA control were purchased from Gemma Genetics (Shanghai, China). For transient transfection, A549 cells were transfected with siRNA using a transfection reagent (Lipofectamine 2000) for 12h, followed by functional assays and subsequent experiments. Total RNA was extracted using an RNA extraction kit (Vazyme, China) according to the manufacturer’s instructions. cDNA was synthesized for real-time PCR using SYBR Green qPCR mix (Vazyme, China).

Immunohistochemistry (IHC)

A total of 6 pairs of paraffin-embedded lung adenocarcinoma (LUAD) and adjacent tissues were collected to assess HAVCR1 expression via immunohistochemistry (IHC). The study was approved by the Ethics Committee of the First Affiliated Hospital of Bengbu Medical University. Paraffin-embedded tissue sections were deparaffinized and incubated with HAVCR1 primary antibody (1:100; Novus). After secondary antibody incubation, sections were processed with an IHC staining kit (MaxVision Biotechnology, China) and visualized under a light microscope (Olympus, Tokyo). Expression was evaluated based on staining intensity and the percentage of positive cells. The IHC scoring criteria were as follows: positive cells were identified by brown-yellow granules in the cytoplasm and nucleus. Staining intensity was scored as 0 for blue, 1 for yellow–brown, and 2 for brown. The percentage of positive cells was scored as 1 for 0–10%, 2 for 11–30%, 3 for 31–50%, and 4 for > 50%. A product of the two scores > 3 was considered positive for HAVCR1. Sections were deparaffinized with a xylene-alcohol gradient. After rinsing in water, endogenous peroxidase activity was quenched with hydrogen peroxide, followed by antigen retrieval in a pressure cooker with EDTA, cooling at room temperature, and rinsing with PBS. Sections were then blocked with serum, incubated with primary and secondary antibodies, visualized with DAB, and counterstained with hematoxylin.

Transwell assay migration invasion and wound healing assays

Transwell assay migration and wound healing assay A549 cells were transfected with HAVCR1 siRNA for 24 h and cultured in 24-well culture plates with 8 mm pore membrane inserts to measure cell migration and invasion capacity. 4 × 10⁴ cells were inoculated in the upper chamber of a transwell assay with 200 μL of serum-free medium, and 800 μL of medium containing 10% FBS was added to the lower chamber. After 48 h of incubation, cells migrating across the membrane were fixed with paraformaldehyde, stained with 1% crystal violet and counted under a light microscope (200 ×). In addition, A549 cells were cultured in 6-well plates and scraped with a 200 μL pipette tip. Cells were cultured in DMEM medium without FBS.

Proliferation and clone formation experiments

Cell proliferation and colony formation assays A549 and H1299 cells were cultured in 96-well plates (3000 cells/well) 24 h after transfection with HAVCR1 siRNA. The proliferative capacity of the treated cells was assayed at 4, 24, 48 and 72 h. 10% Cell Counting Kit-8 (CCK8) reagent (Bio-sharp, Hefei, China) was added to each plate according to the kit instructions, and the OD450 values were analyzed by an enzyme marker (BioTek, United States). Regarding colony formation experiments, 3000 cells were inoculated in cell culture plates and allowed to grow until visible colonies were formed. Then we fixed the clones with paraformaldehyde for 15 min, stained the clones with 1% crystal violet for 20 min and counted the number of clones (> 50 cells).

Statistical analysis

Data processing and statistical analyses were performed using Linux (Ubuntu 18.04.5), python (version 3.1.1), and R (version 4.1.3) programs. Kaplan–Meier survival curves were plotted by the survminer R package (version 0.4.9) for OS analysis. Wilcoxon signed rank test was performed for statistical comparison between the two groups. The chi-square test was used to analyze significant differences between the two categorical variables. p < 0.05 indicated statistical significance.

Results

Identification of 33 hub genes for efferocytosis

167 efferocytosis-related genes (ERGs) were obtained from the literature, and we used one-way Cox regression analysis to obtain 33 prognostic genes associated with OS. We have shown a forest plot of ten of these genes. (Fig. 2A). ER-related gene mutations were found in 9.09% (56/616) of LUAD patients. Among them, HAVCR1 was the gene with the highest mutation rate, followed by CX3CR1 and TLR2 (Fig. 2B). Meanwhile, we observed varying degrees of DNA copy number variation in ERGs, among which, CD300LF, SPHK1, CD300LB, and PECAM1 showed extensive CNV amplification, and CNV depletion was present in some ERGs (Fig. 2C). Figure 2D demonstrated the locations of CNV alterations in ER-related genes on chromosomes; finally, we analyzed the differential expression patterns of these 33 genes between normal and tumor tissues in TCGA, and found that most of these genes showed significant differential expression (Supplementary Fig. 1A).

Fig. 2

Identification of 33 hub genes for efferocytosis and Single-cell validation of Efferocytosis features. A Forest plot of 10 genes. B Distribution of ERGs and mutation frequency in the TCGA-LUAD cohort. C CNV alteration frequency of ERGs in LUAD, the height of the bar represents the mutation frequency. D Location of CNV alterations in ERGs on chromosomes. E t-SNE plot showing the cell types identified by marker genes. F Efferocytosis score (ER) of each cell; G Heatmap showing the 5 most important marker genes in each cell cluster. H Violin plot showing Efferocytosis score for each cell; I GSEA in the high ER group enrichment analysis, including apoptosis, cell charring, iron death and necrotic apoptosis

Single-cell validation of Efferocytosis features

Single-cell RNA sequencing data were collected from LUAD patients using Bischoff et al. Using marker genes for different cell types, we labeled cells into seven major clusters, namely CD8 T cells, CD4 T cells, macrophages, mast cells, endothelial cells, epithelial cells, and B cells (Fig. 2E). To quantify the activity of efferocytosis (ER) in different cell types, we used the "AddModuleScore" function in the Seurat software package to calculate the expression levels of 33 ER-related genes in all cells (Fig. 2F). Among the seven cell types, we observed significantly elevated ER activity in macrophages and endothelial cells (Fig. 2H). The heightened ER activity in macrophages and endothelial cells is indicative of an active role in the clearance of cellular debris and apoptotic bodies, which is a hallmark of efferocytosis. This process is crucial in the tumor microenvironment as it can influence the immune response and the progression of LUAD. The heatmap demonstrated the expression of 5 marker genes in each cell population (Fig. 2G). GSEA revealed the upregulation of 4 cell death-related pathways (apoptosis, cellular pyroptosis, ferroptosis, and necrotic apoptosis) in the high-scoring group, suggesting that multiple forms of cell death may occur within lung adenocarcinoma in this group. (Fig. 2I).

Our observations suggest that the modulation of ER activity in macrophages and endothelial cells could serve as a therapeutic target to enhance the efficacy of immunotherapy. This could be particularly relevant in LUAD, where understanding the interplay between efferocytosis and immune cell function is key to developing more effective treatment strategies.

Consistent clustering classifies LUAD patients into two types

By unsupervised clustering of 33 ERGs, we sought to identify unprecedented subtypes associated with efferocytosis in LUAD. We found that k = 2 was optimal for categorizing the entire cohort into clusters C1 (n = 249) and C2 (n = 251) (Fig. 3A, B). Principal component (PCA) analysis of LUAD patients was well distributed in both clusters (Fig. 3C). Kaplan–Meier survival analysis showed a significantly better prognosis for C2 compared to C1 (p = 0.0072) (Fig. 3D). Additionally we also obtained the same results by GSE31210 and GSE50081 validation sets (Supplementary Fig. 1B–G), where we compared the expression of ERGs in both subtypes. Some ERGs were highly expressed in C2, such as NT5E, CX3CL1, S100B and TLR2 (Fig. 3E). To investigate the biological characteristics of two ER subtypes, we conducted GSVA analysis. C1 subtype showed significant enrichment in oncogenic pathways like MYC targets and G2M checkpoint, and in processes like glycolysis and DNA repair. C2 subtype was enriched in immune pathways, including IL2/STAT5 and IL6/JAK/STAT3 signaling (Fig. 3F). In the GSEA analysis based on the KEGG gene set, we observed a significant enrichment of the C1 subgroup into cell cycle- and oncogenesis-related pathways, and further enrichment in the C2 subgroup into immune-related pathways (Fig. 3G–H). Finally, we used univariate and multivariate Cox regression analyses to confirm cluster grouping as an independent prognostic factor for LUAD (Fig. 3I, J). Functional enrichment analysis of DEGs revealed that ER is associated with the assembly of the MHC protein complex and has molecular functions in binding amide, peptides, and the MHC protein complex itself. KEGG analysis showed that these genes are involved in pathways like Phagosome and the Intestinal immune network for IgA production (Supplementary Fig. 1H, I).

Fig. 3

Consistent clustering classifies LUAD patients into two types. A LUAD patients were clustered into two molecular clusters when k = 2 based on 33 ERGs; B indicates that clustering results were best when k = 2; C PCA analysis showed significant differences between the two subtypes; D Kaplan–Meier survival analysis indicates that C2 has a better prognosis. E Differences in expression levels of ERGs between the two subtypes, red markers represent differentially expressed genes, p < 0.05; F GSVA demonstrates the HALLMARK pathway of different subtypes, red represents promotion and blue represents inhibition; G GSEA demonstrates the cancer pathway of C1; H GSEA shows the immune pathway of C2 subtype. I Unifactorial demonstration of clinicopathologic factors and ER subtypes. J Multifactorial demonstration of clinicopathologic factors and ER subtypes

Functional enrichment analysis and immunoinfiltration landscape of two subtypes

Our analysis also indicated that C2 is more responsive to immunotherapy due to its high expression in immune-related pathways (Supplementary Fig. 2A). We collected TILs scores from different samples and found that TILs scores were higher in the C2 group (Supplementary Fig. 2C). To explore the roles of ER subtypes in the LUAD tumor microenvironment (TME), we initially utilized the ESTIMATE algorithm to assess the overall immune infiltration between the two subgroups, encompassing stromal, estimate, and immunity scores. The C2 subgroup demonstrated elevated scores in all these categories (Supplementary Fig. 2D). Subsequently, leveraging the CIBERSORT algorithm, we identified a higher presence of plasma cells and M0-type macrophages in the C1 subgroup of lung adenocarcinoma patients, whereas the C2 subgroup exhibited increased infiltration of M2-type macrophages (Supplementary Fig. 2F). We then evaluated the relationship between the two subgroups and 28 immune cells using the ssGSEA method. This analysis revealed significantly higher infiltration levels of eosinophils, immature dendritic cells, natural killer T cells, and natural killer cells in the C2 subgroup. In contrast, activated CD4 T cells were notably overexpressed in the C1 subgroup. Additionally, we observed higher levels of MDSC, activated CD8 T cells, and effector memory CD8 T cells in C2 compared to C1 (Supplementary Fig. 2G). Employing the ssGSEA algorithm, we obtained scores for immune cells and associated pathways, indicating heightened activity in immune pathways within the C2 subgroup (Supplementary Fig. 2B). Finally, comparing the immune escape scores of the two subgroups using the TIDE algorithm, we found a higher TIDE score in C1, suggesting a greater likelihood of immune escape in this subgroup and potentially better immunotherapy efficacy for the C2 subgroup (Supplementary Fig. 2E).

In summary, the LUAD patients were classified into two distinct subtypes based on survival-related efferocytosis genes. The C2 subtype showed higher immune and stromal scores, robust immune-related pathway activation, and elevated antitumor activity, while C1 presented with poorer prognosis, oncogenic activation, and a heightened immune escape profile.

Construction and validation of ER-related prognostic features

Next, based on the 33 prognostically relevant ERGs, we used LASSO and multivariate Cox regression to construct ER-related prognostic features. We first performed a tenfold cross-validated LASSO regression analysis on these 33 genes and screened out 7 genes (ENTPD1, HAVCR1, IL1A, MERTK, P2RX1, PANX1, SPHK1) for further analysis. And we found that these 7 genes have good diagnostic efficacy in the early stage of lung adenocarcinoma [37,38,39,40,41,42]. The efferocytosis-associated risk score (ERGRS) for each LUAD patient was based on the following equation calculated as ERGRS = − 0.403*ENTPD1 + 0.248*HAVCR1 + 0.143*IL1A + − 0.221*MERTK + − 0.346*P2RX1 + 0.305*PANX1 + 0.331*SPHK1 (Fig. 4A, B and Supplementary Fig. 3A). We assigned LUAD patients to the high-risk or low-risk group based on the median risk score. Kaplan–Meier analysis showed that patients in the high-risk group had worse OS (p < 0.05); Fig. 4C, E and Supplementary Fig. 3B), and the analysis of subjects' work characteristics (ROC) curves showed that the area under the curve (AUC) of the ERGRS in the TCGA training set at 1, 3, and 5 years respectively reached 0.71, 0.69, and 0.70, indicating good predictive performance, 0.68,0.76,0.74 for validation cohort GSE31210; and 0.59, 0.68, and 0.67 for validation cohort GSE50081 (Fig. 4D, F and Supplementary Fig. 3C). The clinical risk factor heatmap demonstrated the differences in expression of these seven genes in the high and low risk groups (Fig. 4J, K Supplementary Fig. 3D). In addition, the survival curves also showed that the low-risk group had a favorable prognostic performance (Supplementary Fig. 3E–M). We validated prognostic models in GSE30219, GSE37745, and GSE68465 cohorts, finding high-risk patients had worse OS (p < 0.05). AUC scores ranged from 0.53 to 0.68 across cohorts. ERGRS scores varied by tumor origin, with head and neck, and renal papillae tumors showing higher scores, and prostate, adrenocortical, and thymoma tumors showing lower scores (Supplementary Fig. 4A–D). Malignant tumor infiltration and metastasis are key indicators of cancer severity, with EMT being a critical mechanism. Angiogenesis is also vital for tumor growth and spread [36]. We explored the link between ERGRS and these malignancy markers by analyzing z-scores for angiogenesis, EMT, and cell cycle in relation to ERGRS (Supplementary Fig. 5A). We discovered significant correlations between ERGRS and angiogenesis (R = 0.4, p < 0.0001), EMT (R = 0.46, p < 0.0001), and cell cycle (R = − 0.027, p < 0.05) in the TCGA pan-cancer cohort and across most tumor types (Supplementary Fig. 5C, D). High ERGRS tumors tend to have more active angiogenesis and aggressive cells. We also calculated the average ERGRS score per tumor, showing that ERGRS reflects a tumor's malignant potential (Supplementary Fig. 5B).

Fig. 4

Construction and validation of ER-related prognostic features. A Ten-fold cross-validation of parameter selection adjusted by LASSO regression. B Screening of coefficients under LASSO analysis. A vertical line is plotted at the value selected by tenfold cross-validation of overall survival. C, E KM survival curves for high and low risk groups (C) TCGA-LUAD; (E) GSE31210; D, F Time-dependent ROC curve analysis (D) TCGA-LUAD; (F) GSE31210; J, K TCGA- LUAD (J) cohort, and (K) distribution of risk scores and patient survival between low and high risk groups in GSE31210. G Column plots combining age, grade, sex, N-stage, total stage, and risk score; H Constructed 1-, 3-, and 5-year survival column plots. calibration curves. I DCA decision curve analysis. L Venn diagram showing efferocytosis gene characteristics obtained by two machine learning algorithms. M HR and 95% CI of members of the efferocytosis gene signature calculated by multivariate Cox regression analysis. N Neural network architecture of the attention mechanism, the binary classification algorithm

Clinical characteristics of prognostic models and construction of diagnostic models

To make the ERGRS more suitable for clinical applications, we constructed a nomogram based on the ERGRS and clinical features (Fig. 4G). Decision curve analysis (DCA) showed that the nomogram had better clinical benefits than other clinical features (Fig. 4I), and the calibration curve showed that the nomogram predictions were in good agreement with the actual observations (Fig. 4H). To demonstrate the prognostic efficacy of ERGRS, we compared ERGRS with other clinical factors, which showed strong prognostic efficacy of ERGRS. To further illustrate the prognostic efficacy of ERGRS, we compared the prognostic efficacy of ERGRS with the existing LUAD model, and we integrated the results of previous studies that used different biologically significant features, such as immune gene signatures [43], B-cell marker genes [44], collagen signatures [45], immune activation-associated gene signature [46], NK cell signature [47], basement membrane-related genes [48], metabolism-immunity-related genes [49] and mitotic spindle-related features [50] etc. Notably, ERGRS exhibited better C index performance than almost all models in the TCGA-LUAD, GSE31210 and GSE50081 datasets (Supplementary Fig. 6A–C).

To reveal the efferocytosis features in lung adenocarcinoma (LUAD) prognosis and immunotherapy, we strategically employed RSF and LASSO regression algorithms to pinpoint key efferocytosis genes associated with patient outcomes. The intersection of genes selected by these algorithms led us to identify IL1A, HAVCR1, and SPHK1 as central to efferocytosis and significantly correlated with LUAD prognosis. Integrating these genes into our predictive model, we observed a robust correlation with LUAD outcomes, underscoring the model's clinical applicability (Fig. 4L). Furthermore, Fig. 4M illustrates the multifactorial Cox analysis of these genes, reinforcing their prognostic significance.

In an effort to substantiate the predictive power of these efferocytosis gene features, we developed an attention-based DNN classification model using the TCGA-LUAD cohort as our training set (Fig. 4N). With the three efferocytosis genes as inputs and ER risk as the classification label, our model achieved an accuracy of 0.78, a sensitivity of 0.8750, a specificity of 0.88, and an AUC of 0.845 (Supplementary Fig. 8A–C). These results not only demonstrate the genes' potential as biomarkers for ER risk prediction but also highlight their role in the immunotherapeutic response of LUAD patients, aligning with our study's goal of uncovering the intricate links between efferocytosis and LUAD prognosis and treatment efficacy.

ERGRS has excellent predictive power for immunotherapy

In our research, we found that low ERGRS group had elevated immunity and stroma scores and higher ESTIMATE scores compared to the high ERGRS group, showed lower tumor purity (Supplementary Fig. 7A). Immune cell infiltration analysis revealed increased CD8 + T cells, B cells, and DCs in the low-ERGRS group (Supplementary Fig. 7B). Antigen presentation genes and immune checkpoints like IDO2, BTLA, CTLA4, and HAVCR2 were differentially expressed, hinting at the low ERGRS group's potential in immunotherapy (Supplementary Fig. 7C, D). Additionally, the high ERGRS group had higher TMB and TNB, suggesting greater immunogenicity (Supplementary Fig. 7E, F). Therefore, we believe that ERGRS may affect the immunotherapy of lung adenocarcinoma. To comprehensively assess the role of ERGRS in LUAD immunotherapy, we performed a systematic analysis. First, we analyzed a cohort of uroepithelial cancers treated with anti-PD-L1 therapy (IMvigor210), and the low-risk group had a significant survival advantage compared with the high-risk group (Fig. 5A). Meanwhile, patients with low risk scores were more sensitive to immunotherapy (Fig. 5B, C). Next, in the GSE78220 cohort, the low-risk score also had a strong ability to predict prognosis and immunotherapy benefit (Fig. 5D). Immune checkpoint inhibitors (ICBs) are revolutionary immunotherapeutic agents that block inhibitory signals of T-cell activation, enabling tumor-reactive T cells to generate effective anti-tumor responses [51]. To further explore the role of risk scores in immunotherapy, we explored the association between ERGRS and ICB-related positive signaling pathways. The results showed that ERGRS was positively correlated with cell cycle, DNA replication, base excision repair, and viral oncogenic effects (Fig. 5E). The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was used to assess patient response to immunotherapy and showed better responsiveness in the low ERGRS group (P [Fisher exact test] = 3.50e−07; Fig. 5F). Submap analysis of another group of melanoma patients receiving immunotherapy (GSE91061) also showed a better response to PD-1 treatment in the low ERGRS (Bonferroni-corrected p < 0.01, Fig. 5H). Finally, the results showed that low ERGRS tended to represent better immunotherapy outcomes (GSE91061, p = 0.014; Fig. 5G, I). High ERGRS patients had worse outcomes and were enriched in cancer pathways like E2F targets and EMT. Since they didn't respond well to immunotherapy, we looked for drugs in CTRP and PRISM databases. Cisplatin showed promise, especially for patients with low CXCL10, suggesting it could benefit chemotherapy. We also identified potential treatments from studies and databases, including BRD-K92856060, monensin and hexaminolevulinate (Supplementary Fig. 8D–I).

Fig. 5

ERGRS has excellent predictive power for immunotherapy. A Survival curves for the high and low ERGRS groups in the IMvigor210 cohort. B Box plots depicting the difference in risk scores between CR/PR patients and SD/PD patients in the IMvigor210 cohort. C Proportion of CR/PR or SD/PD patients receiving immunotherapy in the high and low risk groups of the IMvigor210 cohort. D Proportion of patients with CR/PR or SD/PD who received immunotherapy in the high and low risk groups of the GSE78220 cohort. E Correlation of ERGRS with ICB response characteristics and each step of the tumor immune cycle. F The TIDE algorithm predicts the response to immunotherapy between the high ERGRS and low ERGRS groups. G Proportion of R or NR patients receiving immunotherapy in the high and low risk groups of the GSE91061 cohort. H Submap algorithm predicting response to immunotherapy between the high and low ERGRS groups. I Distribution of ERGRS in different immunotherapy response groups in GSE91061

HAVCR1 promotes lung adenocarcinoma progression

We analyzed the difference in the expression of HAVCR1 between tumor tissues and normal tissues in TCGA pan-cancer. The results showed that HAVCR1 was highly expressed in most tumors (Fig. 6A). In the above analysis, we explored the expression of HAVCR1 in pan-cancer and has not been fully investigated in LUAD. Therefore, we performed a series of cellular experiments. The knockdown effect of HAVCR1 is shown in Fig. 6B, with interference 1 exhibiting higher knockdown efficiency. CCK8 assay showed that inhibition of HAVCR1 may significantly inhibit the proliferative capacity of LUAD cells, with the inhibitory effect of Interference 1 being even more pronounced (Fig. 6C). Immunohistochemical staining of HAVCR1 was performed on 12 lung adenocarcinoma tissues and surrounding lung tissues from the First Affiliated Hospital of Bengbu Medical University. The results showed that all lung adenocarcinoma tissues were ( +) and surrounding lung tissues were (−) (Fig. 6D). And it was also demonstrated by cell scratch assay that knockdown of HAVCR1 significantly inhibited the migratory ability of LUAD cells (Fig. 6E). The results of the Transwell assay showed that knockdown of HAVCR1 significantly inhibited the migration and invasion ability of LUAD cells (Fig. 6F). All these results demonstrated that HAVCR1 was a prognostic marker for LUAD. The same trend was verified by colony formation assay (Fig. 6G).

Fig. 6

HAVCR1 promotes lung adenocarcinoma progression. A Box line plot demonstrating the expression level of HAVCR1 in pan-cancer. B RT-qPCR screening for suitable interferences. C CCK8 assay experiment. D Immunohistochemical images of tumor tissue and normal tissue; E Cell scratch experiment. F Transwell assay validation of the interference. G Clonal spot experiment. (*p < 0.05. **p < 0.01. ***p < 0.001. ****p < 0.0001. ns, p > 0.05)

Discussion

In recent years, the 5-year survival rate of lung adenocarcinoma has increased from 5 to 15% due to the use of novel chemotherapeutic agents, targeted drugs and immunosuppressants, but it remains the leading cause of life-threatening health problems in cancer. Immune checkpoint inhibitors (ICIs) offer significant advantages in terms of efficacy and safety, bringing new hope for NSCLC treatment [52]. Given these circumstances, generating prognostic models based on key biomarkers that can predict patient survival, immune microenvironmental status, gene mutations, and potential drug sensitivities has great potential for improving the accuracy of LUAD therapy. For example, Xie.et.al studied G protein-coupled receptor-related models [53] and Lian.et.al identified T cell exhaustion related models [54]. Each has provided new insights into the treatment of LUAD patients. However, the identification of biomarkers associated with the role of efferocytosis in LUAD is very scarce.

In this study, we performed a comprehensive bioinformatics analysis of efferocytosis in LUAD using data from TCGA-LUAD and 5 GEO validation cohorts. A profile of efferocytosis-related genes and TME features in LUAD was revealed. Molecular typing using 33 efferocytosis genes showed that C1 had a poorer prognosis and a higher infiltration of plasma cells and M0-type macrophages. The C1 tumor proliferative capacity and malignant phenotype was demonstrated by GSVA results. In contrast, C2 had a better prognosis and was highly correlated with immunity. Next, in we screened ER-related prognostic features and constructed a stable and reliable ERGRS-related prognostic signature including seven ER-related prognostic signatures (PANX1,P2RX1, HAVCR1, SPHK1, IL1A, MERTK and ENTPD1). We compared the ERGRS model to existing LUAD prognostic models, finding it has a high AUC, indicating strong predictive accuracy, and performed well across various datasets. Its generalizability is supported by consistent performance in external cohorts. ERGRS offers unique advantages by predicting survival and providing immune landscape insights, crucial for personalized treatment and immunotherapy. It also stratifies patients into high- and low-risk groups with precision, surpassing traditional models. Our ERGRS was validated by a pan-cancer cohort to have malignant tumor proliferation characteristics. ERGRS not only has strong predictive ability in prognosis, but also effectively predicts immunotherapy response in LUAD patients. The tumor immune microenvironment consists of immune cells, stromal cells, and a series of cytokines that provide the environment for tumor survival [55]. Components of this tumor microenvironment can modulate the immune phenotype of the tumor, thus affecting the survival prognosis of patients [56]. In this study, the immune cell infiltration in the tumor microenvironment of patients with high and low ERGRS was thoroughly analyzed. Using multiple immune infiltration algorithms, a significantly higher abundance of macrophage M0 and M1 infiltration was found in the high-risk population compared to the low-risk population.

In the ERGRS-based risk grouping, we found that the low ERGRS group predominantly showed better OS, high immune infiltration, and lower TNB, TMB. The low ERGRS group had higher immune checkpoint expression, suggesting that patients with low ERGRS may have a more favorable response to ICI therapy. Two widely used prediction tools, the TIDE and Submap algorithms, also demonstrated a better response to immunotherapy in the low ERGRS group, which is consistent with our analysis and suggests that ERGRS may be useful for early identification of immunotherapy-sensitive populations.

To address the current poor response of the high ERGRS group to immunotherapy, we systematically screened them for potential therapeutic agents using a comprehensive screening framework that has been proven to be effective in previous studies [57, 58]. In this study, one CTRP-derived drug (BRD-K92856060) and two PRISM-derived drugs (monensin and hexaminolevulinate) were screened as possible candidates for the high-ERGRS group.

We screened ER genes by RSF and LASSO and constructed a diagnostic model by deep learning, which showed that IL1A, HAVCR1 and SPHK1 were the most central genes [59]. Interestingly, previous studies have shown that these genes play important roles in regulating the progression and prognosis of many cancers. IL1A has been reported to be associated with inflammatory diseases, cancer and immunotherapy resistance [60, 61]. SPHK1 promotes the development of non-small cell lung cancer by activating the STAT3 and EMT pathways [62, 63]. HAVCR1 has been less studied in lung adenocarcinoma, which prompted us to explore the function of HAVCR1 as a potential prognostic biomarker for lung adenocarcinoma. In follow-up experiments, we found that HAVCR1 promoted the proliferation, migration and invasive ability of lung adenocarcinoma cells. Finally, by analyzing the Pan-TCGA dataset, we found that HAVCR1 was significantly up-regulated in multiple cancers and was significantly associated with poor prognosis.

In conclusion, this study introduces a novel efferocytosis-associated risk profile that accurately predicts clinical prognosis and treatment response in patients with lung adenocarcinoma (LUAD). Efferocytosis has been poorly studied in LUAD and provides new insights into the critical role of efferocytosis in cancer progression. However, limitations of our work remain. The current results highlight the need for prospective clinical trials to further validate the clinical applicability of our ERGRS. The mechanisms underlying the poor prognosis of patients at high ERGRS risk that we obtained should be further explored in experiments. In addition, further in vivo experiments are needed to explore the functional role of HAVCR1 in lung adenocarcinoma, which could help to provide stronger clues to guide clinical application.

Integration of ERGRS into clinical practice

As we look to the future, the integration of the ERGRS into existing clinical workflows and decision-making processes is of paramount importance. We propose the following strategies to ensure that ERGRS can be effectively utilized in a clinical setting:

Integration with clinical decision support systems (CDSS): We propose to integrate the ERGRS algorithm into existing CDSS. This will allow clinicians to receive real-time, personalized medical recommendations based on the ERGRS score as part of their diagnostic and treatment decision-making process.

Workflow optimization: We understand that integrating ERGRS into the clinical environment may require adjustments to current workflows. This includes embedding an ERGRS scoring module within the Electronic Health Record (EHR) system to naturally incorporate ERGRS scoring into the daily workflow of physicians.

Quality improvement program: We view the implementation of ERGRS as an ongoing quality improvement process. We will regularly evaluate the effectiveness of ERGRS in clinical decision-making and make adjustments based on feedback.

We believe that through these measures, ERGRS can be effectively integrated into clinical workflows, thereby improving the accuracy of treatment and survival rates for patients with lung adenocarcinoma (LUAD).