Identification of metastasis related miRNAs in hepatoma carcinoma
To construct a model to evaluate the risk of metastasis, we first obtained microarray data retrieved from Gene Expression Omnibus (GEO) database, and 241 samples (GSE=6578) were randomly divided into training group (n=169) and testing group (n=72). The Flow chart is illustrated in Figure 1. In the training group, 41 patients were diagnosed with HCC metastasis, while in the testing group, 19 were diagnosed. There was no big difference in metastasis between these two groups (p> 0.05). Then we used LASSO regression analysis to screen miRNAs highly associated with metastasis. The results of the training group with maximal AUC showed 8 potential metastasis related miRNAs, including miR-30c, miR-124a, miR-125b, miR-185, miR-206, miR-323, miR-325 and miR-1 (Figure 2A and 1B). After multivariate logistic regression analysis, only miR-30c, miR-185 and miR-323 were independent risk factors for HCC metastasis (Table 1).
Nomogram construction of hepatoma carcinoma
We developed a nomogram for predicting HCC metastasis, based on 3 independent risk factors (Figure 2C). In the training group, the nomogram performed well to distinguish HCC patients with metastasis from non-metastasis with a high AUC of 0.869 (95% CI 0.813-0.925); the sensitivity and specificity were as high as 80.49% and 78.91%, respectively (Figure 3A). The optimal cutpoint of 131.59 was calculated by Cutpoint R package, and the patients in testing group were stratified into the high group (total points >131.59) and the low group (total points ≤131.59). Consequently, the AUC in the testing group was 0.821 (95%CI 0.770-0.872) for predicting metastasis with a sensitivity of 48.5% and specificity of 92.3% (Figure 3B, Table 2). Moreover, the calibration plot for metastasis probability confer a good agreement between the prediction by our nomogram and actual observation in the training group and testing group, respectively (Figure 3C and 2D).
We further determined clinical applicability of nomogram using decision curve analysis (DCA) and clinical impact curve analysis (CICA). The results suggested our nomogram could have an optimal clinical net benefit (Figure 3E and 2F).
Association between nomogram and survival outcome
To investigate the correlation between nomogram and survival outcome, HCC patients without follow-up were removed from further analysis in the datasets. We calculated the total points of 211 cases (Alive group, n= 152; dead group, n= 59), and observed that the dead group had higher total points compared with alive group (Figure 4A). Also, the prognostic value of miR-30c, miR-185, miR-323 were obtained using KM plotter (Figure 4B). Kaplan Meier overall survival (OS) curves displayed these results: patients with high miR-30c level had a good prognosis compared with cases with low miR-30c level (n=76, HR= 0.55, p=0.042); no obvious difference was observed between low and high miR-185 levels; unfortunately, missing expression data for miR-323 in the database resulted in an unknown prognostic role.
Functional enrichment analysis of miRNA targets
We used the multiMiR package to predict miRNA-target interactions from 8 external databases, including DIANA-microT, ElMMo, MicroCosm, miRanda, miRDB, PicTar, PITA and TargetScan. For each miRNA, If the miRNA-target interactions were predicted in 3 different database, then these targets were screened out. We then performed functional enrichment analysis. These genes were enrich in the terms of transcription coregulator activity, phosphatidylinositol binding, transcription corepressor activity, kinase regulator activity by Molecular Function (MF) enrichment analysis ; in the terms of glutamatergic synapse, organelle subcompartment, filopodium by Cellular component (CC) enrichment analysis; in the terms of cell-cell adhesion via plasma-membrane adhesion molecules, embryonic organ development by Biological Process (BP) enrichment analysis (Figure 5A-C). These genes were enrich in the terms of Transcriptional Regulation by MECP2, trans−Golgi Network Vesicle Budding by Reactome pathway analysis (Figure 5D).
The PPI network of these miRNA targets is constructed by string (Supplementary Figure 1) , and MCODE (Molecular COmplex Detection ) in Cystoscope was used to generate 3 stable subnetworks (Figure 6A-C). Hub genes for subnetwork 1 included “CBLB, NEDD4, CUCL2”, subnetwork 2 included “PICALM, VAMP2 and NECAP1” and subnetwork 3 included “ACTR1A, DYNLL2 and RAB7A”. The prognostic value of NEDD4, CUCL2, PICALM, VAMP2, NECAP1, DYNLL2 and RAB7A in HCC were observed except for CBLB and ACTR1A using KM plotter (Figure 6D).