Potential therapies for residual hepatoblastoma following incomplete ablation treatment in a nude mouse subcutaneous xenograft model based on lncRNA and mRNA expression profiles.

Tumor recurrence following radiofrequency ablation (RFA) treatment in liver cancer is an important factor affecting patient prognosis. Furthermore, the biological role of long non‑coding RNAs (lncRNAs) in residual hepatoblastoma (HB) tissues after RFA remains largely unknown. By using microarray technology, this study investigated the expression of lncRNAs and mRNAs among four pairs of HB tissues (incomplete ablation treatment and no treatment) in a nude mouse subcutaneous xenograft model. Subsequently, bioinformatics analysis was used to understand the functions and pathways of the identified mRNAs. Finally, a connectivity map (CMap) analysis was conducted to identify potential therapeutic strategies for residual HB tissues. Compared with the untreated nude mouse subcutaneous xenograft model, in the experimental group, a significant difference in the expression of 740 lncRNAs and 663 mRNAs was detected. Subsequently, bioinformatics analysis revealed that the differentially expressed mRNAs were significantly enriched in pathways associated with antigen processing, the presentation of endogenous antigens, the regulation of cellular metabolic processes, MAPK signaling and cell cycle regulation. Additionally, six compounds (valproic acid, metformin, tanespimycin, wortmannin, fulvestrant and MK‑886) were identified by CMap analysis as potential therapeutic agents for the treatment of residual HB tissues. These findings provide a novel insight into the pathogenesis of residual HB and potential therapeutic strategies for aggressive tumor recurrence following RFA treatment in patients with HB.

Introduction

Liver cancer is one of the most common malignancies and is the fourth most common cause of cancer-related death worldwide (1). In 2018, there were estimated to be 840,000 new cases and 780,000 deaths from liver cancer worldwide. Hepatocellular carcinoma (HCC), cholangiocarcinoma and hepatoblastoma (HB) are the main types of liver cancer. Among them, HB is the most common liver tumor in children <4 years old; it has a 5-year survival rate of 60–87.7% (2–4). Generally, liver transplantation, chemotherapy and surgical resection are the preferred treatment strategies for patients with HB (5,6). Although the prognosis of HB has markedly improved in recent decades, the prognosis remains poor for patients with advanced or chemotherapy refractory disease (7).

According to several international guidelines, radiofrequency ablation (RFA) has recently been suggested as a first-line therapeutic strategy for adult patients with early-stage liver cancer (8–11). Studies have reported that patients with liver cancer treated with RFA have a 5-year survival rate of 66–86% (12–14). However, the incidence of recurrence has been reported to be markedly higher in response to RFA compared with hepatic resection, and the recurrence rate may be as high as 60–80% within 5 years of RFA (15–17). Among patients who have undergone hepatic resection, tumor recurrence after RFA is an important factor affecting patient prognosis (18,19). Although several studies have reported the successful RFA treatment of children with liver cancer (4,20), there remains little focus on RFA for the treatment of HB.

Long non-coding RNAs (lncRNAs) are arbitrarily considered to be >200 nucleotides in length and are not translated into proteins (21–23). In the past, lncRNAs were considered transcriptional regulation ‘noise’. However, with the development of techniques for high-throughput sequencing, a large number of newly discovered dysregulated lncRNAs have been reported to serve an important role in tumor growth, invasion and metastasis (24–27). For example, Dong et al reported that upregulated lncRNA-taurin upregulated 1 (TUG1) serves an important role in regulating HB cell function, tumor progression and tumor angiogenesis (28). However, the biological role of lncRNAs in HB tissues remains widely unknown, particularly in residual HB tissues after incomplete RFA treatment.

The connectivity map (CMap) is a wide-ranging drug perturbation database that contains specific gene expression profiles from cultured human cells grown with small bioactive molecules. This database enables the identification of drugs that affect the same molecular pathway as identified differentially expressed genes (29–31). Compounds with a highly positive connectivity score (i.e., ‘score’, 1) may have a highly positive connection with the same molecular pathways associated with the differentially expressed genes. Furthermore, a highly negative connectivity score (i.e., ‘score’, −1) indicates a more negative connection between the compounds and differentially expressed genes.

The aim of this study was to analyze the gene expression profiles of residual HB tissues after incomplete RFA treatment using microarray data combined with the CMap database, in order to improve understanding of the processes underlying rapid proliferation of residual HB tissues and to identify potential therapeutic targets. A flow chart summarizing the present work is presented in Fig. 1.

Figure 1.

Flow chart of the present study. lncRNAs, long non-coding RNAs; RFA, radiofrequency ablation.

Materials and methods

Instruments and reagents

The Cool-tip™ RFA system was purchased from Medtronic. NanoDrop ND-1000 was obtained from NanoDrop Technologies (Thermo Fisher Scientific, Inc.). The Gene Expression Hybridization kit (cat. no. 5188-5242) was purchased from Agilent Technologies, Inc. TRIzol® reagent was obtained from Invitrogen (Thermo Fisher Scientific, Inc.).

Cells

Human HepG2 cells were purchased from Cell Bank of the Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences. Cells were cultured in Dulbecco's modified Eagle's medium (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% fetal bovine serum (Gibco; Thermo Fisher Scientific, Inc.) at 37°C in atmosphere containing 5% CO2.

Animals

Male athymic BALB/c nude mice (age, 6 weeks; weight, 8–22 g) were obtained from Shanghai SLAC Laboratory Animal Co., Ltd. [production permit no. SCXK(Hu)2007-00058]; mice were housed in specific pathogen-free conditions and no infectious diseases were detected in the animals. The animals were housed in class 10,000 clean rooms, at room temperature (20–25°C), with 30–70% humidity and a 12-h artificial light/dark cycle. All experimental protocols were approved by the Animal Care and Experiment Committee of Guangxi Medical University (Nanning, China).

Nude mouse subcutaneous xenograft model

A total of eight immunodeficient nude mice were subcutaneously transplanted with 100 µl HepG2 tumor cells (2×106 cells/ml). After 4 weeks, the tumors had an average diameter of 0.8–1.0 cm and the mice were randomly divided into two groups for further experimentation. One group of nude mice (n=4) underwent incomplete ablation via the ultrasound-guided Cool-tip™ RFA system, after being anesthetized with 1.5% pentobarbital sodium (60 mg/kg body weight; intraperitoneal injection). Referring to previous reports, the parameters of construction of the incomplete ablation model were as follows: Radiofrequency power was set at 30 W and the temperature was 70±5°C; after continuous ablation for 10 sec, the mice were euthanized for further analysis (32,33). Mice in the control group (n=4) did not receive any treatment. After euthanasia, tumors were immediately removed within 10 min, images were captured and the tumor samples were divided equally into two parts (>100 mg). One part was used for pathological experiments, and the other was stored at −80°C or in liquid nitrogen after washing with PBS before being used for array hybridization.

Hematoxylin and eosin (H&E) staining

Tissues were fixed in 4% paraformaldehyde (Sinopharm Chemical Reagent Co., Ltd.) at 20°C for >72 h. Subsequently, the tissues were embedded in paraffin, sliced into 5-µm sections and underwent H&E staining according to a previously published protocol (34). A light microscope (Olympus Corporation) was used to capture the images.

RNA extraction and microarray hybridization

Total RNA was extracted from each sample using TRIzol® reagent, according to the manufacturer's instructions. The quality of RNA was evaluated using a NanoDrop ND-1000. Subsequently, RNA integrity was assessed using standard denaturing agarose gel electrophoresis (0.8–2%). Sample labeling and array hybridization were performed according to the Agilent One-Color Microarray-Based Gene Expression Analysis protocol (Agilent Technologies, Inc.). Firstly, the purified mRNA was obtained by removing ribosomal RNA from total RNA using the mRNA-ONLY™ Eukaryotic mRNA Isolation kit (Epicentre; Illumina, Inc.), according to the manufacturer's protocol. Secondly, the Arraystar Flash RNA Labeling kit (Arraystar, Inc.) was used to amplify the mRNA and transcribe it into fluorescent cRNA, according to the manufacturer's protocol. Subsequently, labeled cRNAs were purified using the RNeasy Mini kit (Qiagen, Inc.), and their consistency and specific activity were calculated using a NanoDrop ND-1000, according to the manufacturer's protocol. The instructions of the Agilent Gene Expression Hybridization kit (cat. no. 5188-5242; Agilent Technologies, Inc.) were followed to perform the chip hybridization. Finally, the Agilent G2505C DNA Microarray Scanner (Agilent Technologies, Inc.) was used to wash, fix and scan the hybridized arrays.

Microarray data analysis

Agilent Feature Extraction software (version 11.0.1.1; Agilent Technologies, Inc.) was applied to extract data for further analysis. Subsequently, the raw data were normalized with log2 transformations (GeneSpring GX v12.1; Agilent Technologies, Inc.). Fold changes >2 and a P<0.05 were considered statistically significant differences in the expression of lncRNA and mRNA.

Co-expression network construction

To further understand the possible role of lncRNAs and mRNAs in the incompletely ablated HB tissue, a lncRNA-mRNA co-expression network was constructed according to the correlation coefficient of the differential expression data. The Pearson conjunction coefficient was used to evaluate the correlation of the differentially expressed genes. Genes with an absolute value of correlation coefficient ≥0.99 and P<0.0001 were selected, and the lncRNA-mRNA co-expression network was drawn using Cytoscape 3.7.2 (https://cytoscape.org/) (35).

Standard gene set enrichment analysis (GSEA)

GSEA 3.0 software (http://www.broad.mit.edu/gsea) was used to further analyze the biological signatures of the gene set in the incompletely ablated HB tissue (36,37). The sets, including hallmark gene sets, chemical and genetic perturbations, canonical pathways, cancer modules, Gene Ontology (GO) gene sets, oncogenic signatures and immunologic signatures were collected from the Molecular Signatures Database. The normalized enrichment score (NES) was used to quantify the magnitude of enrichment. Using 1,000 permutation test runs, the cut-off criteria were false discovery rate q-value ≤0.01 and nominal P≤0.001.

Functional annotation analysis and pathway enrichment analysis of the differentially expressed genes

In addition, bioinformatics analysis, including Kyoto Encyclopedia of Genes and Genomes (KEGG; www.kegg.jp) pathway analysis and GO analysis, was performed to explore the roles of the mRNAs differentially expressed in residual tissue after RFA treatment compared with untreated HB tissue. The GO project (http://www.geneontology.org) provides the properties of genes and/or gene products through three domains: Biological process (BP), cellular component (CC) and molecular function (MF). R package version 2.8.0 (‘GOstats.’ and ‘GeneAnswers’ packages; http://www.r-project.org/) is a free script that synthesizes GO terms and comprehensively visualizes the results. KEGG is an online database of genomes, enzymatic pathways, biochemical reactions, diseases and drugs. The pathway database was used to perform pathway analysis of differentially expressed mRNAs to infer their molecular functions.

Protein-protein interaction (PPI) networks and hub genes

The Search Tool for the Retrieval of Interacting Genes database (http://string-db.org/) was used to search for affiliations among proteins encoded by the genes that were differentially expressed in residual HB tissues after incomplete RFA treatment. The database aims to integrate large amounts of data on protein-protein associations (38). The minimum required interaction score for a connection node to be considered was 0.9. Any gene encoding proteins with >10 connection nodes was regarded as a hub gene.

CMap analysis

To identify drugs that may affect the same molecular pathways that were linked to the differentially expressed genes between the residual HB tissues after incomplete RFA treatment and untreated HB tissues, the observed differentially expressed genes and hub genes were compared with the gene expression levels in the CMap database (http://portals.broadinstitute.org/cmap/). Prior to conducting the CMap search, the probe IDs of the genes were converted into the Human Genome HT U133A Array format using the Batch Query function on the Affymetrix GeneChip website (https://www.affymetrix.com/site/mainPage.affx), and these data were then inputted into the CMap database. A compound with a connectivity score ≤-0.75, which indicated a highly negative association with the expressed genes in the database, was considered to be a potential drug for RFA-treated liver cancer. Only the compounds present in both comparison results of DEGs and hub genes were selected for further research.

Prediction of interactions between chemicals and proteins

The Search Tool for the Interacting Chemicals database (STITCH; http://stitch.embl.de) was used to identify interactions between proteins and the aforementioned compounds, which may help to identify a novel therapy for HB following incomplete RFA (39). In the fifth version of STITCH, the database includes >9,600,000 proteins and 430,000 compounds. The information and SMILES networks of compounds required for the STITCH interaction analysis were collected from the PubChem Compounds website (https://pubchem.ncbi.nlm.nih.gov/) for further analysis. Compounds with clinical application were selected as candidate drugs for further research.

Molecular docking

Docking was carried out using SYBYL-X 2.0 software (Certara, Princeton, NJ, USA) to discover novel ligands between the aforementioned candidate drugs and proteins. The crystal structures of proteins were obtained from the Research Collaboratory for Structural Bioinformatics Protein Data Bank (RCSB PDB) (http://www.rcsb.org/) (40). A consensus score (CScore) ≥4 was used to represent high binding affinities between the proteins and novel drugs (41,42).

Results

Pathological examination of the tumor tissue after RFA

In the untreated group, the tumor was red to the naked eye (Fig. 2A), whereas in the RFA-treated group, the tumor exhibited a white burning appearance (Fig. 2C). Subsequently, the two groups underwent pathological examination using H&E staining (Fig. 2B and D). A small number of fragmented dark blue nuclear structures were detected at the junction of the necrotic zone and the marginal zone in the RFA-treated group (magnification, ×200). Light pink/dark red staining indicates vascular tissue (magnification, ×200).

Figure 2.

Pathological examination of the tumor tissue from the RFA and untreated groups. (A) Gross view of tumor tissue in the untreated group. (C) Gross view of tumor tissue in the RFA-treated group. The black box represents the ablated area, and the blue arrow points to the puncture site of the ablation needle. (B) HE staining of the untreated group, observed under a microscope (×200 magnification). (D) HE staining of the RFA-treated group, observed under a microscope (×200 magnification). HE, hematoxylin and eosin; RFA, radiofrequency ablation.

Expression profile of lncRNAs and mRNAs in incomplete RFA residual tissues of a subcutaneous HepG2 cell-based tumor transplantation nude mouse model

To identify lncRNAs and mRNAs, the expression of which was modified by RFA treatment of HB, a subcutaneous tumor transplantation model of nude mice was generated using HepG2 cells; one group was treated with incomplete RFA and the other received no treatment.

Fold-changes (incomplete RFA-treated nude mice vs. untreated nude mice subcutaneously transplanted with HepG2 cells) and P-values were calculated after normalizing the in vivo expression data (Fig. 3A and D). Using the GeneSpring GX v12.1 software package, the median expression levels of the lncRNAs and mRNAs in the experimental and control groups were ~5.2 and 6.2, respectively (Fig. 3C and F). Overall, 740 lncRNAs and 663 mRNAs were significantly differentially expressed in the four incomplete RFA-treated samples compared with the untreated samples (fold change ≥2.0, P<0.05; Fig. 3B and E). Among these, 161 lncRNAs and 87 mRNAs were significantly enriched in the experimental group, whereas 579 lncRNAs and 576 mRNAs had significantly lower expression in the experimental group. The expression of lncRNA RP11-150O12.5 (fold change, 19.60; P=0.044) was significantly upregulated, but the expression of AK094929 (fold change, 26.38; P=0.008) was significantly downregulated. These differentially expressed lncRNAs may act as the chief regulators in the recurrence and progression of HB after incomplete RFA treatment.

Figure 3.

Expression profile of lncRNAs and mRNAs in incomplete radiofrequency ablation residual tissues obtained from a subcutaneous tumor transplantation model. Heatmaps of (A) lncRNAs and (D) mRNAs. Volcano plots of differentially expressed (B) lncRNAs and (E) mRNAs. Box plots of (C) lncRNAs and (F) mRNAs after quantile normalization. lncRNA, long non-coding RNA.

Construction of the lncRNA-mRNA co-expression network

A co-expression network was constructed to more comprehensively understand the relationship between the lncRNA and mRNA expression levels in the incompletely ablated HB tissue (Fig. 4). Overall, 1,201 pairs of co-expression correlations between lncRNAs and mRNAs were collected; these pairs included 387 mRNAs (43 upregulated and 344 downregulated) and 468 lncRNAs (78 upregulated and 390 downregulated). Within this co-expression network, there were 975 positive relationships and 226 negative relationships. AV30S1-TTI1 and IGKV-HOXA13 were the most positive relationships. XLOC_006929-PTH1R was the most negative relationship.

Figure 4.

Co-expression network of differentially expressed lncRNAs and mRNAs in incomplete radiofrequency ablation residual tissues obtained from a subcutaneous tumor transplantation model. Blue circles represent differentially expressed lncRNAs; red triangles represent differentially expressed mRNAs; red lines indicate positive correlation between lncRNA and mRNA expression; blue lines indicate negative correlation.

GSEA

GSEA is an approach used to identify chemical and genetic perturbations, pathways, cancer modules, GO terms, oncogenic signatures and immunologic signatures related to all expressed genes. In all of the assessed gene sets, there were eight gene sets that exhibited significant enrichment in the incompletely ablated HB tissue compared with the control tissue, including TUMORS PCA2 [Normalized Enrichment Score (NES)=2.13], IL6 DEPRIVATION (NES=2.14), EARLY T LYMPHOCYTE (NES=2.15), ENDOCHONDRAL BONE MORPHOGENESIS (NES=2.09), TRAF6KO EFF CD8 TCELL (NES=2.09), NEOCORTEX BASAL RADIAL GLIA (NES=2.10), IMMATURE B LYMPHOCYTE (NES=2.07) and MYC TARGETS (NES=2.06) (Fig. 5).

Figure 5.

Gene set enrichment analysis of the differentially expressed genes in the incomplete radiofrequency ablation residual tissues. NES, normalized enrichment score.

Functional annotation of the differentially expressed genes

GO and KEGG analyses were conducted to predict the functions and interacting pathways of the differentially expressed genes in the incomplete ablation HB tissue. With regards to GO BP terms, the upregulated genes were most enriched in ‘antigen processing and presentation of endogenous antigen’ (GO:0019883), ‘embryonic morphogenesis’ (GO:0048598), ‘antigen processing and presentation of exogenous peptide antigen via MHC class I, TAP-dependent’ (GO:0002479), ‘protein localization to chromatin’ (GO:0071168) and ‘antigen processing and presentation of exogenous peptide antigen via MHC class I’ (GO:0042590) (Fig. 6A). In addition, the downregulated genes were particularly enriched in ‘regulation of lipid metabolic process’ (GO:0019216), ‘phosphorus metabolic process’ (GO:0006793), ‘phosphate-containing compound metabolic process’ (GO:0006796), ‘regulation of cellular metabolic process’ (GO:0031323) and ‘nucleobase-containing compound metabolic process’ (GO:0006139) (Fig. 6B). For MF, the upregulated genes were most enriched in ‘voltage-gated cation channel activity’ (GO:0022843), ‘voltage-gated ion channel activity’ (GO:0005244), ‘voltage-gated channel activity’ (GO:0022832), ‘cation channel activity’ (GO:0005261) and ‘metal ion transmembrane transporter activity’ (GO:0046873) (Fig. 6C). The downregulated genes were most enriched in ‘SNARE binding’ (GO:0000149), ‘syntaxin binding’ (GO:0019905), ‘tRNA binding’ (GO:0000049), ‘protein methyltransferase activity’ (GO:0008276) and ‘p53 binding’ (GO:0002039) (Fig. 6D). For CC, the upregulated genes were most enriched in ‘9+0 non-motile cilium’ (GO:0097731), ‘photoreceptor cell cilium’ (GO:0097733), ‘endosome lumen’ (GO:0031904), ‘transmembrane transporter complex’ (GO:1902495) and ‘nuclear heterochromatin’ (GO:0005720) (Fig. 6E). The downregulated genes were most enriched in ‘intracellular’ (GO:0005622), ‘intracellular membrane-bounded organelle’ (GO:0043231), ‘organelle inner membrane’ (GO:0019866), ‘nucleus’ (GO:0005634) and ‘obsolete intracellular part’ (GO:0044424) (Fig. 6F). Furthermore, the upregulated genes were enriched in the following pathways: ‘MAPK signaling pathway’, ‘acute myeloid leukemia’, ‘arrhythmogenic right ventricular cardiomyopathy (ARVC)’, ‘thyroid hormone synthesis’ and ‘pancreatic cancer’ (Fig. 7A), and the downregulated genes were enriched in ‘human papillomavirus infection’, ‘regulation of lipolysis in adipocytes’, ‘glycerolipid metabolism’, ‘ECM-receptor interaction’ and ‘cell cycle’ pathways (Fig. 7B).

Figure 6.

GO enrichment analysis of the upregulated mRNAs: (A) Biological process; (C) molecular function; (E) cellular component, and downregulated mRNAs: (B) Biological process; (D) molecular function; (F) cellular component. GO, Gene Ontology.

Figure 7.

KEGG pathway analysis of the (A) upregulated and (B) downregulated mRNAs. KEGG, Kyoto Encyclopedia of Genes and Genomes.

PPI network and hub genes for the differentially expressed genes

Hiding the disconnected nodes, a PPI network comprising 247 nodes and 462 edges was generated to view the interactions among the 663 differentially expressed genes (Fig. 8). Finally, 11 hub genes, including two upregulated (VAMP2 and HELZ2) and nine downregulated genes (HIST1H2BA, CREBBP, HSPA8, GNGT2, NCOA1, DVL2, NUP85, CD36 and PPP2R5A) were identified from the PPI network.

Figure 8.

Protein-protein interaction network of the differentially expressed genes, as constructed by Cytoscape. The size and brightness of the node indicates the degree. Line thickness represents combined connection strength. Green and blue arrows represent the upregulated and downregulated hub genes, respectively.

Repurposed drugs for the treatment of RFA-treated HB

After removing duplicates, 41 compounds with connectivity scores ≤-0.75, which indicated that the compounds may negatively affect the differentially expressed genes, were identified from the CMap database. In addition, 146 compounds with connectivity scores ≤-0.75 were selected from the database using the hub genes imputed. Finally, only six drugs that appeared at the intersection of the aforementioned results were selected as candidate drugs: Valproic acid, metformin, tanespimycin, wortmannin, fulvestrant and MK-886 (Fig. 9).

Figure 9.

Venn diagram of the number of overlapping predicted drugs using all DEGs and hub genes. DEGs, differentially expressed genes.

Construction of a drug-target network for the six candidate drugs

To identify the interactions between proteins and chemicals, a drug-target network was constructed for the six candidate drugs using the STITCH database (Fig. 10). The results identified four drug-target networks containing a total of 37 genes, which may be relevant for HB therapy. Additionally, two drug-target networks were not available: Tanespimycin and wortmannin. These compounds and target proteins may provide potential treatment direction in the future.

Figure 10.

Drug-target interaction networks. (A) Valproic acid, (B) metformin, (C) fulvestrant and (D) MK-886.

Molecular docking study

Two pairs of compounds and target proteins with the highest combined score [solute carrier family 2 member 4 (SLC2A4) and metformin, combined score 0.99; ATP-dependent translocase ABCB1 (ABCB1) and valproic acid, combined score 0.95] were selected to verify the binding affinities using molecular docking analysis. The CScore of the affinities between SLC2A4 (PDB code: 5VNM) and metformin was 4 (Fig. 11A), and the CScore of the affinities between ABCB1 (PDB code: 6QEX) and valproic acid was 5 (Fig. 11B). These findings may provide new insights about the selection of future drugs.

Figure 11.

Molecular docking of predicted compounds and target proteins. (A) Metformin and solute carrier family 2 member 4; (B) valproic acid and ATP-dependent translocase ABCB1.

Discussion

HB is the predominant type of pediatric liver tumor, and has two primary pathological types: Epithelial and mixed (2). Epithelial cells are the major histological cell type in HB and this type of cancer usually has a good prognosis (43). However, the small-cell undifferentiated subtype of HB has a poor prognosis (2,44,45). RFA has been widely used for treating liver cancer, particularly patients with early detection of cancer. However, because of residual tumor cells, microscopic metastasis and vascular invasion, this type of therapy is often unable to offer complete tumor clearance (46). Notably, after sub-lethal treatment of cancer, surviving tumor cells have a significantly enhanced ability to invade and metastasize (19,47,48). Increasing evidence has shown that lncRNAs may be involved in tumorigenesis and disease progression in HB. Dong et al (28,49) reported that knockdown of lncRNA-CRNDE and lncRNA-TUG1 reduces growth and angiogenesis in HB cells, which may be a prospective target in the diagnosis and therapy of HB. In addition, Zhang et al reported that lncRNA OIP5-AS1, as a competing endogenous RNA, may inhibit cell proliferation, metastasis and epithelial to mesenchymal transition in HB (50). However, few studies have described the changes in lncRNA expression in incompletely ablated HB tissue.

In the present study, HepG2 tumor cells were subcutaneously transplanted into eight male nude mice. Subsequently, four nude mice were randomly selected as the experimental group and subjected to incomplete ablation treatment, which served as an in vivo model for residual HB cells following insufficient RFA. After microarray analysis, it was determined that a total of 740 lncRNAs and 663 mRNAs were significantly differentially expressed compared with in the untreated nude mouse subcutaneous xenograft model. A co-expression network was then constructed to predict the target genes and factors regulated by the lncRNAs. Simultaneously, GSEA was performed to further analyze the chemical and genetic perturbations, pathways, cancer modules, GO terms, oncogenic signatures and immunologic signatures related to all expressed genes. GO term analysis, KEGG pathway analysis and a PPI network analysis were then performed to understand the function and pathways of the differentially expressed genes. The results revealed that the genes were enriched in T lymphocytes, CD8 effector T cells, B cell lymphoma tumors, etc. The results from GO term and KEGG pathway analyses indicated that the differentially expressed genes were involved in functions including antigen processing, as well as the presentation of endogenous antigens, regulation of cellular metabolic processes, MAPK signaling and cell cycle regulation. These results indicated that the differentially expressed mRNAs may have an important role in the invasion and metastasis of incompletely ablated HB cells through the regulation of immunologic activity, energy changes and cellular metabolism.

While there are numerous approaches for treating patients with HB, complete surgery is still highly recommended (2). However, the use of RFA therapy in HB is a good alternative (4,51–53). RFA therapy combined with surgery and/or chemotherapy may be a promising and effective way to treat children with HB (52,53). However, research has shown that after RFA in patients with early-stage liver cancer (Barcelona Clinic Liver Cancer stage 0 or A), the frequency of aggressive intrasegmental recurrence is markedly increased to 15% and this is an important factor affecting patient prognosis (47). In this study, six compounds were identified as candidate drugs for the treatment of patients with aggressive intrasegmental recurrence: Valproic acid, metformin, tanespimycin, wortmannin, fulvestrant and MK-886.

Valproic acid is a fatty acid with anticonvulsant and anti-manic properties that has marked therapeutic value for the treatment of bipolar disorder and epilepsy, particularly generalized epilepsies (54). A recent discovery reported a novel ability of valproic acid, in that it can inhibit histone deacetylase and may serve an important role in cancer treatment (54–56). Zhu et al reported that, combined with sorafenib, valproic acid may inhibit growth and induce apoptosis of HCC cells (57). However, valproic acid is also known to cause liver injury. Metformin is a biguanide hypoglycaemic agent that is a first-line therapeutic strategy for the treatment of type 2 diabetes (58). In 2001, Schneider et al first reported the activity of metformin as an antitumor agent in hamsters (59). To date, the value of metformin in the treatment of cancer has been further explored in areas such as liver cancer, ovarian cancer, breast cancer, prostate cancer and colorectal cancer (60–64). Zhang et al (65) reported that metformin and curcumin in combination effectively inhibit the growth, metastasis and angiogenesis of HCC. Tsai et al (66) reported that therapeutic programs combined with metformin and rapamycin could enhance autophagic cell death in HCC. Tanespimycin is a benzoquinone antineoplastic antibiotic that has been widely investigated for the treatment of leukaemia and solid tumors, such as breast cancer (67), gastric cancer (68) and pancreatic cancer (69). Wortmannin is an androstadiene metabolite isolated from Penicillium wortmannii. PX-866 is an analogue of wortmannin and has stable antitumor activity (70). Fulvestrant is a pure estrogen receptor antagonist and has been widely used in the treatment of advanced breast cancer (71,72). MK-886 is a member of the class of indoles, and evidence has suggested that it can cause cell death in several tumor cell types, such as gastric cancer cells, prostate cancer cells and lung tumor cells (73–75).

In conclusion, this study generated a model to imitate residual HepG2 cells following incomplete RFA treatment using a nude mouse subcutaneous xenograft model. Using microarray data analysis, profiles of lncRNAs and mRNAs were obtained in incompletely ablated HB cells. By analyzing those data in relation to data from the CMap database, six compounds were identified as potential drugs for the treatment of aggressive intrasegmental recurrence in patients with HB.