Liquid–Liquid Phase Separation-Related Genes Associated with Tumor Grade and Prognosis in Hepatocellular Carcinoma: A Bioinformatic Study

Abstract

Aim

The aim of the present study was to identify the association between tumor grade and liquid–liquid phase separation (LLPS)-related genes, and to generate a LLPS-related gene-based risk index (LLPSRI) as a prognostic tool for hepatocellular carcinoma (HCC).

Methods

Weighted gene correlation network analysis was performed to test whether the LLPS-related gene modules were associated with tumor grade of HCC. The candidate modules were subjected to functional enrichment analysis. We generated a LLPSRI using the expression profiles of the hub genes among the candidate modules in order to identify patients at high risk. Then, the biological characteristics of the high-risk patients were revealed using gene set enrichment analysis. Additionally, an independent external data set was used to validate the LLPSRI.

Results

Four gene modules showed a significant positive correlation with tumor grade and involved various cancer-related pathways. Among the hub genes, six were selected to generate the LLPSRI, which was significantly associated with prognosis of HCC patients. The LLPSRI could successfully divide patients with HCC into high- and low-risk groups, and patients in the high-risk group showed shorter overall survival than those in the low-risk group. E2F, MYC, and mTORC1 signaling may be important determinants of survival in the high-risk group. The prognostic value of the LLPSRI was validated with the independent external data set.

Conclusion

We identified LLPS-related gene modules that are associated with HCC tumor grade. The LLPSRI may be useful as a prognostic marker of HCC, and it may reliably stratify patients into groups at low or high risk of worse survival. Our analysis also suggests that certain biological characteristics of HCC may be associated with high risk of worse survival.

Keywords:

• hepatocellular carcinoma
• liquid–liquid phase separation
• WGCNA
• prognostic biomarker

Introduction

Hepatocellular carcinoma (HCC) is the most common type of malignant tumor, and it is a global health challenge whose incidence is growing worldwide.^1,2 Infection with hepatitis B virus (HBV) or hepatitis C virus (HCV) as well as cirrhosis are considered risk factors for HCC development.³ Our understanding of the pathophysiology of HCC has improved in recent decades, but this knowledge has yet to be translated into clinical practice. Novel concepts are required to reveal the complex mechanisms underlying this disease. In recent years, liquid–liquid phase separation (LLPS) has emerged as a new concept to explain how cells perform their vital functions.⁴ LLPS refers to the generation of membraneless intracellular condensates where essential processes occur, including chromatin organization, X-chromosome inactivation, transcription, DNA damage response, autophagy, and tumorigenesis.^5–7 Aberrant LLPS is associated with various diseases, including cancer.^8–11 For example, LLPS underlies the formation of nucleoli, whose morphology is considered an indicator of tumor malignancy and differentiation.¹² Numerous genes encoding proteins or RNA have been implicated in LLPS.⁹ The field of LLPS-mediated regulation of biological processes is still in an early stage, and much more needs to be learned about the molecular mechanisms involved and what they mean for cells and organisms.

At the same time, HCC is a highly heterogeneous disease,³ and certain patients with early-stage disease may be at higher risk of worse survival. A “one-size-fits-all” method to predict patients prognosis is unreliable, so novel prognostic models are urgently needed. Considering that LLPS is associated with HCC development, we hypothesized that the expression pattern of LLPS-related genes may be associated with prognosis of HCC. In the present study, we identified modules of LLPS-related genes that are associated with tumor grade, and from a subset of hub genes we created an index to predict prognosis of patients with HCC.

Materials and Methods

Data Processing

The TCGA-LIHC data set, including RNA sequencing (RNA-seq) expression profiles (displayed as read counts) and clinical information of HCC patients in The Cancer Genome Atlas (https://www.cancer.gov/), were downloaded to perform weighted gene correlation network analysis (WGCNA)¹³ and to generate an LLPS-related gene-based risk index (LLPSRI) in HCC. The RNA-seq profiles were normalized using the voom method in the limma package.¹⁴ To be included in our study, both gene expression profiles and clinical information had to be available, and the tumor had to be HCC. Data for healthy liver tissues, samples from patients with combined HCC and cholangiocarcinoma, and HCC samples of unknown tumor grade were excluded from our analysis. Another processed HCC data set (LIRI-JP) from the JP Project of the International Cancer Genome Consortium (https://dcc.icgc.org/) was downloaded from the HCCDB database¹⁵ and used to validate the LLPSRI. The LLPS-related genes were downloaded from PhaSepDB (http://db.phasep.pro/).¹⁶ The workflow of our study is shown in Figure 1. Our institution waived the requirement for ethical approval for our study because it retrospectively meta-analyzed data already deposited in public databases.

Figure 1 Workflow of the present study.

Abbreviations: TCGA-LIHC, the Cancer Genome Atlas-liver hepatocellular carcinoma; HCC, hepatocellular carcinoma; WGCNA, weighted gene correlation network analysis; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; LLPS, liquid–liquid phase separation; LASSO, least absolute shrinkage and selection operator.

Identification of Modules Associated with Tumor Grade Using WGCNA

A total of 2158 LLPS-related genes were found in the RNA-seq expression profiles of the TCGA-LIHC data set. We extracted the expression profiles of LLPS-related genes in order to perform WGCNA using the WGCNA package. First, hierarchical clustering analysis was performed, and the samples with significantly different gene expression patterns were removed. Then, the soft thresholding power value was screened during module construction. A suitable power value was identified if the degree of independence was > 0.9. The minimum module size was set to 30, and the individual module was assigned a unique color. If a module showed a significant positive correlation with tumor grade (P < 0.05), it was considered a candidate module.

Functional Enrichment Analysis of Candidate Modules

The clusterProfiler package¹⁷ was used to perform functional enrichment analysis to reveal the potential biological functions of candidate modules, based on Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Results were considered significant if they were associated with P < 0.01 after adjustment using the Benjamini & Hochberg method.

Hub Genes and the LLPSRI

In WGCNA, module membership (MM) was defined as the association between a gene and its module. In the present study, a gene with MM > 0.8 was considered a hub gene.¹⁸ The expression profiles of the hub genes were used to generate the LLPSRI. Briefly, univariate Cox regression analysis was carried out to identify overall survival (OS)-associated genes. A gene with P < 0.05 was considered significant. The gene expression profiles of the OS-associated genes were used to perform least absolute shrinkage and selection operator (LASSO) Cox regression analysis using the “glmnet” package.¹⁹ The features with non-zero coefficients were selected in 10-fold cross-validation for multivariate stepwise Cox regression, and the LLPSRI was defined as

LLPSRI=Expr_gene1*Coef_gene1+Exprgene2*Coef_gene2+Expr_gene3*Coef_gene3+...

where “Expr” represents the expression value of the candidate genes in the multivariate Cox regression analysis, and “Coef” is the corresponding estimated regression coefficient. Then, patients with HCC were divided into low- or high-risk groups based on the median LLPSRI. OS was compared between the low- and high-risk groups using the log-rank method.

Gene Set Enrichment Analysis (GSEA)

We performed GSEA^20,21 to explore the biological characteristics of high-risk patients. The hallmark gene set and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway gene sets from the Molecular Signatures Database (version 7.2)^20,22 were used as the reference gene set. GSEA was carried out using the GSEA software written in JAVA (http://www.broadinstitute.org/gsea). We set the cutoff criteria as a false discovery rate value < 0.05.

Validation of the LLPSRI

As we did with patients in the TCGA-LIHC dataset, we assigned each patient in the JP-ICGC an LLPSRI according to the above formula. Then, patients with HCC were divided into low- or high-risk groups based on the median PRS. OS was compared between the low- and high-risk groups using the log-rank method.

Results

Multiple LLPS-Related Gene Modules are Associated with HCC Tumor Grade

Four samples of HCC were removed as they showed significantly different expression patterns from the other samples (Figure 2A). In the end, 352 HCC samples were used for WGCNA. The power was 4, which was the lowest value for the scale with an independence degree of up to 0.90 (Figure 2B). WGCNA identified nine modules (Figure 2C), among which the black, turquoise, brown, and yellow modules showed a significant positive correlation with tumor grade and were therefore selected as candidate modules for subsequent analysis (Figure 2D). The green module, in contrast, negatively correlated with tumor grade. Gene significance (GS) was defined by WGCNA as the correlation of a gene with tumor grade. GS and MM correlated significantly with each other in these four candidate modules (Figure 2E-H).

Figure 2 Weighted correlation network analysis. (A) Four samples show significantly different expression patterns. (B) Network topologies at different soft-thresholding powers, indicated by the red numbers. (C) Dendrogram of genes clustered based on a dissimilarity measure (1-Topological Overlap Matrix). Each color represents a module in the gene co-expression network constructed by WGCNA. (D) Heatmap of the correlation between module eigengenes and tumor grade. (E–H) Gene significance and module membership correlated significantly with each other in the brown, black, turquoise, and yellow modules.

Tumor Grade-Related Gene Modules are Involved in Various Pathways

Analysis of GO indicated the genes of these four candidate modules were involved mainly in GOs related to nucleic acid metabolism (Figure 3A). KEGG pathway enrichment analysis showed that the genes were involved tumor-related pathways, such as the p53 signaling pathway, cell cycle, and FoxO signaling pathway (Figure 3B).

Figure 3 Functional enrichment analysis. (A) Gene Ontology analysis and (B) Kyoto Encyclopedia of Genes and Genomes pathway analysis of the four modules.

The LLPSRI May Be a Prognostic Indicator in HCC

Forty-three LLPS-related genes were considered as OS-associated genes in univariate Cox analysis (Table 1), and LASSO analysis showed nine of them to have non-zero coefficients (Figure 4A and Table 1). A subset of six genes (CDCA8, CENPA, KPNA2, NCL, PRC1, and RAD51) was selected by multivariate stepwise Cox regression to generate the LLPSRI (Table 1). The LLPRRI was significantly associated with OS [hazard ratio (HR) 2.718, 95% confidence interval (CI) 2.160–3.420, P < 0.001]. The patients were divided into high- and low-risk groups according to the median LLPSRI (Figure 4B), and patients in the high-risk group had shorter OS than those in the low-risk group (Figure 4C). In fact, the LLPSRI was associated with OS at an even greater significance level than some clinicopathological characteristics, including sex, age, American Joint Committee on Cancer stage, and vascular invasion (Figure 4D).

Table 1 The Liquid-Liquid Phase Separation-Related Overall Survival-Associated Genes and Coefficients

Features	Univariate Cox Analysis				LASSO	Multivariate Cox Analysis
	Beta	HR	95% CI for HR	P value	Coefficients	Beta	HR	95% CI for HR	P value
ANLN	0.304465	1.355899	1.189463–1.545623	0.000005	0.000000
AURKB	0.188270	1.207159	1.06828–1.364093	0.002535	0.000000
BOP1	0.207506	1.230605	1.048903–1.443782	0.010906	0.000000
CCNB1	0.305345	1.357093	1.17198–1.571445	0.000045	0.000000
CDCA8	0.392861	1.481213	1.270048–1.727487	0.000001	0.282526	0.635958	1.888831	1.314818–2.713443	0.000580
CDK1	0.260708	1.297848	1.135288–1.483685	0.000134	0.000000
CDKN3	0.183832	1.201814	1.05629–1.367388	0.005246	0.000000
CENPA	0.298749	1.348171	1.190937–1.526164	0.000002	0.074927	0.274736	1.316184	1.007311–1.719766	0.044078
CPSF6	0.961201	2.614836	1.717433–3.981155	0.000007	0.000000
DCAF13	0.406715	1.501875	1.199044–1.88119	0.000400	0.146485
DHX9	0.384644	1.469091	1.060259–2.035569	0.020799	0.000000
EIF3H	0.278869	1.321635	1.018506–1.714981	0.035913	0.000000
EME1	0.205810	1.228520	1.06892–1.411949	0.003748	0.000000
FANCI	0.226450	1.254139	1.069437–1.470741	0.005339	0.000000
HNRNPU	0.476012	1.609643	1.134128–2.28453	0.007710	0.000000
ILF3	0.411729	1.509426	1.09493–2.080833	0.011949	0.000000
KIF23	0.240368	1.271718	1.116499–1.448515	0.000296	0.000000
KIF2C	0.306683	1.358911	1.197292–1.542346	0.000002	0.000000
KPNA2	0.572342	1.772412	1.450444–2.165852	0.000000	0.294180	0.414368	1.513414	1.110543–2.062433	0.008693
LMNB1	0.260174	1.297156	1.107278–1.519595	0.001273	0.000000
MCM2	0.288300	1.334158	1.151078–1.546356	0.000129	0.000000
MCM3	0.294040	1.341837	1.100714–1.635782	0.003621	0.000000
MCM4	0.312781	1.367222	1.140203–1.639441	0.000735	0.000000
MCM6	0.404000	1.497804	1.240398–1.808626	0.000027	0.000000
MCM7	0.277971	1.320448	1.100858–1.58384	0.002741	0.000000
MKI67	0.246531	1.279579	1.128826–1.450465	0.000116	0.000000
NCL	0.830998	2.295608	1.611022–3.271102	0.000004	0.269699	0.487740	1.628631	1.083987–2.446928	0.018863
NEK2	0.211361	1.235358	1.089473–1.400777	0.000979	0.000000
PABPC1	0.330291	1.391374	1.150525–1.682642	0.000660	0.000000
PCLAF	0.200604	1.222140	1.05863–1.410905	0.006191	−0.031519
POLR2H	0.289243	1.335417	1.026984–1.736479	0.030875	0.117190
PRC1	0.219217	1.245102	1.078638–1.437255	0.002756	−0.021436	−0.646627	0.523810	0.37609–0.72955	0.000131
RAD51	0.208031	1.231252	1.067492–1.420134	0.004278	−0.148117	−0.363251	0.695412	0.531551–0.909785	0.008059
RAD54L	0.228929	1.257253	1.116881–1.415267	0.000151	0.000000
RFC4	0.307928	1.360603	1.117403–1.656736	0.002178	0.000000
RPL37	0.287910	1.333637	1.058401–1.680449	0.014637	0.000000
RPLP2	0.206228	1.229033	1.014312–1.489208	0.035290	0.000000
SNRPA	0.390740	1.478074	1.138796–1.918431	0.003316	0.000000
SNRPD1	0.311754	1.365818	1.078644–1.729449	0.009638	0.000000
SPC25	0.300570	1.350628	1.168214–1.561525	0.000049	0.000000
TOP2A	0.201420	1.223139	1.091744–1.370347	0.000513	0.000000
TTK	0.286449	1.331690	1.170482–1.515102	0.000014	0.000000
ZNF207	0.753354	2.124112	1.379752–3.270048	0.000621	0.000000

Abbreviations: HR, hazard ratio; LASSO, least absolute shrinkage and selection operator; CI, confidence interval.

Figure 4 The liquid–liquid phase separation-related gene-based risk index (LLPSRI) in the TCGA-LIHC data set. (A) Nine genes were identified as optimal features in LASSO analysis. (B) Risk factor correlation diagram. Top panel: risk predictions for each patient. The horizontal dotted line represents the median risk value, which distinguished a group at low risk (blue) and a group at high risk (red). Middle panel: patients ranked by predicted risk value, showing slightly longer survival time among patients at low risk than among those at high risk. Blue dots represent living patients; red dots, deceased patients. Bottom panel: expression of candidate genes. (C) The high-risk group showed worse overall survival than the low-risk group. (D) The LLPSRI is a prognostic factor independent of routine clinicopathological characteristics.

Biological Characteristics of Patients in the High-Risk Group

Five hallmark gene sets, hallmark E2F targets, hallmark G2M checkpoint, hallmark MTORC1 signaling, hallmark MYC targets v1, and hallmark unfold protein response, were significantly enriched in the samples of high-risk group (Figure 5A). Two pathways, cell cycle and spliceosome, were significantly enriched in the samples of high-risk group (Figure 5B).

Figure 5 Gene set enrichment analysis. (A) Five hallmark and (B) two Kyoto Encyclopedia of Genes and Genomes gene sets were significantly enriched in the high-risk samples.

Validation of the LLPSRI Using an Independent Data Set

All HCC patients in the LIRI-JP data set were assigned LLPSRIs using the above formula. The LLPRRI was significantly associated with OS (HR 2.939, 95% CI 1.886–4.578, P < 0.001). As we did with patients in the TCGA-LIHC data set, we divided the LIRI-JP patients into high- and low-risk groups according to the median LLPSRI (Figure 6A). Patients in the high-risk group had shorter OS than those in the low-risk group (Figure 6B).

Figure 6 Validation of the prognostic value of the liquid–liquid phase separation-related gene-based risk index using an external data set. (A) Patients with HCC were divided into high- and low-risk groups, whose risk score, survival and gene expression were compared. (B) The overall survival rate was higher among low-risk patients than among high-risk patients.

Discussion

Many gene-based signatures have been proposed for prognostic stratification of HCC.^23–25 Lin et al proposed a twenty gene-based prognostic stratification in HCC.¹⁸ Lv et al proposed a seven gene-based risk index for predicting early recurrence in HCC.²⁶ Huang et al reported a RNA binding protein-related prognostic signature for HCC.²⁷ These studies have usually not started from subsets of genes with particular biological significance. The present study, in contrast, started from the newly introduced concept of LLPS as important in health and disease. Focusing on genes related to LLPS, we defined several gene modules that were associated with HCC tumor grade. Our LLPSRI was generated using the hub genes in the modules and was found to be associated with OS. Furthermore, its prognostic value was validated with a second, independent data set. These results suggest that LLPS may be a useful concept for explaining the pathogenesis of HCC. The LLPSRI may provide a valuable reference for the personalized management of HCC patients.

According to our GSEA, the target genes of transcription factors E2F and MYC are significantly enriched in high-risk samples, so the networks regulated by these two transcription factors are likely related to LLPS. Previous studies suggested critical roles for E2F and MYC in hepatocarcinogenesis,^28–30 which the present analysis supports in the case of HCC. Our results justify further work to elucidate how E2F and MYC are involved in LLPS. In addition, we found that high-risk patients may show greater mTORC1 signaling, which is considered a major tumor-initiating pathway in HCC.³¹ This suggests that mTOR inhibitors may benefit our HCC patients in the high-risk group, but not those in the low-risk group. Thus, the LLPSRI may be useful for screening patients for potential response when conducting prospective trials of mTOR inhibitors.

Our LLPSRI consists of six LLPS-related genes, some of which have previously been associated with HCC. Aberrant expression of KPNA2 correlates with early recurrence and poor prognosis in patients with small HCC.³² PRC1 promotes early recurrence of HCC in association with the Wnt/β-catenin signalling pathway.³³ Silencing CDCA8 can suppress HCC growth and stemness.³⁴ RAD51 appears to influence the onset, progression and drug resistance of various cancers.³⁵ NCL and CENPA should be explored further to confirm and clarify their involvement in HCC.

The present study may provide new insights into the associations of LLPS with HCC tumor grade and prognosis, but several limitations should be noted. Most importantly, our analysis is based exclusively on bioinformatics, so our findings should be verified and extended in experimental studies, ultimately prospective studies with patients. Second, some of the genes that we identified here have yet to be linked to HCC. Indeed, whether the LLPSRI is linked to prognosis through causality or simple correlation needs to be explored.

Conclusion

We identified LLPS-related gene modules that are associated with HCC tumor grade. The LLPSRI may be useful as a prognostic marker of HCC, and it may reliably stratify patients into groups at low or high risk of worse survival. Our bioinformatic analysis also identified some biological characteristics of HCC that may be linked to worse survival. In these ways, our study generates numerous testable hypothesis to help guide future research into the mechanisms and personalized treatment of HCC.

Data Sharing Statement

The raw analyses from this study can be obtained from the corresponding author upon reasonable request.

Disclosure

Zhao-Shan Fang and Zhi Zhang are co-first authors for this study. 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.

Additional information

Funding

This work was supported by the Senior Talent Launch Fund of Wuming Hospital (2019-17), the Self-financing Research Project of the Guangxi Province Health and Family Planning Commission (Z20180668), the Project of the Nanning Scientific Research and Technology Development Plan (20153012), and the Guangxi Graduate Education Innovation Project (YCBZ2018041).