Screening and biomarker assessment of ferroptosis genes FLT3 and ALOX5 in lung adenocarcinoma

The Cancer Genome Atlas and FerrDb analysis

We retrieved transcriptome data from The Cancer Genome Atlas (TCGA) database, consisting of 522 patients with LUAD and 59 samples of normal lung tissue. FerrDb (http://www.zhounan.org/ferrdb) is an openly accessible database that provides information on the regulators, markers, and links to diseases related to ferroptosis. In this study, 108 drivers were selected for further study. The edgeR package from Bioconductor in R was employed to identify and analyze the differentially expressed genes (DEGs) between normal lung tissues and LUAD samples [14]. A significance threshold of p-value<0.05 was used for identifying differentially DEGs.

Kaplan–Meier survival plot

The Kaplan-Meier (K-M) Plotter (http://www.kmplot.com/analysis/) linked with PAN-CANCER software was employed to generate survival curves for LUAD patients based on their expression levels of ferroptosis genes [15]. The gene expression levels were categorized into two groups (high and low) using the median value, and subsequent survival analyses were conducted. Hazard ratio (HR) and logarithmic rank p values were calculated for 95% confidence intervals.

Gene correlation analysis

To further confirm the correlation between ferroptosis genes, we utilized the corrplot package in R to analyze the TCGA expression data and assess the expression level correlation coefficients using the Spearman method. On the x-axis of the analysis, the selected ferroptosis-related genes with significant differences in expression were represented, while on the y-axis, homologous genes were represented. Tumor and para-cancer tissues from LUAD patients were analyzed. To identify the functional classes of selected genes, STRING (https://string-db.org/) was utilized to create networks of protein-protein interactions (PPIs) via Cytoscape to reveal functional relationships between iron failure genes and the top 10 most relevant genes.

Functional analysis of selected ferroptosis genes

To examine the functions of the selected genes, pathway enrichment analysis was performed by submitting the selected genes to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, systematically analyzed the metabolic pathways and gene product functions in cells [16, 17]. Besides, gene ontology (GO) term enrichment analysis was perform by GOrilla online software [18]. The annotation, visualization, and synthetic discovery databases (DAVID) (https://david.ncifcrf.gov/) were also employed to conduct GO and KEGG pathway enrichment analysis. In line with the DAVID manual, 23 ferroptosis genes were uploaded, and their functions were subsequently mapped. The important molecular functions and pathways of the target genes were extracted. We used histograms and Bubble Diagram Drawing Software to represent the GO and KEGG results of the ferroptosis genes.

TIMER database analysis

To investigate the relationship between the expression level of 23 ferroptosis genes and the tumor immune-microenvironment in LUAD, the tumor immune estimation resource database (TIMER, https://cistrome.shinyapps.io/timer/) was applied to estimate the abundance of different types of immune cells that infiltrate tumors based on TCGA data via deconvolution statistics [19]. In addition, the “Correlation” module was utilized to explore the relationship between ferroptosis-related genes expression and a specific set of immune markers.

Patient selection and evaluation

The patients included in this study had primary LC treated at the Zhejiang Cancer Hospital (Zhejiang, China) from January 2014 to June 2021. The diagnosis of LC was established based on clinical and X-ray examinations and confirmed through histological analysis of tumor biopsies. The staging and TNM classification of the LC cases were determined according to the 8th edition of the American Joint Committee on Cancer (AJCC) guidelines. The inclusion criteria included previous surgery and histological confirmation of LUAD. Participants who had multiple organ dysfunction or a history of organ transplantation, lung squamous cell carcinoma (LUSC), or small cell lung cancer (SCLC) were excluded from the study, as well as those who were unwilling to participate. All participants in the study provided samples of their tumor tissues that had been embedded in paraffin, as well as paired paracancerous specimens, which were obtained during surgical operations. In addition to the tissue specimens, information about the participants’ clinical characteristics was collected in a careful manner.

IHC and H&E staining

Briefly, the paraffin-embedded tumor tissue specimens and paired paracancerous specimens were sectioned into 4-μm sections and subjected to hematoxylin and eosin (H&E) staining using the standard protocol. Additionally, the ALOX5 gene expression level of these samples was quantified by immunohistochemistry (IHC). IHC assays were performed using an automated analytical instrument (BenchMark ULTRA system, Roche, Basel, Switzerland) and stained with 5-LO Polyclonal Antibody (YT0027, Immunoway, Texas, USA) according to the manufacturers’ protocols. Microscopic examination of the tumor cells and IHC analysis was carried out using an optical microscope (BX43, Olympus, Tokyo, Japan). Two experienced pathologists independently assessed and scored the percentage of tumor cells and IHC results retrospectively.

Statistical analysis

Multivariate Cox regression and K-M analysis were conducted using IBM SPSS Statistics 20.0 (IBM, New York, USA). Additional statistical analysis was carried out using R software (version 3.5.2). The correlation between gene expression was evaluated using Spearman correlation, with the strength of correlation judged based on the following absolute values: 0.00–0.19, “very weak”; 0.20–0.39, “weak”; 0.40–0.59, “moderate”; 0.60–0.79, “strong”; 0.80–1.0, “very strong”. Inter-group differences in IHC quantitative results were analyzed using paired samples t-tests, with statistical significance determined by p<0.05. This same significance level was used as the inclusion standard in the analysis.

Expression levels of ferroptosis genes in LUAD

We obtained 108 driver ferroptosis genes through the FerrDb database. The mRNA expression levels of iron ferroptosis genes were analyzed in TCGA using the edgeR package from Bioconductor in order to uncover disparities in their expression between LUAD and normal tissues. Building upon the results, 89 ferroptosis genes exhibited differential expression (p<0.05) between LUAD and normal tissues (Figure 1). Of these, 56 ferroptosis genes were up-regulated, among which 5 genes (MIOX, CDKN2A, TFR2, NOX5 and ALOXE3) had a log2 fold change (FC) greater than 2. Thirty-three ferroptosis genes were down-regulated, among which 3 genes (DUOX1, CDO1 and EPAS1) had a log2 FC less than −2 (Table S1).

Figure 1:

Human ferroptosis gene expression levels in lung adenocarcinoma (LUAD) from TCGA. (A) Up-regulated genes. (B) Down-regulated genes. In the heatmap, the annotation “t” represents “tumor” and “n” represents “normal”, *p<0.05, **p <0.01, ***p<0.001.

Prognostic significance of ferroptosis genes in LUAD

We investigated and explored the potential association between ferroptosis gene expression and OS in LUAD. First, the effect of ferroptosis gene expression on survival rate was evaluated using TCGA. The findings indicated a significant correlation between the expression of ferroptosis genes and the prognosis of LUAD. Further evaluation of the impact of ferroptosis gene expression on the survival probability of LUAD patients was conducted using the K-M Plotter. Twenty-three ferroptosis genes were selected according to a P<0.05 and consistent prognostic trends based on data analysis from RNA-Seq and gene chips. The results of the data analysis based on RNA-Seq (Figure 2A) and RNA microarray (Figure 2B) showed that high expression of ferroptosis genes (CS, NOX5, G6PD, PGD, VDAC2, GOT1, HMOX1, ULK1, CDKN2A and HILPDA) was notably correlated with poor prognosis. On the contrary, patients with LUAD who had high expression levels of various other ferroptosis genes (FLT3, ALOX5, NCOA4, GABARAPL2, WIPL1, ZEB1, PEBP1, CDO1, EPAS1, TAZ, GLS2, ALOX15B, and DPP4) was markedly longer (Figure 2A, B).

Figure 2:

Forest plot of the prognostic survival analysis of 23 ferroptosis genes. (A) Forest plot of the prognostic survival analysis of 23 ferroptosis genes based on RNA-Seq data analysis. (B) Forest plot of the prognostic survival analysis of 23 ferroptosis genes based on gene chip data.

Correlation analysis between ferroptosis genes and constructing protein interaction networks

We use the corrplot package in R to further verify the correlation between ferroptosis genes (Figure 3A) and observed that ALOX5 had a notable positive correlation with CDO1 (0.42) and EPSA1 (0.43). CS was positively correlated with VDAC2 (0.41). ZEB1 was positively correlated with EPAS1 (0.49), CDO1 (0.41) and FLT3 (0.44). Additionally, we observed a relatively stronger positive correlation between G6PD and PGD (0.64), as well as between CDO1 and EPAS1 (0.63). Finally, TAZ was negatively correlated with NCOA4 (−0.47).

Figure 3:

Correlation analysis between ferroptosis genes and their interacting proteins in LUAD, and functional analysis of 23 ferroptosis genes. (A) Correlation analysis between the ferroptosis genes. The darker the color, the higher the correlation coefficient. Red and blue represent positive and negative correlation, respectively. (B) PPI network of ferroptosis genes including the top 10 proteins. (C) GO biological process analysis of 23 ferroptosis genes in LUAD. (D) KEGG pathway analysis of 23 ferroptosis genes in LUAD. The circle’s size and color show the number and −log10 (p-value), respectively. The x-axis represents the average score of genes in the specific KEGG terms. KEGG, Kyoto Encyclopedia of genes and genomes.

Functional interactions between proteins are essential for molecular mechanisms and the metabolism of malignant tumors. To confirm the interaction of ferroptosis proteins in LUAD progression, the STRING tool was utilized for the analysis of the PPI network. The top 10 proteins, along with their corresponding gene names and annotations, were listed in Table 1. The high degree of connectivity between these protein-protein pairs suggests that these ferroptosis proteins play critical roles during LUAD progression (Figure 3B).

Functional annotation of ferroptosis genes

By GO analysis, ferroptosis genes were grouped into biological processes, cell components, and molecular functions. The biological processes included pentose-phosphate shunting and pentose biosynthetic process. The cellular components primarily involved were pre-autophagosomal structure membrane and the molecular function analysis indicated a major involvement of oxidoreductase activity (Figure 3C).

By KEGG pathway analysis, the ferroptosis genes were mapped to 10 high-abundance KEGG pathways, among which four major pathways were cancer central carbon metabolism, cancer microRNAs, carbon metabolism, and antibiotic biosynthesis (Figure 3D). In addition, Tables S2 and S3 list the details of 23 ferroptosis genes in the three GO classes and pathways.

Correlation of ferroptosis gene expression with immune marker sets and immune infiltration levels

We further examine the correlation between the expression of ferroptosis genes and immune marker sets in LUAD. FLT3 expression was positively correlated with the expression of CD3E, CCR7, CD2, CD19, CD3D, CD79a, CD86 and CSF1R. ALOX5 was positively correlated with the expression of CD86, CSF1R, CD68, ITGAM, VISIG4, IRF5, CD163, MS4A4A, HLA-DPB1, HLA-DPA1, ITGAX and STAT5A (Figure 4A). These findings indicate that elevated levels of FLT3 and ALOX5 expression may have a significant impact on immune evasion. Moreover, our results suggests that high expressed FLT3 and ALOX5 is broadly linked to immunity in patients with LUAD. Besides, The relationship between the expression of FLT3 and ALOX5 and the infiltration levels of six main types of immune cells, which contains B cells, CD4+/ CD8+ T cells, dendritic cells (DCs), macrophages, and neutrophils was analyzed. The results displayed that FLT3 expression was positively associated with infiltration levels of B cells and moderately correlated with infiltration levels of CD4+ T cells, neutrophils, and DCs in LUAD (Figure 4B). Besides, ALOX5 expression demonstrated a positive correlation with infiltration levels of macrophages, neutrophils and CD4+ T cells, with the strongest correlation with the infiltration level of DCs in LUAD (Figure 4C).

Figure 4:

Correlation of ferroptosis gene expression with immune marker sets and immune infiltration levels and IHC results. (A) Heat maps show the correlations between ferroptosis genes and immune marker sets correlations. The darker the color, the higher the correlation coefficient. Red and blue represent positive and negative correlation, respectively. (B) FLT3 was positively correlated with infiltrating levels of B cells and was moderately correlated with infiltration levels of CD4+ T cells, neutrophils and DCs in LUAD. (C) ALOX5 expression was positively correlated with the infiltration level of macrophages, neutrophils and CD4+ T cells, with the strongest correlation with the infiltration level of DCs in LUAD. p<0.05 was considered statistically significant. (D) Box-plot of IHC results in tumor tissue and matched paracancerous tissue, p<0.01. (E) H&E and IHC staining results.

Patient characteristics

From 2020 to 2021, we recruited 43 newly diagnosed LC patients to investigate any potential differences in ALOX5 expression between tumor tissue and paracancerous tissue. The study participants had a median age of 61 years, and their tumor size and TNM stage were carefully assessed and recorded (Table S4). Table 2 provides a summary of their characteristics. For each patient, we obtained paraffin-embedded tumor tissue specimens and paired paracancerous specimens through surgical procedures, which were then subjected to H&E and IHC staining.

IHC assay of ALOX5 protein expression

IHC assay was performed on the 43 patient-derived samples to quantify the ALOX5 protein expression level. Based on the positive signal intensity, samples were divided into ++, +∼++, and + levels, corresponding to a score of 1, 2 and 3, respectively. Additionally, the samples were scored 0, 1, 2, 3 and 4, corresponding to the percentage of positive cells for ALOX5 as follows: 0%, ≤25%, 26%–50%, 51%–75% and >75%, respectively. The final score of the ALOX5 expression level is the multiplication of these two scores (Table S4). Our results showed that compared with normal tissue, ALOX5 expression significantly decreased in tumor tissue (p <0.01), indicating that it is a potential biomarker in LC diagnostic (Figure 4D, E).