Identification of potential key molecular biomarkers in lung adenocarcinoma by bioinformatics analysis.

Background: Lung cancer is one of the most common malignant tumors in the world, of which the rate of incidence has continuously increased over recent years. Lung adenocarcinoma (LUAD) is the most frequent pathological type of lung cancer.
Methods: In order to discover the key markers for the occurrence and development of LUAD, we collected messenger RNA (mRNA) expression datasets in the Gene Expression Omnibus (GEO), namely, GSE2514, GSE7670, and GSE40275. The differentially expressed genes (DEGs) were screened using the online interface between GEO and R (GEO2R). Then, DEGs were functionally annotated in the Database for Annotation, Visualization, and Integrated Discovery (DAVID). Next, a protein-protein interaction (PPI) network was drawn by using the Search Tool for the Retrieval of Interacting Genes (STRING) web tool and Cytoscape software. Finally, Kaplan-Meier plotter was utilized to analyze the overall survival (OS) of the hub genes. The correlation between fibroblast growth factor 2 (FGF2) and immune infiltration was studied by TIMER web services.
Results: In this study, we obtained a total of 284 DEGs through the intersection of 3 datasets, and found that DEGs were highly related to biological processes such as "cell adhesion", "cell differentiation", and "cell proliferation". After that, the hub genes were obtained by analyzing the PPI network. Finally, we found that the abnormal expression of hub genes is obviously related to poor prognosis in LUAD patients. The expression level of FGF2 was positively correlated with the immune infiltration in LUAD. Conclusions: In general, the DEGs and hub genes can provide new research targets for the development of LUAD, as well as potential diagnosis and treatment strategies for disease treatment. In particular, FGF2 expression was found to be involved in the immune microenvironment of LUAD. 2022 Translational Cancer Research. All rights reserved.

Introduction

Lung cancer is one of the most common malignant tumors in the world, with recent statistics indicating that about 26% of tumor patients die from this disease (1). Moreover, lung adenocarcinoma (LUAD) is one of the most common pathological types of non-small cell lung cancer (NSCLC). In China, the incidence of lung cancer continues to follow an annually increasing trend (2), and LUAD has become the most common type. Therefore, it is necessary to explore the mechanism of occurrence and development of LUAD, particularly the abnormally expressed genes that play an important role. By discovering genes and studying their corresponding functions, it may be possible to provide new strategies for the treatment of LUAD.

Microarray technology and bioinformatics have been extensively used in tumor research. Many current studies are performed using various analytical tools to compare microarray datasets to obtain DEGs in various tumors. Then, the role of DEGs is explored in molecular functions (MFs), biological processes (BPs), and different pathways to provide new ideas for tumor research (3,4). However, due to various factors such as sample size, tumor stage, grade, and ethnicity, a complete set of DEGs is often not available. Therefore, it is important for us to compare the latest microarray datasets repeatedly in order to obtain more representative DEGs. Herein, 3 LUAD messenger RNA (mRNA) datasets were screened from the Gene Expression Omnibus (GEO) database, and a total of 284 DEGs were discovered by intersection. Subsequently, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were employed for the functional annotation of DEGs, and a protein-protein interaction (PPI) network was established through Cytoscape (https://cytoscape.org/). In summary, we used the above methods to illustrate the molecular mechanism of LUAD. As a result, 284 DEGs and 10 hub genes were screened, which may be used as potential biomarkers for LUAD. In this study, we found low expression of fibroblast growth factor 2 (FGF2) in LUAD and it is correlated with the long-term prognosis among the patients. Specifically, with higher FGF2 expression, patients tend to have a longer life span after diagnosis, which is not consistent with previous studies (5,6). Thus, more effort is needed to clarify the mechanism between FGF2 and LUAD. Herein, we explore the role of FGF2 in the immune infiltration of LUAD, which may be another prognostic factor that provides a novel prospect for the LUAD’s treatment. We present the following article in accordance with the REMARK reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-21-2676/rc).

Methods

Source of data

The datasets were collected from the open GEO (https://www.ncbi.nlm.nih.gov/geo/). We selected 3 lung cancer datasets, GSE2514 (7), GSE7670 (8), and GSE40275 (9), and human LUAD was the main pathological type. Among them, GSE2514 contains 19 tumor samples and 19 normal lung tissue samples, GSE7670 contains 57 tumor samples and 57 normal lung tissue samples, and GSE40275 contains 8 tumor samples and 43 normal lung tissue samples. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

DEGs screened

The GEO2R online tool (https://www.ncbi.nlm.nih.gov/geo/geo2r/) was used to screen tumor tissues and normal human lung tissues, and statistically analyze each dataset to calculate the adjusted P value and |logFC| and define the required DEGs that meet the conditions, that is, adjusted P<0.01, |logFC| ≥1.0. Finally, the intersection was acquired from the Venn diagram online drawing tool (http://bioinformatics.psb.ugent.be/webtools/Venn/).

GO and KEGG pathway analysis of DEGs

The GO analysis is a shared method for large-scale functional enrichment research, in which gene functions can be divided into cellular component (CC), MF, and BP. The KEGG is a common biological information database. The Database for Annotation, Visualization, and Integrated Discovery (DAVID) tool was conducted to perform GO and KEGG pathway analysis on DEGs (https://david.ncifcrf.gov/). Statistical significance was considered when P<0.01 and the gene count was ≥10.

Construction of PPI network and selection of hub genes

The Search Tool for the Retrieval of Interacting Genes (STRING), an open online tool (http://string-db.org/), was employed to analyze PPI information. A combined score was >0.4. Then, the PPI network was constructed by Cytoscape software. Finally, the degree of each protein node was calculated using the cytohubba plugin of Cytoscape. The top 10 genes were defined as hub genes of the network in our study. Clustering of the top 10 genes was constructed using an open multi-omics and clinical/phenotypic data exploration tool, UCSC Xena (http://xena.ucsc.edu/) (10).

Survival analysis of hub genes

The mRNA lung cancer database was used to evaluate the value of hub genes on survival and prognosis of LUAD patients by Kaplan-Meier plotter (http://kmplot.com/analysis/). The Kaplan-Meier plotter is an open online tool that mainly analyzes the impact of 54,675 genes on the survival and prognosis of patients in most tumor types. The value of hub genes for survival and prognosis in LUAD patients can be assessed using the mRNA lung cancer database in Kaplan-Meier plotter. The gene probe was based on “only jetset best probe set”. According to the median of mRNA expression, tumor patients were divided into high expression group and low expression group. Statistical significance was considered when P<0.01. The immunohistochemical data of 9 hub genes in LUAD or normal lung tissue were obtained from the Human Protein Atlas (HPA) (https://www.proteinatlas.org/).

PrognoScan database analysis

PrognoScan database is an online database (http://dna00.bio.kyutech.ac.jp/PrognoScan/), which uses a large number of cancer microarray dataset for the analysis of gene expression and relations between survival prognosis of patients. FGF2 expression and survival in patients with LUAD were analyzed in the PrognoScan database. Cox P<0.05 was statistically significant.

TIMER database analysis

TIMER is a comprehensive online web tool for systematic analysis of immune infiltration in different cancer types (https://cistrome.shinyapps.io/timer/). The DiffExp tool in TIMER was used to study the difference of FGF2 expression between tumor and normal tissues in The Cancer Genome Atlas (TCGA) tumor database. The gene module visualized the correlation between FGF2 expression and immune infiltration level. The results included the correlation analysis with B cells, CD4+ T cells, CD8+ T cells, Neutrophils, macrophages, and dendritic cells. The scatter plots were used to show the correlation between FGF2 and different immune cell infiltration. The gene expression level was determined by log2 RSEM.

Statistical analysis

The GEO2R is used to screen DEGs. The functional enrichment research of DEGs uses GO analysis. The STRING and Cytoscape were used to construct PPI network. The Kaplan-Meier Plotter database and Prognoscan database were used to analyze the survival of hub gene. The correlation between gene expression and immune infiltration was displayed by Spearman correlation and statistical significance. The correlation of variables is established by the following absolute values: 0.00–0.19 “very weak”, 0.20–0.39 “weak”, 0.40–0.59 “medium”. P<0.05 was statistically significant.

Results

The DEGs were retained in LUAD. Through the analysis of GSE2514, GSE7670, and GSE40275, we obtained 840, 1,410, and 5,454 DEGs, respectively. Afterwards, the intersection of the DEGs of each dataset was taken by the Venn diagram, and 284 DEGs were obtained, of which 77 were highly expressed in tumor tissues and 207 were lowly expressed ().

Figure 1

Venn diagram was drawn to obtain DEGs. (A) Up-regulated genes; (B) down-regulated genes. DEGs, differentially expressed genes.

Functional enrichment of DEGs was carried out via GO and KEGG analysis. The GO analysis revealed that DEGs were mainly enriched in BPs such as cell adhesion, cell division, extracellular matrix organization, positive regulation of gene expression, positive regulation of angiogenesis, response to hypoxia, and so on. CC was mainly enriched in extracellular exosome, extracellular region, plasma membrane, membrane, integral component of plasma membrane, extracellular space, and so on. MF was primarily involved protein binding, protein kinase binding, protein homodimerization activity, calcium ion binding, and heparin binding execution. The details are shown in and . The KEGG analysis was mainly related to the cell cycle ().

Table 1

GO and KEGG pathway analysis of DEGs in LUAD samples

Category	Term	Description	Count	P value
BP Term	GO:0007155	Cell adhesion	33	7.79e−12
BP Term	GO:0051301	Cell division	21	1.83e−06
BP Term	GO:0030198	Extracellular matrix organization	17	1.92e−07
BP Term	GO:0010628	Positive regulation of gene expression	17	8.76e−06
BP Term	GO:0045766	Positive regulation of angiogenesis	16	7.94e−10
BP Term	GO:0001666	Response to hypoxia	16	1.97e−07
BP Term	GO:0007179	Transforming growth factor beta receptor signaling pathway	12	3.70e−07
BP Term	GO:0030336	Negative regulation of cell migration	11	4.11e−06
BP Term	GO:0051726	Regulation of cell cycle	11	4.34e−05
CC Term	GO:0005886	Plasma membrane	92	9.65e−05
CC Term	GO:0070062	Extracellular exosome	76	7.16e−07
CC Term	GO:0005576	Extracellular region	59	7.18e−10
CC Term	GO:0016020	Membrane	58	5.25e−05
CC Term	GO:0005887	Integral component of plasma membrane	53	3.36e−09
CC Term	GO:0005615	Extracellular space	52	1.78e−09
CC Term	GO:0009986	Cell surface	30	6.20e−09
CC Term	GO:0005578	Proteinaceous extracellular matrix	21	6.95e−09
CC Term	GO:0031012	Extracellular matrix	19	8.52e−07
CC Term	GO:0045121	Membrane raft	18	2.21e−08
CC Term	GO:0016324	Apical plasma membrane	15	1.80e−04
CC Term	GO:0009897	External side of plasma membrane	12	4.88e−04
CC Term	GO:0043235	Receptor complex	11	2.77e−05
CC Term	GO:0016323	Basolateral plasma membrane	11	4.96e−04
CC Term	GO:0000776	Kinetochore	10	4.29e−06
CC Term	GO:0030496	Midbody	10	1.77e−04
MF Term	GO:0005515	Protein binding	181	9.00e−08
MF Term	GO:0019901	Protein kinase binding	22	5.44e−07
MF Term	GO:0042803	Protein homodimerization activity	30	4.36e−06
MF Term	GO:0005509	Calcium ion binding	28	2.46e−05
MF Term	GO:0008201	Heparin binding	12	4.57e−05
KEGG pathway	hsa04110	Cell cycle	11	6.77e−04

GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; LUAD, lung adenocarcinoma; BP, biological process; CC, cellular component; MF, molecular function.

Figure 2

GO enrichment analysis result of DEGs. (A) BP; (B) CC; (C) MF. GO, Gene Ontology; DEGs, differentially expressed genes; BP, biological process; CC, cellular component; MF, molecular function.

The PPI network was constructed and hub genes were extracted. The PPI network consisted of 262 nodes and 1,363 interaction lines, in which the red nodes were highly expressed genes, and the green nodes were lowly expressed genes ().

Figure 3

Establish a PPI network and extract the top 10 hub genes. (A) The PPI network was established by Cytoscape; (B) the hub genes were extracted by using cytohubba; (C) top 10 hub gene obtained by screening. PPI, protein-protein interaction; EZH2, enhancer of zeste 2 polycomb repressive complex 2 subunit; AURKA, aurora kinase; CCNA2, cyclin A2; CCNB1, cyclin B1; FGF2, fibroblast growth factor 2; FOXM1, forkhead box M1; MKI67, marker of proliferation Ki-67; TYMS, thymidylate synthetase; TOP2A, topoisomerase II alpha; CENPE, centromere protein E.

Next, the PPI networks with connectivity degree were further screened out. The top 10 genes were obtained, namely enhancer of zeste 2 polycomb repressive complex 2 subunit (EZH2), aurora kinase A (AURKA), cyclin A2 (CCNA2), cyclin B1 (CCNB1), FGF2, forkhead box M1 (FOXM1), marker of proliferation Ki-67 (MKI67), thymidylate synthetase (TYMS), topoisomerase II alpha (TOP2A), and centromere protein E (CENPE) (), and a network diagram of their interaction was drawn (). The results of hierarchical clustering showed that the hub genes could basically distinguish LUAD from normal lung tissue ().

Table 2

Top 10 hub genes with higher degree of connectivity

Gene symbol	Degree	Gene description
EZH2	43	Enhancer of zeste 2 polycomb repressive complex 2 subunit
AURKA	41	Aurora kinase A
CCNA2	40	Cyclin A2
CCNB1	40	Cyclin B1
FGF2	39	Fibroblast growth factor 2
FOXM1	39	Forkhead box M1
MKI67	39	Marker of proliferation Ki-67
TYMS	38	Thymidylate synthetase
TOP2A	38	Topoisomerase (DNA) II alpha
CENPE	37	Centromere protein E

Figure 4

Constructed clustering of hub genes via the UCSC Xena platform. In these tissues, red denotes up-regulated genes, and blue denotes down-regulated genes. EZH2, enhancer of zeste 2 polycomb repressive complex 2 subunit; AURKA, aurora kinase; CCNA2, cyclin A2; CCNB1, cyclin B1; FGF2, fibroblast growth factor 2; FOXM1, forkhead box M1; MKI67, marker of proliferation Ki-67; TYMS, thymidylate synthetase; TOP2A, topoisomerase II alpha; CENPE, centromere protein E.

Survival analyses of hub genes were performed. According to the hub genes expression, the survival curves of LUAD patients were drawn via Kaplan-Meier plotter (). It was found that patients with high expression of EZH2, AURKA, CCNA2, CCNB1, FOXM1, MKI67, TYMS, TOP2A, and GEPNE have a poor prognosis, while patients with high FGF2 gene expression have a better prognosis, which is consistent with the previous expression results in the mRNA datasets. The expression of hub genes in LUAD and normal lung tissues was obtained from immunohistochemical data of HPA (). The expression results of the hub genes in different sample types of immunohistochemistry were consistent with the results of the above database.

Figure 5

Survival analysis of hub gene. P<0.05. HR, hazard ratio. EZH2, enhancer of zeste 2 polycomb repressive complex 2 subunit; AURKA, aurora kinase; CCNA2, cyclin A2; CCNB1, cyclin B1; FGF2, fibroblast growth factor 2; FOXM1, forkhead box M1; MKI67, marker of proliferation Ki-67; TYMS, thymidylate synthetase; TOP2A, topoisomerase II alpha; CENPE, centromere protein E.

Figure 6

Immunohistochemical images of lung adenocarcinoma and normal lung tissue from the Human Protein Atlas database. Magnification: ×100. EZH2, enhancer of zeste 2 polycomb repressive complex 2 subunit; AURKA, aurora kinase A; CCNA2, cyclin A2; CCNB1, cyclin B1; FGF2, fibroblast growth factor 2; FOXM1, forkhead box M1; MKI67, marker of proliferation Ki-67; TYMS, thymidylate synthetase; TOP2A, topoisomerase II alpha.

FGF2 and tumor immune microenvironment. We analyzed the expression of FGF2 in TCGA database by TIMER web service. It also found that FGF2 was highly expressed in normal tissues compared with tumor tissues (). We further confirmed that the high expression of FGF2 was associated with better prognosis in LUAD patients in PrognoScan database (). Then, TIMER was used to analyze the relationship between FGF2 and tumor purity in LUAD, and it was found that the expression level of FGF2 was negatively correlated with tumor purity (). The expression level of FGF2 was found to be positively correlated with the immune infiltration of B cells, CD4+ T cells, CD8+ T cells, Neutrophils, macrophages, and dendritic cells. There was moderate correlation with neutrophils infiltration, but weak correlation with CD4+ T cells, CD8+ T cells, macrophages, and dendritic cells infiltration.

Figure 7

The expression of FGF2 in LUAD. (A) The expression of FGF2 in lung adenocarcinoma was analyzed from TCGA database by TIMER. (B) Kaplan-Meier survival curves comparing the FGF2 differential expression in lung adenocarcinoma in the PrognoScan database. *, P<0.05; **, P<0.01; ***, P<0.001. FGF2, fibroblast growth factor 2; LUAD, lung adenocarcinoma; TCGA, The Cancer Genome Atlas; TPM, transcripts per million; HR, hazard ratio; CI, confidence interval.

Figure 8

Correlation between FGF2 expression and immune infiltration in lung adenocarcinoma. FGF2, fibroblast growth factor 2; TPM, transcripts per million.

Discussion

The most common type of lung cancer is NSCLC, accounting for 85% of all lung cancer cases. As the most common subtype of NSCLC, LUAD comprises approximately 40% of the cases (11), and the incidence of LUAD is increasing annually. Although there has been significant improvement in the diagnosis and treatment strategies for LUAD, such as targeted treatments, the long-term survival rate of patients has remained low. The main reason is the dynamic drug resistance of the tumors (12). Therefore, it is urgent that the pathogenesis of LUAD be further studied in order to find new diagnostic and therapeutic targets. With the development of bioinformatics, the screening of DEGs provides a valuable tool for research. We may find new gene targets by comparing and extracting a large-scale LUAD mRNA expression dataset.

In this study, we screened 284 DEGs by comparing 3 datasets, including 77 high expression genes and 207 low expression genes. After GO analysis of 284 DEGs, it was found that their biological functions were mainly concerned cell adhesion and division, while KEGG pathway was enriched in the cell cycle. We further extracted the top 10 hub genes by PPI network, which were EZH2, AURKA, CCNA2, CCNB1, FGF2, FOXM1, MKI67, TYMS, TOP2A, and CENPE. Among them, FGF2 was a low-expression gene, and the others were high expression genes. Finally, survival analysis showed that the expression results of hub genes were consistent in patients, and the prognosis was strongly correlated with the hub gene expression.

The EZH2 gene is a component of polycomb repressive complex 2 (PRC2), a methyltransferase. Increasing DNA methylation and inactivating tumor suppressor genes are important factors by which EZH2 promotes tumor progression (13,14). It can also enhance cell proliferation by promoting cell cycle processes (15). A study found that EZH2 acetylation enhanced the metastasis and invasion ability of tumor cells and that it was correlated with poor prognoses in LUAD patients (16). It has been reported that the vascular endothelial growth factor (VEGF)/vascular endothelial growth factor receptor 2 (VEGFR-2) pathway can regulate the EZH2 expression in LUAD (17).

The AURKA gene is a serine/threonine kinase that plays a key role in regulating the cell cycle and mitosis of normal cells (18), and its abnormal expression has been also correlated with various types of tumors (19-21). Some studies have found that the activation of AURKA in NSCLC can increase the resistance to anti-EGFR targeted therapies (22,23). In addition, the differential expression of AURKA in LUAD has been confirmed by previous studies (24,25). In our study, we also found that AURKA was highly expressed in adenocarcinoma, and the high expression of AURKA in LUAD patients was obviously correlated with poor prognosis. Therefore, the level of AURKA expression may be used as a diagnosis tool and treatment target for LUAD. The expression of AURKA can be further studied for its role in the development and progression of LUAD.

The FGF2 gene is a member of the fibroblast growth factor family. It mainly binds FGFR on the cell surface and is involved in cell proliferation, metastasis, and angiogenesis. In previous studies, it was found that FGF2 expression in NSCLC was negatively correlated with the survival rate (26,27). However, in our study, we found that FGF2 was lowly expressed in LUAD tissues, and survival analysis also showed that patients with low FGF2 expression had better prognoses. This shows that the mechanism of FGF2 in LUAD is not fully clarified, and further study should be undertaken. Therefore, we further investigated the mechanism of FGF2 affects the prognosis of LUAD patients. It is known that the immune microenvironment plays an important role in NSCLC, and immunotherapy (PD-L1 inhibitors) has a significant effect in the treatment of lung tumors (28,29). So, we analyzed the effect of FGF2 expression on the immune microenvironment in LUAD. It is noteworthy that this study found that the high FGF2 expression decreased tumor purity, and the FGF2 expression level was positively correlated with immune cell infiltration. Although the correlation was not very strong, the significance was very high (). Thus, we can appropriately believe that immune infiltration is one of the reasons for the better prognosis of patients with high FGF2 expression.

The FOXM1 gene is a transcription factor of the forkhead box family. Its high expression is associated with various types of tumors. It plays an important role in cell proliferation, metastasis, drug resistance, and the promotion of vascular survival and cell cycling (30-32). In addition, studies have found that FOXM1 is correlated with the development of LUAD, being specifically associated with its proliferation, invasion, and metastasis (33-35).

The CENPE gene is a kinesin-like motor protein, and its overexpression is related to the development and progression of tumors (36). In NSCLC, the high expression of CENPE was correlated with the negative prognoses of patients (37). In addition, studies have found that CENPE can directly act on FOXM1 to regulate the proliferation of LUAD cells (33).

The TOP2A gene mainly acts on DNA transcription and replication enzymes, and its overexpression is associated with many tumors. Also, studies have shown that TOP2A can target CCNB1 and CCNB2 in LUAD to promote tumor proliferation and metastasis (38), leading to poor patient prognosis (39). In addition, CCNA2 and CCNB1 are proteins that can regulate the cell cycle. Their abnormal expression was correlated with the development and progression of NSCLC and poor patient prognosis (40-42). For example, CCNA2 can promote the epithelial-mesenchymal transition of NSCLC cells via integrin (43).

The TYMS gene is a key enzyme that maintains DNA synthesis and repair, and is considered one of the targets of new antifolate drugs such as pemetrexed (44). Its high expression was also linked to the poor prognosis of NSCLC (45).

The expression of MKI67 is related to the proliferation and growth of tumors. It is widely used as a proliferation and prognostic marker for research (46). It is also highly expressed in lung cancer (46-48).

In summary, 284 DEGs and 10 hub genes were obtained from the GEO database in our study. These 10 hub genes were strongly correlated with patient prognosis. Previous studies have demonstrated that the abnormal expression of hub genes was often strongly correlated with the development and prognosis of lung cancer, which is consistent with our results with LUAD. In the future, genes with high expression of EZH2, AURKA, CCNA2, CCNB1, FGF2, FOXM1, MKI67, TYMS, TOP2A, and CENPE may be used as targets for tumor therapy to prolong the survival time of LUAD patients. The relationship between FGF2 and immune microenvironment can also be explored to provide new strategies for the study and treatment of LUAD. However, some detailed mechanisms of hub genes in LUAD remain unclear and require further study in the future. In summary, we screened the GEO database to obtain hub genes that were highly correlated with LUAD. These genes may be potential biomarkers or therapeutic targets for LUAD.

We used some bioinformatics methods to find key genes via this study. However, due to the technical condition limitations, some other bioinformatics hot issues include weighted gene co-expression network analysis, competing endogenous RNA (ceRNA) network, DNA repair gene mutations and increased Neoantigen load and activated T cell infiltration were not carried out in this study, which will be the direction of our future research in LUAD. At the same time, as a result of the limitation of basic experimental conditions, this research was limited to the analysis of open database, and did not conduct functional research on hub genes. It is necessary to explore the mechanism of hub genes in the development of LUAD through functional research in the future.