Integrative analysis of dysregulated lncRNA-associated ceRNA network reveals potential lncRNA biomarkers for human hepatocellular carcinoma

Literature metrology analysis method

All HCC-related literature was obtained from the Science Citation Index Expanded (SCI-E) from the Web of Science (WOS) of Clarivate Analytics on February 1, 2019. The documents were analyzed by two independent authors. The literature data retrieval strategy was as follows: title = (‘hepatocellular carcinoma’) or title = (‘hepatocellular cancer’) or title = (‘liver carcinoma’) or title = (‘liver cancer’) and title = (‘human’) and title = (‘biomarker’). All references were dated between 2009 to 2018 and only research articles and reviews in English were included. The data were obtained from the WOS and did not include animal studies.

HCC biomarker-related literature were collected from the WOS and analyzed using VOSviewer 1.6.5 software (Leiden University, Leiden, Netherlands) and CiteSpace V software (Drexel University, Philadelphia, PA, USA), respectively. VOSviewer 1.6.5 and CiteSpace V software were used to perform literature cluster analysis and key word hotspot analysis.

Patients and samples

We collected data from 349 patients with HCC and RNA sequencing data from TCGA dated up to November 1, 2018. Annotation information for the RNA sequencing datasets were obtained using Affymetrix Human Genome Array platforms. The study was in accordance with the TCGA database portal platform guidelines. The subjects were simplified based on the following exclusion criteria: (a) without completed data information; (b) histologic diagnosis was not HCC; and (c) two or more types of cancers, including HCC. 313 HCC tumors and 44 normal liver tissue samples were included in this study. Of the 313 HCC patients, 233 patients had histopathological stage I/II HCC, and 80 had stage III/IV HCC, according to the 7th American Joint Committee on Cancer (AJCC) Tumor Node Metastasis (TNM) staging system. 249 of these cases had lymphatic metastasis while 64 did not. Details of the RNA sequencing datasets from TCGA, sample descriptions, and clinicopathological features are provided in Table 1. The flow diagram for integrated bioinformatics analysis from TCGA is shown in Fig. 1.

Integration of RNA sequence data and differential expression analysis

RNA level 3 expression data were processed and standardized based on the mRNA expression data of TCGA. The original RNA sequencing raw reads were processed and normalized using the TCGA RNASeqV2 system to fit the analysis. HCC level 3 normalized miRNAs sequencing data (Illumina HiSeq 2000 microRNA sequencing platforms) (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga/using-tcga/types) were also downloaded from TCGA. We analyzed and contrasted the significant dysregulated lncRNAs, mRNAs, and miRNAs in tumor tissues from 313 HCC patients and 41 normal liver tissues using limma R package tool (false discovery rate (FDR)<0.05, fold change>2, P < 0.05). We used overlapping subclass analysis to identify the separate dataset co-differentially expressed genes, including lncRNAs, miRNAs, and mRNAs using Venn 2.1 software (http://bioinfogp.cnb.csic.es/tools/venny/index.html). The lists of integrated HCC tissue dysregulated lncRNAs, miRNAs, and mRNAs were saved for further analysis.

Gene functional enrichment analysis

Significantly dysregulated intersected mRNAs were selected and imported into the Gene Ontology (GO) tool (http://www.geneontology.org) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.kegg.jp/) to identify their molecular function and to find potential regulated signal pathways for these genes. Up-regulated and down-regulated mRNAs from the overlapping subclasses were analyzed. Visualization of the GO and KEGG were plotted using R software.

Construction of the ceRNA network

The lncRNA-miRNA-mRNA ceRNA network was built based on the theory that lncRNA can bind miRNA, acting as so-called miRNA sponges, with miRNA binding to the mRNAs and negatively regulating gene expression. HCC tissues significantly dysregulated lncRNAs, mRNAs, and miRNAs (FDR<0.05, fold change >2, P < 0.05) and were selected to build the ceRNA network, in which the fold changes of genes were rooted in the TCGA database, to determine whether these intersection genes were involved in ceRNA regulation. MiRcode (https://omictools.com/mircode-tool), miRanda (http://www.microrna.org/microrna/home.do) and Targetscan (http://www.targetscan.org/) were used to predict the miRNA target lncRNAs and miRNA-mRNAs interactions in the different databases. Our study combined the significant function regulation genes in GO and KEGG, and miRNAs predicted target genes to further assist in the selection of the intersection mRNAs. The subset of intersection miRNAs were selected to negatively regulate lncRNAs and to assist in the selection of the intersection mRNAs used to build the ceRNA network according to the ceRNA regulation theory. We used Cytoscape software 3.0 (National Institute of General Medical Sciences, Bethesda, MD, USA) for this analysis.

Analysis of the association between ceRNA network key lncRNAs and HCC clinical features from TCGA

Abnormally expressed key lncRNAs may play an important role in HCC progression based on the lncRNAs in the ceRNA network. The key lncRNAs involved in the network were selected as target lncRNAs that may be associated with HCC progression. We explored the potential association between ceRNA network key lncRNAs and the TCGA clinical features of HCC patients, which included race, gender, TNM stage, tumor grade, lymphatic metastasis, and chronic HBV or HCV infection, using multiple linear regression analysis.

Kaplan–Meier survival curve analysis

Kaplan–Meier survival analysis was performed to investigate whether the expression of ceRNA network key lncRNAs was associated with the overall survival of HCC patients. Kaplan–Meier survival analysis parameters were calculated using the publicly available TCGA HCC patient datasets and Gene Expression Profiling Interactive Analysis (GEPIA) tools (http://gepia.cancer-pku.cn/). The survival distributions of patients with HCC in TCGA, and the key lncRNAs expression level changes were analyzed using Kaplan–Meier, log-rank, and hazard ratio (HR). P < 0.05 was the cutoff criterion.

Preparation for human HCC samples and qRT-PCR validation

Samples from tumor tissue and paired non-tumor liver tissue were collected from 20 HCC patients (aged 40–69 years) at Lanzhou University Second Hospital (Lanzhou, China), for qRT-PCR validation. Patients were diagnosed with HCC according to their histopathology. All patients provided informed consent and their clinical information was collected by an investigator using patient interviews and medical records. The collection of the tumor samples from HCC patients was approved by the School of Public Health, Lanzhou University (Lanzhou, China) (Lzuggwsxy-20190806) and conformed to the Helsinki Declaration and current legislation. Samples were collected and stored in RNAlater (Ambion, Foster City, CA, USA) at −80 C.

The total RNA from the tissue samples was isolated using the TRIzol® reagent (Invitrogen; Thermo Fisher Scientific, Inc. Carlsbad, CA, USA). The Reverse Transcription Kit (Promega Corporation, Madison, WI, USA) and GoTaq® qPCR Master Mix of Power SYBR® Green (Promega Corporation) were used to synthesize cDNA and for qRT-PCR detection. qRT-PCR was performed using the Step One PlusTM PCR System (Applied Biosystems; Thermo Fisher Scientific). qRT-PCR relative fold change results were calculated using the 2^−ΔΔCt method.

Statistical analyses

Data were analyzed using SPSS Statistics V21.0 (IBM, Armonk, NY, USA) and expressed as mean ± SD. All analyses were performed three times and represent data from three individual experiments. A two-tailed Student’s t-test was used to measure the significance of differences between subgroups. Kaplan–Meier survival analysis was used to investigate the correlation between the changes in lncRNAs expression levels and the prognostic overall survival times for patients. Statistical significance was P < 0.05.

Literature metrology analysis of HCC research

922 publications from 2008 to 2018 matched the search criteria. These HCC-biomarker-related studies were analyzed by VOSviewer and three primary clusters were identified: pathogenesis related, clinical patients related, and etiologically related. Cluster analysis showed that there were three major focuses in HCC research (Fig. 2A).

Key words and article titles from the 922 papers were analyzed using VOSviewer software. The integrated analysis is shown in Fig. 2B; colors were assigned to the key words by VOSviewer. The different color shades represent the usage frequency of the key words; the colors, ranging from blue to yellow, represent a low to high frequency of occurrence, respectively. Key words with yellow (high frequency) represented the research hotspots in this field. The key words analysis revealed that “prognosis”, “effect”, “expression”, and “concentration” were frequently used.

The high frequency key words were identified by CiteSpace V software analysis as the frontier research fields. One of these frontier research keywords was “lncRNA” (Fig. 3). More recently dated studies including the keyword “lncRNA” with greater frequency. Based on these results, we determined that the objectives of our study were to find the relationships between the expression levels of lncRNA and HCC progression, and to identify the potential diagnostic and prognostic biomarkers for this disease. We utilized the TCGA database HCC-related RNA sequence data mining as the data source for comparisons of gene differences and bioinformatics analysis.

HCC-specific lncRNAs

We found that 323 lncRNAs were significantly dysregulated in HCC tumor tissues from the TCGA database (fold change>2, P < 0.05). The 323 significantly dysregulated lncRNAs in HCC with different tumor stages and lymph node metastasis status were carefully analyzed. 231 lncRNAs were found to be significantly dysregulated in HCC tumor stage I/II (non-lymphatic metastasis) compared with normal liver tissues; 208 lncRNAs were significantly dysregulated in HCC tumor stage I/II (lymphatic metastasis) compared with normal liver tissues; 252 lncRNAs were significantly dysregulated in HCC tumor stage III/IV (non-lymphatic metastasis) compared with normal liver tissues; and 199 lncRNAs were significantly dysregulated in HCC tumor stage III/IV (lymphatic metastasis) compared with normal liver tissues (Fig. 4A). We selected the 128 intersected lncRNAs, which including 85 upregulated and 43 downregulated genes, for further analysis and construction of the ceRNA network (Table S1).

Function analysis of intersected mRNAs

We found that 2,026 HCC tissues had significant differences in mRNA expression and were included in the Venn diagram intersection subset analysis. These differentially expressed genes may play key roles in the progression of HCC. We analyzed the potential biological regulatory functions of these 2,026 mRNA by GO enrichment of functions and KEGG pathway analyses. The most enriched function by GO analysis of upregulated mRNAs was the ‘Mitotic cell cycle’ (Fig. 5A). The most enriched function by GO analysis of downregulated mRNAs was ‘Small molecule metabolic process’ (Fig. 5B). KEGG pathway analysis indicated that 60 signaling pathways were involved in regulation by upregulated mRNAs, and the most enriched pathway was ‘Cell cycle’ (Fig. 6A). In addition, there were 152 signaling pathways involved in regulation by downregulated mRNAs, and the most enriched pathway was ‘Metabolic pathways’ (Fig. 6B). The MAPK signaling pathway has been shown to participate in the progression of HCC (Zhao et al., 2018), and the P53 signaling pathway is a key pathway in HCC cell proliferation and apoptosis (Zhao et al., 2019b). Bladder cancer, small cell lung cancer, pathways in cancer, and the PI3K-AKT signaling pathway may also be involved in the regulation of cancer progression (Tang et al., 2019).

Prediction of miRNA targets and construction of ceRNA network

168 miRNAs were found to be significantly dysregulated between HCC tissues and normal liver tissues (fold change>2, P < 0.05). We selected the intersected subset of 76 miRNAs related to HCC tumor histological type and lymphatic metastasis (Fig. 4C). We predicted the potential relationships between these 76 miRNAs and the above intersected subset of 128 lncRNAs (Fig. 4A) by miRanda software. There were 59 specific miRNAs interacting with 92 specific lncRNAs (Table S2).

The ceRNA network was constructed based on the predicted miRNA-targeted genes. We predicted that miRNAs targeted mRNAs using mRBase targets and Targetscan based on the information from the 59 miRNAs described in Table S2. The intersected mRNAs were chosen from the predicted mRBase and Targetscan mRNAs. Bioinformatics was used to analyze the dysregulated intersection subset of 2026 mRNAs. 59 miRNAs were related to the 164 intersected mRNAs (Table S3).

Tables S2 and S3 were used to construct the ceRNA network. The network was visualized using Cytoscape software 3.0. There were 66 lncRNAs, 33 miRNAs, and 93 mRNAs included in the ceRNA network (Fig. 7). The relationships among the ceRNA network genes are shown in Table S4. The DAVID database (https://david.ncifcrf.gov/mRNAs) was used to analyze genes in the ceRNA network that may be involved in the regulation of signaling pathways. The top 15 KEGG pathways, as determined by analysis of the regulatory signaling pathways, are listed in Table 2. Five cancer-related signaling pathways were enriched and categorized as pathways in cancer, small cell lung cancer, PI3K-Akt signaling pathway, p53 signaling pathway, and microRNAs in cancer. Another 10 non-cancer-related pathways were established and included metabolic pathways, ECM-receptor interaction, and aldosterone synthesis and secretion.

Associations between lncRNA signature and HCC clinical features

The 66 lncRNAs involved in the ceRNA network were selected for further investigation to identify the association between the key lncRNAs and the TCGA database of the clinicopathological features of the 313 patients with HCC. Clinicopathological features included race, gender, tumor grade, TNM stage, lymphatic metastasis, and chronic hepatitis B or C virus infection. 26 lncRNAs (13 upregulated and 13 downregulated) were found to be differentially expressed in HCC patients with different clinical features (P < 0.05). The results suggested that LOC642852, DDX12P, HERC2P2, UCA1, MCM3AP-AS1, CCDC163P, LINC01018, C3P1, AKR7L, AKR1C6P, A1BG-AS1, AQP7P1, UNQ6494, and LINC00261 were associated with tumor grade. CCDC163P, MCM3AP-AS1, HERC2P2, DDX12P, LOC146880, LINC01018, AKR7L, C3P1, A1BG-AS1, AKR1C6P, LINC00261, AQP7P1, and UNQ6494 were associated with TNM stages. LOC642852, CCDC163P, MAFG-AS1, THUMPD3-AS1, SMIM10L2A, AQP7P1, and PRR26 were associated with lymphatic metastasis. UCKL1-AS1, LOC146880, CCDC163P, AKR7L, and LINC01018 were associated with chronic HBV infection. UCA1, GUSBP11, C3P1, LINC00261, and GOLGA2P7 were associated with chronic HCV infection. We also found that some lncRNAs may be associated with gender or race (Table 3).

Prognostic analysis of lncRNA expression and HCC patients’ overall survival

Kaplan–Meier survival analysis was performed based on the RNA sequencing data analysis and clinical features from TCGA in HCC patients. This analysis was conducted to determine the relationships between the 26 key lncRNAs related to the clinicopathological features. The overall expression of these 26 lncRNAs in relation to HCC prognosis in TCGA was calculated using the Cox proportional hazard regression model. Six lncRNAs were associated with HCC overall survival (log-rank P < 0.05), and of these, CCDC163P, LOC146880, LOC642852, MCM3AP-AS1, and THUMPD3-AS1 were negatively correlated with prognosis (P < 0.05); UNQ6494 was positively correlated with prognosis (P < 0.05) (Fig. 8).

qRT-PCR validation

UCKL1-AS1, LOC146880, UCA1, C3P1, LINC00261, and LINC01018 may be important in the progression of HCC and their expressions were detected in 20 patients with newly diagnosed HCC and their paired non-tumor liver tissue samples. qRT-PCR was used to assess the reliability and validity of our bioinformatics analysis results. UCKL1-AS1, LOC146880, and UCA1 were upregulated and C3P1, LINC00261, and LINC01018 were downregulated in HCC tissues (P < 0.05). qRT-PCR validation and bioinformatics analysis gave similar results in 20 newly diagnosed HCC patients (Table S1), suggesting that the bioinformatics analysis used in this study was credible (Fig. 9).