Integrated analysis of a competing endogenous RNA network reveals an 11-lncRNA prognostic signature in ovarian cancer

Identification of significantly differentially expressed lncRNAs and mRNAs

According to the annotation information of Ensembl and miRbase databases, 14,148 lncRNAs, 19,645 mRNAs and 2,588 miRNAs were annotated. After removing the low-abundance genes, 8,860 lncRNAs, 18,438 mRNAs and 1,104 miRNAs were retained. Based on the given threshold, 3,279 differentially expressed lncRNAs (DElncRNAs), of which 2,210 DElncRNAs were upregulated and 1,069 DElncRNAs downregulated, and 4,789 differentially expressed mRNAs (DEmRNAs), of which 3,145 DEmRNAs were upregulated and 1,644 DEmRNAs downregulated, were identified between OVC tissues and normal tissues (Figure 1A and 2A).

Weighted correlation network analysis

The top 5000 most variant lncRNAs, measured by the median absolute deviation (MAD), were selected for the WGCNA. To obtain a higher average connectivity degree in the lncRNA group, the soft threshold power was set to 7, which was the lowest threshold to make the scale-free R² reach 0.85 (Figure 1B). Then, the cluster dendrogram was clustered based on the selected threshold, and 5 color modules were identified in the lncRNA group (Figure 1C). The lncRNAs that could not be included in any module were put into the gray module, which was removed in the subsequent analysis. The other four modules were represented in blue, yellow, turquoise and brown (Figure 1C). The relationship between modules and traits (OVC tissues and normal tissues) is shown in Figure 1D. Among these four modules, the turquoise module (R = 0.87, P < 0.01) and blue module (R = -0.84, P < 0.01) were significantly correlated with OVC (Figure 1D). Intra-modular analysis found that lncRNAs in the blue (R = 0.95, P < 0.01) and turquoise modules (R = 0.99, P < 0.01) were also highly correlated with OVC (Figure 1E, 1F).

Similarly, the top 5000 most variant mRNAs were selected for the WGCNA. As shown in Figure 2B, the soft threshold power was set to 18 to ensure a scale-free network, and we obtained eight modules for the next analysis (Figure 2C). Then, mRNAs in the eight color modules were continuously used to analyze the relationship between modules and traits (OVC tissues and normal tissues) (Figure 2D). The blue module (R = 0.89, P < 0.01) was positively associated with OVC, while the turquoise module (R = -0.84, P < 0.01) was negatively correlated with OVC (Figure 2D). Intra-modular analysis showed that mRNAs in both the turquoise module (R = 0.96, P < 0.01) and blue module (R = 0.99, P < 0.01) were highly related with OVC (Figure 2E, 2F).

Functional enrichment analysis for the overlapped mRNAs

The blue and turquoise modules, which were closely related to OVC, were selected and intersected with the above 4,789 DEmRNAs. In total, 1,827 mRNAs were obtained (Figure 3A, 3B). With a threshold of FDR < 0.01, GO enrichment analysis determined that 87 BP terms, 29 CC terms and 23 MF terms were significantly enriched in the overlapped mRNAs. KEGG pathway analysis revealed that 17 statistically significant signaling pathways were concentrated on the overlapped mRNAs. The top 10 terms of functional enrichment analysis are shown in Figure 3C.

Construction of the lncRNA-miRNA-mRNA regulatory network

Overall, 1,885 lncRNAs, 1,827 mRNAs and 1,104 miRNAs were used to build the ceRNA regulatory network. According to the hypothesis of ceRNA, the candidate ceRNA triplets should meet the following criteria: (1) the lncRNA-mRNA interaction shows significantly positive correlations; (2) the miRNA-mRNA interaction and lncRNA-miRNA interaction show significantly negative correlations; (3) the lncRNA-miRNA-mRNA triplets show significant results in hypergeometric testing. Finally, the lncRNA-miRNA-mRNA ceRNA regulatory network, including 201 lncRNAs, 85 miRNAs and 146 mRNAs, was established based on the interactions of 487 lncRNA-miRNA pairs and 248 miRNA-mRNA pairs (Figure 4A). The expression of these 201 lncRNAs and 146 mRNAs between 371 OVC tissues and 88 normal tissues is shown in Figure 4B.

Functional enrichment and PPI analyses of the mRNAs in the ceRNA regulatory network

Based on the functional enrichment analysis, the 146 mRNAs involved in the identified ceRNA regulatory network were remarkably enriched for 55 GO BP terms and 13 KEGG pathways (FDR < 0.05). The top 10 enriched GO BP terms, mainly including "regulation of phospholipase activity", "extracellular structure organization", "positive regulation of phospholipase activity" and several system development-related terms, are shown in Figure 4C. The top 10 enriched KEGG pathways, mainly including the "Hippo signaling pathway", "AGE-RAGE signaling pathway in diabetic complications", "PI3K-Akt signaling pathway" and several other cancer-related pathways, are shown in Figure 4D.

With the threshold of a minimum required interaction score > 0.4, the STRING database was used to construct a PPI network based on the 146 mRNAs involved in the identified ceRNA regulatory network. After wiping out the isolated nodes in the network, 96 nodes and 180 edges were mapped in the PPI network (Figure 4E). The top 3 mRNAs with the highest degrees were FGF2 (degree = 17), KIT (degree = 13) and NCAM1 (degree = 11).

Prognosis values of lncRNAs and mRNAs in the ceRNA network

To identify the mRNAs and lncRNAs with potential prognostic value in predicting OS, survival analyses of 146 mRNAs and 201 lncRNAs in the ceRNA regulatory network were conducted using the univariate Cox proportional hazards regression model. The results showed that 24 RNAs, including 12 mRNAs and 12 lncRNAs, were significantly associated with OS (Figure 5A, 5B). Among these significant RNAs, six lncRNAs (TYMSOS, LINC01619, LINC00488, CHRM3-AS2, AL391069.3 and AC104667.2) and two mRNAs (ABCG8 and MYCN) were negatively associated with OS (P < 0.05). In comparison, the remaining six lncRNAs (PCOLCE-AS1, MEF2C-AS1, LINC01558, HOXB-AS2, CACNA1G-AS1 and AC026904.1) and ten mRNAs (ZCCHC24, NBL1, SLC22A3, GFPT2, LRRC17, TCF15, PTGIS, FSTL3, PRDM6 and ARHGAP6) were positively associated with OS (P < 0.05). Based on their respective optimal cut-offs, Kaplan-Meier curves of these 12 mRNAs and 12 lncRNAs associated with OS were drawn (Figure 5C, 5D). The results of the Kaplan-Meier curve analysis and log-rank test of these 24 RNAs were proven to be consistent with the results of the univariate Cox regression analysis. Furthermore, we extracted the OS-specific ceRNA sub-network of these 24 RNAs. As shown in Figure 5E, The OS-specific ceRNA sub-network contained 4 RNA molecules, including 6 OS-related lncRNAs, 14 interacting miRNAs and 9 OS-related mRNAs.

Construction of the OVC-specific lncRNA-based prognostic signature

Through the univariate Cox regression analysis, we identified 30 potential OS-associated lncRNAs (Figure 6A). Then, we used the LASSO regression model to further identify an optimal subset of the lncRNA-based signature reliably associated with OS. As a result, 11 lncRNAs, including C4A-AS1, LINC02408, AC087521.1, LINC00488, AC010275.1, CHRM3-AS2, LINC00337, Z98257.1, AC104667.2, SOCS2-AS1 and CACNA1G-AS1, were identified for modeling (Figure 6B, 6C). Furthermore, a lncRNA-based risk score was developed by integrating the expression data for these 11 lncRNAs with corresponding coefficients weighted by the LASSO model as follows: Risk score = (0.147 × LINC02408 expression value) + (-0.066 × LINC00488 expression value) + (-0.161 × LINC00337 expression value) + (0.104 × Z98257.1 expression value) + (-0.129 × CHRM3-AS2 expression value) + (-0.14 × C4A-AS1 expression value) + (-0.156 × AC104667.2 expression value) + (0.088 × AC087521.1 expression value) + (-0.159 × SOCS2-AS1 expression value) + (0.096 × AC010275.1 expression value) + (0.111 × CACNA1G-AS1 expression value).

The risk score was calculated for all 371 OVC patients with the lncRNA-based prognostic signature. The accuracy of this prognostic signature in OVC long-term OS predictions was evaluated by ROC curves. The values of the area under the ROC (AUC) for the 1-year, 3-year and 5-year OS prediction models was 0.67, 0.688 and 0.702, respectively (Figure 6D). After estimating the maximally selected rank statistics, patients were divided into two groups based on the optimal risk score cut-off (0.612). Patients with risk scores higher than 0.612 were classified into the high-risk group (283 patients), whereas those with risk scores less than or equal to the cut-off value were allocated to the low-risk group (88 patients). Survival analysis suggested that high-risk patients were at a markedly increased risk of death in OVC and predicted a shorter OS (Figure 6F and 6H). The heatmap revealed that the expression of these 11 lncRNAs varied by risk score (Figure 6J). This prognostic signature also showed an excellent performance to predict disease-free survival (DFS) (Figure 6E, 6G, 6I, 6K). In the TCGA-DFS dataset (n = 177), the AUC for this prognostic signature for the 1-year, 3-year and 5-year DFS prediction models were 0.658, 0.676 and 0.728, respectively (Figure 6E). Based on the same cut-off (0.612), patients in the low-risk group had better survival time than those in the high-risk group and a decreased risk of death (Figure 6G, 6I).

Expression of the 11 lncRNAs was analyzed among normal tissues, tumor high-risk and tumor low-risk groups. As shown in Figure 7, all 11 lncRNAs among the three groups were significantly differentially expressed. Furthermore, the expression of Z98257.1 in the high-risk group was not significantly higher than that in the low-risk group; the other ten lncRNAs were significantly differentially expressed between every two-group comparison (Figure 7H). Coexpression analysis revealed that these 11 lncRNAs showed a weak or moderate-intensity of coexpression interactions as evaluated by the Pearson correlation coefficient (Figure 7L).

To investigate whether the prognostic value of the 11-lncRNA signature was independent of other clinical variables, we performed univariate and multivariate Cox regression analysis, which included risk score, age, race, TNM classification, venous invasion and lymphatic invasion as covariates. The results showed that either in the TCGA-OS dataset (HR = 1.36, P < 0.0001) or the TCGA-DFS dataset (HR = 1.32, P < 0.0001), the 11-lncRNA risk score was an independent prognostic factor (Table 1). In addition, age (HR = 1.02, P = 0.01) and race (HR = 0.55, P = 0.0052) were also independent prognostic factors of OS (Table 1). Therefore, the 11-lncRNA signature may offer an approach for risk assessment and predict the prognosis of OVC patients.

The analysis of the hallmark pathway gene sets highlighted that the 50 key signaling pathways in OVC samples were quite different between the high- and low-risk groups (Figure 8A, 8B). Of note, OVC samples in the high-risk group were most significantly enriched for TGF-, NOTCH-, WNT-, HEDGEHOG-, KRAS- and EMT (epithelial-mesenchymal transition) signaling pathways. In comparison, OVC samples in the low-risk group were most significantly enriched for MYC-, MTORC1- and some hormonal fluctuations and metabolism pathways.

Construction of a nomogram integrating the 11-lncRNA signature and clinical factors

Considering the clinical relevance and prognostic value of age and race, a nomogram based on the 11-lncRNA signature, age and race was established (Figure 8C and Table 1). Calibration curves for estimating 1-year, 3-year and 5-year survival showed that there were good correlations between the prediction made by the nomograms and the actual observation (Figure 8D). The results suggested that the nomogram could accurately predict 1-year, 3-year and 5-year survival of OVC patients.

Patients and samples

The RNA-Seq data and clinical phenotype information of OVC tissues obtained from the TCGA project and normal tissues obtained from the GTEx project were downloaded from the UCSC Xena database (http://www.genome.ucsc.edu) [46]. All the OVC patients were selected according to the following criteria: (1) patients had a histopathological diagnosis of OVC; (2) patients had no other tumors; (3) patients did not receive preoperative chemoradiation; (4) patients had complete RNA-Seq data and clinical prognosis data. Finally, 459 samples, including 371 OVC tissues and 88 normal tissues, were enrolled in our study. This study was performed following the TCGA and GTEx publication guidelines. As the data were all retrieved from TCGA and GTEx, approval from a local Ethics Committee was not necessary.

Data preprocessing

The mRNAs and lncRNAs were identified and annotated from the RNA-Seq data via the Ensembl database (http://www.ensembl.org), and miRNAs were annotated based on the miRbase database (http://www.mirbase.org) [47, 48]. The RNA-Seq data that could not be annotated in the database were discarded. The expression of lncRNAs, mRNAs and miRNAs with zero counts in more than 20% of samples was filtered out to remove low-abundance genes. [49, 50].

Differential expression analysis

Differentially expressed lncRNAs (DElncRNAs) and mRNAs (DEmRNAs) between OVC tissues and normal tissues were analyzed using the "edgeR" R package [51, 52]. P-values were corrected with a false discovery rate (FDR). Fold changes of expression levels (log2 absolute) ≥ 2 and FDR < 0.05 were set as the thresholds.

WGCNA

To identify the significant lncRNAs and mRNAs associated with the carcinogenesis of OVC, WGCNA analysis was performed according to the protocol and recommendations of the "WGCNA" R package [53] to construct independent signed networks and establish the weighted gene coexpression modules and module-trait relationship from the SOC samples and normal samples. The scale-free topology fitting index (R²) > 0.85 was set as the threshold to construct the weighted gene coexpression network. A minimum cluster size of 30 and a merge threshold function of 0.25 were chosen as the thresholds to identify coexpressed gene modules. A biweight mid-correlation coefficient (r) ≥ 0.8 and P-value < 0.05 were set as the thresholds to determine statistically significant modules with the carcinogenesis of OVC.

Functional enrichment analysis of the overlapped mRNAs

Enrichment analysis was performed to understand the potential mechanism and function of the overlapped mRNAs. Gene Ontology (GO) terms, including molecular function (MF), biological process (BP) and cellular component (CC) categories, and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the "ClusterProfiler" R package, with FDR < 0.01 as the threshold [54].

Construction of the lncRNA-miRNA-mRNA regulatory network

Based on the ceRNA hypothesis, the lncRNA-miRNA-mRNA ceRNA network was constructed for the overlapped mRNAs and lncRNAs through "GDCRNATools" following the following three steps: (1) miRcode (http://www.mircode.org/) database was used to predict the lncRNA-miRNA interactions, and miRTarBase database (http://mirtarbase.mbc.nctu.edu.tw) was used to predict the miRNA-mRNA interactions; (2) according to the common miRNA, the lncRNA-miRNA and miRNA-mRNA interactions were merged into the potential lncRNA-miRNA-mRNA ceRNA triples. Moreover, the hypergeometric test was used (P-value < 0.05 for the hypergeometric tests) to determine the significance of the shared miRNAs by lncRNA and mRNA pairs; (3) according to the ceRNA hypothesis, in the lncRNA-miRNA-mRNA ceRNA triples, the expression of lncRNA and mRNA pairs should be positively correlated, and the expression of lncRNA-miRNA pairs and miRNA-mRNA pairs should be negatively correlated. The Pearson correlation analysis was performed (P-value < 0.05 for the Pearson's correlation tests) to measure the expression correlation of lncRNA-miRNA, miRNA-mRNA and lncRNA-mRNA pairs [55–57]. Cytoscape (version 3.7.1) software was used to construct and visualize the lncRNA-miRNA-mRNA regulatory network [58].

Functional enrichment and protein-protein interaction (PPI) analyses of the mRNAs in the ceRNA network

To understand the underlying pathways and processes of the ceRNA regulatory network, the "ClusterProfiler" R package was used to analyze the functional profiles (GO BP terms and KEGG pathway) of the mRNAs in the ceRNA regulatory network [54]. With a threshold of FDR < 0.05, the results of the GO BP and KEGG enrichment analyses were displayed using the "GO plot" R package. To identify the protein-protein interactions between the mRNAs in the ceRNA regulatory network, STRING (Search Tool for the Retrieval of Interacting Genes) database (http://www.string-db.org/) was used to construct a PPI network (minimum required interaction score > 0.4), which was visualized by Cytoscape (version 3.7.1) software [58, 59]. By ranking the degree of connectivity between the mRNAs in the PPI networks, high-degree genes, the so-called "hub" genes, were identified using the Cytohubba plugin [60].

Survival analysis of lncRNAs and mRNAs in the ceRNA regulatory network

Univariate Cox proportional hazards regression analysis was performed to estimate the expression of lncRNAs and mRNAs in the ceRNA regulatory network with overall survival (OS) through the "survival" package in R software. Based on the gene expression of each patient with optimal cut-off values, which was determined by the "survminer" R package, the patients were divided into high-expression and low-expression groups. Then, using the "survival" R package, Kaplan-Meier curves and log-rank methods (Mantel-Haenszel test) were used to examine the statistically significant OS-related mRNAs and lncRNAs, which were determined by the univariate Cox regression analysis. P-value < 0.05 was considered statistically significant.

Construction of lncRNA-associated prognostic signatures

The potential OS-related lncRNAs chosen from the univariate Cox regression analysis (P-value < 0.1) were inserted into the least absolute shrinkage and selection operator (LASSO) analysis to calculate the coefficients. Ten-fold cross-validation was conducted to tune the optimal value of the penalty parameter λ, which gives the minimum partial likelihood deviance. According to the linear combination of the expression values weighted by the regression coefficient from the LASSO analysis, risk scores for each sample were calculated. All patients were classified into either the high-risk or low-risk group based on the optimal cut-off value of their risk score. The optimal cut-off value was identified using the "survminer" R package. The survival curve in the high-risk and low-risk groups was estimated and compared by the Kaplan-Meier method and the log-rank test. A P-value < 0.05 was considered statistically significant. The receiver operating characteristic (ROC) curves and areas under the ROC curves (AUC values) were performed to evaluate the sensitivity and specificity of the risk score model for survival prediction using the "survival ROC" R package.

Multivariate Cox regression analysis was performed to determine whether the risk score model was independent of other clinical features, including age, race, TNM classification, venous invasion and lymphatic invasion. The hazard ratio and 95% confidence intervals for each variable were calculated. A P-value < 0.05 was considered statistically significant.

Gene Set Variation Analysis (GSVA)

To further explore the biological function of the 11-lncRNA signature in the occurrence and development of OVC, we performed GSVA between high-risk and low-risk groups using the "GSVA" R package [61]. Each gene set associated with a pathway was trimmed to only contain unique genes to reduce pathway overlaps and pathway redundancies [62]. The hallmark gene sets were used as the reference gene set [63]. FDR < 0.05 was set as the cut-off criterion.

Construction of nomogram predictive models

To precisely predict the 1-year, 3-year and 5-year OS of OVC patients, a prognostic nomogram was established using the "rms" R package based on the 11-lncRNA signature, age and race. Calibration curves were used to compare the concordance between nomogram-predicted survival and observed survival.

Statistical method

For the correlation analysis, Pearson's correlation test was performed. For the survival analysis, Cox regression analysis and Kaplan-Meier log-rank test were performed. To analyze differences between groups, the Mann-Whitney test was used for continuous data, whereas the Chi-square test was used for categorical data.