3.1 Competing endogenous network construction and GO enrichment analysis
From the genes downloaded from the exoRbase database, we find the difference genes of mRNA, lncRNA, circRNA between the tumor and non-tumor tissues respectively. And then explored the relationships among mRNA, lncRNA, microRNA, circRNA, pairs of lncRNA and miRNA, miRNA and mRNA and circRNA and miRNA were used to develop ceRNA network by cytoscape. Totally, 398 pairs of lncRNAs-miRNAs, 681 pairs of miRNAs-mRNAs and 200 pairs of circRNAs-miRNAs were confirmed. Finally, 10 DElncRNAs, 72 DEmiRNAs, 48 DEmRNAs and 15DEcircRNAs were involved to construct ceRNA network.(Fig. 1A) To better understand the biological consequences of the exosome-related genes, we carried out GO enrichment analysis. The Go enrichment result illustrated that the exosome- related genes were significantly enriched in BP and MF, which were related to protein localization to nucleus and phosphatase regulator activity (Fig. 1B). All these signaling pathways are reported to take significant roles in the development of cancer. The results may indicate that abnormal expression of exosome-related genes might influence the pathway activation and mitosis and further lead to HCC.
3.2 Construction an exosome-related genes prognostic model in TCGA
Intersect the mRNAs found in ceRNA network with TCGA genes, we got 48 DEGs. Then, we carried out a univariate Cox regression and 18 DEGs were found to be significantly correlated with OS in HCC patients (P < 0.01). And a lasso analysis was carried out to narrow down the gene range. We have selected 6 DEGs, which appeared more than 500 times of 1000 cycles totally. At last, a multivariate Cox regression was carried out, and five genes (XPO1, IFI30, FBXO16, CALM1, MORC3) were chosen to establish the prognostic model as shown below: PI = (0.469 * expression value of XPO1) + (0.141* expression value of IFI30) + (0.218* expression value of FBXO16) + (0.146* expression value of CALM1) + (-0.343 * expression value of MORC3). The best cut-off value of 2.06 for the risk score was got by the software of X-tile. Then HCC patients were divided into low-risk and high-risk groups according to the best cut-off value. The overall survival of the high-risk group was significantly poorer than that of low-risk group according to the Kaplan–Meier curve (P < 0.001, HR = 3.51) (Fig. 2A). The area under the time-dependent ROC curves (AUCs) for 1-, 2‐ and 3‐year overall survival (OS) were 0.74, 0.71, 0.7 and 0.68 respectively, indicating a good predictive performance of this prognostic model (Fig. 2C).
3.3 External validation of exosome-related genes prognostic model in ICGC
Samples from ICGC were further used to validate the predictive power of this prognostic mode. Same way as above, patients were divided into low- and high-risk groups with the same cut-off value of 2.06 and the survival rate of HCC patients in low-risk group was higher than that in the high-risk group (P = 0.009, HR = 2.65) (Fig. 2B). And the AUCs of the five-gene prognostic model were 0.63, 0. 58 and 0.63 for the 1-, 2- and 3‐year survival times (Fig. 2D), which showed that this model is reliable.
3.4 Establishment a predictive nomogram
Univariate Cox regression analysis illustrated that T.N.M. stage and risk-score were both meaningful factors (Fig. 3A). Multivariate Cox regression analysis showed T.N.M. stage and risk-score were both independent prognostic factors (Fig. 3B).So, the nomogram was built by including T.N.M. stage and the risk score (Fig. 3C). Assign a score to the independent factors of each level, and calculate the total score by summing the scores of each person. By transforming the relationship of the total score, the survival rate of 1-, 2- and 3-years can be obtained. The calibration map was used for internal verification of this nomogram, which showed that there was a better conformity between the predicted results and the realistic observation results (Fig. 3D-F).
3.5 Computation of immune cell type fractions
A new tool called CIBERSORT is used to directly analyze the types and expression of immune cells in the tissue, combined with the Leukocyte signature matrix (LM22). And we assessed the differences of twenty-two types of immune cells between the two groups of HCC patients. And 2 immune cells (T cells regulatory, P = 0.043, Macrophages M2, P = 0.048) were different between the two groups (Fig. 4A). From our result, we may conclude that the poor prognosis relates to the high content of Tregs and M2. As shown in Fig. 4B, the twenty-two immune cell proportions between high-risk and low-risk groups were different, and the relationship among immune cells showed in Fig. 4C.
3.6 Evaluation of the immunotherapy effect
Blockade of immune checkpoints has been found to be a promising treatment method of many cancers. Subsequently, we studied the different expression of key immune checkpoint molecular including PD-L1, PD-L2, TIGIT and IDO1. As shown in Figs. 5A-5D, high-risk HCC patients seemed to express higher immune checkpoint molecules than low-risk patients (p < 0.05). From our research, we may conclude that the poor prognosis is related to the high expression of exosomes-related genes.