Four-miRNA signature as a prognostic tool for lung adenocarcinoma

Introduction

Lung cancer is a highly prevalent disease and the most frequently diagnosed cancer. It is also highly fatal, being the most common cause of cancer-related deaths in males and the second most common cause in females.¹ Lung adenocarcinoma (LUAD) is the most often diagnosed histological subtype of non-small-cell lung cancer (NSCLC).² Despite the advent of targeted therapies that have improved outcomes for a subset of patients,^3–6 the overall 5-year survival rate for lung cancer patients remains low at 17.4%.⁷ For patients with NSCLC, the best prognosis is achieved using a complete surgical resection of stage IA disease, yielding a 5-year survival rate up to 70%. However, this requires early and accurate disease prognosis. Lung cancer can be detected early by a computed tomography scan in high-risk individuals, but this method has a high false-positive rate that can lead to unnecessary treatment.⁸ Therefore, more prognostic indicators are urgently needed.

miRNAs are highly conserved, endogenous noncoding RNAs (20–22 bp in length) that function in a variety of biological processes at the transcriptional or posttranscriptional level. Abnormal expression of miRNAs is a common problem in cancer development. Specifically, there is growing evidence for a role of various miRNAs in the progression and biology of NSCLC.⁹ Several miRNAs have been identified as markers in cancers, including a four-miRNA signature in colon cancer¹⁰ and a three-miRNA signature in cervical cancer.¹¹ Therefore, miRNAs are promising indicators for cancer diagnosis, prognosis, and targeted therapies.

The Cancer Genome Atlas Project (TCGA) is a National Cancer Institute effort to profile 29 different tumor types using genomic platforms and to make raw and processed data readily available to all researchers.¹² The TCGA has released massive datasets of miRNA sequencing data from patients with LUAD. For this study, we downloaded high-throughput miRNA data from the TCGA database to identify several miRNAs with differential expression between LUAD and paired healthy lung tissues. We then evaluated the prognostic value of these differentially expressed miRNAs and found a four-miRNA signature that could effectively predict patient survival. Furthermore, we analyzed the pathways and functions of the target genes modulated by these four miRNAs. Ultimately, we hope that this signature can be used as a novel prognostic tool for LUAD, and that our functional analysis may provide more insights into the molecular mechanisms of this prevalent and devastating disease.

Patients and methods

Data processing

The preprocessed mature miRNA expression profiles, displayed as reads per million, and clinical information from the TCGA (TCGA data version 2016_01_28 for LUAD) were downloaded from Broad Institute GDAC Firebrowse (http://firebrowse.org/). The inclusion criteria were as follows: 1) the sample included both miRNA expression profiles and clinical information, and 2) the sample included prognosis information.

Generation of a prognostic signature for LUAD using differentially expressed miRNAs

To create a LUAD prognostic signature, we first removed miRNAs with expression levels of zero in >50% of the patients. We then normalized the miRNA expression profiles by log2 transformation using the R language project. We identified the miRNAs that were differentially expressed between LUAD and paired healthy tissues using the “limma” package in R. Finally, we calculated fold changes in the expression of individual miRNAs between LUAD and healthy tissues and considered miRNAs with a |log2FC| >1 and a P-value <0.05 to be significant.

The four-miRNA signature index and survival

In our univariate Cox proportional hazards regression analysis, we selected the median expression level for each miRNA as the cutoff point to dichotomize patients into high- and low-expression groups. To determine if the expression level was associated with a prognostic outcome, we applied a Kaplan–Meier survival analysis with the log-rank method. We then confirmed the prognostic value of the four statistically significant miRNAs using a multivariate Cox proportional hazards regression analysis. For each of these four miRNAs, we assigned each patient one point for the high-risk expression level. Thus, each patient received a score ranging from 0 to 4, which we called the signature index. We plotted a survival curve based on the patients’ survival status, survival days, and signature index using the Kaplan–Meier survival method (log-rank). We considered patients with a signature index of 3 or 4 to be high risk and those <3 to be low risk.

Prediction of the target genes of the prognostic miRNA signature

We predicted the target genes of the four prognostic miRNAs using the TargetScan (http://www.targetscan.org/), miRDB (http://www.mirdb.org/miRDB/), and miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/index.php) online analysis tools. To further enhance the reliability of the bioinformatic analysis, we created a Venn diagram for each miRNA to identify the overlapping results among the three prediction tools. We performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of the overlapping target genes using the “clusterProfiler” package in R.¹³ The clusterProfiler package implements methods to analyze and visualize the functional profiles of genes and gene clusters. We set the cutoff criteria as a P-value <0.05 and a gene count ≥2.

Statistical analysis

We analyzed the expression levels of the miRNAs in LUAD and matched healthy tissues using unpaired t-tests. We performed Kaplan–Meier survival analysis and univariate/multivariate Cox proportional hazards regression analysis to compare the expression levels (low vs high) and prognostic significance (low risk vs high risk) of each miRNA. We considered P-values <0.05 to be statistically significant. Statistical analysis was performed using IBM SPSS Statistics software program version 22.0 (IBM Corp., Armonk, NY, USA).

Results

Differentially expressed miRNAs in LUAD and healthy tissues

Our analysis included a total of 438 LUAD tissues and 45 matched healthy tissues from the same patients. The detailed clinical characteristics including age at diagnosis, metastasis, lymph node status, pathological stage, and T stage are listed in Table 1. According to the cutoff criteria (P<0.05 and |log2FC| >1), 187 miRNAs were differentially expressed between LUAD and matched healthy tissues. This included 148 miRNAs that were upregulated in the LUAD tissues and 39 miRNAs that were downregulated. The results of the expression analysis are presented as a volcano plot (Figure 1) to demonstrate that the distribution of P-values and |log2FC| was reasonable with respect to each other.

Table 1 Summary of patient cohort information

Figure 1 Volcano plot of the differentially expressed miRNAs between LUAD and paired healthy tissues. Notes: Red indicates high expression in LUAD tissues, while blue indicates low expression. Gray shows the miRNAs with a |log2FC| <1 and/or an adjusted P-value ≥0.05. The volcano plot was constructed using the ggplot2 package of the R language. Abbreviation: LUAD, lung adenocarcinoma.

Identification of four miRNAs associated with survival in LUAD

For each of the 187 differentially expressed miRNAs, we used the median expression level as a cutoff to stratify the 438 patients into two groups: a high-level group and a low-level group. The univariate Cox proportional hazards regression analysis revealed that a total of 22 miRNAs had prognostic value. We subsequently applied a multivariate Cox proportional hazards regression analysis to verify these results and identified four miRNAs, miR-148a-5p (Figure 2A), miR-31-5p (Figure 2B), miR-548v (Figure 2C), and miR-550a-5p (Figure 2D), as independent prognostic indicators. The results of the univariate and multivariate analyses for the 22 miRNAs are displayed in Table 2. All four miRNAs were upregulated in the LUAD tissues. Patients with high expression levels of miR-148a-5p, miR-31-5p, and miR-1548v had longer survival times compared with patients with low expression levels of those miRNAs. In contrast, patients with high expression levels of miR-550a-5p had shorter survival times than those with low expression levels.

Figure 2 Four miRNAs were associated with overall survival in LUAD patients using Kaplan–Meier curves and log-rank tests. Notes: The patients were stratified into high- and low-risk groups according to the median expression of each miRNA. Each plot represents the Kaplan–Meir curves obtained using a different miRNA: (A) miR-148a-5p, (B) miR-31-5p, (C) miR-548v, and (D) miR-550a-5p. Abbreviation: LUAD, lung adenocarcinoma.

Table 2 Univariate and multivariate analyses in LUAD patients Note: *P value of <0.05 was considered statistically significant difference. Abbreviations: LUAD, lung adenocarcinoma; HR, hazard ratio; CI, confidence interval.

A four-miRNA signature index for LUAD prognosis

We scored the four-miRNA signature by value assignment for each miRNA (Table 3). Each patient was assigned one score for a low expression level of miR-148a-5p, miR-31-5p, and miR-548v, respectively, and one score for a high expression level of miR-550a-5p. Then, we summed the scores for each patient, resulting in a signature index ranging from 0 to 4. All patients were divided into high- or low-risk groups according to this index. Specifically, we considered patients with a signature index of 3 or 4 to be high risk and those with an index <3 to be low risk. We compared the two groups in a Kaplan–Meier survival analysis using the log-rank method. This revealed that patients in the high-risk group (n=123) had significantly shorter survival times than those in the low-risk group (n=315) (P=0.002; Figure 3). This indicates that our four-miRNA signature created was associated with survival in LUAD.

Table 3 Patient scoring method using the four-miRNA signature index Note: aA score of 3 or 4 indicates the patients to be high risk, while a score <3 indicates low risk. Abbreviation: RPM, reads per million.

Figure 3 The Kaplan–Meier curves obtained using the four-miRNA signature to separate patients into high- and low-risk groups.

Target prediction and functional analysis

To explore the potential biological functions of the four miRNAs (miR-148a-5p, miR-31-5p, miR-548v, and miR-550a-5p), we predicted their respective target genes using the TargetScan, miRDB, and miRTarBase online analysis tools. Each analysis tool identified unique target genes for the four miRNAs, but all three tools identified 12 overlapping target genes for miR-148a-5p, 17 for miR-31-5p, nine for miR-548v, and 6 for miR-550a-5p (Figure 4A–D). We analyzed the 44 potential target genes using the KEGG signal pathway and GO databases. Those 44 potential target genes are shown in Table 4. The GO annotation results have three parts: a cellular component, a biological process, and a molecular function. The results of the significantly enriched GO term analysis are shown in Figure 5A–C (P<0.05, gene count ≥2). This analysis revealed that the molecular functions primarily associated with these genes were protein serine/threonine kinase activity and protein-phosphorylated amino acid binding. The KEGG pathways identified in the analysis were significantly enriched with the Hippo and sphingolipid signaling pathways (P<0.05, gene count ≥2) (Figure 5D). Together, these results suggest that the differentially expressed miRNAs may modulate protein phosphorylation, particularly in the identified signaling pathways.

Figure 4 The overlapping target genes of the respective miRNAs were predicted using the miRDB, TargetScan, and miRTarBase online analysis tools: (A) miR-148a-5p, (B) miR-31-5p, (C) miR-548v, and (D) miR-550a-5p.

Table 4 The overlapping potential target genes of the four miRNAs

Figure 5 The significantly enriched Gene Ontology annotation of the target genes of the four miRNAs including their (A) cellular component, (B) biological process, and (C) molecular function. (D) Kyoto Encyclopedia of Genes and Genomes pathway analysis of the predicted targets of the four miRNAs.

Discussion

LUAD is a devastating cancer subtype with a low survival rate. Despite its prevalence, prognostic tools for LUAD are limited. The prognosis of patients with LUAD would be improved considerably if tumor behavior could be reliably predicted at the time of the initial diagnosis. This requires the identification of novel biomarkers and an understanding of the molecular mechanisms underlying LUAD development. In this paper, we identified a total of 187 miRNAs that were differentially expressed between LUAD and matched healthy tissues, 4 of which were associated with overall survival in patients with LUAD. Because miRNAs are known to modulate gene expression, we screened the target genes of those four miRNAs and used bioinformatics to predict the pathways and biological functions associated with their targets.

Several miRNA markers have been identified for the prediction and prognosis of cancers such as head and neck squamous cell carcinomas,¹⁴ bile duct cancer,¹⁵ and glioma.¹⁶ In this paper, we sought to identify miRNA markers associated with prognosis for LUAD. We created a four-miRNA signature index that successfully separated patients into low- and high-risk groups. Specifically, patients deemed high risk by our four-miRNA signature index had significantly shorter survival times than those in the low-risk group (P=0.002). We additionally performed a receiver operating characteristic (ROC) curve study after multivariate analysis, but were unable to achieve the suggested area under the curve (AUC) >0.7 (data not shown). However, this discrepancy is similar to the results of previous studies. Several studies have investigated the relationship between miRNA expression patterns and the survival and prognosis of patients with LUAD. Li et al reported that miR-101-1, miR-220a, miR-4661, and miR-450a-2 were associated with overall survival in patients with LUAD.¹⁷ This prognosis model had significant efficacy for Caucasian patients, but as in the current study, their model did not reach an AUC of 0.7 (AUC =0.629). This study also did not show their sensitivity and specificity, so it is difficult to fully evaluate the efficacy of their model. Alternatively, Peng et al performed a comprehensive analysis based on LUAD miRNome profiling studies using the robust rank aggregation method. They identified a panel of six miRNAs (miR-21-5p, miR-210-3p, miR-182-5p, miR-183-5p, miR-126-3p, and miR-218-5p) associated with LUAD and evaluated the expression and prognostic values of those miRNAs in an independent cohort.¹⁸ They did not perform a ROC curve study. Sathipati and Ho proposed an optimized support vector regression method to identify an miRNA signature for estimating LUAD patients survival using their miRNA expression profiles. Their miRNA signature consists of 18 miRNAs associated with LUAD: hsa-let-7f-1, hsa-miR-16-1, hsa-miR-152, hsa-miR-217, hsa-miR-18a, has-miR-193b, hsa-miR-3136, hsa-let-7g, hsa-miR-155, hsa-miR-3199-1, hsa-miR-219-2, hsa-miR-1254, hsa-miR-1291, hsa-miR-192, hsa-miR-3653, hsa-miR-3934, hsa-miR-342, and hsa-miR-141.¹⁹ In our study, we determined that increased expression levels of four mature miRNAs (miR-148a-5p, miR-31-5p, miR-548v, and miR-550a-5p) were associated with clinical outcomes in patients with LUAD. Interestingly, the miRNA signatures generated by these various studies do not overlap, perhaps because of variability in the methods used to identify miRNAs. And, no study can build a predictive tool to achieve the desired AUC. Therefore, predicting the prognosis of LUAD with miRNA signature still needs more in-depth research. However, our study and the previous studies all emphasize new and constructive strategies to identify prognostic markers for LUAD. One strength of our study is that we removed miRNAs with expression levels of zero in >50% of the patients, which would enhance the feasibility of implementation.

To gain insight into the molecular functions of the four miRNAs, we predicted the target genes for the four miRNAs and identified their associated KEGG pathways and GO annotations. The molecular functions of the target genes were mainly associated with protein serine/threonine kinase activity and protein-phosphorylated amino acid binding. Thus, the four miRNAs might be closely related to protein phosphorylation. The target genes were also significantly enriched for the Hippo signaling pathway (six genes, P=7.87×10⁻⁶). Abnormal regulation of the Hippo signaling pathway is involved in the development of a variety of human cancers such as breast cancer, hepatocellular carcinoma, hematological cancer, and lung cancer.²⁰ Thus, abnormal regulation of signaling pathways, including Hippo, may play a crucial role in the pathogenesis and progression of LUAD. The interaction between the four miRNAs and the Hippo pathway merits further exploration, as more research is necessary to determine the exact molecular mechanisms of the role of miRNAs in LUAD.

Our study has a few limitations. Although we included a relatively large sample size of 438 patients, we used the median expression level as the cutoff value to stratify patients. Despite the ability of our model to successfully stratify patients into high- and low-risk groups, it did not yield the suggested AUC. However, this was similar to the results of previous studies examining the relationship between miRNA expression and LUAD prognosis. Thus, the prognostic value of the four-miRNA signature needs further validation with other methods and a larger sample size. Further functional investigation is also required to explore the molecular functions of the four miRNAs in LUAD progression.

Conclusion

Taken together, our results identified a four-miRNA signature that could be a prognostic tool for patients with LUAD. The four miRNAs were expressed in >50% of patients with LUAD and were more highly expressed in LUAD tissues than in matched healthy ones. These miRNAs modulated genes associated with protein phosphorylation and the Hippo signaling pathway, which has previously been associated with cancer progression. Ultimately, we hope that this miRNA signature can help to predict LUAD prognosis and uncover the mechanisms underlying LUAD development.

Acknowledgments

This study was supported by the National Natural Science Foundation of China (Grant No 81360290), the Guangxi Natural Science Foundation (Grant Nos 2017GXNSFAA198249, 2016GXNSFBA380090 and 2015GXNSFAA139128), the Key Research and Development project of Guangxi (Grant No Guike AB17195002), the Basic Ability Enhancement Program for Young and Middle-age Teachers of Guangxi (Grant No 2017KY0120), the Self-raised Scientific Research Funds of Ministry of Health of Guangxi Province (Grant No Z2016480) and the China Scholarship Council (201608455001).

Disclosure

The authors report no conflicts of interest in this work.

---

References

1.		Torre LA, Bray F, Siegel RL, Ferlay J, Lortet-tieulent J, Jemal A. Global cancer statistics, 2012. CA Cancer J Clin. 2015;65(2):87–108.
2.		Travis WD. Pathology of lung cancer. Clin Chest Med. 2011;32(4):669–692.
3.		Zhou C, Wu Y-L, Chen G, et al. Erlotinib versus chemotherapy as first-line treatment for patients with advanced EGFR mutation-positive non-small-cell lung cancer (OPTIMAL, CTONG-0802): a multicentre, open-label, randomised, phase 3 study. Lancet Oncol. 2011;12(8):735–742.
4.		Maemondo M, Inoue A, Kobayashi K, et al. Gefitinib or chemotherapy for non-small-cell lung cancer with mutated EGFR. N Engl J Med. 2010;362(25):2380–2388.
5.		Paez JG, Jänne PA, Lee JC, et al. EGFR mutations in lung cancer: correlation with clinical response to gefitinib therapy. Science. 2004;304(5676):1497–1500.
6.		Kwak EL, Bang Y-J, Camidge DR, et al. Anaplastic lymphoma kinase inhibition in non-small-cell lung cancer. N Engl J Med. 2010;363(18):1693–1703.
7.		Siegel RL, Miller KD, Jemal A. Cancer statistics, 2016. CA Cancer J Clin. 2016;66(1):7–30.
8.		Aberle DR, Abtin F, Brown K. Computed tomography screening for lung cancer: has it finally arrived? Implications of the National Lung Screening Trial. J Clin Oncol. 2013;31(8):1002–1008.
9.		Boeri M, Pastorino U, Sozzi G. Role of microRNAs in lung cancer. Cancer J. 2012;18(3):268–274.
10.		Xu J, Zhao J, Zhang R. Four microRNAs signature for survival prognosis in colon cancer using TCGA data. Sci Rep. 2016;6:38306.
11.		Liang B, Li Y, Wang T. A three miRNAs signature predicts survival in cervical cancer using bioinformatics analysis. Sci Rep. 2017;7(1):1–8.
12.		Chandran UR, Medvedeva OP, Barmada MM, et al. TCGA expedition: a data acquisition and management system for TCGA data. PLoS One. 2016;11(10):e0165395.
13.		Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–287.
14.		Wong N, Khwaja SS, Baker CM, et al. Prognostic microRNA signatures derived from The Cancer Genome Atlas for head and neck squamous cell carcinomas. Cancer Med. 2016;5(7):1619–1628.
15.		Wang M, Wen T-F, He L-H, Li C, Zhu W-J, Trishul NM. A six-microRNA set as prognostic indicators for bile duct cancer. Int J Clin Exp Med. 2015;8(10):17261–17270.
16.		Yan W, Li R, Liu Y, et al. MicroRNA expression patterns in the malignant progression of gliomas and a 5-microRNA signature for prognosis. Oncotarget. 2014;5(24):12908–12915.
17.		Li X, An Z, Li P, Liu H. A prognostic model for lung adenocarcinoma patient survival with a focus on four miRNAs. Oncol Lett. 2017;14(3):2991–2995.
18.		Peng Z, Pan L, Niu Z, et al. Identification of microRNAs as potential biomarkers for lung adenocarcinoma using integrating genomics analysis. Oncotarget. 2017;8(38):64143–64156.
19.		Sathipati SY, Ho S-Y. Identifying the miRNA signature associated with survival time in patients with lung adenocarcinoma using miRNA expression profiles. Sci Rep. 2017;7(1):7507.
20.		Yu F-X, Zhao B, Guan K-L. Hippo pathway in organ size control, tissue homeostasis, and cancer. Cell. 2015;163(4):811–828.