The data of liver cancer patients and RNA-seq expression results were downloaded with RTCGA Toolbox package in R (version 3.5.3) (8). The gene microarray with cancer tissue data (GSE54236) (9,10) was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) (11).

SPSS software 23.0 (IBM Corporation, Armonk, NY, USA) was used for data analysis. Boxplots were used for discrete variables to measure the expression differences, and Chi-square tests were used to examine the correlation between ZNF385D-AS2 expression and clinical data (12). Receiver-operating characteristic curve (ROC) was drawn by ‘p-ROC package’ to evaluate the capability of diagnosis. We divided patients into high and low ZNF385D-AS2 expression groups using the optimal cutoff value of overall survival (OS) as determined by the Youden index. Kaplan-Meier curves were used to compare the differences in the OS and relapse-free survival (RFS) by using survival package in R (13). Univariate Cox analysis was used to select the related variables. Then, multivariate Cox analysis was applied for the influence of ZNF385D-AS2 expression on OS and RFS of the patients (14).

GSEA is a computational method that determines whether an a priori defined set of genes shows statistically significant concordant differences between two biological states. In this study, GSEA was performed by using GSEA software 3.0 from the Broad Institute (UC San Diego, San Diego, CA, USA) (http://software.broadinstitute.org/gsea/index.jsp) (15). The gene expression data were RNA-seq data from TCGA-LIHC and GEO database. The gene set of ‘c2. cp.biocarta.v6.2.symbols.gmt’, ‘c3. cp.biocarta.v6.2.symbols.gmt’, ‘c5.cp.biocarta.v6.2.symbols. gmt’ and ‘h.all.v6.2.symbols.gmt’, which summarizes and represents specific, well-defined biological states or processes, was downloaded from the Molecular Signatures Database (http://software.broadinstitute.org/gsea/msigdb/index.jsp) (15). The normalized enrichment score (NES) was acquired by analyzing with permutations for 1,000 times. A gene set was considered to be significantly enriched at a normal P-value <0.05 and false discovery rate (FDR) <0.25.

The Database for Annotation, Visualization, and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/) and KO-Based Annotation System (KOBAS) (http://kobas.cbi.pku.edu.cn/) were used to conduct relevant pathway analysis (16,17), and Gene Ontology (GO) analysis was performed for the functional annotation of the predicted genes (17). Three GO terms [biological process (BP), cellular component (CC) and molecular function (MF)] were utilized to identify the enrichment of target genes. In addition, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed for the functional annotation of these genes. GO terms and KEGG pathway with P-values <0.05 were considered statistically significant. The enrichment map of annotation analysis was drawn using Cytoscape (version 3.3.1) (http://www.cytoscape.org/cy3.html) and R (version 3.5.3).

In order to analyze the specific functions of lncRNA ZNF385D-AS2, it was necessary to analyze the interaction between lncRNAs and the coding genes as well as the interacting miRNAs (18). Based on data from the TCGA-LIHC and GEO (GSE54236) databases, using the Pearson and Spearman correlation analysis based on the Logistic function (‘WGCNA’ package in R), the co-expression relationship between ZNF385D-AS2 and the coding genes and miRNAs were identified, and the co-expression network between ZNF385D-AS2 and the coding genes and miRNAs were mapped in Cytoscape. At the same time, comparative analysis of several lncRNA databases, such as NONCODE (http://www.noncode.org/), LncRBase (http://starbase.sysu.edu.cn/), Co-LncRNA (http://www.biobigdata.com/Co-LncRNA/), was carried out to identify the miRNAs that exhibited regulatory relationships with ZNF385D-AS2 and mRNAs which these miRNAs may regulate. After comparing the data from the TCGA and GEO databases, the possible ceRNA patterns in HCC were predicted and plotted which were also mapped in Cytoscape. Circos plots was generated based on the remappings from these two predicted results. Circos plots were generated using the Circos visualization tool in R (‘RCircos’ package in R) (19).

Both gene expression and clinical data of patients with liver cancer were downloaded from The Cancer Genome Atlas (TCGA-LIHC) database (9,10). The total number of patients was 427. After initial screening, we omitted 13 normal samples and 74 tumor samples with missing or unclear information, and the remaining 303 tumor samples and 37 normal samples were available. The detailed clinical characteristics, including age, sex, TNM stage, survival status, pathological status, and ethnic compositions are shown in Table I.

Using boxplots, we measured the differences in ZNF385D-AS2 expression in liver cancer patients and control subjects. As shown in Fig. 1A, we evaluated the overall expression trend of ZNF385D-AS2 in liver cancer, then found that ZNF385D-AS2 expression was significantly lower in primary HCC tissues than that in normal liver tissues (P<0.001). Moreover, as shown in Fig. 1A-K, there was differential ZNF385D-AS2 expression in the groups according to sex (P=0.046; Fig 1C), vital status (P=0.037; Fig. 1D), clinical stage (P=0.001; Fig. 1E), T stage (P=0.001; Fig. 1F) and survival time (P=0.033; Fig. 1K). Of note, differences in ZNF385D-AS2 expression were observed according to patient age as well as the TNM stage of cancer, and clinicopathological parameters (Fig. 1). We also collected data on the expression of ZNF385D-AS2 in several common digestive tumors based on TCGA database. After horizontal comparison, we found that except for cholangiocarcinoma and rectal cancer, the expression of ZNF385D-AS2 was reduced in most digestive system tumors (Fig. 2).

According to Chi-square tests, the correlation between the clinical features and the expression of ZNF385D-AS2 was analyzed and is documented in Table II. The expression of ZNF385D-AS2 was highly associated with sex (χ²=3.846, P=0.050), T stage (χ²=3.875, P=0.049), M stage (χ²=4.221, P=0.040), N stage (Fisher's exact text, P<0.001) and clinical stage (χ²=4.365, P=0.037).

We generated Kaplan-Meier curves of overall survival (OS), and log-rank tests showed that ZNF385D-AS2 low-expression was associated with poor OS (P=0.0079; Fig. 3A). Further subgroup analysis (Fig. 3A-O) showed that ZNF385D-AS2 low-expression was associated with the poor OS of patients with sex (male) (P=0.015; Fig. 3D), M0 stage (P=0.0099; Fig. 3H), N0 stage (P=0.028; Fig. 3I), age ≤55 (P=0.023; Fig. 3J), early pathological stage (G1/2) (P=0.044; Fig. 3L) and race (Asian) (P=0.0046; Fig. 3N).

As shown in Fig. 4, the ROC of ZNF385D-AS2 was executed, and the area under the curve (AUC) was 0.594, which represented moderate diagnostic ability.

In ZNF385D-AS2 low-expression patients, we used univariate analysis and selected the critical variables including age, sex clinical stage, pathological grade and TMN classification. Multivariate analysis with the Cox proportional hazards model indicated that clinical stage (HR=1.418, P=0.011) and T classification (HR=1.713, P<0.001) were independent prognostic factors for patients with HCC (Table III).

To identify the biological functions activated in liver cancer, we conducted GSEA between high and low ZNF385D-AS2 expression data sets. GSEA revealed significant differences (FDR <0.25, P<0.05) in the enrichment of ‘MSigDB Collection’, and the specific contents are shown in Table IV. In the nucleus, ‘histone deubiquitination’, ‘V-D-J recombination’ and ‘translation factor activity RNA binding’ were found to be differentially enriched in the ZNF385D-AS2 low-expression phenotype. Meanwhile, Table IV shows that ‘cytoplasmic exosome RNase complex’, ‘nuclear exosome RNase complex’, ‘retrograde transport endosome to Golgi’ and ‘Golgi organization’ exhibited a positive correlation with ZNF385D-AS2. In the cytoplasm, low expression of ZNF385D-AS2 was also related to various cytological behaviors such as ‘protein transport along microtubule’, ‘negative regulation of defense response to virus’ and ‘oxidoreductase activity acting on the CNH group of donors NAD⁺ or NADP⁺ as acceptor’. Moreover, Table IV shows that as the expression of ZNF385D-AS2 is decreased, various intracellular signaling pathways are also affected. These signaling pathways encompass the metabolic synthesis and degradation of substances, such as ‘one carbon pool by folate’, ‘glyoxylate and dicarboxylate metabolism’, ‘fatty acid metabolism’, ‘selenoamino acid metabolism’, ‘fructose and mannose metabolism’ as well as ‘aminoacyl tRNA biosynthesis’, ‘glycosyl-phosphatidylinositol GPI anchor biosynthesis’, ‘RNA degradation’ and ‘lysine degradation’. Moreover, it also affects signal transfer and material transport within HCC cells, including ‘ABC transporters’, ‘Snare interactions in vesicular transport’ and ‘basal transcription factors’, in addition to the progression of various diseases, including ‘acute myeloid leukemia’, ‘hepatitis’ and ‘endometrial cancer’.

During the progression of liver cancer, the expression of ZNF385D-AS2 was found to gradually decrease. In this process, expression levels of a certain number of mRNAs were found to be altered as ZNF385D-AS2 decreased, and these results are summarized in Table V. The expression levels of MEK, EGFR, ERB2, mTOR initially increased and then decreased. A positive regulation relationship was found between the expression levels of ZNF385D-AS2 and GLI1, CAMP, cyclin D1, Wnt. However, there still exists a negative regulatory relationship between the expression level of JAK2, PDGF-ERK and BRCA1. It is worth noting that P53, a gene that inhibits the development of cancer, was found to first decrease and then increase as the amount of ZNF385D-AS2 expression decreases.

We next used R's ‘edgr’ package to calculate the difference (log fold-change >1, P<0.05) in expression between mRNAs, miRNAs and lncRNAs in the TCGA-LIHC and GEO (GSE54236) databases. Then we conducted comparative analysis of several lncRNA databases, such as NONCODE, LncRBase, Co-LncRNA to identify miRNAs that have regulatory relationships with ZNF385D-AS2 and mRNAs that these miRNAs may regulate. After comparing the data from the TCGA and GEO databases, the possible ceRNA patterns in HCC were predicted and plotted in Cytoscape (Fig. 5A). At the same time, we used Pearson and Spearman correlation analysis based on the Logistic function (‘WGCNA’ package in R), of the co-expression relationship between ZNF385D-AS2. The mRNAs and miRNAs were identified, and the co-expression network between ZNF385D-AS2 and the mRNAs and miRNAs was also mapped (Fig. 5B). The result strongly suggested that there is a regulatory signaling axis that exists between ZNF385D-AS2 and miR-96 and miR-182.

We then integrated the regulatory networks obtained by these two different methods, and again performed GO terms and KEGG pathway analyses for these mRNAs after removal of the duplicates. Cytoscape was used to conduct an analysis map of genes enriched in the GO terms and to construct an interaction network for the related genes (Fig. 6A). As shown in Fig. 6B and C, these specific mRNAs were most highly enriched in the following GO terms: Molecular function (MF) (growth factor activity, transforming growth factor β receptor binding), biological process (BP) (axon guidance, cell development, regulation of MAPK cascade, positive regulation of cell division, positive regulation of pathway-restricted SMAD protein phosphorylation, cell growth, SMAD protein signal transduction) and cellular component (CC) (extracellular space). Based on the results in Fig. 5, we used R to perform a dot-enrichment of the KEGG pathway obtained from KOBAS. These signaling pathways are able to affect material metabolism within HCC cells, signal transduction, progression of various diseases, and other biological behaviors (Fig. 7).

Alignments of gene co-expression maps to GRCh38.95 identified reference genomic regions contributing to the composition of 59 genome sets. According to the co-expressing genes and ceRNA regulatory network predictive analysis results, we comprehensively analyzed the gene co-expression result and ceRNA regulatory network in GEO (GSE54236) and TCGA-LICH with ZNF385D-AS2. In order to find the most represented reference fragments, all GRCh38.95 loci present in the gene co-expressed maps were deduced and merged, resulting in 59 reference donor fragments, which settled in the outermost track. The numbers of gene co-expressed maps containing each of these fragment labels were then marked in the second track, excluding duplicate counts. We found 96 pairs of genes which existed in a co-expressed phenomenon. In the inner sector, we linked these pairs of gene sets that had co-expression relationships with lines. The sum of gene co-expressed map alignments across the whole genome was used as the links for the 59 gene co-expressed maps in the Circos plot (Fig. 8).