Identification and validation of a ferroptosis-related gene to predict survival outcomes and the immune microenvironment in lung adenocarcinoma

Obtaining DEGs of LUAD

A total of 382 FRGs were downloaded from the public database FerrDb, and noncoding genes were deleted from the data set. Difference analysis showed that a total of 70 FRGs were differentially expressed between tumors and normal tissues; the characteristics of these FRGs in the TCGA-LUAD cohort are shown in Fig. 1A. In addition, all LUAD samples were clustered by consensus clustering based on 70 differentially expressed FRGs, and the best clustering effect was achieved when K = 4 (Additional file 1: Figure S1A–D). Surprisingly, significant differences were observed in OS between the four clusters (Additional file 1: Fig. S1E).

Fig. 1figure 1

Construction of 7–gene Risk Signature. A The heatmap of FRDEGs between normal and tumor samples. B The forest plot of univariate Cox regression. C Differential expressed volcanic maps. The positions of seven key FRGs in the volcanic map. D, E Forest plots show independent prognostic roles of RiskScore and other clinical features

Construction and validation of prognostic model

As shown in the forest map, 15 DEGs associated with LUAD prognosis were identified by the R software package ‘survival’ in the TCGA cohort (Fig. 1B). Subsequently, seven FRGs with independent prognostic significance for patients with LUAD were obtained by multivariate Cox regression analysis on the 15 appealed DEGs; the locations of these FRGs in the differential volcanic map are shown in Fig. 1C. The Kaplan–Meier curve revealed the survival differences between the seven FRGs at different levels (Additional file 1: Fig. S1F–L). Ultimately, we established a seven-gene risk signature using the multivariate Cox regression coefficient, modeled as follows: Risk score of each patient = (0.262297470303534 * expression level of ALOX12B) + (0.0055627086031226 * expression level of DDIT4) + (0.00646442374221167 * expression level of SLC7A5) + (− 0.0131510931059613 * expression level of TRIB3) + (− 0.0300191156801198 * expression level of IL33) + (0.015384279700975 * expression level of RRM2) + (0.0046423481143834 * expression level of CAV1). According to the risk index ranking of 504 patients in TCGA-LUAD, patients were divided into a high-risk group and a low-risk group, with the median risk index as the cutoff between high and low risk. An independent prognostic analysis of the clinical characteristics of LUAD (age, gender, stage, and risk score) was conducted (Fig. 1D, E) after screening out the clinical information. The risk score and clinical stage were independent factors affecting the prognosis of LUAD. The clinical correlation chart in Additional file 3: Figure S3A shows the proportion of patients with high and low risk scores at each stage. Further analysis showed that the OS rate of TCGA-LUAD participants in the high-risk group was significantly lower than the OS rate of those in the low-risk group (P < 0.001) (Fig. 2A). Figure 2B shows that a difference in the OS rate was also present between the low-risk and high-risk groups in the GSE30219 cohort; the OS rate of the high-risk group was significantly lower than that of the low-risk group (P < 0.001). The ROC curve (Fig. 2C) shows that our risk signature had high sensitivity and specificity in predicting the OS rate of patients with LUAD. The AUCs were 0.707, 0.685, and 0.673 for 1-, 2-, and 3-year OS, respectively. We then added clinical features, including TIDE, tumor inflammation signature (TIS), age, and stage, to compare the accuracy of OS prediction. As shown in Fig. 2D, the risk index was the best predictor of OS. The DCA (Fig. 2E) showed that the net benefit of the risk index was the largest. Finally, a nomogram based on the two independent prognostic factors was constructed to quantify the survival probability of patients with LUAD in 1, 3, and 5 years (Fig. 2F).

Fig. 2figure 2

Verification of 7-gene Risk Signature and Construction of Nomo Graph. A, B Survival analysis between the high- and low-risk groups based on OS. C ROC curve of the 7-gene signature model. D 1-year survival rate predicted by 7-gene signature and other clinical features was compared in ROC curve. E Decision curve. F Nomogram for predicting 1,3,5-year survival rates in LUAD patients

Enrichment analysis

We conducted an enrichment analysis to explore the potential biological characteristics of 70 DEGs. The 25 pathways shown in Additional file 2: Figure S2A were enriched in the KEGG enrichment analysis. The more prominent pathways, represented by red bubbles, included fluid shear stress and atherosclerosis, HIF-1 signaling pathways, lipid and atherosclerosis, and ferroptosis. GO enrichment analysis showed that the biological processes (including responses to oxidative stress and metal ions) and the cell components (including the apical plasma membrane, cell apex, and oxidoreductase complex) were significantly enriched (Additional file 2: Figure S2B). Furthermore, GSEA enrichment analysis showed that the high- and low-risk groups were enriched incompletely different functions or pathways. Among high-risk patients, small cell pulmonary carcinoma, melanoma, cancer pathways, thep53 signaling pathway, and pathogenic Escherichia coli infection were significantly enriched (Additional file 2: Figure S2C).

Difference of immune related indexes between LUAD patients

We calculate the StromalScore, ImmuneScore, and EstimateScore for each LUAD sample usingthe estimate algorithm. The t-tests revealed significant differences in StromalScore, ImmuneScore, and EstimateScore between high- and low-risk groups (Fig. 3A–C). To explore the influence of different levels of immune-related scores on the prognosis of patients, the StromalScore, ImmuneScore, and EstimateScore of all patients with LUAD were divided into high- and low-score groups by average. The OS rates of the two groups were compared using the Kaplan–Meier method; Fig. 3D shows that a higher EstimateScore represented a better prognosis. The same was true for ImmuneScore and StromalScore (Additional file 5: Figure S5A, B). Significant differences in risk scores were also observed among the four clusters (Additional file 3: Figure S3C).

Fig. 3figure 3

Difference of immune related indexes between LUAD patients. A, B, C Violin Chart. Differences in StromalScore, ImmuneScore and EstimateScore between high- and low-risk patients. D Different OS based on different levels of ESTIMATEScore. E, F, G Differences in StromalScore, ImmuneScore and EstimateScore among the four FRGs clusters. H Proportion of high- and low-infiltration of immunocytes in four FRGs clusters

LUAD immune subtype based on FRG clustering

In addition, the t-test revealed significant differences in immune-related scores among the four clusters obtained by the consensus clustering analysis (Fig. 3E–G), suggesting that the four clusters obtained for70 differentially expressed FRGs may serve as immune subtypes. To explore the relationship between the aforementioned four clusters and the two groups clustered according to the level of immunocyte infiltration, we drew a clinical correlation heat map, shown in Fig. 3H. Significant differences were observed in immunocyte infiltration among the four clusters. For example, the samples with high immunocyte infiltration in the composition of C3 were significantly more than those with low immunocyte infiltration, whereas the samples with low immunocyte infiltration in the C2 and C4 clusters were more. The mulberry diagram in Fig. 4A shows the corresponding relationship between immune infiltration level, FRG cluster, risk level, and other clinical features. The above results show that the four clusters identified by the clustering of 70 FRGs may represent the immune subtypes of LUAD.

Fig. 4figure 4

Difference and correlation of immunocyte infiltration in LUAD patients. A Alluvial diagram of the relationship between RiskScore and several features. B Immune cell infiltration heatmap, 7 algorithms to assess immune cell infiltration level among the high-and low-risk groups. C Correlation heatmap, correlation of infiltration levels of immune cells in LUAD

Further analysis of immune cell infiltration in LUAD tumor microenvironment

The CIBERSORT algorithm was used to quantify the levels of 22 immune cells in each LUAD sample, and at-test was used to compare the levels of immune cells between the high-risk and low-risk patients. As shown in Fig. 4B, significant differences were observed in the levels of B cell memory, CD4 memory T cells, activated CD4 memory T cells, regulatory T cells (Tregs), resting NK cells, monocytes, M0 macrophages, resting dendritic cells, resting mast cells, and activated mast cells between high-risk and low-risk patients. The relative proportion of 22 immune cells in all LUAD samples is shown by a barplot (Additional file 3: Figure S3D). When Kaplan–Meier curves were used to compare the effects of immunocytes with different infiltration levels on the OS rates, we found that different infiltration levels of 12 typesof cells, such as activated CD4 memory T cells and resting mast cells, were associated with different OS rates (Additional file 4: Figure S4A–L). To learn more about the differences between immune cells in high- and low-risk patients, we analyzed the data using six other algorithms (TIMER, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC); Fig. 4B shows the different infiltration levels of immunocytes in each database.

Relationship between immunocyte infiltration and clinical characteristics

The correlation heat map shown in Fig. 4C reveals the relationship between the immune cells; most modules are blue, suggesting that the immune cell infiltration in the LUAD microenvironment is negatively correlated. To understand the correlation between immune cell infiltration and clinical characteristics, we calculated the immune cell infiltration level of all LUAD samples in TCGA using the ssGSEA algorithm and then downloaded the clinical information for all LUAD samples from TCGA using the R package ‘TCGAbiolinks’. The 580 LUAD samples with clinical information were clustered into two groups according to the level of immune cell infiltration and visualized by a heat map (Fig. 5A). The left portion, dominated by blue, represents the samples with less immunocyte infiltration, whereas the right portion, dominated by red, represents the samples with higher immunocyte infiltration. The clinical information, such as survival status, stage, gender, tumor location, and related driving gene mutation state, corresponded to the heat map below. In addition, to better understand the immune differences between high-risk and low-risk groups, we quantified 29 immune-related functions using the ssGSEA algorithm. As shown in Fig. 5B, 16 immune-related functions or pathways, such as Adcs and CCR, showed significant differences. When the Kaplan–Meier curve was used to compare the effects of different levels of immune-related functions on OS rates, we found that various levels of multiple immune-related functions or pathways predicted different OS rates. Figure 5C shows that varying levels of mast cell-related immune functions or pathways corresponded to a different prognosis for patients with LUAD. After univariate Cox regression analysis, we identified activated mast cells and activated CD4 memory T cells as prognostic immune cells of LUAD; these two cell types, which had high levels of infiltration, were adverse prognostic factors for LUAD (Fig. 6A). Further analysis confirmed that the expression levels of seven risk genes significantly affected the infiltration level of immunocytes in the LUAD microenvironment, and the infiltration level of activated mast cells, which was most related to the prognosis of LUAD, was affected by the expression levels of RRM2 and IL33 (Fig. 6B). The expression of RRM2 not only regulates the infiltration of activated mast cells but is also related to the infiltration of 19 immunocytes, such as CD4 memory T cells, dendritic cells, and macrophages (Fig. 6C).

Fig. 5figure 5

Analysis of Immunocytes Clustering and immune-related functions in LUAD Patients. A Clustering of immunocytes infiltration. TCGA-LUAD cohort was divided into two groups with high-and low-infiltration levels according to 24 immunocytes. B The difference of immune-related function scores between high- and low-risk groups. C Different infiltration levels of Mast-cell represent different LUAD OS

Fig. 6figure 6

Analysis of the relationship between FRGs and prognosis-related immunocytes and immunotherapy-related indicators. A Immunocytes which play independent prognostic role in LUAD patients. B Correlation between the expression levels of 7 hug FRGs and the infiltration levels of immunocytes. C Immunocytes whose infiltration level correlated with RRM2 expression. D Difference in TMB among high- and low-risk groups; E, F, G, H Violin Chart. Differences in Four Immunotherapy Related scores between high- and low-risk Groups

Risk signature and immunotherapy

To explore the role of the risk signature in guiding immunotherapy, we identified16 patients with pulmonary carcinoma receiving anti-PDL1 immunotherapy in the GSE126044 cohort. Additional file 5: Figure S5G shows the differences in immunocyte infiltration and response to anti-PD-L1 immunotherapy between high-risk and low-risk patients. High-risk patients had lower immunocyte infiltration levels, but their response rates to anti-PD-L1 immunotherapy were higher. In addition, we found that the expression of RRM2 was higher in patients who had no response to anti-PDL1 treatment (P = 0.052) (Additional file 5: Figure S5F). We then downloaded the TIDE and MSI scores of the LUAD cohort from the TIDE database, and we downloaded the TMB data of patients with LUAD using the R packages ‘TCGAbiolinks’ and ‘maftools’. The Wilcoxon rank sum test revealed significant differences between the five immune therapy-related scores described above in high-risk and low-risk patients (Fig. 6D–H). The patients with LUAD with high TMB had worse OS (Additional file 5: Figure S5D). Subsequent correlation analysis showed that the risk score was positively correlated with TMB (R = 0.23, P = 0.000018) (Additional file 5: Figure S5C). The correlation heat map in Additional file 3: Figure S3B shows the distribution of clinical features in the high-risk and low-risk groups. The distribution of TMB, gender, stage, T, and N was significantly different in the high-risk and low-risk groups.

To study the relationship between the seven-gene risk signature and immunotherapy, we counted the mutations of seven key genes in patients with LUAD. As shown in the waterfall diagram in Additional file 5: Figure S5E, ALOX12B had the highest the mutation probability (18%), followed by RRM2 and IL33, and missense mutation was the most common form. To further explore the relationship between the expression of the seven risk genes and immunocyte infiltration, we conducted a correlation analysis, as shown in Fig. 7A, B. The expression level of RRM2 positively regulated the infiltration abundance of activated mast cells (R = 0.12) and activated CD4 memory T cells (R = 0.41), and the infiltration level of resting mast cells was significantly negatively correlated with the expression level of ALOX12B (R =  −0.13) (Additional file 3: Figure S3F). When the ssGSEA algorithm was used to analyze the scores of immune-related functions or pathways, the scores of the checkpoints between high- and low-risk LUAD samples were significantly different. The checkpoints were divided into high and low groups using the average, and the Kaplan–Meier curve showed significant differences in survival (Fig. 7C). Next, we analyzed the expression levels of 49 common immune checkpoints in high-risk versus low-risk patients, as shown in the box plot in Fig. 7D. The expression levels of ALK, ROS1, CD44, and other 27 checkpoints were different. As shown in Fig. 7E, F, the expression levels of ALK and ROS1 decreased as the risk score increased. Furthermore, the correlation analysis showed that ALK and ROS1 were significantly associated with key risk genes. The expression of ROS1 was positively correlated with that ofIL33 and negatively correlated with that ofSLC7A5, TRIB3, and RRM2 (Fig. 7G). The expression of ALK was positively correlated with the expression of IL33 and CAV1and negatively correlated with the expression of ALOX12B, SLC7A5, and TRIB3 (Fig. 7H). The expression of two additional checkpoints, PD-1 and PD-L1, were also significantly correlated with multiple prognostic FRGs (Additional file 5: Figure S5H, I).

Fig. 7figure 7

FRGs regulate the level of immunocytes infiltration and immune checkpoint expression in LUAD patients. A, B The change trend of the infiltration level of Mast-cells and T cells CD4 memory activated with the expression of RRM2. C Different Check-point levels represent different LUAD OS. D 27 common checkpoint genes are differentially expressed between high and low risk groups. E, F The expression of ROS1 and ALK change with risk score. G, H Correlation of the expression of ROS1 and ALK with 7 key FRGs

Overexpression of RRM2 in tumor tissues

To detect the expression of RRM2 in normal lung tissues and lung cancer tissues, we performed immunohistochemical experiments for comparison. The results showed that the expression of RRM2 in tumor tissues was higher than that in normal lung tissues (Fig. 8).

Fig. 8figure 8

The expression of RRM2 in normal lung tissue and lung cancer tissue was detected by immunohistochemistry. 674, 677, 712, 728, 733 and 734 are patient numbers. Expression of RRM2 in LUAD tissues was higher than that in normal tissues

Silencing RRM2 induced ferroptosis in lung cancer cells

Ferroptosis is a recently discovered form of cell death characterized by lipid peroxidation accumulation and iron dependence. Further study on the specific mechanism of ferroptosis is expected to bring new prospects for cancer treatment. To determine the clinical relevance of RRM2 expression, we detected RRM2 expression in a normal lung cell line (BEAS-2B) and seven lung cancer cell lines (H1299, A549, H460, H23, H838, PC-9, and H1975). Western blotting showed that RRM2 protein was highly expressed in H1975 cells. In addition, we detected the expression of 4HNE and ACSL4 in different lung cancer cell lines. The results showed that 4HNE and ACSL4 were highly expressed in H1975 cells (Fig. 9A), indicating that H1975 cells may be more sensitive to ferroptosis. Therefore, we used H1975 cells for subsequent experiments. After H1975 cells were transfected with RRM2siRNA for 48 h, cell death was detected by flow cytometry. The results showed that silencing RRM2 could induce cell death (Fig. 9B). To further study the effect of silencing RRM2 on ferroptosis, Western blotting was used again to detect the expression of several proteins in H1975 cells with RRM2 siRNA transfected. The results showed that the expression of RRM2 in the silencing group was significantly decreased, indicating that silencing RRM2 was successfully. In the meantime, the expression of ferroptosis-related indicators like SLC7A11 and GPX4 were decreased, while ACSL4 were increased in H1975 cells after silencing RRM2 (Fig. 9C). Then, we continued to detect the levels of ferrous ion and lipid peroxidation in lung cancer cells after silencing RRM2. The results showed that silencing RRM2 could induce an increase in ferrous ion level (Fig. 9D) and lipid peroxidation accumulation (Fig. 9E) in H1975 cells. Ultimately, we demonstrated that silencing RRM2 induced ferroptosis in H1975 cells.

Fig. 9figure 9

Effects of silencing RRM2 on ferroptosis. A Expression level of RRM2,4HNE and ACSL4 in different lung cancer cell lines by Western blotting. B The results of cell death level after silencing RRM2 by flow cytometry. C The expression of ferroptosis-related indicators was changed in H1975 cells after silencing RRM2. D The fluorescence intensity of cells after silencing RRM2 was evaluated using a rotating disk super-resolution laser confocal microscope. E Detection of cell lipid peroxidation level after silencing RRM2 by flow cytometry

Screening drugs for high-risk FRGs

To predict effective therapeutic drugs for patients with high-risk LUAD, we constructed a ridge regression model to predict the drug IC50using the GDSC cell line expression profile and TCGA gene expression profile; we screened six antitumor drugs with significantly lower IC50 in high-risk LUAD: bosutinib, dasatinib, gefitinib, tipifarnib, docetaxel, and JNK inhibitor VIII (Additional file 6: Figure S6A–F). These six antitumor drugs could kill tumor cells in patients with high-risk LUAD. In addition, we compared the risk genes with the reference data set of the CMap database and obtained the correlation of 1309 drugs according to the enrichment of risk genes in the reference gene expression profile. With P < 0.05 as the standard, 158 drugs were found to be significantly enriched in high-risk genes; of these, 86 drugs promoted the expression of high-risk genes, and 72 drugs inhibited the expression of high-risk genes. Medrysone, phenoxybenzamine, vorinostat, thioguanosine, apigenin, and chrysin were selected as the six drug candidates enriched with more prominent inhibitors of high-risk gene expression, and their chemical structures were found in the PubChem database (Additional file 6: Figure S6G-L). The detailed information on these six drugs is shown in Table 2.

Table 2 Six small molecule drugs for high-risk LUAD patients slected from CMAP database

留言 (0)

沒有登入
gif