Introduction
Liver cancer is one of the most commonly diagnosed cancers and the fourth leading cause of cancer death worldwide, with about 841,000 new cases and 782,000 deaths annually.1–3 Its prognosis remains dismal, although much progress has been made in both surgical and nonsurgical interventions over the past several decades. Due to the impact of liver cancer, a scoring system that can predict patient survival and inform accurate drug administration is urgently needed. Supported by recent progress in genetic profiling, including genomic microarrays and high-throughput sequencing technology in combination with bioinformatic analyses, we can now more readily identify targets for clinical treatment and prognostic prediction.
Liver cirrhosis caused by persistent viral infections (i.e. the hepatitis B and C viruses), which leads to cellular damage and altered tissue regeneration and inflammatory microenvironments, is currently the most significant risk factor for developing hepatocellular carcinoma (HCC).4–7 As early as 1986, the interesting phenomenon was noted that many of the same signal transduction pathways were involved in both tumorigenesis and wound healing.8 Subsequent studies confirmed that the similarities between tumorigenesis and wound healing include fibroblast recruitment and activation,9–13 extracellular matrix component deposition,14–16 infiltration of immune cells,17–20 neovascularization,21–23 and cellular lineage plasticity.24–26 Although there are some specific links between tumorigenesis and wound healing, the underlying mechanisms remain unclear. Fortunately, high-throughput bioinformatics may open a new avenue by which we may better understand these processes. In the present study, we used The Cancer Genome Atlas (TCGA), the International Cancer Genome Consortium (ICGC), and the Gene Expression Omnibus (GEO) database to generate a nomogram for evaluating the prognosis of HCC patients.
Methods
HCC data sources and preprocessing
Public gene expression data and clinical annotations were collected from the TCGA and ICGC databases. Patients without survival information were excluded. A TCGA-LIHI cohort (369 samples) and an ICGC-LIRI-JP cohort (225 samples) were enrolled in this study. FPKM-normalized, log2-transformed RNA expression data and clinical demographic data were downloaded from the respective websites (https://icgc.org/ , https://portal.gdc.cancer.gov/ ). These data were then merged to form a single metacohort after batch effect removal using the “sva” R package.
Unsupervised clustering of wound healing-related genes
A total of 565 wound healing-related genes were identified by mining Molecular Signature Database (MSigDB). Clustering was performed on the metacohort generated by merging the TCGA-LIHC and ICGC-LIRI-JP cohorts. The R package “Nbcluster” was used to determine the optimum number of clusters. The R package “Kmeans” was used to perform k-means clustering and to assign the clusters.
Gene set variation analysis (GSVA) and functional annotation
To investigate the differences in the biological processes between clusters with different healing patterns, we performed GSVA enrichment analysis using the “GSVA” R package. GSVA is a non-parametric and unsupervised method that is commonly used to estimate variation in pathways and biological processes in samples from expression datasets. The KEGG, GOBP, and HALLMARK gene sets were downloaded from the MSigDB database for use in the GSVA analysis. Adjusted p values less than 0.05 were considered statistically significant. The “cluster Profiler” R package was used to perform functional annotation for the m6A gene signature or other genes with an FDR cutoff less than 0.05.
Identification of differentially expressed genes (DEGs) between clusters with different healing patterns
To identify relevant wound healing-related genes, the empirical Bayesian approach included in the “limma” R package was used to identify DEGs between the different healing pattern clusters. The significance criteria for a DEG were set as an adjusted p value less than 0.05 and a log2-fold change greater than 1 or less than -1. Functional DEG annotation was performed with the “cluster Profiler” R package.
Generation of the heal.immune gene signature and heal.immune risk scores
We defined the heal.immune score for prognosis prediction based on the metacohort (merged TCGA-LIHC and ICGC-LIRI-JP cohorts) consisting of 594 samples with complete overall survival (OS) data. The meta cohort was randomly divided into a training set (70%) and a validation set (30%). Sixty wound healing-related genes with fold-changes greater than 1 or less than -1 between the two healing pattern clusters and twenty-five immune checkpoint-related genes (ICGs) were included for in silico signature development. A stepwise Cox estimation of high/low-stratified gene expression data was used to assess the prognostic value of the signature genes and to determine which genes to include in the heal.immune gene signature. Univariate Cox estimation excluded 25 genes with p values greater than 0.05, leaving 61 genes for subsequent multivariate Cox estimation. Fifty-five additional genes with p values greater than 0.02 were excluded via multivariate Cox estimation, leaving FCER1G, PLAT, ITGA5, CCNB1, CD86 and CD40 to be included in the heal.immune gene signature. Heal.immune risk scores were calculated as follows:
heal.immune risk score = 0.7701082 * f(FCER1G) − 0.6348783 * f(PLAT) + 1.1085626 * f(ITGA5) + 0.8197798 * f(CCNB1) -0.7765288 * f(CD86) + 0.5306283 * f(CD40)
f(gene) = 0 (gene ≤ cutpoint)
f(gene) = 1 (gene > cutpoint)
The cutpoint of each gene was determined via the “Surv_cutpoint” function from survival package in R.
For the calculation of the heal.immune risk score in the validation group, a similar algorithm was applied, with 2−ΔCT , instead of 0 or 1 as f (gene) value for each risk gene. The CT value of actin in each sample was used as a reference.
Validation the prognostic value of heal.immune gene signature in clinical samples
A total of 74 patients, as an outer validation group, were enrolled in this study that were from the Huashan Hospital of Fudan University between 2015 and 2017. All clinical samples were collected with informed consent from patients, and the research was approved by the Ethics Committee of Huashan Hospital, Fudan University (Shanghai, China).
RNA extraction and cDNA synthesis
Frozen tumor specimens were manually ground in liquid nitrogen using a mortar and pestle, instantly transferred into the lysis buffer, and homogenized using a needle and syringe. Total RNA was extracted using AllPrep DNA/RNA Mini kit (Qiagen, Hilden, Germany), according to manufacturer’s instructions. The quantity and quality of isolated RNA samples were determined by microliter spectrophotometer (ThermoFisher, Waltham, MA, USA). Afterwards, 1 µg of total RNA was converted into 200 µL of cDNA using the PrimeScript™ RT Reagent Kit (TaKaRa Bio, Shiga, Japan) according to the manufacturer’s instruction.
Real-time qPCR
The 10 µL reaction contained 5 µL TB Green Mix (TB Green Premix Ex TaqII; TaKaRa Bio), 0.2 µL ROX, 2 µL ddH2O, 2 µL cDNA template, and 0.4 µL of each primer (forward and reverse). Amplification and detection were performed using the ABI PRISM 7900 Sequence Detection System (Applied Biosystems Inc., Foster City, CA, USA). The real-time qPCR primer sequences used are provided in Table 1. Conditions for amplification were 95°C for 30 s, followed by 40 cycles of 95°C for 5 s, and 60°C for 30 s.
Table 1Primer sequences of the heal.immune gene set
| Gene | Forward primer | Reverse primer |
|---|
| FCER1G | AGCAGTGGTCTTGCTCTTACT | TGCCTTTCGCACTTGGATCTT |
| PLAT | AGCGAGCCAAGGTGTTTCAA | CTTCCCAGCAAATCCTTCGGG |
| ITGA5 | GGCTTCAACTTAGACGCGGAG | TGGCTGGTATTAGCCTTGGGT |
| CCNB1 | TTGGGGACATTGGTAACAAAGTC | ATAGGCTCAGGCGAAAGTTTTT |
| CD86 | CTGCTCATCTATACACGGTTACC | GGAAACGTCGTACAGTTCTGTG |
| CD40 | ACTGAAACGGAATGCCTTCCT | CCTCACTCGTACAGTGCCA |
| β-actin | GGACCTGACTGACTACCTCAT | CGTAGCACAGCTTCTCCTTAAT |
Immunohistochemical staining
For immunohistochemistry (IHC), 5-µm paraffin-embedded sections of patient/mice tumors were baked at 60°C for 1 h, deparaffinized in xylene, and rehydrated in a graded series of ethanol solutions. Antigens were unmasked by microwave heating of the samples in 10 mM sodium citrate buffer (pH 6.0) for 15 m (5 m each for 3 times), and the reaction was quenched using hydrogen peroxide (3%). After washing with phosphate-buffered saline, samples were incubated with the following primary antibodies, respectively, at 4°C overnight: anti-FCER1G (A12889; ABclonal, Woburn, MA, USA), anti-CD86 (A1199; ABclonal), anti-CyclinB1 (A19037; ABclonal), anti-TPA (10147-1-AP; Proteintech, Rosemont, IL, USA), anti-ITGA5 (27224-1-AP; Proteintech), or anti-CD40 (66965-1-Ig; Proteintech). DAB (3,3′-diaminobenzidine) was used as a detection system.
Nomogram construction
Of the 594 samples in the metacohort, 524 had complete clinical annotations. Table 2 summarizes the baseline clinical demographics. Samples were randomly divided into a training set (70%) and a validation set (30%). A nomogram was constructed to predict 1-, 3- and 5-year survival via a multivariate Cox regression model including high-low-stratified heal.immune risk score, TNM stage, sex, and age as variates in the “rms” R package. The accuracy of the nomogram’s predictive power was assessed via calibration and receiver operating characteristic (ROC) curves.
Table 2Clinical traits of the meta-cohort used for nomogram construction
| TCGA | ICGC | Metacohort |
|---|
| Age in years | | | |
| >60 | 195 | 148 | 326 |
| <60 | 165 | 39 | 198 |
| Unknown | 3 | 0 | |
| Sex | | | |
| Male | 245 | 140 | 277 |
| Female | 118 | 47 | 247 |
| Unknown | 0 | 0 | |
| Stage | | | |
| I | 170 | 29 | 198 |
| II | 84 | 81 | 164 |
| III | 81 | 65 | 146 |
| IV | 4 | 12 | 16 |
| Unknown | 24 | 0 | |
| Event | | | |
| Alive | 233 | 152 | 375 |
| Dead | 130 | 35 | 149 |
| Unknown | 0 | 0 | |
| Risk score | | | |
| High | 174 | 98 | 263 |
| Low | 189 | 89 | 261 |
Results
Identification of two healing patterns in HCC
A total of 565 wound healing-related genes reported in the Gene Ontology (GO) biological process database from the MSigDB were identified and used in this study. Genetic variation analysis showed that wound healing-related genes were rarely mutated, except for TP53 (Supplementary Fig. S1). We first compared the healing patterns between tumor and nontumor liver tissue to investigate healing dysregulation in HCC. We found a significant difference in the gene expression patterns in these tissue types (Fig. 1A). To better visualize this effect, we used principle component analysis to reduce the dimensionality of the data such that tumor and normal tissue could be categorized into two ellipses (Fig. 1C).
After batch effect removal, two large cohorts of RNA sequencing expression data, i.e. the TCGA_LIHC and ICGC_LIRI_JP cohorts, were combined into one metacohort containing 594 samples. Based on the expression levels of wound healing-related genes and k-means clustering, we identified two optimized patterns, which we named “healing cluster 1” (C1) and “healing cluster 2” (C2). C1 and C2 comprised 501 and 93 HCC samples, respectively (Fig. 1B). The expression levels of wound healing-related genes in C1 were much higher than those of C2. Kaplan-Meier curves showed that C1 patients had a relatively better OS compared with C2 patients (log-rank p test=0.016; Fig. 1D). These observations suggest that the expression patterns of wound healing-related genes can be used in HCC classification and prognostication.
C1 and C2 patterns differ in their associated biological processes
Next, we checked for differences in the biological processes of C1 and C2 by comparing the GSVA enrichment scores of gene sets from MSigDB. We found that the biological processes associated with C2 were enriched for immune cell accumulation and immune response (Fig. 2A). Functional annotation of the genes upregulated in C2 reinforced this finding (Fig. 2C). In C1, the enriched biological processes and functional annotation pointed to metabolism (Fig. 2B, D). Based on the above functional annotation, we compared the stromal scores, immune scores, tumor purity, and the inflammatory response to wounding enrichment scores of C1 and C2, and we found statistically significant differences in all cases (Fig. 2E–H). These results indicate that wound healing processes, especially a wound healing-associated inflammatory response, are differentially regulated in the two clusters.
HCC healing clusters are correlated with immune infiltration
To explore the differences between the C1 and C2 tumor microenvironments, we investigated immune cell infiltration and immune-related gene expression to characterize the immunological landscape. First, the numbers of eight types of infiltrated stromal cells were calculated using a microenvironment cell population counter (also known as an ‘MCP’ counter). Significantly higher numbers of stromal cells and immune cells, i.e. T cells, lymphocytes, natural killer cells, monocytes, dendritic cells, neutrophils and endothelial cells, were observed in C2 than in C1 (Fig. 3B). Notably, fibroblasts robustly accumulated in C2, consistent with the higher expression levels of wound healing-related genes. Furthermore, C2 exhibited higher expression levels of ICGs (Fig. 3D).
Among 565 healing-related genes, 60 were differentially expressed in the two clusters. Hierarchical clustering of the DEGs and ICGs partitioned them into two sets with distinct expression patterns (Fig. 3A). We analyzed the correlation between the two sets and found a significant positive correlation between the DEGs and ICGs (Fig. 3C). Collectively, these results suggest that the two HCC healing clusters have distinct immune microenvironments.
Development of the heal.immune score and validation of its prognostic value
As described above, we demonstrated that wound healing-related genes were associated with patient outcomes. Next, we aimed to devise a scoring system based on the wound healing- and immune checkpoint genes to predict HCC prognosis. The expression values of the genes were replaced with a stratification of high/low via a customized cutoff established using the “survival” R package. First, 565 wound healing-related genes (based on GO) and 25 ICGs were considered. Of these, 505 genes with fold-changes less than 2 between the clusters were discarded. A stepwise Cox estimation analysis was further used to narrow the gene list. Ultimately, six genes with p values less than 0.02 in the multivariate Cox estimation, including FCER1G, PLAT, ITGA5, CCNB1, CD86 and CD40, were used in the heal.immune gene set and established a scoring system (Fig. 4A). Based on univariate Cox estimation, risk scores, referred as heal.immune score, were calculated by adding the log (HR, Hazard ratio) of each gene multiply its value (1 for high, 0 for low), and the results were categorized as high risk or low risk using the median value as a cutoff (Fig. 4B, D). These genes have high expression heterogeneity in HCC and low correlation among themselves (Fig. 4C). We divided the 594 patients of the metacohort into a training group (70%) and a validation group (30%). Cox estimation verified the six-gene set as an OS risk factor in validation group (Fig. 4E). Moreover, the heal.immune score was highly related to healing cluster and patient OS. Kaplan-Meier curves showed that patients in the heal.immune score high group had much shorter OS compared with those in the low group of the validation group (Fig. 4G). Moreover, we validated the heal.immune gene set in the independent cohort from Huashan Hospital, Fudan University. The Fudan validation cohort verified the six-gene set as an OS risk factor (Fig. 4F). As shown in Figure 4H, the heal.immune gene set can also effectively discriminate the risk of OS. The IHC staining of heal.immune gene set of two patients with highest or lowest heal.immune score revealed a distinct expression between the two groups, which is consistent to the qPCR result (Fig. 4I).
Establishment of a nomogram using the heal.immune score in combination with clinical traits
The 524 (out of 594) patients with full clinical data were enrolled for nomogram construction. The baseline characteristics of the patients are shown in Table 2. The patients were divided into a training cohort (70%) and a validation cohort (30%). Multivariate and univariate Cox estimation showed that age, sex, TNM stage, and heal.immune score were independent risk factors for HCC (Table 2). A prognostic model constructed by fitting previous risk factors into a multivariate Cox model (Fig. 5C) confirmed that the heal.immune score is an independent risk factor for OS in HCC patients. Based on the strength of the multivariate Cox model, a nomogram integrating heal.immune score, age, sex, and TNM stage was constructed (Fig. 5A). A gross score was calculated by adding up all of the points. The calibration curve for predicting 3- and 5-year OS indicated that the nomogram-predicted survival closely corresponded with actual survival outcomes (Fig. 5B). For predicting 1-, 3- and 5-year survival, the area under the curve (AUC) values of the ROC curves were 0.82, 0.76 and 0.73 (Fig. 5D) in the training group and 0.84, 0.76 and 0.72 in the test group (Fig. 5E). This new nomogram based on age, sex, TNM stage, and heal.immune score provides a relatively accurate tool for predicting survival in HCC patients.
Discussion
Although HCC has been studied in great detail, methods for early diagnosis as well as treatment effects and prognosis have not been well characterized. For diagnosis and treatment, it is necessary to further understand its molecular mechanism of occurrence and progression. Thanks to advances in high-throughput sequencing technology, the relationships between genetic changes and immune infiltration and the occurrence and progression of diseases can be analyzed to obtain important bioinformatics-derived data that can be used for diagnosis, treatment, and prognosis.
It is widely recognized that there are biological similarities between wound healing and cancer features and progression; for example: lymphocyte infiltration in wound healing and tumor microenvironments; dermal cell migration and crawling in wound healing and epithelial-mesenchymal transition; wound healing- and cancer-associated fibroblast activation; and extracellular matrix (ECM) remodeling, which occurs during both processes. Unfortunately, these biological similarities have not attracted the attention of most researchers, especially in the HCC field. These two processes are distinctive to a degree. A key difference between wound healing and cancer is that wound healing is a self-limiting process, whereas tumors continue to expand, evolve, and spread. This difference is related to differences in the composition of the microenvironment. During wound healing, once re-epithelialization is complete, inflammation disappears, but this is not the case during tumorigenesis.
Here, we describe the analysis of a sample of 594 patients in TCGA and ICGC cohorts. By comparing the expression levels of 565 wound healing-related genes between HCC and normal tissue, we found that most genes were dysregulated in HCC, i.e. they differed in healing pattern, albeit with etiological similarity. We clustered patients into the C1 and C2 healing clusters based on the expression levels of wound healing-related genes. Interestingly, there were clear correlations between healing pattern and the degree of immune infiltration as well as the proportions of immune cells. Defying our expectation, the degree of immune cell infiltration was negatively correlated with OS.
Normally, the MHCI/II complexes of antigen-presenting cells display exogenous or endogenous antigens for detection by T cell surface receptors, thereby activating an immune response that kills target cells.27–29 However, T cells are also regulated by co-inhibitory signals that have negative effects on T cell activity, proliferation and survival, thus preventing excessive immune activation.30–32 Such opposing pathways represent a sophisticated self-protective mechanism. These immunosuppressive molecules, called immune checkpoints, are profoundly intertwined with tumor occurrence and progression. These factors help tumors evade the immune system and induce immune tolerance,33–35 thereby rendering large numbers of tumor-infiltrated immune cells ineffective at combatting tumors. Our current findings reflect this phenomenon. Twenty-five ICGs were found to be highly expressed in C2 patients. From an alternative perspective, C2 patients may benefit from immunotherapy, and immune checkpoint inhibitors may improve their prognosis.
Via stepwise Cox estimation, we established FCER1G, PLAT, CCNB1, CD40, CD86 and ITGA5 as a gene signature to calculate the heal.immune risk score. FCER1G encodes a high-affinity IgE receptor associated with leukocyte infiltration in dermatitis, and it plays a key role in allergic reactions.36,37 In neoplasms, the function of FCER1G remains elusive; however, in renal cell carcinoma and multiple myeloma, FCER1G upregulation leads to cancer generation and development via interactions with immune cells.38,39PLAT encodes a fibrinolytic enzyme that plays a role in cell migration and tissue remodeling.40–42CCNB1 encodes cyclin B1, which promotes cell proliferation and tumor growth in various human cancers.43–45ITGA5 encodes an integrin family member. Extensive research has revealed that integrins mediate multiple pathological processes, including thrombotic disease, infectious diseases, inflammation, fibrosis, and cancer, by mediating cellular adhesion to the ECM and engaging and activating downstream signaling pathways.46–50 Collectively, these risk genes are critically positioned to regulate immune cell infiltration, ECM remodeling, cell cycle regulation and cell-ECM interactions. Further exploration of their functions in HCC is urgently needed.
Finally, we established age, sex, TNM stage, and heal.immune score as independent prognostic risk factors based on multivariate and univariate Cox models, and these factors were incorporated into a nomogram. We confirmed the reliability of this scoring system by evaluating prognosis in the validation cohort. In this study, we integrated, for the first time, wound healing-related genes and ICGs to devise a metric to predict HCC prognosis, and we showed that this approach has good accuracy. The nomogram described here might serve as a new tool for clinical diagnosis and treatment.