A novel five-lncRNA signature panel improves high-risk survival prediction in patients with cholangiocarcinoma

Establishment of a five-lncRNA signature predictive model from the TCGA cohort

Based on the screening criteria, we obtained1192 differentially expressed lncRNA, including 744 up-regulated and 448 down-regulated. Among them, 33 lncRNA showed >4-fold decreased expression including HULC, and 51 lncRNA exhibited >4 fold increased expression (Figure 1A). Unsupervised hierarchical cluster analysis showed that the expression of differentially expressed lncRNA distinguished CCA samples from the normal samples (Figure 1B).

Univariate Cox regression analysis assessing the association between 1192 differentially expressed lncRNA and CCA patients’ overall survival in TCGA cohort, 10 lncRNA (LINC01336; AP000943.4; AC006504.8; AC090114.2; AC004921.1; AC134682.1; AL449106.1; AL359715.5; AC016876.1 and HULC) were selected for further multivariate Cox regression analysis. In this multivariate Cox regression analysis, five lncRNA (AC006504.8, AC090114.2, AP000943.4, AL359715.5, and HULC) were identified as independent predictors for CCA survival outcomes (Table 1). Among them, the expression of AC006504.8, AC090114.2, and AP000943.4 was up-regulated (Figure 1C–1E); expression of AL359715.5 and HULC was down-regulated (Figure 1F, 1G). Similarly, these 5 lncRNAs also showed similar behaviors in paired carcinoma and paracancerous tissues (Supplementary Figure 1). Based on these five lncRNA expression levels and their corresponding coefficient in multivariate Cox regression analysis, a five-lncRNA signature predictive model for CCA patients prognosis was established as follows: Risk Score= (0.6542 × expression level of HULC) + (-1.2388 × expression level of AL359715.5) + (1.3769 × expression level of AC006504.8) + (-3.6697 × expression level of AC090114.2) + (-1.5165 × expression level of AP000943.4).

Performance evaluation of the five-lncRNA signature model for CCA prognosis in the discovery cohort

We calculated the risk score for each CCA patient in the TCGA cohort based on the five-lncRNA signature model. According to the median of log(Risk Score), -0.0995 was set as the cutoff level to stratify CCA patients into a high-risk or low-risk group(sensitivity is 82.01%, specificity is 86.02%). Finally, 18 patients whose log (Risk Score) >-0.0995 were assigned to the high-risk group, and 18 patients whose log(Risk Score) <-0.0995 were assigned to the low-risk group (Figure 2A–2C). The KM curves show that patients in the high-risk group had a worse prognosis than those in the low-risk group (median OS is 18.5 months vs. 60 months; log-rank p=0.002) (Figure 2D). The univariate Cox regression analysis between the lncRNA-based risk score and CCA patient survival score showed that high-risk patients have a worse prognosis (HR=6.760; 95%CI, 1.572-29.068; p=0.008). The multivariate Cox regression analysis revealed that high-risk score is an independent predictor for CCA patient survival after adjusting the clinical covariate including recurrence status and residual tumor status (HR=4.879; 95%CI, 1.587-14.996; p=0.006) (Table 2).

The AUC value and F1 score were calculated to assess the performance of the five-lncRNA signature model. Among the down-regulated lncRNA, the AUC for AL359715.5 and HULC were 0.777 and 0.821, respectively. The AUC for up-regulated lncRNA AC090114.2, AC006504.8, and AP000943.4 were 0.727, 0.728 and 0.77, respectively. For the whole model, the AUC was 0.938 (Figure 2E), and the F1 score was 0.7222. These results indicate that the five-lncRNA signature model has a high predictive value for CCA prognosis.

Performance validation of the five-lncRNA signature model in the validation cohort

We validated the prediction ability of the five-lncRNA signature model in the WMU cohort to assess the robustness of the model for survival prediction in patients with CCA. Based on the five-lncRNA model and cutoff point derived from the TCGA cohort, patients in the WMU cohort were divided into high-risk group (n=54) and low-risk group(n=36). We compared the KM curve between these two groups and found that the OS of patients in the high-risk group was worse than those in the low-risk group (p<0.001, Figure 3B). The AUC value was 0.816 in the WMU cohort (Figure 3C).

Stratified survival analysis of five-lncRNA model in the TCGA cohort

Further stratified analysis showed that the five-lncRNA signature model could accurately distinguish the prognosis between the high-risk group and the low-risk group in young patients <60 years (n=14, p=0.032, Figure 4E) and older patients ≥60 years (n=22, p=0.037, Figure 4A). Similarly, stratifying the patients based on the stage of disease, revealed that the five-lncRNA signature model has good discriminatory ability for earlier-stage patients (n=28, p=0.007, Figure 4F) and advanced-stage patients (n=8, p=0.028, Figure 4B). For patients with or without recurrence, the five-lncRNA model could divide patients into the high-risk or low-risk groups in those with recurrence (n=19, p=0.005, Figure 4C) and without recurrence (n=17, p=0.02, Figure 4G). Moreover, the five-lncRNA signature model could separate the high-risk group and low-risk group for patients with tumors (n=19, p=0.008, Figure 4D) and those who were tumor-free (n=15, p=0.18, Figure 4H). Multivariate Cox regression analysis combined with stratification analysis showed that there was no significant difference in OS between the high-risk and low-risk groups with five-LncRNA markers in tumor-free patients, and it this suggests that patients in the early stages of tumor development may benefit significantly from these prognostic biomarkers.

Identifying the functions of the five-lncRNA signature model

Co-expression analysis showed significant co-expression of 1429 DPCGs, 1440 DPCGs, 300 DPCGs, 495 DPCGs, and 552 PCGs with HULC, AL359715.5, AP000943.4, AC006504.8, AC090114.2, respectively. Functional enrichment analysis indicated that 72 GO biological processes (BP) terms, 21 GO cellular components (CC) terms, and 35 GO molecular functions (MF) terms were enriched for HULC-related DPCGs. Biological processes were mainly involved in the oxidation-reduction process, xenobiotic metabolic process, metabolic process; cellular components were mainly involved in extracellular exosome, mitochondrial matrix, blood microparticle; molecular functions were mainly involved in oxidoreductase activity, electron carrier activity, monooxygenase activity (Supplementary Figure 2A, 2C). There was significant enrichment of 60 KEGG pathways in HULC-associated DPCGs, including leucine, isoleucine and valine degradation, complement and coagulation cascades, fatty acid degradation, carbon metabolism and chemical carcinogenesis (Supplementary Figure 2B, 2C).

47 GO BP terms, 11 GO CC terms, and 33 GO MF terms were enriched for AL359715.5-related DPCGs, whose biological processes were mainly associated with drug metabolic process, lipid metabolic process, lipoprotein metabolic process; cellular components were mainly associated with organelle membrane, mitochondrion, peroxisome; molecular functions were mainly associated with iron ion binding, heme binding, cholesterol transporter activity (Supplementary Figure 3A, 3C).

56 KEGG pathways were enriched for AL359715.5-related DPCGs, which were primarily linked to Drug metabolism - cytochrome P450, Metabolism of xenobiotics by cytochrome P450, Retinol metabolism, Peroxisome, and Cholesterol metabolism (Supplementary Figure 3B, 3C).

16 GO BP terms, 18 GO CC terms and 4 GO MF terms were enriched for AC006504.8-related DPCGs, whose biological processes were mainly involved in cell division, sister chromatid cohesion, mitotic nuclear division; cellular components were mainly involved in condensed chromosome kinetochore, midbody, nucleoplasm; molecular functions were mainly involved in protein binding, ATP binding, and cadherin binding involved in cell-cell adhesion (Supplementary Figure 4A, 4C). It was a significant enrichment of 7 KEGG pathways in AC006504.8-related DPCGs, including DNA replication, cell cycle, the p53 signaling pathway, Fanconi anemia pathway and progesterone-mediated oocyte maturation (Supplementary Figure 4B, 4C).

7 GO BP terms, 2 GO CC terms and 1 GO MF terms were raised in AC090114.2-associated DPCGs, whose biological processes were mainly related to sister chromatid cohesion, DNA replication initiation, G1/S transition of mitotic cell cycle; cellular components were mainly associated with cytosol and MCM complex; molecular functions were mainly associated with protein binding (Supplementary Figure 5A, 5C). 4 KEGG pathways were raised in AC090114.2-related DPCGs, which were mainly related to DNA replication, cell cycle, cellular senescence and oocyte meiosis (Supplementary Figure 5B, 5C).

1 GO CC terms and 1 GO MF terms could be found for DPCGs related to AP000943.4, whose cellular components were associated with cytoskeleton; molecular functions were associated with extracellular matrix organization and involved a KEGG pathway for Human papillomavirus infection (Data not shown).

Functional evaluation of common DPCGs for five-lncRNA signature model

The intersection of the DPCGs corresponding to the five-LncRNA signature model showed that 171 DPCGs were shared by this five-LncRNA signature (Figure 5A). The common KEGG pathway is cell cycle, DNA replication, oocyte meiosis, Fanconi anemia pathway, and progesterone-mediated oocyte maturation (Figure 5B and Supplementary Table 1).

GSEA between the high-risk group and low-risk group

Through GSEA analysis, we clarified the significant difference in survival between the high-risk and low-risk groups. The results showed significant enrichment of markers including the "complement pathway" in the high-risk group. Pathways including IL-2 Receptor Beta Chain in T cell Activation, Keratinocyte Differentiation, T cell receptor pathway, and Neurotrophin signaling pathway were enriched in the low-risk group [18], many of these being closely related to the occurrence and development of cancer [19] (Figure 6 and Supplementary Table 2).

Screening differentially expressed lncRNA in CCA patients in the discovery cohort

Pre-processed level 3 RNA sequencing count data and relative clinical information for CCA patients were obtained from the Genomic Data Commons Data Portal database (TCGA, https://portal.gdc.cancer.gov/projects/TCGA-CHOL) [43]. Ultimately, 60,483 RNA-ENSG_ID expression profiles in 36 CCA patients were included for further analysis as the discovery cohort (TCGA cohort).

The Gencode.v27.long_noncoding_RNAs.gtf compressed file was downloaded from the GENCODE database (release 27) (http://www.gencodegenes.org/) and the transformed data (antisense, lincRNA, and sense_intronic) was determined as lncRNA [44]. It was then filtered by removing the exon-expressing lncRNA from any known coding gene by GENCODE-based gene annotation. A total of 59,264 lncRNA-ENSG_IDs, and RNA-ENSG_ID and lncRNA-ENSG_ID were intersected, 13,126 lncRNA-Gene_names were obtained (Supplementary Figure 6). Samples in which the lncRNA with RPKM expression value is 0 more than 20% were removed. Finally, 3651 lncRNA were keptfor further analysis. Furthermore, we downloaded Homo_sapiens.GRCh38.91.chr.gtf zip file from Ensembl genome browser 91 (ftp://ftp.ensembl.org/pub/release-91/gtf/homo_sapiens) to obtain 57,000 RNA-Gene_symbol expression profiles, then used "mygene" package in R to find the RNA with human entrezgene symbol (GeneID) in the 57,000 RNAs.

Differentially expressed genes (DEGs) between 36 CCA tissues (30 iCCA tissues, 6 pCCA/dCCA tissues) and 9 normal tissues with the |log₂ fold change (log₂ FC)| ≥1 and false discovery rate(FDR) <0.01 were considered as selection criteria for subsequent analysis. At the same time, differentially expressed analysis was performed on the RNAs identified with GeneID.

Establishment and performance evaluation of a lncRNA-based prediction model in the discovery cohort

Firstly, we assessed the association between differentially expressed lncRNA and the overall survival (OS) of CCA patients in the TCGA cohort by univariate Cox proportional hazards regression analysis. LncRNA that were statistically significant were included in further multivariate Cox regression analysis. Then, using the multivariate Cox regression analysis, we identified independent lncRNA predictors. Furthermore, the corresponding coefficients of each lncRNA in the Cox regression model were obtained to calculate the risk score [45]. Finally, a lncRNA-based prediction model to evaluate the CCA patient survival outcomes was established as below:

$\text{LncRNA-based Risk Score}={\sum }_{i=1}^N(Coei\times EVi)$

In this formula, N represents the number of prognostic lncRNA, Coe_i is the coefficient of the lncRNA_i in the multivariate Cox regression analysis, EV_i represents the expression level of the lncRNA_i.

We used the lncRNA-based model to calculate the risk score for each CCA patients in the TCGA cohort. Setting the median value of log (Risk Score) as a threshold, CCA patients were divided into high-risk and low-risk groups and the difference of OS between these two groups was compared. Then univariate and multivariate Cox regression analysis was performed to determine whether the lncRNA signature was an independent predictor variable of other clinicopathologic features for survival outcomes. Further stratification analysis was conducted on clinicopathologic characteristics, which were statistically significant in a multivariate Cox regression model to determine the lncRNA signature model’s predictive capacity within the same clinical features. We calculated the area under the time-dependent receiver operating characteristic (ROC) curve (AUC) within a 3-year survival period to evaluate the sensitivity and specificity of the lncRNA model to predict survival outcomes.

Verification of lncRNA signature model for survival prediction in the validation cohort

We performed validation of discovered lncRNA in fresh frozen tissues from 90 CCA patients who underwent surgery between November 2012 and December 2015 at the First Affiliated Hospital of Wenzhou Medical University (WMU cohort). The patient inclusion criteria were as follows: (1) pathological diagnosed with primary CCA; (2) had completed clinicopathological and follow-up monitoring; (3) had no anti-tumor treatment before this surgical resection. Exclusion criteria were as follows: (1) previous radiofrequency or other anti-tumor treatment before surgery; (2) patients who were lost to follow-up after surgery. The study was approved by the institutional review boards of First Affiliated Hospital of Wenzhou Medical University and written informed consent was obtained from all patients. Patients demographics and clinicopathological characteristics are shown in Supplementary Table 3.

The lncRNA expression of primary CCA tumor fresh frozen samples was assessed by real-time quantitative PCR. RNAeasy mini kit (Qiagen, CA, USA) was used to extract the total RNA. High-Capacity cDNA Reverse Transcription kit from Applied Biosystems (Grand Island, NY, USA) was used to synthesize cDNA from 2 μg of total RNA. Semi-quantitative detection of mRNA was performed using an ABI 7300 RT-PCR system. Relative quantification of mRNA levels was performed with 18S ribosomal RNA as an internal reference gene and data from the ΔΔCt method. The statistical mean and standard error were determined by the ΔCt value. All data were independently inputed three times. The primer sequences used in the present study are shown in Figure 4A.

Through the same lncRNA signature model and cutoff level derived from the discovery cohort, patients in the WMU cohort were divided into high-risk and low-risk groups. Then, we investigated the performance of the lncRNA signature model in the WMU cohort.

Co-expression and functional enrichment analysis

The Spearman correlation coefficient was calculated between lncRNA expression level and the differentially expressed protein-coding genes (DPCGs). The DPCGs with Spearman correlation coefficient greater than 0.50 were considered to be lncRNA-related DPCGs. Gene Ontology (GO), Cell Component (CC), Biological Process (BP), Molecular Function (MF), Cell Component (CC) and Kyoto Encyclopedia Gene and Genome (KEGG) pathway enrichment assays were performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID, version 6.8, https://david.ncifcrf.gov/) [46]. Gene Set Enrichment Analysis (GSEA 3.0) analyses were carried out to elucidate the survival difference between high-risk and low-risk groups. GSEA analyses were implemented with java software GSEA (http://software.broadinstitute.org/gsea/index.jsp) [18]. GO (BP, MF, CC) terms, KEGG pathways, and GSEA analyses with adjusted P value or a false discovery rate (FDR) less than 0.05 were considered statistically significant.

Statistical analysis

The continuous variables were expressed as mean ± standard deviation or median (quartile range), and categorical variables were presented as frequencies (percentages). A Chi-square test was used to compare the differences between independent groups. Cox proportional hazard regression analysis was conducted to evaluate the association of lncRNA signature in predicting overall survival in CCA patients. Kaplan-Meier (KM) analysis was used to determine survival outcomes. The median values were used as a cut-off level to plot the KM curves, and the log-rank test was performed to evaluate the statistical significance. The results of the stepwise multivariate Cox regression analysis of the AIC (Akaike Information Criterion, assessing the goodness of fit of a statistical model) test yielded a predictive model with optimal interpretation and information effectiveness. A linear correlation model was performed to evaluate the relationships between the variables and the Pearson correlation coefficient or Spearman rank correlation coefficient was used to present the result. Unless otherwise indicated, all statistical tests were two-sided and a p-value<0.05 was considered as statistically significant. All data analysis was performed with R (version 3.3.3; http://www.r-project.org/). The differential expression of the lncRNA profile was estimated by the R “edgeR” package. Unsupervised hierarchical clustering analysis was accomplished by the R “pheatmap” package and represented as a volcano plot. KM survival analysis and Cox proportional hazard regression analysis was performed by the R “survival” package. The AUC value was calculated by the R “Survival ROC” package. The R “clusterProfiler”, “pathview” and “venn” package were used to find the common DPCGs of lncRNAs for KEGG pathways.

Data sharing and data accessibility

The data used to support the findings of this study are available from the corresponding author upon request.