Bioinformatics analysis of molecular genetic targets and key pathways for hepatocellular carcinoma

Introduction

Hepatocellular carcinoma (HCC), the second most common cancer in the world, is a major contributor to cancer mortality and incidence.¹ More than 580,000 new cases are increased every year in Asia, with China accounting for more than 50%.² Several methods have been used for decreasing the incidence rate, containing early detection and diagnosis. These methods are well efficacious and could reduce the mortality and incidence. However, the 10-year survival rate remains unsatisfactory in patients.³ As is well known, HCC is a heterogeneous disease. The occurrence and development of HCC are associated with the cellular context, environmental influences, and gene aberrations.¹ Recently, several studies suggested that cellular pathways and multiple genes may participate in the progression of HCC.^4–6 However, the underlying molecular mechanisms remain unknown. Thus, it is urgent and important to dig the hub molecules and to uncover the key molecular mechanisms.

The miRNAs are a series of small, non-coding regulatory RNAs, which regulate gene expression and perform a crucial role in the regulation of proliferation, cellular differentiation and cell death. A part of miRNAs was regarded as tumor suppressor genes, or oncogenes that contribute to several cancer progressions. Notably, some key miRNAs provide a novel approach for HCC therapy. For instance, microRNA-21, microRNA-26a and microRNA-29a-3p have been reported association with HCC and could serve as biomarkers;^7,8 Shigoka et al also showed that plasma miR-92a could effectively discriminate HCC from the control subjects.⁹

In recent years, a large amount of data about transcriptome of cells have been established using numerous of high‑throughput technologies. In particular, the microarray is a novel high-throughput platform for the analysis of gene expression. During tumorigenesis, the microarray technique, as a useful tool for the analysis of general genetic alteration, has been extensively used. Through accurate microarray analyses, the key factors processed by biochip data extraction, biological data clustering and sequence alignment can be identified. These analytical techniques provide a new and quick method to explore the molecular pathogenesis mechanism of various types of cancers.

In this study, 3 gene expression profiles and 2 miRNA expression profiles were downloaded from the Gene Expression Omnibus (GEO) database to obtain differentially expressed genes (DEGs) and differentially expressed (DE) miRNAs between the HCC tissues and normal liver tissues samples. Furthermore, functional enrichment and the PPI network analyses were used to identify the DEGs. Meanwhile, with the mRNA–microRNA interaction analysis, the ten hub genes in DEGs were tested and four key genes and miRNAs as well as their potential molecular mechanisms were studied.

Materials and methods

Microarray data

Three gene expression profiles (GSE84006, GSE14323, GSE14811) and two miRNA expression profiles (GSE40744, GSE36915) were retrieved from the GEO database.¹⁰

Identification of DEGs and DE miRNAs

Limma is a Bioconductor package of R that offers a comprehensive solution for the analysis of gene expression data. limma has been widely used for gene detection.¹¹ This package contains sophisticated facilities for reading, normalizing and exploring for large datasets. Moreover, it also applies an adjusted P-value (adj. P) to help correct false-positives. Herein, the adj. P for selection of DEGs and DE miRNAs was set as <0.05. Funrich webserver (http://funrich.org/index.html) was utilized to analyze the overlap of DEGs in the three datasets and the overlapping DEMs in the two miRNA expression profiles.¹²

Functional and pathway enrichment analysis

Gene ontology (GO) is widely used in annotating genes, gene products and sequences. The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a comprehensive database for biological interpretation of genome sequences and other high-throughput data. Both analyses were available in the Database for annotation, visualization and integrated discovery database (DAVID database), which is a bioinformatics data resource composed of an integrated biology knowledge base and analysis tools to extract meaningful biological information from large quantities of genes and protein collections.¹³ Herein, GO and KEGG analyses were applied by using the DAVID database to identification of DEGs. The cut-off criterion was set (P-value <0.05).

Construction of PPI network and analysis of modules

The online database of STRING (https://string-db.org/) was applied to assess the PPI containing direct (physical) and indirect (functional) associations.¹⁴ Then, the Cytoscape program, an open-source tool, was used for the visual exploration of the interaction networks among different biomolecules, which contained protein, gene and other types of interactions.¹⁵ Herein, the DEGs were mapped to STRING database to assess the PPI and visualized by Cytoscape program. The cut-off criterion of combined score and node degree was set to >0.15 and ≥10, respectively. Recent study suggested that the newly proposed method MCC performs better than other 10 previous reported methods.¹⁶ Therefore, the plug-in of Molecular Complex Detection (MCODE) was utilized to screen hub gene modules with the degree cut-off, haircut on, k-core, node score cut-off, max depth set as 10, 0.2, 2, 0.2 and 100. In addition, the functional and pathway enrichment analyses of DEGs were carried out by DAVID database in each module and the cut-off criterion of P-value was set as <0.05.

Prediction of miRNA targets and prognosis analysis

The target genes of the DEMs from GSE40744 and GSE36915 were predicted by miRNet webserver (https://www.mirnet.ca/), which is commonly used for predicting microRNA targets.¹⁷ The target genes were superposed with the DEGs to gain an intersection for further analysis. To identification of the hub genes, the results from MCODE, CytoHubba and miRNA-gene networks were combined. The Kaplan–Meier plotter webserver was applied to analyze the prognostic significance of the identified hub genes.¹⁷

Results

Identification of DEGs in expression profiles

GSE84006 datasets contained 38 tumor tissues and 38 paired adjacent non-tumor tissues. GSE14811 datasets contained 56 HCC samples and 56 liver samples. GSE14323 included 19 normal samples, 17 cirrhosis HCC samples, 41 cirrhosis samples and 38 HCC samples. The miRNA expression profile of GSE40744 included 9 normal samples and 67 HCV-associated HCC tumor samples. GSE36915 included 21 non-tumor samples and 68 tumor samples. A total of 8,655, 11,544, and 1,375 DEGs were identified after the analyses of the GSE84006, GSE14323, and GSE14811 datasets, respectively. Among them, 592 genes were found in all these datasets (Figure 1A).

Figure 1 Identification of overlapping DEGs and DE miRNAs. (A) Identification of overlapping upregulated and downregulated DEGs in GSE84006, GSE14323 and GSE14811. (B) Identification of overlapping upregulated and downregulated DE miRNAs in GSE36915 and GSE40744.

Functional and pathway enrichment analyses of DEGs

The GO analysis indicated 366 overlapping genes that are involved in numerous biological processes (BP), such as substance metabolic, single organisms and macromolecules (Figure 2A). In terms of cellular components (CC), DEGS were mainly enriched in the membrane-bounded organelles and intracellular membrane-bounded organelles (Figure 2B). The overlapping DEGs were primarily associated to the structural organization of the extracellular matrix, organic cyclic compounds binding, protein binding and heterocyclic compound binding in terms of molecular functions (Figure 2C). Additionally, the DEGs were enriched in three Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, which included the metabolic pathways, protein processing in the endoplasmic reticulum and the thyroid cancer pathway (Figure 2D).

Figure 2 GO method analyzed the overlapping DEGs. The number of DEGs were showed of black bars. (A) Biological processes (BP) of top 10; (B) Cellular components (CC) of top 10; (C) Molecular function (MF) of top 10; (D) KEGG pathway of top 10.

Construction of PPI network and analysis of modules

A total of 532 nodes and 3,370 edges were mapped in the PPI network of the identified DEGs. As shown in Figure 3, these DEGs-overlapped showed significant interactions and networks. The combined scores higher than 0.4 in PPIs were used for constructing the PPI networks. The MCODE plugin was used to analyze the entire PPI network and four modules were selected. As shown in Figure 4, enrichment analysis of the module genes of the KEGG pathway reveals enrichment in pathogenic Escherichia coli infection and the regulation of actin cytoskeleton. The first 50 genes provided by MCC method were selected and sequentially ordered as follows: SNRPD2, SF3B4, HNRNPA1, SNRPC, HNRNPU, RNPS1, POLR2H, SF3A2, PRPF3, POLR2I and so on (Table 1 ranked by MCC method).

Figure 3 PPI network construction and module analysis. DEGs used blue nodes to represent. Orange colors of nodes represented the correlation to modules with other DEGs; meanwhile, the lines represented the relationship of two nodes.

Figure 4 KEGG pathway analysis. The most important four modules analyzed by KEGG pathway. Bars signified the number of DEGs. Different colors represented different modules.

Table 1 Top 50 in PPI network ranked by MCC method

miRNA-DEG pairs

As shown in Figure 1B, a total of 22 DE miRNAs were identified by evaluating the dataset of miRNA expression dataset (GSE36915 and GSE40744). Then, 180 nodes and 224 edges were mapped in the miRNA-gene network of the identified DEGs and miRNAs (Figure 5). Moreover, the genes predicted by the miRNet webserver were considered as DE miRNAs target genes. For further identification of reliable hub genes, the identified target genes were compared to the DEGs. In this process, only the overlapping genes were considered as the hub genes. After combination of the results provided by MCODE, CytoHubba and the miRNA-gene network, ten hub genes were selected, including PLK1, PRCC, PRPF4, PSMA7, CHERP, GPKOW, HNRNPA1, PSMB5, PSMC4 and RNPS1.

Figure 5 miRNA-gene network. Regulation of DEGs in miRNA-gene network. Light blue nodes stand for DEGs, light orange nodes represent DE miRNAs. The lines represent the regulation of relationship between two nodes.

Kaplan–Meier survival analysis

To verify the results from bioinformatics analysis, the Kaplan–Meier plotter database was performed to predict the prognostic value of these hub genes. As shown in Figure 6, the OS for 364 patients with HCC was analyzed using the Kaplan–Meier survival plot. As a result, PLK1, PRCC, PRPF4 and PSMA7 were significantly correlated with worse OS for HCC patients.

Figure 6 Survival rates of four target genes. Using the Kaplan–Meier plotter database to analyze prognostic significance of the target genes in individuals of HCC. (A) PLK1. (B) PRCC. (C) PRPF4. (D) PSMA7. The red lines signified individuals with high expression of target gene and black lines with low expression.

Discussion

Although HCC has been extensively studied, HCC is still a seriously difficult problem in the early diagnosis. The reason for that is the lack of depth understanding of the molecular mechanisms in the occurrence, progression and development. Therefore, a depth research, relating to the key genes and mechanisms of HCC development and progression, is important for HCC diagnosis and pre-treatment. It is convenient and easy to discover the general genetic alterations in the development and progression of diseases by the novel microarray technology, which could investigate the underlying hub genes and mechanisms for the diagnosis and therapies.^18–20

In this study, we analyzed 592 overlapping DEGs, both up-regulated and down-regulated genes, in the GSE84006, GSE14323 and GSE14811 expression profile datasets (Figure 1). The GO analysis indicated that the 592 overlapping DEGs were correlated to substance metabolic, single organisms and macromolecules at the level of BPs. Besides, we tested DEGs in the KEGG pathway enrichment analysis. The results showed that the overlapping DEGs were mostly enriched in the metabolic pathway, protein processing in the endoplasmic reticulum and the thyroid cancer pathway. The metabolic pathway was identified as one of the essential pathways of the major modules of the overlapping DEGs (Figure 2). These enriched pathways may provide potential strategies for the development of new therapeutic agents.

Carcinogenesis is a complicated and complex process influencing with some specific genes and involving multiple signaling cascades. Recently, bioinformatics analyses have been used for identifying novel diagnosis markers and therapeutic targets for many cancers.^21–23 Zou et al²⁴ identified NRAGE as a novel biomarker for HCC by integrated bioinformatics analysis. Furthermore, in early HCC, the G2/M checkpoint participated in the progression.²⁵ Several similar researches have been shown for HCC now. Using bioinformatics analysis, Wang et al²⁶ have found the miR-192-3p-XIAP axis involved in HCC. However, compared with this study, previous studies only chose a signal data or file alone, as well as analyzed the one module method, with a high degree of connectivity, to select those hub genes. Notably, the genes provided by our study that not only combining the results from CytoHubba, MCODE method and network tools for the identification of the key genes, but also using the Kaplan–Meier database. This is forceful increasing the reliability of these results (Figures 3–5).

We predicted 10 hub genes containing PLK1, PRCC, PRPF4, PSMA7, CHERP, GPKOW, HNRNPA1, PSMB5, PSMC4 and RNPS1. A few studies have reported these genes in other cancers. To ensure the results obtained from the bioinformatics analysis, the Kaplan–Meier plotter database was been chosen to predict the prognostic value of these key genes. Only PLK1, PRCC, PRPF4 and PSMA7 significantly correlated with worse OS for the HCC patients (Figure 6).

PLK1, a serine/threonine-protein kinase that belongs to the polo-like kinase family, participates in various biological processes about embryonic development, cell cycle and RNA processing.²⁷ Xu et al²⁸ tested high expression of PLK1 in the HCC tissues showed significantly worse effect in the histological type. PLK1 and HOTAIR accelerate proteasome degradation of SUZ12 and ZNF198 during hepatitis B virus-induced liver carcinogenesis.²⁹ PLK1 phosphorylation of PTEN causes a tumor-promoting metabolic state.³⁰

PRCC, the most common TFE3 fusion partner in papillary renal carcinoma, is associated with pre-mRNA splicing factors.³¹ In a subset of papillary renal cell carcinomas, PRCC has been repeatedly reported and presumably causes cancer.^32,33 By RNA sequencing and FISH experiment, the clinicopathological analysis and detection of the gene fusion verified that the PRCC-MITF gene fusion defines a novel member of the MiT family translocation renal cell carcinoma.³³

PRPF4, a key component of the spliceosomes, plays a critical role in pre-mRNA splicing and its mutations result in retinitis pigmentosa due to photoreceptor defects.^34,35 PRPF4 is necessary part for cell survival, and in zebrafish, PRPF4 participate in posterior lateral line primordium migration.³⁶ However, PRCC and PRPF4 are not reported with HCC, and studying these important genes can therefore increase the understanding of not HCC but other cancers. Notably, it is not coincident in the PRPF4 mRNA level and its corresponding miRNA level. That result potential showed that miRNA may be not directly control PRPF4.

The gene PSMA7 encodes a member of the peptidase T1A family, which is a 20S core alpha subunit.³⁷ PSMA7 is often over-expressed in several human common cancers. In addition, overexpression of PSMA7 correlates with reduced survival rate in several types of cancers. Holmila et al³⁸ tested a diagnostic use of PSMA7 to screen for early HCC. Tan et al³⁹ investigated that PSMA7 inhibited the tumorigenicity of A549 human lung adenocarcinoma cells. The effects of sh-PSMA7 on cell proliferation and vascular endothelial growth factor expression via the ubiquitin-proteasome pathway in cervical cancer.⁴⁰ Salivary exosomal PSMA7 is a promising biomarker of inflammatory bowel disease.⁴¹

Conclusion

In conclusion, a total of four genes including PLK1, PRCC, PRPF4 and PSMA7 were identified as HCC biomarkers, and metabolic pathway, protein processing in the endoplasmic reticulum and thyroid cancer pathway were confirmed to be the essential and most mechanisms in the development of HCC. Our study, however, had some obvious limitations such as no analysis of clinic liver tissue samples, or only using the Kaplan–Meier plotter database to analyze the prognostic value of these molecular target genes. Further studies should explore these markers in the HCC patients to ensure the prognostic effect and need to dig the underlying mechanisms or related pathways. In the next step, we would try to use computer to predict which compound could target these hub genes to inhibit HCC development; meanwhile, using some biological method we could investigate these potential functions in vivo and in vitro. Overall, our results showed that the interactive network of miRNAs and mRNA is highly complex. However, we need more experiments to investigate these novel, key and hub genes.