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.).

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% CO₂.

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).

A total of eight immunodeficient nude mice were subcutaneously transplanted with 100 µl HepG2 tumor cells (2×10⁶ 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.

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.

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.

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.

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).

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.

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.

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.

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.

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.

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).

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).

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.

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.

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).

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).

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.

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).

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.

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.