Screening of prognostic biomarkers for endometrial carcinoma based on a ceRNA network

Data sources and pretreatment

EC RNA-seq and miRNA-seq data were downloaded from the TCGA database (https://cancergenome.nih.gov/) (Tomczak, Czerwińska & Wiznerowicz, 2015). RNA-seq data included 551 EC and 35 paracancerous samples. miRNA-seq data included 546 EC and 33 paracancerous samples. As the data was obtained from the TCGA website, further approval by an ethics committee was not required. This study meets the publication guidelines provided by TCGA (http://cancergenome.nih.gov/publications/publicationguidelines). RNA-seq and miRNA-seq matrix files describing samples were downloaded from TCGA. After that, mRNA and lncRNA expression profiles were extracted from the RNA-seq matrix information. Thus, we obtained three matrix files: mRNA, lncRNA and miRNA expression profiles.

Screening of differentially expressed genes

Statistical analysis was performed using the exactTest function in the edge.R software package (Robinson, McCarthy & Smyth, 2010). and differentially expressed genes were screened out. Next, the P value was FDR-corrected. mRNA, lncRNA and miRNA with FDR < 0.05 were screened out. RNAs with fold change >4 (|logFC| > 2) were further selected as differentially expressed mRNA, lncRNA and miRNAs. A volcano plot of differentially expressed genes was generated.

ceRNA topological network construction

Topological and stability analysis of the ceRNA network

The NetWorkAnalyzer toolkit (Shannon et al., 2003) from Cytoscape (Doncheva et al., 2012) was used to analyze the characteristics of the ceRNA topological network. NetWorkAnalyzer is mainly used to analyze network diameter, connection numbers and average clustering coefficient. In this study, we calculated network connections, the path length and the closest centrality of nodes.

Functional and enrichment analysis of ceRNA

Biological function characterization of DEmRNAs in the ceRNA network was performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) version 6.7 (https://david-d.ncifcrf.gov/home.jsp), and pathway enrichment analysis was performed using the KO-Based Annotation System (KOBAS) 3.0 (http://kobas.cbi.pku.edu.cn/) (Huang da, Sherman & Lempicki, 2009; Xie et al., 2011).

Prognostic analysis of the ceRNA module

Clinical information of samples was downloaded from the TCGA database, and survival data were extracted. Combined with expression profile data, a Kaplan–Meier (K–M) survival curve was generated for each node in the ceRNA topological network using survival (Therneau, 2015) in the R package. Survival differences for genes, lncRNAs and miRNAs in the ceRNA module were also analyzed. Patients were dichotomized for survival analysis using a log-rank test with optimal cutoff values determined by the “surv_cutpoint” function of the “survminer” R package. P < 0.05 was considered statistically significant.

Validation of DEmRNAs and DEmiRNAs in the ceRNA network

The GSE17025 (Allard et al., 2008) and GSE35794 datasets were used to validate the expression of candidate mRNAs and miRNAs, respectively. Transcriptome profiling data in the GSE17025 (Platform: GPL570) dataset contained 91 EC and 12 non-tumor samples. The microRNA expression profile in the GSE35794 (Platform: GPL10850) dataset contained 18 EEC and four normal samples.

Screening of differentially expressed RNAs

The mRNA expression profiles of 35 normal samples and 551 EC samples were compared. After statistical validation, 2632 differentially expressed mRNAs were screened out including 1672 upregulated and 960 downregulated mRNAs (Figs. 1A and 1D; Table S1). We obtained 1178 differentially expressed lncRNAs after screening, including 867 that were upregulated and 311 that were downregulated. Results are shown in Figs. 1B and 1E; Table S1.

The miRNA expression profiles of 33 normal samples and 546 EC samples were compared, and 189 differentially expressed miRNAs were obtained after statistical validation. One hundred and forty of these 189 miRNAs were upregulated and 49 were downregulated. Results are shown in Figs. 1C and 1F; Table S1.

Construction of the ceRNA network and analysis of topology and stability

Based on the differentially expressed lncRNA and miRNA profiles, lncRNA and corresponding miRNAs were paired using the miRcode database. We established 556 pairs, with 99 lncRNAs and 27 miRNAs (Table S2). Twenty-five of these 27 DEmiRNAs were found in the starBase database. The target genes of these 25 miRNAs were predicted using three public databases: miRDB, miRTarBase and TargetScan. Nine hundred and fifty-six target genes were found in total. These target genes were compared with 2632 differentially expressed genes that were screened out, and genes in common were selected out. Seventy-three DEmRNAs that could interact with 20 of the 25 DEmiRNAs in all three datasets were selected (Table S3). After removing the remaining seven DEmiRNAs and corresponding lncRNAs, 97 DElncRNAs, 20 DEmiRNAs and 73 DEmRNAs were ultimately used to construct a ceRNA network used the “ggalluvial” R package for visualization (Fig. 2A; Table S4).

The topology of the ceRNA network was analyzed to verify the network’s reliability. The distribution of the numbers of nodes is shown in Fig. 2B. Node number decreased with increases in distribution degree. The highest degree was 51; the lowest was zero. Most nodes in the network were nodes with fewer interactions. The closeness centrality (CC) of a given node can be used calculate the connection steps from the node to another. More centralized nodes show higher scores. Closeness centrality indicates the shortest path, and we demonstrate in Fig. 2C that there were more nodes displaying the same numbers of connections, and that these nodes are relatively centralized in the network. Those nodes with more connections are relatively sparsely distributed. Path reflects the combination of all nodes in the network. Figure 2D shows the shortest pathway length distribution in the network, demonstrating that path length is mainly focused in the center. There were fewer extreme values in this analysis. The upper threshold was four and the lower threshold was one, suggesting that the majority of nodes in the network can be connected by a relatively short path.

Node degree distribution density is shown in Fig. 2E. We found that the number of nodes decreased as node degree increased, indicating that most nodes in the network are isolated. It is possible that during diseases occurrence, some key nodes change and interact with nearby nodes; therefore, co-expression occur and downstream biological processes are affected. In addition, we calculated the connection degree of each gene by topology to clarify its importance in the ceRNA network (Table S5). MEG3 (connection degree = 21), hsa-mir-195 (connection degree = 51) and ZEB1 (connection degree = 6) were considered the most important genes among the lncRNAs, miRNAs, and mRNAs, respectively. Because hsa-mir-195 was the gene with the highest connection degree in the ceRNA network, we conclude that it might exert a strong influence on EC tumorigenesis.

Positive correlations among the ceRNA expression levels

According to ceRNA theory, lncRNAs can positively regulate mRNAs by competitively competing with miRNAs (Fig. 2A). To verify this finding, regression analysis was performed between the DElncRNAs and DEmRNAs targeted by hsa-mir-195. Strong positive correlations were found between the DElncRNAs and DEmRNAs targeted by hsa-miR-195 (minimum requirement of an interaction score >0.3; Fig. 3; Table S6). Figures 3A–3F shows the top six interactions in which C2orf48 interacts with PSAT1, C2orf48 interacts with KIF23, C2orf48 interacts with CCNE1, C2orf48 interacts with CEP55, C2orf48 interacts with CBX2, and C2orf48 interacts with CDC25A under hsa-mir-195 regulation.

Biological processes and pathway enrichment analysis of the ceRNA network

Seventy-three mRNAs in the ceRNA network were subjected to Gene Ontology (GO) and pathway enrichment analysis using the DAVID 6.7 database. The first 20 biological functions that demonstrated statistical significance were analyzed by “GOplot”.R package. Pathways were visualized with Cytoscape software. GO analysis demonstrated that the biological processes significantly enriched by mRNAs in the ceRNA network include GO:0000187 ∼activation of MAPK activity, GO:0060412 ∼ventricular septum morphogenesis and GO:0045599 ∼negative regulation of fat cell differentiation (Fig. 4A; Table S7). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis demonstrated that these genes are involved in regulating microRNAs in cancer, p53 signaling pathways, the cell cycle, the Rap1 signaling pathway and other signaling pathways closely associated with tumorigenesis (Fig. 4B; Table S8).

Survival and gene expression analysis

To clarify how DElncRNAs, DEmiRNAs, and DEmRNAs in the ceRNA network affect the prognoses of patients suffering from EC, K-M survival analysis was performed for the 97 lncRNAs, 20 miRNAs and 73 mRNAs, separately. The results demonstrate that 92 of the 97 lncRNAs were significantly associated with prognosis based on their respective optimal cutoffs (Kassambara et al., 2014) and 59 DElncRNAs exhibited positive correlations with overall survival; in contrast, the remaining 33 DElncRNAs exhibited negative associations with overall survival (log-rank P < 0.05; Table S9). Sixteen of the 20 differentially expressed miRNAs correlated to prognosis. Ten DEmiRNAs positively correlated to overall survival, and the remaining 6 DEmiRNAs were negatively associated with overall survival (log-rank P < 0.05; Table S10). Sixty-one of the 73 DEmRNAs were highly relevant for overall survival based on their respective optimal cutoffs. Forty-four DEmRNAs positively correlated to overall survival, and the remaining 17 DEmRNAs were negatively correlated to overall survival (log-rank P < 0.05; Table S11). Figures 5A–5D shows the survival curves of the top four DElncRNAs (FAM41C, ADARB2.AS1, C8orf49 and MIR7.3HG), Figs. 5E–5H shows the survival curves of the top four DEmiRNAs (has-mir-195, has-mir-140, has-mir-205 and has-mir-200a) and Figs. 5I–5L shows the survival curves of the top four DEmRNAs (PSAT1, KIF23, MCM4 and CCNE1). Then, we verified the expression of four mRNAs and miRNAs through the TCGA endometrial cancer dataset, and found that hsa-mir-195 and hsa-mir-140 are poorly expressed in endometrial cancer (Figs. 5Q–5R). The survival curves showed that as expression decreases, patient prognoses worsen, suggesting that these two miRNAs promote endometrial cancer development by negatively regulating the expression of cancer-promoting genes. Hsa-mir-205 and hsa-mir-200a were highly expressed in endometrial carcinoma tissues (Figs. 5S–5T), and their survival curves showed that as gene expression increases, survival prognosis worsens. We conclude that both miRNAs negatively regulate the expression of tumor suppressor genes and thereby promote the occurrence and development of tumors. The four genes that are highly expressed in endometrial cancer tissues (PSAT1, KIF23, CCNE1 and MCM4) showed a pattern whereby increasing expression was associated with worsening five-year survival prognosis (Figs. 5M–5P).

Analysis of correlations between gene expression and clinicopathological parameters

The UALCA online database (Chandrashekar et al., 2017) was used to analyze correlations between mRNA expression of the above genes and patient clinicopathological parameters. We found that PSAT1, KIF23, CCNE1 and MCM4 expression was higher at different clinical stages than in the control group, and expression was higher in different pathological subtypes than in the normal group. At the same time, we found that the expression levels of all four genes was higher in the serous subtype than in endometroid type. All differences were statistically significant, suggesting that these four markers could also be used as biomarkers for molecular subtypes of EC (Figs. 6A–6H).

Validation of DEmRNAs and DEmiRNAs in the ceRNA network

The top four DEmRNAs and DEmiRNAs correlated with overall survival were selected for validation. Consistent with our earlier results, mean expression levels of all four mRNAs were significantly higher in EC tissues than in non-cancerous tissues in the GSE17025 dataset (Figs. 7A–7D). Mean expression levels of hsa-mir-195, hsa-mir-140 were significantly lower in EC tissues than in non-cancerous tissues in the GSE35794 dataset (Figs. 7E–7F), and hsa-mir-205 expression was significantly higher in EC tissues. However, due to the small number of patients, no significant differences were found in hsa-mir-200a expression (Figs. 7G–7H).

Flow chart of construction and analysis of ceRNA network

The flow chart of construction of lncRNAs-miRNAs-mRNA regulation network in endometrial carcinoma of is shown in Fig. 8.