RNA sequencing (RNA-seq) transcriptome data (count and transcript per million values) were obtained from The Cancer Genome Atlas (TCGA) database (TCGA-LIHC; http://portal.gdc.cancer.gov) and assessed using the R package ‘TCGAbiolinks’ (version 2.24.3). Data from 424 samples were downloaded. Similarly, the clinical information of 377 patients with HCC was obtained through TCGAbiolinks (BiocManager, version 1.30.18). The gene expression of patients with multiple locus samples was replaced by the average of the multiple samples. The genes detected more than twice had the expression value replaced by the average value. A total of 98 patients with virus-related HCC and the corresponding RNA-seq expression profiles were selected according to the following inclusion criteria: i) HBV, HCV or HBV + HVC infection; ii) no risk factors for alcohol, tobacco, non-virus induced cirrhosis and fatty liver disease; and iii) complete clinicopathological information including age, sex, tumor-node-metastasis (TNM) stage (20), survival time and status.

The validation set was obtained from the International Cancer Genome Consortium (ICGC) database (ICGC-LIRI-JP; http://dcc.icgc.org/). The RNA-seq transcriptome data and clinical information of the ICGC-LIRI-JP cohort were downloaded and processed in the same manner as the TCGA data. According to the aforementioned inclusion criteria, 53 patients with virus-related HCC and the corresponding RNA expression profiles were chosen. In addition, 174 patients with HCC with virus infection (HBV, HVC or HBV + HVC) were selected for extended validation of prognosis. The clinical data of the aforementioned patients with virus-related HCC are presented in Tables SI and SII.

The specific genes of platelets have been reported in previous experimental studies (21,22). The genes specifically expressed among the top 50 genes in the study by Gnatenko et al (21) and genes with a relative expression index >2 in the study by Raghavachari et al (22) were combined as the final platelet signature in the present study. The PSS of each sample was calculated using the R package ‘GSVA’ (version 1.44.5). The patients were then divided into high and low PSS groups according to the cut-off value (0.73) obtained using X-tile software (version 3.6.1; http://x-tile.software.informer.com).

Platelet-related function and hallmark signatures were obtained from the Molecular Signature Database (M15669; http://www.gsea-msigdb.org/gsea/msigdb/index.jsp) of the Gene Set Enrichment Analysis database (http://www.gsea-msigdb.org/gsea/index.jsp). Subsequently, gene set variation analysis (GSVA) was used to calculate the score of each sample in platelet-related function and hallmark pathways. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) database (https://david.ncifcrf.gov/), and results with P<0.05 were considered statistically significant.

The stromal, immune and tumor purity scores of each sample were obtained using the R package ‘ESTIMATE’ (version 1.0.13) (23). CIBERSORT (24) and xCell (25) were used to assess the extent of immune cell infiltration. The immune checkpoints were selected based on a previous study (26).

DEGs (foldchange >2; P<0.05) between the high and low PSS groups were identified using the R packages ‘DESeq2’ (version 1.36.0) (27) and ‘edgeR’ (version 3.38.4) (28). Mutual DEGs identified were used as the final DEGs for subsequent analysis. The PPI network was assessed using the Protein-Protein Interaction Networks Functional Enrichment Analysis (STRING) database (https://string-db.org/) and visualized by Cytoscape software (version 3.9.1; http://cytoscape.org/).

First, univariate Cox analysis of the aforementioned DEGs was used to select the prognostic genes with a significance of P<0.05. Subsequently, least absolute shrinkage and selection operator (LASSO) regression analysis was used to identify the final gene signature according to the minimum λ value (29). The RS of patients with virus-related HCC was evaluated according to the following formula: Risk score=∑ coef (gene) × expr (gene). Patients were then divided into high and low RS groups according to the median value. Based on the RS, Kaplan-Meier analysis was performed and time-dependent receiver operating characteristic (ROC) curves were generated to assess whether the RS could accurately predict the prognosis of patients with virus-related HCC. Furthermore, univariate and multivariate Cox analyses were performed, and a nomogram, calibration curve and ROC curve were generated to evaluate whether the RS can be used as an independent clinical prognostic predictor. Expression and survival of genes for model in HCC was obtained through the University of ALabama at Birmingham CANcer data analysis Portal (UALCAN) database (https://ualcan.pa-th.uab.edu/) and the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/).

The half-maximal inhibitory concentration (IC₅₀) values of chemotherapy and targeted drugs for patients with virus-related HCC were assessed using the R package ‘pRRophetic’ (version 0.5). The Tumor Immune Dysfunction and Exclusion (TIDE; http://tide.dfci.harvard.edu/) algorithm (30) was used to predict the anti-programmed cell death protein 1 and anti-cytotoxic T-lymphocyte associated protein 4 (CTLA4) response of each patient with virus-related HCC. Patients whose TIDE scores were <0 were considered non-responders, whilst patients whose TIDE scores were >0 were considered responders. Moreover, microsatellite instability (MSI), dysfunction and exclusion of T cells were also obtained through the TIDE algorithm.

A total of five paired HBV-related HCC and adjacent tissues from HCC patients who received hepatectomy from May 2019 to June 2019 were collected at the Affiliated Hospital of North Sichuan Medical College (Nanchong, China). The use of samples was approved by Research Ethics Committee of Affiliated Hospital of North Sichuan Medical College (Nanchong, China; approval no. 2023059). Total RNA was extracted using the total RNA isolation kit (cat. no. RE-03011; FOREGENE; Chengdu Fuji Biotechnology Co., Ltd.) from five paired HBV-related HCC and adjacent tissues following the manufacturer's instructions. The primescript^™ RT reagent kit (cat. no. RR037A; Takara Biotechnology Co., Ltd.) was used for reverse transcription. The reverse transcription volume was 20 µl and the temperature protocol was 37°C for 15 min and 85°C for 5 sec. Subsequently, quantitative PCR was performed using the TB Green^® Premix Ex Taq™ II (cat. no. RR820A; Takara Biotechnology Co., Ltd.). The reaction volume was 20 µl and the thermocycling conditions were 95°C for 30 sec, 95°C for 5 sec and 60°C for 30 sec (40 repeats), and a default melt curve procedure. Primer sequences were synthesized by Beijing Tsingke Biotech Co., Ltd. The specific sequences used are as follows: Cyclin-J-like protein (CCNJL), (forward) 5′-CCTGCGCGAGAAGGAACTG-3′ and (reverse) 5′-CGTTGTAGCGATCCATGAAGTG-3′; nuclear receptor subfamily 0 group B member 1 (NR0B1), (forward) 5′-CTCACTAGCTCAAAGCAAACGC-3′ and (reverse) 5′-GCGCTTGATTTGTGCTCGT-3′; tripartite motif containing 54 (TRIM54), (forward) 5′-AGGAGGTGTGCCAGACTATC-3′ and (reverse) 5′-GGTCGCCATACTGACGGATG-3′; GAPDH, (forward) 5′-GGAGCGAGATCCCTCCAAAAT-3′ and (reverse) 5′-GGCTGTTGTCATACTTCTCATGG-3′. RNA expression levels were normalized using GAPDH (31).

Statistical analyses were performed using R software (version 4.2.1; The R Foundation). Comparisons and pairwise comparisons between two groups were assessed using the Wilcoxon rank-sum test or the Wilcoxon signed-rank test, respectively. Comparisons between ≥3 groups were evaluated using the Kruskal-Wallis test, followed by Dunn's or Steel-Dwass correction.

Comparison of clinical features were assessed using the χ² or Fisher's exact test. Kaplan-Meier analysis was used to assess survival using the ‘survival’ package in R (version 3.4.0) and the cutoff values were obtained using X-tile software (version 3.6.1; http://x-tile.software.informer.com) or the median. DEGs (foldchange >2; P<0.05) between the high and low PSS groups were identified using the ‘DESeq2’ (version 1.36.0) and ‘edgeR’ (version 3.38.4) packages obtained by BiocManager (version 1.30.18). Mutual DEGs in DESeq2 and edgeR were used as the final DEGs for subsequent analysis. Correlation analysis was performed using Spearman's correlation analysis. P<0.05 was considered to indicate a statistically significant difference.

A schematic of the study design is presented in Fig. 1A. First, the transcriptome and clinical profiles of 98 patients with virus-related HCC were obtained from TCGA, as well as the gene set of the platelet signature from previous studies (21,22). A PSS was then calculated for each sample using GSVA. Based on the PSS, X-tile was used to divide patients into high and low PSS groups (Fig. 1B). Kaplan-Meier survival analysis indicated that patients in the high PSS group had a worse prognosis (P=0.026) than patients in the low PSS group (Fig. 1C), indicating that PSS is associated with patient prognosis in virus-related HCC. Notably, no significant difference in platelet counts was demonstrated between the high and low PSS groups (Fig. S1A).

Previous studies reported that platelet count and activation are associated with poor prognosis in patients with HCC (13,14). Furthermore, the accumulation of platelets in the TME promotes HCC progression by accelerating inflammation (16). Therefore, the present study assessed the difference in platelet activity and inflammation in the TME between the high and low PSS groups using GSVA. Notably, it was demonstrated that the main functions of platelets (including activation, aggregation and adhesion) were significantly enriched (P<0.05) in the high PSS group; however, there was no significant difference in the platelet-derived growth factor receptor α or β signaling pathway between the high and low PSS groups (Fig. 1D). Furthermore, GSVA of hallmark pathways indicated that inflammation-related pathways (including IL6-JAK-STAT3, TNF-α via NFκB, inflammatory response and IL2-STAT5 signaling pathways) were significantly enriched (P<0.05) in the high PSS group (Fig. 1E). The P53 and KRAS_UP signaling pathways were also significantly enriched (P<0.05) in the high PSS group, but there was no significant difference in other protumor pathways (including the Wnt-β-catenin, PI3K-AKT-MTOR and MTORC1 signaling pathways) between the high and low PSS groups (Fig. S1E). Subsequently, the expression of several inflammatory factors were compared between the high and low PSS groups. The expression levels of the proinflammatory factors IL1B and IL6 were significantly higher (P=0.04 and P=0.0073, respectively) in the high PSS group than in the low PSS group (Fig. 1F and G). However, the expression levels of IL10, TGFB1 and IFNG were not significantly different between the high and low PSS groups (Fig. S1B-D).

Inflammation, as a double-edged sword in tumor progression, recruits a large number of immune cells to accumulate in the TME to form an inflammatory immune microenvironment (32). Given the high association between inflammation and immunity, the present study assessed the immune infiltration and expression of immune checkpoints in the TME between the high and low PSS groups. First, the ESTIMATE algorithm was used and the results revealed that the high PSS group had significantly higher stromal (P=0.0057) and immune (P=0.024) scores than the low PSS group (Fig. 2A and B). Moreover, the low PSS group had a significantly higher tumor purity (P=0.0056) than the high PSS group (Fig. 2C). To further assess these results, the present study used two algorithms (xCell and CIBERSORT) to evaluate the infiltration of immune cells in the TME in the high and low PSS groups. The xCell results demonstrated that the compositions of activated myeloid dendritic cells, macrophages (including M1 and M2) and neutrophils were significantly higher (P<0.05) in the high PSS group than in the low PSS group. Meanwhile, the high PSS group had significantly lower proportions (P<0.05) of CD4+ central memory T cells, naïve CD8+ T cells, eosinophils, mast cells and B-cell plasma than the low PSS group (Fig. 2D). Similarly, the CIBERSORT results revealed that the high PSS group had significantly higher levels of infiltration (P<0.05) of M2 macrophages and neutrophils than the low PSS group (Fig. S2A). These results indicate that the high PSS group may have infiltration of more immunosuppressive cells than the low PSS group, which contributes to the formation of an immunosuppressive microenvironment. Furthermore, as the expression of immune checkpoints is another crucial factor in the formation of an immunosuppressive microenvironment (33), the present study also compared the mRNA expression levels of 11 immune checkpoints between the high and low PSS groups. The expression levels of Hepatitis A virus cellular receptor 2 (HAVCR2; TIM3), CD80 and CD86 were significantly higher (P<0.05) in the high PSS group than the low PSS group, whilst the expression of other immune checkpoints did not significantly differ between the two groups (Fig. 2E). The aforementioned results suggest that the high PSS group may have a stronger immunosuppressive microenvironment than the low PSS group.

Intrinsic molecular differences are known to determine the hallmarks of malignant biological behavior in tumor cells (34). Based on the aforementioned results, the present study compared the transcriptome differences between the high and low PSS groups to assess the underlying key molecules and mechanisms. DESeq2 and edgeR (two differentially expressed analysis R packages for sequencing data) were used to identify 137 mutual DEGs between the high and low PSS groups (Fig. 3A). Subsequently, functional enrichment analysis of the 137 DEGs was performed using DAVID to evaluate their functional characteristics. The biological processes results revealed that these DEGs were enriched in cell-cell signaling, chemokine-mediated signaling pathways and inflammatory responses. In molecular function, these DEGs were also enriched in structural molecule activity, receptor binding, CXCR chemokine receptor binding and chemokine activity. Furthermore, the cellular component demonstrated that these DEGs were mainly enriched in extracellular components (Fig. 3B). The results of GO analysis indicated that these DEGs were fully involved in cellular interactions (such as tumor and immune cells) and inflammatory immune responses in the TME. KEGG analysis then revealed that the 137 DEGs were mainly enriched in mineral absorption, Staphylococcus aureus infection, pancreatic secretion, the estrogen signaling pathway, the IL-17 signaling pathway and viral protein interactions with cytokines and cytokine receptors (Fig. 3C). Notably, the IL-17 signaling pathway and viral protein interactions with cytokines and cytokine receptors also suggest that these DEGs might serve crucial roles in inflammatory immune responses and virus-related tumoral regulation. In addition, a PPI network was constructed using STRING to show the interactions of 137 DEGs (Fig. 3D).

Tumoral prognostic models are currently widely used in several clinical and basic studies to identify novel molecular subtypes and aid in clinical prognostic assessment (26). After identifying 137 DEGs, the present study performed univariate Cox analysis for these DEGs. The genes with P<0.05 were then chosen for LASSO regression analysis (Fig. 4A and B), and 3 genes (CCNJL, NR0B1 and TRIM54) were identified to construct the prognostic model. Based on the regression coefficients of the LASSO regression analysis and the expression values of the 3 genes in every sample, an RS was calculated for every patient with virus-related HCC and the median of the RS was used as the cutoff value to divide patients into high and low RS groups (Fig. 4C).

The 3 genes (CCNJL, NR0B1 and TRIM54) were mainly highly expressed in the high RS group (Fig. 4C). Kaplan-Meier survival analysis indicated that the patients with a high RS had a worse prognosis (P=0.0079) than the patients with a low RS (Fig. 4D). Subsequently, a time-dependent ROC curve was used to evaluate the prognostic model and it was demonstrated that the area under the curve (AUC) values at 1, 3 and 5 years were 0.80, 0.74 and 0.74, respectively (Fig. 4E). To assess the robustness of the prognostic model, 53 patients with virus-related HCC were chose from the ICGC database. Patients with death or high expression values of the 3 genes were mainly concentrated in the high RS group (Fig. 4F). Kaplan-Meier survival analysis also indicated that the patients with a high RS had a worse prognosis (P=0.0039) than the patients with a low RS (Fig. 4G). Furthermore, the AUC values at 1, 2 and 3 years were 0.83, 0.76 and 0.79, respectively (Fig. 4H). Survival analysis was also performed for 174 patients with HCC with virus infection and other risk factors, and consistent results were obtained (Fig. S3A).

Additionally, the mRNA expression levels of the 3 genes were assessed in 5 paired HBV-related HCC tissues (Fig. 4I), and a Sankey diagram was constructed to show the distribution of patients between the RS and PSS groups (Fig. S3B). Meanwhile, data from the UALCAN database indicated that all 3 genes (CCNJL, NR0B1 and TRIM54) were significantly highly expressed in tumors, and, in comparison with low/medium expression, the high expression of these genes was significantly associated with poor prognosis in patients with HCC (Fig. S3C-H). Considering that genes are not exactly consistent at the mRNA and protein expression levels (35), the present study used the HPA database and revealed that the staining intensity of CCNJL and TRIM54 was notably stronger in tumor tissues than in normal tissues (Fig. S3I and J). Unfortunately, immunohistochemical staining results for NR0B1 could not be obtained from the HPA database.

After validating the robustness of the model in the present study, the clinical value of the RS was evaluated. First, univariate and multivariate Cox analyses of several clinical factors (including age, sex and TNM staging) and RS were performed. The results of both analyses demonstrated that RS could be an independent prognostic factor for patients with virus-related HCC (Fig. 5A and C). Subsequently, a nomogram was constructed to predict the 1-, 3- and 5-year survival rates of the patients (Fig. 5E). The total score was determined based on individual factor scores and contributed to the survival prediction. Furthermore, the calibration curves were used to describe the accuracy of the model: The calibration curves revealed the highest degree of agreement between the predicted and observed survival probabilities at 3 years in the training set, in comparison with the validation set (Fig. 5G). Notably, the 3-year ROC curves of age, sex, TNM staging and RS demonstrated that RS had the highest AUC value (AUC=0.74; Fig. 5H). Subsequently, the results of Cox analysis from the validation set suggested that RS was a potential independent prognostic factor in virus-related HCC (Fig. 5B and D). The calibration curve also revealed similar results to the training set, and the RS acquired the highest AUC value (AUC=0.79) in the 3-year ROC curves (Fig. 5F, I and J). These results indicate that RS is a potential clinical prognostic indicator in virus-related HCC.

Increasing evidence has demonstrated that patients with HCC benefit from systemic therapy (such as chemotherapy, targeted therapy and immunotherapy), although certain patients are resistant to systemic therapy (36,37). However, there is insufficient evidence to assist in judging the efficacy of systemic therapy in patients with virus-related HCC. The aforementioned results of the present study indicate the importance of RS in prognostic evaluation. Therefore, the potential of RS in predicting the response to systematic treatments was subsequently assessed. First, the IC₅₀ of several chemical or targeted drugs indicated that the patients in the high RS group were more likely to be sensitive to chemotherapy and targeted therapy than those in the low RS group. Among these drugs, the IC₅₀ values of AKT inhibitor III, gemcitabine, paclitaxel, rapamycin, sorafenib and sunitinib were significantly lower (P<0.05) in the high RS group than in the low RS group (Fig. 6A-H). The TIDE, MSI, dysfunction and exclusion scores were then obtained for each sample using the TIDE database. Notably, it was demonstrated that the high RS group had significantly higher TIDE (P=0.00035) and exclusion (P<0.0001) scores, as well as significantly lower MSI scores (P=0.016) than the low RS group, whilst there was no significant difference in the dysfunction score between the high and low RS groups (Fig. 6I-L).

Furthermore, it was demonstrated that the low RS group had a significantly higher response rate to immunotherapy than the high RS group (P=0.0002; Fig. 6M). To further support this finding, a correlation analysis was performed between the RS and several immune checkpoints. The results revealed that RS was significantly positively correlated (P<0.05) with the expression of CD80, CD86, CTLA4, HAVCR2, galectin 9, TNFRS4 and TNFRS9 (Fig. 6N-T), although there was no significant correlation between RS and the expression of programmed cell death protein 1, CD274, lymphocyte-activation gene 3 and T cell immunoreceptor with Ig and ITIM domains (Fig. S4). Therefore, the results suggest that patients in the high RS group were better suited for chemotherapy or targeted therapy, whilst patients in the low RS group were better suited for immunotherapy.

Given the excellent performance of PSS in predicting platelet activity and inflammation, RS was subsequently combined with PSS to assess the performance in predicting platelet activity, inflammation and pro-tumor pathways. The high RS combined with high PSS group had significantly higher scores for platelet activation, platelet aggregation, platelet-derived growth factor receptor β signaling pathway and platelet activation signaling and aggregation (P<0.05) than the low RS combined with low PSS group (Fig. 7A). These results indicated that the high RS combined with high PSS group may have the highest platelet activity among the four groups. Compared with the low RS combined with low PSS group, it was demonstrated that inflammatory-related pathways (including IL6-JAK-STAT3, TNFA via NFκB, inflammatory response and the IL2-STAT5 signaling pathways) were significantly enriched (P<0.05) in the high RS combined with high PSS group, which was consistent with previous results based on the PSS. Notable, it was also demonstrated that tumor-related pathways (including the PI3K-AKT-MTOR, MTORC1 and P53 signaling pathways and genes up-regulated by KRAS activation) were significantly enriched (P<0.05) in the high RS combined with high PSS group, which was not demonstrated in the separate PSS groupings (Figs. 7B-D and S5A-I). Notably, mutated P53 exerts oncogenic effects through other tumor-related pathways and activates downstream feedback to inhibit the death of tumor cells (38).

Furthermore, it was demonstrated that the high RS combined with high PSS group had a significantly higher (P<0.05) expression of inflammatory factors including IL1B, IFNG, TGFB1 and IL10 than the low RS combined with low PSS group (Figs. 7E and F and S5J-L). Additionally, Kaplan-Meier survival analysis indicated that the high RS combined with high PSS group had the worst prognosis (P=0.002) among the four groups (Fig. 7G). In summary, these results indicate that RS combined with PSS could predict platelet activity, inflammation and protumor pathways in the TME.