The long noncoding RNA expression profile of HCC
The research flow diagram of this study was shown in Fig. 1. To identify lncRNAs that were differentially expressed in HCC, 6 GEO datasets were integrated into analysis (39 normal samples vs 44 HCC samples), including GSE58043, GSE67260, GSE112613, GSE89186, GSE64631 and GSE70880. We obtained 66 differentially expressed lncRNAs, among which 37 were up-regulated and 29 were down-regulated (Fig. 2a). Then we downloaded the relevant expression profiles of the TCGA LIHC dataset (50 normal samples vs 374 HCC samples) and obtained 2685 differentially expressed lncRNAs, among which 2323 were up-regulated and 362 were down-regulated (Fig. 2b). There were 26 lncRNAs that both differentially expressed in GEO joint dataset and TCGA LIHC dataset, among which 15 were up-regulated and 11 were down-regulated (Fig. 2c). An additional file shows these in more detail [see Additional file 3]. The expression heatmaps of GEO joint dataset and TCGA LIHC were displayed in Fig. 2d and Fig. 2e, respectively.
The potential interactive miRNAs of 26 differentially expressed lncRNAs were screened by R software based on the highly conserved microRNA families file downloaded from the miRcode V11 database. R software predicted the target genes of corresponding miRNAs that were simultaneously supported in miRDB, miRTarBase and TargetScan. Only 7 lncRNAs (LINC00221, FAM99B, LINC00355, MAGI2-AS3, CRNDE, PWRN1, FBXL19-AS1) were predicted to have highly conserved targeted miRNAs [see Additional file 4]. We used GEPIA2 to carry out survival analyses for these 7 lncRNAs, and found that only FBXL19-AS1 (Fig. 3a), FAM99B (Fig. 3b) and CRNDE (Fig. 3c) were associated with the prognosis of HCC patients. Among the 3 lncRNAs, the function of FBXL19-AS1 in HCC has not yet been studied and FBXL19-AS1 might serve as a novel HCC biomarker. Therefore, FBXL19-AS1 was selected for further study.
FBXL19-AS1 was significantly up-regulated in both GEO joint dataset and TCGA LIHC datasets and mainly enriched on cell cycle, cancer and inflammation-related pathways by GSEA (Fig. 3d). We speculated that FBXL19-AS1 might play an important role in the occurrence, development and prognosis of HCC. In addition, we found that FBXL19-AS1 was located in cytosol according to the LncLocator database [see Additional file 5], which was the basis of establishing a more reliable ceRNA network. Interestingly, FBXL19, the complementary gene of FBXL19-AS1, was found to be positively correlated with the tumor infiltrating immune cells (TIICs) in HCC, including B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils and dendritic cells [see Additional file 6].
Expression of FBXL19-AS1 in HCC tissues and its prognostic value
In order to verify the expression status of FBXL19-AS1 in HCC patients and to study its clinical significance, we collected 57 pairs of fresh tissues including HCC and adjacent non-tumor tissues. The qPCR results showed that the expression of FBXL19-AS1 in HCC tissue was significantly higher than adjacent non-tumor tissue (Fig. 4a).
The results on the study of the correlation between FBXL19-AS1 expression and clinicopathological features showed in Table1 that the over expression of FBXL19-AS1 was significantly correlated with higher GGT (P = 0.046), higher AFP (P = 0.020) and worse TNM stage (P = 0.043), which was consistent with what we predicted in GEPIA2 [see Additional file 4].
To further verify the prognostic value of FBXL19-AS1 in HCC, we performed Kaplan-Meier survival analysis and log-rank test in 57 HCC patients with intact prognostic information. The results showed that HCC patients with elevated FBXL19-AS1 had shorter overall survival (P = 0.030) (Fig. 4b), hinting FBXL19-AS1 might be an important prognostic factor for HCC.
Expression of FBXL19-AS1 in plasma and its diagnostic value
FBXL19-AS1 expression level was also checked through qPCR in plasma samples, which were taken from 79 healthy subjects, 77 patients with hepatitis B, 80 patients with cirrhosis, and 92 patients with HCC (Fig. 4c). The results showed that the expression of FBXL19-AS1 in plasma was significantly higher in hepatitis B, cirrhosis, and HCC patients than in healthy subjects (P < 0.001). While the expression of FBXL19-AS1 in HCC patients was higher than that of hepatitis B patients (P = 0.016) and cirrhosis patients (P = 0.004). The main clinical features of the plasma samples were shown in Table 2. From which the expression of FBXL19-AS1 in each group was correlated with the plasma alanine aminotransferase (ALT) (P < 0.001), aspartate aminotransferase (AST) (P < 0.001), albumin (ALB) (P < 0.001), alpha fetoprotein (AFP) (P < 0.001), but not with gender, age, CEA or other biochemical indicators.
ROCs of FBXL19-AS1 in HCC, drawn to evaluate the diagnostic value, indicated FBXL19-AS1 was with moderate diagnostic ability to distinguish HCC patients from healthy people (AUC = 0.875, P < 0.001). Intergation of AFP, the most commonly detected biomarker for the diagnosis of HCC whose predictive value alone was rather low (AUC = 0.769, P < 0.001), with FBXL19-AS1, its preditive validity was significantly improved (AUC = 0.931, P < 0.001). (Table3, Fig. 4d-i).
78 miRNAs were obtained as the potential target miRNAs of FBXL19-AS1 through binding site search from miRcode V11 database. 24 miRNAs were selected after comparing with the 467 differentially expressed miRNAs screened from TCGA (P < 0.05). Subsequently, miRNA target genes were predicted and only 7 of the above 24 miRNAs met our requirements (has-miR-142-3p, hsa-miR-125a-5p, hsa-miR-216b-5p, hsa-miR-107, hsa-miR-17-5p, has-miR-20b-5p, has-miR-22-3p). Afterwards, we explored signaling pathways related with the 7 miRNAs through DIANA-miRPath. As shown in Fig. 5, the 7 miRNAs were essentially involved in the initiation and progression of many types of cancers. Among them, pathways that might be related to HCC were mTOR signaling pathway, Hedgehog signaling pathway, hepatitis B, hepatitis C, p53 signaling pathway, Wnt signaling pathway, GnRH signaling pathway, ErbB signaling pathway, PI3K-Akt signaling pathway, MAPK signaling pathway, etc.
Validation of miRNA expression based on GEO and TCGA database
To evaluate the expression of these 7 miRNAs in HCC, 9 GSE microarrays from GEO and dataset from TCGA were selected for verification [see Additional file 8]. Compared with normal liver tissues, hsa-miR-216b-5p (SMD = 0.593, P = 0.000, Fig. 6a), hsa-miR-107 (SMD = 0.729, P = 0.000, Fig. 6b), hsa-miR-17-5p (SMD = 0.502, P = 0.001, Fig. 6c) were up-regulated in HCC tissues, while has-miR-125a-5p (SMD = -0.947, P = 0.000, Fig. 6d), hsa-miR-22-3p (SMD = -0.563, P = 0.000, Fig. 6e) were down-regulated in HCC tissues. In addition, hsa-miR-20b-5p (SMD = 0.194, P = 0.217, Fig. 6f) and hsa-miR-142-3p (SMD = -0.425, P = 0.056, Fig. 6g) were not differentially expressed between normal tissues and HCC tissues, more studies with larger sample size were still needed.
We analyzed the correlation between the expression of FBXL19-AS1 and these 7 miRNAs based on the TCGA LIHC dataset. Significant correlations were not found in hsa-miR-216b-5p (r = -0.0010, P = 0.9830, Fig. 7a), hsa-miR-107 (r = -0.0504, P = 0.3027, Fig. 7b), hsa-miR-17-5p (r = 0.0715, P = 0.1438, Fig. 7c) or hsa-miR-125a-5p (r = 0.0303, P = 0.5357, Fig. 7d), but the relations were prominent in hsa-miR-22-3p (r = -0.2861, P < 0.001, Fig. 7e), hsa-miR-20b-5p (r = 0.0993, P = 0.0420, Fig. 7f) and hsa-miR-142-3p (r = -0.1435, P = 0.0032, Fig. 7g). Hence, we took hsa-miR-22-3p, hsa-miR-20b-5p and hsa-miR-142-3p for subsequent analyses.
LncRNA-miRNA-mRNA network construction
To further explore the potential downstream targets of these 3 miRNAs, three online bioinformatics servers (miRDB, miRTarBase and TargetScan) were used. There were 399 target genes of these 3 miRNAs simultaneously supported by all three databases. Then, we screened out 12,841 differentially expressed mRNAs in TCGA LIHC (P < 0.05). In addition, 12,194 mRNAs were predicted to be co-expressed with FBXL19-AS1 in HCC by cBioportal database (P < 0.05). Finally, 205 mRNAs were selected as targets through the intersection of the above three gene sets.
A new ceRNA network was formed among lncRNA (FBXL19-AS1), three miRNAs (hsa-miR-22-3p, hsa-miR-20b-5p, hsa-miR-142-3p) and 205 mRNAs (Fig. 8a). The diamond in the middle represented FBXL19-AS1, the gray triangles were miRNAs, and the circles were mRNAs. The circles in red meant the corresponding mRNAs were elevated in HCC, blue circles represented the decreased expression of related mRNA in HCC, deeper color indicated increased logFC, and larger size indicated smaller P value.
PPI network construction and screening of hub genes
We constructed a PPI network of these 205 mRNAs based on the STRING database and then visualized by Cytoscape 3.7.2. After removing the free nodes, the PPI network containing 158 nodes and 272 edges (Fig. 8b). Thereafter, hub genes identified by 12 algorithms (MCC、DMNC、MNC、Degree、EPC、BottleNeck、EcCentricity、Closeness、Radiality、Betweenness、Stress、ClusteringCoefficient) constituted a subnetwork with 9 nodes and 9 edges (Fig. 8c), which revealed the 9 hub genes (STAT3, CNOT7 BTG3, E2F1, TRIM37, YWHAZ, RBBP7, KIF23, ESR1) played important roles in the pathogenesis of HCC.
Functional analysis of 9 hub genes
GO and KEGG enrichment analyses were performed on the 9 hub genes (P < 0.05). Top 15 terms and pathways were selected for demonstration by P value. GO functional enrichment analysis revealed hub genes mainly enriched in transcription factor activity and transcriptional activator activity (Fig. 9a, Table 4). KEGG pathway analysis indicated the 9 hub genes might influence the occurrence and progression of HCC by participating in pathways such as hepatitis C, hepatitis B, microRNAs in cancer, cell cycle, viral carcinogenesis, and proteoglycans in cancer (Fig. 9b, Table 5). In addition, the 9 hub genes were also involved in pathways associated with non-small cell lung cancer, pancreatic cancer, breast cancer and other diseases.
Verification of ceRNA network
In order to establish a more reliable ceRNA network, we performed correlation analyses among FBXL19-AS1, 2 miRNAs and 9 mRNAs based on 370 HCC tissues and 50 normal tissues from TCGA LIHC dataset. Considering none of the 9 hub genes was correlated with hsa-miR-142-3p, it was not enrolled into the subsequent analyses. Significant correlations were found in FBXL19-AS1, hsa-miR-22-3p, hsa-miR-20b-5p, and 7 mRNAs, except STAT3 and CNOT9 (Fig. 10). Survival analyses of OS and DFS of the remaining 7 hub genes were performed by the Kaplan-Meier method in the GEPIA2 database. Remarkable survival differences existed in nearly all the 7 hub genes, except the DFS of YWHAZ (Fig. 11), which indicated all the 7 hub genes were associated with the prognosis of HCC. Therefore, we constructed a ceRNA network consisting of 1 lncRNA, 2 miRNAs, 7 hub genes, and seven lncRNA-miRNA-mRNA regulatory pathways (FBXL19-AS1/miR-22-3p/YWHAZ axis, FBXL19-AS1/miR-22-3p/ESR1 axis, FBXL19-AS1/miR-20b-5p/E2F1 axis, FBXL19-AS1/miR-20b-5p/BTG3 axis, FBXL19-AS1/miR-20b-5p/KIF23 axis, FBXL19-AS1/miR-20b-5p/KIF23 axis, FBXL19-AS1/miR-20b-5p/TRIM37 axis, FBXL19-AS1/miR-20b-5p/RBBP7 axis) (Fig. 12A). To further verify the reliability of the network, we analyzed the expression status of the elements in the network based on TCGA LIHC dataset. As shown in Fig. 12B, FBXL19-AS1, has-miR-20b-5p, has-miR-22-3p and 7 hub genes were all differentially expressed between normal tissues and HCC tissues.