Skip to main content
Download PDF

Interaction, immune infiltration characteristics and prognostic modeling of efferocytosis-related subtypes in glioblastoma

BMC Medical Genomics volume 16, Article number: 248 (2023) Cite this article

Abstract

Background

Efferocytosis is a biological process in which phagocytes remove apoptotic cells and vesicles from tissues. This process is initiated by the release of inflammatory mediators from apoptotic cells and plays a crucial role in resolving inflammation. The signals associated with efferocytosis have been found to regulate the inflammatory response and the tumor microenvironment (TME), which promotes the immune escape of tumor cells. However, the role of efferocytosis in glioblastoma multiforme (GBM) is not well understood and requires further investigation.

Methods

In this study, we conducted a comprehensive analysis of 22 efferocytosis-related genes (ERGs) by searching for studies related to efferocytosis. Using bulk RNA-Seq and single-cell sequencing data, we analyzed the expression and mutational characteristics of these ERGs. By using an unsupervised clustering algorithm, we obtained ERG clusters from 549 GBM patients and evaluated the immune infiltration characteristics of each cluster. We then identified differential genes (DEGs) in the two ERG clusters and classified GBM patients into different gene clusters using univariate cox analysis and unsupervised clustering algorithms. Finally, we utilized the Boruta algorithm to screen for prognostic genes and reduce dimensionality, and the PCA algorithm was applied to create a novel efferocytosis-related scoring system.

Results

Differential expression of ERGs in glioma cell lines and normal cells was analyzed by rt-PCR. Cell function experiments, on the other hand, validated TIMD4 as a tumor risk factor in GBM. We found that different ERG clusters and gene clusters have distinct prognostic and immune infiltration profiles. The ERG signature we developed provides insight into the tumor microenvironment of GBM. Patients with lower ERG scores have a better survival rate and a higher likelihood of benefiting from immunotherapy.

Conclusions

Our novel efferocytosis-related signature has the potential to be used in clinical practice for risk stratification of GBM patients and for selecting individuals who are likely to respond to immunotherapy. This can help clinicians design appropriate targeted therapies before initiating clinical treatment.

Peer Review reports

Introduction

Glioblastoma (GBM) is the most common primary brain malignancy, accounting for approximately 12%-15% of all brain tumors [1, 2]. Despite significant advancements in chemotherapy, radiation therapy, and surgical treatment, GBM patients’ 5-year survival rate remains less than 5% [3]. Although epidemiological studies suggest that ionizing radiation increases glioblastoma incidence [4], most GBM patients have no clear pathogenetic cause. The widespread heterogeneity within and between individuals is the root cause of GBM treatment failure, making it one of the most aggressive and treatment-resistant malignancies [5]. Therefore, discovering new biomarkers and establishing effective molecular staging systems to select appropriate treatments for GBM patients is crucial, as molecular alterations are increasingly important in glioma classification and grading [6, 7].

Phagocytes, such as macrophages and immature dendritic cells, play a vital role in efferocytosis, the process of recognizing and engulfing dying cells during apoptosis [8]. Unlike regular cytokinesis, efferocytosis preserves the membrane integrity of dead cells, preventing exposure to immunogenic substances and avoiding secondary cell damage caused by inflammatory responses [9, 10]. Although numerous genes that promote efferocytosis are involved in tumor development and metastasis and are frequently overexpressed in various cancers, including lung cancer, breast cancer, and leukemia [11, 12]. Current theories suggest that efferocytosis may contribute to tumor progression due to the unaccompanied release of inflammatory factors and production of killer effector T cells, resulting in a suppressive immune microenvironment and the immune escape of tumor cells [12, 13]. In contrast, uncleared apoptotic cells and secondary necrosis promote a proinflammatory environment and antitumor immunity [14]. The receptors for phagocytes, such as tumor-associated macrophages (TAM), have been extensively studied. For example, MerTK is involved in epidermal growth factor receptor (EGFR) inhibitor resistance in non-small-cell lung cancer [15], and blocking phagocytic receptors with the membrane-linked protein V can effectively slow tumor progression in prostate cancer [16]. Additionally, TIM-4 acts as a PS receptor on the surface of phagocytes and promotes angiogenesis in colorectal cancer by upregulating vascular endothelial growth factor (VEGF) [17].

The crucial role of efferocytosis in cancer development and progression is attributed to its effect on tumor cell growth, metastasis, EMT, and angiogenesis [18]. Although traditional oncology therapies such as chemotherapy and radiotherapy trigger apoptosis and efferocytosis, they also lead to tumor inflammation and limit antitumor immunity [19]. Studies have indicated that solely blocking efferocytosis cannot completely inhibit the production of tumor immunosuppressive cells and mediators [20]. However, combined inhibition of tumor cell apoptosis and efferocytosis has been observed to effectively suppress metastatic recurrence of tumors. Notably, the immunosuppressive microenvironment in GBM patients forms a comprehensive and self-sufficient system [21]. Additionally, the efferocytosis process may function as an immune checkpoint similar to PD-1/PD-L1, which could be targeted for therapeutic interventions [22]. Therefore, developing a combination therapy that targets both conventional oncology therapy and efferocytosis presents a significant technical challenge.

In this study, we categorized 20 genes associated with efferocytosis into different clusters and assessed their impact on the prognosis of GBM patients. To achieve this, we used an unsupervised consensus clustering method and combined three GBM cohorts. The univariate Cox analysis and Boruta algorithm were used to identify the differentially expressed genes among the ERG clusters. Then, a scoring system was established based on the PCA algorithm. The primary objective of the study was to determine whether this novel efferocytosis characteristic could accurately predict the prognosis of GBM patients and assist medical professionals in identifying potentially responsive patients for the development of effective immunotherapies.

Materials and methods

Acquisition of raw data

We collected RNA-seq data and clinical information of GBM patients from two databases, UCSC Xena (https://xena.ucsc.edu/) and CGGA (http://www.cgga.org.cn/). In addition, we downloaded mutation data of TCGA-GBM patients and removed 8 duplicates, leaving us with a TCGA cohort comprising 161 GBM tissue samples and 5 normal samples. The CGGA-325 cohort contained 139 GBM samples, while the CGGA693 cohort had 249 GBM samples. Gene expression profiles were measured using the transcript per million estimation and log2-based transformation. After combining the mRNA expression data of GBM patients from these three cohorts, we used the “sva” package [23] to perform batch correction. We identified 103 efferocytosis-related genes from the genecards portal (http://www.genecards.org/) using the keyword “efferocytosis”. After comparing these genes on the Pubmed website (https://pubmed.ncbi.nlm.nih.gov/) for their research applications in efferocytosis, we selected the most plausible 22 genes. These 22 efferocytosis-related genes are listed in Supplementary Table 1, and we have provided a brief illustration of the efferocytosis process in Fig. 1.

Fig. 1
figure 1

Brief process of efferocytosis

Full size image

Consensus unsupervised clustering

To identify distinct ERG clusters based on the expression of 20 ERGs, we employed consensus unsupervised cluster analysis using the “ConsensusClusterPlus” software program [24]. We used the “Pam” algorithm with “Euclid” as the distance measure, and resampled items with a rate of 80% for 1000 replications to determine the optimal k value based on the proportion of ambiguous clustering (PAC). Next, we used the “limma” package to identify differentially expressed genes among the various ERG clusters, where expression levels were considered significant if |log2 FC| exceeded 1 and the adjusted P-value was below 0.05 [25, 26].

Functional enrichment analysis

To explore the biological functions of the differentially expressed genes in various ERG clusters, we conducted GO enrichment and KEGG pathway analyses. We used the clusterProfiler package and applied the BH method to adjust the P-value [27]. It is important to note that the GO enrichment functions include cellular components, biological processes, and molecular functions [28, 29].

Establishment of the efferocytosis-related signature

After conducting one univariate Cox analysis (with a significance level of P < 0.05), we proceeded to a second clustering analysis on differentially expressed genes (DEGs) identified in the various ERG clusters. Among the DEGs associated with the ERG gene cluster, we defined genes positively associated with the ERG gene cluster as gene signature A and genes negatively associated with the ERG gene cluster as gene signature B. We utilized the “Boruta” software package to further screen significant genes among the candidate genes, applying the Boruta algorithm with a maxRuns value of 500. By using the Boruta algorithm, we retained the genes identified as “confirmed”. Principal component analysis (PCA) was applied to reduce the dimensionality of the ERG gene cluster. Subsequently, an ERG score was assigned to each patient by computing the score for each GBM sample using the following formula: score = ∑PCA A—∑PCA B. To classify GBM patients into high-risk (HR) and low-risk (LR) groups, we used the “surv_cutpoint” function from the “survminer” package to determine the optimal cutoff values [30].

Prediction of immunotherapy

We aimed to evaluate the ability of the ERG signature to predict the response to immunotherapy in GBM patients, by utilizing the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm that integrates both tumor immune dysfunction and exclusion factors [31]. We obtained a subset of genes associated with cancer and immunity by using the website developed by Xu et al. [32] (http://biocc.hrbmu.edu.cn/tip/) and selecting genes that were positively associated with anti-PD-L1 drug response based on Mariathasan’s study features [33]. We then applied the GSVA method to calculate the enrichment scores for gene signatures associated with the cancer immune cycle, considering p-values less than 0.05 to indicate statistically significant differences between the two groups. To assess the relationship between risk scores and the two genetic features mentioned above, we used the R package “ggcor”.

Immune microenvironment-related analysis

We utilized various algorithms for immune cell infiltration analysis, including XCELL [34], TIMER [35], QUANTISEQ, MCPCOUNT, EPIC [36], CIBERSORT [37] and CIBERSORT- ABS. The results from these algorithms were compared and further analyzed using the “ComplexHeatmap” R package [38, 39]. We examined the correlation between immune cells and risk scores using Spearman correlation analysis. To differentiate GBM patients with low ERG scores from those with high scores based on their immune cell features, we employed the single-sample gene set enrichment analysis (ssGSEA) technique. We also estimated the immune and stromal scores of each glioma sample using the R program “ESTIMATE,” which provides an estimate of the quantities of immune and stromal components present in vivo [40]. To comprehensively investigate the tumor microenvironment’s heterogeneity in different datasets and cell types, we utilized the Tumor Immune Single-Cell Hub (TISCH) database. TISCH is an extensive single-cell RNA-seq database dedicated to the tumor microenvironment (http://tisch.comp-genomics.org).

Drug sensitivity

The semi-inhibitory concentration, also known as IC50, represents the drug concentration that corresponds to a 50% ratio of apoptotic cells to the total number of cells and is often used to assess a drug’s ability to induce apoptosis. A lower IC50 value indicates a higher induction of apoptosis, whereas a higher value suggests that the cells are more tolerant to the medication. To assess the efficacy of our ERG signature in targeted chemotherapy, we used the “pRRophetic” tool to compute the IC50 values for various chemotherapeutic agents typically used in GBM treatment [41]. We compared the IC50 values between the high and low-scoring groups and evaluated the patients’ sensitivity to each drug.

Transfection of cells and real-time PCR

U251MG, LN229, and SW1783 human glioma cells and human astrocytes (NHA) were cultured in Dulbecco’s Modified Eagle’s Medium (DMEM, Gibco, C11995500BT, Canada) supplemented with 10% fetal bovine serum (FBS, Gibco, 10091148, Canada) and 1 × penicillin/streptomycin (Gibco, 15,140–122, Canada). All cultures were maintained in a CO2 incubator (TFS3111, USA) at 37 °C with 5% CO2. TIMD4 gene knockdown was achieved using small interfering RNA (siRNA). The specific TIMD4 siRNA sequences can be found in Supplementary Table 2. In brief, cells were seeded at 50% confluency in 6-well plates and transfected with negative control (NC) and siBARD1 using Lipofectamine 3000 (Invitrogen, USA).

Total RNA was extracted from cell lines and tissues using TRIzol (Sigma-Aldrich, T9424, America) according to the manufacturer’s instructions. cDNA was synthesized using the PrimeScriptTM RT Reagent Kit (Takara, RR047, Japan). Real-time polymerase chain reaction (RT-PCR) was performed using SYBR Green Master Mix (Q111-02, Vazyme) to quantify mRNA expression levels normalized to GAPDH mRNA levels. The 2 − ΔΔCt method was used to calculate the expression levels. All primers were provided by Qingdao BioScience (Beijing, China), and the primer sequences can be found in Supplementary Table 2.

Cell counting Kit8 assay and transwell assay

First, cells (1000 cells per well) were seeded into a 96-well plate and incubated at 37 °C for 4 h with CCK-8 reagent (10 μL) (Dojindo, CK18, Japan). The absorbance was measured at a wavelength of 450 nm using an ELx800 plate reader (Thermo, Multiskan Spectrum, USA) to count the cells. Cell growth was represented as fold change from day 0 to day 4 and presented in a graph.

Cell invasion and migration studies were performed using a transwell assay. The upper chambers of a 24-well plate were filled with treated SW1783 cells (2 × 10^5 cells) and incubated for 48 h. To evaluate the invasive and migratory abilities of the cells, the top surface of the plate was pre-coated with a matrix gel solution (BD Biosciences, USA) or left uncoated. The remaining cells at the bottom layer were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet (Solarbio, China) after removing the surface cells.

Statistical analysis

All analyses were performed using R version 4.1.1, 64-bit, and its support package. To assess prognostic value and compare patient survival in various subgroups within each data set, Kaplan–Meier survival analysis and log-rank tests were performed. For significance tests comparing the various groups, Kruskal–Wallis and Wilcoxon’s tests were performed. By displaying univariate and multivariate forest plots, we examined if this created ERG signature is an independent predictive factor in comparison to other clinical features. The “stats” package and “prcomp” function were used to perform principal component analysis. Spearman correlation analysis was used to investigate the correlation coefficients. In all statistical investigations, P < 0.05 was considered statistically significant.

Results

Genetic and transcriptional alterations of ERGs in GBM

In this study, we collected a total of 22 ERGs, 21 of which were identified in the TCGA and CGGA cohorts. Initially, we compared the expression levels of ERGs between GBM and normal tissues using TCGA, CGGA, and GTEx expression profiles. Our analysis showed that all ERGs, except RAB17, had different expression levels between tumor and normal tissues (Fig. 2A). We also examined the somatic mutation frequencies of the 20 differentially expressed ERGs and found that they had low mutation frequencies, with only 31 out of 390 GBM samples (7.95%) having mutations in these genes (Fig. 2B). In addition, somatic copy number variation analysis showed a general decrease in copy number variation (CNV) for genes such as FPR2, TYRO3, IGF2R and AXL, while GAS6 showed a gain in CNV (Fig. 2C). The chromosomal localization of these ERGs is shown in Fig. 2D. We constructed an efferocytosis-related network to demonstrate the comprehensive landscape of ERG interactions, regulator connections, and their prognostic value in patients with GBM (Fig. 2E). Additionally, we compared the expression levels of TIMD4 in HA cells and three GBM cell lines by cell line experiments and found that BARD1 was significantly highly expressed in tumor cells, especially SW1783 cells (Fig. 2F). We then examined the expression level of TIMD4 5 days after transfection by qRT-PCR to test the effectiveness of siRNA knockdown of TIMD4 in SW1783 cell lines (Fig. 2G). Subsequent CCK-8 cell assays showed that knockdown of BARD1 significantly reduced the proliferative capacity of the SW1783 cell line (Fig. 2H). In addition, GBM cells transfected with si-TIMD4 exhibited weaker migratory invasion ability in transwell assays (Fig. 2I, J). Thus, TIMD4 is a pro-carcinogenic factor in GBM. These findings suggest that there are significant differences in the genetic profiles and expression levels of ERGs between GBM and control samples, and that ERGs may play a crucial role in the development of GBM.

Fig. 2
figure 2

Genetic and transcriptional alterations of ERGs in GBM. A Distribution of expression of ERGs between normal and GBM. B Mutation frequency of 20 ERGs in the TCGA cohort of 390 GBM patients. C Copy number variation (CNV) of 20 ERGs in TCGA-GBM. D Localization of 20ERGs in chromosomal regions. E Network plot showing the correlation between the 20 ERGs. Red connecting lines indicate positive correlations, while blue indicates negative correlations. F TIMD4 was highly expressed in GBM cell lines compared to normal human astrocyte NHA cell lines. G RT-qPCR was performed to detect the relative expression of TIMD4 in GBM cells transfected with si-RNAs or negative control (NC). H CCK8 assay showed that SW1783 cells with reduced TIMD4 expression had significantly reduced proliferative capacity compared to the NC group. I, J Transwell assay showed that down-regulation of TIMD4 expression inhibited the migration and invasion ability of SW1783 cells. All data are expressed as mean ± SD of three independent experiments. * p < 0.05, ** p < 0.01, *** p < 0.001

Full size image

Validation of single-cell sequencing data

To investigate the expression of 20 ERGs within the tumor microenvironment (TME), we analyzed the GBM single-cell dataset GSE141982 obtained from the TISCH database. Among the 20 ERGs, MPO expression was not detected. The GSE141982 dataset consisted of 16 cell clusters and 4 major cell types, which were distributed and counted as shown in Fig. 3A and B. Our analysis revealed that IGF2R, NCF1, and FPR2 were predominantly expressed in CD8T cells, with lower expression levels observed in malignant cells. Conversely, FN1 and GAS6 were primarily expressed in endothelial cells. Notably, our findings showed that almost all ERGs were associated with immune cell infiltration, indicating that efferocytosis plays a crucial role in the GBM immune microenvironment (Fig. 3C, D).

Fig. 3
figure 3

20 ERGs in single-cell RNA sequencing. A, B Annotation of all cell types in GSE141982 and the percentage of each cell type. C, D The expression of 20 ERGs in each cell type

Full size image

Identification of ERG clusters

We utilized the TCGA-GBM, CGG325, and CGGA-693 cohorts and performed PCA analysis to demonstrate a significant reduction in the corrected batch effect (Fig. 4A, B). The optimal number of clusters was determined to be k = 2, based on our findings, and we divided the 549 GBM patients into two clusters (Fig. 4C, Supplementary Figure S1, and Supplementary Table 3). Patients in cluster A had a worse prognosis than those in cluster B, as evidenced by Kaplan–Meier survival analysis (Fig. 4D, P < 0.001). Significant differences in transcriptional profiles between the two clusters were also observed using PCA analysis (Fig. 4E). We conducted the “ssGSEA” algorithm to explore the tumor microenvironment in both clusters, and the results indicated that cluster A had a higher abundance of immune cells than cluster B, except for CD56dim NK cells and type 2 helper T cells (Fig. 4F).

Fig. 4
figure 4

Identification of two ERG clusters. A Principal component analysis of common gene profiles before the combination of the TCGA-GBM, CGGA-325, and CGGA-693 cohorts. B Principal component analysis of common gene profiles after the combination of TCGA-GBM, CGGA-325, and CGGA-693 cohorts. C Heat map of the consensus matrix defining two clusters (k = 2) and their associated regions. D Kaplan–Meier survival analysis of OS in 2 ERG clusters. E PCA analysis shows significant differences in the transcriptome between the two clusters. F The abundance of tumor-infiltrating immune cells between two ERG clusters was calculated by ssGSEA. ns no significance, *** p < 0.001

Full size image

We then examined ERG expression and clinicopathological features in both clusters and found significant differences, with cluster A exhibiting significantly higher ERG expression (Fig. 5A). GSVA enrichment analysis revealed that immune activation-related pathways were significantly enriched in cluster A, including the JAK/STAT signaling pathway, programmed cell death, cytokine receptor interactions, chemokine signaling pathway activation, NOD-like, and Toll-like receptor signaling pathway (Fig. 5B). Differential expression analysis between the two clusters revealed 641 DEGs (Supplementary Table 4), which were enriched in functions related to efferocytosis and immune-related pathways such as neutrophil-mediated immunity, immune receptor activity, cytokine receptor activity, and leukocyte activation (Fig. 5C). KEEG enrichment analysis showed that these DEGs were associated with the progression of certain autoimmune diseases (Fig. 5D).

Fig. 5
figure 5

Clinicopathological and biological characteristics of two ERG clusters. A Differences in clinicopathological features and expression levels of ERG between the two different subtypes. B Differential enrichment of the KEGG pathway between the two ERG clusters based on GSVA analysis. C, D GO and KEGG enrichment analyses of DEGs among two ERG clusters

Full size image

Identification of ERG gene clusters

To investigate the association between gene expression and GBM prognosis, we conducted univariate Cox regression analysis on the 641 DEGs and identified 296 genes that were prognostic for GBM (Supplementary Table 5). We then performed clustering analysis on these 296 DEGs, and found that the optimal number of clusters was 2 based on the slope of the cumulative distribution function curve (Fig. 6A, Supplementary Figure S2, and Supplementary Table 6). The heat map in Fig. 6B shows the expression of the 296 DEGs in the two gene clusters, as well as the clinicopathological characteristics of each sample. Remarkably, all 20 ERGs were more highly expressed in gene cluster A (Fig. 6C). Furthermore, Kaplan–Meier survival analysis revealed that patients in gene cluster A had a worse prognosis than those in gene cluster B (Fig. 6D, p < 0.001). To explore the tumor microenvironment of both gene clusters, we used the “ssGSEA” algorithm. Consistent with previous ERG clusters, differential enrichment of immune cells showed a higher abundance of immune cells in cluster A than in cluster B, except for CD56dim NK cells and type 2 helper T cells (Fig. 6E).

Fig. 6
figure 6

Identification and characterization of two ERG gene clusters. A Heat map of the consensus matrix defining the two gene clusters (k = 2) and their associated regions. B Relationship between clinicopathological features and the two gene clusters. C Expression differences of 20 ERGs between two ERG gene clusters. D Kaplan–Meier survival analysis of OS between two ERG gene clusters. E The abundance of tumor-infiltrating immune cells between two ERG gene clusters was calculated by ssGSEA. ns no significance,* p < 0.05, *** p < 0.001

Full size image

Development of an efferocytosis-related scoring system

We identified 296 differentially expressed genes (DEGs) and categorized 54 genes that were positively correlated with the ERG gene cluster signature as gene signature A, while 242 genes that were negatively correlated were assigned to gene signature B (Figure S6). To identify candidate genes from these important genes, we employed the “Boruta” package, resulting in 53 genes for signature A and 202 genes for signature B (Supplementary Table 7). We then calculated ERG scores for 549 GBM patients using the PCA formula described earlier. After determining the best cut-off value for the PCA score, we divided the TCGA and CGGA cohorts into high-risk (HR) and low-risk (LR) groups. Significant differences in overall survival (OS) were observed between the two groups in the whole cohort and in three separate cohorts, with patients in the LR group often having a better prognosis (Fig. 7A, D, E, F). The risk curves and survival status plots emphasized the strong discriminatory power of this ERG signature (Fig. 7B, C). The Sankey plots showed that the majority of ERG gene cluster A with poorer prognosis belonged to the HR group, indicating poorer survival outcomes (Fig. 7G). Both ERG cluster A (Fig. 7H) and ERG gene cluster A (Fig. 7I) with poorer prognosis had higher ERG scores, further highlighting the prognostic value of this signature.

Fig. 7
figure 7

Construction and validation of the ERG scoring system. A Kaplan–Meier survival analysis of OS between high- and low-achieving subgroups in the entire cohort. B Distribution of ERG scores in the whole cohort. C Relationship between ERG characteristics and survival status in the whole cohort. D Kaplan–Meier survival analysis of OS between high and low-scoring subgroups in the TCGA-GBM cohort. E Kaplan–Meier survival analysis of OS between high- and low-achieving subgroups in the CGGA-325 cohort. F Kaplan–Meier survival analysis of OS between high-scoring and low-scoring subgroups in the CGGA-693 cohort. G Sankey diagram demonstrating the relationship between patient survival status and ERG group. H The difference in scores between the two ERG clusters. I The difference in scores of two ERG gene clusters

Full size image

Validation of an efferocytosis-related scoring system

After conducting univariate Cox regression analysis, we observed that age, IDH mutation status, and ERG score were significantly associated with OS in all datasets of GBM patients (Fig. 8A). Subsequently, we performed multivariate Cox regression analysis and identified that age and ERG score were independent prognostic indicators for GBM patients (Fig. 8B). Additionally, using the chi-square test, we found that our ERG grouping was associated with the gender and IDH mutation status of the patients (Fig. 8C), and that lower ERG scores were linked to IDH mutations (Fig. 8D, E). These findings suggest that our developed ERG signature is a reliable predictor of OS for GBM patients, irrespective of other clinical characteristics.

Fig. 8
figure 8

Prognostic value of ERG scores in patients with GBM. A Univariate and B multivariate COX analysis to assess ERG signature and clinical characteristics (including age, gender, and IDH mutation status). C Histogram of clinical characteristics associated with ERG scores. D, E Correlation of ERG scores with IDH mutation status. ***p < 0.001

Full size image

The correlation with tumor microenvironment

In the tumor microenvironment (TME), immune cell infiltration plays a crucial role in immune response. Initially, we conducted Spearman correlation analysis in the TCGA cohort to examine the relationship between ERG scores and immune cell abundance in the GBM TME using different algorithms. Figure 9A illustrates the immune cell infiltration landscape in different risk groups, where we found that most immune cells were positively correlated with ERG scores. Interestingly, macrophage abundance linked to efferocytosis showed a positive correlation with ERG scores across all algorithms (Fig. 9B). To further investigate the association between ERG grouping and immune cells and functions, we quantified the ssGSEA enrichment scores for different immune cell subpopulations, related functions, or pathways. The results revealed that the high-scoring subgroups had more infiltration of B cells, CD8 + T cells, dendritic cells (DCs), immature DCs (IDCs), macrophages, neutrophils, plasmacytoid DCs (pDCs), helper T cells, type 1 T helper cells (Th1), type 2 T helper cells (Th2), tumor-infiltrating lymphocytes (TILs), and regulatory T cells (Tregs) (Fig. 9D). Similarly, all immune pathways, including APC_co_inhibition, Check-point, HLA, Inflammation-promoting, and T_cell_co-inhibition, were higher in the high-scoring subgroup (Fig. 9D). Furthermore, stromal scores, immune scores, and estimation scores were higher in the high-risk (HR) group (Fig. 9F). We also observed that some common immune checkpoints (ICs) were more highly expressed in the HR group, and we present the TME scores, immune checkpoint expression, and immune cell infiltration landscape between the two groups in the form of a heat map (Fig. 9C). Our findings in the CGGA-325 and CGGA-693 cohorts were similar to those discussed above (Supplementary Figure S3). Additionally, all 20 ERGs were more highly expressed in the HR group, suggesting that the higher overall immune level and immunogenicity of TME in the HR group were likely triggered by the effects of efferocytosis (Fig. 9E).

Fig. 9
figure 9

Analysis of the immune microenvironment in different risk groups. A Differences in immune infiltration status between different risk groups were evaluated by seven algorithms. B Bubble plot of the correlation between ERG score and immune cell abundance. C Heatmap showing differences in TME score, immune checkpoint expression, and immune cell infiltration calculated by ssGSEA among different risk subgroups. D Differences in ssGSEA scores of immune cells and immune function in the two score subgroups. E Differences in the expression of immune checkpoint genes in the two score groups. F Comparison of the differences in StromalScore, ImmuneScore, and ESTIMATEScore between the two score subgroups. ns no significance,* p < 0.05, ** p < 0.01, *** p < 0.001

Full size image

Prediction and validation of immunotherapy efficacy, prediction of targeted chemotherapeutic drugs

To determine the suitability of patients for immunotherapy, we utilized the TIDE score to assess potential immune dysfunctions in the tumor and regional lymph nodes. We found that patients in the low-risk subgroup were more likely to respond positively to immunotherapy (Fig. 10A, B, C). We also examined the sensitivity of three classical chemotherapeutic agents, Lapatinib, Bortezomib, and Elesclomol. In the low-risk group, patients treated with Lapatinib and Bortezomib had higher IC50 values, indicating greater sensitivity to these drugs. Conversely, patients treated with Elesclomol showed higher sensitivity in the high-risk group (Fig. 11A-C).

Fig. 10
figure 10

Prediction and validation of the effect of immunotherapy. A Distribution of ERG scores between responders and non-responders. B ERG score predicts ROC curve for immunotherapy response. C Distribution of TIDE scores between high- and low-risk groups in the TCGA-GBM dataset. D The relationship between ERG scores, ICB response traits, and each stage of the tumor immune cycle. E The plot of the difference in enrichment scores between the high-risk and low-risk groups on the immunotherapy prediction pathway. F The plot of differences between the high-risk and low-risk groups on each step of the cancer-immune cycle. * p < 0.05, ** p < 0.01, *** p < 0.001

Full size image
Fig. 11
figure 11

ERG signature predicts chemotherapy sensitivity. A Lapatinib, B Bortezomib, C Elesclomol

Full size image

Furthermore, we investigated the differences between the two subgroups in predicting immune checkpoint blockade (ICB) response characteristics. We observed that the LR group had higher scores in DNA replication, cell cycle, viral oncogenesis, base excision repair, and p53 signaling pathway, compared to the HR group (Fig. 10E). We also evaluated the relationship between ERG scores and ICB-related positive signals and found a negative correlation between ERG scores and signals such as DNA replication, cell cycle, depletion pathway, mismatch repair, base excision repair, and microRNAs in cancer (Fig. 10D). To assess the biological function of the chemokine system and immunomodulators, we compared the differences in the activity of tumor immune steps between the high- and low-risk groups. We observed that the HR group exhibited upregulated activity in most steps of the tumor immune cycle, including the release of cancer cell antigens (step 1), the presentation of cancer antigens (step 2), priming and activation (step 3), and the entry of immune cells into the tumor (step 4), such as T cell recruitment, CD8 T cell recruitment, Th1 recruitment, DC cell recruitment, and Th22 cell recruitment (Fig. 10F). Additionally, we found a positive correlation between each of these steps in the tumor immune cycle and ERG scores (Fig. 10D).

Discussion

GBMs can be categorized into different subgroups based on their gene expression profiles, which include mutations in isocitrate dehydrogenase (IDH), promoter methylation of O6-methylguanine-DNA methyltransferase (MGMT), and amplification of epidermal growth factor receptor (EGFR), reflecting their histological and morphological characteristics [42, 43]. However, relying solely on tumor size, histologic grade, or individual genetic features to predict prognosis and determine treatment options for GBM patients is insufficient due to the complex and multifactorial nature of GBM development. Thus, there is an urgent need for more accurate models for preclinical selection [29, 44].

Potential targets for cancer therapy include efferocytosis-related genes and pathways, such as phosphatidylserine, TYRO3, MerTK, indoleamine-2,3-dioxygenase 1, membrane-linked protein V, CD 47, TGF-β, and IL-10 [45]. These signals are often upregulated in GBM and associated with poor prognosis [46, 47]. However, treating glioma is complicated by tumor heterogeneity, changes in immune checkpoints, and extensive immunosuppression in the hypoxic microenvironment [48,49,50]. In order to address these challenges, Wu et al. proposed MerTK, an efferocytosis-related receptor, as a potential therapeutic target for glioblastoma [51]. In light of the crucial role that efferocytosis plays in GBM progression, along with immunosuppressive medication and promotion of tumor growth, a new efferocytosis-related scoring system was developed to evaluate the risk and predict personalized therapy [52].

Due to the high heterogeneity of glioblastoma (GBM), identifying various subtypes is often the best prognostic approach for patient intervention (with diverse phenotypes associated with efferocytosis). Despite efferocytosis and tumor development being popular topics in medical research, there is still insufficient literature on the combined effect of efferocytosis-related phenotypes in GBM. In this study, we categorized 549 GBM patients from three cohorts into two ERG clusters. Patients in cluster A had worse overall survival (OS) than those in cluster B, indicating that these efferocytosis-associated genes might affect GBM prognosis. To assess the difference in the tumor immune microenvironment, we compared the enrichment scores of tumor-infiltrating immune cells (TIICs) in the ERG clusters using ssGSEA. Despite having a worse prognosis, cluster A showed higher immune cell infiltration levels. The unique brain immunology contributes to GBM’s distinct tumor microenvironment, where multiple peripheral immune components, including various types of monocytes and lymphocytes, are present in the tumor immune microenvironment (TIME). However, their infiltration rate is significantly lower in gliomas than other tumors. GBM’s tumor-infiltrating lymphocytes (TILs) are low, while the content of CD4 + T cells and CD8 + T cells increases with tumor malignancy [53]. Treg cells also play an important role in the immunosuppressive microenvironment, as a component of the glioma microenvironment [54]. However, in advanced gliomas, TAMs are mainly characterized by an “M2” phenotype, which induces immunosuppressive responses and tumor immune escape [55, 56]. Although natural killer cells (NK cells) are potent innate cytotoxic lymphocytes, in the context of immunosuppression, tumor-associated neutrophils (TANs) play a crucial role. TANs, myeloid-derived suppressor cells (MDSCs), and the combined negative regulation of Treg and NK cells infiltrating in the TIME of GBM are generally considered functionally incompetent [57]. Our findings are consistent with GBM-derived cytokines and chemokines reprogramming TIICs. The distinct immune profiles of the two ERG clusters suggest that some underlying genes need to be identified. Therefore, we extracted differentially expressed genes (DEGs) from both and found that these genes were enriched in cytophagy and immune-related functions, indicating that they could be targeted in immunotherapy.

Using a bioinformatics approach, we developed an ERG signature for GBM based on the PCA algorithm and key DEGs to investigate the impact of efferocytosis-related phenotypes on prognosis. Our ERG signature effectively stratifies GBM patients based on risk and serves as an independent predictor of survival, with lower scores indicating better OS compared to higher scores. In line with previous findings, higher ERG scores are associated with increased tumor-infiltrating immune cells (TIICs) and poorer prognosis. Additionally, higher tumor microenvironment (TME) scores are linked to higher ERG scores. Several bioinformatics studies have shown that high mesenchymal and immune scores are associated with malignancy progression and a very poor prognosis [58]. Given that immunosuppressive cells within the TME can render immunotherapy ineffective, a high TME score is considered a red flag for GBM patients [59]. In addition to the suppressive role of TME, hypoxic conditions can also protect tumors from immune responses by inhibiting natural killer cell and connective tissue cell activity, and promoting immunosuppressive cytokine release and cell function enhancement. Our study shows that higher ERG scores are associated with an immunosuppressive microenvironment, highlighting the OS advantage of patients in lower-scoring groups.

The use of immune checkpoint inhibitors (ICIs) has revolutionized treatment for multiple types of cancer, including melanoma, lung cancer, and kidney cancer, resulting in a significant increase in overall survival for oncology patients [60]. ICIs targeting CTLA4 and PD-1/PDL-1 pathways have improved immune activation, paving the way for new therapies. Although the efficacy of immune checkpoint inhibition therapy for GBM is currently insufficient, ICIs are still the most clinically established form of immunotherapy [61]. Therefore, identifying patients who are likely to benefit from immunotherapy early on is crucial. Recently, a study demonstrated that administering PD-1 inhibitors two weeks before surgery, as a neoadjuvant regimen, improved overall survival in patients with recurrent GBM, compared to postoperative adjuvant therapy. This success supports the theory that PD-1 inhibitors can enhance antitumor immune responses [62, 63]. Nonetheless, current clinical practice lacks specific biomarkers for GBM immunotherapy.

We examined the expression of immune checkpoints (ICs) in high-risk (HR) and low-risk (LR) subgroups. In the HR group, most ICs (PD-1, CTLA-4, IDO, LAG-3, and TIM-3) were highly expressed. The interaction between PD-1 and PD-L1 generates an immune regulatory axis that promotes GBM cell invasion in brain tissue [63]. Elevated PD-L1 in glioma cells binds to PD-1 on tumor-associated macrophages (TAMs) and tumor-infiltrating lymphocytes (TILs), creating a suppressive immune microenvironment and resulting in a poor prognosis for GBM patients [64, 65]. Tumor-derived antigens increase LAG-3 expression, leading to CTL deficiency [66]. TIM-3 controls T-cell depletion by interacting with the ligand Gal-9 and contributes to tumor immune evasion [67]. Overexpression of TIM-3 in GBM is associated with worse prognosis, lower quality of life, and increased malignancy [68, 69]. The cancer immune cycle reflects the immunological response of the human immune system to cancer [70]. In our study of GBM tumor immune cycle and immune checkpoint blockade (ICB) response, we observed a significant positive correlation between ERG scores and ICB-related negative signals, as well as a positive correlation with the suppressive tumor immune cycle. These results further support the presence of immunosuppression and an inflammatory tumor microenvironment in the HR group. With the advent of bioinformatics, various algorithms have been successfully utilized to predict immunotherapy outcomes in tumors [31]. Using the TIDE algorithm, we explored the immunotherapeutic potential of our ERG signature and found significantly higher TIDE scores in the HR group. Lower ERG scores were linked to a better prognosis and higher response rates to immunotherapy.

Unlike inflammatory diseases, tumors employ efferocytosis to foster an immunosuppressive milieu that polarizes macrophages towards the M2 phenotype. This phenotype inhibits anti-tumor immunity, facilitates tissue repair, and stimulates vascular growth, ultimately leading to a dismal prognosis. Augmented expression of efferocytosis-related positive molecules boosts tumor cell survival, migration, invasion, and metastasis. Upon binding to its receptor, phosphatidylserine (PS) impedes the generation of NF-κB and type I IFN, thereby constraining the antitumor immune response.

Conclusions

Despite offering valuable insights, our bioinformatics-based investigation has limitations that must be acknowledged. To confirm our findings, larger prospective studies and more in vivo and in vitro experiments are necessary, particularly for validating the efferocytosis-related signature in a genuine and larger cohort. Nevertheless, our study indicates that the ERG cluster signature may be associated with the prognosis and response to immunotherapy in GBM patients, and can direct future research on efferocytosis. Although the role of efferocytosis in tumor cell apoptosis is still in its preliminary stage of exploration, our study has successfully developed a signature linked to efferocytosis in GBM. In the future, our computational scoring method could help clinicians precisely assess the prognosis and immune status of GBM patients and recognize particular subgroups who may gain from tailored immunotherapy and chemotherapy treatments.

Availability of data and materials

The datasets analyzed for this study were obtained from the UCSC Xena website (https://xenabrowser.net/datapages/) and CGGA dataset (http://www.cgga.org.cn/).

References

  1. Scott JG, Berglund A, Schell MJ, Mihaylov I, Fulp WJ, Yue B, et al. A genome-based model for adjusting radiotherapy dose (GARD): a retrospective, cohort-based study. Lancet Oncol. 2017;18(2):202–11. https://doi.org/10.1016/S1470-2045(16)30648-9. Epub 2016/12/21.

    Article  PubMed  Google Scholar 

  2. Molinaro AM, Taylor JW, Wiencke JK, Wrensch MR. Genetic and molecular epidemiology of adult diffuse glioma. Nat Rev Neurol. 2019;15(7):405–17. https://doi.org/10.1038/s41582-019-0220-2. Epub 2019/06/23.

    Article  PubMed  PubMed Central  Google Scholar 

  3. Ostrom QT, Gittleman H, Fulop J, Liu M, Blanda R, Kromer C, et al. CBTRUS statistical report: primary brain and central nervous system tumors diagnosed in the United States in 2008–2012. Neuro Oncol. 2015;17 Suppl 4(Suppl 4):iv1–62. https://doi.org/10.1093/neuonc/nov189. Epub 2015/10/30.

    Article  PubMed  Google Scholar 

  4. Ostrom QT, Bauchet L, Davis FG, Deltour I, Fisher JL, Langer CE, et al. The epidemiology of glioma in adults: a “state of the science” review. Neuro Oncol. 2014;16(7):896–913. https://doi.org/10.1093/neuonc/nou087. Epub 2014/05/21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Richards LM, Whitley OKN, MacLeod G, Cavalli FMG, Coutinho FJ, Jaramillo JE, et al. Gradient of developmental and injury response transcriptional states defines functional vulnerabilities underpinning glioblastoma heterogeneity. Nat Cancer. 2021;2(2):157–73. https://doi.org/10.1038/s43018-020-00154-9. Epub 2022/02/06.

    Article  CAS  PubMed  Google Scholar 

  6. Niu X, Sun J, Meng L, Fang T, Zhang T, Jiang J, et al. A five-IncRNAs signature-derived risk score based on TCGA and CGGA for glioblastoma: potential prospects for treatment evaluation and prognostic prediction. Front Oncol. 2020;10:590352. https://doi.org/10.3389/fonc.2020.590352. Epub 2021/01/05.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Wang Z, Gao L, Guo X, Feng C, Lian W, Deng K, et al. Development of a nomogram with alternative splicing signatures for predicting the prognosis of glioblastoma: a study based on large-scale sequencing data. Front Oncol. 2020;10:1257. https://doi.org/10.3389/fonc.2020.01257. Epub 2020/08/15.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Banerjee HN, Bartlett V, Krauss C, Aurelius C, Johnston K, Hedley J, et al. Efferocytosis and the story of “Find Me,” “Eat Me,” and “Don’t Eat Me” signaling in the tumor microenvironment. Adv Exp Med Biol. 2021;1329:153–62. https://doi.org/10.1007/978-3-030-73119-9_8. Epub 2021/10/20.

    Article  CAS  PubMed  Google Scholar 

  9. Sachet M, Liang YY, Oehler R. The immune response to secondary necrotic cells. Apoptosis. 2017;22(10):1189–204. https://doi.org/10.1007/s10495-017-1413-z. Epub 2017/09/02.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Gardai SJ, McPhillips KA, Frasch SC, Janssen WJ, Starefeldt A, Murphy-Ullrich JE, et al. Cell-surface calreticulin initiates clearance of viable or apoptotic cells through trans-activation of LRP on the phagocyte. Cell. 2005;123(2):321–34. https://doi.org/10.1016/j.cell.2005.08.032. Epub 2005/10/22.

    Article  CAS  PubMed  Google Scholar 

  11. Zhou Y, Yao Y, Deng Y, Shao A. Regulation of efferocytosis as a novel cancer therapy. Cell Commun Signal. 2020;18(1):71. https://doi.org/10.1186/s12964-020-00542-9. Epub 2020/05/07.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Lantz C, Radmanesh B, Liu E, Thorp EB, Lin J. Single-cell RNA sequencing uncovers heterogenous transcriptional signatures in macrophages during efferocytosis. Sci Rep. 2020;10(1):14333. https://doi.org/10.1038/s41598-020-70353-y. Epub 2020/09/02.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Voll RE, Herrmann M, Roth EA, Stach C, Kalden JR, Girkontaite I. Immunosuppressive effects of apoptotic cells. Nature. 1997;390(6658):350–1. https://doi.org/10.1038/37022. Epub 1997/12/06.

    Article  CAS  PubMed  Google Scholar 

  14. Liang YY, Schwarzinger I, Simonitsch-Klupp I, Agis H, Oehler R. Impaired efferocytosis by monocytes in multiple myeloma. Oncol Lett. 2018;16(1):409–16. https://doi.org/10.3892/ol.2018.8620. Epub 2018/06/22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Yan D, Parker RE, Wang X, Frye SV, Earp HS 3rd, DeRyckere D, et al. MerTK promotes resistance to irreversible EGFR tyrosine kinase inhibitors in non-small cell lung cancers expressing wild-type EGFR family members. Clin Cancer Res. 2018;24(24):6523–35. https://doi.org/10.1158/1078-0432.CCR-18-0040. Epub 2018/09/09.

    Article  CAS  PubMed  Google Scholar 

  16. Pujol-Autonell I, Ampudia RM, Planas R, Marin-Gallen S, Carrascal J, Sanchez A, et al. Efferocytosis promotes suppressive effects on dendritic cells through prostaglandin E2 production in the context of autoimmunity. PLoS ONE. 2013;8(5):e63296. https://doi.org/10.1371/journal.pone.0063296. Epub 2013/05/22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Tan X, Zhang Z, Yao H, Shen L. Tim-4 promotes the growth of colorectal cancer by activating angiogenesis and recruiting tumor-associated macrophages via the PI3K/AKT/mTOR signaling pathway. Cancer Lett. 2018;436:119–28. https://doi.org/10.1016/j.canlet.2018.08.012. Epub 2018/08/18.

    Article  CAS  PubMed  Google Scholar 

  18. Fadeel B. Plasma membrane alterations during apoptosis: role in corpse clearance. Antioxid Redox Signal. 2004;6(2):269–75. https://doi.org/10.1089/152308604322899332. Epub 2004/03/18.

    Article  CAS  PubMed  Google Scholar 

  19. Gheibi Hayat SM, Bianconi V, Pirro M, Sahebkar A. Efferocytosis: molecular mechanisms and pathophysiological perspectives. Immunol Cell Biol. 2019;97(2):124–33. https://doi.org/10.1111/imcb.12206. Epub 2018/09/20.

    Article  PubMed  Google Scholar 

  20. Werfel TA, Elion DL, Rahman B, Hicks DJ, Sanchez V, Gonzales-Ericsson PI, et al. Treatment-induced tumor cell apoptosis and secondary necrosis drive tumor progression in the residual tumor microenvironment through MerTK and IDO1. Cancer Res. 2019;79(1):171–82. https://doi.org/10.1158/0008-5472.CAN-18-1106. Epub 2018/11/11.

    Article  CAS  PubMed  Google Scholar 

  21. Turner WA, Fadel HE, Krauss JS Jr. Detection and quantitation of fetomaternal hemorrhage. South Med J. 1986;79(5):571–5. https://doi.org/10.1097/00007611-198605000-00012. Epub 1986/05/01.

    Article  CAS  PubMed  Google Scholar 

  22. Liang X, Luo M, Shao B, Yang JY, Tong A, Wang RB, et al. Phosphatidylserine released from apoptotic cells in tumor induces M2-Like macrophage polarization through the PSR-STAT3-JMJD3 axis. Cancer Commun (Lond). 2022;42(3):205–22. https://doi.org/10.1002/cac2.12272. Epub 2022/02/23.

    Article  PubMed  Google Scholar 

  23. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012;28(6):882–3. https://doi.org/10.1093/bioinformatics/bts034. Epub 2012/01/20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Wilkerson MD, Hayes DN. Consensusclusterplus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26(12):1572–3. https://doi.org/10.1093/bioinformatics/btq170. Epub 2010/04/30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Zhang H, Zhang N, Wu W, Zhou R, Li S, Wang Z, et al. Machine learning-based tumor-infiltrating immune cell-associated IncRNAs for predicting prognosis and immunotherapy response in patients with glioblastoma. Brief Bioinform. 2022;23(6):386. https://doi.org/10.1093/bib/bbac386. Epub 2022/09/23.

    Article  CAS  Google Scholar 

  26. Zhang N, Zhang H, Wu W, Zhou R, Li S, Wang Z, et al. Machine learning-based identification of tumor-infiltrating immune cell-associated IncRNAs for improving outcomes and immunotherapy responses in patients with low-grade glioma. Theranostics. 2022;12(13):5931–48. https://doi.org/10.7150/thno.74281. Epub 2022/08/16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. Clusterprofiler 4.0: a universal enrichment tool for interpreting omics data. Innovation (Camb). 2021;2(3):100141. https://doi.org/10.1016/j.xinn.2021.100141. Epub 2021/09/25.

    Article  CAS  PubMed  Google Scholar 

  28. Zhao S, Zhang X, Gao F, Chi H, Zhang J, Xia Z, et al. Identification of copper metabolism-related subtypes and establishment of the prognostic model in ovarian cancer. Front Endocrinol (Lausanne). 2023;14:1145797. https://doi.org/10.3389/fendo.2023.1145797. Epub 2023/03/24.

    Article  PubMed  Google Scholar 

  29. Zhao S, Chi H, Yang Q, Chen S, Wu C, Lai G, et al. Identification and validation of neurotrophic factor-related gene signatures in glioblastoma and Parkinson’s disease. Front Immunol. 2023;14:1090040. https://doi.org/10.3389/fimmu.2023.1090040. Epub 2023/02/25.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Zhao H, Chen Y, Shen P, Gong L. Identification of immune cell infiltration landscape and their prognostic significance in uveal melanoma. Front Cell Dev Biol. 2021;9:713569. https://doi.org/10.3389/fcell.2021.713569. Epub 2021/09/14.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018;24(10):1550–8. https://doi.org/10.1038/s41591-018-0136-1. Epub 2018/08/22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Xu L, Deng C, Pang B, Zhang X, Liu W, Liao G, et al. Tip: a web server for resolving tumor immunophenotype profiling. Cancer Res. 2018;78(23):6575–80. https://doi.org/10.1158/0008-5472.CAN-18-0689. Epub 2018/08/30.

    Article  CAS  PubMed  Google Scholar 

  33. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. Tgfbeta attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018;554(7693):544–8. https://doi.org/10.1038/nature25501. Epub 2018/02/15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Aran D, Hu Z, Butte AJ. Xcell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. 2017;18(1):220. https://doi.org/10.1186/s13059-017-1349-1. Epub 2017/11/17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, et al. Timer2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 2020;48(W1):W509–14. https://doi.org/10.1093/nar/gkaa407. Epub 2020/05/23.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Racle J, de Jonge K, Baumgaertner P, Speiser DE, Gfeller D. Simultaneous enumeration of cancer and immune cell types from bulk tumor gene expression data. Elife. 2017;6:e26476. https://doi.org/10.7554/eLife.26476. Epub 2017/11/14.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol Biol. 2018;1711:243–59. https://doi.org/10.1007/978-1-4939-7493-1_12. Epub 2018/01/19.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Xie J, Zhang J, Tian W, Zou Y, Tang Y, Zheng S, et al. The pan-cancer multi-omics landscape of FOXO family relevant to clinical outcome and drug resistance. Int J Mol Sci. 2022;23(24):15647. https://doi.org/10.3390/ijms232415647. Epub 2022/12/24.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Chen X, Jiang X, Wang H, Wang C, Wang C, Pan C, et al. DNA methylation-regulated SNX20 overexpression correlates with poor prognosis, immune cell infiltration, and low-grade glioma progression. Aging (Albany NY). 2022;14(12):5211–22. https://doi.org/10.18632/aging.204144. Epub 2022/07/01.

    Article  CAS  PubMed  Google Scholar 

  40. Yao Z, Zhu G, Too J, Duan M, Wang Z. Feature selection of omic data by ensemble swarm intelligence based approaches. Front Genet. 2021;12:793629. https://doi.org/10.3389/fgene.2021.793629. Epub 2022/03/31.

    Article  PubMed  Google Scholar 

  41. Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One. 2014;9(9):e107468. https://doi.org/10.1371/journal.pone.0107468. Epub 2014/09/18.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Aldape K, Zadeh G, Mansouri S, Reifenberger G, von Deimling A. Glioblastoma: pathology, molecular mechanisms and markers. Acta Neuropathol. 2015;129(6):829–48. https://doi.org/10.1007/s00401-015-1432-1. Epub 2015/05/07.

    Article  CAS  PubMed  Google Scholar 

  43. Zhao S, Ji W, Shen Y, Fan Y, Huang H, Huang J, et al. Expression of hub genes of endothelial cells in glioblastoma-A prognostic model for GBM patients integrating single-cell RNA sequencing and bulk RNA sequencing. BMC Cancer. 2022;22(1):1274. https://doi.org/10.1186/s12885-022-10305-z. Epub 2022/12/07.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Zhao S, Chi H, Ji W, He Q, Lai G, Peng G, et al. A bioinformatics-based analysis of an anoikis-related gene signature predicts the prognosis of patients with low-grade gliomas. Brain Sci. 2022;12(10):1349. https://doi.org/10.3390/brainsci12101349. Epub 2022/10/28.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Tajbakhsh A, Gheibi Hayat SM, Movahedpour A, Savardashtaki A, Loveless R, Barreto GE, et al. The complex roles of efferocytosis in cancer development, metastasis, and treatment. Biomed Pharmacother. 2021;140:111776. https://doi.org/10.1016/j.biopha.2021.111776. Epub 2021/06/02.

    Article  CAS  PubMed  Google Scholar 

  46. Wium M, Paccez JD, Zerbini LF. The dual role of tam receptors in autoimmune diseases and cancer: an overview. Cells. 2018;7(10):166. https://doi.org/10.3390/cells7100166. Epub 2018/10/17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Onken J, Vajkoczy P, Torka R, Hempt C, Patsouris V, Heppner FL, et al. Phospho-AXL is widely expressed in glioblastoma and associated with significant shorter overall survival. Oncotarget. 2017;8(31):50403–14. https://doi.org/10.18632/oncotarget.18468. Epub 2017/09/09.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Malta TM, de Souza CF, Sabedot TS, Silva TC, Mosella MS, Kalkanis SN, et al. Glioma CpG island methylator phenotype (G-CIMP): biological and clinical implications. Neuro Oncol. 2018;20(5):608–20. https://doi.org/10.1093/neuonc/nox183. Epub 2017/10/17.

    Article  CAS  PubMed  Google Scholar 

  49. Chen R, Cohen AL, Colman H. Targeted therapeutics in patients with high-grade gliomas: past, present, and future. Curr Treat Options Oncol. 2016;17(8):42. https://doi.org/10.1007/s11864-016-0418-0. Epub 2016/06/24.

    Article  PubMed  Google Scholar 

  50. Colwell N, Larion M, Giles AJ, Seldomridge AN, Sizdahkhani S, Gilbert MR, et al. Hypoxia in the glioblastoma microenvironment: shaping the phenotype of cancer stem-like cells. Neuro Oncol. 2017;19(7):887–96. https://doi.org/10.1093/neuonc/now258. Epub 2017/03/25.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Wu J, Frady LN, Bash RE, Cohen SM, Schorzman AN, Su YT, et al. Mertk as a therapeutic target in glioblastoma. Neuro Oncol. 2018;20(1):92–102. https://doi.org/10.1093/neuonc/nox111. Epub 2017/06/13.

    Article  CAS  PubMed  Google Scholar 

  52. Che Mat MF, Abdul Murad NA, Ibrahim K, Mohd Mokhtar N, Wan Ngah WZ, Harun R, et al. Silencing of PROS1 induces apoptosis and inhibits migration and invasion of glioblastoma multiforme cells. Int J Oncol. 2016;49(6):2359–66. https://doi.org/10.3892/ijo.2016.3755. Epub 2016/11/15.

    Article  CAS  PubMed  Google Scholar 

  53. Gieryng A, Pszczolkowska D, Walentynowicz KA, Rajan WD, Kaminska B. Immune microenvironment of gliomas. Lab Invest. 2017;97(5):498–518. https://doi.org/10.1038/labinvest.2017.19. Epub 2017/03/14.

    Article  CAS  PubMed  Google Scholar 

  54. Wainwright DA, Dey M, Chang A, Lesniak MS. Targeting Tregs in malignant brain cancer: overcoming IDO. Front Immunol. 2013;4:116. https://doi.org/10.3389/fimmu.2013.00116. Epub 2013/05/31.

    Article  PubMed  PubMed Central  Google Scholar 

  55. de Groot AE, Pienta KJ. Epigenetic control of macrophage polarization: implications for targeting tumor-associated macrophages. Oncotarget. 2018;9(29):20908–27. https://doi.org/10.18632/oncotarget.24556. Epub 2018/05/15.

    Article  PubMed  PubMed Central  Google Scholar 

  56. Liu Y, Cao X. The origin and function of tumor-associated macrophages. Cell Mol Immunol. 2015;12(1):1–4. https://doi.org/10.1038/cmi.2014.83. Epub 2014/09/16.

    Article  CAS  PubMed  Google Scholar 

  57. Poli A, Wang J, Domingues O, Planaguma J, Yan T, Rygh CB, et al. Targeting glioblastoma with NK cells and mAb against NG2/CSPG4 prolongs animal survival. Oncotarget. 2013;4(9):1527–46. https://doi.org/10.18632/oncotarget.1291. Epub 2013/10/16.

    Article  PubMed  PubMed Central  Google Scholar 

  58. Zhang C, Cheng W, Ren X, Wang Z, Liu X, Li G, et al. Tumor purity as an underlying key factor in glioma. Clin Cancer Res. 2017;23(20):6279–91. https://doi.org/10.1158/1078-0432.CCR-16-2598. Epub 2017/07/30.

    Article  CAS  PubMed  Google Scholar 

  59. Alghamri MS, Banerjee K, Mujeeb AA, Mauser A, Taher A, Thalla R, et al. Systemic delivery of an adjuvant CXCR4-CXCL12 signaling inhibitor encapsulated in synthetic protein nanoparticles for glioma immunotherapy. ACS Nano. 2022;16(6):8729–50. https://doi.org/10.1021/acsnano.1c07492. Epub 2022/05/27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Cheng F, Liang H, Butte AJ, Eng C, Nussinov R. Personal mutanomes meet modern oncology drug discovery and precision health. Pharmacol Rev. 2019;71(1):1–19. https://doi.org/10.1124/pr.118.016253. Epub 2018/12/14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Qi Y, Liu B, Sun Q, Xiong X, Chen Q. Immune checkpoint targeted therapy in glioma: status and hopes. Front Immunol. 2020;11:578877. https://doi.org/10.3389/fimmu.2020.578877. Epub 2020/12/18.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Cloughesy TF, Mochizuki AY, Orpilla JR, Hugo W, Lee AH, Davidson TB, et al. Neoadjuvant anti-PD-1 immunotherapy promotes a survival benefit with intratumoral and systemic immune responses in recurrent glioblastoma. Nat Med. 2019;25(3):477–86. https://doi.org/10.1038/s41591-018-0337-7. Epub 2019/02/12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Berghoff AS, Kiesel B, Widhalm G, Rajky O, Ricken G, Wohrer A, et al. Programmed death ligand 1 expression and tumor-infiltrating lymphocytes in glioblastoma. Neuro Oncol. 2015;17(8):1064–75. https://doi.org/10.1093/neuonc/nou307. Epub 2014/10/31.

    Article  CAS  PubMed  Google Scholar 

  64. Zhou J, Pei X, Yang Y, Wang Z, Gao W, Ye R, et al. Orphan nuclear receptor TLX promotes immunosuppression via its transcriptional activation of PD-L1 in glioma. J Immunother Cancer. 2021;9(4):e001937. https://doi.org/10.1136/jitc-2020-001937. Epub 2021/04/17.

    Article  PubMed  PubMed Central  Google Scholar 

  65. Gordon SR, Maute RL, Dulken BW, Hutter G, George BM, McCracken MN, et al. PD-1 expression by tumour-associated macrophages inhibits phagocytosis and tumour immunity. Nature. 2017;545(7655):495–9. https://doi.org/10.1038/nature22396. Epub 2017/05/18.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Andrews LP, Marciscano AE, Drake CG, Vignali DA. LAG3 (CD223) as a cancer immunotherapy target. Immunol Rev. 2017;276(1):80–96. https://doi.org/10.1111/imr.12519. Epub 2017/03/05.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Das M, Zhu C, Kuchroo VK. Tim-3 and its role in regulating anti-tumor immunity. Immunol Rev. 2017;276(1):97–111. https://doi.org/10.1111/imr.12520. Epub 2017/03/05.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Han S, Feng S, Xu L, Shi W, Wang X, Wang H, et al. Tim-3 on peripheral CD4(+) and CD8(+) T cells is involved in the development of glioma. DNA Cell Biol. 2014;33(4):245–50. https://doi.org/10.1089/dna.2013.2306. Epub 2014/02/12.

    Article  CAS  PubMed  Google Scholar 

  69. Li G, Wang Z, Zhang C, Liu X, Cai J, Wang Z, et al. Molecular and clinical characterization of Tim-3 in glioma through 1,024 samples. Oncoimmunology. 2017;6(8):e1328339. https://doi.org/10.1080/2162402X.2017.1328339. Epub 2017/09/19.

    Article  PubMed  PubMed Central  Google Scholar 

  70. Chen DS, Mellman I. Oncology meets immunology: the cancer-immunity cycle. Immunity. 2013;39(1):1–10. https://doi.org/10.1016/j.immuni.2013.07.012. Epub 2013/07/31.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We thank all the participants involved in this study.

Funding

This work was supported by Wuxi Taihu Lake Talent Plan, Supports for Leading Talents in Medical and Health Profession(2020THRC-DJ-SNW).

Author information

Author notes

  1. Songyun Zhao, Qi Wang and Yuankun Liu contributed equally to this work.

Authors and Affiliations

  1. Department of Neurosurgery, Affiliated Wuxi People’s Hospital of Nanjing Medical University, Wuxi, China

    Songyun Zhao, Yuankun Liu, Wei Ji & Chao Cheng

  2. Department of Gastroenterology, Affiliated Hospital of Jiangsu University, Zhenjiang, China

    Qi Wang

  3. Department of Lung Cancer Surgery, Tianjin Medical University Cancer Institute and Hospital, Tianjin, China

    Pengpeng Zhang

  4. The First Affiliated Hospital of Nanjing Medical University, Nanjing, China

    Pengpeng Zhang & Jiaheng Xie

  5. Department of Plastic Surgery, Xiangya Hospital, Central South University, Changsha, China

    Jiaheng Xie

Authors
  1. Songyun Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Qi Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Yuankun Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Pengpeng Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Wei Ji
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Jiaheng Xie
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Chao Cheng
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

SZ and QW conceived the study. SZ, QW and YL drafted the manuscript. PZ performed the literature search and collected the data. WJ and PZ analyzed and visualized the data. CC and JX helped with the final revision of this manuscript. All authors contributed to the article and approved the submitted version.

Corresponding authors

Correspondence to Jiaheng Xie or Chao Cheng.

Ethics declarations

Ethics approval and consent to participate

This article does not contain any studies with human participants or animals.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Supplementary Figure 1.

Consensus unsupervised clustering of 549 GBM patients based on 8 ERGs. Supplementary Figure 2. Consensus unsupervised clustering of 549 GBM patients based on 473 DEGs. Supplementary Figure 3. Differences in tumor immune microenvironment between low and high-scoring subgroups in the CGGA cohort. Supplementary Table 1. 22 tumor-related ERGs were obtained from literature search results. Supplementary Table 2. siRNA and primer sequence information. Supplementary Table 3. Queue information was obtained for the two ERG clusters. Supplementary Table 4. Differential genes between two ERG clusters. Supplementary Table 5. Results after univariate Cox analysis of differential genes between two ERG clusters. Supplementary Table 6. Detailed cohort information for the 2 ERG gene clusters by the consensus clustering algorithm. Supplementary Table 7. Applying the “Boruta” package to screen for important genes.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhao, S., Wang, Q., Liu, Y. et al. Interaction, immune infiltration characteristics and prognostic modeling of efferocytosis-related subtypes in glioblastoma. BMC Med Genomics 16, 248 (2023). https://doi.org/10.1186/s12920-023-01688-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12920-023-01688-4

Keywords

BMC Medical Genomics

ISSN: 1755-8794

Contact us