Establishment and validation of the ERS model
By comparing the gene expression profiles of 374 tumor tissues and 50 normal liver tissues from the TCGA database, we identified 35 differentially expressed ERS-related genes, including 24 upregulated genes and 11 downregulated genes (Fig. 2a) in tumor tissues. Using univariate Cox regression, nine ERS-related genes (p<0.05) associated with OS were found among the 35 DEGs (Additional Table 1). Following the application of multivariate Cox regression analysis based on these ERS-related genes, a prognostic model was developed by defining four sERS-related genes (Fig. 2b). The format for calculating the risk score was as follows: ERSrisk=0.24550*STC2+0.27670* DDX11+0.11410* MAGEA3+0.38278* BRSK2. The ERS risk was selected as the cutoff value for the ROC curve, and the patients were separated into high-risk and low-risk groups. Calibration curves were used to evaluate the discrimination and calibration of our model (Additional Fig. 1). The ERS model C-indexes were 0.646 in the TCGA-LIHC cohort and 0.653 in the ICGC-JP-LIRI cohort. K‒M survival analysis revealed that patients in the high-risk group had a poor prognosis, while those in the low-risk group had a relatively good prognosis (Fig. 2c).
To further assess the predictive ability of the model for HCC, a stratified analysis was conducted based on pathological and clinical characteristics such as TNM stage, age, sex, and tumor stage (Additional Fig. 2a). Patients with a high ERSrisk exhibited a greater likelihood of early mortality than did those with a low ERSrisk (Additional Fig. 2b). The area under the curve (AUC) values for 1-, 3-, and 5-year OS predictions were superior to those of other clinical characteristics (Fig. 2d). Variations in ERSrisk were observed between T1 and T2~T4, as well as between stage I and stage II~III, but not in advanced HCC (Fig. 2e).
Next, we performed a comparative analysis of mutations in the high-risk and low-risk groups (Additional Fig. 3a and b). We found that the most significant difference in ERSrisk was TP53. The ERSrisk of the mutant type was significantly greater than that of the wild type (Additional Fig. 3c). Moreover, the expression levels of MAGEA3 and DDX11 were increased (Additional Fig. 3d), which suggested that ERS may be closely related to TP53 mutation.
In addition, we performed an intergroup comparative analysis of 198 drugs. After cross-validation
with the TCGA and ICGC databases, we found that patients in the high-risk group were more
sensitive to five drugs, namely, ubrosertib, taselisib, pictilisib, OSI.027, and doramapimod (Additional Fig. 4a), but more resistant to 12 drugs, namely, BI.2536, ERK_2440, ERK_6604, mitoxantrone,
PD0325901, RO.3306, SB216763, SCH772984, selumetinib, topotecan, trametinib, and VE.822 (Additional Fig. 4b). Interestingly, four of the five sensitive drugs targeted the PI3K/mTOR signaling pathway; however, most of the resistant drugs targeted the ERK/MAPK signaling pathway (Table 1). This meant that even if they were all proliferation-related drugs, patients with high ERSriskERS risk were more likely to receive drugs targeting the PI3K/mTOR signaling pathway rather than the ERK/MAPK signaling pathway.
ERS model correlated with macrophages and endothelial cells
Furthermore, we analyzed the differences in infiltrating immune cells between the high-risk and low-risk groups. The results revealed significant differences in a variety of immune cells, including T cells, M0 macrophages, M1 macrophages, dendritic cells, neutrophils, and mast cells. The number of M0 macrophages in the high-risk group was significantly greater than that in the low-risk group, while the number of M1 macrophages was less than that in the low-risk group. Similar results were also observed in the ICGC-JP-LIRI cohort (Fig. 3a). Moreover, ssGSEA was used to assess the correlation between immune cells and ERSrisk, and it was found that ERSrisk was significantly related to macrophages (Additional Fig. 5); therefore, we believe that ERS is closely related to macrophages in the tumor immune microenvironment of HCC.
Additionally, we evaluated the relationship between immune checkpoints and ERS. The response to PD-1/CTLA4 immunotherapy in patients with different immunophenotypes was compared between the high-risk and low-risk groups in the TCGA-LIHC cohort. The low-risk group had a greater non-response rate to anti-CTLA4/PD-1 immunotherapy, suggesting that the high-risk group may have a better response to anti-CTLA4/PD-1 therapy (Fig. 3b). Moreover, the TIDE score of the high-risk group was lower than that of the low-risk group, while the PD-L1 score was higher than that of the low-risk group (Fig. 3c), which indicated that the high-risk group had a lower risk of evading immunotherapy and may be more effective in immune checkpoint treatment.
We integrated the single-cell data from GEO and removed batch effects. After a standardized procedure, all cells were divided into 23 clusters. Cell type annotations were verified through the “singleR” R package and cellmarker2 database (Additional Fig. 6). The 23 clusters were labeled with eight cell types that were assigned based on the expression of cell type-specific markers: T cells, NK cells, endothelial cells, malignant cells, myeloid-derived cells, HSCs, B cells, and epithelial cells (Fig. 4a). Subsequently, the gene tag was strategically positioned, leading to the identification of significant expression of STC2 within the endothelial cell population (Fig. 4b). Consequently, the endothelial cell population was subsequently extracted, followed by in-depth subpopulation analysis. The endothelial cells were partitioned into ten distinct clusters and subsequently classified into five subtypes: PLVAP+ECs, CD9+ECs, IGFBP3+ECs, PLPP3+ECs, and TFF3+ECs (Fig. 4c; Additional Fig. 7), with predominant expression of STC2 observed in PLVAP+ ECs (Fig. 4d).
Since we detected significant differences in macrophages between the high-risk and low-risk groups, we reclassified myeloid-derived cells and divided them into six cell types: pDCs, DC1s, DC2s, TAM1s, TAM2s, and monocytes (Fig. 4e; Additional Fig. 8). However, the gene was not significantly expressed in the macrophages.
Considering the differences in the expression of gene tags, we visualized the signaling pathways involving PLVAP+ ECs and found that they were closely related to each cell population (Fig. 4f). Interestingly, the THY1 signaling pathway only existed in PLVAP+ ECs and was not found in other cell populations, especially those closely related to TAM1 (Fig. 4g), which may indicate that there was a unique interaction between PLVAP+ ECs and TAM1. This may indicate that STC2 is involved in the reprogramming of endothelial cells and may regulate the reprogramming and function of macrophages through PLVAP+ECs.
The above results suggested that ERS plays a role in endothelial cells, thereby affecting the function of macrophages.
Establishment and validation of the radiogenomic signature
Although we already understand the important role of ERS in the occurrence and development of HCC and the TME, we hope to make early judgments on patients through non-invasive methods to better guide subsequent treatment. Utilizing resources from the TCIA databases, we embarked on a comprehensive analysis to discern specific radiological features that could predict ERSrisk. 1,766 original image features were extracted from the ROIs. After intragroup consistency analysis and intergroup consistency analysis, a total of 1,162 features whose ICC> 0.9 were retained. Thirty features were identified from the training cohort after Spearman’s correlation analysis and mRMR. The LASSO Cox regression algorithm further narrowed down a fusion signature that retained nine archetypal features (Fig. 5a and b). The resulting signature was formulated as follows (Fig. 5c):
radscore=
AP-wavelet-LHL-firstorder-Skewness*0.034045
-VP-wavelet-HHH-firstorder-TotalEnergy*0.06068
+VP-wavelet-HLL-glcm-Correlation*0.045428
-AP-wavelet-LLH-glcm-Idmn*0.038591
+AP-wavelet-HLL-gldm-DependenceVariance*0.032164
-AP-wavelet-LLH-glrlm-RunVariance*0.061446
-AP-wavelet-LHL-glcm-Imc2*0.012521
+VP-wavelet-LLL-firstorder-RootMeanSquared*0.052823
-VP-wavelet-LLH-glcm-Correlation*0.058891
The C-index of the radiogenomic signature was 0.647 in the TCG-LIHC cohort. The radiogenomic signature was significantly associated with OS in the TCGA-LIHC cohort (Fig. 5d).
Multi-scale validation of the radiogenomic signature
We used a surgical therapy cohort to validate the radiogenomic signature and obtained good results, with a C-index of 0.607. The prognosis of the RDS-low group was significantly better than that of the RDS-high group (Fig. 6a). Our study subsequently assessed the relationships between the radiogenomic signature and systemic combination therapy response(Fig. 6b). Interestingly, as shown in Fig. 6c and d, we found that patients in the RDS-low group (45%) had a greater objective response rate than did those in the RDS-high group (5%). In addition, the PD rate was greater in the RDS-high group than in the RDS-low group, which was also confirmed in the entire cohort (Table 2). Next, we explored the ability of the radiogenomic signature to predict systemic combination therapy response. We found that the signature had excellent performance in predicting the ORR, with an AUC of 0.785 and the DCR, with an AUC of 0.768 (Fig. 6e). Finally, Kaplan–Meier plots of progression-free survival (PFS) confirmed that patients in the RDS-low group had a longer PFS than those in the RDS-high group (Fig. 6f).