Bioinformatics analysis identifies key genes of prognostic value in lung cancer

2.1 PCA, volcano plot and Venn diagram analysis

PCA (principal component analysis) was performed as previously depicted [7]. The GEO dataset GSE1037, which consists data on cDNA microarray containing 40,386 elements, was used to analyze the gene expression profiles of 38 surgically resected samples of lung neuroendocrine tumors and 11 SCLC (small cell lung cancer) cell lines (accessed via: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE1037) with the R package (limma analysis), after applying a batching correction. PCA was used to visualize the differences between the normal group (primary normal lung tissue, n = 19) and the cancer group, including PTC (primary typical carcinoid, n = 12), AC (adenocarcinoma, n = 12), PLCNC (primary large cell neuroendocrine carcinoma, n = 8), PSCLC (primary small cell lung carcinoma, n = 15) and PLCC (primary large cell carcinoma, n = 14).

In the volcano plot, the vertical lines represent genes with an absolute value of ∣log2(fold change)∣≥ 2, indicating upregulated or down-regulated DEGs, while the horizontal line indicates a significance level of p < 0.05. Significantly upregulated genes are denoted by red dots, while significantly downregulated genes are represented by blue dots. By using a Venn diagram, we identified 28 upregulated genes and 125 downregulated genes that were differentially expressed based on the intersection of the gene sets.

2.2 GO/KEGG pathway enrichment analysis and protein-protein interaction analysis

The expression of 14 key genes in lung tumor and non-tumor tissues was analyzed using the FactoMineR package in R, and the results were visualized using heat maps. GO-KEGG analysis (https://www.genome.jp/kegg/) was performed as previously described [8]. The online platform Metascape (https://metascape.org/) and R language version 4.3.1 (Beagle Scouts, Murray Hill, NJ, USA) were used for enrichment analysis of those above-mentioned intersecting differential genes. PPI (protein-protein interaction) analysis was performed on the STRING online platform (https://string-db.org/).

2.3 Correlation and cluster analysis of differential gene expression

The correlation of key gene expression between normal lung tissue and lung carcinoma tissue was analyzed using the clusterProfiler package of R language [9, 10].

2.4 GEO, GEPIA analysis and Kaplan-Meier curve analysis

The expression of 14 DEGs was analyzed by box plot from GSE1037. The GEPIA online software (http://gepia.cancer-pku.cn/) was used to analyze the expression of MFAP4, PDZD2, FBLN1, FBLN5, EFEMP1, KDR, S1PR1, CAV1, GRK5, EDNRA, EDNRB, CALCRL, PTGER4 and ADRB1 in different types of lung cancer. The prognosis of patients was assessed using the Kaplan-Meier curve tool to analyze and visualize the survival data for the investigated patients (http://kmplot.com/analysis/).

2.5 ROC curve analysis

ROC (Receiver Operating Characteristic) curves of differential genes were analyzed using the R language based on the TCGA-GTEX (Genotype-Tissue Expression) data. The type of data from TCGA comprised RNAseq data, while data from GEO were derived from cDNA microarrays.

2.6 Cell lines

A549 and SK-MES-1 cells were purchased from Procell (Wuhan, China), cultured in DMEM (Dulbecco’s Modified Eagle Medium) with 10% FBS (fetal bovine serum) and penicillin-streptomycin, then incubated at 37 ℃ in a 5% CO₂ (carbon dioxide) environment, following the manufacturer’s instructions.

2.7 Western blotting assay

Western blotting analysis was performed as previously depicted [11]. Briefly, proteins were extracted using RIPA lysis buffer, then denatured, electrophorized, transferred, immunoblotted and visualized using a chemiluminescent ECL reagent, followed by immunoblotting with the following antibodies: anti-mouse MFAP4 (1:1000, Santa, sc-398438, USA), anti-mouse p-ERK1/2 (1:1000, Abcam, ab278538, England), anti-mouse ERK1/2 (1:1000, Abcam, ab17942, England), anti-mouse p38 (1:1000, Santa, sc-7972, USA), anti-mouse p-p38 (1:1000, Santa, sc-166182, USA), anti-mouse p-JNK (1:1000, Santa, sc-6254, USA), anti-mouse JNK (1:1000, Santa, sc-7345, USA) and anti-mouse GAPDH (1:1000, Santa, sc-47724, USA).

2.8 Cell proliferation, migration and angiogenesis

Cell proliferation and migration assays were conducted following previously described methods [12]. The indicated cell lines were seeded into a 96-well plate to assess cell viability at various time points (0 h, 24 h, 48 h and 72 h). Additionally, the same cell lines were seeded into a transwell chamber after incubating for 24 h, following which the number of migrated cells was counted. For angiogenesis, human umbilical vein endothelial cells (HUVECs) were cultured in a medium obtained from lung cancer cells transfected with the designated plasmids for 48 h to conduct the tubule formation experiment. The cells were washed, detached, seeded and imaged as previously described [13].

2.9 Statistical analysis

The SPSS v22.0 software (IBM, Armonk, NY, USA) was used for statistical analyses. The Student’s t-test and one-way ANOVA (Analysis of variance) were performed, and data are presented as means ± SEM (standard error of the mean) of three independent experiments. A p-value of 0.05 or less was considered significant.

3.1 Analysis of differential genes DEGs in lung cancer based on GEO data

Bioinformatics methods were applied to analyze DEGs from GSE1037 to identify key genes related to the prognosis of lung cancer patients. After batch correction, the data were analyzed by PCA (Fig. 1A). The results showed significant differences between diverse kinds of lung cancer, including AC, PLCC, PLCNC, PSCLC and PTC, and the normal tissue group. We analyzed the DEGs between the lung cancer and normal group by volcano plot and identified 743 upregulated and 769 downregulated genes in PTC, 313 upregulated and 199 downregulated genes in AC, 454 upregulated and 501 downregulated expression genes in PLCNC, 325 upregulated and 1014 downregulated expression genes in PSCLC, and 537 upregulated and 416 downregulated genes in PLCC (Fig. 1B). From these, 28 upregulated and 125 downregulated intersection DEGs were identified among the AC, PTC, PLCC, PLCNC and PSCLC samples via Venn analysis (Fig. 1C).

Fig. 1.Analysis of differential genes DEGs in lung cancer based on GEO data. (A) PCA was used to analyze the difference between normal and lung cancer, including AC (adenocarcinoma), PTC (primary typical carcinoid), PLCC (primary large cell carcinoma), PLCNC (primary large cell lung carcinoma) and PSCLC (primary small cell lung carcinoma). (B) Volcano plot was used to assess DEGs between normal and different types of lung cancer. Blue: decreased expression genes, red: increased expression genes, grey: no significant change expression genes. (C) Venn plot analyzed the intersection genes among AC (adenocarcinoma), PTC (primary typical carcinoid), PLCC (primary large cell carcinoma), PLCNC (primary large cell lung carcinoma) and PSCLC (primary small cell lung carcinoma).

3.2 Differential gene enrichment analysis, cluster analysis and protein interaction network PPI analysis

Enrichment analysis was performed on the above intersection DEGs (Fig. 2A), and the results indicated that the screened DEGs were associated with angiogenesis, cell-substrate adhesion, circulatory system and locomotion process. KEGG enrichment analysis showed that the DEGs were mainly involved in the MAPK pathway. Additionally, the PPI network analysis identified 14 genes among the DEGs, namely MFAP4, PDZD2, FBLN1, FBLN5, EFEMP1, KDR, S1PR1, CAV1, GRK5, EDNRA, EDNRB, CALCRL, PTGER4 and ADRB1 (Fig. 2B).

Fig. 2.Differential gene enrichment analysis and protein interaction network PPI analysis. (A) GO-KEGG enrichment analysis of the intersection DEGs (differentially expressioned genes). For GO database analysis, the functions of genes were divided into three parts: cellular component (CC), molecular function (MF), and biological process (BP). (B) PPI (protein-protein intreaction) network analysis of the intersection DEGs was analyzed on the Metascape online platform (https://metascape.org/gp/index.html#/main/step1).

Next, we analyzed the expression of the selected 14 DEGs in tumor tissues and normal tissues of lung cancer patients and found that all were downregulated expression in lung carcinoma tissues compared to normal lung tissues (Fig. 3A). Moreover, the PPI analysis (Fig. 2B) demonstrated a consistent interaction pattern among the 14 DEGs (Fig. 3B). Notably, our analysis using R revealed that nearly all of the 14 DEGs displayed significant correlations with each other (Fig. 3C).

Fig. 3.Differential gene expression correlation and cluster analysis. (A) Heatmap analysis of the expression of MFAP4 (microfibril-associated protein 4), PDZD2 (PDZ domain containing 2), FBLN1 (fibulin 1), FBLN5 (fibulin 5), EFEMP1 (EGF containing fibulin extracellular matrix protein 1), KDR (kinase insert domain receptor), S1PR1 (sphingosine-1-phosphate receptor 1), CAV1 (caveolin 1), GRK5 (G protein-coupled receptor kinase 5), EDNRA (endothelin receptor type A), EDNRB (endothelin receptor type B), CALCRL (calcitonin receptor-like receptor), PTGER4 (prostaglandin E receptor 4), and ADRB1 (adrenoceptor beta 1) between normal and carcinoma tissues. (B) PPI among MFAP4, PDZD2, FBLN1, FBLN5, EFEMP1, KDR, S1PR1, CAV1, GRK5, EDNRA, EDNRB, CALCRL, PTGER4 and ADRB1 was analyzed on the STRING online platform. (C) Correlation of MFAP4, PDZD2, FBLN1, FBLN5, EFEMP1, KDR, S1PR1, CAV1, GRK5, EDNRA, EDNRB, CALCRL, PTGER4 and ADRB1 determined by R language analysis.

3.3 Analysis of differential gene expression in different subtypes of lung cancer from GEO and TCGA datasets

Based on the GEO chip data GSE1037, we analyzed the expression levels of the 14 DEGs in AC, PTC, PLCC, PLCNC and PSCLC by box plot (Fig. 4) and found no significant difference in PDZD2 between the different subtypes of lung cancer and the normal group. The expression of MFAP4 was downregulated in all kinds of lung cancer except for PTC, and PTGER4 was downregulated in different subtypes of lung cancer, excluding PTC and PLCNC. Additionally, FBLN1, FBLN5, EFEMP1, KDR, S1PR1, CAV1, GRK5, EDNRA, EDNRB, CALCRL and ADRB1 were all significantly decreased in AC, PTC, PLCC, PLCNC and PSCLC compared to the normal group (Fig. 4).

Fig. 4.Expression of key genes in different subtypes of lung cancer. Box plot was used to analyze the expression of MFAP4 (microfibril-associated protein 4), PDZD2 (PDZ domain containing 2), FBLN1 (fibulin 1), FBLN5 (fibulin 5), EFEMP1 (EGF containing fibulin extracellular matrix protein 1), KDR (kinase insert domain receptor), S1PR1 (sphingosine-1-phosphate receptor 1), CAV1 (caveolin 1), GRK5 (G protein-coupled receptor kinase 5), EDNRA (endothelin receptor type A), EDNRB (endothelin receptor type B), CALCRL (calcitonin receptor-like receptor), PTGER4 (prostaglandin E receptor 4), and ADRB1 (adrenoceptor beta 1) in different types of lung cancer from GSE1037.

Similarly, we used the GEPIA online platform to analyze the expression of 14 DEGs in LUAD (T = 483, N = 347) and LUSC (T = 486, N = 338). As shown in Fig. 5, for LUAD, which included 483 tumor cases and 347 normal cases, and LUSC, which comprised 486 tumor tissues and 338 normal tissues, all of the 14 DEGs were downregulated in the LUSC tissues compared to the normal group. Further, in LUAD tissues, all 14 DEGs showed significantly lower expression levels compared to the normal group, except for EDNRA.

Collectively, our data indicate that these 14 DEGs were downregulated in different subtypes of lung cancer.

Fig. 5.Analysis of differential gene expression in lung adenocarcinoma and lung squamous cell carcinoma based on the TCGA database. GEPIA (gene expression profiling interactive analysis) online software was used to analyze the expression of MFAP4 (microfibril-associated protein 4), PDZD2 (PDZ domain containing 2), FBLN1 (fibulin 1), FBLN5 (fibulin 5), EFEMP1 (EGF containing fibulin extracellular matrix protein 1), KDR (kinase insert domain receptor), S1PR1 (sphingosine-1-phosphate receptor 1), CAV1 (caveolin 1), GRK5 (G protein-coupled receptor kinase 5), EDNRA (endothelin receptor type A), EDNRB (endothelin receptor type B), CALCRL (calcitonin receptor-like receptor), PTGER4 (prostaglandin E receptor 4), and ADRB1 (adrenoceptor beta 1) in LUAD (lung adenocarcinoma) and LUSC (lung squamous cell carcinoma).

3.4 Influence of DEGs on lung adenocarcinoma and lung squamous cell carcinoma survival

Next, we analyzed the relationship between the expression of the 14 DEGs and survival rates. As shown in Fig. 6, for LUAD, a low expression of MFAP4, PDZD2, FBLN1, FBLN5, EFEMP1, KDR, S1PR1, CAV1, GRK5, EDNRA, EDNRB, CALCRL, PTGER4 and ADRB1 were related with poor survival. Comparatively, for LUSC, only MFAP4 and PTGER4 were positively correlated with survival rate, and there was no difference between other gene expressions and survival (Fig. 6 and Supplementary Fig. 1).

Fig. 6.The influence of differential genes on the survival and prognosis of lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC). Kaplan-Meier plot analysis was performed to determine the relationship between key gene expression and survival in LUAD and LUSC. p < 0.05.

3.5 Based on the TCGA database to analyze the value of differential genes in lung cancer diagnosis

The area under the ROC curve (AUC) of the nomogram of the training data set was AUC ≥0.9, with a CI >0.9, indicating good predictive ability in the TCGA cohort test, as a high AUC (≥0.9) represents high prediction accuracy. For LUAD diagnosis, MFAP4, FBLN1, FBLN5, EFEMP1, S1PR1, CAV1, GRK5, EDNRB, CALCRL, PTGER4 and ADRB1 exerted good predictive ability in the TCGA cohort test (Fig. 7). Comparatively, for LUSC diagnosis, only MFAP4, FBLN5, EFEMP1, KDR, S1PR1, CAV1, GRK5, EDNRB, CALCRL and ADRB1 exerted good predictive ability in the TCGA cohort test (Fig. 7). Further, the ROC curves of the nomogram of the training data set for PDZD2, EDNRA and KDR in LUAD, PDZD2, PTGER4, FBLN1 and EDNRA were AUC <0.9 (Supplementary Fig. 2).

Fig. 7.Significance of differential genes from the TCGA database in lung cancer diagnosis. ROC (receiver operating characteristic) curve analysis was used to assess the accuracy of each prognostic gene in different types of lung cancer. AUC (Area Under the ROC Curve) ≥0.9, CI (confidence interval): 0.8–0.9.

3.6 Effects of upregulated MFAP4 on angiogenesis and MAPK pathway in lung cancer cells

According to the GO/KEGG enrichment analysis, differential genes were significantly enriched in angiogenesis, cell adhesion movement and MAPK pathway, further confirming the role of MFAP4 at the cell level. From the above-screened genes, MFAP4 had the most differential expression and was the most significant marker associated with patient survival and was thus selected for validation. As shown in Fig. 8A, MFAP4 expression was decreased in A549 and SK-MES-1 cells compared with BEAS-2B cells. MFAP4 overexpression inhibited cell proliferation in A549 and SK-MES-1 cells compared to the control cell lines. The knockdown of MFAP4 induced cell proliferation in the indicated cell lines (Fig. 8B) [14]. Moreover, cell migration was decreased in A549 and SK-MES-1 cells following transfection with MFAP4, while it was increased in MFAP4 deficient A549 and SK-MES-1 cell lines (Fig. 8C). In addition, conditioned medium from MFAP4 overexpression reduced angiogenesis in HUVECs cell lines, and the knockdown of MFAP4 facilitated angiogenesis in HUVECs cell lines (Fig. 8D). Mechanistically, the overexpression of MFAP4 suppressed the phosphorylation of ERK1/2, p38 and JNK (c-Jun N-terminal kinases) in A549 and SK-MES-1 cells. MFAP4 deficiency induced the expression of p-ERK1/2, p-p38 and p-JNK in A549 and SK-MES-1 cell lines (Fig. 8E). Collectively, these data suggest that MFAP4 negatively regulates cell proliferation, migration and angiogenesis.

Fig. 8.Effects of upregulated MFAP4 on angiogenesis and MAPK pathway in lung cancer cells. (A) The expression of MFAP4 was assessed in BEAS-2B, A549 and SK-MES-1 via western blotting. (B) A549 and SK-MES-1 were divided into four groups (Control, MFAP4, shNC and shMFAP4) depending on different plasmids transfection. Cell proliferation of the cell line was assessed by MTT assay. (C) A549 and SK-MES-1 were transfected with the indicated plasmids, and transwell assay was performed to assess cell migration with quantitative data attachment. (D) Human umbilical vein endothelial cells (HUVECs) were cultured with a transfected medium of lung cancer cells to conduct tubule formation experiments with quantitative data attachment. (E) The related protein expression was measured by western blotting. Quantitative protein expression (MFAP4, p-ERK1/2, p-p38, and p-JNK (c-Jun N-terminal kinases)) data is attached. Error bars represent data from three independent experiments (mean ± SD). ##p < 0.01. **p < 0.01. #compared to the control group. *compared to the shNC group.