A six-microRNA signature can better predict overall survival of patients with esophagus adenocarcinoma

RNA-seq and clinicopathological data of EAC patients

From TCGA data portal (https://portal.gdc.cancer.gov/), RNA-seq data and associated clinical information were downloaded in January 2019. The annotation information was provided by GENCODE datasets (www.gencodegenes.org). Given that some miRNAs and mRNAs display little or no expression in some tissues or do not vary sufficiently, only those with raw count value >20 in more than 80% of samples were retained for further analysis. Once normalization by edgeR was completed, this was followed by conversion of the expression patterns of miRNAs and mRNAs to log2 (normalized value +1) in preparation for the subsequent processing. Samples with a <1-month censor time are removed, because they cannot be representative samples for analyzing prognostic factors. A total of 84 EAC subjects with the corresponding clinical data including age, gender, height, weight, race, alcohol history, Barrett’s disease history, tumor size, lymph node status, metastasis status, and TNM stage were collected in this study (Table 1). The EAC patients’ dataset contained 96 samples (84 EAC and 12 normal tissues) and 272 miRNAs. Since the data came from the TCGA database, no further approval was required from the Ethics Committee.

Construction and validation of the miRNA risk score

Eighty-four patients were stratified to two categories in a random manner: training set = 42, test set = 42. Training set was analyzed to build a miRNA model that was later confirmed in test and entire sets. In the training set, we screen out miRNAs with a significant p-value less than 0.1 by using univariate survival analysis based on Cox proportional hazards of each miRNA. The least absolute shrinkage and selection operator (LASSO) is a generalized linear regression algorithm capable of variable selection and regularization simultaneously (Gao, Kwan & Shi, 2010). We determined the lambda by using the cross-validation routine cv.glmnet with an n-fold equal to 10. Least absolute shrinkage and selection operator was performed to reduce above selected prognostic miRNAs further and to construct the risk score system.

For determination of survival risks, a prognostic model was created on the basis of miRNA data as follows: $$\text{Risk score}={\displaystyle \sum _{i=1}^n\beta i\ast genei}$$β stands for the coefficient of the miRNA, and gene refers to miRNA expression value.

Using the median score in training set as the cutoff, we stratified the subjects to low-risk and high-risk groups. The Kaplan–Meier (KM) and log-rank methods were applied to compare the survival rate between the groups by using the R “survival” package. The time-dependent receiver operating characteristic (ROC) curve was plotted by using the R “timeROC” package to evaluate specificity and sensitivity of the miRNA expression-based prognostic signature. Thereafter, this signature was validated in test set and entire set. ROC and KM curves were also carried out to validate accuracy and feasibility of the miRNA model. Then stratified analysis based on clinical parameters was performed in the entire set. All ROC and KM curves were plotted with R (version 3.5.2), and p < 0.05 represented statistical significance.

Gene set enrichment analysis

Subjects were stratified to two groups (high and low) based on the risk score of the 6-miRNA signature. We used gene set enrichment analysis (GSEA, http://software.broadinstitute.org/gsea) (Subramanian et al., 2005) to figure out potential functional annotations in the two groups. The BioCarta dataset (c2.cp.biocarta.v6.2.symbols.gmt) served as the reference gene set. False discovery rate < 0.05, enrichment score > 0.5 were set as the significance threshold.

Functional enrichment analysis

Using the miRNA target prediction tool starBase (http://starbase.sysu.edu.cn/index.php), the target genes of the 6-miRNA signature were predicted based on five datasets, including TargetScan, PITA, miRmap, microT, and miRanda. Metascape is a free online platform having a large-scale set of functional annotation tools to understand biological mechanisms behind a large pool of genes (http://metascape.org/gp/index.html#/main/step1). We used Metascape to analyze functional enrichment of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and Gene Ontology (GO) based on the prognostic target genes of miRNAs and visualized by R “gglot2” package.

The predictive 6-miRNA signature for the training set

The overall design and workflow of this study were presented in Fig. 1. According to the results of the univariate Cox regression analyses, 64 miRNAs associated with survival data were selected for patients with EAC (Table S1). The lambda value was set by using the lambda.min, which is the value of lambda giving minimum mean cross-validated error, and then six miRNAs with nonzero coefficients were defined (Fig. S1). Based on the LASSO Cox regression models, a risk score was determined for each subject according to 6-miRNA status: Risk score = (−0.6089 × hsa-let-7b) + (−0.1974 × hsa-mir-23a) + (0.3369 × hsa-mir-3074) + (0.0294 × hsa-mir-424) + (0.2421 × hsa-mir-425) + (0.2435 × hsa-mir-505).

In the training set, the patients with EAC were divided into a high-risk group and a low-risk group. The risk scores of patients were ranked, and the dotplot was developed for the survival status of each patient. Compared with the mortality of patients in the high-risk group, those in the low-risk group was much lower (Figs. 2A and 2C). Moreover, based on a heatmap of the 6-miRNA profile, the levels of hsa-mir-3074, hsa-mir-424, hsa-mir-425, and hsa-mir-505 were lower in the low-risk group than those of the high-risk group. The level of hsa-let-7b and hsa-mir-23a were higher in the low-risk group than those of the high-risk group (Fig. 2D). The KM curve indicated that the survival time of patients in the high-risk group was shorter than those in the low-risk group (Fig. 2E). We described the predictive value of the 6-miRNA signature by using a time-dependent ROC curve. The area under curve (AUC) at 1, 2, and 3 years of the signature was 0.860, 0.962, 0.959, respectively (Fig. 2B).

The predictive power of 6-miRNA signature in test set and entire set

The 6-miRNA signature was applied to the test set and the entire set for evaluation of its prognostic value. The distribution of risk scores, the expression values of six miRNAs and the survival status of patients ranked according to the risk scores were presented in test set (Figs. 3A, 3C and 3E) and entire set (Figs. 3B, 3D and 3F). In test set and entire set, patients with the low-risk scores exhibited better overall survival than those with the high-risk scores based on the KM curve (Figs. 3G and 3H). The 3-year AUC of the 6-miRNA-based signature was 0.840 and 0.868, respectively, for the test set and the entire set (Figs. 3I and 3J).

To assess the independent prognostic value of the 6-miRNA signature, various clinicopathological factors were subjected to univariate Cox regression and multivariate Cox regression analysis. The result indicated that the 6-miRNA signature was an independent prognostic factor after adjustment for other clinicopathological factors (HR = 2.95, CI [1.43–6.07], p = 0.00338, Table 2). When stratified by clinical factors (age, gender, caucasian, height, weight, alcohol consumption history, Barrett’s disease, TNM stage), a nearly universal result was obtained for all subgroups (Fig. 4), showing that high-risk score was strongly associated with poor prognosis and vice versa. Regardless of height, weight, TNM stage, alcohol consumption history and Barrett’s disease, the 6-miRNA signature is significantly effective. Therefore, the present results suggest that the 6-miRNA signature can predict the clinical prognosis of EAC.

Functional analysis of the 6-miRNA signature

BioCarta pathway enrichment was conducted through GSEA in high-risk group in the entire set. It revealed that high-risk patients were associated with some pathways, including “proteasome pathway,” “MCM pathway,” “G2 pathway,” and “cell cycle pathway” (Figs. 5A–5D). Through a miRNA prediction tool, starBase, 179 target mRNAs for hsa-let-7b, 147 for hsa-mir-23a, 382 for hsa-mir-424, 37 for hsa-mir-425, and 11 for hsa-mir-505 were obtained. Unfortunately, no target gene for hsa-mir-3074 was predicted. We conducted functional enrichment of these target genes by GO and KEGG categories. Cellular component, molecular function and biological process of these target genes based on p-values were showed (Figs. 5E–5G). The top 20 KEGG pathways of these target genes were plotted (Fig. 5H). Among these pathways, MAPK signaling pathway, hippo signaling pathway, foxo signaling pathway, and TGF-beta signaling pathway were reported to be related to metastasis of cancer (Blum et al., 2019; Janse van Rensburg & Yang, 2016; Kim et al., 2018; Sun et al., 2018). Some other pathways are also known to be associated with cancers, such as pathways in cancer, miRNAs in cancer, cell cycle, autophagy.