Identification of a novel DNA repair-related prognostic signature predicting survival of patients with hepatocellular carcinoma

Introduction

Hepatocellular carcinoma (HCC), which accounts for approximately 75% of liver cancers, is the sixth most common cancer and the third most common cause of cancer-related death worldwide.^1,2 Recently, epidemiological study has indicated that the worldwide incidence of HCC is highly heterogeneous due to variable risk factors, and that most HCC cases occur in South-eastern Asia and Saharan Africa,³ where hepatitis B virus (HBV) infection is the main risk factor.⁴ Diagnosis and treatment of HCC have greatly improved as medical technology rapidly develops, but prognosis remains very poor, with 5-year survival rates under 40%.⁵ Although multiple prognostic biomarkers, including alpha‐fetoprotein (AFP),⁶ glypican‐3,⁷ des‐γ‐carboxyprothrombin,⁸ and cytokeratin‐19,⁹ are widely used in clinical practice, their effectiveness can be influenced by various factors and varies from person to person; thus, the biomarkers are unreliable in clinical implementations.^10,11 Meanwhile, a large body of evidence has confirmed that molecular biomarkers play key roles in tumor diagnosis, prognosis assessments, and in identifying potentially high-risk patients. For example, CDX2 is under-expressed in patients with high-risk stage II colon cancer, suggesting its potential as an adjuvant to chemotherapy and a prognostic biomarker for stage II and III colon cancer.¹² Sulfite oxidase (SUOX) is a favorable diagnostic and prognostic biomarker of HCC and may forecast the performance and tumor recurrence risk when combined with serum AFP.¹³ However, the predictive power of a single indicator is insufficient. Conversely, gene signatures making use of several genes provide better predictive performance,^14,15 and multi-gene prognostic features can guide clinicians in selecting appropriate treatment. Identifying core genes and pathways in models can enhance clinical applications, identify new therapeutic targets, and provide new insights into cancer progression. Therefore, identifying effective and reliable biomarkers for prognostic prediction is urgently needed for HCC.^16,17

In this study, gene set enrichment analysis (GSEA) was used to explore bioinformatics information and to screen several genes for further analysis. Previous differential analyses have often compared gene expression in two groups, specifically seeking genes that were significantly up‐ or down-regulated, and have not examined biological information and related gene regulatory networks.^18,19 Meanwhile, GSEA does not require a differential gene threshold; it only examines whether genes in a set are randomly distributed or up- or down-regulated in some phenotypes relative to others, and identifies biological associations (related to functional information, processes, or disease states).

In the present study, we used mRNA expression data for enrichment analysis. GSEA showed that the hallmark DNA repair gene set was significantly up-regulated in the tumor phenotype, and we selected 103 enriched DNA repair genes to construct a seven‐gene signature that can predict the prognosis of HCC.

Materials and methods

Data collection and processing

Whole mRNA expression data and clinical information were obtained from the publicly available Cancer Genome Atlas database (https://cancergenome.nih.gov/). Patients meeting the following criteria were excluded: 1) patient had no follow‐up or survival status, 2) patient had clinical information but no mRNA expression data. Based on these requirements, 365 HCC samples and 50 adjacent noncancerous tissue samples were included in our study. Clinical information was extracted, including age, TNM classification, pathologic stage, histologic grade, history risk factor, and vascular invasion and neoplasm status, among others. The detailed clinical characteristics are summarized in Table 1.

Table 1 Clinical characteristics of patients with HCC in our study

Gene set enrichment analysis (GSEA)

Gene set enrichment analysis was used to determine whether the hallmark gene sets downloaded from the Molecular Signature database predicted significant differences between the HCC tumors and adjacent non-tumor tissue in the training set.²⁰ mRNA expression levels in the two groups were analyzed. Gene sets at nominal p<0.05 and false discovery rate (FDR) <0.25 were considered to be significantly enriched and to identify biological processes warranting further investigation.

Statistical analysis

mRNA expression profiles were extracted and log₂-transformed. Univariate Cox regression analysis evaluated mRNAs significantly associated with overall survival (OS) with p-values <0.05. Further screening was then performed using multivariate Cox proportional hazards regression analysis, and seven mRNAs were ultimately included. Hazard ratios (HRs) were calculated for each gene. Afterwards, a prognostic risk score model was constructed based on a linear combination of expression levels weighted by regression coefficients obtained from the multivariate Cox regression analysis according to the following equation:

where β is the HR for each gene. Three hundred sixty-five patients were classified into high‐ and low‐risk groups using a median risk score. A Kaplan‐Meier curve was used to compare the prognoses of the two groups, and differences were assessed using the log‐rank test. Student’s t-test was used to examine the differential expression of included genes in HCC tumor and adjacent non-tumor tissues. The chi-squared test analyzed mortality in the two groups. A ROC curve was drawn to analyze the sensitivity and specificity in predicting patient overall survival, and the optimal diagnostic threshold was evaluated using Youden’s Index. The area under the curve of the ROC curve (AUC) was examined with the Z test. P<0.05 was considered statistically significant. Univariate and multivariate Cox regression analyses were performed to investigate the relationship between signature risk score and clinicopathological features. All statistical analyses were conducted using R 3.3.1 (www.r-project.org) and SPSS 19.0 (IMB, Armonk, NY, USA).

Results

Gene set enrichment analysis and gene screening

Clinical information and matched mRNA expression data from 365 patients with HCC were obtained from TCGA. The hallmark gene sets in the Molecular Signatures Database comprise 50 gene sets, each of which summarizes and clearly defines a biological state or process and displays coherent expression. GSEA was used to search for significant differences between tumor and normal tissue in these gene sets. The results showed that eight gene sets were significantly enriched, with p-values <0.05 and FDRs <0.25, including the DNA repair, G2/M checkpoint, E2F targets, MYC targets V2, Mitotic spindle, MYC targets V1, Unfolded protein response, and Wnt/β-catenin signaling sets (Table 2, Figure 1). We selected the lowest p-value DNA repair process, which contained 103 core enriched genes for further investigation.

Table 2 Gene sets enriched in the HCC tumor phenotype (365 samples)

Figure 1 Eight gene sets significantly enriched in the HCC tumor phenotype using GSEA. Including DNA repair, E2F targets, G2/M checkpoint, Mitotic spindle, MYC targets V1, MYC targets V2, Unfolded protein response, Wnt/β-catenin signaling.

Identification of DNA repair-related genes related to survival

To identify genes associated with survival in HCC, univariate Cox regression analysis was applied. Twenty-six genes were significantly correlated with OS (P<0.00001), and multivariate Cox regression analysis was used to further pare down the panel. As a result, a total of seven genes, namely, ADA, FEN1, POLR2G, SAC3D1, UPF3B, SF3A3, and SEC61A1, were confirmed to be significantly associated with survival. These were used to construct a gene signature. The hazard ratio (HR) of each gene was calculated: ADA, SAC3D1, UPF3B, SF3A3, and SEC61A1 had HRs >0, and were considered high-risk factors associated with short survival; FEN1 and POLR2G had HRs <0, and were considered protective factors associated with long survival (Table 3). The mRNA expression levels of these genes in adjacent normal and tumor tissues were also investigated, and we found that the seven genes were all significantly up-regulated in tumor tissues (p<0.05, Figure 2A). In addition, we explored somatic mutations, including copy-number variations and point mutations in the seven genes of the prognostic signature though the cBioPortal website (www.cbioportal.org). The results showed that total mutation frequency was 6.53%, and all genes in the prognostic signature had missense mutations or copy-number amplifications in patients with HCC (Figure 2B and C), which suggests that these gene play an role in HCC.

Table 3 Seven mRNAs significantly associated with OS of patients with HCC

Figure 2 Seven screened mRNA which associated with OS. Notes: (A) Different expression of seven mRNAs in tumor and adjacent non-tumor tissues, ****represents for p-value <0.0001). (B) The total genetic alteration of 7 genes in the HCC sample. (C) Genetic alterations in each of the seven genes in the prognosis signature. Red and blue represent amplification and deletion, respectively.

Development and validation of a seven-gene signature for survival prediction

The seven screened genes were selected to formulate a gene signature. The risk score of every patient was calculated according to the formula:

Risk score = (0.167) x ADA value + (0.54) × SAC3D1 value + (0.254) × SEC61A1 value + (0.667) × SF3A3 value + (0.233) × UPF3B value + (−0.276) × FEN1 value + (−0.414) × POLR2G value. The patients were sorted according to risk score and divided into high- and low-risk groups. Risk score distribution and survival status of each patient is shown in Figure 3A. The high-risk group had a higher mortality rate. The mortality rate of the high-risk group was 44.3%, while that of the low-risk group was 26.0%. A chi-squared test showed significant differences in mortality between the two groups (p=0.0000). The Kaplan–Meier curve revealed that patients in the high-risk group had poorer prognoses (log-rank test p<0.0001; Figure 3B). The median survival for the high-risk group was 2.75 years, and that of the low-risk group was 6.94 years. A ROC curve was generated, and the signature had 68.9% prediction accuracy and 82.8% specificity with a Youden’s index of 0.517, showing good sensitivity and specificity with an average AUC of 0.746 (P<0.01) (Figure 3C). Furthermore, we compared the predictive value for survival of the mRNA-based signature and the other factors using ROC curve analysis. Our results show that T stage (AUC =0.638) and pathological stage (AUC=0.602) have higher predictive value for survival than do other routine factors. However, the 7-mRNA signature (AUC=0.746) may be better for predicting survival in HCC than these routine clinicopathological factors (7-mRNA signature versus T stage P-value =0.031, Figure 3D).

Figure 3 Evaluation and validation the survival prediction of the risk score signature. Notes: (A) Risk score distribution and survival status of each patient; (B) Kaplan-Meier survival curves shows a correlation between the expression of the seven-mRNA signature model and overall survival of patients. (log‐rank test P<0.0001); (C) ROC curve show the sensitivity and specificity of the seven-mRNA signature model in predicting the patient overall survival (AUC=0.746, P<0.01); (D) The ROC curves compare predictive value for survival of 7 mRNA signature with routine clinicopathologic.Abbreviations: ROC, Receiver operating characteristic; AUC, Area Under Curve.

Relationship between seven-gene signature risk score and clinicopathological features

We further investigated the relationship between risk score distribution and clinical parameter stratification and found that only tumor grade, pathological stage, and plasma AFP level are associated with our risk score (p<0.05, Figure 4). Other clinical features, such as age, gender, HBV or HCV virus infection, and status of liver fibrosis, had no effect on risk score. Thus, the score was related to three parameters, but independent of other clinical information. Furthermore, we used univariate and multivariate Cox regression to examine the combined ability of the signature risk score and other clinicopathological features to predict survival. Univariate Cox regression showed that tumor histological grade, tumor size, neoplasm status, new tumor events after initial treatment, distant metastasis, and signature risk score were significantly related to survival. The Kaplan–Meier curve suggested that the above results and conclusions were consistent, showing that patients with tumors, new tumor events during follow-up, large tumor size, high tumor grade, distant metastasis, or high risk scores have worse prognosis. Notably, the risk scores had greater prognostic value because they were significantly different between univariate and multivariate analyses, with p-values <0.05 (HR =2.38, 95% CI (confidence interval) =1.355–4.184), implying that the risk score is an independent prognostic indicator for HCC (Figure 5, Table 4).

Table 4 Univariate and multivariate Cox regression analyses for each clinical feature

Figure 4 The relationship between risk score distribution and clinical parameters stratification. Including age, gender, tumor pathologic stage, tumor grade, plasma AFP level and liver fibrosis status. Notes: *Represent for p-value<0.05, NS represent no statistical difference.

Figure 5 Kaplan–Meier suggested that some clinical features could effectively predict the survival of patients with HCC. Notes: (A) tumor size; (B) metastatic diagnosis; (C) pathologic stage; (D) new event after initial treatment; (E) neoplasm cancer status.

The patients were then stratified into different subgroups according to clinical parameters for further analysis of the relationship between signature risk score and patient prognosis. The Kaplan–Meier curve shows that only in the female subgroup, there was no significant difference between prognoses in the high- and low-risk groups, implying that HCC may be a disease that requires further explanation. There were significant differences between prognoses in the high- and low-risk groups in subgroups stratified by age, pathological stage, new tumor events after initial treatment, tumor status, alcohol consumption, or viral infection. The prognoses of patients in the high-risk groups were always poor. These results indicate that the risk score signature is a stable prognostic marker for patients with HCC and is not influenced by factors that often cause variations in the efficacy of traditional biomarkers (Figure 6).

Figure 6 Kaplan–Meier curve illustrates the prognostic value of risk score signature based on stratification of different clinical features. Notes: (A) Gender; (B) Age; (C) Pathological stage; (D) New event after initial treatment; (E) Person neoplasm cancer status; (F) Virus infection; (G) Alcohol consumption.

Discussion

HCC is a highly heterogeneous disease,^21,22 which increases the difficulty of treatment and prognostic prediction. There is evidence that commonly used clinicopathological parameters, such as TNM staging, age, gender, viral infection, and serum AFP levels, are not sufficient to accurately predict patient outcomes.^23–28 Therefore, a large number of mRNAs have been examined and identified as biomarkers in the process of cancer development and prognosis. For example, Zhang et al confirmed that ATPase family AAA domain containing 2 (ATAD2) is significantly overexpressed in gastric cancer and might serve as an independent marker of poor prognostic in patients.²⁹ However, such biomarkers are limited because they can be influenced by many factors, making their predictive ability unstable. Therefore, data mining techniques as well as univariate and multivariate Cox regression analysis are utilized to develop gene signatures containing various related genes and combining their predictive effects to enhance prediction efficiency. Such signatures are widely used in molecular diagnosis, individualized treatment, and accurate survival prediction,^30,31 and the prediction value is superior to those of single biomarkers.^18,32–34

Over the years, with the development of high-throughput sequencing technology, the combination of microarray data and bioinformatics has been widely used to find differentially expressed genes in various cancers and to find differential diagnostic and prognostic markers, which usually focus on significant differences. However, GSEA classifies genes based on a series of correlation calculations to determine whether a gene set can reveal the phenotype distribution of related genes. It only emphasizes a specific function rather than exploring a wide range of genes. In addition, the method can find biomarkers in just a few steps. In this study, we performed GSEA analysis using mRNA expression profiles from 365 HCC patients, and we found that eight biological processes or pathways were significantly enriched, with p<0.05 and FDR <0.25. The DNA repair process had the lowest p-value, and we chose it for further analysis. Next, univariate Cox regression analysis was used to identify mRNAs associated with survival in HCC. We found that 80 mRNAs were associated with survival, with p<0.05. We selected the 26 mRNAs most significantly associated with survival (p<0.0001) and further analyzed and identified 7 mRNAs that predicted survival in patients with a multivariate Cox proportional hazards model. These together were used to derive a signature risk score. According to the risk score, patients were separated into high- and low-risk groups. A Kaplan-Meier curve showed that the high-risk group had poor prognosis. Furthermore, we performed univariate and multivariate Cox regression analysis to explore the combined ability of the signature risk score and other clinicopathological parameters to predict survival. The results confirmed that the risk score is a stable, independent prognostic indicator and suggests that the risk score could be of important significance for patients with HCC as an effective clinical classification tool.

DNA damage stimulated by carcinogens or metabolic processes can be translated into genetic mutations, further leading to genomic instability and malignant transformation. To prevent this, DNA damage response and repair interactions counteract threats to genomic integrity; thus, changes in DNA repair capacity caused by genetic polymorphisms may be associated with cancer propensity.³⁵ Studies have shown that DNA repair genes are responsible for repairing different types of DNA damage,^36,37 and genetic polymorphisms have been found in almost all human DNA repair genes.³⁸ Meanwhile, polymorphisms in DNA repair genes are associated with various cancers, such as lung,³⁹ epithelial ovarian,⁴⁰ and pancreatic cancer.⁴¹ Several studies have predicted HCC survival using DNA repair-related genes. For example, Long, et al shows that the DNA repair genes XRCC4 and XXCC5 may be genetic markers for aflatoxin (AFB₁) exposure-induced hepatocellular carcinoma risk and outcome.⁴² Within the identified gene signature in our study, SAC3D1 is highly expressed in hepatocellular carcinoma tissues and is associated with poor prognosis. It can be used as a biomarker for prognosis of HCC.⁴³ Flap endonuclease 1 (FEN1) is a structurally specific nuclease that takes an essential role in a variety of DNA metabolic pathways and in the maintenance of genomic stability. Human FEN1 mutations that impair its exonuclease activity are associated with cancer development.⁴⁴ Meanwhile, Liu et al reported that functional FEN1 genetic variation increases the risk of colorectal cancer, HCC, gastric cancer, and esophageal cancer.⁴⁴ SF3A3 is a nucleic acid metabolism-related gene involved in DNA replication, repair, and transcriptional regulation.⁴⁵ SF3A3 is also an activity-suppressing genes of the tumor suppressor P53. Silencing of this gene increases the expression level of p53, thereby inducing tumor cell cycle arrest and cell death.⁴⁶ ADA primarily regulates immune responses by degrading adenosine.⁴⁷ The POLR2G, SEC61A1, and UPF3B genes have been rarely studied in relation to tumors. POLR2G is mainly involved in biological processes, and can mediate chromosomal imbalances that lead to abnormal regulation of genes in normal and ectopic tissues of the ovary.⁴⁸ SEC61A1 is the major subunit of the SEC61 complex and is the major polypeptide conduction channel in the endoplasmic reticulum membrane. Its missense mutation can cause genetic immune-related diseases, such as plasma cell deficiency.⁴⁹ UPF3B is located on the X chromosome, and its mutation can cause neurodevelopmental disorders.⁵⁰ In our study, all seven genes were significantly associated with the prognosis of HCC. However, no DNA repair gene signature has previously been developed for the prediction of HCC prognosis. Therefore, we used bioinformatics methods to identify a signature composed of ADA, FEN1, POLR2G, SAC3D1, UPF3B, SF3A3, and SEC61A1, which are associated with DNA repair, and demonstrated its prognostic value in HCC.

Conclusion

In summary, we developed a seven‐gene risk signature related to the DNA repair process that can forecast survival in HCC, with higher risk scores indicating poor prognosis. Additionally, its prognostic value remained stable despite variations in age, clinical pathological stage, historical risk factors and neoplasm status. This signature could serve as an effective and stable classification tool for clinical practice.