Construction of a genome instability-derived lncRNA-based risk scoring system for the prognosis of hepatocellular carcinoma

Identification of GILncRNAs in HCC patients

The cumulative somatic mutations for each sample were first computed and ranked in the decreasing order. The top 25% samples and the last 25% ones in the ranking list were defined as genomic unstable (GU)-like group and genomic stable (GS)-like group, respectively (Figure 1A). Each group was consisted of 91 samples. Next, significant expressed lncRNAs were identified by comparing the lncRNA expression profiles between GU-like group and GS-like group. 88 lncRNAs were obtained based on the criteria of |fold change| > 1.5 and FDR adjusted P<0.05 and served as GILncRNAs. Of these, 32 lncRNAs were upregulated and 56 were downregulated in GU-like group (Supplementary Data 2). Then, we performed unsupervised hierarchical clustering analysis of all 364 samples using the filtered expression profile of differently expressed GILncRNAs. All samples were hence re-grouped based on the cluster result (Figure 1B). The group with higher cumulative somatic mutations was defined as GU-like group, and the other group was named as GS-like group. As shown in Figure 1B, the somatic mutation pattern was significantly different between the two groups. The median value of somatic cumulative mutations was 156.7 in the GU-like group while 124.0 in the GS-like group (P<0.001) (Figure 1C).

To validate the potential function of identified 88 lncRNAs and uncover their association with genomic instability, Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were carried out using clusterProfiler software in R-version 3.5.2. Firstly, protein coding genes (PCGs) closely related to the expression of lncRNAs were screened out using Pearson Correlation Coefficients. The top 10 PCGs most correlated with each lncRNAs were retained for lncRNAs–mRNA co-expression network construction. In the co-expression network, the nodes represented lncRNAs and mRNAs, and the lines represented the relationships between lncRNAs and mRNAs (Figure 1D). Go analysis of the PCGs in the network indicated that genomic instability might impact various metabolism process including fatty acid / fatty acid derivative metabolic process, pentose metabolic process, cellular aldehyde metabolic process, small molecule catabolic process and so on (Figure 1E). In terms of KEGG pathway analysis, PCGs in the network were found to enriched in a variety of metabolism pathway, HIF-1 signal pathway, extracellular matrix-receptor interaction, biosynthesis of amino acids and so on (Figure 1F). Enrichment results further reveal genomic instability-related lncRNAs were involved in various biological process of cancer. Expression of lncRNA could break the regulatory balance between lncRNA and PCGs, subsequently interfere with a variety of metabolism pathway, leading to gene damage repair process and exacerbation of genomic instability.

Construction of GILncSig for outcome prediction in the training set

To explore the role of GILncRNAs in the prognosis of HCC patients, 343 samples downloaded from TCGA were randomly divided into 2 groups named the training set (n=172) and testing set (n=171) respectively. As shown in Supplementary Table 1, there was no significant difference observed in the common clinical features characteristic between two groups (P > 0.05, Chis-square test). Then training set were used to establish GILncSig of HCC patients. We performed univariate Cox proportional hazard regression analysis to investigate the association of GILncRNAs and overall survival (OS) of HCC patients in the training set and found 9 GILncRNAs were closely related to the prognosis of HCC patients (P<0.05; Supplementary Table 2). Furthermore, multivariate Cox regression analysis was performed to evaluate the independent prognostic value of 9 GILncRNAs. 3 of 9 candidate lncRNAs including AC145343.1, AC004862.1 and ZFPM2-AS1 were obtained with prognostic significance in multivariate Cox analysis (P<0.05) (Table 1). Finally, GILncSig was established to predict outcome of HCC patients in the training set according to the equation mentioned above. The GILncSig was constructed as follow: GILncSig score = (0.3804 × expression level of AC145343.1) + (0.1253 × expression level of ZFPM2-AS1) + (−0.2344 × expression level of AC004862.1). A positive/negative regression coefficient demonstrated a positive/negative association between risk score and the expression level of lncRNA. Higher GILncSig score means higher risk of poor prognosis. Of the GILncSig, AC145343.1 and ZFPM2-AS1 tended to be risky factors while AC004862.1 was more likely to be a protective factor for the survival of HCC patients.

The GILncSig of each sample in the training set were calculated and then these patients were equally divided into high risk group and low risk group according to computed risk score with a decreasing order. Next we investigated the survival time of patients in the two groups using Kaplan–Meier analysis. The result showed that longer survival of patients in the low-risk group compared that in the high-risk.

Group (P<0.001, log rank test; Figure 2A). The time-dependent ROC curve was illustrated in Figure 2B and demonstrated an AUC of 0.781 for the GILncSig. We also plotted the expression levels of lncRNAs in the GILncSig and the count of somatic mutations with the increasing score in patients of the training set. As shown in Figure 2C, the risk lncRNA AC145343.1 and ZFPM2-AS1 showed up-regulated expression while the protective lncRNA AC004862.1 showed opposing expression pattern in the samples with high risk scores. Comparison analysis showed the number of somatic mutation between two groups has no significance (P = 0.16) but the average count of somatic mutation in patients of high-risk group was higher than that of low-risk group (Figure 2D).

Independent examination of GILncSig for HCC patients

The RNA-seq data of testing group including 171 HCC samples was analyzed to assess the stability, validity and prediction ability of the GILncSig. Patients in the testing group were also evenly separated into the high-risk group and low-risk group based on their GILncSig score. Kaplan–Meier survival analysis showed that patients in the high-risk group had a shorter survival time compared to those in the low-risk group (P<0.05) (Figure 3A; Left panel). The time-dependent ROC curves analysis of the GILncSig in the testing set yielded an AUC of 0.665 (Figure 3B; Left panel). Then we sorted the samples in the testing group according to GILncSig score with the increasing order to further demonstrate the alteration of the expression levels of the GILncSig (Figure 3C; Left panel). Concurring with the expectation, the risk lncRNAs AC145343.1 and ZFPM2-AS1 were more highly expressed and the protective lncRNA AC004862.1 showed lower expression pattern in the patients of testing set with higher risk scores. The number of the distribution of somatic mutation in patients with high scores was significantly higher than in patients with low scores (P<0.05, Figure 3D; Left panel).

Similar results have been obtained in the whole TCGA samples in terms of the prognostic value of the GILncSig. Kaplan–Meier plot analysis demonstrated that patients in the low-risk group had a more favorable survival time compared to those in the high-risk group (P<0.001) (Figure 3A; Right panel). As for ROC curve, the TCGA set yielded an AUC of 0.728 (Figure 3B; Right panel). The expression patterns of lncRNAs AC145343.1, ZFPM2-AS1 and AC004862.1 in TCGA set were consistent with that in the training and testing group (Figure 3C; Right panel). The high-risk group also showed increased somatic mutation counts compared to the low-risk group (P<0.05) (Figure 3D; Right panel).

Independence of the GILncSig from common clinical variables

To access the prognostic value of other clinical factors, we first conducted univariate Cox regression analyses on age, gender, tumor grade, tumor stage and GILncSig of each set. The result showed that only tumor stage and GILncSig were closely related to the survival time of HCC patients in each set (Table 2). Multivariate Cox regression was further analyzed to explore the independence performance of the GILncSig. Tumor stage and GILncSig also exhibited significant difference in each group based on the multivariate analysis result (Table 2). Then we continued to examine whether the prognosis performance of GILncSig was independent of the tumor stage. We excluded the patients with unknown stage and stratified remaining samples into early-stage group (stage I-II, n = 238) and late-stage group (stage III-IV, n = 83) for stratification analysis. Samples in each group were further separated into high- risk and low-risk group. A significant difference in survival outcome between high-risk and low-risk group was observed in both early-stage group and a late-stage group (P<0.01) (Figure 4). It was clear that the GILncSig could serve as an independent prognostic factor for the overall survival of HCC patients.

Comparison of the GILncSig with existing lncRNA-related signatures in prognosis value

To further illustrate the performance of GILncSig in our study, we recruited two recently published lncRNA signatures for survival prediction of HCC patient. Li’s study included 12-lncRNA signature (LilncSig) and Ma’s study included 4-lncRNA signature (MalncSig) [28, 29]. For LilncSig, Li et al. analyzed 12 pairs of HCC and adjacent normal mucosal tissues and identified 3900 differentially expressed lncRNAs as candidate biomarkers for the prognosis of HCC [28]. The 12-lncRNA signature was constructed using the least absolute shrinkage and selection operator (LASSO) cox regression method [8]. For MalncSig, the probe expression profiles of 225 HCC samples and 220 paired non-tumor tissue samples were derived from Gene Expression Omnibus (GEO)-GSE14520 [29]. Univariate coxregression and LASSO model were applied to screen lncRNAs linking to the overall survival. Then the multivariate Cox regression model was implemented to construct the prognostic score model [20]. Comparison analysis was performed between GILncSig and two recruited lncRNA signatures. As shown in Figure 5, the AUC of overall survival (OS) for the GILncSig is 0.728, which is significantly higher than that of LilncSig (AUC = 0.619) and MalncSig (AUC = 0.575). In addition, GILncSig consisted of 3 lncRNAs while MalncSig/LilncSig included 4/12 lncRNA for outcome prediction. Based on the AUC and lncRNA number, our optimized GILncSig significantly outperformed the two recently published lncRNA signatures in the sense of prognostic performance.

Comparison of the GILncSig with TP53 mutation status in prognosis value.

It is reported that mutation of TP53 gene is associated with worse survival in cancer patients [30]. Statistic results suggested that patients in the high-risk group have a significantly higher percentage of TP53 mutations than patients in the low-risk group among the training set, testing set and all TCGA samples (Figure 6A) (P<0.001). The proportions of TP53 mutation in the high-risk group of the training set, testing set and TCGA set were 49%, 42% and 45% respectively. In the low-risk group of the training set, testing set and TCGA set, TP53 mutation proportions were 20%, 12% and 16% respectively. TP53 transcriptionally targets hundreds of genes and regulates the expression of gene contributed to cell cycle, apoptosis, DNA repair proteins and metabolic [31, 32]. Recent studies had demonstrated that mutation of TP53 increased genomic instability and served as an independent prognostic marker [33–37]. Therefore, we continued to compare the performance of the GILncSig and TP53 mutation status in prognosis value. According to the GILncSig and TP53 mutation status, we classified all samples into TP53 Mutation/GS−like group, TP53 Mutation/GU−like group, TP53 Wild/GS−like group and TP53 Wild/GU−like group. Figure 6B revealed the survival curve of four risk groups. With regard to GU−like patients, the survival of TP53 Mutation group was more closely resembles that of TP53 Wild group. However, for patients with TP53 mutation, the survival of GU−like group was not similar to the GS−like group. Moreover, the survival outcome of patients in both GU−like group and TP53 mutation group were significantly worse than that in both GS−like group and TP53 wild group, indicating that the GILncSig and TP53 mutation status exhibited better prognostic performance than TP53 mutation status alone.

Unfavorable impact of AC145343.1 on HCC

Of the GILncSig, AC145343.1 served as a risky factor and the most important lncRNA for prognostic prediction according to the regression coefficient. Hence, we further assessed the function of AC145343.1 with regards to HCC. Firstly, we sought to characterize the potential prognostic ability of AC145343.1 indicated by survival analysis using The Encyclopedia of RNA Interactomes (ENCORI) as described in our previous study [38]. A significant decline of survival time in high AC145343.1 set were observed (Figure 7A). To evaluate the phenotype effect of AC145343.1 in vitro, we transfected siRNA into HepG2 cell to establish AC145343.1-downregulated cell line. To assess the proliferation inhibitory potential of AC145343.1 in HCC, we employed CCK-8 and colony formation assay in HepG2 with/without AC145343.1 downregulation. After AC145343.1 silencing, HepG2 exhibited obviously lower cell viability and significantly decreased colony area compared to the control group (Figure 7B, 7C). Invasion and migration roles of AC145343.1 were further observed, manifested by transwell assays and wound healing array. Transwell assay indicated silencing AC145343.1 remarkably decreased the number of HepG2 that migrated across the transwell chamber (Figure 7D). Wound healing array also revealed that AC145343.1-downregulated HepG2 exerted a significant delay in wound healing when compared with the control group (Figure 7E). Taken together, these results supported AC145343.1 knockdown repressed the proliferative, migratory and invasive abilities of HepG2 cells.

Data collection

Clinical characteristics, RNA sequencing (RNA-Seq) data and somatic mutation variation information of HCC patients were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). 377 samples with RNA expression profiles, survival information and common clinical characteristics were extracted for next analysis. Clinical and pathological characteristics of all samples were shown in Supplementary Data 1. Somatic mutation data of 364 HCC patients were also downloaded from TCGA. After integration of 377 samples with their somatic mutation data, 343 patients remained. These patients were randomly separated into two groups according to previous study, named training set and testing set respectively [63]. The training set with 172 HCC patients was used to recognize clinical outcome-related lncRNA signature and establish prognostic risk model. The testing set with 171 HCC patients was used to evaluate the performance of prognostic risk model of training group.

Genome instability-associated lncRNAs analysis

Somatic acquired genomic instability is one of the hallmarks of malignancy cancer [64]. Aberrant lncRNA levels are contributed to abnormal mutation and expression of genes involved in both tumor initiation and progression [65]. To recognize genome instability-associated lncRNAs (GILncRNAs), we constructed a computational frameworks integrating somatic mutator information and lncRNA expression profiles of tumor genome (Figure 8): (i) total cumulative somatic mutations of each sample was calculated; (ii) samples were ranked based on cumulative number of somatic mutations from high to low; (iii) The top 25% samples in the ranking list were set as genomic unstable (GU)-like group while the last 25% were set as genomically stable (GS)-like group; (iv) lncRNA expression profiles between the GU group and GS group were filtered; (v) GILncRNAs were identified according to the criteria: |fold change| > 1.5 and false discovery rate (FDR) adjusted P<0.05.

Statistical analysis

To identify GU group and GS group using profiles of differentially expressed lncRNAs, hierarchical cluster analyses were carried out based on Euclidean distances and Ward’s linkage method. Univariate and multivariate Cox proportional hazard regression analysis was performed to evaluate the prognostic value of the expression level of GILncRNAs. To obtain outcome prediction model using GILncRNAs, we combined the regression coefficients from the multivariate regression analysis with the expression data of identified prognostic lncRNAs and constructed the following equation of genome instability-derived lncRNA signature (GILncSig) according to previous study [63]:

$$\text{GILncSig}(\text{sample})={\displaystyle \sum _{\text{i}=1}^{\text{n}}\text{coef}({\text{lncRNA}}_{\text{i}})}\text{*}\text{expr}{\text{(lncRNA}}_{\text{i}}\text{)}$$

GILncSig (sample) in the formula is a risk score of the prognosis of HCC patient. lncRNAi indicates the ith prognostic lncRNA and expr (lncRNAi) represents the expression level of lncRNAi of HCC patient. coef (lncRNAi) is the contribution index of lncRNAi to prognostic risk score which was acquired from the regression coefficient of multivariate Cox analysis. Patients were further equally divided into high-risk group and low-risk group according the risk score. High-risk group of unfavorable prognosis with high GILncSig and low-risk group with low GILncSig were obtained to verified the performance of outcome prediction model. To demonstrate the survival rate and median survival of different group, we performed the Kaplan–Meier analysis and P<0.05 was considered to be meaningful. Independence validation among GILncSig and other important clinical factors were evaluated using multivariate Cox regression and stratified analysis. Hazard ratio (HR) and 95% confidence interval (CI) were obtained through Cox analysis. R-version 3.5.2 was used to compute ROC curves and ROC AUC to validate the performance of the constructed GILncSig.

Functional enrichment analysis

To unreal the co-expressed lncRNA-mRNA pairs, Pearson Correlation Coefficients were calculated based on the expression profile between every differentially expressed lncRNA and mRNA. The top 10 items were defined as significant co-expressed mRNAs of lncRNAs. To disclose the biological functional properties associated with the proposed GILncRNAs, functional enrichment analysis of significant co-expressed mRNAs including Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway was applied. The enrichment analysis was performed using clusterProfiler software in R-version 3.5.2 [63].

Cell proliferation and colony formation assays

Human liver cancer cell lines HepG2 were purchased from the Cell Bank of the Type Culture Collection of the Chinese Academy of Sciences, Shanghai Institute of Biochemistry and Cell Biology. Cells were cultured in Dulbecco’ s Modified Eagle Medium (DMEM) supplemented with 10% fetal bovine serum (FBS), 1% penicillin and 1% streptomycin (Gibco Life Technologies, Lofer, Austria), then incubated in a humidified incubator with settled parameters (37° C, 5% CO₂). AC145343.1 specific siRNA and negative control siRNA were obtained from Vigene Biosciences and transfected into HepG2 using X-tremegene siRNA transfection reagent (Roche Diagnostics, Shanghai, China) according to the standard guidelines.

Cell proliferation ability of HepG2 with/without AC145343.1 downregulation was evaluated by the Cell Counting KIT-8 (CCK-8, KeyGEN BioTECH, Nanjing, China) according to manufacturer's guidelines. Briefly, 3×10³ cells suspended were seeded in each well of 96-well plates and incubated overnight for cell attachment. The CCK-8 cell proliferation reagent (10 μl) was added to each well at 24 h, 48 h and 72 h. 4 h after CCK-8 administration, cell proliferation ability is detected. For colony formation assay, 1×10³ cells suspended were seeded in each well of 6-well plate and incubated for 2 weeks. Cell colonies were fixed with 4% formaldehyde solution and stained with crystal violet for image visualization.

Wound healing assay and transwell assay

Migration ability of HepG2 with/without AC145343.1 downregulation was observed by wound healing assay. 5×10⁵ cells suspended were seeded into a 6-well plate. When the cells permeated 90% of the plate, a ‘wound’ in cell monolayers was scratched using a 1 ml pipette tip. Then images of wounds were captured by an inverted microscope at 0 h, 24 h, 48 h and 72 h. Migration ability was calculated by calculating the wound confluence parameter. Invasive ability of cells was further accessed by transwell assay. The transwell chambers were placed in a 24-well plate which contained complete cultured medium (10% FBS), then added a layer of Matrigel before cell seeding. 8×10⁴ cells were seeded into upper transwell chambers and inserted into 300 μl serum-free medium. The cells of the upper surface were removed using a cotton swab 24 h later. Then the cells on the bottom surface were fixed with 4% formaldehyde solution and stained with 0.5% hematoxylin solution for 20 min. Images of invaded cells were collected by an inverted microscope.

Availability of data and materials

The datasets used during the current study are available from TCGA database (https://cancergenome.nih.gov/), and Supplementary Materials.