3.1 Identification and validation of TEXSRGs-related clinical phenotypes

The 342 HCC samples in the TCGA dataset were split into high and low categories based on the median value after the abundance of TEX was determined using ImmuCellA. As shown in Fig. 1A, 892 genes, including 447 upregulated genes and 491 downregulated genes, were found as TEX-related genes using cut-off parameters of p-value less than 0.05 and |logFC|≥ 0.5. These stemness-related genes were crossed with the TEX-related genes to create a total of 146 TEXSRGs, which were then used in the NMF clustering study. When k = 2, the categorization of TEXSRG-related subgroups was the most reliable (Fig. 1B, C). The PCA findings demonstrated that Cluster 2 patients could be easily differentiated from Cluster 1 patients (Fig. 1D). Furthermore, there were substantial variations in Grade, TNM stage, recurrence rate, and patient survival status between the two clinical categories (Fig. 1E). Compared to Cluster 1, HCC patients in Cluster 2 had superior mortality rates (Fig. 1F). Additionally, TEX and mRNAsi were more abundant in Cluster 1 HCC patients than in Cluster 2 patients (Fig. 1G). Genetic mutation research revealed that the two groups mutation rates were considerably different from one another (Fig. 1H). The sample with the highest mutation frequency in Cluster 1 was TP53, and the specific genes were AXIN1, BAP1, MUC4, and RYR2, while the sample with the highest mutation frequency in Cluster 2 was CTNNB1, and the specific genes were APOB, OBSCN, ABCA13, and LRP1B. Lastly, we discovered through GSVA analysis that Cluster 1 was primarily enriched in processes related to the cell cycle, whereas Cluster 2 was primarily enriched in processes related to metabolism (Fig. 1I).

Fig. 1

Identification of TEXSRGs-related clinical phenotypes in the TCGA cohort. A 892 genes were found as TEX-related genes. B, C When k = 2, the categorization of TEXSRG-related subgroups was the most reliable. D The PCA findings demonstrated that Cluster 2 patients could be easily differentiated from Cluster 1 patients. E There were substantial variations in Grade, TNM stage, recurrence rate, and patient survival status between the two clinical categories. F Compared to Cluster 1, HCC patients in Cluster 2 had superior mortality rates. G TEX and mRNAsi were more abundant in Cluster 1 HCC patients than in Cluster 2 patients. H Genetic mutation analysis. I GSVA analysis

Subsequently, we validated the TEXSRG-related clinical phenotypes in the ICGC dataset. Similar to the TCGA dataset, the HCC sample was divided into two subgroups (Fig. 2A) and patients in Cluster 2 could be easily distinguished from those in Cluster 1 (Fig. 2B). There were also significant differences between the two clinical categories in terms of TNM staging and patient survival (Fig. 2C). Mortality was higher in patients with HCC in Cluster 2 compared to Cluster 1 (Fig. 2D). In addition, TEX and mRNAsi were more abundant in HCC patients in cluster 1 than in cluster 2 (Fig. 2E). Finally, we found by GSVA analysis that Cluster 1 was mainly enriched in cell cycle-related processes, while Cluster 2 was mainly enriched in metabolism-related processes (Fig. 2F).

Fig. 2

Validation of TEXSRGs-related clinical phenotypes in the ICGC cohort. A The HCC sample was divided into two subgroups. B Patients in Cluster 2 could be easily distinguished from those in Cluster 1. C There were also significant differences between the two clinical categories in terms of TNM staging and patient survival. D Mortality was higher in patients with HCC in Cluster 2 compared to Cluster 1. E TEX and mRNAsi were more abundant in HCC patients in cluster 1 than in cluster 2. F GSVA analysis

3.2 TME heterogeneity among different clinical phenotypes

The results of the analysis according to the ESTIMATE algorithm showed that although there was no difference in stromal score and tumor purity between the two clinical subgroups, Cluster 1 had a significantly higher immune score than Cluster 2 (Fig. 3A). Subsequently, we further characterized their immunological profile in a variety of immune-related cell types using CIBERSORT, MCPcounter, TIMER, and xCELL. According to TIMER, samples in Cluster 1 had higher concentrations of B cells, CD4 T cells, neutrophils, macrophages, and myeloid dendritic cells than samples in Cluster 2, as shown in Fig. 3B. Figure 3C shows that samples in Cluster 1 in the MCPcounter database had higher abundance levels of T cells, CD8 T cells, B cells, monocytes, macrophage monocytes, neutrophils, endothelium cells, cancer-related fibroblasts, and myeloid dendritic cells than samples in Cluster 2. Comparing samples from Cluster 1 and Cluster 2, CIBERSORT found that Cluster 1 samples contained more naïve B cells, memory B cells, CD8 T cells, resting memory CD4 T cells, follicular helper T cells, Tregs, resting NK cells, activated NK cells, M0 macrophage, and resting myeloid dendritic cells (Fig. 3D). As determined by the xCELL database (Fig. 3E), samples in Cluster 1 had lower numbers of B cells, naïve CD4 T cells, CD8 T cells, class-switched memory B cells, common lymphoid progenitor, common myeloid progenitor, myeloid dendritic cells, M1 macrophage, memory B cell, mast cell, monocyte, NK T cells, Th1 CD4 T cell, Th2 CD4 T cell, Tregs, and activated myeloid dendritic cell when contrasted with those in Cluster 2. When the four aforementioned methods were merged, samples in Cluster 1 showed higher immune cell infiltrates in their TME. Considering the high level of immune cell infiltration but the poor prognosis in Cluster 1, we speculate that this may be related to TEX infiltration.

Fig. 3

TME heterogeneity among different clinical phenotypes. A The results of the analysis are according to the ESTIMATE algorithm. Immunological profiles in a variety of immune-related cell types were characterized using TIMER (B), MCPcounter (C), CIBERSORT (D), and xCELL (E). ns not significant; *p < 0.05; **p < 0.01; ***p < 0.001

3.3 The role of clinical phenotypes in predicting the efficacy of immunotherapy

Even though the impacts of dysfunction scores were reversed, TIDE analysis revealed that cluster 1 had substantially greater TIDE and exclusion scores than cluster 2 (Fig. 4A). Cluster 1 had a reduced “response” proportion when compared to the anticipated immunotherapy response rate (Fig. 4B). Last but not least, we discovered that patients in cluster 1 had reduced levels of PD-L1 and higher levels of CD276, CTLA4, CXCR4, IL1A, LAG3, TGFB1, TNFRSF4, TNFRSF9, and PD1 compared to patients in cluster 2 (Fig. 4C).

Fig. 4

The role of clinical phenotypes in predicting the efficacy of immunotherapy. A TIDE analysis. B Cluster 1 had a reduced “response” proportion when compared to the anticipated immunotherapy response rate. C patients in cluster 1 had reduced levels of PD-L1 and higher levels of CD276, CTLA4, CXCR4, IL1A, LAG3, TGFB1, TNFRSF4, TNFRSF9, and PD1 compared to patients in cluster 2. ns not significant; *p < 0.05; **p < 0.01; ***p < 0.001

3.4 Identification and validation of TEXSRGs-related prognostic model

To further quantify the prognostic value of TEXSRGs in HCC, we did a univariate Cox regression analysis of TEXSRGs combined with survival status (Fig. 5A), and prognostically relevant TEXSRGs were obtained to be included in the LASSO Cox analysis (Fig. 5B). Eight genes (ZIC2, ESR1, PAFAH1B3, TNNT1, CDCA7, HMGA2, MYRIP, and FCER1G) were included to construct a TEXSRGs-score that could accurately predict the prognosis of HCC according to the following formula: TEXSRGs-score = (0.1298 × ZIC2) + (0.0349 × TNNT1) + (0.0371 × CDCA7) + (0.00487 × HMGA2) + (0.0691 × PAFAH1B3)—(0.01806 × MYRIP) + (0.1883 × FCER1G)—(0.0395 × ESR1). According to the median TEXSRGs-score, we determined the TEXSRGs-score for each sample and divided them into two groups: high and low. Patients with high TEXSRGs-scores had considerably lower survival rates than patients with low TEXSRGs-scores, according to the Kaplan–Meier curve (Fig. 5C). The later grade, later TNM stage, recurrence, and survival status were associated with the high TEXSRGs-score (Fig. 5D). By measuring the AUC, the prognostic strength of the TEXSRG-related signature was assessed. According to the findings, the TEXSRGs-related signature performed well, with an AUC of 0.702, 0.682, and 0.692 for predicting the 1-, 2-, and 3-year survival of HCC patients, respectively (Fig. 5E). Additionally, TEX and mRNAsi were more abundant in patients with high TEXSRGs-score than those with low TEXSRGs-score (Fig. 5F). Furthermore, we validated the TEXSRGs-related prognostic signature in the ICGC cohort. We found that patients with high TEXSRGs-scores had considerably lower OS than patients with low TEXSRGs-scores (Fig. 5G) and the TEXSRGs-related signature performed well, with an AUC of 0.741, 0.693, and 0.709 for predicting the 1-, 2-, and 3-year survival of HCC patients, respectively (Fig. 5H). The TEXSRGs-related signature was an independent influence on prognosis in both datasets (Fig. 5I). According to the GSEA results, the differential genes between the two groups were mainly enriched in biological processes related to energy metabolism (Fig. 6). Lastly, we developed a nomogram to predict the survival rates for HCC patients by integrating the TEXSRGs-score and TNM stages (Fig. 7A), which performed well (Fig. 7B).

Fig. 5

Identification and validation of TEXSRGs-related prognostic model. A Univariate Cox regression analysis. B LASSO Cox analysis. C Patients with high TEXSRGs-scores had considerably lower survival rates than patients with low TEXSRGs-scores. D The later grade, later TNM stage, recurrence, and survival status were associated with the high TEXSRGs-score. E The prognostic strength of the TEXSRG-related signature was assessed. F TEX and mRNAsi were more abundant in patients with high TEXSRGs-score than those with low TEXSRGs-score. G Patients with high TEXSRGs-scores had considerably lower OS than patients with low TEXSRGs-scores. H TEXSRGs-related signature performed well, with an AUC of 0.741, 0.693, and 0.709 for predicting the 1-, 2-, and 3-year survival of HCC patients, respectively. I The TEXSRGs-related signature was an independent influence on prognosis in both datasets

Fig. 6

Identification of HALLMARK, GO, and KEGG enrichment between high- and low-TEXSRGs scores subgroups

Fig. 7

The predictive significance of the TEXSRGs-score was verified in the nomogram. A A nomogram model was built. B The likelihood of 1-, 2-, and 3-year survival rates predicted and actual overlap on the calibration curves showed good agreement

3.5 Role of TEXSRGs-score in predicting immunotherapy efficacy and target drug screening

Even though the impacts of dysfunction scores were reversed, TIDE analysis revealed that samples with high TEXSRGs-score had substantially greater TIDE and exclusion scores than those with low TEXSRGs-score (Fig. 8A). Samples with high TEXSRGs-score had a reduced “response” proportion when compared to the anticipated immunotherapy response rate (Fig. 8B). Moreover, patients with high TEXSRGs-score had reduced levels of PD-L1 and higher levels of CD276, CTLA4, CD4, CXCR4, IL1A, LAG3, TGFB1, TNFRSF4, TNFRSF9, and PD1 compared to patients with low TEXSRGs-score (Fig. 8C). Last but not least, we discovered five tumor-sensitive medications (Fig. 8D) and discovered that patients with higher TSRGs-scores had lower IC50 values for the medications kahalide F and ARV-825, indicating that they were more receptive to these medications (Fig. 8E).

Fig. 8

The role of TEXSRGs-score in predicting the efficacy of immunotherapy. A TIDE analysis. B Samples with high TEXSRGs-score had a reduced “response” proportion when compared to the anticipated immunotherapy response rate. C Difference analysis of ICIs genes. D The five tumor-sensitive medications. E Difference analysis of IC50 values. ns not significant; *p < 0.05; **p < 0.01; ***p < 0.001

3.6 Relationship between TEXSRGs-score and TEX and tumor stemness in clinical samples

According to the GEPIA database [22], increased expression of PAFAH1B3, ZIC2, and lower expression of ESR1 were found in HCC samples when compared to those in normal samples (Fig. 9A). We subsequently validated the expression levels of the eight TEXSRGs in clinical samples using qRT-PCR and found that the results were consistent with the GEPIA analysis (Fig. 10A). The expression of TEXSRGs in each sample was log-transformed, and the TEXSRGs-score of each sample was calculated by substituting the formula, and the clinical samples were divided into two groups of high and low TEXSRGs-score according to the median value. We then stained the tumor samples using IHC to detect the infiltration of CD8+ T cells in the tissues and found that the percentage of infiltration of CD8+ T cells was higher in the tumor tissues with high TEXSRGs-score than in those with low TEXSRGs-score (Fig. 10B). We also compared the differences in expression levels of TEX marker genes (HAVCR2, TIGIT, CTLA4, LAG3, and PD1) and the putative tumor stem cell marker genes (CD44 and PROM1) between high and low TEXSRGs-score groups in clinical tumor samples using qRT-PCR. We found that these genes were significantly upregulated in the high TEXSRGs-score group, suggesting that these samples may be in a state of high TEX infiltration and tumor stem cell activity (Fig. 10C). Furthermore, apart from CDCA7, FCER1G, and HMGA2, all genes were strongly associated with patient prognosis (Fig. 9B). Next, we focused on PAFAH1B3, ZIC2, and ESR1 because of their differential expression and prognostic potential. The human protein atlas (HPA) [23] database was performed to explore the three TEXSRGs’ protein expression in healthy and HCC tissues. Only PAFAH1B3 had differential expression among all of them (Fig. 9C). Lastly, we used immunohistochemistry to confirm its protein level and discovered that PAFAH1B3’s protein level was up-regulated in HCC tissues compared to normal tissues, which was compatible with the bioinformatics analysis discussed above (Fig. 10D).

Fig. 9

The expression of the eight TEXSRGs in normal and HCC tissues. A Expression levels and (B) prognostic values were explored in the GEPIA database. C Protein levels of ESR1 and PAFAH1B3 were explored in the HPA database. ns not significant; *p < 0.05

Fig. 10

Relationship between TEXSRGs-score and TEX and tumor stemness in clinical samples. A The expression levels of the eight TEXSRGs in clinical samples were validated using qRT-PCR. B The percentage of infiltration of CD8 T cells was higher in the tumor tissues with high TEXSRGs-score than in those with low TEXSRGs-score. C The differences in expression levels of TEX marker genes (HAVCR2, TIGIT, CTLA4, LAG3, and PD1) and the putative tumor stem cell marker genes (CD44 and PROM1) between high and low TEXSRGs-score groups in clinical tumor samples were compared using qRT-PCR. D PAFAH1B3’s protein level. ns not significant; *p < 0.05; **p < 0.01; ***p < 0.001