Integrated Proteomics and Bioinformatics to Identify Potential Prognostic Biomarkers in Hepatocellular Carcinoma

Introduction

Liver hepatocellular carcinoma (HCC) has been extensively studied as one of the most common and devastating malignancies.^1–3 Approximately 75–85% of primary liver cancers are caused by HCC and it is the third leading cause of death by cancer globally.⁴ The unclear molecular mechanism led to the poor progression of the HCC. Although numerous studies have a focus on the molecular understanding and treatment strategy of the HCC, patient survival remains at a dismally low rate which less than 8 months.^5,6 It is important to uncover novel genes with prognostic significance that shed light on improving the HCC patient’s outcome. Therefore, exploring new targets for early diagnosis and treatment of HCC is of great meaning.

A single gene or protein or any biomolecule is not responsible for any disorder or disease. However, a collection of genes coordinating with each other leads to disease occurrence. The weighted gene co-expression network analysis (WGCNA) was widely performed as an analysis method for large-scale data sets and find highly relevant gene modules, which are used to construct a free-scale gene co-expression network.^7,8 Related module detection helps in integrating various networks with varied molecular profiles to identify candidate biomarkers or therapeutic targets in many diseases,⁹ such as depressive order¹⁰ and chronic kidney disease.¹¹ Lau et al established an unfolding gene–gene relationship network based on differentially regulated genes of HCC.¹² In the current, network-based analysis has been widely employed to characterize specific physiological functions, signaling and metabolic networks, and genes with clinical significance in HCC.¹³

In this study, we aimed to profile a collection of novel genes as potential biomarkers for the diagnosis of HCC. We identified eight pivotal genes of HCC by the WGCNA, Protein-Protein Interaction (PPI), proteomic DEGs, TCGA cohort DEGs and TCGA cohort KM of HCC. The proposed algorithm was used in the development of module detection study specifically for co-expression network in HCC. An integrated methodological description highlights protein correlated to identify key genes for HCC, strengthening the experimental validation and drug therapeutics to diagnose and prognosis of HCC.

Methods

Patients and Tissue Samples

A total of 24 pairs of HCC and corresponding adjacent non-tumor specimens used for proteomics examination were simultaneously obtained from 24 HCC patients in Nanfang Hospital, Southern Medical University. All clinical samples in vitro were stored at -80 °C. Patients were not subjected to any neoadjuvant therapy before surgery. And the patients’ information was obtained from medical charts and follow-up. Informed consent was obtained from each patient. The study was approved by the Protection of Human Subjects Committee of Nanfang Hospital (approval no. NFEC-201801-K4). We obeyed the principles of the 1983 Declaration of Helsinki. All experiments in this paper obeyed this principle.

Data Independent Acquisition Proteomics

All the samples were lysed in 8 M urea-containing protease and phosphatase inhibitors (ThermoFisher Scientific, MA, USA). Protein concentration was determined by Bradford protein assay. Trypsin (Promega, Madison, WI) was added to digest the proteins at 37°C. DIA was performed with 45 isolation windows, and each window overlapped 1 m/z, total cycle time was 3.98 s. The DIA data were processed and analyzed by Spectronaut X (Biognosys, Schlieren, Switzerland) with default settings. The retention time prediction type was set to dynamic iRT. Spectronaut X determined the data extraction and the ideal extraction window dynamically depending on iRT calibration and gradient stability. FDR cutoff on precursor and protein level was lower than 1%. Decoy generation was set to mutated apply a random number of amino acid position swamps (min=2, max=length/2). Otherwise, all the selected fragment ions passing the filters were used for quantification.

Data Acquisition

Data were obtained from the following approaches: RNA sequencing data and clinical information (including gender, age, stage, clinical stage, pathological grade, tumor pathology subtype, overall survival time) for HCC data type and corresponding clinical features of 40 normal samples and 372 cancer samples were obtained from TCGA database (http://portal.gdc.cancer.gov).¹⁴

Differentially Expressed Genes (DEGs) Screening

In R 3.6.1 (R Foundation for Statistical Computing, Vienna, Austria), the optimal cutoff values for overall survival analysis of stromal score by maxstat package were used as the grouping criterion.¹⁵ Sequencing data were normalized and analyzed for differences. In order to reduce the false-positive rate, log2 |fold change (FC)| was corrected by the “ashr” method.¹⁶ |FC| > 2 and P < 0.05 were used as the screening criteria for DEGs.

Construction of Protein-Protein Interaction (PPI) Network

The genes related to clinical progression were imported into the STRING online database (https://string-db.org),¹⁷ and a PPI network was constructed with a mutual score greater than 0.4 as a threshold. The plug-in in Cytoscape v3.6.1 was used to screen the top 30 genes in the network as the central genes using five algorithms: betweenness, closeness, degree, Edge Percolated Component (EPC), and Maximal Clique Centrality (MCC).^18,19 The genes obtained by various algorithms were intersected to obtain the hub genes.

Functional Enrichment Analysis

Gene Ontology (GO) analysis is to describe the biological and molecular functions of genes and to describe different levels and dimensions of cellular components.^20,21 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis can bring together genes that are highly similar in sequence and perform the same function.²² DAVID (https://david.ncifcrf.gov) provides a comprehensive set of functional annotation tools for investigators to understand the biological meaning behind a large list of genes.²³ DEGs in the key module were put into DAVID for GO and KEGG analysis. P < 0.05 was set as the screening criteria for analysis results.

Weighted Gene Co-Expression Networks and Their Modules

A co-expression network targeting DEGs was established by using the R package Weighted gene co-expression networks and their modules (WGCNA) (version 1.68).²⁴ WGCNA is a freely accessible R software package (version R 3.6.0), developed for the construction of weighted gene co-expression networks. Rather than focus only on differential gene expression, WGCNA uses information from the genome to identify a set of genes of interest and converts the association of thousands of genes with phenotypes into associations between several gene sets and phenotypes, eliminating the problem of multiple hypothesis test correction. The parameter β is a soft-thresholding power parameter that strengthens strong correlations and penalizes weak correlations between genes. A hierarchical clustering tree was constructed with different branches of the tree representing different gene modules. The adjacency matrix was transformed into a topological overlap matrix (TOM). Genes were divided into different gene modules based on the TOM-based dissimilarity measure.

Gene Set Enrichment Analysis (GSEA)

The GSEA was used to screen the significantly changed pathways.²⁵ Pre-ranked GSEA was performed with 1000 permutations. The P-value was calculated by family-wise error rate (FWER), which is a robust method for multiples testing.²⁶ The GSEA plots were visualized by limma R package.²⁷

Statistical Analysis

Unpaired Student’s t-test was employed to perform the two-group comparison. t-test was utilized to analyze the difference between the T and P groups. A p-value <0.05 as well as |fold change| >2 was used to filter DEGs.

Results

Identification of DEGs Associated with the HCC and GO/KEGG Analysis

To identify the DEGs of HCC, we subjected sample from 24 pairs of samples including coupled HCC and adjacent non-tumorous tissues. As seen in Figure 1A, the heatmaps showed that 185 proteomic DEGs were significantly differentially expressed in carcinoma compared with adjacent non-tumor tissues (|FC|>2, p<0.05). What's more, 77 proteomic DEGs were significantly upregulated and 108 were downregulated. To analyze the biological functions of the highly significant module of 161 DEGs systematically, we constructed a PPI network based on STRING database and was visualized by Cytoscape 3.4.0 (Figure 1B). The PPI complex was filtered to obtain 152 nodes and 586 pairs of PPI relationships. Thereafter, to detect the fundamental function of DEGs, GO term and KEGG pathway enrichment analysis were performed. The results revealed that the red module was mainly enriched in the small molecule catabolic (GO: 0044282, P = 6.73e–19) and participated in the retinol metabolism pathway (hsa00830, P = 2.82e–15) (Figure 2).

Figure 1 DEGs identification and Construction of the PPI network. (A) The heatmap of DEGs between the LIHC and adjacent non-tumorous tissues. (B) Protein-protein interaction (PPI) network of 185 common DEGs. Upregulated DEGs were labeled in red and down-regulated were labeled in blue; p-value and fold change were log-transformed, and applied as the bolder with and node fill color; topological degree was used as the criteria for node size.

Figure 2 Functional analysis of common DEGs. (A) The top 10 significantly enriched GO-biological processes of the 185 DEGs; (B) The top 10 significantly enriched GO-cellular component of the 185 DEGs; (C) The top 10 significantly enriched GO-molecular function of the 185 DEGs; (D) The top 10 significantly enriched KEGG pathway enrichment analysis of the 185 DEGs.

Construction of a Weighted Gene Co-Expression Network

To decipher the biological network of HCC, co-expression modules of genes and the relationship between modules and pathogenesis of LIHC were performed by WGCNA analysis. As seen in Figure 3A, sample clustering was performed to detect variation and outliers of all data sets and no outlier was identified, indicating that all of the samples were used for the next step analysis. In the next step, we choose β = 7 as the soft threshold to construct a scale-free network (R² = 0.950) (Figure 3B). The densely interconnected gene modules were identified based on average linkage hierarchical clustering on all data sets using the matrix above (Figure 3C). Visualization of the eigengene network representing the relationships among the modules and the clinical trait weight are shown in Figure 3D and E. The results revealed that DEGs could be divided into six modules (blue, brown, green, turquoise, yellow, grey) based on their overall functions. Among these six modules, the turquoise, yellow, and green modules showed positive correlations and blue and brown modules showed negative correlations. One module labelled with turquoise was significantly co-expressed (p = 4e–12) (Figure 3F), suggesting its important function of genes for HCC process in the module.

Figure 3 WGCNA analysis of the LIHC proteome profiling. (A) Samples clustering were conducted to detect outliers between the LIHC and adjacent non-tumorous tissues; (B) Analysis of the scale-free topology model fit index for soft threshold powers (β); (C) Cluster dendrogram was generated by hierarchical clustering to show the modules of highly interconnected groups of genes between LIHC and adjacent non-tumorous tissues; (D) Interactive relationship analysis of co-expression genes. The light color indicates topological overlap, while the darker color indicates a high topological overlap. (E) The heatmap shows the eigengene adjacency. (F) Heatmap of the correlation between module eigengenes and clinical traits. P-value is shown in each color cell coded by the correlation between modules and traits (red indicates positive correlation; blue indicates negative correlation).

Identification of Pivotal Genes and Correlation Analysis with Clinical Characteristics

In the next step, we compared the proteomic DEGs, proteomic turquoise module, TCGA cohort DEGs and TCGA cohort KM of HCC. The results showed that eight common pivotal genes were identified (Figure 4A), including ACADS, HSD17B13, PON3, AMDHD1, CYP2C8, CYP4A11, SLC27A5, CYP2E1. We selected these eight genes with a hazard ratio (HR) <1 as the hub genes. We used these genes to build a prognostic prediction model. In addition, the module signature with other independent clinical variables integrated to establish a comprehensive model for monitoring progression in HCC (Table 1). We then choose seven independent risk features into our model consisting of age, tumor grade, pathological T stage, platelet count, weight, gender, stage and risk score for univariate (Figure 4B). The multivariate analysis showed that the HP of risk score, T stage and pathological T stage were greater than 1 (Figure 4C). The results suggesting that their higher expression was associated with a lower stage/risk score.

Table 1 Clinical Features of the Patients with LIHC

Figure 4 Construction and assessment of pivotal genes for predicting progression. (A) The common genes between our DEGs, proteomic turquoise module, TCGA cohort DEGs and TCGA cohort KM of LIHC. (B) Univariate-Cox regression analysis for screening appropriate and significant features into final nomogram model (C) Multivariate-Cox regression analysis for screening appropriate and significant features into final nomogram model.

Risk Score Analysis and KEGG Analysis for the Validation of Prognostic Model

Based on the median risk score, HCC samples were divided into low-risk (n =136) and high-risk groups (n=136). The distribution of risk score and survival status are shown in Figure 5A and B. The transcriptome data of 372 HCC patients were selected for the gene set enrichment analysis (GSEA) procedure using the model nomogram scores as the reference phenotype. We found that DNA replication, spliceosome, the cell cycle, the cellular senescence and microRNAs in cancer were upregulated in the high-risk group. However, the carbon metabolism, neuroactive ligand–receptor interaction, peroxisome, PPAR signal pathways and complement and coagulation cascades were downregulated in the low-risk group (Figure 5C). All of these aberrant pathways were enriched for hallmarks of malignant tumors with a false discovery rate (FDR) of <0.05.

Figure 5 Risk score analysis, Kaplan–Meier analysis and time-dependent KEGG analysis for the validation of prognostic model. (A) Rank of risk score and distribution of groups. Patients with EC were divided into low- and high-risk subgroups based on the median value of the risk score calculated. (B) The survival status and survival time of patients with EC ranked by risk score. In (A and B), green dots represent for patients with a low level of risk score and red dots represent for patients with a high level of risk score. (C) KEGG pathway showed five positively correlated groups and five negatively correlated groups.

Discussion

Cancer is increasingly a global problem, and primary liver carcinoma is the second top leading cancer with high mortality worldwide.^28,29 As one of the main histological types of primary liver carcinoma, HCC accounts for approximately 90% of pathogenesis³⁰ and leading to a much lower 5-year survival.³¹ Hence, it is significant to develop uncovering novel diagnostic and prognostic biomarkers of HCC, which would shed light on improving the HCC patient’s outcome. Up to now, some biomarkers have been verified to provide insights into disease outcomes. Lin et al elevated α-Fetoprotein (AFP) levels are useful prognostic factors for patients with hepatocellular carcinoma.³² Lee et al Identified that gamma-carboxy prothrombin (DCP) levels can predict prognosis in patients with HCC.³³ However, due to the high heterogeneity of HCC, more prognosis biomarkers need to be discovered.

In our study, a comprehensive analysis of DEGs was conducted in 24 pairs of HCC tissues and adjacent tissues using WGCNA. Firstly, we established protein-protein networks to find 77 upregulated proteomic DEGs and 108 downregulated proteomic DEGs by calculating the degree in Cytoscape. Then, GO term and KEGG pathway enrichment analysis were performed and found that red module was mainly enriched in the small molecule catabolic and participated in the retinol metabolism pathway, which has been reported that was closely related to liver disease.^34–36 Then the WGCNA analysis was performed to identify co-expression gene modules, verified the relationship between the gene network and phenotype, and found out key genes in the network. A module is defined as a group of genes that have separate function than other groups for a better understanding of molecular mechanisms of HCC. The expression levels of these key genes showed obvious differences among the significant modules. In the blue modules, these genes were significantly downregulated in HCC tissues, while upregulated in the green, turquoise and yellow module. Furthermore, the module labelled with turquoise was significantly co-expressed (P = 4e–12), suggesting its important function of genes for HCC process in the module.

Accordingly, we trained and validated eight pivotal genes related to prognosis, including ACADS, HSD17B13, PON3, AMDHD1, CYP2C8, CYP4A11, SLC27A5, CYP2E1. These genes were identified by the proteomic DEGs, WGCNA, PPI networks, TCGA cohort DEGs and TCGA cohort KM of HCC. It was reported that DNA methylation plays a key role in regulating ACADS which made it a potential therapeutic target for treating HCC.³⁷ HSD17B13 is a kind of genetic variants in retinol metabolism, which is associated with liver injury progression.³⁸ Lu et al found that HSD17B13 was a tumor suppressive gene that inhibits HCC cell glycolysis as revealed by glucose utilization, lactate production, and extracellular acidification ratio.³⁹ PON3 is a member of Paraoxonase (PON) proteins primarily expressed in the liver, its expression and specific activities were found to negatively correlate with several inflammatory disorders and served as a prognostic predictor in HCC.⁴⁰ Song et al identified AMDHD1 was as a novel tissue-specific gene in human and mouse liver tissues.⁴¹ CYP2C8, CYP4A11 and CYP2E1 are members of cytochrome P450 family which related iron ion binding and oxidoreductase activity and functioned in drug metabolism in the liver.^42–44 Elevating cytochrome P450 gene expression has been linked to the aggravation of various cancers and affects various regulated metabolites.⁴⁵ Zhu et al found that the expression levels of SLC27A5 responsible for the rate-limiting steps of fatty acid uptake.⁴⁶ Gao et al firstly found that SLC27A5 expression is downregulated in HCC by DNA hypermethylation, and reduced SCL27A5 expression contributes to tumor progression and poor prognosis.⁴⁷ Taken together, our findings indicate that the eight common genes in the turquoise modules are closely associated with the development of HCC and contribute to a common biological function of deciphering the biological network of HCC, including the growth, development and progression of HCC.

Furthermore, the prognostic signature with other independent clinical variables was integrated to establish a comprehensive model for monitoring progression in HCC. The univariate and multivariate analysis has shown that the HP of risk score, T stage and pathological T stage were greater than 1, indicating that these genes (ACADS, HSD17B13, PON3, AMDHD1, CYP2C8, CYP4A11, SLC27A5, CYP2E1) may be potentially new biomarkers to identify at-risk HCC patients when further validated. The higher expression was associated with a lower stage/risk score. A GSEA-based assessment of co-expression gene modules was conducted given the clinical significance of HCC, we sought to investigate the potential relationships among the co-expression gene modules with other clinical features. We found that DNA replication, spliceosome, the cell cycle, the cellular senescence and microRNAs in cancer were upregulated in the high-risk group, which are related to HCC tumor procession.⁴⁸ The above results revealed that these genes might be the novel potential prognostic biomarkers for HCC.

Conclusion

In conclusion, we found novel biomarker genes (ACADS, HSD17B13, PON3, AMDHD1, CYP2C8, CYP4A11, SLC27A5, CYP2E1) to predict the prognosis of HCC using integrated proteomics and bioinformatics. Functional analysis of these key genes showed enrichment in the HCC regulatory pathway. Furthermore, the higher expression of these genes was associated with a lower stage/risk score of HCC. Therefore, these genes could be potential prognostic biomarkers for HCC. Further biological and basic studies are needed to validate our findings.