Identification of four prognostic LncRNAs for survival prediction of patients with hepatocellular carcinoma

Collection of TCGA dataset and related clinical information

The RNA-Seq data of patients in TCGA-LIHC project were downloaded from Genomic Data Commons Data Portal, as well as the corresponding clinical information (Cancer Genome Atlas Research, 2011). To obtain lncRNA expression profile, we first got ensembl ID of lncRNAs which had been filtered by The Atlas of Noncoding RNAs in Cancer (TANRIC) (Li et al., 2015), then 12,309 lncRNA ID were used to retrieve the required expression information from above RNA-Seq data. Similarly, 19,817 mRNAs ID were downloaded from GENCODE (release 26) (Harrow et al., 2012), then we got the expression profiles of mRNAs from the RNA-Seq data.

Identification of potential prognostic signatures

After removing lncRNAs whose expression is zero in more than 50% of the samples, a total of 4,683 lncRNAs were utilized to do the differential expression analysis among normal and LIHC. To find the prognostic lncRNAs, these patients were analyzed by univariate and multivariate Cox proportional hazards regression methods. Briefly, univariate survival analysis was performed with expression value of lncRNAs as independent variable, while multivariate analysis was conducted with expression value of lncRNAs and other factors (including clinical features or other lncRNAs) as the independent variables. To construct the predictive risk score model, we only used those lncRNAs which were still independent of other factors after the multivariate Cox regression analysis. The coefficient of each prognostic lncRNA in risk score model was derived from corresponding multivariate analysis of these identified lncRNA markers (Zhou et al., 2015). All patients were divided into two groups according to the risk score. Then, for mRNAs differentially expressed between these two group, prognostic mRNAs were recognized with the same method of lncRNA, as well as the mRNA-based model. The combined transcripts based model was fitted with both above identified lncRNA and mRNA biomarkers.

Statistical analysis

All bioinformatics analysis and statistical tests were performed with R language (version 3.3.2): (i) limma package for differential expression analysis; (ii) survival package for univariate and multivariate Cox proportional hazards regression survival analysis; (iii) pROC package for ROC curve and the area under the ROC curves (AUC) analysis; (iv) Gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), Disease Ontology (DO) enrichment analysis were performed with clusterProfiler package (Yu et al., 2012). The co-expression correlation between CTC-297N7.9 and TMEM220 were evaluated with Pearson correlation coefficient. The expression profiles of TMEM220 in normal tissues and cancer cell lines were retrieved from Genotype-Tissue Expression (GTEx) project and Cancer Cell Line Encyclopedia (CCLE) project, respectively (Barretina et al., 2012; Consortium GT, 2013).

Validation for the expression values of lncRNAs

First, nucleotide sequences of each identified lncRNA were downloaded from GENCODE (release 26) (Harrow et al., 2012). The probe sequences of microarray were retrieved from matching platform annotation files deposited in Gene Expression Omnibus (GEO) database with R package GEOquery. Next, probes were re-mapped to lncRNA sequences using Magic-BLAST tool (version 1.2.0, from NCBI), and probes which uniquely mapped to the sequence with no mismatch were retained to represent the corresponding lncRNA. For each GEO dataset which focused on HCC, we used GEO2R web application supplied by GEO (Barrett et al., 2013) to validate differential lncRNA expression.

Differential expression profile and identification of prognostic lncRNAs

From TCGA, RNA-Seq expression profiles of LIHC were downloaded, which covered 371 primary tumor samples and 50 (a subset of 371 HCC patients) paired adjacent normal liver tissue samples. After the initial analysis of all lncRNAs, a total of 1,074 lncRNAs were considered to be differentially expressed between tumor and normal (cutoff of absolute log₂FC > 1 and adjusted p-value <0.01), the expression patterns of which were displayed in Fig. 1A. Then we conducted a univariate Cox proportional hazards regression analysis of the above 1,074 lncRNAs, which revealed a subset of 108 lncRNAs significantly associated with patients’ OS (overall survival, p-value <0.01, Fig. 1B). Given that OS can be affected by many factors, a multivariate Cox proportional hazards regression analysis was performed with lncRNAs and other clinical features (such as age of initial diagnosis, pathologic stage and TP53 mutation status) to further evaluate the independent prognostic value of these lncRNAs. Finally, four of the 108 lncRNAs were demonstrated as independent predictors of overall survival (Fig. 1C). The detailed description of these four identified lncRNAs was shown in Table 1.

Risk score model based on expression level of independent prognostic lncRNAs

Four lncRNAs screened from above analysis were used to construct a risk score model, risk score = (0.2567 * expression level of RP11-322E11.5) + (0.1307 * expression level of RP11-150O12.3) + (−0.2320 * expression level of AC093609.1) + (−0.1857 * expression level of CTC-297N7.9), then the risk score of each patient was calculated according to this model. With the median value of risk score as cutoff, patients in TCGA datasets were divided into high risk group (HRG, n = 186) and low risk group (LRG, n = 185). As shown in Fig. 2, the four lncRNAs displayed significantly different expression patterns between normal and tumor, LRG and HRG, which indicated possible protective roles of AC093609.1 and CTC-297N7.9, carcinogenic effects of RP11-322E11.5 and RP11-150O12.3 in HCC tumorigenesis and prognosis.

Performance evaluation of risk score model for survival prediction

The Kaplan–Meier curves of the above four lncRNAs were presented in Fig. S1. The distribution of risk score, following-up time and survival status were shown in Figs. 3A and 3B. In the survival analysis, K–M plots exhibited a significant difference in patients’ OS between LRG and HRG (log rank p < 0.0001, Fig. 3C). Patients in LRG had significantly longer OS (mean 2.64 years) when compared with the HRG (mean 1.76 years). The univariate Cox proportional hazards regression analysis implied the significant association between signature lncRNAs based risk score and patients’ OS (HR = 2.718, 95% CI [2.103–3.514], p = 2.32e−14, Table 2). The comparison between prognostic lncRNAs, mRNAs and combined transcripts based risk score model were shown in Table S1.

Furthermore, ROC analysis were also performed to evaluate specificity and sensitivity of the risk score models for five-year survival prediction. The AUC derived from lncRNA risk score model was 0.701, suggesting that our risk score was a good indicator for survival prediction (Fig. 3D). Additionally, risk model developed with the combination of lncRNAs and mRNAs exhibited the best prediction power when compared with only lncRNAs or mRNAs (AUC = 0.73 vs 0.70 or 0.68, respectively). To evaluate the independent predictive performance of our model when considering the traditional clinical prognostic factors, a multivariate Cox regression analysis was conducted with risk score and other clinical features as independent variables in this TCGA dataset. The results revealed a significant correlation with survival time after adjusting for multiple clinical variates (HR = 2.513, 95% CI [1.886–3.348], p = 3.12e−10, Table 2).

Stratified Cox proportional hazards regression analysis

Moreover, the results also indicated HBV infection status and pathological stage iii as independent prognostic factors in the above multivariate analysis (p = 0.0129 and p = 0.0141, respectively), so that we performed stratification analysis for HBV infection status and pathological stage next. All patients were first stratified into group with (n = 140) or without HBV infection (n = 230). For patients with HBV infection, significant longer OS can still be observed in LRG than HRG (Fig. 4A, p < 0.0001, mean 3.04 years vs. 1.31 years), but not for patients without HBV infection (Fig. 4B, p = 0.18, mean 2.45 years vs. 2.11 years). Similarly, when 370 patients stratified according to different pathological stage, the Cox regression analysis also exhibited a strong classification power between LRG and HRG both in stage i–ii (Fig. 4C, n = 256, p = 0.00013, mean 2.75 years vs. 1.88 years) and stage iii–iv subgroup (Fig. 4D, n = 90, p = 0.0079, mean 2.39 years vs. 1.50 years).

Expression profiling of identified prognostic lncRNAs

After the examinations of HCC associated datasets in GEO, variations of lncRNA expression were summarized in Table 3. For AC093609.1, CTC-297N7.9 and RP11-150O12.3, many other studies manifested the same significant difference as our study in liver tumor tissue (p < 0.05). Despite no evidence of significant difference was seen for RP11-322E11.5 (p > 0.05), the expression was also up-regulated in GSE98269, which was the same as our result (lgFC > 0).

Functional implication of prognostic lncRNA signatures

To further explore the possible role of these four lncRNA in HCC tumorigenesis and prognosis, we performed enrichment analysis with GO, KEGG and DO databases. First we carried out the expression analysis of mRNA between LRG and HRG, 2,510 in total were recognized as differentially expressed genes. Then we conducted functional enrichment analysis for up-regulated and down-regulated mRNA respectively. The GO, KEGG and DO enrichment showed that 1,189 up-regulated mRNAs are mainly involved in chromosome function, cell cycle, p53 signal pathway and associated with a variety of cancers (Figs. 5A–5B, Fig. S2A). For 1,321 down-regulated mRNAs, it revealed that these mRNAs are participant in many pathways or diseases (Figs. 5C–5D, Fig. S2B), such as blood microparticle, biosynthesis of basic nutrients, drug metabolism and multiple metabolic diseases.

Biological significance of the prognostic lncRNA CTC-297N7.9

We further examined whether there were any protein-coding genes which are neighborhoods of these lncRNAs in genomic location coordinates. Among these four prognostic lncRNAs, CTC-297N7.9 locates upstream of TMEM220, which is a transmembrane protein-coding gene (Fig. 6A). Then expression correlation of CTC-297N7.9 with TMEM220 were evaluated by Pearson correlation analysis, which implied a strong positive correlation between this paired lncRNA and mRNA (Fig. 6B, correlation coefficient = 0.795, p < 0.0001). Moreover, to explore the possible role of TMEM220 in HCC, we retrieved expression profile of this gene from several databases. From Fig. 6C, we can see a liver-specific expression of TMEM220 among different types of normal tissue. However, this tissue-specific expression pattern disappeared in liver cancer cell when compared with multiple tissue-derived cancer cell lines.