Genome-wide analysis of lncRNAs, miRNAs, and mRNAs forming a prognostic scoring system in esophageal squamous cell carcinoma

Data collection and pretreatment

The sequencing data (RNA-sequencing and miRNA-sequencing) and clinical information of ESCC patients were obtained from the TCGA database (https://portal.gdc.cancer.gov/). Based on the annotation file (Homo_sapiens.GRCh38.95.chr.gtf) downloaded from the Ensembl database (http://asia.ensembl.org/info/data/ftp/index.html), we identified 19876 lncRNAs and 19645 protein-coding genes. At the same time, we identified 2069 miRNA according to the annotation file (mature.fa) downloaded from miRBase database (http://mirbase.org/ftp.shtml). LncRNAs, mRNAs and miRNAs expressing raw count value >1 were screened for subsequent operation. This study was in line with the published guidelines provided by TCGA (https://cancergenome.nih.gov/publications/publicationguidelines). Since our data was obtained from the TCGA database, no ethics committee approval was required.

Differentially expressed analysis

The analysis and extraction of differentially expressed lncRNAs and mRNAs between 81 tumor tissues and 11 normal tissues were conducted by using the edgeR package of R language (Robinson, McCarthy & Smyth, 2010; R Core Team, 2013). Similarly, the differentially expressed miRNAs between 95 tumor tissues and 13 normal tissues were analyzed and extracted using edgeR package. |log2FC| > 2 and FDR < 0.05 (FC, fold change; FDR, false discovery rate) were considered to be significant. After edgeR normalization, log2 (normalized value +1) transformation was performed on the expression profiles of miRNAs, mRNAs and lncRNAs for subsequent manipulation.

Survival analysis

The clinical datasets of the ESCC cohort were downloaded from TCGA. Samples with a survival time of t = 0 days were removed to avoid introducing more mixed factors, and the remaining 80 samples were retained for the survival analysis. The clinical and pathological characteristics of the remaining 80 samples are summarized in Table 1. Univariate Cox regression analysis was used in R software to evaluate whether lncRNA, miRNA and mRNA were correlated with OS. RNAs with P < 0.05 were screened as prognostic biomarkers. RNAs with hazard ratio (HR) <1 were defined as protective signature, while RNAs with HR for death >1 were defined as risky RNAs.

LncRNA-mRNA co-expression network

The correlation between prognostic lncRNA and mRNA expression profiles was analyzed by Spearman method, and the lncRNA-mRNAs pairs that the absolute value of correlation coefficients > =0.4 and p < 0.05 were selected to construct the co-expression network. The co-expression network result was displayed using Cytoscape software version 3.6.0 (https://cytoscape.org/) (Shannon, 2003). CytoHubba, a plugin in the Cytoscape software, was adopted to calculate the degree of each node and select modules of hub genes from the network (Chin et al., 2014).

Functional enrichment analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional and pathway enrichment analysis were performed for mRNAs in prognosis-related co-expression RNA network using the Database for Annotation, Visualization and Integrated Discovery bioinformatics resources (DAVID) (https://david-d.ncifcrf.gov/), and P < 0.05 was considered as the cut-off criterion to screen the Enriched terms and pathways (Altermann & Klaenhammer, 2005; Ashburner et al., 2000; Huang, Sherman & Lempicki, 2009a; Huang, Sherman & Lempicki, 2009b).

Prognostic scoring system

RNAs with univariate Cox regression P < 0.01 were selected for the stepwise Cox regression procedures. Akaike information criterion (AIC) was used to evaluate the relative goodness of fitted model. Furthermore, Multivariate Cox regression coefficient was multiplied by the expression level of independent biomarkers (P < 0.05) to construct prognostic Cox models of lncRNA, miRNA and mRNA, respectively. Finally, a prognostic scoring system in ESCC was constructed, based on above-described multiple types of RNA. Receiver operating characteristics curves (ROC) and area under ROC curves (AUCs) were applied to evaluate the efficiency of each model. Statistical computing was performed using R software version 3.5.2. A flow diagram of the prognostic scoring system is presented in Fig. 1.

Statistical analysis

The statistical analyses in the present study were conducted by SPSS Statistics 18.0 and R software version 3.5.2. P value <0.05 was defined as statistically significance. Univariate Cox and multivariate Cox regression analyses were used to identify prognostic biomarkers. Survival curves were plotted by Kaplan–Meier (K-M) analysis, and differences in survival rates were assessed using a log-rank test.

Differentially expressed lncRNA, miRNA and mRNA

Analysis of expression profiles in ESCC compared with normal esophageal tissues identified a total of 1662 lncRNAs, 79 miRNAs and 2063 mRNAs (Table S1). Among them, 818 and 844 lncRNAs were respectively up-regulated and down-regulated (Fig. 2A); 52 miRNAs were up-regulated, and 27 were down-regulated (Fig. 2B); 869 up-regulated mRNAs and 1196 down-regulated mRNAs were obtained (Fig. 2C). Expression heatmaps were constructed by the top 50 up-regulation and the top 50 down-regulation to visualize the most significant lncRNAs, miRNAs and mRNAs (Fig. S1). The heatmap of the lncRNAs (Fig. S1A), miRNAs (Fig. S1B) and mRNAs (Fig. S1C) showed that the tumors clustered separately from the normal tissues.

Prognostic lncRNAs, miRNAs, mRNAs and co-expression network

Using univariate Cox regression analysis on the remaining 80 samples with survival times >0, 62 prognostic lncRNAs, eight prognostic miRNAs, and 66 prognostic mRNAs were identified in ESCC (P-value <0.05) (Table S2). The prognostic lncRNAs and mRNAs in ESCC were used to generate the co-expression network consisting of 22 lncRNAs, 40 mRNAs, and 77 interaction pairs (Fig. 3) (Table S3). Cytoscape analysis of the co-expression network revealed the top five prognostic RNAs (CDCA2, MTBP, CENPE, PBK, AL033384.1) (Table 2). Based on the median expression of each top 5 RNAs, 80 ESCC patients were divided into two groups (high expression vs low expression). The prognostic value of RNA was demonstrated by K-M plots (Fig. 4).

We then performed GO functional enrichment analysis of mRNAs in prognosis-related co-expression RNA network (Fig. 5). The results showed that the prognostic mRNAs mainly enriched in biological process (BP) including cell cycle, mitotic cell cycle and nuclear division. Cellular component (CC) analysis indicated enrichment in intracellular non-membrane-bounded organelle, non-membrane-bounded organelle and cytoskeletal part. Besides, in the molecular function (MF), the mRNAs were significantly clustered into purine nucleotide binding, ribonucleotide binding and ATP binding terms (Table S4). No pathways were significantly enriched in the KEGG enrichment analysis of prognostic mRNAs.

Prognostic scoring system

To create prognostic scoring system, RNAs with univariate Cox regression P < 0.01 were selected for the stepwise Cox regression procedures. Next, based on 5 lncRNAs, 2 miRNAs and 3 mRNAs respectively, we constructed three prediction models with single type of RNA to calculated the risk scores for predicted survival (Table S5). The formulas for the three prognostic models were as follows: lncRNA-based prognostic score = (0.447 × expression level of LINC01068) + (0.3677 × expression level of LINC00601) + (0.3075 × expression level of TTTY14) + (−0.8750 × expression level of AC084262.1) + (−0.4744 × expression level of LINC01415); miRNA-based prognostic score = (1.2932 × expression level of miR-5699-3p) + (0.7202 × expression level of miR-552-5p); mRNA-based prognostic score = (0.5139 × expression level of MLIP) + (0.5746 × expression level of TNFSF10) + (−1.0069 × expression level of SIK2). Of three prognostic models, seven RNAs were shown to be risky RNAs (LINC01068, LINC00601, TTTY14, miR-5699-3p, miR-552-5p, MLIP, TNFSF10, HR >1) and three RNAs were the protective RNAs (AC084262.1, LINC01415, SIK2, HR <1) (Figs. 6A–6C).

Using these three formulas, we calculated the prognostic score for each of the 80 patients separately and ranked them according to the increased prognostic scores. we divided the ESCC patients into two group (high-risk or low-risk) using the median prognostic score as a cutoff. As shown in Figs. 6D–6F, patients in the high-risk group had a worse prognosis than the low-risk group in all three models (P < 0.0001). We also used ROC curves to estimate the specificity and sensitivity of these prognostic models. All three prognostic models showed moderate prognostic evaluation ability, with AUC of 1 year values of 0.855, 0.859 and 0.785, separately, and AUC of 3 year values of 0.909, 0.709 and 0.762, separately (Fig. 7). Figure 8 shows the distribution of patient prognostic scores, the survival status and tumor RNAs expression of all 80 ESCC patients. Patients in the high-risk group had more deaths than those in the low-risk group in all three models (Figs. 8D–8F). Moreover, patients in the low-risk group tend to express protective RNA, while patients in the high-risk group tend to express risky RNA (Figs. 8G–8I).

In order to improve the prediction accuracy and understand the potential molecular mechanism of prognostic markers, we constructed a prognostic Cox model with multiple types of RNA for ESCC using the ten RNAs provided above (Table S6). The formula was as follows: RNA-based prognostic score = (0.42895 × expression level of LINC01068) + (0.34829 × expression level of LINC00601) + (0.2185 × expression level of TTTY14) + (−1.393 × expression level of AC084262.1) + (−0.33364 × expression level of LINC01415) + (1.06024 × expression level of miR-5699-3p) + (0.34784 × expression level of miR-552-5p) + (0.3418 × expression level of MLIP) + (0.05437 × expression level of TNFSF10) + (−1.38365 × expression level of SIK2). The ten RNAs forest plots of RNA-based prognostic model was presented in Fig. 9. K-M analysis showed patients in the low-risk group had better long-term survival than those in the high-risk group (P-value <0.0001; Fig. 10A). Furthermore, the AUC of 1 year value was 0.916 and 3 year value was 0.917 (Fig. 10B), indicating that the combination of different types of RNA patterns is a more accurate prognostic model than single type of RNA prediction model. For tumor staging, we also generated K-M plots and corresponding ROC curves. ESCC patients were also divided into two groups by tumor, node and metastasis (TNM) stage, and the prognosis of the two groups was different (P-value <0.05; Fig. 10C). The AUC of 1 year value and 3 year value based on TNM staging were 0.612 and 0.548, respectively (Fig. 10D). Although TNM staging is often used in clinical prognostic prediction, its prognostic AUC value is limited. Besides, combining multiple clinical parameters, we performed cox regression analysis of the prognostic score. As shown in Table 3, in both univariate Cox and multivariate Cox regression analysis, prognostic score was significantly correlated with survival (P < 0.001), that is, prognostic score was an independent prognostic factor in ESCC patients.