Identification of potential prognostic long non-coding RNA for predicting survival in intrahepatic cholangiocarcinoma.

Intrahepatic cholangiocarcinoma (ICC) is an aggressive biliary epithelial tumor with poor prognosis. There are increasing evidences that long non-coding RNAs (lncRNAs) are dysregulated in multifarious tumors, revealing potential significant role of lncRNAs in tumorigenesis.We used the ICC dataset retrieved from The Cancer Genome Atlas and the Gene Expression Omnibus database to obtain the lncRNAs expression profiles and identify potential prognostic lncRNAs for predicting the prognosis in ICC. Univariate and multivariate Cox regression analyses were performed to construct a prognostic index (PI). Furthermore, coexpression analysis and functional assessment were performed to initially investigate the function of these prognostic lncRNAs.A total of 255 differentially expressed lncRNAs (DElncRNAs) were identified among two RNA sequencing dataset of a total 63 ICC patients with 98 samples using R platform. Thirteen of 255 DElncRNAs were identified as prognostic lncRNAs and used for a PI. Patients with high PI were associated with poor prognostic (P = .0064), and the Cox regression showed consistent result (P = .042). The time-dependent receiver operating characteristic analysis showed the PI performed well in ICC survival prediction with an area under curve of 0.921, 0.801, and 0.717 for 1-, 3-, and 5-year survival, respectively.In conclusion, we included 13 identified prognostic DElncRNAs and constructed a prognostic signature/PI. ICC patient with higher PI was associated with poorer prognosis. However, the clinical role as well as biological functions of constructed PI and these prognostic DElncRNAs need to be verified in future study.

Introduction

Intrahepatic cholangiocarcinoma (ICC), which is derived above the second-order bile ducts, is an aggressive biliary epithelial tumor with poor prognosis. The incidence and patient mortality of ICC have increased globally over the past few decades.[ Multiple factors, including genetic and environmental, participate in hepatocarcinogenesis of ICC. Chronic hepatitis B and C, hepatobiliary flukes, primary sclerosing cholangitis, biliary tract cysts, hepatolithiasis and toxins are considered as the main environmental factors associated with ICC.[ With the development of high-throughput technologies enabling to study through genome-wide screening, the primary genetic factors of ICC carcinogenesis were found including IDH1, IDH2, FGFR1, FGFR2, FGFR3 and KRAS/MAPK.[ However, it is little known that which and how these biomarkers affecting prognosis of ICC. Meanwhile, more comprehensive and new prognostic signatures, as well as potential treatment targets, are urged to be discovered toward ICC.

With advances in technologies and years of study, non-coding RNAs including long non-coding RNAs (lncRNAs) are gradually found being greatly significant in many biological and pathological processes though regulating transcriptional and translational output instead of protein-coding capacity like mRNAs.[ Furthermore, there are increasing evidences that some lncRNAs are dysregulated in multifarious tumors, revealing potential significant role of lncRNAs in tumorigenesis.[ Recent studies have shown that lncRNAs played an important role as prognostic factors and biomarkers in hepatocellular caricinoma (HCC).[ And the mechanisms and potential therapeutics targets were well elucidated.[ However, the relative lncRNAs and associated mRNAs of ICC are still unknown and discovery for new lncRNAs involved in ICC is urgently needed.

In the present study, we intended to use the ICC dataset retrieved from The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) database to obtain the lncRNAs expression profiles and identify potential prognostic lncRNAs for predicting the prognosis in ICC. These findings could provide new insights into the molecular mechanisms for ICC.

Materials and methods

Data source

Two RNA sequencing (RNA-Seq) datasets of a total 63 ICC patients with 98 samples were retrieved from TCGA data portal and GEO database (GSE107943) in March 2019, including 33 ICC tissues with 8 normal liver tissues in the former and 30 with 27 in the latter, respectively. The relative clinical data of ICC patients were also retrieved from TCGA while these information from GSE107943 were not accessible. Patients with distal cholangiocarcinoma were excluded in the study. In addition, extra ethics approval was not essential for the data were all obtained from TCGA and GEO database. The authors cannot access to information that could identify individual participants during or after data collection.

Screening of differentially expressed lncRNAs (DElncRNAs) and differentially expressed mRNAs (DEmRNAs)

The DElncRNAs and DEmRNAs were calculated by R platform using both edgeR[ and DESeq2[ packages. In the comparison between tumor tissues and normal control, false discover rate or adjust value < .05, and |log2 foldchange (log2FC)| > = 1 were identified as differentially expressed. The overlapping DElncRNAs and DEmRNAs between edgeR and DESeq2 in both datasets were used for further analyses. Hierarchial clustering analyses were conducted by R platform.

Identification of prognostic DElncRNAs and construction of prognostic index based on prognostic DElncRNAs

Toward 33 ICC patients containing survival data, univariate Cox regression analysis for overall survival (OS) was performed to identify prognostic DElncRNAs. DElncRNAs with P value < .05 were considered as prognostic and used for further analyses. For constructing a prognostic signature, a multivariate Cox analysis with OS was performed to determine regression coefficient (β) as the weight for each prognostic DElncRNA. The formula of prognostic index (PI) was as followed: prognostic index = expression of DElncRNA1 × β1DElncRNA1 + expression of DElncRNA2 × β2DElncRNA2 +…expression of DElncRNAn × βnDElncRNAn.[ According to PI for each ICC patient, they would be divided into low and high risk group by median PI. In the meantime, survival analysis and time-dependent receiver operating characteristic (ROC) curve were created in R platform with survivalROC package for evaluating the effectiveness of the PI in predicting OS among ICC patients.[ Besides, univariate and multivariate Cox regression analyses were also performed among all the accessible patient characteristics to reduce bias and thus confirm the effect of prognostic signature.

Coexpression analysis and functional assessment

Correlated pairs of DElncRNAs and DEmRNAs were identified by Pearson correlation coefficients (PPCs) based on their expression levels. DEmRNAs with Pearson correlation coefficient > 0.9 were considered as prognostic DElncRNA-correlated DEmRNAs, and used for construction of coexpression network and following functional assessment. Cytoscape software was applied to establish the DElncRNAs and DEmRNAs coexpression network.[

The function of lncRNAs is usually associated with relative mRNAs instead of encoding proteins themselves.[ Thus, functional enrichment analysis of Gene Ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways for co-expressed DEmRNAs was performed to infer biological function of prognostic DElncRNAs through the Database for Annotation, Visualization, and Integrated Discovery version 6.8 (DAVID v6.8).[P value < .05 was considered as significantly enriched functional annotations in GO terms and KEGG pathways.

Statistical analysis

All statistical analyses were performed using the Statistical Package for Social Sciences 22.0 for Windows (SPSS Inc., Chicago, IL) and R 3.3.0. Associated factors of OS were identified using univariate Cox proportional hazard regression while those patient characteristics with P value < .05 were included in the multivariate Cox proportional hazards regression model for adjustment. All statistical assessments were two-tailed, and P value < .05 was considered statistically significant.

Results

Identification of DElncRNAs

For the expression data of 33 ICC tissues and 8 normal liver tissues in TCGA, 2438 lncRNAs were identified as DElncRNAs using DESeq2 (Fig. 1A), while 2730 lncRNAs using edgeR (Fig. 1B). For the expression data of 30 ICC tissues and 27 normal liver tissues in GSE107943, 3084 lncRNAs were identified as DElncRNAs using DESeq2 (Fig. 1C), while 4544 using edgeR (Fig. 1D). Among them, 255 DElncRNAs overlapped in the 4 DElncRNAs series were used for further analyses (Fig. 1E). The heatmap of the 255 DElncRNAs in two databases were shown in Figure 2.

Figure 1

Analysis of differentially expressed lncRNAs (DElncRNAs). (A) DElncRNAs of patients in TCGA identified using DESeq2 package; (B) DElncRNAs of patients in TCGA identified using edgeR package; (C) DElncRNAs of patients in GSE107943 identified using DESeq2 package; (D) DElncRNAs of patients in GSE107943 identified using edgeR package; (E) Overlapping DElncRNAs.

Figure 2

(A) Heatmap of 255 DElncRNAs in ICC patients in TCGA; (B) Heatmap of 255 DElncRNAs in ICC patients in GSE107943.

Construction of prognostic DElncRNAs expression based prognostic signature

Univariate Cox regression analyses were performed between DElncRNAs and the OS of 33 patients in TCGA to identify prognostic DElncRNAs. Among the 255 obtained DElncRNAs, 13 of them were identified and the information of 13 prognostic DElncRNAs was shown in Table 1. A multivariate Cox regression analysis was subsequently performed to determine the contribution for each prognostic DElncRNA to the signature and the results were also shown in Table 1. The prognostic signature was finally expressed as follow: prognostic index = expression of SFTA1P × (0.707716) + expression of LINC01117 × (0.016705) + expression of LINC00702 × (−0.058047) + expression of LINC01116 × (0.006796) + expression of LINC00632 × (0.031144) + expression of DIRC3 × (0.058592) + expression of UPP2-IT1 × (0.397472) + expression of LOXL1-AS1 × (0.006703)+ expression of LINC00944 × (−0.022045) + expression of LDLRAD4-AS1 × (0.007969) + expression of LINC00519 × (0.057222) + expression of SERTAD4-AS1 × (0.008442) + expression of AC012613.2 × (−0.26531). Patients were then divided into low risk group and high risk group according to medium PI (Fig. 3). ROC curve showed the PI could helpfully predict 1-year, 3-year and 5-year survival of ICC patients, and the area under curve of survival ROC reached 0.921, 0.801, and 0.717, respectively (Fig. 4A). In the meantime, survival analysis showed patients in low risk group had significantly longer OS than patients in high risk group (Log-rank test, P = .0064), indicating low PI could predict a significantly better prognosis (Fig. 4B). Furthermore, univariate and multivariate regression analyses including varies patient characteristics were performed. Although pathological stage was not showed as a risk factor of OS in univariate analysis, we considered that it generally affected the prognosis of ICC and still put it in multivariate analysis. The results were showed in Table 2 and perineural invasion (P = .014) and PI (P = .042) were the independent risk factors of OS in ICC patients.

Table 1

The information of 13 prognostic DElncRNAs and prognostic analysis in patients with ICC.

Figure 3

PI analysis of ICC patients. (A) The low and high score group for the prognostic DElncRNAs signature in ICC patients; (B) The survival status and duration of ICC cases; (C) Heatmap of 13 prognostic DElncRNAs expression in ICC. The color from blue to red shows a trend from low expression to high expression.

Figure 4

ROC and Kaplan-Meier curves for PI in TCGA ICC cohort. (A) Time-dependent ROC curves analysis for survival prediction by PI. (B) Kaplan-Meier survival curves showing overall survival outcomes according to relative high-risk and low-risk patients (Log-rank test, P = .0064).

Table 2

Univariate and multivariate analysis for overall survival.

A total of 4053 DEmRNAs were identify using the same methods for DElncRNAs (Fig. 5). With Pearson correlation coefficient > 0.9, a total of 220 coexpression pairs including 133 DEmRNAs and 6 DElncRNAs for ICC were obtained. UPP2-IT1, LINC01116, LINC00702, LINC00632, DIRC3, and LINC01117 were coexpressed with 102, 31, 30, 27, 18, and 12 DEmRNAs, respectively (Fig. 6). In addition, all the coexpressed DEmRNAs included in network were found to be upregulated.

Figure 5

Venn plot of overlapped DEmRNAs. A. identified using DESeq2 package for patients in TCGA; B. identified using edgeR package for patients in TCGA; C. identified using DESeq2 package for patients in GSE107943; D. identified using edgeR package for patients in GSE107943.

Figure 6

DElncRNA-DEmRNA coexpression network. The blue ellipses and red rhombuses were represented the DEmRNAs and DElncRNAs, respectively.

Based on the coexpressed DEmRNAs, we performed functional enrichment analyses investigating potential functions of the 6 optimal DElncRNAs. GO terms and KEGG pathways with ten minimized P value were respectively showed in Figure 7. Extracellular space (P value = 1.18E−10), Extracellular region (P value = 1.48E−07) and Blood microparticle (P value = 1.67E−05) were three most significantly enriched GO terms while Complement and coagulation cascades (P value = 4.47E−08) with 10 DEmRNAs, Biosynthesis of amino acids (P value = 1.15E-05) with 8 DEmRNAs, Metabolic pathways (P value = 4.05E−05) with 29 DEmRNAs were the 3 most significantly enriched KEGG pathways.

Figure 7

Function assessment of GO terms and KEGG pathway.

Discussion

ICC, which is a primary liver cancer with incidence only second to the HCC, usually has a poor prognosis even after undergoing radical therapy comparing to HCC.[ Completeness of resection, max size of tumor, number of tumors, vascular invasion, visceral invasion, and lymph node metastases are the widely accepted prognostic factors of ICC.[ Currently, MTHFR, TS, GST01, MRP2/ABC2, FIC1, NKG2D, XRCC1, PTGS2, and COX-2 were considered to be the molecular pathogenesis of ICC involving multiple signal pathways.[ However, no specific biomarker has been reported to effectively predict prognosis of ICC. Thus, it is a great need with significant importance to investigate prognostic biomarkers and signatures of ICC.

In this study, we comprehensively performed an overall analysis of lncRNA and mRNA expression of ICC. Moreover, we identified 13 differentially expressed lncRNA associated with carcinogenesis and development of ICC. Then we constructed a prognosis signature/PI of 13 DElncRNAs that was able to effectively predict survival of ICC patients and was verified as a significantly independent risk factor of OS. For now, we were the first that comprehensively analyzed lncRNA profile using both TCGA and GEO database and we hoped that our study would improve the understanding in carcinogenesis and development of ICC.

Comparing with mRNAs and micro RNAs (miRNAs), it seems like lncRNAs are expressed in a more cell-type- and tissue-specific way, thus having greater advantages as diagnostic and prognostic biomarkers.[ Wang et al[ conducted a study involving 77 ICC patients to investigate mRNAS associated with carcinogenesis and prognosis of ICC and the results showed that CYP2D6, CYP2D7 and PCSK6 were associated with overall survival in ICC. Ma et al[ identified CPS1 and its lncRNA CPS1-IT1 as risk factors associated with poor liver function and reduced survival rates in ICC. Moreover, Zhang et al[ identified CCAT1 functions as an oncogenic lncRNA in ICC through inhibiting miR-152. Recently, the results of study by Lv et al[ showed that EMP1–008, ATF3–008 and RCOR3–13 were observed significantly downregulated in ICC with tumor metastasis. These findings suggest lncRNAs may play an important role in carcinogenesis, development and even metastasis of ICC, indicating the potential value of lncRNAs to be a diagnostic and prognostic biomarker for ICC patients. However, none of them used comprehensive databases with expression profiling by high throughput sequencing.

For the 13 prognostic DElncRNAs finally identified, all of them were upregulated in ICC. And three of them were considered as tumor suppressors while the rest as tumor promoters according to their coefficient in PI. Moreover, 5 of 13 were focused on by previous studies while we knew little about the rest. A study involving 68 pairs of gastric cancer and normal tissues revealed SFTA1P functioned as a tumor suppressor in gastric cancer,[ which was inconsistent with our result where SFTA1P functioned as a tumor promoter. However, a bioinformatic analysis study suggested SFTA1P played an important role in the regulation of both oncogene and tumor supressor gene during the carcinogenesis of lung squamous cell carcinoma,[ revealing complicated function of SFTA1P. Zhang et al[ reported that LINC01116 promoted the proliferation and migration of osteosarcoma cells, which was consistent with our result. There was a controversy about the function of LINC00702, as a tumor promoter[ or a tumor suppressor[. Apparently, our results supported the latter. In addition, DIRC3 was reported once as a prognosis biomarker in laryngeal squamous cell carcinoma.[ For FOXL1-AS1, previous studies revealed that LOXL-AS1 promoted proliferation, thus promoting tumor development in both medulloblastoma and prostate cancer.[

There still were some limitations in our research. Firstly, although comparatively large cohorts came from two databases were used to identified DElncRNAs and DEmRNAs, the limited patients from single TCGA database were included in clinical analyses. Secondly, we did not validate the signature that we developed with extra cohorts. Lastly, we initially predicted the function of prognostic DElncRNAs but exploring further mechanisms and pathways should be focused on in future studies.

Conclusion

In short, we assessed lncRNA data based on RNA-Seq from both TCGA and GEO database, and constructed a prognostic signature including 13 prognostic DElncRNAs. ICC patient with higher PI was associated with poorer prognosis. The effectiveness of signature was proved and the functions of 13 prognostic DElncRNAs were initially investigated. However, the clinical role and the specific mechanisms of constructed PI and the prognostic DElncRNAs should be verified and explored in future studies.

Author contributions

All authors made substantive intellectual contributions to this study to qualify as authors. YH conceived of the design of the study. ZYZ performed the study, collected the data, and contributed to the design of the study. ZYZ prepared the manuscript. ZMW edited the manuscript. All authors read and approved the final manuscript.