Identification of metastasis related miRNAs in hepatoma carcinoma
A flowchart detailing the overall study is shown in Fig. 1. 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). 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 (Fig. 2A and 2B). 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 (Fig. 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 (Fig. 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% (Fig. 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 (Fig. 3C and 3D).
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 (Fig. 3E and 3F).
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 (Fig. 4A). Also, the prognostic value of miR-30c, miR-185, miR-323 were obtained using KM plotter (Fig. 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 (Fig. 5A-C). These genes were enrich in the terms of Transcriptional Regulation by MECP2, trans−Golgi Network Vesicle Budding by Reactome pathway analysis (Fig. 5D).
The PPI network of these miRNA targets is constructed by string (Supplementary Fig . 1), and MCODE (Molecular Complex Detection ) in Cystoscope was used to generate 3 stable subnetworks (Fig. 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 (Fig. 6D).