Introduction
============

Hepatocellular carcinoma, accounting for 90% of primary liver cancers, is one of the most common malignancies worldwide particularly in East Asia and sub-Saharan Africa. More than 250,000 new cases and approximately 500,000--600,000 deaths occur due to this disease annually ([@ref-11]). Despite the recent advances in early diagnosis, imaging techniques and therapeutic intervention for HCC, low overall 5-year survival rate of HCC patients remains unsatisfactory ([@ref-27]).

LncRNAs, defined as ncRNAs longer than 200 nucleotides, play essential roles in almost every aspect of biological processes. Over the past decades, emerging evidences have associated the misexpression of lncRNAs with tumor initiation, metastasis, and prognosis ([@ref-26]; [@ref-31]; [@ref-34]). Many studies, up to now, have highlighted the molecular mechanism and biological characters of lncRNA in HCC occurrence and progression ([@ref-24]; [@ref-29]). Moreover, some lncRNAs can also be served as valuable prognostic predictor for HCC patients ([@ref-7]; [@ref-21]; [@ref-32]).

With the development of high-throughput technologies, such as RNA-seq and microarray, tremendous expression datasets of gene transcripts have been created for multiple cancer types. TCGA (The Cancer Genome Atlas), a project comprised of 33 cancer types, aiming to improve the prevention, diagnosis, and treatment of cancer ([@ref-6]). In our present study, expression data and related clinical information of 371 LIHC (liver hepatocellular carcinoma) patients were obtained from TCGA, then four lncRNAs were recognized as independent prognostic biomarkers in the survival analysis. In addition, we also investigated the possible biological mechanism uncovering the association of these prognostic lncRNAs with cancer occurrence and progression.

Materials and Methods
=====================

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 ([@ref-6]). 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) ([@ref-20]), 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) ([@ref-16]), 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 ([@ref-35]). 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 ([@ref-30]). 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 ([@ref-1]; [@ref-9]).

Validation for the expression values of lncRNAs
-----------------------------------------------

First, nucleotide sequences of each identified lncRNA were downloaded from GENCODE (release 26) ([@ref-16]). 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 ([@ref-2]) to validate differential lncRNA expression.

Results
=======

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~2~FC \> 1 and adjusted *p*-value \<0.01), the expression patterns of which were displayed in [Fig. 1A](#fig-1){ref-type="fig"}. 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](#fig-1){ref-type="fig"}). 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](#fig-1){ref-type="fig"}). The detailed description of these four identified lncRNAs was shown in [Table 1](#table-1){ref-type="table"}.

![Identification of prognostic lncRNA for patients with HCC.\
(A) Heatmap of 1,074 differentially expressed lncRNAs between normal and tumor tissues. (B) *Z* score of 1,074 deregulated lncRNAs from univariate Cox aggression analysis (brown dots indicate lncRNAs with *p* value \<0.01). (C) *Z* score of four identified prognostic lncRNAs from univariate and multivariate Cox aggression analysis.](peerj-05-3575-g001){#fig-1}

10.7717/peerj.3575/table-1

###### Information of the four identified prognostic lncRNAs for HCC.

![](peerj-05-3575-g007)

  Ensembl ID        Gene name         Chr.   Coordinate               *Z* score[^a^](#table-1fn1){ref-type="fn"}   *P* value[^a^](#table-1fn1){ref-type="fn"}
  ----------------- ----------------- ------ ------------------------ -------------------------------------------- --------------------------------------------
  ENSG00000267583   *RP11-322E11.5*   18     35,443,869--35,467,088   3.386                                        6.935e−04
  ENSG00000254290   *RP11-150O12.3*   8      37,597,480--37,599,858   3.313                                        8.284e−04
  ENSG00000230587   *AC093609.1*      2      43,097,746--43,102,691   −3.595                                       3.143e−04
  ENSG00000264016   *CTC-297N7.9*     17     10,741,267--10,769,016   −5.151                                       2.376e−07

**Notes.**

Statistics derived from univariate Cox proportional hazards regression analysis.

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](#fig-2){ref-type="fig"}, 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.

![Expression patterns of four prognostic lncRNAs.\
(A) Expression pattern of *RP11-322E11.5*. (B) Expression pattern of *RP11-150O12.3*. (C) Expression pattern of *AC093609.1*. (D) Expression pattern of *CTC-297N7.9*. Normal indicates paired adjacent tissue from a subset of 371 primary tumors; Low/High risk subgroups were divided by our risk score model; statistical significance between different subgroups were performed with the two-tailed Wilcoxon test.](peerj-05-3575-g002){#fig-2}

Performance evaluation of risk score model for survival prediction
------------------------------------------------------------------

The Kaplan--Meier curves of the above four lncRNAs were presented in [Fig. S1](#supp-1){ref-type="supplementary-material"}. The distribution of risk score, following-up time and survival status were shown in [Figs. 3A](#fig-3){ref-type="fig"} and [3B](#fig-3){ref-type="fig"}. 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](#fig-3){ref-type="fig"}). 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](#table-2){ref-type="table"}). The comparison between prognostic lncRNAs, mRNAs and combined transcripts based risk score model were shown in [Table S1](#supp-3){ref-type="supplementary-material"}.

![Performance evaluation of risk score model for overall survival prediction.\
(A) Bar plot of risk score for every patient and ordered by the value of risk score. (B) Distribution of following-up time and risk score for every patient. (C) Kaplan--Meier plots of overall survival for patients with high or low risk score. *P* value was calculated with two-sided log-rank test. (D) ROC curves of the risk score based on prognostic lncRNAs, mRNAs and combined transcripts for five-year survival prediction.](peerj-05-3575-g003){#fig-3}

10.7717/peerj.3575/table-2

###### Univariate and multivariate Cox proportional hazards regression analysis.

![](peerj-05-3575-g008)

                                                Univariate analysis     Multivariate analysis                          
  --------------------------------------------- ----------------------- ----------------------- ---------------------- ----------
  Risk score[^a^](#table-2fn1){ref-type="fn"}   2.718 (2.103--3.514)    2.32E−14                2.513 (1.886--3.348)   3.12E−10
  Age                                           1.01 (0.997--1.024)     0.139                   1.009 (0.993--1.025)   0.2602
  Gender (Male/Female)                          0.816 (0.573--1.163)    0.26                    0.968 (0.646--1.450)   0.8748
  TP53 mutation (Y/N)                           1.563 (1.078--2.266)    0.0184                  1.395 (0.915--2.128)   0.1217
  HBV infection (Y/N)                           1.697 (1.199--2.402)    0.0029                  1.656 (1.113--2.463)   0.0129
  Pathologic stage (ii/i)                       1.423 (0.872--2.323)    0.1583                  0.864 (0.513--1.455)   0.5827
  Pathologic stage (iii/i)                      2.676 (1.754--4.083)    4.92E−06                1.794 (1.125--2.861)   0.0141
  Pathologic stage (iv/i)                       5.496 (1.695--17.821)   0.0045                  2.836 (0.818--9.833)   0.1003

**Notes.**

Derived from risk score model based on the expression of prognostic lncRNAs.

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](#fig-3){ref-type="fig"}). 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](#table-2){ref-type="table"}).

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](#fig-4){ref-type="fig"}, *p* \< 0.0001, mean 3.04 years vs. 1.31 years), but not for patients without HBV infection ([Fig. 4B](#fig-4){ref-type="fig"}, *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](#fig-4){ref-type="fig"}, *n* = 256, *p* = 0.00013, mean 2.75 years vs. 1.88 years) and stage iii--iv subgroup ([Fig. 4D](#fig-4){ref-type="fig"}, *n* = 90, *p* = 0.0079, mean 2.39 years vs. 1.50 years).

![Stratification survival analysis for patients with available HBV infection and pathological information.\
(A--B) Stratification analysis of patients with or without HBV infection. (C) Stratification analysis of patients with pathological stage i and ii. (D) Stratification analysis of patients with pathological stage iii and iv. *P* value were all calculated with two-sided log-rank test.](peerj-05-3575-g004){#fig-4}

Expression profiling of identified prognostic lncRNAs
-----------------------------------------------------

After the examinations of HCC associated datasets in GEO, variations of lncRNA expression were summarized in [Table 3](#table-3){ref-type="table"}. 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](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE98269), which was the same as our result (lgFC \> 0).

10.7717/peerj.3575/table-3

###### Expression profiles of 4 lncRNAs in GEO datasets.

![](peerj-05-3575-g009)

  LncRNA          Year   Accession                                                                 Normal   Tumor   lgFC    *P* value
  --------------- ------ ------------------------------------------------------------------------- -------- ------- ------- -----------
  RP11-322E11.5   2017   [GSE98269](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE98269)   3        3       0.41    4.50E--01
  RP11-150O12.3   2008   [GSE6222](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE6222)     2        10      1.41    1.60E--01
                  2013   [GSE50579](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50579)   10       67      1.93    6.05E--03
                  2014   [GSE55092](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE55092)   91       49      0.28    3.78E--06
                  2014   [GSE58043](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE58043)   7        7       1.29    7.16E--02
                  2015   [GSE67260](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE67260)   5        5       3.17    3.09E--03
                  2016   [GSE75271](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE75271)   5        50      0.47    2.30E--02
  AC093609.1      2008   [GSE6222](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=)            2        10      −1.88   4.15E--02
                  2011   [GSE29721](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE29721)   10       10      −0.49   1.33E--03
                  2013   [GSE17548](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE17548)   20       17      −0.33   4.69E--04
                  2013   [GSE47197](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE47197)   63       61      −0.6    1.22E--12
                  2013   [GSE50579](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50579)   10       67      −1.17   5.68E--05
                  2014   [GSE62232](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62232)   10       81      −0.48   2.74E--06
                  2014   [GSE55092](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE55092)   91       49      −0.54   4.79E--20
                  2014   [GSE58043](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE58043)   7        7       −1.48   4.48E--02
                  2015   [GSE67260](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE67260)   5        5       −2.06   6.63E--03
                  2017   [GSE70880](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE70880)   16       16      −1.48   2.71E--06
                  2017   [GSE98269](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE98269)   3        3       −1.72   1.11E--02
  CTC-297N7.9     2017   [GSE98269](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE98269)   3        3       −1.4    2.48E−02
                  2015   [GSE67260](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE67260)   5        5       −2.24   4.44E--04

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](#fig-5){ref-type="fig"}--[5B](#fig-5){ref-type="fig"}, [Fig. S2A](#supp-2){ref-type="supplementary-material"}). For 1,321 down-regulated mRNAs, it revealed that these mRNAs are participant in many pathways or diseases ([Figs. 5C](#fig-5){ref-type="fig"}--[5D](#fig-5){ref-type="fig"}, [Fig. S2B](#supp-2){ref-type="supplementary-material"}), such as blood microparticle, biosynthesis of basic nutrients, drug metabolism and multiple metabolic diseases.

![KEGG and DO enrichment analysis of mRNAs deregulated in high risk group compared with low risk group.\
(A--B) KEGG and DO analysis of up-regulated mRNAs. (C--D) KEGG and DO analysis of down-regulated mRNAs.](peerj-05-3575-g005){#fig-5}

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](#fig-6){ref-type="fig"}). 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](#fig-6){ref-type="fig"}, 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](#fig-6){ref-type="fig"}, 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.

![Biological features of prognostic lncRNA *CTC-297N7.9*.\
(A) Relative genomic location of *CTC-297N7.9* and *TMEM220*. (B) Expression correlation analysis of *CTC-297N7.9* and *TMEM220* in patients with HCC. (C) Expression patterns of *TMEM220* in normal tissues and cancer cell lines. Relative expression was calculated with average expression level of all normal tissues or cancer cell lines as one, respectively.](peerj-05-3575-g006){#fig-6}

Discussion
==========

Nowadays, much is known about the cellular changes and etiological agents responsible for HCC, while the underlying molecular pathogenesis of HCC is still unclear. Although much progress has been made, many issues still remain unresolved, such as the prognosis evaluation of HCC cases ([@ref-22]). Traditional prognostic staging models, including BCLC, Okuda and CLIP scoring systems, have been used to stage HCC ([@ref-3]; [@ref-15]). Mounts of novel prognostic biomarkers based on gene expression also have been carried out to predict HCC prognosis ([@ref-4]; [@ref-14]; [@ref-17]). However, neither of which can distinguish every situation. Up to now, mounts of reports have regarded lncRNAs as oncogenes or tumor suppressors. Previous studies on lncRNAs in HCC mainly focused on well-known cancer related transcripts. For example, overexpression of *HOTAIR* and *MALAT-1* were found to be associated with tumor recurrence of HCC after liver transplantation ([@ref-19]; [@ref-29]). In addition, suppression of these two lncRNAs in HepG2 (liver cancer cell line) could reduce cell viability, invasiveness, and increase the sensitivity to apoptosis.

In the present study, we found 1,074 dysregulated lncRNAs in LIHC, some of which have been recognized to be associated with cancer. For example, *LINC01419* was reported to be significantly overexpressed in HBV-related and HCV-related HCC ([@ref-33]), which was consistent with our result. A single nucleotide polymorphism lies in *HAGLR* was identified as a risk site for ovarian cancer ([@ref-5]). Expression of *HOTTIP* is associated with HCC patients' clinical progression and survival outcome ([@ref-24]). *LINC01093* and *FAM99A*, both specific expression in liver tissue, was significantly down-regulated in HCC compared with normal livers, as well as lower in cirrhotic tissues relative to normal tissues ([@ref-10]). *LINC00844* expression was also down-regulated in HCC and breast cancer ([@ref-12]; [@ref-18]).

Through univariate and multivariate survival analysis, four lncRNAs were identified to have the ability to predict patients' overall survival independently. We then developed a risk score model based on expression of the four prognostic lncRNAs. With this risk model, patients were divided into LRG or HRG, which exhibited strong predictive power for patients' outcome. Patients in LRG are more likely to survival longer than those in HRG. Moreover, prediction power could be enhanced when fitted lncRNA in the model based on mRNA signatures, which also indicate the significance of lncRNA as a valuable supplement to the regular mRNA markers. In addition, other prognostic factors for HCC, like HBV infection status and pathological stage, also showed significant correlation with survival expectations in our multivariate Cox regression analysis. Therefore, we performed stratification analysis for above clinical factors; the results suggested the ability of our risk score model to distinguish patients with significantly different prognostic outcomes (except for patients without HBV infection subgroup).

To confirm our results, we first explored the expression profiles of above lncRNAs in other datasets, which are similar with our findings. We then search some databases to find appropriate cohort to make a cross-validation of our survival prediction model. Considering the fact that many microarrays are not designed for lncRNA, we didn't find suitable dataset with both expression profiles of the four lncRNAs and detailed clinical data for survival analysis in GEO database. In [@ref-13]'s [@ref-13] study, RNAseq analysis were performed on 254 liver tumors, which may be helpful for our validation. However, the original sequencing data deposited in European Gnome-phenome Archive (EGA) database were not publicly accessible, and the processed expression matrix were also not supplied by the author. As a result, we only make a partial cross-validation of the expression profiles of these four prognostic lncRNAs.

As the risk score model was constructed on the basis of lncRNAs expression value and its powerful discrimination of expected survival, differentially expressed mRNAs between LRG and HRG were considered to share similar expression patterns with these signature lncRNAs. On the other hand, these differentially expressed mRNAs can also reflect the possible cause of the significantly different prognosis between these two groups. Enrichment analysis revealed that up-regulated genes are mainly involved in cell cycle pathway and associated with multiple cancers. Those down-regulated genes, linked with a variety of metabolic or cardiovascular diseases, are mainly enriched in drug metabolism pathway, which are consistent with the poor prognosis in HRG.

Although accumulating studies have demonstrated that lncRNAs involve in diverse cancers, the underlying mechanisms are still not well understood. In [@ref-23]'s [@ref-23] study, lncRNA *HULC*, acting as a post-transcriptional regulator of gene expression, was detected to be highly specific up-regulation in HCC. Moreover, subsequent research has shown that HULC can trigger protective autophagy and attenuate the sensitivity of HCC cells to chemotherapeutic agents ([@ref-28]). Among the four prognostic lncRNAs, *RP11-150O12.3* was recognized as an independent predictor of gastric cancer (GC) prognosis ([@ref-25]), which also related to survival time of patients with colorectal adenocarcinoma ([@ref-31]). *CTC-297N7.9*, locating upstream of protein coding gene *TMEM220*, showed its protective role in tumorigenesis and progression. As a transmembrane protein, *TMEM220* is specifically highly expressed in liver tissue compared with other normal tissues. Recently, a study find that *TMEM220* was down-regulated and highly methylated in GC tissues compared to normal tissues. The demethylation experiment also implied that the methylation of *TMEM220* was inversely correlated with its expression ([@ref-8]). Therefore, we suggest that *CTC-297N7.9* may be able to regulate the methylation of *TMEM220*, or be involved in cell autophagy through the interaction with functional proteins, and then affect the prognosis of HCC patients.

Conclusions
===========

In summary, we recognized four lncRNAs associated with OS of HCC patients, and constructed a risk score model based on lncRNAs expression profiles to make a survival prediction. The four lncRNAs could independently predict OS both in univariate and multivariate survival analysis. Our risk score model had a strong power to distinguish patients with different survival outcomes. These results shown the possibility of lncRNAs as novel biomarkers and therapy targets for HCC. Moreover, this study also revealed the possible role of lncRNAs in occurrence and progression of liver cancer.

Supplemental Information
========================

10.7717/peerj.3575/supp-1

###### Survival analysis of four independent prognostic lncRNA markers

Kaplan-Meier survival curve for (A) *RP11-322E11.5* (B) *RP11-150O12.3* (C) *AC093609.1* and (D) *CTC-297N7.9.* Patients were divided into high or low expression group with median expression value as the cutoff.

###### 

Click here for additional data file.

10.7717/peerj.3575/supp-2

###### GO enrichment analysis of mRNAs deregulated in high risk group compared with low risk group

\(A\) GO analysis of upregulated mRNAs. (B) GO analysis of downregulated mRNAs.

###### 

Click here for additional data file.

10.7717/peerj.3575/supp-3

###### Comparison between lncRNAs, mRNAs and combined transcripts based risk score model in the overall and stratification survival analysis

###### 

Click here for additional data file.

10.7717/peerj.3575/supp-4

###### The corresponding probe name of four lncRNAs among different microarray platforms

###### 

Click here for additional data file.

HCC

:   Hepatocellular carcinoma

ROC

:   Receiver operating characteristic curve

HBV

:   Hepatitis B virus

HCV

:   Hepatitis C virus

lncRNA

:   Long non-coding RNA

TCGA

:   The Cancer Genome Atlas

LIHC

:   Liver hepatocellular carcinoma

TANRIC

:   The Atlas of Noncoding RNAs in Cancer

GO

:   Gene ontology

KEGG

:   Kyoto Encyclopedia of Genes and Genomes

DO

:   Disease ontology

AUC

:   Area under the ROC curve

GTEx

:   Genotype-Tissue Expression

CCLE

:   Cancer Cell Line Encyclopedia

OS

:   Overall survival

HRG

:   High risk group

LRG

:   Low risk group

BCLC

:   Barcelona Clinic Liver Cancer

CLIP

:   The Cancer of the Liver Italian Program

GC

:   Gastric cancer

Additional Information and Declarations
=======================================

The authors declare there are no competing interests.

[Zhonghao Wang](#author-1){ref-type="contrib"} conceived and designed the experiments, performed the experiments, analyzed the data, contributed reagents/materials/analysis tools, prepared figures and/or tables.

[Qian Wu](#author-2){ref-type="contrib"} performed the experiments, analyzed the data, wrote the paper, reviewed drafts of the paper.

[Shu Feng](#author-3){ref-type="contrib"} performed the experiments, analyzed the data, wrote the paper.

[Yanhua Zhao](#author-4){ref-type="contrib"} analyzed the data, contributed reagents/materials/analysis tools.

[Chuanmin Tao](#author-5){ref-type="contrib"} conceived and designed the experiments, reviewed drafts of the paper.

The following information was supplied regarding data availability:

The origin RNAseq data used in our study were all downloaded from the TCGA data portal: <https://portal.gdc.cancer.gov/>.
