Research Paper Advance Articles
Key microRNAs and hub genes associated with poor prognosis in lung adenocarcinoma
- 1 Department of Thoracic Surgery, The First Affiliated Hospital of Zhengzhou University, Zhengzhou 450052, China
- 2 Biotherapy Center, The First Affiliated Hospital of Zhengzhou University, Zhengzhou 450052, China
Received: April 25, 2020 Accepted: October 8, 2020 Published: January 10, 2021https://doi.org/10.18632/aging.202337
How to Cite
Copyright: © 2021 Ye et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
In the study, we obtained 36 pairs of lung adenocarcinoma (LUAD) tissues and adjacent non-tumorous tissues. Then, we chose a specific hub-target gene of miRNA and used qRT-PCR to evaluate the expression of PECAM1. We found that the expression level of PECAM1 mRNA in LUAD was significantly lower than that in adjacent nontumor tissues (P<0.0001). Univariate and multivariate analyses were conducted on 481 LUAD patients from The Cancer Genome Atlas (TCGA) according to the Cox proportional hazard regression model to evaluate the impact of PECAM1 expression and other clinicopathological factors on survival. The results showed that the low expression of PECAM1 was an important independent predictor of poor overall survival (HR, 0.704; 95% CI, 0.518-0.957; P = 0.025). Based on the Tumor Immune Estimation Resource (TIMER) database, the relationship between PECAM1 expression and B cell, CD8+ T cell, CD4+ T cell, macrophage, neutrophil, and dendritic cell infiltration was weak in LUAD (P<0.01). In particular, a more significant positive correlation between PECAM1 expression and HLA-complex members, CD1C, NRP1, and ITGAX expression in dendritic cell was detected in LUAD. The mechanism which PECAM1 involved in the development of LUAD may be closely related to changes in the immune microenvironment.
The increasing incidence of LUAD may be due to increased smoking and air pollution . In addition, LUAD is commonly diagnosed at a late stage, when the cancer has already spread to nearby tissues and metastasized . Late diagnosis and inaccurate prognostication lead to delayed treatment and low survival rates . Therefore, there is an urgent requirement for effective biomarkers to accurately predict the outcomes of LUAD patients in order to provide appropriate adjuvant treatment.
Currently, many studies have found that a class of miRNAs could inhibit the stability and translation of mRNA by binding to specific gene sequences, which be used as potential biomarkers for different types of cancer. Due to the challenges in early diagnosis and prognosis of LUAD, it is important to identify the genetic factors involved in this disease. A variety of bioinformatics tools have been used to evaluate gene expression level, identify unique genes using RNA sequencing and next-generation sequencing data and explore potential significance in the development of tumors . Microarray analysis were widely used to recognize genetic changes in the initiation and development of tumors. In our study, we used mRNA and miRNA microarray data from the Gene Expression Omnibus (GEO) genomics datasets to analyze the differentially expressed genes (DEGs) as well as miRNAs (DEMs). In addition, we hope to provide a potential therapy targets and prognostic markers for LUAD.
Our analysis found that low PECAM1 expression in LUAD is an independent risk factor for poor overall survival (OS). The analysis of TIMER database indicated that the expression of PECAM1 significantly correlated with CD8+ T cells, B cells, macrophages, CD4+ T cells, neutrophils, and dendritic cells infiltration in lung adenocarcinoma. Our study provides potentially prognostic biomarkers with LUAD patients and highlights PECAM1 as a valuable prognostic biomarker in LUAD.
Analysis of DEGs and DEMs
For the GSE74190 dataset, a total of 55 miRNAs (28 upregulated and 27 downregulated miRNAs) were extracted (Figure 1A and Supplementary Table 1), and a total of 711, 597, and 896 genes were extracted from the GSE31908, GSE10072, and GSE43458 datasets, respectively (Figure 1B–1D). These miRNAs and genes had the differently expression between LUAD samples and normal samples. Among the DEGs, 254 overlapped; 49 of these were upregulated (Figure 1E), and 205 were downregulated (Figure 1F). Based on the log fold change (logFC) value, we show the top 25 upregulated and downregulated DEGs as a heat map in Figure 1G.
Figure 1. Volcano plot of gene expression profile data in LUAD and normal samples and the heatmap of the overlapping DEGs. (A) Volcano plot of GSE74190. (B) Volcano plot of GSE31908. (C) Volcano plot of GSE10072. (D) Volcano plot of GSE43458. (E) Venn diagram of the upregulated overlapping DEGs. (F) Venn diagram of the downregulated overlapping DEGs. (G) Heatmap of the overlapping DEGs. Green represents a low fold change (FC) value, and red represents a high fold change FC value. Each column represents one dataset, and each row represents one gene. The number in each rectangle represents the FC in LUAD samples compared with normal samples. The gradual color change from red to green represents the changing process from upregulation to downregulation.
The enrichment analysis in overlapping DEGs and DEMs
The biological processes (BP) results showed that overlapping DEGs were dramatically concentrated in cell adhesion, negative regulation of angiogenesis, and others in Figure 2A. Regarding molecular function (MF), the overlapping differentially expressed genes were significantly concentrated in heparin binding, extracellular matrix (ECM) binding, and others in Figure 2B. Regarding cellular components (CC), the overlapping DEGs were concentrated in proteinaceous ECM, membrane raft, and others in Figure 2C. In addition, the Kyoto Encyclopedia of Genes and Genomes (KEGG) results indicated that the overlapping DEGs were concentrated in pathways involved in ECM-receptor interactions, the PI3K-Akt signaling pathway, and the PPAR signaling pathway in Figure 2D.
Figure 2. Functional and pathway enrichment analyses of the overlapping DEGs in LUAD. (A) The BP analysis of DEGs. (B) The MF analysis of DEGs. (C) The CC analysis of DEGs. (D) The KEGG pathway analyses of DEGs. The x-axis represents the q value (−log10), and the y-axis represents the GO term. The GO terms were measured by the rich factor, q value, and number of genes enriched. The greater the rich factor, the greater is the degree of enrichment and the greater the P value [0, 1]. The brighter the color of red, the more significant is the term.
To understand thoroughly and fully the biological characteristics of the 55 DEMs in LUAD, we also fulfilled Gene Ontology (GO) annotation and biological pathway analysis. Regarding BP, the DEMs were concentrated in nucleoside, nucleotide, and others in Figure 3A. Regarding MF, the DEMs were significantly enriched in transcription factor activity in Figure 3B. Moreover, regarding CC were nucleus, cytoplasm, heterogeneous nuclear ribonucleoprotein, and actin cytoskeleton in Figure 3C. The most enriched biological pathway terms for the DEMs were related to the ERBB receptor signaling network, nectin adhesion pathway, proteoglycan syndecan-mediated signaling events and others in Figure 3D.
Figure 3. Functional and pathway enrichment analyses of the DEMs in LUAD. (A) The BP analysis of DEMs. (B) the MF analysis of DEMs. (C) The CC analysis of DEMs. (D) The biological pathway analysis of DEMs. The upper x-axis represents the P value (−log10), and the lower x-axis represents the percentage of genes (blue). The y-axis represents the GO term. Yellow represents a P = 0.05 as reference, and red represents the specific P value. The longer the rectangular zone, the smaller is the P value.
Construct the protein-protein interaction (PPI) network of the overlapping DEGs
To assess the interactions of the overlapping DEGs, a PPI network was established in the STRING online website. In total, the network of the DEGs comprises 221 nodes and 607 edges (Figure 4). The more characteristic properties of node degree, the more significant influence the gene has on maintaining the stability of the PPI network. We used Molecular Complex Detection (MCODE) to evaluate functional modules. Eight hub genes were identified with degrees > 14: SPP1 (degree = 29), VWF (degree = 29), PECAM1 (degree = 26), COL1A1 (degree = 24), CDH5 (degree = 23), END1 (degree = 21), COL3A1 (degree = 19), and TEK (degree = 14) (Table 1).
Figure 4. PPI network analysis of the overlapping DEGs. Red nodes, upregulated genes; green nodes, downregulated genes.
Table 1. Hub genes for hypomethylated, highly expressed genes ranked in cytoHubba.
|Catelogy||Rank methods in cytoHubba|
|Gene symbol top 15||SPP1||SPP1||SPP1||SPP1||SPP1||PECAM1|
We used Kaplan-Meier (KM) plotter to evaluate the 55 DEMs and 8 hub genes. Among the DEMs, our results showed that low expression level of hsa-miR-126, hsa-miR-218, hsa-miR-30a, hsa-miR-145, hsa-miR-1, hsa-miR-195, hsa-miR-551b, hsa-miR-497, and hsa-miR-101 (P < 0.05) was associated with poor OS in LUAD (Figure 5A–5I). In contrast, high expression level of hsa-miR-215, hsa-miR-31, hsa-miR-196a, and hsa-miR-21 (P < 0.05) was also associated with poor OS in LUAD (Figure 5J–5M). Regarding the hub genes, high expression level of SPP1, COL1A1, and COL3A1 and low expression level of VWF, PECAM1, EDN1, CDH5, and TEK were associated with the poor OS in LUAD (P < 0.05) in Figure 6.
Figure 5. Overall survival analyses of DEMs. KM curves depicting OS for LUAD patients with high and low expression of (A) hsa-miR-126, (B) hsa-miR-218, (C) hsa-miR-30a, (D) hsa-miR-145, (E) hsa-miR-1, (F) hsa-miR-195, (G) hsa-miR-551b, (H) hsa-miR-497, (I) hsa-miR-101, (J) hsa-miR-215, (K) hsa-miR-31, (L) hsa-miR-21, and (M) hsa-miR-198a.
Figure 6. Prognostic value of eight DEGs in LUAD patients. Prognostic value of (A) SPP1 (log-rank P = 0.0015), (B) COL3A1 (log-rank P = 0.0018), (C) OL1A1 (log-rank P = 4e-08), (D) VWF (log-rank P = 0.045), (E) CDH5 (log-rank P = 0.0034), (F) TEK (log-rank P = 8.9e-10), (G) PECAM1 (log-rank P = 0.0036), and (H) EDN1 (log-rank P = 0.018).
Transcriptional expression of the hub genes and correlation analysis
As shown in Figure 7, based on the Gene Expression Profiling Interactive Analysis (GEPIA) website, mRNA levels were upregulated for three genes and notably downregulated for five genes in LUAD. Moreover, using the Human Protein Atlas (HPA) website, we examined the protein expression level of the hub genes and observed that protein expression of SPP1, COL1A1, and COL3A1 was noticeably upregulated in tumor samples, whereas the protein expression of VWF, PECAM1, EDN1, CDH5, and TEK was noticeably downregulated (Figure 8). In the upregulated hub genes, increased expression of SPP1, COL1A1, and COL3A1 was strongly associated each other in LUAD. In the downregulated hub genes, increased expression in PECAM1, VWF, CDH5, and TEK was strongly associated each other in LUAD (Figure 9).
Figure 7. Expression levels of eight hub genes in human lung adenocarcinoma. (A) SPP1; (B) COL3A1; (C) COL1A1; (D) VWF; (E) CDH5; (F) TEK; (G) PECAM1; and (H) EDN1. The gray and red boxes represent normal and cancer tissues, respectively. The expression data are first log2(TPM+1) transformed for differential analysis, and the log2FC is defined as median (Tumor, T)-median (Normal, N). *P < 0.05, T vs N. Abbreviations: FC, fold change; TPM, transcript per kilobase million.
Figure 8. Validation of eight hub genes using The Human Protein Atlas database. Expression of (A) SPP1, (B) COL3A1, (C) COL1A1, (D) VWF, (E) CDH5, (F) TEK, (G) PECAM1, and (H) EDN1.
Figure 9. Correlation analysis of upregulated hub genes (SPP1, COL3A1, and COL1A1) and downregulated hub genes (PECAM1, VWF, CDH5, and TEK) in LUAD. (A) SPP1 and COL1A1, (B) SPP1 and COL3A1, (C) COL1A1 and COL3A1, (D) PECAM1 and VWF, (E) PECAM1 and CDH5, (F) PECAM1 and TEK, (G) VWF and CDH5, (H) VWF and TEK, and (I) CDH5 and TEK. TPM, transcript per million. The x-axis represents the TPM of the hub gene MYC (log2). The expression levels of upregulated hub genes (SPP1, COL3A1, and COL1A1) were positively correlated with each other. The expression levels of downregulated hub genes (PECAM1, VWF, CDH5, and TEK) were positively correlated with each other.
Identification of hub genes by ROC analysis
We used the receiver operating characteristic (ROC) curve to predict accuracy of hub genes (Figure 10). The upregulated gene SPP1 and downregulated genes VWF, CDH5, TEK, and PECAM1 had AUC values > 0.97. The ROC curves of the genes had excellent predictive performance and were able to discriminate LUAD from healthy tissues.
The network of DEMs and predicted targets is presented in Figure 11. The interactive relationship between hub genes and miRNAs is presented in Supplementary Figure 1 and Supplementary Table 2. Notably, hsa-miR-31 targeted COL1A1, VWF, and TEK; hsa-miR-196a targeted COL1A1 and COL3A1; hsa-miR-215 targeted TEK; hsa-miR-218 targeted COL1A1; hsa-miR-1 targeted EDN1; hsa-miR-145 targeted PECAM1; hsa-miR-195 targeted VWF, COL1A1, and COL3A1; hsa-miR-101 targeted CDH5; and hsa-miR-497 targeted VWF and COL1A1.
PECAM1 mRNA expression in LUAD tumor and normal samples
In four analyses, the mRNA levels of PECAM1 in tumor samples were lower than those in normal samples (P < 0.05; Figure 12A–12D). In addition, the qRT-PCR was conducted to demonstrate the transcript level of PECAM1 which was lower in our LUAD tumor samples (P <0.0001) in Figure 12E.
Figure 12. Comparison of PECAM1 mRNA levels in normal lung tissues and LUAD tissues across four analyses of LUAD. (A–D) Each study showed that the PECAM1 mRNA expression levels in normal lung tissues were significantly higher than those in LUAD tissues. (E) qRT-PCR analysis of PECAM1 mRNA expression in 47 pairs of LUAD tissues and adjacent nontumor tissues. (F) Kaplan-Meier curve of the relationship between PECAM1 mRNA expression and the prognosis of LUAD patients based on the TCGA database.
Univariate and multivariate analyses of the impact of PECAM1 expression on survival outcomes
The expression of PECAM1 was significantly related to survival, univariate and multivariate Cox regression model were analyzed to investigate the effect of PECAM1 expression of 481 LUAD patients and other pathologic and clinical variables on survival. Based on univariate analysis, we found that pathological stage (Hazard Ratio (HR), 2.446; 95% Confidence Interval (CI), 1.778-3.363; P < 0.001), T stage (HR, 2.295; 95% CI, 1.554-3.388; P < 0.001), N stage (HR, 2.468; 95% CI, 1.830-3.328; P < 0.001), M stage (HR, 1.926; 95% CI, 1.086-3.416; P = 0.025), and PECAM1 expression (HR, 0.634; 95% CI, 0.473-0.874; P = 0.005) were the significant factors of survival (Table 2). Moreover, we also analyzed the expression of PECAM1 and other clinicopathological variables (T stage, N stage, and M stage, pathological stage, age) by using multivariate Cox analysis. Above results suggested that the low expression of PECAM1was a significant independent predictor of poor overall survival (HR, 0.704; 95% CI, 0.518-0.957; P = 0.025; Table 2). According to the previously clinical follow-up data provided by the TCGA database, patients with the low expression of PECAM1 had poor OS (P = 0.0035, Figure 12F).
Table 2. Univariate analysis and multivariate analysis of correlation of PECAM1 expression with LUAD patients.
|Clinicopathological features||Univariate analysis||Multivariate analysis|
|HR||95% CI||P||HR||95% CI||P|
|Bold values indicate P < 0.05. HR, hazard ratio; CI, confidence interval.|
The PECAM1 expression was correlated to immune cell infiltration in LUAD
The analysis of TIMER website observed that the expression of PECAM1 was associated with infiltration of several immune cells in LUAD (P < 0.01; Figure 13). Moreover, we investigated the relationship between PECAM1 expression and gene markers of immune cells. The results mentioned above demonstrated that the correlations between the expression of PECAM1 and gene markers of B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and DCs were determined (Figure 14A–14F). In particular, we detected a positive correlation between the expression of PECAM1 and gene markers of DCs in LUAD (HLA-DPA1: cor = 0.481, P = 0e+∞; HLA-DPB1: cor = 0.445, P = 2.1e-26; HLA-DRA: cor = 0.407, P = 0e+∞; HLA-DQB1: cor = 0.326, P = 2.95e-14; NRP1: cor = 0.259, P = 2.74e-09; ITGAX: cor = 0.475, P = 2.34e-30; CD1C: cor = 0.353, P = 1.38e-16; Figure 14F). Thus, we conjectured that the low expression of PECAM1 might exert influence on anti-tumor immunity in the LUAD microenvironment.
Figure 13. Correlation of PECAM1 expression with immune cell infiltration levels in LUAD. Tumor-infiltrating immune cells included B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and DCs. Gene expression levels against tumor purity are displayed in the left-most panel.
Figure 14. Correlation of PECAM1 expression with gene markers of tumor-infiltrating immune cells in LUAD. (A) Correlation with gene markers of B cells in LUAD. (B) Correlation with gene markers of CD8+ T cells in LUAD. (C) Correlation with gene markers of CD4+ T cells in LUAD. (D) Correlation with gene markers of neutrophils in LUAD. (E) Correlation with gene markers of macrophages in LUAD. (F) Correlation with gene markers of DCs in LUAD.
LUAD is a highly malignant tumor that is particularly prone to brain and bone metastasis, which negatively affects patient survival . LUAD is the main cause of cancer death in large part because of the challenges faced in detecting the cancer at an early stage. Lung cancer disproportionately affects the elderly, and management of this complex and lethal disease in older patients is challenging . The senescence and frailty of the elderly adversely affect the management of cancer and patient outcome . Earlier diagnosis and treatment of LUAD may reduce deaths among older patients. Therefore, identifying prognostic markers with high specificity and sensitivity is crucial. Many studies indicated that miRNA and RNA can regulate the expression of genes and involved in many biological processes of human malignant tumors . In particular, recent studies have showed that distinct miRNA and mRNA expression profiles are related to the development of LUAD .
In present research, the expression of mRNA and miRNA data downloaded from the GEO database yielded 254 DEGs and 55 DEMs. The enrichment analysis of DEGs were involved in ECM-receptor interaction pathways, the PI3K-Akt signaling pathway, and the PPAR signaling pathway. The DEMs were associated with ERBB receptor signaling network and the IGF1 pathway. Our study and others have shown that lung cancer progression depended not only on the presence of driver mutations within the cancer cells, but also on their interactions with the cellular and ECM environment . Recent evidence suggested that altered metabolic pathways play a crucial role in both cancer occurrence and development, and established link between PPAR signaling, metabolism, and cancer currently represents one of the most active research fields in the literature [11, 12].
We identified eight hub genes 13 and DEMs that were related to the prognosis of patients in LUAD. High expression of four DEMs (miR-215, miR-31, miR-196a, and miR-21) and three hub genes (SPP1, COL1A1, and COL3A1) was associated with worse OS, whereas low expression of nine DEMs (hsa-miR-126, hsa-miR-30a, hsa-miR-218, hsa-miR-1, hsa-miR-145, hsa-miR-195, hsa-miR-551b, hsa-miR-101, and hsa-miR-497) and five hub genes (VWF, PECAM1, EDN1, CDH5, and TEK) was associated with worse OS. The analysis of ROC showed that hub gene levels had excellent predictive performance and were able to discriminate LUAD from healthy tissue. The up-regulated SPP1, COL1A1, and COL3A1 in LUAD were positively correlated with each other, while the down-regulated PECAM1, VWF, EDN1, CDH5, and TEK were in LUAD also positively correlated with each other. Guo et al studied the expression and clinical significance of SPP1 in LUAD and found that SPP1 influenced the initiation and progression of LUAD, which was considered as a new target for molecular treatment . Similarly, we tried to determine the mechanisms of PECAM1 in LUAD. We found that PECAM1 mRNA expression was lower in LUAD. Through using qRT-PCR, it demonstrated that the expression of PECAM1 were down-regulated in tumor tissues. Plus, low PECAM1 expression was significantly associated with clinicopathological features related to development of LUAD. Univariate and multivariate survival analyses identified the low PECAM1 expression was an important independent predictor of unfavorable LUAD prognosis, which demonstrate its potential as a biomarker for LUAD.
In the microenvironment of lung adenocarcinoma, tumor-infiltrating immune cells plays a significant role in cancer progression, which not only inhibit tumor progression by attacking and killing lung cancer cells, but also promote tumor progression by changing the tumor microenvironment . Thus, we evaluated the relationship between PECAM1 expression and infiltration of immune cells. The positive correlation was seen between PECAM1 expression and CD4+ T cells, B cells, CD8+ T cells, neutrophils, macrophages, and DCs. We further analyzed the relationship between PECAM1 expression and molecular markers of tumor-infiltrating immune cells. The present study indicated that the expression of PECAM1 was related to molecular markers of B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells. Especially, PECAM1 expression was significantly correlated with HLA-complex member (ITGAX, CD1C, NRP1) in dendritic cells. Due to DCs are vital antigen-presenting cells (APCs) triggering T-cell–mediated immune response , Thus, the functional impairment of dendritic cells could influence the body’s anti-tumor immunity reaction. While fundamental studies should be established causally, the result indicated that low expression of PECAM1 could compromise DC-mediated antitumor immune response in LUAD.
Dysregulated miRNAs have shown promise as noninvasive diagnostic and prognostic biomarkers for LUAD. Numerous studies indicated that altered expression of miRNAs is related to various cancers, including lung cancer . In our study, high expression of four DEMs (has-miR-31, has-miR-21, has-miR-196a, and has-miR-215) and low expression of nine DEMs (hsa-miR-126, hsa-miR-218, hsa-miR-30a, hsa-miR-145, hsa-miR-1, hsa-miR-195, hsa-miR-551b, hsa-miR-497, and hsa-miR-101) were associated with worse OS in LUAD. Liu et al found that miR-126-3p was crucial in the invasion of NSCLC by targeting CCR1 . Shi et al determined that downregulation of miR-218 was noted to serve an essential role in epithelial-mesenchymal transition (EMT) and cancer metastasis of lung cancer by targeting Slug/ZEB2 signaling . MiR-30a played a crucial role in the occurrence of NSCLC . Liu et al reported that miR-145 was downregulated in tumor tissues and that miR-145 overexpression resulted in downregulation of N-cadherin, vimentin, and E-cadherin, suggesting decreased EMT activity . Korde et al found that levels of miR-1 were lower in NSCLC than normal tissue and that overexpression of miR-1 reduced tumor growth and angiogenesis . Liang et al found that downregulating miR-195 simultaneously stimulated cell proliferation and inhibited cell apoptosis in LUAD . Lin et al identified miR-551b as a cancer-specific miRNA that may predict the progression and prognosis of LUAD . Cui et al indicated that inhibition of miR-101 promotes progression of NSCLC through activation of the Wnt/β-catenin signaling pathway . Yin et al reported that miR-145 and miR-497 were downregulated in lung cancer cell lines and that miR-145 and miR-497 overexpression repressed TGF-β-induced EMT and the migration and invasion of cancer cell . Recent researches have showed that miR-31 modulates the NSCLC cell cycle by targeting hMLH1 . Bao et al revealed that serum miR-196a-5p may be a valuable noninvasive biomarker for the clinical diagnosis of NSCLC . In our study, all miRNAs which associated with prognosis have been reported previously in lung cancer except miR-215. And maybe they could be the useful biomarker predictive of prognosis in LUAD.
PECAM1 can code for CD31, belonging to the immunoglobulin superfamily of adhesion molecule, which can maintain and restore vascular integrity by the wingless-related integration site signaling pathway . Yu et al analyzed the related microarray datasets (GSE9804) and identified that PECAM1 played a crucial role in the tumorigenesis of LUAD via regulation of vascular endothelial growth factor(VEGF) expression . The mRNA transcription and the protein expression levels, we found that PECAM1 expression was significant down-regulated in LUAD tumor tissues. Importantly, we found that low expression of PECAM1 was a significant independent predictor of poor OS by univariate and multivariate survival analyses in LUAD. SPP1 has also been demonstrated in mice and humans to affect tumor cell migration, adhesion, and invasion, indicating its possible use as a diagnostic biomarker and potential therapeutic target in cancer . Our enrichment analysis indicated that SPP1 is enriched in cell adhesion, the PI3K-Akt signaling pathway, and ECM-receptor interactions, which play an important role in LUAD. In addition, we confirmed that COL1A1 expression was higher in LUAD tissues. COL1A1 was reported to be enriched in various biological processes such as cell proliferation, metastasis, invasion, and angiogenesis . Wei et al analyzed the GSE75037 microarray dataset and determined that COL1A1, CTGF, COL1A2, COL3A1, and BGN were involved in cell adhesion, the PI3K-Akt signaling pathway, and ECM-receptor interaction . Above study was consistent with our findings. Besides its essential role in hemostasis, there is growing evidence that VWF has antitumor effects, mainly related to angiogenesis and apoptosis . EDN1 were implicated in cancer progression, apoptosis, EMT, stromal reaction, and tumor invasion, spread, and drug resistance . CDH5 is an adhesion molecule and member of the selectin family that contributes significantly to tumorigenesis and tumor progression . TEK is a receptor tyrosine kinase that is expressed in endothelial cells. Dysregulated TEK expression has also been observed in several cancer, indicating that TEK expression may have prognostic value in these cancers .
In sum, these data supported that DEMs and hub genes may have clinical utility as prognostic markers in LUAD. Although there are some limitations in the current research, the combination of multichip analysis and bioinformatics analysis has significant advantages in identifying potential diagnostic biomarkers and therapeutic targets of cancers.
Therefore, the present research provided the potential prognosis biomarkers as well as future therapeutic targets for LUAD and indicated that PECAM1, in particular, may be a novel biomarker of survival that provided a novel diagnostic biomarkers and therapeutic targets for the treatment of LUAD and may help to improve the therapeutic efforts in the future. However, further studies were needed to explore the molecular mechanisms and biological functions of these microRNAs and hub genes and to evaluate whether they can serve as biomarkers or therapeutic targets in LUAD patients.
Materials and Methods
LUAD dataset preprocessing
We downloaded GSE31908, GSE10072, and GSE43458 datasets and GSE74190 miRNA expression data from GEO by searching using the following key terms: (“microRNA” OR “miRNA” OR “mRNA”) AND (“Lung adenocarcinoma” OR “LUAD”). The GSE31908 dataset included 20 normal lung tissues and 30 LUAD tissues. GSE10072 comprised 107 samples, including 49 normal samples and 58 tumors. GSE43458 contained 110 samples, including 80 LUAD tissues and 30 normal tissues. GSE74190 comprised 80 samples, including 36 tumors and 44 normal tissues.
The GEO2R could compare different groups in GEO datasets and identify DEGs and DEMs. Adjusted P < 0.05 and an absolute logFC > 1 for the DEGs and the DEMs were used as the cutoff criteria. The Venn plots were used to show the overlapping DEGs that were up or down-regulated in the GSE31908, GSE10072, and GSE43458 datasets.
Clinical specimen collection
In total, we obtained 36 patients of LUAD tissues and corresponding normal tissues from the First Affiliated Hospital of Zhengzhou University between April 2019 and January 2020. These patient samples had never received chemoradiotherapy prior to surgical resection. Written informed consents were received from all patients, and this study was approved by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University.
The enrichment analysis in overlapping DEGs and DEMs
Database for Annotation, Visualization, and Integrated Discovery (DAVID) has integrated many biological data analytical tools. We visualized the key BP, MF, CC and KEGG pathways of DEGs by using the DAVID. Similarly, we also performed a functional enrichment analysis of the DEMs by FunRich, which is also primarily used for the functional enrichment analysis of genes and proteins. Visualization of GO terms and KEGG results was performed by bubble diagrams and bar charts.
Construct the PPI network of the overlapping DEGs
With setting an interaction score ≥0.4, we analyzed the PPI network of DEGs by STRING website. Then, we downloaded the network relationship file, and the top 15 genes were identified in accordance with Cytoscape 3.7.1 and MCODE, an app of Cytoscape software (MCC, MNC, Degree, EPC, Closeness, Radiality ranking of MCODE). We used the intersection of these sets to obtain the final hub genes.
The survival analysis of DEMs and hub genes
Respectively, KM plotter was used to perform an OS analysis for DEMs and hub genes in 513 and 866 LUAD patient samples. Based on the gene transcriptional expression of a given gene, we separated patients into high and low groups and performed KM plots by the KM plotter. In addition, the chart showed the HR with the 95% CI and the log-rank P value, and the number at risk is displayed below the curves.
The expression and correlation analysis of the hub genes
GEPIA, a well-developed interactive online server applied to a standard processing pipeline, including 9736 tumor and 8587 normal tissues mRNA data. We used GEPIA to appraise the mRNA expression levels of hub genes in LUAD and normal tissues. A correlation analysis presents pairwise genes using the Pearson, Spearman, and Kendall correlation statistics based on GTEx and/or TCGA expression data in GEPIA.
The ROC curves were carried out to assess the specificity of the DEGs for LUAD diagnosis. The ROC curves were plotter and the AUC values were calculated using the ROC package for R.
Construct the miRNA–mRNA network
We used the miRecords database to predict the target genes of DEMs using 11 established algorithms. Only target genes of DEMs predicted by at least three prediction databases were accepted. Next, the regulatory network of the predicted genes and miRNAs that targeted them were visualized in Cytoscape.
Verification of PECAM1 by the GEO database
By using the Mann-Whitney U test, we validated mRNA expression levels of PECAM1 between LUAD tissues and normal tissues in GEO databases. The details of GEO databases are as follows: GSE40791, GSE32863, GSE75037, and GSE7670. Above datasets were obtained from the GEO datasets. The preprocessing of each dataset adopts the preprocessing results of the original author.
Verification of PECAM1 expression by qRT-PCR
Total RNAs were isolated from LUAD tissues using TriZol reagent (TakaRa, Japan). Then total RNAs was reverse-transcribed using the PrimeScript RT reagent Kit with gDNA Eraser (TakaRa, Japan). PECAM1 expression was validated by qRT-PCR using SYBR Premix Ex Taq II (TakaRa, Japan) on the LightCycler 480 II Real-Time PCR System (Roche, Basel, Switzerland). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was amplified as an internal control. The primer sequences (sunya Biotech, zhengzhou, China) were as follows:
GAPDH forward: GGAGCGAGATCCCTCCAAAAT
GAPDH reverse: GGCTGTTGTCATACTTCTCATGG
Univariate and multivariate cox regression analyses
The Cox regression models were used to calculate univariate and multivariate analysis. The factors that had a significant effect on survival were identified. We calculated the value of the clinicopathological factors and PECAM1 expression on survival by the univariate and multivariate analysis. Firstly, we downloaded mRNA expression and clinical data from TCGA. Then we screened the whole survival and gene expression data. Finally, the data of 481 patients were obtained, which were analyzed by Cox regression models. The median expression was used as a cutoff value, the patients at the top 50% expression level of PECAM1 were assigned to the high group and the other patients were assigned to the low group. We used R software’s survival package to analyze and visualize the data.
The online tool TIMER is a resource for comprehensive analysis of immune infiltration across multiple cancer types. We used TIMER on LUAD sample data to assess these infiltrates. Then, the relationship between expression of PECAM1 and immune infiltration was demonstrated. Furthermore, the database was used to explore the association between PECAM1 expression and gene markers of immune cells in LUAD.
We used the Mann-Whitney U test to test the differential expression of PECAM1 in tumor and corresponding normal samples. Univariate and multivariate survival analysis were carried out by using the Cox proportional hazards regression models. All statistical analyses were used by IBM SPSS statistics (version 26.0) and R software (version 4.0.2). P < 0.05 was considered statistically significant.
LUAD: lung adenocarcinoma; miRNA: microRNA; GEO: Gene Expression Omnibus; DEGs: differentially expressed genes; DEMs: differentially expressed miRNAs; DAVID: Database for Annotation, Visualization, and Integrated Discovery; PPI: protein-protein interaction; OS: overall survival; NSCLCs: non–small-cell lung cancers; mRNA: messenger RNA; TCGA: The Cancer Genome Atlas; DCs: dendritic cells; logFC: log fold change; BP: biological processes; MF: molecular function; ECM: extracellular matrix; CCs: cellular components; KEGG: Kyoto Encyclopedia of Genes and Genomes; GO: Gene Ontology; MCODE: Molecular Complex Detection; KM: Kaplan-Meier; GEPIA: Gene Expression Profiling Interactive Analysis; HPA: Human Protein Atlas; ROC: receiver operating characteristic; AUC: area under curve; HR: hazard ratio; TIMER: Tumor Immune Estimation Resource; qRT-PCR: quantitative reverse transcription PCR.
Yu Qi substantially contributed to the concept and design of the study. Guanchao Ye, Lan Huang, and Yafei Liu collected data on lung adenocarcinoma from several different databases. Chunyang Zhang, Yinliang Sheng, and Bin Wu used bioinformatics to analyze the data. Lu Han, Bo Dong, and Chunli Wu co-wrote the manuscript. Yu Qi examined all the data and modified the manuscript. All authors read and approved the final manuscript.
Conflicts of Interest
The authors declare no conflicts of interest.
This study was supported by the National Natural Science Foundation of China (Grant No. 81773045) and the Joint Construction Project of Henan Medical Science and Technology Project, China (Grant No. LHGJ20190301).
- 1. Li H, Li QD, Wang MS, Li FJ, Li QH, Ma XJ, Wang N. Smoking and air pollution exposure and lung cancer mortality in Zhaoyuan County. Int J Hyg Environ Health. 2013; 216:63–70. https://doi.org/10.1016/j.ijheh.2012.06.003 [PubMed]
- 2. Du R, Shen W, Liu Y, Gao W, Zhou W, Li J, Zhao S, Chen C, Chen Y, Liu Y, Sun P, Xiang R, Shi Y, Luo Y. TGIF2 promotes the progression of lung adenocarcinoma by bridging EGFR/RAS/ERK signaling to cancer cell stemness. Signal Transduct Target Ther. 2019; 4:60. https://doi.org/10.1038/s41392-019-0098-x [PubMed]
- 3. Nasim F, Sabath BF, Eapen GA. Lung cancer. Med Clin North Am. 2019; 103:463–73. https://doi.org/10.1016/j.mcna.2018.12.006 [PubMed]
- 4. Alshamsan A, Khan S, Imran A, Aljuffali IA, Alsaleh K. Prediction of Chlamydia pneumoniae protein localization in host mitochondria and cytoplasm and possible involvements in lung cancer etiology: a computational approach. Saudi Pharm J. 2017; 25:1151–57. https://doi.org/10.1016/j.jsps.2017.05.007 [PubMed]
- 5. Denisenko TV, Budkevich IN, Zhivotovsky B. Cell death-based treatment of lung adenocarcinoma. Cell Death Dis. 2018; 9:117. https://doi.org/10.1038/s41419-017-0063-y [PubMed]
- 6. Vitale SG, Capriglione S, Zito G, Lopez S, Gulino FA, Di Guardo F, Vitagliano A, Noventa M, La Rosa VL, Sapia F, Valenti G, Rapisarda AM, Peterlunger I, et al. Management of endometrial, ovarian and cervical cancer in the elderly: current approach to a challenging condition. Arch Gynecol Obstet. 2019; 299:299–315. https://doi.org/10.1007/s00404-018-5006-z [PubMed]
- 7. Gajra A, Akbar SA, Din NU. Management of lung cancer in the elderly. Clin Geriatr Med. 2016; 32:81–95. https://doi.org/10.1016/j.cger.2015.08.008 [PubMed]
- 8. Chen X, Xu X, Pan B, Zeng K, Xu M, Liu X, He B, Pan Y, Sun H, Wang S. miR-150-5p suppresses tumor progression by targeting VEGFA in colorectal cancer. Aging (Albany NY). 2018; 10:3421–37. https://doi.org/10.18632/aging.101656 [PubMed]
- 9. Li Y, Zhang H, Fan L, Mou J, Yin Y, Peng C, Chen Y, Lu H, Zhao L, Tao Z, Chen J, Wang Y, Qi X, et al. MiR-629-5p promotes the invasion of lung adenocarcinoma via increasing both tumor cell invasion and endothelial cell permeability. Oncogene. 2020; 39:3473–88. https://doi.org/10.1038/s41388-020-1228-1 [PubMed]
- 10. Götte M, Kovalszky I. Extracellular matrix functions in lung cancer. Matrix Biol. 2018; 73:105–21. https://doi.org/10.1016/j.matbio.2018.02.018 [PubMed]
- 11. Vitale SG, Laganà AS, Nigro A, La Rosa VL, Rossetti P, Rapisarda AM, La Vignera S, Condorelli RA, Corrado F, Buscema M, D’Anna R. Peroxisome proliferator-activated receptor modulation during metabolic diseases and cancers: master and minions. PPAR Res. 2016; 2016:6517313. https://doi.org/10.1155/2016/6517313 [PubMed]
- 12. Laganà AS, Vitale SG, Nigro A, Sofo V, Salmeri FM, Rossetti P, Rapisarda AM, La Vignera S, Condorelli RA, Rizzo G, Buscema M. Pleiotropic actions of peroxisome proliferator-activated receptors (PPARs) in dysregulated metabolic homeostasis, inflammation and cancer: current evidence and future perspectives. Int J Mol Sci. 2016; 17:999. https://doi.org/10.3390/ijms17070999 [PubMed]
- 13. Guo Z, Huang J, Wang Y, Liu XP, Li W, Yao J, Li S, Hu W. Analysis of expression and its clinical significance of the secreted phosphoprotein 1 in lung adenocarcinoma. Front Genet. 2020; 11:547. https://doi.org/10.3389/fgene.2020.00547 [PubMed]
- 14. Merlo A, Dalla Santa S, Dolcetti R, Zanovello P, Rosato A. Reverse immunoediting: when immunity is edited by antigen. Immunol Lett. 2016; 175:16–20. https://doi.org/10.1016/j.imlet.2016.04.015 [PubMed]
- 15. Stankovic B, Bjørhovde HA, Skarshaug R, Aamodt H, Frafjord A, Müller E, Hammarström C, Beraki K, Bækkevold ES, Woldbæk PR, Helland Å, Brustugun OT, Øynebråten I, Corthay A. Immune cell composition in human non-small cell lung cancer. Front Immunol. 2019; 9:3101. https://doi.org/10.3389/fimmu.2018.03101 [PubMed]
- 16. Wang K, Chen R, Feng Z, Zhu YM, Sun XX, Huang W, Chen ZN. Identification of differentially expressed genes in non-small cell lung cancer. Aging (Albany NY). 2019; 11:11170–85. https://doi.org/10.18632/aging.102521 [PubMed]
- 17. Liu R, Zhang YS, Zhang S, Cheng ZM, Yu JL, Zhou S, Song J. MiR-126-3p suppresses the growth, migration and invasion of NSCLC via targeting CCR1. Eur Rev Med Pharmacol Sci. 2019; 23:679–89. https://doi.org/10.26355/eurrev_201901_16881 [PubMed]
- 18. Shi ZM, Wang L, Shen H, Jiang CF, Ge X, Li DM, Wen YY, Sun HR, Pan MH, Li W, Shu YQ, Liu LZ, Peiper SC, et al. Downregulation of miR-218 contributes to epithelial-mesenchymal transition and tumor metastasis in lung cancer by targeting Slug/ZEB2 signaling. Oncogene. 2017; 36:2577–88. https://doi.org/10.1038/onc.2016.414 [PubMed]
- 19. Song F, Xuan Z, Yang X, Ye X, Pan Z, Fang Q. Identification of key microRNAs and hub genes in non-small-cell lung cancer using integrative bioinformatics and functional analyses. J Cell Biochem. 2020; 121:2690–703. https://doi.org/10.1002/jcb.29489 [PubMed]
- 20. Liu Q, Chen J, Wang B, Zheng Y, Wan Y, Wang Y, Zhou L, Liu S, Li G, Yan Y. miR-145 modulates epithelial-mesenchymal transition and invasion by targeting ZEB2 in non-small cell lung cancer cell lines. J Cell Biochem. 2019; 120:8409–18. https://doi.org/10.1002/jcb.28126 [PubMed]
- 21. Korde A, Jin L, Zhang JG, Ramaswamy A, Hu B, Kolahian S, Guardela BJ, Herazo-Maya J, Siegfried JM, Stabile L, Pisani MA, Herbst RS, Kaminski N, et al. Lung endothelial MicroRNA-1 regulates tumor growth and angiogenesis. Am J Respir Crit Care Med. 2017; 196:1443–55. https://doi.org/10.1164/rccm.201610-2157OC [PubMed]
- 22. Liang Y, Rong X, Luo Y, Li P, Han Q, Wei L, Wang E. A novel long non-coding RNA LINC00355 promotes proliferation of lung adenocarcinoma cells by down-regulating miR-195 and up-regulating the expression of CCNE1. Cell Signal. 2020; 66: 109462. https://doi.org/10.1016/j.cellsig.2019.109462 [PubMed]
- 23. Lin K, Xu T, He BS, Pan YQ, Sun HL, Peng HX, Hu XX, Wang SK. MicroRNA expression profiles predict progression and clinical outcome in lung adenocarcinoma. Onco Targets Ther. 2016; 9:5679–92. https://doi.org/10.2147/OTT.S111241 [PubMed]
- 24. Cui Y, Zhang F, Zhu C, Geng L, Tian T, Liu H. Upregulated lncRNA SNHG1 contributes to progression of non-small cell lung cancer through inhibition of miR-101-3p and activation of Wnt/β-catenin signaling pathway. Oncotarget. 2017; 8:17785–94. https://doi.org/10.18632/oncotarget.14854 [PubMed]
- 25. Yin Q, Han Y, Zhu D, Li Z, Shan S, Jin W, Lu Q, Ren T. miR-145 and miR-497 suppress TGF-β-induced epithelial-mesenchymal transition of non-small cell lung cancer by targeting MTDH. Cancer Cell Int. 2018; 18:105. https://doi.org/10.1186/s12935-018-0601-4 [PubMed]
- 26. Zhong Z, Dong Z, Yang L, Chen X, Gong Z. MicroRNA-31-5p modulates cell cycle by targeting human mutL homolog 1 in human cancer cells. Tumour Biol. 2013; 34:1959–65. https://doi.org/10.1007/s13277-013-0741-z [PubMed]
- 27. Bao M, Pan S, Yang W, Chen S, Shan Y, Shi H. Serum miR-10a-5p and miR-196a-5p as non-invasive biomarkers in non-small cell lung cancer. Int J Clin Exp Pathol. 2018; 11:773–80. [PubMed]
- 28. Villar J, Zhang H, Slutsky AS. Lung repair and regeneration in ARDS: role of PECAM1 and Wnt signaling. Chest. 2019; 155:587–94. https://doi.org/10.1016/j.chest.2018.10.022 [PubMed]
- 29. Yu DH, Huang JY, Liu XP, Ruan XL, Chen C, Hu WD, Li S. Effects of hub genes on the clinicopathological and prognostic features of lung adenocarcinoma. Oncol Lett. 2020; 19:1203–14. https://doi.org/10.3892/ol.2019.11193 [PubMed]
- 30. Hao C, Cui Y, Owen S, Li W, Cheng S, Jiang WG. Human osteopontin: Potential clinical applications in cancer (Review). Int J Mol Med. 2017; 39:1327–337. https://doi.org/10.3892/ijmm.2017.2964 [PubMed]
- 31. Sun S, Wang Y, Wu Y, Gao Y, Li Q, Abdulrahman AA, Liu XF, Ji GQ, Gao J, Li L, Wan FP, Li YQ, Gao DS. Identification of COL1A1 as an invasion-related gene in Malignant astrocytoma. Int J Oncol. 2018; 53:2542–54. https://doi.org/10.3892/ijo.2018.4568 [PubMed]
- 32. Wei Z, Zhongqiu T, Lu S, Zhang F, Xie W, Wang Y. Gene coexpression analysis offers important modules and pathway of human lung adenocarcinomas. J Cell Physiol. 2020; 235:454–64. https://doi.org/10.1002/jcp.28985 [PubMed]
- 33. Franchini M, Frattini F, Crestani S, Bonfanti C, Lippi G. von Willebrand factor and cancer: a renewed interest. Thromb Res. 2013; 131:290–92. https://doi.org/10.1016/j.thromres.2013.01.015 [PubMed]
- 34. Bagnato A, Rosanò L. The endothelin axis in cancer. Int J Biochem Cell Biol. 2008; 40:1443–51. https://doi.org/10.1016/j.biocel.2008.01.022 [PubMed]
- 35. Berx G, van Roy F. Involvement of members of the cadherin superfamily in cancer. Cold Spring Harb Perspect Biol. 2009; 1:a003129. https://doi.org/10.1101/cshperspect.a003129 [PubMed]
- 36. Ha M, Son YR, Kim J, Park SM, Hong CM, Choi D, Kang W, Kim JH, Lee KJ, Park D, Han ME, Oh SO, Lee D, Kim YH. TEK is a novel prognostic marker for clear cell renal cell carcinoma. Eur Rev Med Pharmacol Sci. 2019; 23:1451–58. https://doi.org/10.26355/eurrev_201902_17102 [PubMed]