Development and Validation of Genome Instability-Associated lncRNAs to Predict Prognosis and Immunotherapy of Patients With Hepatocellular Carcinoma

The pathophysiology of hepatocellular carcinoma (HCC) is prevalently related to genomic instability. However, research on the association of extensive genome instability lncRNA (GILnc) with the prognosis and immunotherapy of HCC remains scarce. We placed the top 25% of somatic mutations into the genetically unstable group and placed the bottom 25% of somatic mutations into the genetically stable group, and then to identify different expression of GILnc between the two groups. Then, LASSO was used to identify the most powerful prognostic GILnc, and a risk score for each patient was calculated according to the formula. Based on a computational frame, 245 different GILncs in HCC were identified. An eight GILnc model was successfully established to predict overall survival in HCC patients based on LASSO, then we divided HCC patients into high-risk and low-risk groups, and a significantly shorter overall survival in the high-risk group was observed compared to those in the low-risk group, and this was validated in GSE76427 and Tongji cohorts. GSEA revealed that the high-risk group was more likely to be enriched in cancer-specific pathways. Besides, the GILnc signature has greater prognostic significance than TP53 mutation status alone, and it is capable of identifying intermediate subtype groups existing with partial TP53 functionality in TP53 wild-type patients. Importantly, the high-risk group was associated with the therapeutic efficacy of PD-L1 blockade, suggesting that the development of potential drugs targeting these GILnc could aid the clinical benefits of immunotherapy. Finally, the GILnc signature model is better than the prediction performance of two recently published lncRNA signatures. In summary, we applied bioinformatics approaches to suggest that an eight GILnc model could serve as prognostic biomarkers to provide a novel direction to explore the pathogenesis of HCC.

Introduction

Hepatocellular carcinoma (HCC) is the most common cause of death in patients with chronic liver disease, and the fourth most common cause of death from cancer (Siegel et al., 2020; Eloranta et al., 2021). Its recurrence rate is up to 70%, and the 5-year survival rate is less than 60% with successful surgical resection or liver transplantation. Despite achieving great progress in the prevention, diagnosis, and treatment of this disease in recent years, the clinical outcome is still unsatisfactory (Forner et al., 2018; Villanueva, 2019). Hence, the identification of a novel and reliable prognostic molecular signature to screen appropriated therapeutic strategies or improve unfavorable prognosis is urgently needed for HCC cohorts.

Molecular biology studies have shown that there is genomic instability or genetic instability in HCC (Dore et al., 2001). Loss of heterozygosity (LOH) and microsatellite instability (MSI) caused by DNA mismatch repair (MMR) gene repair errors are considered to be the two main phenotypic features of genome instability (Tubbs and Nussenzweig, 2017). Genomic stability can be a hallmark of both human genetic disease and cancer (Andor et al., 2017). The relative stability of the genome is the basic prerequisite for cells to be faithfully passaged. When there are genetic defects or exposure to adverse environmental factors such as biological, physical, and chemical tests, it will lead to genomic instability (Gerlach and Herranz, 2020). Studies have constructed 10-miRNA signatures related to DNA damage response, and have shown that 10-miRNA signatures are associated with poor prognoses of ovarian cancer (Wang et al., 2017). More and more evidence indicates that lncRNA plays an important role in tumors, and abnormal lncRNA expression may affect tumor cell proliferation, tumor progression or metastasis (Zhang et al., 1993; Luo et al., 2006; Sanchez Calle et al., 2018). Therefore, constructing lncRNAs related to genomic instability may be a prognostic factor for HCC.

In this study, we download HCC patients from TCGA and GSE76427, and a systematic analysis was conducted on genome instability lncRNA signature (GILncSig) to screen different expression genes (DEGs) in HCC. After univariate and multivariate Cox regression analysis, eight GILncSigs were screened to establish the prognostic risk score model to evaluate the prognosis of HCC. Predicting genes were validated in testing set. Finally, GSE76427 was used to compare the expression level of eight GILncSigs between HCC tissues and the paired adjacent normal tissue (PANT).

Materials and Methods

Data Acquisition

RNA-seq datasets, as well as corresponding clinical information, were downloaded from the TCGA database (https://cancergenome.nih.gov/) and processed through the R package “TCGA-Assembler.” A total of 424 HCC patients (normal = 50, tumor = 374) were enrolled in our study from the TCGA dataset. GSE76427 (normal = 52, tumor = 115) was download from the GEO database (http://www.ncbi.nih.gov/geo) (Grinchuk et al., 2018). We divided all TCGA HCC samples into a training set and a test set. The training set included 184 samples for the creation of a clinical outcome lncRNA risk model. The test set included 181 patients and was used to validate the predictive ability of the prognostic risk model. Meanwhile, the somatic mutation information data was downloaded from TCGA (Supplementary Table S1). The detailed clinical characteristics of these patients are summarized in Table 1.

TABLE 1

TABLE 1. Clinical characteristics of HCC patient datasets in this study.

Selection of GILncSig in HCC

To identify genome instability-associated lncRNAs, a mutator hypothesis-derived computational frame combining lncRNA expression profiles and somatic mutation profiles in a tumor genome was developed: 1) the cumulative number of somatic mutations for each patient was computed; 2) patients were ranked in decreasing order of the cumulative number of somatic mutations; 3) the top 25% of patients were defined as genomic unstable (GU)-like group, and the last 25% were defined genomically stable (GS)-like group; 4) expression profiles of lncRNAs between the GU group and GS group were compared using significance analysis of microarrays (SAM) method; 5) differentially expressed lncRNAs (false discovery rate adjusted p < 0.05) were defined as genome instability-associated lncRNAs.

Limma packages in R language were used to screen the differentially expressed GILncSig in HCC. The default Benjamini-Hochberg FDR was employed to adjust the p-value to remove false positive results. Differently expressed GILncSig were processed by hierarchical clustering analysis with R package “pheatmap” (Robinson and Oshlack, 2010).

Functional Enrichment Analysis

Gene set enrichment analysis (GSEA) was processed by GSEA software version 3.0. used to enrich key signaling pathways. We computed the Pearson correlation coefficients to measure the correlation between the paired expression of lncRNAs and mRNAs, and the top 10 mRNAs were considered as co-expressed lncRNA-associated partners. To predict the potential functions of lncRNAs, we performed functional enrichment analysis of co-expressed lncRNA-associated mRNA partners to determine significantly enriched GSEA Gene Ontology (GO) terms and GSEA Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. The functional enrichment analysis was performed using clusterProfiler software in R-version 3.5.2. FDR <0.25 and |ES| >0.8 were the cutoff criteria to decide a statistical significance (Carlevaro-Fita et al., 2019).

Construction of a Risk Model With Prognostic Value in HCC

Hierarchical cluster analyses were performed using Euclidean distances and Ward’s linkage method. Univariate and multivariate Cox proportional hazard regression analysis was used to evaluate the association between the expression level of genome instability-associated lncRNA and overall survival. Based on the coefficients from the multivariate regression analysis and the expression levels of prognostic genome instability-associated lncRNAs, we constructed a genome instability-derived lncRNA signature (GILncSig) for outcome prediction as follows:

$\mathbf{G}\mathbf{I}\mathbf{L}\mathbf{n}\mathbf{c}\mathbf{S}\mathbf{i}\mathbf{g}\mathbf{R}\mathbf{i}\mathbf{s}\mathbf{k}\mathbf{s}\mathbf{c}\mathbf{o}\mathbf{r}\mathbf{e}={\displaystyle {\sum }_{\mathbf{i}=1}^{\mathbf{n}}\mathbf{c}\mathbf{o}\mathbf{e}\mathbf{f}}\left(\mathbf{I}\mathbf{n}\mathbf{c}\mathbf{R}\mathbf{N}\mathbf{A}\mathbf{i}\right)\ast \mathbf{exp}\mathbf{r}\left(\mathbf{I}\mathbf{n}\mathbf{c}\mathbf{R}\mathbf{N}\mathbf{A}\mathbf{i}\right)$

where GILncSig (patient) is a prognostic risk score for the breast cancer patient, lncRNAi represents the ith prognostic lncRNA, expr (lncRNAi) is the expression level of lncRNAi for the patient, and coef (lncRNAi) represents the contribution of lncRNAi to prognostic risk scores that were obtained from the regression coefficient of multivariate Cox analysis. The median score of the patients in the training set was used as a risk cutoff to classify patients into the high-risk group with high GILncSig or low-risk group with low GILncSig (Tibshirani, 1997; Simon et al., 2011).

Assessment of Immune Cells Infiltration

We performed the ssGSEA method to calculate the enrichment scores on the basis of metagenes. The reason why we considered the metagene a robust approach was the two main characteristics: 1) the use of a set of genes instead of single genes that represent one immune subpopulation because the use of single genes as markers for immune subpopulations can be misleading as many genes are expressed in different cell types; and 2) the assessment of relative expression changes of a set of genes in relation to the expression of all other genes in a sample. Referring to the Bindea et al. study, the author incorporate 535 metagenes that represented diverse 24 immune subpopulations: innate immune cells (dendritic cells [DCs], immature DCs [iDCs], activated DCs [aDCs], plasmacytoid DCs [pDCs], eosinophils, mast cells, macrophages, natural killer cells [NKs], NK CD56dim cells, NK CD56bright cells, and neutrophils) and adaptive immune cells (B cells, T cells, T helper cells, T gamma delta [Tgd] cells, T helper 1 [Th1] cells, Th2 cells, Th17 cells, regulatory T [Treg] cells, CD8⁺ T cells, T central memory [Tcm] cells, T effector memory [Tem] cells, T follicular helper [Tfh] cells, and cytotoxic cells).

Quantitative Reverse Transcription Polymerase Chain Reaction (qRT-PCR) Assays

A total of 150 HCC patients’ tissues and corresponding adjacent tissues were collected to explore the expression of eight GILncs in the tissue samples by using qRT-PCR, which was performed according to the manufacturer’s instructions.

Statistical Analyses

The normality of the variables was evaluated using the Shapiro-Wilk normality test. Kaplan-Meier analysis was used to generate survival curves using the “survival” and “survminer” packages, and the cut-off values were determined through the “surv_cutpoint” function in the packages. To calculate the hazard ratios and identify the independent prognostic factors, univariate and multivariate Cox regression analyses were performed using the “survival” package. Subsequently, we employed the area under the ROC curve (AUC) to measure the prediction accuracy. All statistical analyses were two-sided and considered p < 0.05 as the threshold for statistical significance. The statistical results were all analyzed by R (version3.6.2).

Result

Identification of GILncSigs in HCC Patients

The workflow of our study was shown in Figure 1. To identify the GILncSigs in HCC, the top 25% somatic mutations per patient (n = 93) and the least 25% somatic mutations per patient (n = 90) of the patients were assigned to GU-like (genomic unstable) group and GS-like (genomic stable) group. As shown in Supplementary Table S2, 245 GILncSigs were significantly different in expression between the GU-like group and the GS-like group (Figure 2A). Using 245 differentially expressed GILncSigs, all 374 patients were arranged into two clusters, the GS-like cluster and the GU-like cluster (Figure 2B). Next, we further analyzed the difference in the cumulative somatic mutation value between the GU group and the GS group, and the results showed that the cumulative somatic mutation value of the GU group was significantly higher than that of the GS group (Figure 2C). Researchers have discovered UBQLN4 is a newly identified driver of genomic instability, and harmful UBQLN4 mutations occur in families with autosomal recessive syndromes. Loss of UBQLN4 can lead to increased sensitivity to genotoxic stress and delayed DNA double-strand break (DSB) repair (Jachimowicz et al., 2019). Therefore, we compared the expression level of UBQLN4 in the GU group and the GS group, and the results showed that the expression level of UBQLN4 in the GU group was significantly higher than that in the GS group (Figure 2D).

FIGURE 1

FIGURE 1. Flow chart of the steps in the performed analyses.

FIGURE 2

FIGURE 2. Identification and functional annotations of genomic instability-related lncRNAs in patients with breast cancer. (A) Heatmap of the top 20 genome instability-associated lncRNAs expressing the most upregulation and downregulation. (B) Unsupervised clustering of 374 HCC patients based on the expression pattern of 245 candidate genomic instability-related lncRNAs. The left blue cluster is GS-like group, and the right red cluster is GU-like group. (C) Boxplots of somatic mutations in the GU-like group and GS-like group. Somatic cumulative mutations in the GU-like group are significantly higher than those in the GS-like group. (D) Boxplots of UBQLN4 expression level in the GU-like group and GS-like group. The expression level of UBQLN4 in the GU-like group is significantly higher than that in the GS-like group. Horizontal lines: median values. Statistical analysis was performed using the Mann–Whitney U test. (E) Coexpression network of genomic instability-related lncRNAs and mRNAs based on the Pearson correlation coefficient. The red circles represent lncRNAs, and the blue circles represent mRNAs. (F) Functional enrichment analysis of GESA GO for mRNAs co-expressed lncRNAs. (G) Functional enrichment analysis of GESA KEGG for mRNAs co-expressed lncRNAs.

To better understand the functions of these GILncSigs, we constructed a lncRNA-mRNA co-expression network, where the nodes are lncRNA and mRNA, if they are related to each other, lncRNA and mRNA will be linked together (Figure 2E). In addition, functional enrichment analysis of GESA for lncRNA-correlated mRNA, GSEA showed that lncRNA-correlated mRNA are mainly cancer-specific pathways (Figures 2F,G).

Development of a GILncSig for Prognosis of Patients With HCC in the Training Set

Next, 374 HCC patients were divided into a training set and a testing set according to best batches. In order to predict the clinical outcomes of HCC with GILncSig, we applied the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm to the 245 GILncSigs in the training set. Eight GILncSigs were selected to build the risk signature based on the minimum criteria, and the coefficients obtained from the LASSO algorithm were used to calculate the risk score for training set (Figure 3A). Next, to inspect whether the eight GILncSigs were related to prognosis in HCC, multivariate Cox regression was employed to analyze the hazard ratio (HR) of the eight GILncSigs in HCC. Forest plotting showed that high expression levels of six genes including AC026803.2, RHPN1-AS1, LINC00221, AL031058.1, ZFPM2-AS1, and THORLNC were significantly related to poor overall survival (OS) of HCC patients. Meanwhile, high expression levels of CR936218.2 and AL359915.1 were closely related to relative better OS in HCC patients (Figure 3B; Table 2). Compared with the GS group, the AC026803.2, RHPN1-AS1, LINC00221, AL031058.1, ZFPM2-AS1, and THORLNC level were significantly higher in the HCC group, however, the CR936218.2 and AL359915.1 levels were significantly lower in the HCC group (Figure 3C). The training set patients were assigned to low-risk and high-risk groups based on the median value of risk scores, and patients in the high-risk group had poor survival compared with the low-risk group (Figure 4A). To test the efficiency of GILncSig, ROC curve was constructed. The risk score’s AUC was 0.791, indicating that its efficiency to predict prognosis was accurate (Figure 4B). Simultaneously, as increasing GILncSig score, the change of expression in the training set (Figure 4E). We further analyzed the difference between the cumulative somatic mutation value between the two groups, and the cumulative somatic mutation value of the high-risk group was significantly higher than that of the low-risk group (Figure 4F). We also compared the expression level of UBQLN4 in the two groups, and the expression level of UBQLN4 in the high-risk group was significantly higher than that in the low-risk group (Figure 4G).

FIGURE 3

FIGURE 3. Identification of the genomic instability-derived lncRNA signature (GILncSig) for outcome prediction in the TCGA set. (A) The coefficients calculated by multivariate Cox regression using LASSO are shown. (B) Forest plots showing the results of the multivariate Cox regression between GILncSig expression and OS. (C) The expression level of 8 GILncs in GS and GU groups. (D) Kaplan–Meier estimates of overall survival of patients with low or high risk predicted by the GILncSig in the TCGA set. Statistical analysis was performed using the log-rank test and univariate Cox analysis. (E) Time-dependent ROC curves analysis of the GILncSig at 1 year (F) LncRNA expression patterns and the distribution of somatic mutation and UBQLN4 expression with increasing GILncSig score. The distribution of somatic cumulative mutations (G) and UBQLN4 expression in the (H) in the high- and low-risk groups for HCC patients. The red represents the high-risk group, and the blue represents the low-risk group. Horizontal lines: median values. Statistical analysis was performed using the Mann–Whitney U test.

TABLE 2

TABLE 2. Multivariate Cox regression analysis of genome instability-related lncRNAs associated with overall survival in HCC.

FIGURE 4

FIGURE 4. Performance evaluation of the GILncSig in the training set and testing set. Kaplan–Meier estimates of overall survival of patients with low or high risk predicted by the GILncSig in the training set (A) and testing set (C). Statistical analysis was performed using the log-rank test and univariate Cox analysis. Time-dependent ROC curves analysis of the GILncSig at 1 year in the training set (B) and testing set (D). LncRNA expression patterns and the distribution of somatic mutation count distribution and UBQLN4 expression for patients in high- and low-risk groups in the training set (E) and testing set (H). The distribution of somatic mutation in patients of high- and low-risk groups in the training set (F) and testing set (I). The distribution of UBQLN4 expression in patients of high- and low-risk groups in the training set (G) and testing set (J). Horizontal lines: median values. Statistical analysis was performed using the Mann–Whitney U test.

Validation of a GILncSig for Outcome Prediction in the Testing Set and TCGA Set

To examine the robustness of the GILncSig in the testing set and TCGA set. The results showed that patients in the high-risk group had poor survival compared with the low-risk group in the testing set and TCGA set (Figures 3D, 4C). The risk score’s AUC was 0.727 in the TCGA set and 0.671 in the testing set, indicating that its efficiency to predict prognosis was accurate (Figures 3E, 4D). Simultaneously, as increasing GILncSig score, the change of expression in the testing set and TCGA set (Figures 3F, 4H). We further analyzed the difference in the cumulative somatic mutation value between the two groups, and the results showed that the cumulative somatic mutation value of the high-risk group was significantly higher than that of the low-risk group in the testing set and TCGA set (Figures 3G, 4J). We also compared the expression level of UBQLN4 in the two groups in the testing set and TCGA set, and the expression level of UBQLN4 in the high-risk group was significantly higher than that in the low-risk group (Figures 3H, 4J). To explore the functional role of this prognostic signature in HCC progression, we performed GSEA analysis, and the results indicated that the high-risk group was more likely to be enriched in cancer-specific pathways (Figures 5A,B). The above-mentioned evidence demonstrated that the prognostic signature based on eight GILncSigs is reliable in HCC prognosis.

FIGURE 5

FIGURE 5. GSEA on the TCGA cohort to explore mechanisms underlying the 8- GILncSig. (A) GSEA GO identifies high- and low-risk related signaling pathways in HCC. (B) GSEA KEGG identifies high- and low-risk related signaling pathways in HCC.

Independent Prognostic Analysis of GILncSig

Next, multivariate Cox regression analysis was performed on the clinical factors and GIlncSig risk score to determine the independent prognostic value of GilncSig (Table 3). As shown in Figures 6A,B, in both the young and old groups, the survival rate of the high-risk group was significantly lower than that of the low-risk group (p < 0.05). Subsequently, the survival rate of the high-risk group was significantly lower than that of the low-risk group in the pathological stage (Figures 6E–J). However, there were different results in male and female patients (Figures 6C,D). Finally, we compared the resulting GIlncSig to the latest published signatures related to lncRNAs; the first signature is the 4-lncRNA signature (BailncSig), and the second signature is the 3-lncRNA signature (LilncSig). The results showed that the AUCs for BailncSig (Bai et al., 2021), LilncSig (Li et al., 2021), and GILncSig OS were 0.724, 0.671, and 0.727 (Figure 7C). Based on the above results, GlncSig has independent research value in HCC.

TABLE 3

TABLE 3. Univariate and multivariate Cox regression to identify independent prognosis predictor in both the TCGA cohort and the GSE76427 cohort.

FIGURE 6

FIGURE 6. Kaplan–Meier curves were performed for patients stratified by clinicopathological features in the TCGA set. Impact of prognostic risk on overall survival for patients younger than 65 years old (A) and older than 65 years old (B); for female (C) and male patients (D); for patients in G1-2 (E) and G3-4 (F); for patients in stage I-II (G) and stage III-IV (H); and for patients in stage T1-2 (I) and stage T3-4 (J).

FIGURE 7

FIGURE 7. Relationship between the GILncSig and TP53 somatic mutation. (A) The proportion of TP53 mutation in high- and low-risk groups in the training set, testing set, and the TCGA set. (B) Kaplan–Meier curve analysis of overall survival is shown for patients classified according to TP53 mutation status and the GILncSig. Statistical analysis was performed using the log-rank test. (C) The ROC analysis at 3 years of overall survival for the GILncSig, LilncSig, and BailncSig. (D) The abundance of each TME infiltrating cell in high-risk and low-risk groups. The upper and lower ends of the boxes represented the interquartile range of values. The lines in the boxes represented median value, and black dots showed outliers. The asterisks represented the statistical p value. (E) PD-L1 expression patterns for patients in high-risk and low-risk groups in the TCGA cohort. The correlation of GILncSig score with clinical response to anti-PD-1 immunotherapy. (F) CTLA4_negative + PD-1_negative, (G) CTLA4_negative + PD-1_positive, (H) CTLA4_positive + PD-1_positive, (I) CTLA4_positive + PD-1_negative.

GILncSig is Better Than TP53 Mutation Status in Predicting HCC Prognosis

It has been confirmed that the abnormality of the TP53 gene is closely related to the occurrence and development of HCC (Hill et al., 2018). As shown in Figure 7A, the proportion of patients with TP53 mutations in the high-risk group was significantly higher than that in the low-risk group (p < 0.05). Next, to evaluate whether GILncSig is better than TP53 mutation status predicting HCC prognosis, we divided all patients into TP53 mutation/GS, TP53 mutation/GU, TP53 wild/GS, and TP53 wild/GU groups. The results showed that patients in the TP53 mutation/GU group had poor survival as compared with the TP53 mutation/GS group in the TCGA set (Figure 7B).

The GILncSig in the Role of PD-1/L1 Immunotherapy

HCC is immunogenic and immunosuppressed, thus, the application of immunotherapy is of great significance for the development of THE treatment of HCC. It has been confirmed that the tumor microenvironment (TME) is correlated with the sensitivity to immunotherapy and HCC prognosis, and the ssGSEA method was used to calculate the enrichment scores on the basis of metagenes, and was employed in the gene expression data from three datasets to calculate an immune score. As shown in Figure 7D, the proportion of immune cells in the high-risk group was significantly lower than that in the low-risk group (p < 0.05), suggesting that the immune ability of the high-risk group is suppressed. Then, we analyzed the expression of PD-L1 between the high-risk and low-risk groups, and significantly higher PD-L1 expression in the high-risk group was observed relative to the low-risk group (Figure 7E). Next, we download the immunotherapy data of TCGA-LIHC patients, and the significant therapeutic advantages and clinical response to PD-1/L1 immunotherapy in patients in the high-risk group compared to those in the low-risk group were confirmed (Figures 7F–I).

Validation of the Expression and Prognosis of GILncSig in GSE76427 and Tongji Cohorts

To further validate GILncSig expression in HCC, GSE76427 was used to measure the expression level of GILncSig, and the results showed that compared with PANT group, the AC026803.2, RHPN1-AS1, LINC00221, AL031058.1, ZFPM2-AS1, and THORLNC levels were significantly higher in the HCC group, however, the CR936218.2 and AL359915.1 levels were significantly lower in the HCC group (Figures 8A,C). In addition, base on the 8-GILncSig, patients in the high-risk group had poor survival as compared with the low-risk group (Figures 8B,D).

FIGURE 8

FIGURE 8. Validation of the prognostic performance of GIlncSig on GSE76427 and Tongji datasets. The expression level of eight GILncs for the HCC and paired adjacent normal tissue (PANT) in GSE76427 (A) and Tongji datasets (C). Kaplan–Meier curves illustrated that patients with high risk had worse overall survival than those with low risk in GSE76427 (B) and Tongji datasets (D).

Discussion

The occurrence of a malignant tumor is a process of multi-gene participation and gradual evolution (Hanahan and Weinberg, 2011). The progression of a normal cell to a malignant cell is actually a long process that involves various genetic mutations that lead to a precancerous lesion and then to a malignant tumor. This series of genomic evolution often requires DNA damage or replication abnormalities to chromosomal instability and even the emergence of a “hyperploidy” phenotype. As a result, most tumors often present a complex genetic map at the time of diagnosis, which is very different from that of normal controls, suggesting a high degree of genomic instability in the body at the time of tumor development (Sieber et al., 2003; Pikor et al., 2013). Genomic instability is an important molecular feature of malignancy (Abbas et al., 2013). The relative stability of the genome is the basic prerequisite for faithful cell passage (Ovejero et al., 2020). Detection of genomic instability is now thought to be an early warning of tumorigenesis (Sieber et al., 2005). Current studies have also confirmed that genes that cause genomic instability are important clues to the causes of tumors (Petropoulos et al., 2019). Moreover, with the continuous improvement of modern molecular biology methods, more and more evidence shows that cancer patients can be timely treated by detecting their genomic instability (Ben-David et al., 2019). Clinically, there have been some targeted drugs targeting genomic instability-related genes, which have brought good news to tumor patients (O'Shaughnessy et al., 2011). It has been established that downregulation of these genomic instability-related genes has been clinically found to significantly enhance the sensitivity of cancer patients to platinum-based chemotherapy (Ledermann et al., 2012; Simmons et al., 2016). To sum up, it is of great significance to elucidate the relevant mechanisms of genomic instability in tumor cells and to conduct relevant assays to reduce the incidence of tumor, delay the progression, and improve the disease condition.

LncRNA is widely defined as a class of RNA molecules that have a transcriptional length is greater than 200 nucleotides and lack an open reading frame (Esteller, 2011). LncRNAs have a potential role in regulating the function of tumor cells (Bhan et al., 2017). LncRNA regulates gene expression at different levels including chromatin assembly, transcriptional, and posttranscriptional (Mercer et al., 2009; Schmitt and Chang, 2016). Recent studies have shown that NORAD and GUARDIN are essential for genomic stability (Hu et al., 2018; Munschauer et al., 2018). Currently, a variety of tumor genomic instability detection technologies have sprung up rapidly, and the understanding of the role of genomic instability in the development of tumors is gradually deepening. However, the genom-wide identification of genomic instability-related lncRNAs and the systematic exploration of their clinical significance in cancer are still in their infancy. Therefore, it is of great significance to identify lncRNAs associated with genomic instability.

First, we downloaded the TCGA-LiHC expression profile data and mutation data. We defined the 25% with the highest frequency as the high mutation group, and the 25% with the lowest frequency as the low mutation group. A total of 245 different lncRNAs were obtained compared with the low mutation group. These lncRNAs were characterized as genomic instability-related lncRNAs. Functional enrichment analysis of mRNAs co-expressed with 245 lncRNAs indicated that these lncRNAs may play an important role in the pathogenesis, DNA damage, and apoptosis pathways of HCC, which is consistent with other studies (Gillman et al., 2021). Abnormal repair of DNA damage is directly related to genomic stability (Tubbs and Nussenzweig, 2017). If the mechanism of repairing DNA damage is defective, it will directly lead to the persistence of DNA damage and the harmful changes of cells, until the tumor is triggered. DNA damage is mainly exogenous, such as chemical exposure, UV irradiation, biological hazards, and endogenous, such as in vivo spontaneous DNA damage events, cell cycle process and DNA replication process block. These damages, if not repaired in time, can induce genomic oxidation, alkylation, and even DNA crosslinking, dimer formation and even DNA breakage. Therefore, whether DNA damage can be repaired in time and correctly directly affects the maintenance of genome stability (Pommier et al., 2016; De Magis et al., 2019; Fang and Ding, 2020). Next, we constructed 8-GILncSig to predict the course of HCC patients. GILncSig divided patients into low-risk and high-risk groups. The mutation frequency of the high-risk group was higher than that of the low-risk group in both the training set and the test set. In addition, similar results were found in the GSE76427 dataset. Furthermore, we compared the resulting GIlncSig to the latest published signatures related to lncRNAs; the first signature is the 4-lncRNA signature (BailncSig), and the second signature is the 3-lncRNA signature (LilncSig). The results showed that GILncSig is superior to the other two models.

Among the eight GILncs, we found that three lncRNAs (RHPN1-AS1, LINC00221, and ZFPM2-AS1) have been reported to be related to HCC. RHPN1-AS1 was upregulated in HCC tissue, and silencing lncRNA RHPN1-AS1 also inhibited the activation of PI3K/AKT/mTOR pathway (Song et al., 2020). Human HCC samples had increased the expression of LINC00221. LINC00221 knockdown repressed HCC cell growth, migration, and invasion and enhanced their apoptosis (Yang et al., 2021). ZFPM2-AS1 was observed to be distinctly upregulated in HCC tissues and associated with shorter overall survival. Inhibition of ZFPM2-AS1 suppressed cell proliferation, metastasis, cell cycle progression while accelerated cell apoptosis (Zhang et al., 2021).

There are several limitations to this study. First, there is a lack of biological verification. Future molecular studies are needed to identify the interactions between the eight GIlncSigs from our prognostic model. Moreover, all data in this study is obtained from public databases, so future studies with prospective validation are still warranted.

Conclusion

In conclusion, this work identified and validated an 8-GIlncSig model that can predict the outcome of HCC. This prognostic model can be a promising tool for predicting the clinical prognosis and guiding strategies for HCC, which also provides a wider insight for understanding the HCC.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Ethics Statement

The studies involving human participants were reviewed and approved by Tongji Medical College, Huazhong University of Science and Technology. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author Contributions

YY and LR designed the study, performed statistical analysis, and drafted the manuscript. YL helped to draft the manuscript. LL conceived the study, participated in its design and coordination, and helped to draft the manuscript. All authors read and approved the final manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.763281/full#supplementary-material

References

Abbas, T., Keaton, M. A., and Dutta, A. (2013). Genomic Instability in Cancer. Cold Spring Harbor Perspect. Biol. 5 (3), a012914. doi:10.1101/cshperspect.a012914

PubMed Abstract | CrossRef Full Text | Google Scholar

Andor, N., Maley, C. C., and Ji, H. P. (2017). Genomic Instability in Cancer: Teetering on the Limit of Tolerance. Cancer Res. 77 (9), 2179–2185. doi:10.1158/0008-5472.can-16-1553

PubMed Abstract | CrossRef Full Text | Google Scholar

Bai, Y., Lin, H., Chen, J., Wu, Y., and Yu, S. A. (2021). Identification of Prognostic Glycolysis-Related lncRNA Signature in Tumor Immune Microenvironment of Hepatocellular Carcinoma. Front. Mol. Biosci. 8, 645084. doi:10.3389/fmolb.2021.645084

PubMed Abstract | CrossRef Full Text | Google Scholar

Ben-David, U., Beroukhim, R., and Golub, T. R. (2019). Genomic Evolution of Cancer Models: Perils and Opportunities. Nat. Rev. Cancer 19 (2), 97–109. doi:10.1038/s41568-018-0095-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhan, A., Soleimani, M., and Mandal, S. S. (2017). Long Noncoding RNA and Cancer: A New Paradigm. Cancer Res. 77 (15), 3965–3981. doi:10.1158/0008-5472.can-16-2634

PubMed Abstract | CrossRef Full Text | Google Scholar

Carlevaro-Fita, J., Liu, L., Zhou, Y., Zhang, S., Chouvardas, P., Johnson, R., et al. (2019). LnCompare: Gene Set Feature Analysis for Human Long Non-coding RNAs. Nucleic Acids Res. 47 (W1), W523–w529. doi:10.1093/nar/gkz410

PubMed Abstract | CrossRef Full Text | Google Scholar

De Magis, A., Manzo, S. G., Russo, M., Marinello, J., Morigi, R., Sordet, O., et al. (2019). DNA Damage and Genome Instability by G-Quadruplex Ligands Are Mediated by R Loops in Human Cancer Cells. Proc. Natl. Acad. Sci. USA 116 (3), 816–825. doi:10.1073/pnas.1810409116

PubMed Abstract | CrossRef Full Text | Google Scholar

Dore, M. P., Realdi, G., Mura, D., Onida, A., Massarelli, G., Dettori, G., et al. (2001). Genomic Instability in Chronic Viral Hepatitis and Hepatocellular Carcinoma. Hum. Pathol. 32 (7), 698–703. doi:10.1053/hupa.2001.25593

PubMed Abstract | CrossRef Full Text | Google Scholar

Eloranta, S., Smedby, K. E., Dickman, P. W., and Andersson, T. M. (2021). Cancer Survival Statistics for Patients and Healthcare Professionals - a Tutorial of Real‐world Data Analysis. J. Intern. Med. 289 (1), 12–28. doi:10.1111/joim.13139

CrossRef Full Text | Google Scholar

Esteller, M. (2011). Non-coding RNAs in Human Disease. Nat. Rev. Genet. 12 (12), 861–874. doi:10.1038/nrg3074

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, J., and Ding, Z. (2020). SNAI1 Is a Prognostic Biomarker and Correlated with Immune Infiltrates in Gastrointestinal Cancers. Aging 12 (17), 17167–17208. doi:10.18632/aging.103667

PubMed Abstract | CrossRef Full Text | Google Scholar

Forner, A., Reig, M., and Bruix, J. (2018). Hepatocellular Carcinoma. The Lancet 391 (10127), 1301–1314. doi:10.1016/s0140-6736(18)30010-2

CrossRef Full Text | Google Scholar

Gerlach, S. U., and Herranz, H. (2020). Genomic Instability and Cancer: Lessons from Drosophila. Open Biol. 10 (6), 200060. doi:10.1098/rsob.200060

PubMed Abstract | CrossRef Full Text | Google Scholar

Gillman, R., Lopes Floro, K., Wankell, M., and Hebbard, L. (2021). The Role of DNA Damage and Repair in Liver Cancer. Biochim. Biophys. Acta (Bba) - Rev. Cancer 1875 (1), 188493. doi:10.1016/j.bbcan.2020.188493

CrossRef Full Text | Google Scholar

Grinchuk, O. V., Yenamandra, S. P., Iyer, R., Singh, M., Lee, H. K., Lim, K. H., et al. (2018). Tumor‐adjacent Tissue Co‐expression Profile Analysis Reveals Pro‐oncogenic Ribosomal Gene Signature for Prognosis of Resectable Hepatocellular Carcinoma. Mol. Oncol. 12 (1), 89–113. doi:10.1002/1878-0261.12153

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanahan, D., and Weinberg, R. A. (2011). Hallmarks of Cancer: the Next Generation. Cell 144 (5), 646–674. doi:10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Hill, M. A., Alexander, W. B., Guo, B., Kato, Y., Patra, K., O'Dell, M. R., et al. (2018). Kras and Tp53 Mutations Cause Cholangiocyte- and Hepatocyte-Derived Cholangiocarcinoma. Cancer Res. 78 (16), 4445–4451. doi:10.1158/0008-5472.can-17-1123

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, W. L., Jin, L., Xu, A., Wang, Y. F., Thorne, R. F., Zhang, X. D., et al. (2018). GUARDIN Is a P53-Responsive Long Non-coding RNA that Is Essential for Genomic Stability. Nat. Cel Biol. 20 (4), 492–502. doi:10.1038/s41556-018-0066-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Jachimowicz, R. D., Beleggia, F., Isensee, J., Velpula, B. B., Goergens, J., Bustos, M. A., et al. (2019). UBQLN4 Represses Homologous Recombination and Is Overexpressed in Aggressive Tumors. Cell 176 (3), 505–519. doi:10.1016/j.cell.2018.11.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Ledermann, J., Harter, P., Gourley, C., Friedlander, M., Vergote, I., Rustin, G., et al. (2012). Olaparib Maintenance Therapy in Platinum-Sensitive Relapsed Ovarian Cancer. N. Engl. J. Med. 366 (15), 1382–1392. doi:10.1056/nejmoa1105535

CrossRef Full Text | Google Scholar

Li, L., Xie, R., and Lu, G. (2021). Identification of m6A Methyltransferase-Related lncRNA Signature for Predicting Immunotherapy and Prognosis in Patients with Hepatocellular Carcinoma. Biosci. Rep. 41 (6), BSR20210760. doi:10.1042/BSR20210760

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, J.-H., Ren, B., Keryanov, S., Tseng, G. C., Rao, U. N. M., Monga, S. P., et al. (2006). Transcriptomic and Genomic Analysis of Human Hepatocellular Carcinomas and Hepatoblastomas. Hepatology 44 (4), 1012–1024. doi:10.1002/hep.21328

PubMed Abstract | CrossRef Full Text | Google Scholar

Mercer, T. R., Dinger, M. E., and Mattick, J. S. (2009). Long Non-coding RNAs: Insights into Functions. Nat. Rev. Genet. 10 (3), 155–159. doi:10.1038/nrg2521

PubMed Abstract | CrossRef Full Text | Google Scholar

Munschauer, M., Nguyen, C. T., Sirokman, K., Hartigan, C. R., Hogstrom, L., Engreitz, J. M., et al. (2018). The NORAD lncRNA Assembles a Topoisomerase Complex Critical for Genome Stability. Nature 561 (7721), 132–136. doi:10.1038/s41586-018-0453-z

PubMed Abstract | CrossRef Full Text | Google Scholar

O'Shaughnessy, J., Osborne, C., Pippen, J. E., Yoffe, M., Patt, D., Rocha, C., et al. (2011). Iniparib Plus Chemotherapy in Metastatic Triple-Negative Breast Cancer. N. Engl. J. Med. 364 (3), 205–214. doi:10.1056/nejmoa1011418

CrossRef Full Text | Google Scholar

Ovejero, S., Bueno, A., and Sacristán, M. P. (2020). Working on Genomic Stability: From the S-phase to Mitosis. Genes (Basel) 11 (2), 225. doi:10.3390/genes11020225

PubMed Abstract | CrossRef Full Text | Google Scholar

Petropoulos, M., Champeris Tsaniras, S., Taraviras, S., and Lygerou, Z. (2019). Replication Licensing Aberrations, Replication Stress, and Genomic Instability. Trends Biochemical Sciences 44 (9), 752–764. doi:10.1016/j.tibs.2019.03.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Pikor, L., Thu, K., Vucic, E., and Lam, W. (2013). The Detection and Implication of Genome Instability in Cancer. Cancer Metastasis Rev. 32 (3-4), 341–352. doi:10.1007/s10555-013-9429-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Pommier, Y., Sun, Y., Huang, S.-y. N., and Nitiss, J. L. (2016). Roles of Eukaryotic Topoisomerases in Transcription, Replication and Genomic Stability. Nat. Rev. Mol. Cel Biol. 17 (11), 703–721. doi:10.1038/nrm.2016.111

CrossRef Full Text | Google Scholar

Robinson, M. D., and Oshlack, A. (2010). A Scaling Normalization Method for Differential Expression Analysis of RNA-Seq Data. Genome Biol. 11 (3), R25. doi:10.1186/gb-2010-11-3-r25

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanchez Calle, A., Kawamura, Y., Yamamoto, Y., Takeshita, F., and Ochiya, T. (2018). Emerging Roles of Long Non-coding RNA in Cancer. Cancer Sci. 109 (7), 2093–2100. doi:10.1111/cas.13642

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmitt, A. M., and Chang, H. Y. (2016). Long Noncoding RNAs in Cancer Pathways. Cancer cell 29 (4), 452–463. doi:10.1016/j.ccell.2016.03.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Sieber, O., Heinimann, K., and Tomlinson, I. (2005). Genomic Stability and Tumorigenesis. Semin. Cancer Biol. 15 (1), 61–66. doi:10.1016/j.semcancer.2004.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Sieber, O. M., Heinimann, K., and Tomlinson, I. P. M. (2003). Genomic Instability - the Engine of Tumorigenesis?. Nat. Rev. Cancer 3 (9), 701–708. doi:10.1038/nrc1170

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, R. L., Miller, K. D., and Jemal, A. (2020). Cancer Statistics, 2020. CA A. Cancer J. Clin. 70 (1), 7–30. doi:10.3322/caac.21590

CrossRef Full Text | Google Scholar

Simmons, A. R., Clarke, C. H., Badgwell, D. B., Lu, Z., Sokoll, L. J., Lu, K. H., et al. (2016). Validation of a Biomarker Panel and Longitudinal Biomarker Performance for Early Detection of Ovarian Cancer. Int. J. Gynecol. Cancer 26 (6), 1070–1077. doi:10.1097/igc.0000000000000737

PubMed Abstract | CrossRef Full Text | Google Scholar

Simon, N., Friedman, J., Hastie, T., and Tibshirani, R. (2011). Regularization Paths for Cox's Proportional Hazards Model via Coordinate Descent. J. Stat. Softw. 39 (5), 1–13. doi:10.18637/jss.v039.i05

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, X. Z., Ren, X. N., Xu, X. J., Ruan, X. X., Wang, Y. L., and Yao, T. T. (2020). LncRNA RHPN1-AS1 Promotes Cell Proliferation, Migration and Invasion through Targeting miR-7-5p and Activating PI3K/AKT/mTOR Pathway in Hepatocellular Carcinoma. Technol. Cancer Res. Treat. 19, 1533033820957023. doi:10.1177/1533033820957023

PubMed Abstract | CrossRef Full Text | Google Scholar

Tibshirani, R. (1997). The Lasso Method for Variable Selection in the Cox Model. Statist. Med. 16 (4), 385–395. doi:10.1002/(sici)1097-0258(19970228)16:4<385:aid-sim380>3.0.co;2-3

CrossRef Full Text | Google Scholar

Tubbs, A., and Nussenzweig, A. (2017). Endogenous DNA Damage as a Source of Genomic Instability in Cancer. Cell 168 (4), 644–656. doi:10.1016/j.cell.2017.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Villanueva, A. (2019). Hepatocellular Carcinoma. N. Engl. J. Med. 380 (15), 1450–1462. doi:10.1056/nejmra1713263

CrossRef Full Text | Google Scholar

Wang, T., Wang, G., Zhang, X., Wu, D., Yang, L., Wang, G., et al. (2017). The Expression of miRNAs Is Associated with Tumour Genome Instability and Predicts the Outcome of Ovarian Cancer Patients Treated with Platinum Agents. Sci. Rep. 7 (1), 14736. doi:10.1038/s41598-017-12259-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, L., Si, H., Ma, M., Fang, Y., Jiang, Y., Wang, J., et al. (2021). LINC00221 Silencing Prevents the Progression of Hepatocellular Carcinoma through Let-7a-5p-Targeted Inhibition of MMP11. Cancer Cel Int. 21 (1), 202. doi:10.1186/s12935-021-01819-w

CrossRef Full Text | Google Scholar

Zhang, X.-w., Li, Q.-h., Xu, Z.-d., and Dou, J.-j. (2021). STAT1-induced Regulation of lncRNA ZFPM2-AS1 Predicts Poor Prognosis and Contributes to Hepatocellular Carcinoma Progression via the miR-653/GOLM1 axis. Cell Death Dis. 12 (1), 31. doi:10.1038/s41419-020-03300-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Shields, T., Crenshaw, T., Hao, Y., Moulton, T., and Tycko, B. (1993). Imprinting of Human H19: Allele-specific CpG Methylation, Loss of the Active Allele in Wilms Tumor, and Potential for Somatic Allele Switching. Am. J. Hum. Genet. 53 (1), 113–124.

Google Scholar