Identification of a novel four-lncRNA signature as a prognostic indicator in cirrhotic hepatocellular carcinoma

Data availability statement

The following information was supplied regarding data availability: Raw data were generated in the study. All original data in this study were downloaded from the public databases GEO (https://www.ncbi.nlm.nih.gov/geo/). For the four-lncRNA signature identification, we used GSE14520.

Probe reannotation and preparation of gene expression profiles

The GSE14520 expression profiles and clinical information analyzed in this study were derived from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). To consider the expression and function of lncRNAs in the probes of HCC, the probes of the Affymetrix HT Human Genome U133A Array were reannotated to obtain the lncRNA and mRNA co-expression matrix using the following steps. Firstly, the Affymetrix HT Human Genome U133A Array probe was mapped to the gencode annotation file (Gencode.v29. transcripts, FASTA format,09/27/2018). Information about specific probe set annotation included the probe set ID, Ensemble ID, probe sequence, among other data. Secondly, the probe sets that were allocated to the Ensemble gene IDs in the gencode annotations were acquired. Human coding and noncoding gene annotation files (Homo_sapiens.GRCh38.95.chr.gtf.gz,GTF format,11/25/2018) from the Ensemble database were used to extract the matching information of gene ID and gene symbol. The Ensemble IDs of the probes were then assigned gene IDs, and only those annotated as “protein coding,” “antisense,” “sense intronic,” “bidirectional promoter lncRNA,” “lincRNA,” “non-coding,” “macro lncRNA,” “processed transcript,” “3′ overlapping ncRNA,” and “sense overlapping” were retained. In addition, probes that corresponded to multiple Ensemble IDs were removed. Finally, 12,096 protein-coding transcripts and 610 lncRNA transcripts (different probe IDs may have corresponded to the same transcript) were reannotated.

The GSE14520 probe expression profile from the GEO database was annotated by the reannotated GPL3921 (Affymetrix HT Human Genome U133A Array), including 225 tumor samples and 220 paired non-tumor tissue samples. Among them, the survival information of 203 patients was obtained from the dataset of GSE14520. If multiple probes corresponded to the same gene, the final expression of the gene was determined by the arithmetic mean of multiple probes. The gene expression profile of GSE14520 was normalized by the limma software package. Normalization was performed using the normalizeBetweenArrays function in the limma package (Bolstad et al., 2003). The inclusion criteria were as follows: pathologically confirmed HCC tissues, and patients with liver cirrhosis, and complete follow-up survival information, including overall survival (OS) and recurrence-free survival (RFS). Exclusion criteria included non-tumor tissues, absence of histological examination results, pathological results indicating cholangiocarcinoma, hepatocholangiocarcinoma, or secondary liver cancer, and missing OS or RFS information.

It contained 225 HCC tissue samples and 220 adjacent non-tumor tissue samples. Out of 225 HCC samples, 221 with complete follow-up information (Survival status, Survival time, recurrence status and recurrence—free time) were included in the microarray. Of these 221 patients, 203 have cirrhosis. All 203 participants were patients who underwent radical resection of HCC at the Liver Cancer Institute and Zhongshan Hospital (Fudan University, Shanghai, China) between 2002 and 2003. Of the 203 participants, 177 were male and 26 were female. The mean age of all patients in the dataset was 51.3 years. Regarding the TNM stage of HCC, 153 of the 203 patients were in stage I or II, and 48 were in stage III. Only 2 patients had no documented cancer stage.

Data and clinical characteristics

All 204 patients with liver cancer from the GSE14520 of the GEO database had complete survival information and were included in the study. These patients were randomly divided into the training set (n = 103) and testing set (n = 100) using the R caret software package. The clinical parameters of the two groups of patients are presented in Table 1.

Construction of lncRNA expression-based prognostic signature

Univariate Cox regression analysis was used to screen for significant prognostic lncRNAs in the training set. A P value <0.05 was selected as the threshold for candidate lncRNAs. Candidate prognostic lncRNAs were further reduced by least absolute shrinkage and selection operator (LASSO) regression. The principle of LASSO regression is to eliminate some variables through penalty rules, and ultimately leave all potential predictors with a non-zero coefficient (Gao, Kwan & Shi, 2010). The penalty argument lambda was determined by the cross-validation method using the R glmnet software package. Lambda.min, the lambda value corresponding to the minimum value of the cross-validation error mean, was identified to determine the potential prognostic lncRNAs (Tibshirani, 1997). Each “significant” lncRNA obtained in the above steps were then fitted into a multivariate Cox regression model, and lncRNAs associated with independent prognostic criteria were selected. The prognostic risk formulas were formed, based on the expression level of those remaining eligible lncRNAs multiplied by the multivariate Cox regression coefficient. Risk scores for each sample were calculated using the risk score formula, and patients were divided into high-risk and low-risk groups based on the median risk score of the training cohort as the cutoff value. The same formulas and cutoff values were used to calculate the risk score and risk grouping for patients in the testing set (Zhao et al., 2018). A Kaplan–Meier survival curve and log-rank test were used to analyze the overall survival (OS) and recurrence-free survival (RFS) of patients in the high-risk and low-risk groups in both the training and testing set (Li et al., 2018). A receiver operating characteristic (ROC) curve plotted with the survivalROC software package (version 1.0.3) was adopted to evaluate the specificity and sensitivity of the survival prediction (Heagerty, Lumley & Pepe, 2000). Both ROC and Kaplan–Meier survival curves were constructed with the R studio software (version 3.5.1). A P value < 0.05 was deemed statistically significant.

Functional enrichment analyses

In order to obtain the lncRNA-mRNA co-expression pair, we conducted Pearson correlation analyses between the four lncRNAs identified and the expression profile of the protein-coding genes of the discovery set. The protein-coding genes with a correlation coefficient >0.4 and P < 0.001 were recognized as the four lncRNA-associated genes. The correlated mRNAs were evaluated by functional enrichment analysis to explore the functions of the four prognostic lncRNAs using the clusterProfiler (version 3.10.1) and org.Hs.eg.db (version 3.7.0) package (Yu et al., 2012). Significantly enriched Gene Ontology (GO) terms and Kyoko Encyclopedia of Genes and Genomes (KEGG) pathways with a P value <0.05 and q value <0.05 were visualized using the R studio (version 3.5.1) software.

Statistical analysis

The prognostic ability of the four-lncRNA signature under different clinical parameters was analyzed by Kaplan–Meier survival analysis, to determine which population is most suited for this prognostic model. The survival curve was drawn using the R survminer software package. Statistical analyses were conducted using the SPSS software (version 24.0) and R studio (version 3.5.1).

Construction of risk prognostic scoring system in the training set

The exploration process of this study is shown in Fig. 1. Firstly, 610 lncRNAs were initially screened in the training set (Table S1), and 35 potential OS-related lncRNAs were obtained by a univariate Cox risk regression model (Table S2) (P < 0.05). The remaining lncRNAs were further selected using LASSO regression analysis, and cross-validation was used to select the penalty parameters (Figs. 2A–2B). Five lncRNAs were identified by lambda.min value (Table S3). The lncRNAs obtained in the above steps were inserted into the multivariate Cox regression model. The expression values of four independent prognostic lncRNAs (P value <0.05) and their correlation coefficients in a multivariate regression model were then used to construct prognostic signatures. Detailed information and the significance of survival prediction by the four lncRNAs are presented in Table 2.

Risk score = (expression quantity of AC093797. 1 × -0.4818) + (expression quantity of AL121748. 1 × 0.4404) + (expression quantity of AL162231. 4 × 1.2845) + (expression quantity of POLR2J 4 ×  − 1.5170).

The risk score formula in the training set was used to calculate the risk score of each patient. The median score for all patients in the training set was deemed the cutoff value, and facilitated the division of patients into a high-risk group (n = 51) and low-risk group (n = 52) (Table S4). Patients of the training group were ranked in ascending order of risk score (Fig. 3A). Moreover, the expression profiles of the four lncRNAs were plotted on a heatmap (Fig. 3B). In addition, a scatter plot was constructed to show the OS status and recurrence status of patients (Figs. 3C–3D). The high-risk group had a worse prognosis than the low-risk group in terms of OS and RFS. The prognosis of the high-risk group was compared with that of the low-risk group using Kaplan–Meier curve analysis. The results showed that the OS (HR: 3.650, 95% CI [1.761–7.566], log-rank P = 0.0002) and RFS (HR: 2.392, 95% CI [1.374–4.164], log-rank P = 0.0015) of high-risk patients were significantly lower than those of low-risk patients (Figs. 4A–4B). Furthermore, according to the four-lncRNA signature, the area under the curve of time-dependent ROC analysis of the 5-year OS and RFS predicted in the training set reached 0.839 and 0.715, respectively (Figs. 4C–4D).

Prognostic performance of the four-lncRNA signature in the testing set

To further assess the prognostic significance of the four-lncRNA signature, internal validation was performed in a group of 100 patients from the testing set. The established formula for the training set was used to calculate the patient risk score for the testing set. The same cutoff value used for the training set was also used to divide the testing set into high-risk (n = 44) and low-risk cohorts (n = 56) (Table S5). Kaplan–Meier survival analysis was also performed in the testing set. As shown in Figs. 5A–5B, the four-lncRNA signature showed favorable prognostic value in differentiating the risk stratification of death (HR: 2.475, 95% CI [1.347–4.547], log-rank P = 0.0015) and recurrence (HR: 2.245, 95% CI [1.349–3.736], log-rank P = 0.0014) in patients. The ROC curve was used to analyze the predictive power of risk models in the validation group for 5-year OS (Fig. 5C) and RFS (Fig. 5D). The OS and RFS of patients in the testing group were plotted (Figs. 5E–5F). These results suggested that the four-lncRNA signature in patients with cirrhotic HCC has considerable potential in predicting OS and RFS.

Stratification analysis

To explore the different clinical characteristics for which the four-lncRNA model is applicable, we conducted subgroup analyses of OS and RFS on patients of both the training set and testing set. The risk grouping principle was based on the median four-lncRNA risk score of the training set. Subgroup analyses combining the training and validation sets showed that the risk score model based on the four-lncRNA signature was more appropriate for patients with the following characteristics: BCLC stages 0–A, solitary tumors, and Tumor, Node, Metastasis system (TNM) stages I–II. The results of subgroup analyses of the OS and RFS are shown in Tables 3 and 4, respectively. Kaplan–Meier curves were constructed for both high-risk and low-risk patients with TNM stages I–II, solitary tumors, and BCLC stages 0–A, to analyze the prognostic value of risk models in OS (Figs. 6A–6F) and RFS of cirrhotic HCC (Figs. 7A–7F). Among the patients with the three aforementioned clinical characteristics, the OS (Fig. 6A (P = 0.00021), Fig. 6B (P = 0.0015), Fig. 6C (P = 0.0012), Fig. 6D (P = 0.0043), Fig. 6E (P = 0.00033), and Fig. 6F (P = 0.0019)) and RFS (Fig. 7A (P = 0.0018), Fig. 7B (P = 0.0019), Fig. 7C (P = 0.0028), Fig. 7D (P = 0.00061), Fig. 7E (P = 0.0023), and Fig. 7F (P = 0.0018)) of high-risk patients were worse than those of low-risk patients.

Functional enrichment analyses

Firstly, the co-expression matrix of the four lncRNAs and 12,100 protein-coding genes were extracted from 204 patients with tumor tissue gene expression profile data from GSE14520. The correlation between lncRNAs and protein-coding genes was analyzed. The significantly correlated protein-coding genes (Pearson coefficient >0.4, P < 0.001) were used for GO analysis (Table S6) and KEGG enrichment analysis (Table S7) to determine the potential mechanism of the four lncRNAs in regulating HCC. Functional enrichment analysis was conducted with the R org.Hs.eg.db (version 3.7.0) and clusterProfiler (version 3.10.1) package (Yu et al., 2012). The first 25 significant GO terms (P < 0.05, q < 0.05) (Fig. 8A) and first 25 significant KEGG pathways (P < 0.05, q < 0.05) (Fig. 8B) were charted using the R software (version 3.5.1). The results of GO analysis revealed that four-lncRNA-related functions showed enrichment in protein-coding genes involved in a large group of GO terms, including those associated with transmembrane transport, redox reactions, and fatty acid metabolism. Furthermore, KEGG analysis showed significant associations with a series of enriched pathways, including the peroxisome proliferator-activated receptor (PPAR) signaling pathway, amino acid metabolism, fatty acid metabolism, and chemical carcinogenesis.