Expression of m6A modification related genes of HCC patients and their associations with clinical-pathologic parameters
First, mRNA expression of 13 m6A modification related genes were download from TCGA and compared between HCC patients and normal controls. As was shown at figure 2A and 2B, significantly higher expression of FTO, YTHDC1, YTHDC2, ALKBH5, KIAA1429, HNRNPC, METTL3, RBM15, YTHDF2, YTHDF1 and WTAP were found at HCC compared to normal tissues (all p<0.001). Interestingly, we also found that expression of most of the 13 m6A modification related genes seemed to be lower than that of other 32 kinds of tumors. Besides, most of the 13 m6A modification related genes were positively correlated with each other (figure 2C). Moreover, genetic change, such as missense mutation, truncating mutation, amplification, deep deletion, diploid, gain were observed in near 80% of the HCC patients (figure 2D), Specifically, each HCC patient may have one or more kinds of genetic changes. The genetic rate of WTAP, KIAA1429, RBM15, METTL3, METTL14, ALKBH5, YTHDC1, YTHDC2, HNRNPC, YTHDF1, YTHDF2, FTO and ZC3H13 were 7%, 4%, 17%, 40%, 5%, 5%, 7%, 8%, 18%, 11%, 9%, 13%, 17%, respectively, suggesting that higher expression of m6A modification related genes may be the result of genetic changes in related genes. Taken together, these results indicated that m6A modification related genes played important roles in HCC.
Prognostic value of m6A modification related genes in HCC cases
After discovering the expression of m6A modification related genes were associated with tumor histologic grade and TNM stage, we further analyzed their prognostic values. Univariate analysis showed that higher expression of YTHDF1, WTAP, HNRNPC, RBM15, METTL3, KIAA1429, YTHDC1, YTHDF2 and lower expression of ZC3H13 were statistically related to poorer OS of HCC patients (all p<0.05, supplementary figure 1A); multivariate analysis showed that expression of YTHDF1, WTAP, HNRNPC, RBM15, METTL3, KIAA1429 and YTHDF2 still remained significantly related with OS after adjusting for gender, age, histologic grade, T stage, N stage, M stage and TNM stage (all p<0.05, supplementary figure 1B-1J).Then, the prognostic values of m6A modification related genes for recurrence of HCC patients were also analyzed. Univariate analysis indicated that overexpression of YTHDF1, WTAP, HNRNPC, RBM15, METTL3, YTHDC1 and YTHDF2 were statistically related with shorter DFS (all p<0.05, supplementary figure 2A); and, multivariate analysis showed that expression of YTHDF1, HNRNPC, RBM15, METTL3 and YTHDF2 were still statistically related with DFS after adjusting for gender, age, histologic grade, T stage, N stage, M stage and TNM stage (all p<0.05, supplementary figure 2B-2H). These results strongly confirmed the important roles played by m6A modification related genes in HCC.
Development of risk signature with 5 m6A modification related genes and its association with clinical-pathologic parameters
To better explore the prognostic value of m6A modification related genes, a risk signature was developed. Based on the results of univariate analysis (figure 3A), ZC3H13, YTHDF1, WTAP, HNRNPC, RBM15, METTL3, KIAA1429, YTHDC1 and YTHDF2 were associated OS and were considered as prognostic-related genes. Then, LASSO analysis was used to further screen the prognostic-related genes. In the end, 5 genes, including YTHDF2, YTHDF1, METTL3, KIAA1429 and ZC3H13, were used to develop the risk signature (figure 3A, 3B). The risk score was then constructed based on the coefficients weighted by LASSO analysis and calculated as follows: risk score= (0.07*YTHDF2) + (0.02*YTHDF1) + (0.11*METTL3) + (0.04*KIAA1429) - (0.1*ZC3H13). We calculated risk score for every HCC cases and assigned them into high risk group and low risk group on the basis of the median risk score. Expression of YTHDF2, YTHDF1, METTL3, KIAA1429 tend to be higher in patients with high-risk score; and expression of ZC3H13 seemed to be higher in patients with low-risk scores (figure 3C). Distribution of histologic grade, T stage, TNM stage were significantly different between high risk subgroup and low risk subgroup (all p<0.05, figure 3C). High risk subgroup contained more patients with advanced histologic grade, T stage and TNM stage compared to patients of the low risk subgroup. Lastly, patients in the high risk subgroup had poorer OS (median OS time: 2.46 vs 5.79 years, HR=1.98, 95%CI: 1.39-2.83, P<0.001, figure 3D) and shorter DFS (median DFS: 1.07 vs 2.97 years, HR=3.83, 95%CI: 2.56-5.90, P<0.001, figure 3E) than that of patients of the low risk subgroup, which was consistent with previous result.
Prognostic value of risk signature for OS and DFS of HCC cases
The risk signature was found to be associated with clinical-pathologic parameters. We next performed univariate and multivariate analysis to analyze its prognostic value. Based on the univariate analysis, T stage, M stage, TNM stage and risk signature were statistically related with OS of HCC patients (all p<0.05, figure 4A). The risk signature still remained statistically related with OS after adjusting for T stage, M stage and TNM stage by multivariate analysis. In multivariate analysis, after adjusting for TNM stage, the risk signature was still significantly related with OS (p<0.01, figure 4B). Similarly, univariate analysis also showed that T stage, TNM stage and risk signature were statistically related with DFS of HCC patients. In univariate analysis, T stage, TNM stage and the risk signature were also significantly associated with DFS in HCC patients (all p<0.001, figure 4C). By incorporating these factors into multivariate analysis, the result suggested only the risk signature was statistically related with DFS (p<0.001, figure 4D). To conclude, these results indicated that risk signature was independent prognostic factor for OS and DFS of HCC patients.
Next, we used time-dependent ROC cure analysis to analyze the predictive value of risk signature for HCC patients. As were shown at figure 5, the AUC of risk signature for predicting 1-, 3- and 5-year OS were 0.765, 0.73 and 0.678, respectively, which exhibited better predictive efficiency compared to TNM stage, YTHDF2, YTHDF1, METTL3, KIAA1429 and ZC3H13 (figure 5A, C, E). Likewise, the AUC of risk signature for predicting 1-, 3- and 5-year DFS were 0.695, 0.643 and 0.68, respectively, which also showed better predictive accuracy than TNM stage, YTHDF2, YTHDF1, METTL3, KIAA1429 and ZC3H13 (figure 5B, D, F).
Validation of risk signature
To independently test the applicability of the signature, 232 HCC patients with available OS information from the ICGC portal (https://dcc.icgc.org/projects/LIRI-JP) was further used to examine the applicability of the signature. Risk score for every patient was computed. Similarly, the signature could effectively stratify high risk HCC patients with poorer OS and low risk patients with better OS (HR=2.309, 95%CI: 1.302-4.369, P=0.006, figure 6A). Moreover, the AUC of risk signature for predicting 1-, 3- and 5-year OS were 0.7, 0.74 and 0.714 (figure 6B), respectively, which convincingly suggested the good discrimination and prediction of our signature.
Correlation of risk signature with tumor-infiltrating immune cells in TCGA and ICGC HCC cohort
CIBERSOR was used to calculate 22 kinds of infiltrating immune cells in patients with different risk scores. In TCGA HCC cohort, significantly higher proportions of macrophages M0 cells, memory B cells, follicular helper T cells and neutrophils were found to be enriched in HCC patients with high risk score, while significantly higher proportions of resting memory CD4 T cells and monocytes were found to be enriched in HCC patients with low risk score (all p<0.05, figure 7A). In ICGC HCC cohort, significantly higher proportions of macrophages M0 cells and Treg cells were found to be enriched in HCC patients with high risk score, while significantly higher proportions of naive B cells and gamma delta T cells were found to be enriched in HCC patients with low risk score (all p<0.05, figure 7B). These results suggested that the risk signature was significantly associated with tumor-infiltrating immune cells, and different kinds of infiltrating immune cells in patients with different risk score may contribute to their different prognosis.
Risk signature as indicator in sorafinib treatment response for HCC patients
To investigate the association between risk signature and sorafenib treatment response, we calculated risk score for each HCC patients treated with sorafenib of GSE109211, which contained 21 sorafenib treatment responders and 46 non-responders. Significantly lower risk scores were found at sorafenib treatment responders compared to that in non-responders (P < 0.001, figure 8A). Moreover, the AUC for predicting sorafinib treatment response was 0.794 (figure 8B). Taken together, the risk signature may be served as an indicator for sorafinib treatment response in HCC patients.
Correlation of risk signature with anti-PD-1 immunotherapy
As a major breakthrough in cancer therapy, immunotherapies represented by immunological checkpoint blockade (PD-1/L1 and CTLA-4) has proved promising clinical efficacy and previous study have proved that combination treatment with anti-PD-1 antibodies and Sorafenib exhibited a more potent antitumor effect, but only a small number of patients could achieved durable responses (31, 32), so in the present study, we also explored whether the risk signature could predict patients’ response to immune checkpoint blockade therapy in an anti-PD-1 cohort of GSE78220. Encouragingly, patients with lower risk score had better OS than patients with higher risk score (HR=3.81, 95%CI: 1.13-11.08, p=0.03, figure 9A). Besides, despite there was no statistical difference, lower risk score was found at patients with complete immunotherapeutic response compared to that in patients with partial response and patients with no response, and lower risk score was also found in alive patients treated with anti-PD-1 than that in patients of death, which may due to the limitation number of patients in the cohort (figure 9B-9C). Moreover, the AUC of the risk signature for predicting 1 year-, 1.5-year and 2-year OS of patients with anti-PD-1 immunotherapies was 0.669, 0.725 and 0.639 (figure 9D). In a word, the aboved results strongly indicated that risk signature was significantly correlated with response to anti-PD-1 immunotherapy, which may be used as a new biomarker for predicting the response to anti-PD-1/L1 immunotherapy.