Identification of crucial genes based on expression profiles of hepatocellular carcinomas by bioinformatics analysis

Microarray data

The raw microarray data of GSE121248, including 70 HCC tissues and 37 adjacent normal tissues, were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/), which was an openly available database. Another two datasets of GSE25097 and GSE22058 were used as the testing sets to verify our results.

Data processing and identification of DEGs

The bioconductor package (version3.7, http://www.bioconductor.org/), affy package (version1.58.0, http://bioconductor.org/packages/affy/) and the RMA algorithm (Irizarry et al., 2003) were used for data preprocessing and normalization in R (version 3.5.1, R Core Team, 2018). The t test methods of the limma package (version 3.38.2, http://bioconductor.org/packages/limma/) (Ritchie et al., 2015) were used to identify DEGs in HCC tissues compared to adjacent normal tissues. Values of log Fold Change >2.0 or log Fold Change <−2.0 and P–value <0.05 served as the cut-off criteria.

Functional annotation of DEGs

Gene Ontology (GO) analysis was applied using the Database for Annotation, Visualization and Integrated Discovery (DAVID, version 6.8, https://david.ncifcrf.gov/) to explore the biological processes, cellular components and molecular function in which the DEGs were involved (Ashburner et al., 2000). GeneAnswers (version 2.24.0, http://bioconductor.org/packages/GeneAnswers/) was a Kyoto Encyclopedia of Genes and Genomes (KEGG)-based R package that not only provided the annotation of biochemistry pathways associated with diseases and drugs but also facilitated visualization of the signaling pathways enriched by DEGs (Feng et al., 2010).

Gene set enrichment analysis of HCC

Gene Set Enrichment Analysis (GSEA) (version 3.0, the broad institute of MIT and Harvard, http://software.broadinstitute.org/gsea/downloads.jsp) was performed between HCC and adjacent normal tissues to investigate the biological characteristic of HCC (Subramanian et al., 2005). In detail, the ‘collapse data set to gene symbols’ was set to true, the ‘permutation type’ was set to phenotype, the ‘enrichment statistic’ was set to weighted, and the Signal2Noise metric was used for ranking genes. GSEA calculated a gene set Enrichment Score (ES) that analyzed genes were enriched in the biological signal conduction on the MsigDB (Molecular Signatures Database, http://software.broadinstitute.org/gsea/msigdb). And the meandiv normalization method was used for enrichment scores in the gene sets. In this study, 1,000 gene of permutations were set to generate a null distribution for enrichment score in the hallmark gene sets and functional annotation gene sets. The gene sets database used for enrichment analysis were mainly including ‘h.all.v6.2.symbols.gmt’, ‘c2.cp.kegg.v6.2.symbols.gmt’, ‘c2.cp.reactome.v6.2.symbols’ and ‘c5.bp.v6.2.symbols’. Nominal P–value <0.05, FDR <0.25 and gene set size >100 were defined as the cut-off criteria.

PPI network and module analysis

The PPI network was analyzed using the STRING database for Homo sapiens (Szklarczyk et al., 2017). STRING was used to calculate PPI networks of DEGs with a combined score >0.4 as the cut-off criteria. Cytoscape (version3.7.0, https://cytoscape.org/) was used to visualize the network (Shannon et al., 2003). Mcode (version1.5.1, Bader Lab, University of Toronto, http://www.baderlab.org/Software/MCODE/) was a Cytoscape plugin for constructing the protein modulecomplex with a degree cut-off = 2, node score cut-off = 0.2, max.depth = 100, and k-core = 2. CytoNCA (version 2.1.6, http://apps.cytoscape.org/download/stats/cytonca/) was used to analyze centrality of protein interaction networks that were representative of potentially crucial proteins in the network (Tang et al., 2015). The crucial proteins were identified based on four different centrality parameters (degree centrality (DC), betweenness centrality (BC), eigenvector centrality (EC) and closeness centrality (CC)). The values of centrality parameters were sorted to select the top 10 ranked proteins, and the proteins were input as the crucial candidates. ClueGO (version 2.5.3, http://apps.cytoscape.org/apps/cluego) and CluePedia (version 1.5.3, http://apps.cytoscape.org/apps/cluepedia) were employed to analyze the KEGG and Reactome pathways.

Transcriptional expression level and survival analysis of crucial genes on HCC

The mRNA expression of crucial genes of HCC and adjacent normal tissues were evaluated based on the shinyGEO database, which included two HCC datasets GSE25097 (containing of 268 HCC and 243 adjacent normal tissues) and GSE22058 (containing of 96 HCC and 96 adjacent normal tissues) (Dumas, Gargano & Dancik, 2016). The Human Protein Atlas (http://www.proteinatlas.org) was used for protein detection by immunohistochemistry. We performed TCGA methylation data of MEXPRESS database (https://mexpress.be/) to explain why oncogene mRNAs were highly expressed in HCC. The patients overall survival was investigated using TCGA_LIHC data and the clinicopathological features was download by TCGA database (https://www.cbioportal.org/). Statistical significance was considered at P < 0.05.

HCC sample collection and revalidation of mRNA expression of crucial genes

Thirty pairs of HCC and paired adjacent tissues were obtained from Xijing Hospital, Fourth Military Medical University of China. Our study was approved by the Ethical Committee and Institutional Review Board of Fourth Military Medical University (Ethical Application Ref: KY20183288-1). The RNA was extracted using RNA extraction kit II (Omega Bio-tek, GA, USA) and cDNA was synthesized using PrimeScript™ RT reagent Kit (TaKaRaBio, Otsu, Japan). The mRNA expression was assessed in 384-well plates using the CFX connect™ Real-Time System (BIO-RAD, Hercules, CA, USA) with the SYBR Premix Ex Taq II (TaKaRaBio, Otsu, Japan) and primers, and qPCR data were analyzed using ^△△Ct method with β-actin as the reference gene. Primer sequences were listed in Table S1.

Statistical analysis

Statistical analysis was performed with SPSS 19.0 software (IBM, Armonk, NY, USA). All the data were presented as mean ± SEM. Overall survival time was calculated by the Kaplan–Meier method and analyzed with the log-rank test. The univariate and multivariate analyses were calculated based on the Cox regression model. The statistical analysis was performed as appropriate by Student’s t test and χ² test. The P-value < 0.05 were considered statistically significant.

Identification of DEGs

A total of 164 DEGs that met the cut-off criteria were identified, including 72 upregulated and 92 downregulated genes (Fig. 1A). The mRNA expression profiles of 30 representative DEGs in 70 HCC and 37 adjacent normal tissues were presented in Fig. 1B.

Functional annotation analysis of DEGs

GO analysis of DEGs revealed that the following GO terms were significantly enriched in condensed chromosome kinetochore (P = 1.33E-10), mitotic nuclear division (P = 7.45E-10), sister chromatid cohesion (P = 1.97E-09), and cell division (P = 2.34E-09) (Table S2). KEGG annotation analysis of DEGs in HCC revealed that the pathway mainly including retinol metabolism (P = 4.21E-05), drug metabolism-cytochrome P450 (P = 5.95E-05), cell cycle (P = 2.44E-04), and the p53 signaling pathway (P = 5.59E-04) (Fig. 2A).

Protein–protein interaction network of DEGs

The PPI network of DEGs with 113 nodes (genes) and 933 edges (interactions) was built by STRING. The CytoNCA plugin was used to analyze the centrality of nodes. The top 10 proteins, which were also defined as crucial proteins based on four different centrality parameters, were presented in Table 1, and they included TOP2A, NDC80, ESR1, ZWINT, CENPF, NCAPG, ENO3 and CCNB1. The majority of these crucial proteins participated in mitotic prometaphase (P = 1.5E-14) and resolution of sister chromatid cohesion (P = 6E-13) (Fig. 2B).

The MCODE plugin was employed to identify the modules in the PPI network. The top 2 significant modules were displayed in Fig. 2C, and the seed of two separate modules are RAD51AP1 and ESR1. KEGG and Reactome pathway analyses revealed that the two modules were mainly related to cell cycle checkpoints (P = 7.94E-15), RHO GTPase effectors (P = 6.87E-13), and Cytochrome P450 (P = 1.37E-12) (Table S3).

Gene Set Enrichment Analysis of HCC

To provide further insight into the whole gene enrichment annotation of HCC and consider the potential role of undifferentiated genes or genes with smaller differences in expression in HCC and adjacent normal tissues, GSEA was performed to explore the gene expression profile based on the GO terms and pathway databases of KEGG, Reactome and Hallmark gene sets. The result revealed that enrichment for three of six gene sets were significant at FDR <0.25, including genes involved in cell cycle and chromosome segregation (Fig. 3).

Validation of mRNA expression of crucial genes

Six significantly crucial genes, NDC80, ESR1, ZWINT, NCAPG, ENO3 and CENPF, were validated in two other HCC datasets (GSE25097 and GSE22058). The results also revealed that NDC80, ZWINT, NCAPG and CENPF were significantly upregulated in HCC, and ESR1 and ENO3 were downregulated in HCC (Fig. 4). Then, immunohistochemistry staining validated from the Human Protein Atlas database showed that ZWINT, NCAPG and CENPF protein expression were strongly upregulated in liver cancer tissues compared with normal tissues (Figs. 5A–5C). ENO3 and ESR1 protein were downregulated in liver cancer tissues compared with normal tissues (Figs. 5D & 5E). The protein expression of NDC80 was absent in the Human Protein Atlas database. Consistent with database analyses, the six crucial gene were successfully validated by qPCR in 30 paired human HCC and adjacent normal tissues, and the results revealed no difference with the analysis results of microarray profiling with the exception of ENO3 (Fig. 6).

The relationship between the methylation and mRNA expression of crucial genes were respectively validated by TCGA data to identify oncogenes with highly expressed mRNA in HCC. The results showed that mRNA expression of crucial genes was negatively correlated with methylation levels. The genes highly expressed or hypomethylated were mostly involved in HCC. Hypermethylated/poorly expressed crucial genes in HCC tended to exhibit high expression in adjacent normal tissues (Table 2).

Kaplan–Meier survival analysis of crucial genes in HCC patients

TCGA data were used to explore the relationship between the crucial genes and the survival of HCC patients, and high expression of NDC80 (P = 1.10E-3), NCAPG (P = 2.80E-4), CENPF (P = 5.72E-4) and ZWINT (P = 2.47E-4) significantly correlated with worse survival probability for HCC patients. ESR1 (P = 2.02E-4) exhibited the opposite trend, and ENO3 was not a significant marker of survival prognosis (P = 0.186) (Fig. 7). There we focused a crucial tumor suppressor gene, ESR1. The RNA expression data from the TCGA_LIHC cohorts revealed that increased ESR1 mRNA expression correlated significantly with gender, age, serum AFP, TNM stage, tumor recurrence and tumor differentiation (Table S4). Univariate analysis also revealed that the crucial genes ESR1 expression was significantly correlated with overall survival and recurrence (Table S5). Multivariate analysis was performed using all of the variables that were identified as significant by univariate analysis. The result showed that ESR1 mRNA expression was an independent prognostic indicator for overall survival and recurrence, and other features of patients with HCC were also analyzed (Table S5).

The mRNA expression of ER related receptor and regulatory genes in HCC

Since estrogen plays a protective role in liver cancer by binding to estrogen receptors, we also further investigated the expression of ER (estrogen receptor) regulon genes in HCC tissues compared with adjacent tissues and the correlation between ER regulon genes and ESR1 using the public GEO database (GSE22058). The mRNA expression levels of sex hormone- related genes (ESR2, ESRRB, ESRRG, PGRMC1, PGRMC2, GPER, SHBG) and ER regulatory gene, TFF1 in HCC tissues were significantly lower than those of adjacent tissues, whereas the expression of estrogen-related receptor gene, ESRRA and ER regulatory genes (MTA1 and FOXA1) were higher in tumor tissues (Fig. 8). In HCC tissues, ESR1 mRNA level was positively correlated with ESRRB (r = 0.3441, P = 0.0006), PGRMC1 (r = 0.5124, P < 0.0001) and GPER (r = 0.2274, P = 0.0259) (Fig. 9). Moreover, ESR1 mRNA level was negatively correlated with ER regulatory genes MTA1 (r =  − 0.2714, P = 0.0079) and TFF1 (r =  − 0.3557, P = 0.0004) (Fig. 9).