Competing endogenous RNA network identifies mRNA biomarkers for overall survival of lung adenocarcinoma: two novel on-line precision medicine predictive tools

Protocol approval

Data were collected as described in our previous studies (Zhang et al., 2019a; Zhang et al., 2019b). The study data in the current study were downloaded from The Cancer Genome Atlas (TCGA) database (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga), UCSC Xena database (https://xena.ucsc.edu/), cBioPortal database (http://www.cbioportal.org/), and Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The current study did not require ethics approval, as all study datasets were downloaded and analysed in accordance with the corresponding data policies of previous databases.

The gene expression information of model cohort

Data were collected as described in our previous studies (Zhang et al., 2020; Zhang et al., 2019b). The current study downloaded the original RNA expression information from TCGA database (https://tcga-data.nci.nih.gov/docs/publications/tcga/) on the Illumina HiSeq 2000 RNA Sequencing platform. The RNA symbols were determined according to GENCODE Version 29 (https://www.gencodegenes.org/human/). The RNA expression dataset contained 8,848 lncRNAs and 21,138 mRNAs from 535 lung adenocarcinoma samples and 59 normal samples. The miRNA expression dataset was obtained from the UCSC Xena database (https://xena.ucsc.edu/). The miRNA expression dataset involved 1,881 miRNAs from 518 lung adenocarcinoma and 46 normal samples. According to the median of the original gene expression value, the original expression values were converted into “0” for lower expression and “1” for higher expression.

Differentially expressed analyses

Differentially expressed analyses between tumour samples and normal samples were performed by the edgeR package with R software (version 3.5.2) as described in our previous studies (Zhang et al., 2019b; Zhang et al., 2018). For differentially expressed analyses, a false discovery rate (FDR) less than 0.05 and a ratio value of 1.5 between tumour tissues and normal tissues were defined as the thresholds.

The clinical information of model cohort

Data were collected as described in our previous studies (Zhang et al., 2019b; Zhang et al., 2018). The clinical information was obtained from the cBioPortal database. Patients with LUAD with inadequate survival information or OS times less than one month (for living patients) were removed (n = 19). Therefore, 503 LUAD patients were identified with adequate OS information. There were 494 patients with adequate gene detection information and survival information when performing intersection between clinical datasets on the gene detection dataset.

The corresponding information of validation cohort

Data were collected as described in our previous studies (Zhang et al., 2019b; Zhang et al., 2018). The current study downloaded GSE37745 and GSE50081 as external validation datasets from the GEO database. Gene expression information was detected on the GPL570 platform (Affymetrix Genome U133 Plus 2.0 Array). The verification group contained 233 LUAD patients and 22,850 RNA expression count values.

Statistical analysis

The statistical analyses were carried out as described in our previous studies (Zhang et al., 2020; Zhang et al., 2019b): continuous variables are presented as the mean ± standard deviation (SD) or median (first percentile, third percentile) as appropriate. Continuous variables were compared by t-test or Mann–Whitney U test as appropriate. Categorical variables were compared through the chi-squared test or Fisher’s exact test as appropriate. Time-dependent receiver operating characteristic (ROC) curves were used to assess the predictive performance of prognostic models. Decision curve analysis (DCA) is an assessment method to compare the predictive performance of prognostic models (Localio & Goodman, 2012; Vickers et al., 2008; Vickers & Elkin, 2006). The statistical analyses were carried out using SPSS Statistics 19.0 (SPSS Inc., USA) and R software (version 3.5.2) with the following packages: “GOplot”, “timeROC”, “rms”, “pROC”, “survival”, “clusterProfler” and “glmnet”, as described in our previous studies (Zhang et al., 2020; Zhang et al., 2019b). A P value < 0.05 was the threshold value for statistical significance.

Study cohorts

In the model group (Doc S1), 182 patients (36.8%) out of 494 patients died during the follow-up period, while 127 patients (54.5%) out of 233 patients died in the verification group (Doc S2). The basic features of the two groups are compared in Table 1.

Differentially expressed analyses

Based on the given threshold, 3,310 upregulated lncRNAs, 675 downregulated lncRNAs, 95 upregulated miRNAs, 30 downregulated miRNAs, 4,913 upregulated mRNAs, and 1,921 downregulated mRNAs were identified as differentially expressed RNAs between LUAD tissues and normal tissues. Based on previous differentially expressed mRNAs, 2,982 mRNAs were identified as prognostic mRNA indicators by univariate Cox analyses.

Construction of a competing endogenous RNA network

To obtain the lncRNA-miRNA pairs, the miRcode database (http://www.mircode.org/) was searched according to the differentially expressed lncRNAs. To obtain the miRNA-mRNA pairs, TargetScan (http://www.targetscan.org/), miRDB (http://mirdb.org/), and miRTarBase (https://bio.tools/mirtarbase) were searched according to the previous lncRNA-targeted miRNAs. The miRNA-predicted mRNAs that could be searched in three databases were defined as miRNA-targeted mRNAs. Then, we determined the interaction between miRNA-targeted mRNAs and prognostic mRNAs to identify prognosis-associated miRNA-targeted mRNAs. Finally, twenty-two lncRNAs, twenty-nine miRNAs, and seventy-three mRNAs were used to build the ceRNA network for the overall survival of Patients with LUAD. The ceRNA regulatory network was depicted by Cytoscape software (Fig. 1).

Functional enrichment analyses

To explore the biological functions of prognostic mRNAs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed on prognostic mRNAs in the ceRNA network. The bar plot and dot plot of prognostic mRNAs are depicted in Figs. 2A and 2B. GO terms and KEGG pathways of prognostic mRNAs are depicted in Figs. 3A and 3B.

Development of a prognostic nomogram

The previous prognostic mRNAs were used to construct a prognostic model for overall survival. The information on these mRNAs is summarized in Table 2. The risk score of the prognostic model was calculated by using the following formula: risk score = (−0.7305*DNAJC27) + (0.4776*NPAS2) + (0.3941*PHKA1) + (−0.5537*CDADC1) + (0.4792*PTGFRN) + (0.3980*DDIT4) + (−0.4133*SCAMP5) + (−0.4367*SPRY2) + (0.3850*CSE1L) + (−0.3940*ELAVL4) + (0.3293*TRIM29) + (0.3685*LPGAT1) + (−0.3279*TRPC3) + (0.3257*DCBLD2). A prognostic nomogram chart is presented in Fig. 4. Kaplan–Meier survival curves (Fig. S1 ) demonstrated that these 14 prognostic mRNAs were significantly correlated with OS (P < 0.05). The results of differential expression analysis of 14 enrolled mRNAs are listed in Table 3.

Predictive performance of the prognostic model

Based on the median value, LUAD patients in the model group could be classified into a high-risk subgroup and a low-risk subgroup. The overall survival rate (Fig. 5A) in the high-risk group was significantly poorer than that in the low-risk group (P < 0.001). Harrell’s concordance indexes (C-indexes) of the prognostic signature for overall survival in the model group were 0.786 (95% CI [0.744–0.828]), 0.736 (95% CI [0.694–0.778]), and 0.766 (95% CI [0.724–0.808]) for 1-year, 2-year, and 3-year OS, respectively (Fig. 5B). The scatter plot and the interaction distribution scatter plot are presented in Fig. 5C and Fig. 5D. Calibration curves are depicted in Fig. S2A for 1 year, Fig. S2B for 2 years, and Fig. S2C for 3-year overall survival.

External validation of the prognostic model

Kaplan–Meier analysis (Fig. 6A) indicated that there was a significant difference in OS between the low-risk group and the high-risk group in the validation dataset (P < 0.001). The C-indexes of the prognostic model for overall survival in the validation dataset were 0.948 (95% CI [0.912–0.984]), 0.936 (95% CI [0.900–0.972]), and 0.890 (95% CI [0.854–0.926]) for 1-year, 2-year, and 3-year overall survival, respectively (Fig. 6B). Calibration curves are depicted in Fig. S3A for 1 year, Fig. S3B for 2 years, and Fig. S3C for 3-year overall survival.

Independence assessment of the prognostic model

Table 4 shows that the prognostic model was an independent influencing factor for OS after adjustment for confounding effects in the model dataset. In the validation dataset, multivariate Cox regression analyses demonstrated that the prognostic model was an independent influencing factor for OS. DCA is shown in Fig. S4A for 1-year OS, Fig. S4B for 2-year OS, and Fig. S4C for 3-year OS. The clinical impact curve is presented in Fig. S4D.

Smart Cancer Survival Predictive System

A precision medicine predictive tool named the Smart Cancer Survival Predictive System was developed for predicting the prognosis of LUAD patients. The Smart Cancer Survival Predictive System (Fig. 7) is available at https://zhangzhiqiao2.shinyapps.io/Smart_cancer_predictive_system_9_LUAD_E1002/.

The Smart Cancer Survival Predictive System could predict full-time mortality risk prediction curves for one specific patient (Fig. 7A). Figure 7B, demonstrates the mortality rate predicted percentage and 95% confidence interval at different user-defined time points.

Gene Survival Analysis Screen System

A second precision medicine predictive tool named the Gene Survival Analysis Screen System (Fig. 8) was developed to explore the survival features of prognostic mRNAs. The Gene Survival Analysis Screen System is available at https://zhangzhiqiao7.shinyapps.io/Gene_Survival_Analysis_Screening_System_9_LUAD_E1002/. Figure 8A depicts and compares the survival curves between two defined subgroups. Figure 8B displays the results of univariate survival analysis for selected variables.