A pyrimidine metabolism-related signature for prognostic and immunotherapeutic response prediction in hepatocellular carcinoma by integrating analyses

Data acquisition

We retrieved the data of expression and clinical information of 374 HCC cases, along with 50 normal tissue samples from the TCGA data portal (https://portal.gdc.cancer.gov/projects/TCGA). The GEO database (http://www.ncbi.nlm.nih.gov/geo/) provided the microarray dataset GSE14520 (N=425). To evaluate the prognostic performance of PMRS, an independent validation cohort was utilized, consisting of 223 tumorous samples and 202 normal tissue samples obtained from the ICGC database (https://dcc.icgc.org/releases/current). For the prediction of immunotherapeutic response, we downloaded and included the IMvigor210 cohort in our study. Gene sets related to pyrimidine metabolism were searched for in the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/index.jsp), which were subsequently utilized in the investigation. The flow chart displayed the flow of this study.

Identification and enrichment analysis of the differentially expressed PMRGs

Differentially expressed pyrimidine metabolism-related genes between tumor and normal tissues were screened out with the criteria through the ‘limma’ R package [20]. The ‘clusterProfiler’ R package was utilized for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of those differentially expressed genes. (https://www.liebertpub.com/doi/10.1089/omi.2011.0118).

Consensus clustering of the PMRGs

Based on the differentially expressed PMRGs, the R package ‘ConsensusClusterPlus’ was utilized to execute consensus clustering and determine discrete clusters related to pyrimidine metabolism [21]. To ensure the robustness and stability of our results, we replicated this process 1000 times. Principal component analysis (PCA) was employed to classify PMRGs into clusters. According to the pertinent R package, the OS possibility of PMRG clusters was compared.

Construction and validation of the PMRS

We sequentially completed the Cox and Lasso regression to seek out the optimal independent predictive signature. Furthermore, we computed each sample’s risk score using the subsequent formula: Risk Score = (Coef₁ * Exp₁ +……+ Coef_n* Exp_n). Patients were separated into two categories based on the median value of the anticipated risk scores, which served as the threshold. The Kaplan-Meier analysis was employed to compare the probability of survival within these two categories. Time-dependent ROC analysis was performed through the R package “timeROC”. Survival curves of entire GEO and ICGC cohorts were used to examine the capacity of PMRS to distinguish between groups with different risks in terms of prognosis. Furthermore, stratified analysis was performed to evaluate the prognostic significance of PMRS within segments distinguished through clinical features.

Establishment and assessment of the pyrimidine metabolism-related nomogram

The independent prognostic factors were determined by the accomplishment of Cox regression upon PMRS along with its related clinical parameters. The nomogram was subsequently devised to calculate the survival probabilities of HCC patients. Utilizing calibration and time-dependent ROC, we assessed the reliability and quality of the nomogram. Additionally, employing DCA, we contrasted the cumulative effect of the entire nomogram with that of this model incorporating solely clinical measures.

Functional enrichment analysis

With the annotated gene set obtained from MSigDB, the gene set enrichment analysis (GSEA) was carried out to identify the differences in linked pathways across various risk categories [22].

Evaluation of immunotherapy efficacy and investigation of the immune profile

The ESTIMATE algorithm is a computational method used to estimate the presence and proportion of immune and stromal cells within tumor tissue samples [23]. Based on the theory of linear support vector regression, the method known as CIBERSORT evaluates 22 immune cell type fractions from mass tumor sample data regarding expression using an array of reference gene expression matrices [24, 25]. The comparison of the immune checkpoint, human leukocyte antigen (HLA), and immune function between distinct risk groups was consecutively conducted. The IPS is determined by four main gene categories that decide the degree of immunogenicity, and it is acquired through unbiased analysis utilizing machine learning techniques. In this work, we used IPS data to assess the prospective immunotherapy response. The Submap method was also used to forecast how various risk groups would respond to anti-PD-1 and anti-CTLA-4 treatment [23]. Furthermore, we also utilized clinical and immunotherapy data derived from the ‘IMvigor210CoreBiologies’ to assess the PMRS’s susceptibility to anti-PD-1 treatment in our analysis.

Identification of drug sensitivity

We downloaded the “DrugPredict.tsv” file and estimated the sensitivity score for common drugs and small chemical molecules in the HCC therapy. Additionally, we mined the correlation between PMRS genes and drug sensitivity using the CellMiner database [26] (https://discover.nci.nih.gov/cellminer/home.do).

Validation of the PMRS gene by single-cell RNA analysis

The Tumor Immune Single-cell Hub (TISCH) database (http://tisch.comp-genomics.org/home/) is an open-access platform that focuses on the single-cell sequencing data of the tumor microenvironment [27]. It provides comprehensive and detailed information to help researchers better understand the immune responses in tumors and corresponding single-cell transcriptome data [28]. The establishment of this database provides an important resource for further exploring tumor immunotherapy. We exploited this collection of data to figure out the association between the levels of expression of four genes responsible for establishing the PMRS and the infiltration of diverse immune cell types.

Cell culture and transfection

Hepatocellular carcinoma cell lines MHCC-97H and HepG2 were obtained from the American Type Culture Collection (USA). These cells were cultured at 37° C, 5% CO2 with 10% FBS and 1% penicillin/streptomycin in RPMI-1640. Lipofectamine 3000 was applied to transfect cells with siRNAs. The three siRNA sequences for DCK are as follows: (siRNA-1) 5’- GCCUGUCUCAGUCGAAUAATT – 3’, (siRNA-2) 5’- GCCUUGAAUUGGAUGGAAUTT – 3’ and (siRNA-3) 5’- GCUCAAAGAUGCAGAGAAATT – 3’.

Transfection efficiency

The efficiency of siRNA transfection was evaluated by Western blot assay. Cellular total proteins were isolated with RIPA lysate buffer and quantified with the bicinchoninic acid assay. Thirty micrograms of protein were loaded onto a polyacrylamide gel and separated by electrophoresis, then transferred to PVDF membranes. Following blocking with 5% BSA, the membranes were incubated overnight at 4° C with specific primary antibodies against DCK (1:1000, Proteintech, China) and GAPDH (1:1000, Proteintech). On the following day, the membranes were washed with TBST three times (10 min per time) and then incubated with a secondary antibody (1:3000) for 2 hours at room temperature. Signal intensities were quantified by the fully automatic chemiluminescence image analysis system.

CCK-8 assay

When their siRNA transfection level exceeded 90%, MHCC-97H and HepG2 cells were digested followed by inoculation onto 96-well culture plates with 2103 per well, and 3 duplicate wells were created for each group. The cells were then cultivated in a 5% CO2 incubator at 37° C and examined using the CCK-8 kit at various intervals of 0h, 24h, 48h, 72h, and 96h. (Beyotime Biotechnology, Shanghai, China).

Transwell assay

The lower chamber was covered with 700 μL complete medium. Then, 5×10⁴ cells of each group in 200 μL serum-free medium were seeded in the upper chamber, and then cultured in an incubator (37° C, 5% CO2) for 24h. Next, the cells underwent fixation with 4% paraformaldehyde for a period of 20 minutes, after which they were subjected to staining with crystal violet for 30 minutes. Use sterile cotton to wipe off the cells remaining in the upper chamber slightly. The transwells were dried at room temperature and photographed under a microscope.

Wound healing assay

After the siRNA-transfected cells in the 6-well plate grew to almost 100%, the cells were scratched in the plate with the 200 μL pipetting head and replaced with the serum-free medium to avoid cell growth effects. Subsequently, the image was captured with a microscope at 0h, 24h, and 48h after scratching.

Statistical analyses

All statistical analyses were performed using R software. Spearman analysis was used to analyze correlation coefficients. For normally distributed continuous data, we used independent sample tests, and for irregularly distributed continuous variables, Mann-Whitney U tests. The data comparison between the three groups was done using one-way ANOVA. For survival analysis, the Kaplan-Meier technique with a two-sided log-rank test was used. Statistical significance was defined as a p-value < 0.05. (*P<0.05, **P<0.01, ***P<0.001, ****P<0.0001).

Data availability statement

This article features the original research contributions presented by our study. For any additional inquiries, interested parties may contact the corresponding author directly.

Flow chart

The flow chart displayed the flow of this study (Figure 1).

Figure 1. The flow chart of this research.

Genetic and transcriptional landscape of PMRGs

From the MSigDB database, we initially collected and cataloged 92 genes associated with pyrimidine metabolism. Subsequently, 54 differentially expressed genes were filtered out by criteria of the absolute value of log₂(FC) > 1 and an adjusted P-value < 0.05. As displayed in the genetic variation analysis, missense mutations were the top common type among the 100 (23.2%) mutation samples out of a total of 431 samples (Figure 2A). RNA polymerase I subunit A (POLR1A), Carbamoyl-phosphate synthetase 2-aspartate transcarbamylase-dihydroorotase (CAD), and RNA polymerase III subunit A (POLR3A) ranked the top three mutant PMRGs. According to the results of differential analysis, there were 53 upregulated genes and 1 downregulated gene (Figure 2B). Additionally, we explored the CNV patterns of 54 PMRGs within HCC. (Figure 2C). Figure 2D depicts the locations of the CNV variations of 54 PMRGs on 23 chromosomes.

Figure 2. The genetic and transcriptional landscape of PMRGs. (A) The mutation profile of differentially expressed PMRGs in the TCGA-LIHC cohort. (B) The differential expression of PMRGs in HCC between tumor and normal tissues. (C) The CNV patterns of 54 PMRGs within HCC. (D) The locations of the CNV variations of 54 PMRGs on 23 chromosomes.

Development of consensus cluster on PMRGs

We used a consensus clustering approach to categorize the pyrimidine metabolism-related clusters to examine the expression properties and probable biological features of 54 PMRGs in HCC. According to the analysis of the cumulative distribution function, k=2 demonstrates exceptional clustering elasticity. The intra-cluster correlation is relatively strong while the inter-cluster correlation is low (Figure 3A–3C).

Figure 3. Identification of potential pyrimidine metabolism-related clusters in HCC patients. (A) The consensus clustering analysis of PMRGs. (B) Consensus CDF. (C) Delta area. (D) Kaplan-Meier survival analysis for different clusters. (E, F) The PCA plots of the PMRG clusters. (G, H) GO and KEGG enrichment analysis of differentially expressed PMRGs.

We employed the Kaplan-Meier survival analysis to show a more prolonged overall survival of HCC patients in cluster B (Figure 3D, P <0.001). The PCA graphs revealed a strikingly distinct pattern between clusters for all PMRGs and differentiated PMRGs. (Figure 3E, 3F). Then, functional enrichment analysis on the differentially expressed PMRGs was performed. As the visualization results demonstrated, biological processes (BPs) were predominantly enriched in the pyrimidine−containing compound metabolic process; cellular components (CCs) were significantly enriched in DNA and RNA polymerase complex; molecular functions (MFs) were mainly enriched in nucleotidyltransferase activity and catalytic activity. KEGG enrichment analysis exhibited that these genes were enriched in pyrimidine metabolism and nucleotide metabolism pathway (Figure 3G, 3H).

Development and verification of the prognostic model related to pyrimidine metabolism

We employed univariate Cox analysis to ascertain genes that exhibit a significant association with the survival of HCC patients (Figure 4A). Moreover, a pyrimidine metabolism-related 4-gene model was constructed to predict prognosis through Lasso and multivariate Cox regression analysis (Figure 4B, 4C). Risk score = (0.5961×POLR2L) + (0.3923×POLR3G) + (0.3057×DCK) +(0.3794×UCK2). On the basis of the cutoff, patients in the LIHC cohort were divided into two categories: high- and low-risk. To evaluate and validate the reliability and accuracy of the model, the TCGA cohort was separated into test and train groups. According to the survival curves, the patients in the high-risk group suffered from poorer prognosis in TCGA cohort and GSE14520 cohort (Figure 4D–4G). The area under curve (AUC) values for the predicted survival rates at 1-, 3-and 5-year were 0.797, 0.696, and 0.667 respectively in TCGA cohort; 0.901, 0.676, and 0.688 in TCGA test group; 0.768, 0.708 and 0.688 in TCGA train group; 0.604, 0.610, and 0.641 in GSE14520 cohort, which demonstrate the extraordinary prognostic predictability of the model (Figure 4H–4K). Figure 4L displays the risk score, clinical factors, and expression of four signature genes between the two risk categories. Furthermore, we developed external validation to evaluate the predictive power of the model utilizing the HCC cohort from the ICGC database. The ICGC cohort was divided into high- (n=116) and low-risk (n=116) groups. The KM survival curve (Figure 4M, P=0.011) revealed that the high-risk group had inferior outcomes than the low-risk group. In addition, the 1-year OS AUC value in the ICGC cohort was 0.703 (Figure 4N).

Figure 4. Identification and validation of the PMRS. (A) Forest plot of the univariate Cox regression. (B, C) Construction of PMRS. (D–G) The OS KM curves between different groups respectively in the whole TCGA cohort, test cohort, train cohort, and GSE14520. (H–K) The time-dependent ROC curves of the PMRS respectively in the whole TCGA cohort, test cohort, train cohort, and GSE14520. (L) The risk score and clinical event between the different groups in TCGA cohort. (M) The OS KM curves between high- and low-risk groups in the ICGC cohort. (N) The time-dependent ROC of the PMRS in the ICGC cohort.

Stratified analysis and development of the pyrimidine metabolism-related nomogram

We segregated individuals into different categories and calculated their overall survival for the stratified analysis to further validate the predictive value of the signature in subgroups with distinct clinical factors. The results of the distribution characteristics of the high- and low-risk categories remained comparable with the overall finding that the low-risk group obtained relatively favorable clinical results. (Figure 5A–5O). Subsequently, both major clinicopathological factors and the gene signature were analyzed employing univariate and multivariate Cox regression models to identify their prognostic significance. As the results displayed, Stage [(HR:1.661, 1.355−2.037, P < 0.001) in univariate Cox; (1.561, 1.264−1.928, P < 0.001) in multivariate Cox] and risk score [(1.318, 1.226−1.417, P < 0.001) in univariate Cox; (1.296, 1.198−1.401, P < 0.001) in multivariate Cox] were significantly associated with HCC prognosis and potential to be independent predictors (Figure 6A, 6B). We constructed a comprehensive nomogram with the identified independent prognostic factors to provide quantitative predictions of 1-, 3-, and 5-year overall survival (OS) probabilities in patients with HCC (Figure 6C). The AUC values were 0.802, 0.803, and 0.537 for risk, nomogram, and age respectively, which revealed the accuracy of the predictions (Figure 6D). The calibration curves indicated strong agreement between the actual and predicted overall survival (OS) rates (Figure 6E).

Figure 5. The distribution features and the OS KM curve of various clinicopathological factors in different risk groups. (A–C) Age, (D–F) Gender, (G–I) Grade, (J–L) Stage, (M–O) T.

Figure 6. Establishment and assessment of the nomogram. (A, B) Univariate and multivariate regression. (C) Nomogram for the prediction of HCC patients’ survival. (D) ROC curve of clinical factors and nomogram. (E) Calibration curves of the nomogram at 1-, 3-, and 5-years. (F) DCA curves of clinical factors and this nomogram.

To evaluate the clinical usefulness of the nomogram, we assessed its performance via decision curve analysis (DCA) [29]. Our results manifested that, in comparison to the model relying solely on clinical characteristics, the comprehensive nomogram could lead to increased net benefits and potentially improve the clinical management of the disease (Figure 6F).

Detection of the immune landscape

To further investigate the biological progress of the two risk groups by conducting GSEA. Our observations suggest that the high-risk group exhibited a significant enrichment of cancer-associated processes, such as cell cycle and ECM receptor interaction. In contrast, the low-risk group demonstrated significant enrichment of complement and coagulation cascades, which are known to function in the immune system (Figure 7A, 7B). The exploration of the immune microenvironment was performed by ESTIMATE and CIBERSORT algorithms. Following the outcomes of the ESTIMATE analysis, the group with the low risk had a significantly higher stromal score. The CIBERSORT analysis revealed that the high-risk group owned a significantly higher proportion of M0 macrophages and activated mast cells, while a higher proportion of plasma cells in the low-risk group (Figure 7C, 7D). These observational findings imply that the tumor growth and prognosis of the different groups may vary in terms of immune cell infiltration. For example, a higher proportion of M0 macrophages may affect antigen presentation and T cell activation, thereby accelerating tumor deterioration; and an increased number of plasma cells may indicate enhanced tumor immune response function, which can help stop the progression of tumors. Figure 7E, 7F displayed that 47 common immune checkpoint genes and 24 human leukocyte antigen (HLA) associated molecules were highly expressed in the high-risk group. B cells, mast cells, neutrophils, NK cells, type-II IFN response, and cytolytic activity were more stimulated in the low-risk group, whereas aDCs, IDCs, macrophage, Treg, and MHC class I were more stimulated in the high-risk group (Figure 7G). Tumor cells have developed various mechanisms to evade immune surveillance, including enhancing TMB, which can disguise the tumor cells from being recognized by the immune system. However, on the other hand, TMB has also been found to be a reliable predictor of the likelihood of an immunological response to the tumor, meaning that higher TMB levels may indicate a better likelihood of a successful immune attack against the tumor. The mutation landscape revealed that 159 (86.41%) of 184 samples in the high-risk group had mutations, with TP53 and CTNNB1 occupying the top two mutation frequencies, while 150 (84.75%) of 177 samples in the low-risk group did (Figure 7H, 7I). TMB was consistently considerably greater in the high-risk group compared to the low-risk group (Figure 7J).

Figure 7. Detection of the immune profile between different groups. The GSEA results for the high-risk (A) and low-risk (B) groups. (C) Difference between two risk groups’ stromal, immune, and ESTIMATE scores. (D) Comparisons of the fractions of immunocyte infiltration between distinct groups. (E) Differential checkpoint gene expression between two groups. (F) The differential HLA gene expression between two groups. (G) The levels of immune function in the two at-risk categories. The mutational status of the leading 20 genes in the high-risk category (H). The low-risk cohort (I). (J) The TMB concentration in the various categories. (*P<0.05, **P<0.01, ***P<0.001, ****P<0.0001).

Analysis of tumor immune dysfunction and exclusion (TIDE) includes the TIDE score, the Dysfunction score, the MSI score, and the Exclusion score. In our study, the TIDE score, Dysfunction score, and MSI score were higher in the low-risk group than in the high-risk group, whereas the Exclusion score was lower in the low-risk group than in the high-risk group, indicating the patient’s immune system in the low-risk group triggered a weaker reaction to the tumor, possibly due to the strong suppression of the immune cell (Figure 8A–8D). TIMER also investigated the relationship between immune cell infiltration and the expression levels of four model genes. Figure 8E–8H demonstrate that DCK and UCK2 were primarily associated with the infiltration of B cells, CD4+/CD8+ T cells, macrophages, neutrophils, and dendritic cells. In addition, we investigated the association between the level of infiltration of these immune cells and the mutational status of four model genes (Figure 8I–8L).

Figure 8. Further exploration of the immune landscape in PMRS. (A–D) TIDE analysis of the two risk groups. (E–H) The relation between the expression of four model genes and immune cell infiltration. (I–L) The relation between the mutation landscape of four model genes and immune cell infiltration. (*P<0.05, **P<0.01, ***P<0.001, ****P<0.0001).

Immunotherapy efficacy prediction

The immunotherapeutic response exhibited by the cohorts categorized as high-risk and low-risk was assessed using a combination of the IPS, Submap algorithms, and an external cohort undergoing analogous immunotherapy treatments. The fact that the percentage of ‘ips_ctla4_neg_pd1_neg’ is greater in the low-risk group than in the high-risk group suggests that patients in the low-risk group have a lesser expression of CTLA-4 and PD-1 immune checkpoint proteins in their tumor microenvironment. This often leads to reduced effectiveness of tumor immune surveillance. Therefore, patients in the low-risk group may respond better to immunotherapy and have a better prognosis compared to those in the high-risk group (Figure 9A–9D). Submap analysis yielded results uncovering the therapeutic response to anti-CTLA4 and anti-PD-1 immunotherapy among patients with HCC. The nominal P-value (P = 0.007) provided evidence supporting a greater likelihood of positive response to anti-PD-1 immunotherapy in the low-risk cohort relative to the high-risk group (Figure 9E). iMvigor210 cohort is designed to evaluate the efficacy and safety of atezolizumab, an immunotherapeutic drug used for the treatment of metastatic urothelial carcinoma (mUC). DCK and UCK2 in four model genes were highly expressed in the response group which referred to potentially serve as a predictive factor related to the response of tumor immune therapy (Figure 9F–9I). As the results in Figure 9J, 9K displayed, the Riskscore exhibits higher values in the response group and lower values in the no response group, as well as higher values in the SD/PD group and lower values in the CR/PR group. This observation commonly implies that the Riskscore can function as a stable and dependable indicator, thereby enhancing its capacity to effectively predict patient responsiveness to immunotherapy drugs.

Figure 9. Immunotherapy efficacy prediction. (A–D) IPS comparisons between the two risk categories. (E) Submap analysis between the different groups. (F–I) The differential expression of four model genes between response and non-response groups. (J) Comparisons of the risk score in different responses to anti-PD-1/L1 immunotherapy in IMvigor210 cohort. (K) Differences in risk score between the response and non-response groups in the IMvigor210 cohort.

Estimation of drug sensitivity

To further examine the clinical applicability of PMRS for the precise treatment of HCC, we evaluated the therapeutic effectiveness of commonly prescribed chemotherapeutic agents in various risk categories. Based on the findings, low-risk individuals were more sensitive to Afatinib, Cediranib, Dasatinib, GNE-317, Gefitinib, Ipatasertib, and Osimertinib, whereas less responsive to Axitinib, Gemcitabine, Selumetinib, Sorafenib, and Oxaliplatin (Figure 10A–10L). Figure 10M demonstrates the selection and analysis of drugs that are sensitive to the expression of four model genes, which were screened through the CellMiner database.

Figure 10. Drug sensitivity analysis. (A) Afatinib. (B) Axitinib. (C) Cediranib. (D) Dasatinib. (E) Gemcitabine. (F) GNE-317. (G) Gefitinib. (H) Ipatasertib. (I) Selumetinib. (J) Sorafenib. (K) Osimertinib. (L) Oxaliplatin. (M) The correlation between the drug sensitivity and expression levels of the four model genes.

Analysis of PMRS in single cell level

Tumor Immune Single Cell Hub (TISCH) is a massive managed database that incorporates nearly 2 million single-cell transcriptomic profiles. All the data were uniformly processed using a standardized methodology, which included quality control, eradication of batch effect, clustering, cell-type annotation, classification of malignant cells, variation in expression analysis, and functional enrichment analysis [28]. The scatter plots illustrated that DCK and POLR2L in the signature were significantly associated with tumor immune microenvironment through single cell analysis based on GSE140228, GSE98638, GSE146115, GSE146409, GSE166635, GSE179795 (Figure 11A–11F).

Figure 11. (A–F) Single-cell analysis of four model genes (DCK, POLR2L, POLR3G, UCK2) via TISCH based on GSE140228, GSE98638, GSE146115, GSE146409, GSE166635, GSE179795 datasets.

The role of key gene DCK in hepatocellular carcinoma was verified in vitro

In light of the previous findings, we conducted cell experiments to validate the use of PMRS. Firstly, the siRNAs targeting the expression of DCK were evaluated using Western blot assay. The results indicated that siRNA-3 exhibited the highest efficiency in downregulating DCK protein levels, and therefore was selected for further analysis conducted (Figure 12A). The activity of hepatocellular carcinoma cells was dramatically inhibited following DCK knockdown in MHCC-97H and HepG2 cell lines (Figure 12B, 12C). Afterward, the colony formation assays showed that the capacity of these two cell lines to generate colonies was significantly decreased after DCK knockdown (Figure 12D, 12E). The results of wound healing and transwell assays demonstrated a significant reduction in the migration and invasion capacity of hepatocellular carcinoma cells following DCK knockdown through transwell assay (Figure 12F, 12G) and wound healing assay (Figure 12H, 12I).

Figure 12. Suppression of DCK inhibited the proliferative and migrative ability of LIHC cells. (A) Transfected efficiency of three siRNAs targeting DCK. (B, C) The activity of MHCC-97H and HepG2 cell lines was dramatically inhibited following DCK knockdown. (D, E) Colony formation assay revealed that the capacity of these two cell lines to generate colonies was significantly decreased after DCK knockdown. (F–I) The migrative and invasive ability of MHCC-97H and HepG2 cell lines was considerably suppressed in the transwell assay and wound healing experiment. (*P<0.05, **P<0.01, ***P<0.001, ****P<0.0001).