Co-expression analysis and ceRNA network reveal eight novel potential lncRNA biomarkers in hepatocellular carcinoma

Data preprocessing and differential gene selection

The clinical data, RNAseq data and microRNA data of LIHC were downloaded from the The Cancer Genome Atlas (TCGA) database (Akbani et al., 2014) (Table S1). There were 49 pairs of microRNA samples and 50 pairs of RNAseq samples in total. The variation between the RNA and microRNAs data was calculated using EdgeR package (Reimers & Carey, 2006; Varet et al., 2016). Only microRNAs recorded with first 10 and last 10 FC values were selected for subsequent analysis (Li et al., 2018; Shao & Li, 2019). Because these biggest expression changed miRNA may have real major function in HCC. For RNAseq data, the EdgeR package was used to calculate the differential correlationship, and the threshold value for FDR was set at <0.01, — logFC —>1 (Table S2). Clinical data were used to calculate the correlation matrix of clinical information for integrated analysis. We used the “heatmap.2” function in the “gplots” package to create the heatmap.

Determination of ceRNA

Using microRNAs and differential expression genes as input data for target prediction, the RNA22 program was used to predict binding sites of microRNAs based on their sequence characteristics (Loher & Rigoutsos, 2012). Based on these ceRNA interactions, we obtained all the mRNA-miRNA pairs with sharing the number of common miRNAs. Because the number of the ceRNA pairs is obey hypergeometric distribution. We estimated their statistical significance by a hypergeometric test. The potential top 20 different expression miRNA with hypergeometric test in ceRNA network P value less than 0.05 was obtained (Li et al., 2018; Shao & Li, 2019). Specific formulas, such as the differential expression matrix of lncRNA, were used to get the Pearson correlation coefficient (PCC) (Table S3). Based on ceRNA’s mechanism, the expression matrix (EM) for lncRNA to bind between themselves resulted in a PCC of EM >0. Co-expression is one of the features of ceRNA network on account of their interactions. The final ceRNA pair was obtained by intersecting the ceRNA threshold with cutoff p-value <0.05. Cytoscape v3.0 was used to construct the ceRNA network. The overall Kaplan–Meier (KM) survival analysis in each subtype was performed using GEPIA2 database (Tang, 2017).

WGCNA model related computing

 Weighted gene co-expression network analysis (WGCNA), a bioinformatic method used to find correlation patterns among genes, was utilized in this research. WGCNA assumes that gene expression networks is scale-free, it uses a ‘soft’ threshold to determine the weights of the edges connecting genes and merge individual genes to a module. An appropriate soft threshold will make the co-expression network closer to a scale-free network. Then we constructed a signed weighted co-expression network using WGCNA based on the gene expression value across the TCGA samples. We obtained 60 co-expression modules according to the correlations of fpkm value among samples. Each module is represented by an value belongs to the ’eigengene’. This value is identified from the principal component analysis (PCA) of all the gene expression value in the module. Then, we find the relationship between modules and the trait, Eventually, unsupervised WGCNA identified major lncRNAs expression modules with different degrees of correlation to TNM staging using WGCNA R software (Langfelder & Horvath, 2008).

Prediction of interrelationship between lncRNA-related mRNAs

Using RNAseq data from TCGA, the correlation prediction of the lncRNA-lncRNA network and lncRNA-mRNA network was constructed. The cutoff of PCC was 0.7. We retained the intersection with the lncRNA using a cutoff of the top 20% degree. Ultimately, we found eight potentially important lncRNAs with high degrees of connection in the ceRNA network and high correlation at the expression level.

Gene ontology and KEGG enrichment analysis

Gene ontology (GO) annotation analysis was performed using DAVID software (Huang da, Sherman & Lempicki, 2009a). Gene functions for these important indirectly regulated mRNA genes elucidate that this mRNA regulated by lncRNA may have some key biological functions (Huang da, Sherman & Lempicki, 2009a; Huang da, Sherman & Lempicki, 2009b). String database was used for protein-protein analysis (Szklarczyk et al., 2017).

Analysis of differential miRNAs and differential lncRNAs

The differential expression of miRNAs and lncRNAs between normal samples and cancer samples was calculated separately using the EdgeR package (Chen et al., 2017; Law et al., 2016; Maza, 2016). There were 1962 significant differential expression genes in lncRNAs, and 310 differentially expressed miRNA genes found between healthy and cancer-treated samples. Almost half of the mRNAs are upregulated and half downregulated (Fig. 1A). Most of the gene expression for lncRNAs and miRNAs was upregulated in cancer-treated samples (Figs. 1B–1C).

Weighted gene co-expression network analysis of lncRNA

To determine if any of the identified coexpression modules were associated with TNM stage, we calculated the PCC between the MEs and TNM stage. All lncRNAs were merged to 60 modules according the degree of coexpression across the data set in WGCNA. As in the previous study, we assigned each coexpression module an arbitrary color for reference (Di et al., 2019; Zhussupbekova et al., 2016). The hierarchical clustering dendrogram of the eigengenes shows the module size (the number of genes per module) and relationships among these modules (Fig. 2A). Most modules had minimal relationships with each other.

Comparison of module-characteristic eigengenes showed TNM was best correlated with the module MElightpink4 (p = 4E–04) and MElightcyan (P < 0.006), composed of 125 lncRNAs (Fig. 2B). The PCC values ranged from −1 to +1 depending on the power of the relationship. A positive value indicated that the lncRNA within a particular co-expression module increased as the TNM increased, whereas the opposite occurred if the sign of the PCC was negative. We learned that the correlationship between the module and TNM stage with PCC value was accompanied by the corresponding P-value in brackets. These modules included genes that were co-expressed in a particular TNM stage can be used to represent the TNM stage of HCC development (Fig. 2B). These gene may be the risk factors and therapeutic targets in the treatment of HCC.

ceRNA network of lncRNA reveals potential biomarkers in liver cancer

As shown in Fig. 3A, ceRNA network was constructed by high degree lncRNAs. The topological characteristics of ceRNA network were analyzed with degree, betweenness, and closeness (Table S4). Many genes in ceRNA network are with high degree, betweenness, or closeness like AC016773.1, AC145285.2, LINC01569 and DANCR, and so on. This implied ceRNA network may regulate many gene expression through these lncRNAs to have effect on progression of liver cancer. To verify the function of these lncRNAs, Kaplan–Meier survival analysis was perform ed for expression level of these lncRNAs. The results of the survival analysis presented in Figs. 3B and 3C show patients with high expression level of DANCR or AL671710.1 have poor prognosis. It is a validation of our analysis and two of the lncRNAs may affect prognosis.

The interactions between lncRNAs and mRNAs

Combining WGCNA with the results of ceRNA prediction analysis, lncRNA was combined with correlation prediction to find eight important lncRNAs including AL671710.1, TRIM52-AS1, C1orf220, AC022762.2, DANCR, LINC01569, AC084018.1 and MIR194-2HG. Among these lncRNAs, downregulation of TRIM52-AS1 play key role in renal cell carcinoma (Liu et al., 2016). DANCR increases stemness features of hepatocellular carcinoma (Yuan et al., 2016). DANCR is also associated with various cancer (Lu et al., 2018a; Lu et al., 2018b; Mao et al., 2017; Sha et al., 2017; Xu et al., 2018; Yuan et al., 2016). However, there is little known about the function of other six lncRNAs in cancer. There were 124 genes targeted by these eight lncRNAs (Figs. 4A–4H), including EIF3 which functions during the initiation phase of translation. TRIM52-AS1 may influence cancer behavior and function through interactions with regulator EIF3. EIF3 plays a key role in human diseases (Gomes-Duarte et al., 2018; Valasek et al., 2017). Also, there were 30 genes targeted by lncRNA ac084018.1, including m6A reader methyltransferase like 3 (METTL3). As reported in previous research, METTL3 promotes liver cancer progression through YTHDF2 (Balacco & Soller, 2019; Berlivet et al., 2019; Chen et al., 2018; Weng et al., 2018). DANCR and AL671710.1 also have crucial roles through certain key genes (Figs. 4C, 4E).

GO and KEGG pathway enrichment analyses of lncRNA-targeted gene

Additionally, we performed GO and KEGG enrichment analysis of the mRNAs in the network (Figs. 5A–5C). We analysed target genes of the lncRNA based on their enrichment scores for associated GO terms and KEGG pathways using David tools (Huang da, Sherman & Lempicki, 2009a; Huang da, Sherman & Lempicki, 2009b) . The GO and KEGG enrichment analysis concerning the target genes of lncRNAs indicated that the top regulated pathways of lncRNAs were Huntington’s disease, RNA polymerase and Pyrimidine metabolism, and the top regulated functions of lncRNAs were SRP-dependent co-translational protein targeting to membrane, translational initiation, viral transcription, nuclear-transcribed mRNA catabolic process, and nonsense-mediated decay. Further-more, those essential mRNA may interact with each other and function in HCC (Fig. 6A). We can conclude that those lncRNAs affect HCC through the functions and pathways listed above (the flow chart was depicted in Fig. 6B).