A four-lncRNA signature for predicting prognosis of recurrence patients with gastric cancer

2.1 Data and preprocessing

We obtained RNA sequencing data (including mRNA and lncRNA) of 407 GC samples and miRNA sequencing data of 477 samples with the corresponding clinical information from the publicly accessible TCGA database (https://gdc-portal.nci.nih.gov/) based on Illumina HiSeq 2000 RNA Sequencing platform. A total of 287 samples with paired RNA and miRNA data as well as the corresponding overall survival (OS) information were selected as the training set. The detailed clinical information of samples in TCGA data sets is shown in Table S1.

Meanwhile, we launched a search in GEO database (https://www.ncbi.nlm.nih.gov/geo/) using “gastric cancer,” “stomach cancer,” and “homo sapiens” as the key words. Potential data sets were selected when the following criteria were met: the number of total samples ≥250 and the number of GC samples with corresponding clinical data ≥200. As a result, two data sets including GSE26253 and GSE62254 were chosen. GSE26253 [12] (GPL8432 Illumina HumanRef-8 WG-DASL v3.0 platform) contained gene expression data of 432 GC tissue samples with clinical information available (validation set 1), while GSE62254 [2] (Affymetrix Human Genome U133 Plus 2.0 Array platform) included gene profiling data of 300 GC tissue samples, among which 282 samples had the corresponding clinical data (validation set 2). The detailed clinical information of samples in GSE26253 and GSE62254 data sets is shown in Tables S2 and S3.

In order to supplement the information of targeted molecular therapy, we collected information about the targeted drugs from databases TCGA, xena.ucsc, and cbioporta. Unfortunately, we only got the information on whether the patients received “targeted treatment” or not; however, there were no specific information of targeted drug molecules. Meanwhile, the detailed clinical characteristics of the patients were analyzed and shown in Table S4. In order to exclude the effect of therapeutic schedule on the screening of the present lncRNAs, correlation analysis was used to eliminate the doubt in this aspect, and the results showed that four important lncRNAs were not related to radiotherapy, chemotherapy, and targeted therapy. Their detailed information is shown in Table S5.

Raw data were standardized as previously described by Chaudhary et al. [13]. HUGO Gene Nomenclature Committee (HGNC) [14] repository (http://www.genenames.org/) enrolled 4,313 lncRNAs, 19,197 protein-coding genes, and 1,914 miRNAs. Using the HGNC database, we annotated 13,105 mRNAs, 1,051 lncRNAs, and 413 miRNAs from the abovementioned data sets, according to RefSeq ID information.

2.2 Differential expression analysis between recurrence and nonrecurrence GC samples

We screened differentially expressed RNAs (DERs) including lncRNAs, mRNAs, and miRNAs between recurrence and nonrecurrence samples in the training set, using FDR <0.05 and |log₂ FC| > 0.263 as a selection threshold. The identified significant lncRNAs, miRNAs, and mRNAs consequently underwent two-way hierarchical clustering analysis based on centered Pearson correlation algorithm [15].

2.3 Prognostic model building and validation

In order to identify recurrence-related feature lncRNAs from the preselected differentially expressed lncRNAs, we employed and compared the support vector machine–recursive feature elimination (SVM–RFE) and random forest–out-of-bag (RF–OOB) algorithms for performance. With regard to SVM–RFE [16] algorithm (100-fold cross validation), the lncRNAs subset with the best accuracy was chosen to be the signature lncRNAs. Using RF–OOB [17] algorithm, the subset of lncRNAs with the minimal value of OOB error was selected to be the optimal feature lncRNAs. Principal component analysis (PCA) method was applied to compare the performances of the two algorithms. The feature lncRNAs identified by the superior approach was applied in further analysis.

Using training set, we did univariable Cox regression analysis on the identified recurrence-related signature lncRNAs. The significant lncRNAs were further included in multivariable Cox regression analysis. log-rank p value <0.05 defined the significance. The lncRNAs that were independent prognostic indicators were selected to construct a prognostic score formula as following:

where here β _lncRNA suggests multivariable Cox regression coefficient of lncRNAn; Exp_lncRNAn indicates the expression value of lncRNAn.

We calculated the risk score for each patient in the training set based on the formula. Patients in the training set were categorized by the median risk score into a high-risk group and a low-risk group. Similarly, the lncRNAs-based risk score formula was applied to distinguish patients in two validation sets. Recurrence-free survival (RFS) time of two risk groups was compared using Kaplan–Meier method and log-rank test. Sensitivity and specificity of the risk score model were assessed using the ROC curve.

Using the training set data, we conducted uni- and multivariable Cox regression analysis to evaluate the association of clinical features and risk score model status with RFS time of patients. Integrating independent prognostic clinical features with risk score status, nomogram was built; and the calibration curves were plotted to ascertain its predictive performance.

2.4 Statistical analysis

A wide range of packages in R software (version 3.4.1) were utilized for bioinformatics and statistical analyses of our study: limma package (version 3.34.7) for differential expression analysis; pheatmap package (version 1.0.8) for two-way hierarchical clustering analysis; e1071 [18] (version1.7-1, https://cran.r-project.org/web/packages/e1071); caret [19] (version 6.0-76, https://cran.r-project.org/web/packages/caret) packages for SVM–RFE method; bootstrap algorithm of randomForest [20] package (https://cran.r-project.org/web/packages/randomForest/) for RF–OOB method; psych [21] package (version 1.8.12) for PCA; survival package (http://bioconductor.org/packages/survivalr/) for uni- and multivariable Cox regression analysis; Kaplan–Meier curves; pROC package (https://cran.r-project.org/web/packages/pROC/index.html) for ROC curve analysis; and rms package (version 5.1-2) for nomogram building.

2.5 Construction of lncRNA-associated ceRNA network

For dissecting the underlying molecular mechanisms of the prognostic signature lncRNAs in GC biology, we designed a three-phase study. Initially, the relations between the signature lncRNAs with the differentially expressed miRNAs were predicted based on DIANA-LncBasev2 [22] database. The negatively correlated lncRNA–miRNA pairs were selected to build an lncRNA–miRNA network.

Second, starBase database (version 2.0, http://starbase.sysu.edu.cn/) was used for the prediction of potential target mRNAs of the miRNAs included in the lncRNA–miRNA network. starBase program compiles prediction results from TargetScan, PicTar, RNA22, PITA, and miRanda [23]. The differentially expressed mRNAs were mapped to the predicted target mRNAs that were enrolled in at least three of the five programs. The miRNA–mRNA pairs with negative correlation were selected to develop an miRNA–mRNA network.

Third, an lncRNA–miRNA–mRNA network was built with predicted lncRNA–miRNA and miRNA–mRNA pairs with negative correlation. Networks were visualized with Cytoscape [24] software (version 3.6.1). All mRNAs in the networks were subject to gene ontology (GO) function and kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analysis by Enrichr [25] tool (http://amp.pharm.mssm.edu/Enrichr/). P value <0.05 was considered significant.

3.1 Identification of differentially expressed mRNAs, lncRNAs, and miRNAs between recurrence and nonrecurrence GC patients

The training set (N = 287) consisted of 61 recurrence patients and 226 nonrecurrence patients. Totally 363 DERs between the recurrent and nonrecurrent patients were identified by performing differential gene expression analysis, which were consisted of 317 mRNAs (115 downregulated and 202 upregulated mRNAs; Table S6), 18 lncRNAs (8 downregulated and 10 upregulated lncRNAs; Table S7), and 28 miRNAs (9 downregulated and 19 upregulated miRNAs; Table S8; Figure 1a and b).

Figure 1

Differentially expressed lncRNAs, miRNAs, and mRNAs of recurrence and nonrecurrence patients. (a) Left image is a scatterplot for differentially expressed RNAs. Lateral axis represents nonrecurrence samples; vertical axis represents recurrence samples. Upregulated, downregulated, and nonregulated RNAs are shown in red, green, and blue spots, respectively. Right image shows ratios of up- and downregulated RNAs. (b) Heatmap displays differentially expressed RNAs in recurrence and nonrecurrence samples.

3.2 Construction of a four-lncRNA prognostic score risk model

We used SVM–RFE and RF–OOB algorithms to select feature lncRNAs informative of recurrence based on the abovementioned differentially expressed lncRNAs. Using SVF-RFM algorithm, we obtained a combination of eight-feature lncRNAs when achieving the highest accuracy = 0.895 (Figure 2a, Table S9). A set of seven lncRNAs was determined by RF–OOB algorithm (minimal OOB error = 0.173; Figure 2b; Table S10). The two groups of feature lncRNAs were compared by PCA. Regarding the eight-feature lncRNAs identified by SVM–RFE, the top four principal components in PCA approximately account for 80% of the total observed variances (Figure 3a). In contrast, the top three principal components for the seven-feature lncRNAs selected by RF–OOB algorithm could describe almost 80% of the total observed variances (Figure 3b). These results suggested that the seven-feature lncRNAs were more indicative of recurrence than the eight-feature lncRNAs and were thus used in further analysis.

Figure 2

Feature lncRNAs selected by SVM–RFE (a) and RF–OOB (b). The red dashed line shows the number of lncRNAs achieving the highest accurracy (a) or the smallest OOB error (b).

Figure 3

Comparative analysis of SVM–RFE (a) and RF–OOB (b) by PCA approach. Upper images illustrate 3D PCA plot of PC1, PC2, and PC3. Recurrence and nonrecurrence samples are labeled in red triangles and blue balls, respectively. Below images exhibit cumulative contribution of different PCs. The red dashed line indicates cumulative proportion of 80%.

By subjecting the seven-feature lncRNAs to a univariable Cox regression analysis, we found six prognosis-related lncRNAs (p value <0.05), which further underwent multivariable Cox regression analysis. Four lncRNAs (LINC00843, SNHG3, C21orf62-AS1, and MIR99AHG), which were independent predictors of prognosis, were depicted in Table 1. A risk score formula was created based on the expression of the four signature lncRNAs for prognosis prediction as follows:

We calculated the risk score for each patient in the training set and ranked all the patients based on four-lncRNA signature. With the median risk score as the cutoff, the training set was categorized into a high-risk group (N = 144) and a low-risk group (N = 143). Significantly longer RFS time was observed in the low-risk patients relative to the high-risk patients (p value = 1.448 × 10⁻³, HR = 2.329[1.363–3.979]; Figure 4a). We validated the four-lncRNA signature risk score in two validation sets. As shown in Figure 4b and c, either validation set was divided by the risk score into two risk groups with statistical significance in RFS time (validation set 1: p value = 3.733 × 10⁻², HR = 1.368[1.017–1.839]; validation set 2: p value = 1.444 × 10⁻², HR = 1.555[1.089–2.221]). The risk score of each sample in the data set TCGA, GSE26253, and GSE62254 is shown in Table S11(1–3), respectively. The AUC values of ROC curves were 0.936, 0.827, and 0.822 for training set, validation set 1, and validation set 2, separately (Figure 4d). These observations demonstrated robust predictive performance of the four-lncRNA signature risk score in GC.

Figure 4

Kaplan–Meier estimates of recurrence-free survival and ROC curves for patients in training set and two validation sets by the four-lncRNA signature. (a–c) Kaplan–Meier curves for training set, validation set 1, and validation set 2. (d) ROC curves and AUC values of three data sets.

3.3 Building nomogram based on gender, histologic grade, and risk score model status

By performing univariable Cox regression analyses to assess the relationship of clinical variables and risk score model status with RFS time of patients in the training set, gender (HR = 1.761, 95% CI = 0.967–1.014, p value = 4.92 × 10⁻²), neoplasm histologic grade (HR = 2.181, 95% CI = 1.232–3.856, p value = 6.19 × 10⁻³), and the four-lncRNA risk score model status (HR = 2.329, 95% CI = 1.363–3.979, p value = 1.45 × 10⁻³) were statistically significant (Table 2). Furthermore, multivariable Cox regression analysis incorporating risk score model status with gender and neoplasm histologic grade was performed. Gender (HR = 1.735, 95% CI = 1.179–3.078, p value = 4.93 × 10⁻²), neoplasm histologic grade (HR = 2.063, 95% CI = 1.156–3.682, p value = 1.43 × 10⁻²), and risk score model status (HR = 2.059, 95% CI = 1.198–3.537, p value = 8.91 × 10⁻³) were found to be independent prognostic predictors (Table 2), indicating that prognostic value of the four-lncRNA risk score is independent of other clinical features.

A composite nomogram on the basis of gender, neoplasm histologic grade, and risk score status was established to predict the probability of 3-year and 5-year RFS (Figure 5a). Calibration plot for goodness of fit of the model exhibited good consistence between the predicted and actual RFS time (Figure 5b).

Figure 5

Nomogram to predict 3-year and 5-year probability of recurrence-free survival using TCGA set. (a) nomogram based on gender, neoplasm grade, and risk score model status. RS model status = four-lncRNA risk score model status. (b) calibration curves of actual and predicted probabilities.

3.4 Characterization of ceRNA network of four signature lncRNAs and functional annotation

Using DIANA-LncBasev2, we analyzed the correlations between the four signature lncRNAs and the differentially expressed miRNAs. As a result, we procured 10 negatively correlated lncRNA–miRNA pairs and constructed an lncRNA–miRNA network. There were 3 lncRNAs (upregulated MIR99AHG, downregulated LINC00843, and SNHG3) and 10 miRNAs (five upregulated and five downregulated miRNAs) in the network (Figure 6). As to the 10 miRNAs (miR-7, miR-552, miR-4676, miR-1304, miR-2110, miR-216a, miR-205, miR-487a, miR-551b, and miR-34b), we predicted potential target mRNAs by using starBase, which were further mapped by the 321 differentially expressed mRNAs. The resulting 351 miRNA–mRNA pairs with negative correlation were used to develop a miRNA–mRNA network. Figure 7 exhibited this network of 10 miRNAs and 178 mRNAs (54 downregulated and 124 upregulated mRNAs). Finally, the preselected negatively correlated lncRNA–miRNA pairs and miRNA–mRNA pairs were integrated into a ceRNA regulatory network. Three signature lncRNAs, 10 miRNAs, and 178 mRNAs were observed in the network (Figure 8). Their detailed information are shown in Tables S12 and S13.

Figure 6

Signature lncRNA–miRNA network. Diamonds and triangles stand for lncRNAs and miRNAs, respectively.

Figure 7

The miRNA–mRNA network. Balls and triangles denote mRNAs and miRNAs, respectively.

Figure 8

A ceRNA network of signature lncRNAs, miRNAs, and mRNAs. Diamonds, triangles, and balls represent signature lncRNAs, miRNAs, and mRNAs, separately. Red and black links suggest lncRNA–miRNA and miRNA–mRNA interactions, respectively.

The GO function and KEGG pathway-enrichment analysis for the mRNAs in the ceRNA network was conducted. These mRNAs were significantly enriched in 19 GO terms predominately related to Wnt-signaling pathway, cell proliferation, migration, and other biological processes (Table 3). Five significant KEGG-signaling pathways were noted, including Wnt-signaling pathway, cell adhesion molecules, GC, regulation of lipolysis in adipocytes, and steroid hormone biosynthesis pathways (Table 3). These results revealed that four signature lncRNAs-related ceRNA regulation participated in a variety of biological functions and signaling pathways in GC.