3.1 Identification of ferroptosis-related subtypes in HCC
To mine the DEFRGs between HCC and normal samples, we first screened the DEGs between HCC and non-tumorous samples in the TCGA-LIHC cohort. As shown in Fig. 1A and Supplementary Table 1, 5622 DEGs were filtrated on the basis of the criteria of |log2(FC)| > 1 along with adjusted p-value < 0.05, including 2462 up-regulated genes, as well as 3160 down-regulated genes. The TOP100 DEGs were displayed in Fig. 1B. We then took the intersection of 5622 DEGs and 167 FRGs to obtain 36 DEFRGs (Fig. 1C). The above 36 genes were further included in the univariate Cox assessment and 15 genes relevant to overall survival (OS) in individuals with HCC (p < 0.05) were identified, namely STMN1, MAFG, TXNIP, SLC7A11, RRM2, PCK2, MT3, IL33, HMOX1, GPT2, CBS, CAPG, AURKA, ASNS and ALB (Fig. 1D). We analyzed the Pearson correlations between these 15 genes. The results showed a positive correlation between STMN1 and RRM2, STMN1 and AURKA, SLC7A11 and MAFG, RRM2 and AURKA, PCK2 and GPT2, and PCK2 and ALB (Fig. 1E). To further understand the interactions between the proteins encoded by these 15 genes, we constructed a PPI network containing 13 nodes with 16 edges (Fig. 1F).
To explore the ferroptosis-related subtypes in HCC, 370 HCC samples were clustered using unsupervised cluster analysis based on the expression of 15 DEFRGs. Two subtypes, cluster 1 along with cluster 2, were identified from the results (Fig. 2A). Cluster 1 included 140 cases and cluster 2 included 230 cases. Survival analysis exhibited remarkable differences in prognosis between the two ferroptosis-related subtypes, and cluster 2 had a decent survival advantage (Fig. 2B). To further explore possible mechanisms of the survival differences between the two clusters, GSVA enrichment assessment was conducted on these two subtypes to assess the differences with regard to function along with pathways. The results manifested that cluster 1 was abundant in oncogenic cascades, consisting of the P53 signaling cascade, and some metabolic cascades (Fig. 2C). Cluster 2 was primarily enriched in glucose metabolic cascades and immune-related cascades, consisting of PPAR signaling cascade and adipocytokine signaling cascade (Fig. 2C).
3.2 FIMRGs in HCC
We then compared the infiltrating immune cells and immune features in Cluster 2 and Cluster 1 using ssGSEA. The results manifested that the enrichment fraction of most immune cells and immune features was greater in cluster 1 compared to cluster 2 (Fig.2D). The correlation between the different immune cells and immune features was shown in Fig.2E. To assess further the genes linked with the immune micro-environment, we included the above infiltrating immune cells and immune features into the WGCNA analysis. No obvious outliers were excluded by cluster analysis (Fig.3A). 5 was selected as the optimal soft threshold with an R2 = 0.85 (Fig.3B). On the basis of the optimal soft threshold, we stratified the genes into different modules according to the dynamic tree-cutting algorithm. After merging, a total of 7 modules were generated (Fig.3C). Correlations of the modules with traits were computed and brown modules were associated with most traits as key modules (Fig.3D). The 1436 genes in key module were defined as immune microenvironment-linked genes in HCC. To further identify FIMRGs, we assessed the relationships of the 1436 immune microenvironment-linked genes with the 15 prognosis-linked DEFRGs. After selecting the relationship pairs with correlation coefficient > 0.5 along with p-values < 0.05, the corresponding genes were combined and de-duplicated to finally obtain 55 genes defined as FIMRGs in HCC (Supplementary Table 2).
To investigate further the biological functions of the above 55 genes and the involved pathways, GO and KEGG enrichment analysis was carried out. A total of 291 GO items (including 279 BP items, 7 CC items and 5 MF items) and 10 KEGG pathways were enriched (Supplementary Table 3). Top5 GO entries under each category were listed in Fig.3E. We established that these 55 genes were most linked with ‘immune response-modulating signaling cascade’, ‘immune response-modulating cell surface receptor signaling cascade’, ‘activation of immune response’, ‘immune response-activating cell surface receptor signaling cascade’, and ‘immune response-activating signal transduction’. In addition, for KEGG pathways, these 55 genes were enriched in ‘Osteoclast differentiation’, ‘Neutrophil extracellular trap formation’, ‘C-type lectin receptor signaling cascade’, ‘HIF-1 signaling cascade’, ‘Ferroptosis’, ‘Phagosome’, ‘Chemokine signaling cascade’, ‘Rap1 signaling cascade’, and ‘B cell receptor signaling cascade’ (Fig.3F).
3.3 Creation and verification of the risk score model based on FIMRGs in HCC
In the following, the univariate Cox assessment was fulfilled to identify prognosis-linked genes on the basis of 55 FIMRGs in the training set. Nine out of the 55genes were identified as closely related to patients' OS in the training set (Fig.4A). Then, the nine genes were submitted further to stepwise multivariate Cox regression assessment. These 9 FIMRGs (PTPRO, CRHBP, HMOX1, FABP5, TREM2, ADAMTS14, LILRA6, CAPG, and RGS2) were identified as prognostic genes and each regression coefficient was determined. Then we generated a risk score model. To assess the availability of the model, the patients were split into high-HCC and low-HCC risk groups on the basis of the median of the risk score in the training data set. Then K-M survival analysis was performed. The results suggested that patients in high-HCC risk group had remarkably worse survival rates versus those in low-HCC risk group (Fig.4B). The AUC value of the ROC curve in the training cohort was 0.646, exhibiting a decent accuracy (Fig.4C). Fig.4D displayed the ranked risk scores along with the survival status of individual patients in the training cohort. The distribution map of survival status revealed that patients suffered a greater mortality risk with increasing risk scores. The expression of nine prognostic genes in low-HCC and high-HCC risk groups was exhibited in Fig.4E. The heatmap results suggested PTPRO, HMOX1, FABP5, TREM2, ADAMTS14, LILRA6, CAPG, and RGS2 were highly expressed in the high-HCC risk group. CRHBP was highly expressed in the low-HCC risk group.
We utilized an independent data set with 221 HCC subjects from the GSE14520 to prove the accuracy of the model externally. Consistent with the TCGA-LIHC cohort, subjects with high-HCC risk exhibited worse OS in the validation set (Fig.4F). Correspondingly, the AUC value of the ROC curve in predicting prognosis was 0.624 (Fig.4G). Fig.4H exhibited the survival status on the basis of the patients with increased risk scores and the heatmap result suggested that PTPRO, HMOX1, FABP5, TREM2, ADAMTS14, LILRA6, CAPG, and RGS2 were highly expressed in the high-HCC risk group (Fig.4I). CRHBP was highly expressed in the low-HCC risk group (Fig.4I). These results further affirmed the validity of our risk model in forecasting the survival of HCC patients. Taken together, these results demonstrated that FIMRGs were associated with the prognosisi of HCC.
3.4 Independent prognostic worthiness of the risk score based on FIMRGs in HCC
To clarify the role of the gene signature in HCC progression, we researched the correlation of the risk score with clinical characteristics in the training set. As indicated in Fig.S1A, the risk score was linked with stage. The expression of prognostic genes under different clinical factors and our risk score distribution of patients with different clinical features in the training set, as well as the verification set were shown in Fig.S1B-C. Subsequently, we conducted univariate along with multivariate Cox regression assessments to assess if the risk score was independent of other clinical parameters (age, stage, gender, and AFP) as a predictive factor for individuals with HCC in the training set and verification set respectively. The univariate assessment exhibited that the stage and risk score was relevant to patients' prognosis (Fig.5A, 5C). The multivariate Cox regression assessment demonstrated that the risk score could be used independently to estimate the survival of patients (Fig.5B, 5D). Furthermore, to predict the survival outcome of HCC, a prognostic nomogram was established through Cox regression model analysis according to two significant indepentent indicators (stage and risk score) (Fig.S3A-D). The C-index value was 0.682 in the training set and 0.636 in the validation set. The time‐dependent AUC was > 0.6 for the prediction of OS within 5 years in both the training and validation cohorts, indicated relatively good discrimination by the nomogram. Then, We plotted the ROC curves (Fig.S3E-F)
3.5 Differential analysis of TME and sensitivity to chemotherapy in high-HCC and low-HCC risk groups
We detected that nine FIMRGs were positively correlated with many immunomodulators (Fig.S2A). Expression of MHC was strongly correlated with antigen presentation and processing capacity. Chemokines along with their receptors facilitate the mobilization of effector TIICs, consisting of CD8+ T cells, antigen-presenting cells, and TH17 cells. Exploring the relevance of nine FIMRGs to immune modulators aims to elucidate the immune role of nine genes in TME.
To explore possible mechanisms for the disparity in prognosis between the high-HCC and low-HCC risk groups, we first performed an analysis of infiltrating immune cells (Fig.6A-B). The data revealed that the relative abundance of memory B cells, CD4 Memory T cells, follicular helper T cells, Macrophages (M0), as well as regulatory T cells (Tregs), was remarkably greater in the high-HCC risk group compared to the low-HCC risk group (Fig.6C). Meanwhile, the relative abundance of naive B cells, Monocytes, resting memory CD4+T cells, Macrophages (M1), activated NK cells, Macrophages (M2), and resting Mast cells was remarkably lower in the high-HCC risk group compared to the low-HCC risk group (Fig.6C). Hence, we speculated that there was a significant distinction in the immune microenvironment between the high-HCC and low-HCC risk groups.
The activity of the cancer-immune cycle is a straightforward and integrated expression of the function of the chemokine system as well as other immunoregulators[24, 27], and the activity of these steps determines the fate of tumor cells. In the high-HCC risk group, the activity of all steps of the immune cycle was found to be upregulated (Fig.6D), including seven steps: cancer cell antigen release (step 1), cancer antigen presentation (step 2), initiation and activation (step 3), immune cell moving toward the tumor (step 4) (CD8 T-cell mobilization, macrophage mobilization, Th1 cell mobilization, NK cell mobilization, DC mobilization, and recruitment of TH17), immune cell invasion into the tumor (step 5), T cell recognizing of cancer cells (step 6) along with killing cancer cells (step 7). All of these showed a higher immune infiltration status in the high-HCC risk group. Theoretically, patients in the high-risk group should have a higher response to immunotherapy because subjects in the high-HCC risk group have an inflammatory TME. Subsquently, we used Tumor Immune Dysfunction and Exclusion (TIDE) (http://tide.dfci.harvard.edu/) to predict immunotherapy response in the high-risk and low-risk groups, and additionally used the submap module in GenePattern Modules (https://www.genepattern.org/modules) using data from melanoma as a template . The similarity between the high-risk and low-risk groups and the two immune checkpoint inhibitor treatment data was predicted. It was found that the predicted TIDE scores were significantly higher in the high-risk group than in the low-risk group, and the data from the high-risk group were significantly similar to the data from the PD1 treatment response group (Fig.S4).
Next, we contrasted the expression of 20 important immune checkpoints between the high-HCC and low-HCC risk groups, the expression of 17 immune checkpoints were remarkably different. As exhibited in Fig.6E, the expression of PD-L1, TIM-3, PD-1, CTLA-4, as well as LAG-3 were remarkably higher in the high-HCC risk group, and only the expression of ADORA2A was decreased, indicating that the high-HCC risk group might have a better response to immune checkpoint inhibitors.
As expected, the enrichment scores for all immunotherapy-positive genetic characteristics were higher in the high-HCC risk group (Fig. 6F). In addition, we further analyzed the sensitivity of high-HCC and low-HCC risk groups to drugs. As revealed in Fig.6G, the low-HCC risk group was more sensitive to five chemotherapy drugs (Axitinib, Rapamycin, Temsirolimus, Docetaxel, and Metformin), suggesting that patients in the low-HCC risk group would benefit more from chemotherapy.
3.6 Multi-Omics analysis of prognostic genes
As exhibited in Supplementary Table 1, the expression of PTPRO, CRHBP, HMOX1, LILRA6, and RGS2 was reduced in HCC tissues, and the expression of TREM2, FABP5 and ADAMTS14 was up-regulated in HCC tissues. To assess further the biological mechanisms underlying the aberrant expression of prognostic genes in HCC tissues, we performed a multi-omics analysis of these genes. 8 genes were mutated in HCC tissues, mainly missense mutations (Fig.7A). As these eight genes presented differential expression in most HCC samples, we hypothesized that the aberrant expression of prognostic genes was not mutation-related. We further investigated the abnormal methylation levels. As shown in Fig.7B. TREM2 was hypermethylated in HCC tissues. Aberrant methylation is known to play an instrumental role in the modulation of expression levels, so we then analyzed the relevance between methylation levels and expression levels. The result indicated that methylation of TREM2 was not correlated with the expression (Fig.7C). Next, we contrasted copy number variation (CNV) in terms of both heterozygosity and homozygosity. As revealed in Fig.7D-E, we detected the presence of amplified variants in the RGS2, FABP5, and TREM2. Nevertheless, there was no remarkable correlation between CNV and expression level variation in RGS2, FABP5, and TREM2 (Fig.7F). Therefore, we hypothesize that the aberrant expression regulation of prognostic genes in HCC does not occur at the genomic level and requires further investigation.
3.7 Establishment of a prognosis-related miRNA-mRNA regulatory network
We then explored the post-transcriptional regulatory mechanisms of prognostic genes. First, we analyzed the DE-miRNAs between HCC samples and normal samples. A total of 399 DE-miRNAs, containing 156 up-regulated miRNAs and 234 down-regulated miRNAs, were identified (Supplementary Table 4, Fig.8A). Top100 DE-miRNAs were shown in the heatmap (Fig.8B). Through the miRbase database, we obtained 67 miRNAs associated with prognostic genes. 13 DE-miRNAs associated with prognostic genes were acquired after intersecting with DE-miRNAs (Fig.8C). 5 miRNAs linked with the overall survival of individuals with HCC were further authenticated by univariate Cox analysis (Fig.8D). We finally constructed a miRNA-mRNA regulatory network comprising four prognostic genes and five prognostic miRNAs (Fig.8E). In the network, the hsa-miR-642a-5p regulated the expression of HMOX1 and PTPRO. The hsa-miR-1304-5p regulated the expression of HMOX1. The hsa-miR-133b regulated the expression of PTPRO. The hsa-miR-125a-5p regulated the expression of ADAMTS14. The hsa-miR-3162-5p regulated the expression of RGS2.
3.8 The potential drugs targeting the prognostic genes
To further explore potential drugs targeting prognostic genes, we extracted four drugs targeting HMOX1 from the DGIdb database, including ASPIRIN, STANNSOPORFIN, SUNITINIB, and SORAFENIB. In addition, one drug targeting RGS2, HALOPERIDOL, and one drug targeting CAPG, VINCRISTINE were also identified. We also discovered that ALCOHOL targeted CRHBP (Fig.8F).
3.9 FIMRGs differentially express in protein level
In HPA, the histological expression of FIMRGs in normal and tumor tissues was shown in Fig.9. There were significant differences in the expression of three genes at the protein level between normal and liver cancer tissues. The protein expressions of FABP5, TREM2 and CAPG were significantly up-regulated. However, no significant expression differences were found in PTPRO, CRHBP and HMOX1. The protein expressions of PTPRO, CRHBP and HMOX1 were lower in normal and tumor tissues. ADAMTS14, LILRA6 and RGS2 failed to find data in the database.