Bioinformatics analysis to reveal key biomarkers for early and late hepatocellular carcinoma.

BACKGROUND: The aim of this study was to find the long non-coding RNAs (lncRNAs) and mRNAs acting as biomarker of early and late hepatocellular carcinoma (HCC).
METHODS: The Cancer Genome Atlas (TCGA) dataset was used to identify shared differentially expressed lncRNAs (DElncRNAs) and mRNAs (DEmRNAs) between early and late HCC and normal tissue. Functional annotation and protein-protein interaction network of shared DEmRNAs were performed. Furthermore, DElncRNAs-DEmRNAs co-expression network of early and late HCC were also performed. The expression of selected candidate genes were validated by the quantitative real time polymerase chain reaction (qRT-PCR).
RESULTS: A total of 1,201 shared DEmRNAs and 162 shared DElncRNAs were identified in both early and late HCC compared with normal controls. Cell cycle, p53 signaling pathway, retinol metabolism and metabolism of xenobiotics by cytochrome P450 were four significantly enriched pathways. Base on the protein-protein interaction network, CDK1, AURKA, CDC20, PLK1, AURKB, HIST1H2BG, BUB1B, CCNA2, CCNB1 and CDT1 were key protein. CTD-2510F5.4 and HAND2-AS1 were hub lncRNAs in both early and late HCC. Overall, the confirmation results of qRT-PCR were generally consistent with our integrated analysis.
CONCLUSIONS: A total of 4 DEmRNAs (CDK1, KIFC1, CENPF and ECM1) and 2 DElncRNAs (CTD-2510F5.4 and HAND2-AS1) were identified as key biomarkers for HCC. This study may contribute to reveal the pathogenesis of early and late HCC and provide new and accurate therapeutic targets for HCC. 2020 Translational Cancer Research. All rights reserved.

Introduction

Hepatocellular carcinoma (HCC) is a usual malignant tumor with an increasing incidence and one of the leading causes of cancer-associated deaths globally (1,2). Viral infections are often associated with liver cirrhosis and HCC (3). The therapeutic strategies of HCC depend on disease stage, operative tolerance, and possibility of liver transplant (4). HCC remains a challenging clinical problem with high incidence, limited therapeutic strategies, and poor prognosis. Therefore, it is a pressing task to screen potential biomarkers and uncover the detailed molecular mechanisms underlying of HCC.

With the advent of transcriptome analysis, bioinformatics have become most common ways to identify potential biomarkers of various diseases (5,6). Increasing evidences have indicated that lncRNAs plays a key roles in tumor-related biological processes, such as proliferation, invasion and metastasis (7,8). Recently, many lncRNAs are considered to be novel potential biomarkers and therapeutic targets for various cancers (9,10). The emerging roles of lncRNAs in HCC have also been explored in previous studies (11-13). However, there are few studies on biomarkers of lncRNA in early and late HCC.

In the current study, we identify shared differentially expressed lncRNAs (DElncRNAs) and mRNAs (DEmRNAs) between early and late HCC and normal tissue base on the Cancer Genome Atlas (TCGA) data. For the shared DEmRNAs, functional annotation and protein-protein interaction network construction were performed. Furthermore, DElncRNAs-DEmRNAs co-expression network of early and late HCC were performed. The expression levels of candidate genes were verified using quantitative real time polymerase chain reaction (qRT-PCR). Our study identified potential key DEmRNAs and DElncRNAs in early and late HCC, and provided novel insights into the mechanisms of lncRNAs in HCC.

Methods

Data collection

The lncRNA and mRNA gene expression profiles and corresponding clinical data of HCC were download from TCGA (http://tcga data.nci.nih.gov/). Here, 257 early HCC tissues, 65 late HCC tissues and 50 normal adjacent samples from patients with HCC were enrolled.

Differential expression analysis

The undetectable lncRNAs and mRNAs (with read count value =0 in more than 20% HCC case) were filtered and deleted. The DElncRNAs and DEmRNAs in HCC compared to adjacent normal tissues were identified using the R package DESeq2. The thresholds for differential expression analysis were defined as P value <0.01 and |Log2fold change|>1.5. We used the R package v3.3.3 (https://www.r-project.org/) to perform hierarchical clustering analysis of DElncRNAs and DEmRNAs.

Functional annotation

To discover the biological functions and the potential pathways of DEmRNAs in HCC, Gene Ontology (GO) classification and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were retrieved using DAVID 6.8 (https://david.ncifcrf.gov/). The terms with P value < 0.05 were significant results.

Protein-protein interaction network construction

The protein-protein interaction network was predicted using the Biological General Repository for Interaction Datasets (http://thebiogrid.org/), and then the protein-protein interaction network was presented by the Cytoscape (3.7.1) (http://www.cytoscape.org/). The nodes represent proteins and edges connect the nodes to show their relationship.

DEmRNA-DElncRNA interaction analysis

The correlation between the lncRNAs and DEmRNAs were analyzed by the Pearson correlation coefficient. The threshold for DElncRNA-DEmRNA co-expression pairs was P value <0.05 and |r|>0.80. We used the Cytoscape (3.7.1) to build the DElncRNA-DEmRNA co-expression network.

qRT-PCR verification

Three gene (CENPF, KIFC1 and ECM1) and three lncRNAs (CTD-2510F5.4, AC092171.4 and HAND2-AS1) were selected as candidate genes for qRT-PCR verification. Here, 6 HCC patients were included in our study from the Department of Hepatobiliary Surgery, The Third Hospital of Hebei Medical University. Twelve tissue samples including 6 HCC tumor samples and 6 cancer adjacent samples were obtained. This study was approved by the institutional ethics committee of The Third Hospital of Hebei Medical University (No. 2018-025-1) and conducted in accordance with the Declaration of Helsinki. The signed informed consents of all the participants were obtained.

Total RNA was isolated using TRizol reagent kit (Invitrogen, China). SuperScript III Reverse Transcription Kit (Invitrogen, China) was utilized to obtain the complementary DNA. With Super Real PreMix Plus SYBR Green (Invitrogen, China), quantitative real-time PCR were generated using the ABI 7500 system (Applied Biosystems, USA). ACTB and GADPH acted as the internal criterion for the targeted genes. The 2−ΔΔCt method was used to address the data.

Results

Identification of DEmRNAs and DElncRNAs

After data processing, we screened 1,306 DEmRNAs (934 up-regulated DEmRNAs and 372 down-regulated DEmRNAs) and 193 DElncRNAs (149 up-regulated DElncRNAs and 44 down-regulated DElncRNAs) in early HCC samples compared with normal samples. A total of 1,740 DEmRNAs (1,222 up-regulated DEmRNAs and 518 down-regulated DEmRNAs) and 265 DElncRNAs (192 up-regulated DElncRNAs and 73 down-regulated DElncRNAs) were obtained between late HCC and normal samples. demonstrated hierarchical clustering analysis of the top 100 DEmRNAs and DElncRNAs between early HCC and normal tissue, respectively. showed hierarchical clustering analysis of the top 100 DEmRNAs and DElncRNAs between late HCC and normal tissue, respectively. Moreover, a total of 1,201 shared DEmRNAs and 162 shared DElncRNAs were identified in both early and late HCC compared with normal controls. was Venn diagram of DEmRNAs and DElncRNAs in early and late HCC.

Figure 1

Hierarchical clustering analysis of top 100 differentially expressed lncRNAs (DElncRNAs) and mRNAs (DEmRNAs) between hepatocellular carcinoma (HCC) and normal tissues. (A) DEmRNAs in early HCC. (B) DElncRNAs in early HCC. (C) DEmRNAs in late HCC. (D) DElncRNAs in late HCC. Row and column represented DElncRNAs/DEmRNAs and tissue samples, respectively. Orange and light blue color mean the (early and late) HCC and normal control, respectively. The color scale represented the expression levels. Red color represents the relative expression level of genes was higher than mean, and green color represents the relative expression of genes was lower than mean.

Figure 2

Venny analysis of DEmRNAs and DEmRNAs in early HCC vs. normal control and late HCC vs. normal control. (A) DEmRNAs; (B) DElncRNAs. DEmRNAs, differentially expressed mRNAs; DElncRNAs, differentially expressed long non-coding RNAs; HCC, hepatocellular carcinoma.

Functional annotation of shared DEmRNAs in early and late HCC

A total of 1,201 shared DEmRNAs in early and late HCC were used to perform the the GO and KEGG pathway enrichment analysis (). GO enrichment analysis displayed that M phaseon (P=8.55E-25), cell cycle phasee (P=2.01E-24), chromosome, centromeric region (P=7.50E-13) and calcium ion binding (P=8.65E-06) were significantly enriched GO terms. KEEG enrichment analysis displayed that cell cycle (P=1.92E-09), p53 signaling pathway (P=7.50E-06), Retinol metabolism (P=4.20E-05) and metabolism of xenobiotics by cytochrome P450 (P=5.35E-04) were four significantly enriched pathways.

Figure 3

Top 15 significantly enriched Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways for shared differentially expressed mRNAs (DEmRNAs) in early and late hepatocellular carcinoma (HCC). The x-axis shows −log P value, and the y-axis shows GO terms or KEGG pathways. (A) Biological process; (B) cellular component; (C) molecular function; (D) KEGG pathways. The change color from red to green indicates that the significance is high to low. The larger the dot size, the greater the proportion of this enrichment item.

Protein-protein interaction network of shared DEmRNAs in early and late HCC

The protein-protein interaction network was consisted of 292 nodes and 504 edges (). Top 10 proteins with higher degrees, including CDK1 (degree=35), AURKA (degree=29), CDC20 (degree=26), PLK1 (degree=23), AURKB (degree=21), HIST1H2BG (degree=17), BUB1B (degree=17), CCNA2 (degree=17), CCNB1 (degree=16) and CDT1 (degree=15).

Figure 4

Protein-protein interaction network. Ellipses are used to represent nodes, and lines are used to represent edges. Red and green represent up- and down-ward adjustments, respectively. The blue border indicates the top 20 up- and down-regulated proteins.

DElncRNAs-DE mRNAs co-expression network of early HCC

A total of 87 DElncRNA-DEmRNA co-expression pairs including 18 DElncRNAs and 85 DEmRNAs were identified with absolute value of P value <0.05 and |r|>0.80. The co expressed DElncRNA-DEmRNA network () was consisted of 103 nodes and 87 edges, and its hub lncRNAs were CTD-2510F5.4 (degree=58) and HAND2-AS1 (degree=11).

Figure 5

DElncRNAs-DEmRNAs co-expression network of early HCC. Ellipses and rhombus represent DEmRNAs and DElncRNAs, respectively. Red and green colors represent up- and down-regulation, respectively. The blue border indicates the top 20 up- and down-regulated DElncRNAs and DEmRNAs. DEmRNAs, differentially expressed mRNAs; DElncRNAs, differentially expressed long non-coding RNAs; HCC, hepatocellular carcinoma.

DElncRNAs-DEmRNAs co-expression network of late HCC

A total of 2,339 DElncRNA-DEmRNA co-expression pairs including 92 DElncRNAs and 389 DEmRNAs were identified with absolute value of P value <0.05 and |r|>0.80. The co-expressed DElncRNA-DEmRNA network () was consisted of 481 nodes and 2,339 edges, and its hub lncRNAs were CTD-2510F5.4 (degree=225), RP11-196G18.22 (degree=149), SNHG1 (degree=146), RP5-821D11.7 (degree=146), TMPO-AS1 (degree=144) and HAND2-AS1 (degree=14).

Figure 6

DElncRNAs-DEmRNAs co-expression network of late HCC. Ellipses and rhombus represent DEmRNAs and DElncRNAs, respectively. Red and green colors represent up- and down-regulation, respectively. The blue border indicates the top 20 up- and down-regulated DElncRNAs and DEmRNAs. DEmRNAs, differentially expressed mRNAs; DElncRNAs, differentially expressed long non-coding RNAs; HCC, hepatocellular carcinoma.

As displayed, CTD-2510F5.4, HAND2-AS1 and ECM1 were down-regulated and AC092171.4, CENPF and KIFC1 were up-regulated in HCC compared to adjacent tissues. Base on the TCGA results, HAND2-AS1 and ECM1 were down-regulated while the other lncRNAs and DEmRNAs were up-regulated in HCC compared to adjacent tissues. In generally, the confirmation results of qRT-PCR were consistent with our integrated analysis.

Figure 7

Validation optimal DEmRNAs and DElncRNAs in HCC tissue by qRT-PCR. All of the assays were performed three times independently at least. The x-axis shows DEmRNAs or DElncRNAs and y-axis shows log2(Foldchange). The log2(Foldchange) >0 and log2(Foldchange) <0 indicates up- and down-regulation, respectively. *P<0.05, **P<0.01, ***P<0.001. DEmRNAs, differentially expressed mRNAs; DElncRNAs, differentially expressed long non-coding RNAs; HCC, hepatocellular carcinoma.

Discussion

HCC is a prevalent malignancy that is considered the leading cause of cancer-related deaths worldwide due to the lack of effective treatment (3). Thence, it is urgent to discover potential biomarkers and therapeutic targets for HCC. In this study, a total of 1201 shared DEmRNAs and 162 shared DElncRNAs were identified in both early and late HCC compared with normal controls. For the shared DEmRNAs, functional annotation and protein-protein interaction network construction were performed. Among which, CDK1 (degree=35) was key nodes. In addition, the expressions of 6 candidate genes were verified by the qRT-PCR experiments.

Cyclin dependent kinase 1 (CDK1), belongs to the Ser/Thr protein kinase family, is associated with cell cycle. A recent study has been reported that metformin inhibits the proliferation of HCC cells through inducing G2/M arrest and reduces the expression of CDK1 (14). Wu et al. found that high expression levels of CDK1 in HCC were associated with poor overall survival (15). CDK1 might be involved in progression of HCC via interaction with PLK1 and FOS (16). In this study, CDK1 was also up-regulated in HCC and significantly enriched in pathways of p53 signaling pathway. Paired-related homeobox 1 involves in metastasis and apoptosis of HCC cells by regulating p53-dependent signaling pathway (17). Therefore, we speculated that CDK1 might play a role in pathogenesis of HCC by regulating p53-dependent signaling pathway.

Kinesin family member C1 (KIFC1), a member of the kinesin superfamily, is to act as a motor protein and move towards the minus-end of the microtubule to further affect its dynamics and functions (18). KIFC1 regulates the formation of centrosomes, and it has been found that upregulation of KIFC1 results in the formation of unipolar spindles in cancer cells (19). It has been reported that KIFC1 is a potential marker that predicts an aggressive disease course in ovarian cancer (19). KIFC1 is overexpressed in HCC and that its expression is closely related to the progression and metastasis of HCC (20). High expression of KIFC1 is associated with the progression of HCC, indicating that KIFC1 is a putative marker for prognostic and therapeutic in HCC (21). In this study, we found that KIFC1 was up-regulated in both TCGA integration analysis and qRT-PCR validation, indicating that it might be a potential therapeutic target for HCC.

Centromere protein F (CENPF) is a cell cycle-associated nuclear antigen associated with the centromere-kinetochore complex whose expression and subcellular localization depend on the cell cycle and degrade rapidly at the end of cell division (22). CENPF is defined as a biomarker of cell proliferation in various malignant tumors, including breast cancer (23), HCC (22), and other tumors (24). Researchers have reported that CENPF may be a biomarkers in diagnosis and prognosis of HCC (22,25). In the present study, we found that CENPF was also up-regulated in both TCGA integration analysis and qRT-PCR validation. CENPF was co-expressed with CTD-2510F5.4 in both early and late HCC. CTD-2510F5.4 expression is positive correlation with prognosis in lung adenocarcinoma (26). CTD-2510F5.4 is a malignant phenotype associated lncRNAs and a biomarker for the diagnosis and prognosis of gastric cancer (27). Therefore, we hypothesized that CTD-2510F5.4 might be involved in progression of HCC by regulating CENPF.

Extracellular matrix protein 1 (ECM1), an 85 kDa glycoprotein, is involved in a number of biological processes including mineralization, cell proliferation and angiogenesis. ECM1 modulates the epithelial–mesenchymal transition progression and cell proliferation of breast cancer via regulating the association between β-catenin and the MUC1 (28). Chen et al. found that ECM1 accelerates migration and invasion of HCC via mediating epithelial-mesenchymal transition (29). In this study, we found that ECM1 was down-regulated in both TCGA integration analysis and qRT-PCR validation, and ECM1 was co-expressed with HAND2-AS1 in both early and late HCC. A study has suggested that HAND2-AS1 is down-regulated in HCC tissues, and HAND2-AS1 is associated with invasion of HCC cells (30). Here, HAND2-AS1 was reduced in HCC, which provided further evidence for previous study. Therefore, we speculated that HAND2-AS1 might be involved with the pathogenesis of HCC through regulating ECM1.

In summary, a total of 1,201 shared DEmRNAs and 162 shared DElncRNAs were identified in both early and late HCC compared with normal controls. Among which, CDK1, KIFC1, CENPF, ECM1, CTD-2510F5.4 and HAND2-AS1 were identified as potential biomarkers for HCC. Overall, the confirmation results of qRT-PCR were generally consistent with our integrated analysis. Several limitations of our study should be mentioned. More samples are needed to validate expression of DEmRNAs and DElncRNAs. Besides, in vivo and in vitro experiments are necessary to reveal biological functions of DEmRNAs and DElncRNAs in HCC in future study.