RNA sequencing data and corresponding clinical information of 179 virus-relevant patients with liver cancer were downloaded from TCGA website (https://portal.gdc.cancer.gov/repository). RNA sequencing data and clinical information of another 260 virus-related HCC samples were obtained from ICGA website (https://dcc.icgc.org/projects/LIRI-JP). Besides, patients from GSE84337 (n=75) in the GEO repository was screened to acquire clinical parameters and normalized gene expression data. Clinical information of virus-related liver cancer patients was shown in Table S2 . Samples lacking significant clinicopathological or survival information were excluded from further analysis.

200 inflammatory response-related genes were found in the Molecular Signatures database and listed in the Supplementary Table 1 . Furthermore, t-distributed Stochastic Neighbor Embedding (t-SNE), a nonparametric and unsupervised algorithm, was employed to sort or condense patients into diverse clusters, based on given signatures or hallmarks by using an R package Seurat (28). According to the OS data, two groups were singled out for comparison to determine the “inflammation^high” and “inflammation^low” clusters. The limma algorithm was applied to filtrate DEGs between the above two groups, generating genes with false discovery rate (FDR) adjusted p-value<0.05 and absolute value of log2 (fold change)>1 were regarded as inflammation-related DEGs.

The STRING database (https://string-db.org/) was used to establish the protein–protein interaction (PPI) network among sufferers with co-expression coefficients >0.4. Also, cytoscape software (version 3.7.2) was exploited to visualize the network. Moreover, the hub genes were screened with the MCC algorithm of the cytoHubba plugin. The correlation between the expression of inflammation-related genes was identified by the “reshape2” R package.

To explore the potential mechanisms and pathways about inflammation-related genes, the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis were conducted among IRGs using the R packages “clusterProfiler,” “enrichplot,” “ggplot2,” and “org.Hs.eg.db.”

Based on the expression of inflammation-related genes (IRGs), we classified distinct inflammation-regulated groups through consensus clustering with the k-means method. The number of patterns and corresponding homologous stability were defined by consensus clustering algorithm using the R package ConsensusClusterPlus with 1,000 repetitions (29).

The immune infiltration characteristics (the immune and stromal scores) of virus-related HCC, based on the RNA-seq dataset of TCGA database LIHC, were evaluated by ESTIMATE algorithm (30). Then, CIBERSORTx was applied to quantify the percentages of 22 immune cell subtypes of each patient in the TME (31). Also, the correlation between the subsets on PD-1, PD-L1, and CTLA-4 expression was assessed.

To define preliminary inflammation-related DEGs that were significantly associated with OS in the training cohort, univariate Cox regression analyses using the R package “survival” were further implemented among favorable and risk DEGs, of which p<0.05 were regarded as positive. Also, Least Absolute Shrinkage and Selection Operator (LASSO) regression with 10-fold cross-validation was explored to narrow down the prognosis-related DEGs applying the R package “glmnet” (32). Meanwhile, based on the Gaussian finite mixture model (GMM), classification was conducted with model-based hierarchical agglomerative clustering with the R package “mclust” (33). Afterwards, the clusters made up of DEGs were classified by GMM and logistic regression analysis was utilized to construct combined models to predict the OS status for patients. Besides, to calculate the predictive value of models, receiver operating characteristic (ROC) curves were established by assessing the area under curves (AUCs). Subsequently, the risk scores of patients were estimated according to the expression level of each inflammatory response-related gene and its relevant regression coefficient. The formula was established as follows: risk score = ∑iCoefficient (mRNAi) × Expression(mRNAi). On the basis of the median risk scores, patients were divided into high- and low-risk clusters among training and validation cohorts. The Kaplan-Meier analysis was applied to compare the OS between the high- and low-risk groups. The predictive value of the prognostic model was assessed on account of ROC analysis. The principal component analysis (PCA), acquiring a low-dimensional cluster distribution from high-dimensional gene sets, was utilized for validating the sectionalization results.

The correlation between IRG score and clinical factors was explored. Univariate Cox and multivariate Cox regression analysis were firstly implemented to prove whether IRG score was an independent prognostic predictor. Ulteriorly, a grouping analysis was conducted to explore whether the IRG score sustained its predictive reliable in disparate subgroups on the basis of multifarious clinical variables. Furthermore, the infiltrating levels of immune cells and immune checkpoint (ICP) were analyzed between the distinct risk subgroups and the relevance between IRG score and tumor mutation burden (TMB) score, microsatellite instability (MSI) score, and cancer stem cells (CSC) score was examined.

Nomogram was constructed by the rms R package. Calibration curves and decision curve analysis (DCA) were utilized to quantify the consistency between the predicted and the observed results for 3-, and 5-years survival rates (34).

On the basis of the cBioPortal database, genetic alteration data was acquired. And the number and quality of mutations between the two IRG clusters were analyzed using the R “Maftools” package (35). Subsequently, the online database TIDE (Tumor Immune Dysfunction and Exclusion, http://TIDE.dfci.harvard.edu/) and immunophenotype score (IPS) were calculated to execute immune checkpoint inhibitor response of each virus-related HCC patient in the two groups to assess the value of the IRG in terms of prognostic immunotherapy response.

The pRRophetic R package was used to predict the half-maximal inhibitory concentration (IC50) value of cancer drugs in diverse risk subgroups, which represented the availability of a substance in inhibiting a particular biological or biochemical process (36, 37).

The samples contained 58 blood samples from virus-related HCC patients from West China Guangan Hospital, Sichuan University, between March and November in 2019. The diagnoses of HCC were confirmed by senior pathologist. None of the patients experienced radiotherapy or chemotherapy treatment before samples collection. Also, 50 blood samples from healthy people were considered as the control cluster. Informed consent was obtained from all participants for the use of their blood samples in this study. This project was approved by the Clinical Research Ethics Committee of Chengdu medical college.

The human cell lines (WLR68, LO2, Huh-7, SMMC7721, HepG2, and HCCLM3) were obtained from the School of Bioscience and Technology, Chengdu Medical College (ChengDu, China). All of them were cultured in DMEM (Gibco) medium, which were supplemented with 10% fetal bovine serum (FBS) at 37 °C with 5% CO₂. In addition, the cells were photographed after treatment with paraformaldehyde.

Approximately 8 ml of whole blood from participants was gathered in EDTA tube. After centrifuged at 1,2000g at 4°C to spin down the blood cells for 10 min, the supernatant was shifted into microcentrifuge tubes. Afterwards, plasma was aliquoted or stored at −80°C. RNA was isolated from 400 μL plasma with the mirVana PARIS kit (Ambion, USA) abided by the manufacturer’s protocol. The PrimeScript™ RT reagent kit (TaKaRa) was further applied for reverse transcriptase reaction. Reverse transcription−quantitative PCR (RTqPCR) were implemented to attest the expression levels of the six hub genes in plasma samples and cell lines. The mRNA expression level of MEP1A, CCL20, ADORA2B, TNFSF9, ICAM4, and SLC7A2, was normalized by GAPDH. Fold differences were calculated for each group with normalized CT values.

Full-length SLC7A2 cDNA was synthesized and cloned into the pCS-CG vector (Addgene, Cambridge, MA, USA). shRNA sequences specifically against CCL20 (shCCL20) and control-shRNA against luciferase (shCtrl) were expressed from pLKO.1-puro (Addgene, Cambridge, MA, USA). Production of lentiviral particles using HEK-293T cells and subsequent infection of Huh7 and SMMC7721 cells were performed according to the manual instructions.

14 male nude mice (5-week-old) were purchased from Beijing Vital River Laboratory Animal Technology Co., Ltd. (Beijing, China). SMMC7721 (1×10⁷) stably transfected with shCtrl or shCCL20/overSLC7A2 was subcutaneously injected into the right gluteal region of each nide mice (n=7). After tumor formed, the tumor volume was calculated every 3 days on the basis of the formula: volume(mm³)=width2 (mm²) *length (mm)/2. All the mice were euthanized and the formed tumors were weighted after 30 days. The animal experiment was approved by the Animal Care Committee of Chengdu Medical College.

All analyses were completed on the strength of R language (Version 4.2.1). Student’s t-test, chi-squared test, or Wilcoxon test was applied to compare the differences between groups. Spearman’s correlation test was performed to evaluate the correlation between variables. p-value of <0.05 was deemed as statistically positive.

The study flowchart is revealed in Figure 1 , which precise procedure is as follows: First, RNA sequencing from the TCGA database for 179 virus-related HCC sufferers was obtained, as well as 200 IRGs from the Molecular Signatures database. t-SNE, was applied to sort or condense patients into diverse clusters, based on 200 IRGs and IRG DEGs were identified from survival analysis. Next, consenus clustering was classified inflammation subgroups to analyze immune infiltration. Furthermore, a prognostic inflammation-associated model was established, and its corresponding stability was verified with various methods. Ultimately, immunological characteristics and drug sensitivity analysis extended on the idea of clinical application, while the correlation between IRGs and tumor microenvironment in virus-related HCC was attested.

The expression matrix of 200 IRGs was adopted to compute the euclidean distance between any two patients from 179 virus-related HCCs, and t-SNE algorithm was further condensed the euclidean distance into two-dimensional points. Subsequently, three clusters with virus-related HCC patients were generated and each patient was assigned to its closest ( Figure 2A ), namely 81, 57, and 41 patients in distinct clusters (Cluster I, Cluster II, and Cluster III), respectively. OS analysis displayed that the most significant differences consisted between cluster I and cluster II. Thus, patients in Cluster II yield the best OS while those in Cluster I had the worst prognosis outcome ( Figure 2B ), indicating that Cluster II and Cluster I might represent the lowest and highest status of inflammation. Accordingly, sufferers in Cluster I and Cluster II were classified into “inflammation^high” and “inflammation^low” groups, separately. To obtain inflammation-related DEGs, expression profiles were compared between the inflammation^high and inflammation^low groups, leading to a total of 47 inflammation-related DEGs identified ( Figure 2C ). Next, A PPI network was constructed, composed of 84 nodes. Among all nodes, 10 hub genes, including CCL20, IL1B, CCL2, CCL22, TIMP1, LIF, TLR3, F3, LTA, and PLAUR were distinguished ( Figure 2D ). Further research found that the DEGs were mostly enriched in immune response, integral component of membrane, and signaling receptor activity ( Figure 2E ; Table S3 ). KEGG analysis also demonstrated that these DEGs were closely related to pathways in inflammation, such as IL-17 signaling pathway, TNF signaling pathway, and NF-kappa B signaling pathway ( Figure 2F ; Table S3 ).

The relevance network of IRGs interactions, regulator relationships, and corresponding survival status in virus-related HCC patients was presented in Figure 3A and Table S4 . To further conclude the relation between expression profiles of IRGs and HCC subtypes, a consensus clustering analysis was conducted to separate patients into different gene clusters based on the expression levels of the IRGs ( Figure 3B ). Three discrepant patterns were determined: 98 cases in Cluster 1, 53 cases in Cluster 2, and 28 cases in Cluster 3 ( Figure 3C ). Afterwards, OS status of the three patterns was revealed, contributing to a consequential difference observed ( Figure 3D ). Additionally, the genomic expression and clinicopathological features of three clusters were displayed in Figure 3E , identifying a substantial difference between IRGs expression and clinical characters.

The CIBERSORT and ESTIMATE algorithms were implemented to confirm the activity or enrichment levels for immune cells in virus-related HCCs ( Table S5 ). The heatmap of three independent immune cell infiltration (ICI) subtypes was presented in view of 179 tumor samples with matched ICI profiles from TCGA-LIHC ( Figure 4A ). The expression of three vital ICPs (PD-1, PD-L1, and CTLA-4) was obviously distinct among three clusters. In virtue of the role of TME scores for evaluating the abundance of immune and stromal elements in TME, the ESTIMATE algorithm was further executed to estimate the TME scores, involving stromal score, immune score, and estimate score, in three clusters, the results of which turned out sufferers in cluster 3 yielded superior TME scores ( Figures 4B, C ). In addition, we explored if the three subclasses generated various tumor immune microenvironments (TIME) ( Figure 4D ). Indeed, the immune-high subgroup had high infiltration levels of Eosinophils, Macrophages M0, Macrophages M1, and Neutrophils, while the cluster 1 had remarkable enrichment of resting mast cells.

Firstly, Univariate Cox regression analysis was utilized on the virus-associated HCC groups, demonstrating that 13 prognosis-related IRGs were correlated with OS ( Figure 5A ; Table S6 ). To prevent model overfitting, LASSO penalized Cox regression modeling and GMM model were simultaneously conducted to filter the vital DEGs, which were positively associated with the prognosis of HCC patients. With the joint method, a novel prognostic gene model with six hub genes (ADORA2B, CCL20, ICAM4, MEP1A, SLC7A2, and TNFSF9) was constructed ( Figures 5B–E ). Then, we computed risk score using the following formula: risk score =∑iCoefficient (mRNAi) × Expression(mRNAi), where i, stands for the expression of six key IRGs. In line with the median risk score, samples were clustered into low- and high-risk subgroups. The distribution patterns from PCA analysis illustrated that patients could be distinguished into high- and low-risk classes ( Figure 5F ). Also, the risk plot of IRG score proved that OS time decreased while mortality rise, as IRG score increased. And survival analysis testified that samples in the low-risk cluster produced significantly longer OS time in comparison with that of the high-risk patients ( Figure 5I , P<0.01, log-rank test) ( Figures 5G–I ; Figure S1 ). Moreover, the expressive relationship among them and heatmap of selected genes were displayed in Figure 5J and Figure 5K , respectively. To comprehend the relationship between immune subtypes and IRG score, an alluvial diagram was drawn for clusters with distinct risk-subgroups, and accompanying survival status ( Figure 5L ). The outcomes demonstrated that cluster 3 with higher IRG score was most likely associated with death. Whereas the cluster 1 exhibited a lower IRG score and best prognosis status. A time-dependent ROC curve was further performed and the area under the curve (AUC) reached 0.805, 0.7, and 0.718 at 1, 3, and 5 years, respectively ( Figure 5M ). Besides, the ROC curve explained that the predictive OS accuracy of IRG score was superior to other clinical parameters (Age, gender, Alcohol consumption, Neoplasm histologic grade, and TNM stage and age) ( Figure 5N ). Tremendous differences in the IRG score of three clusters were discovered ( Figure 5O ), implying a higher IRG score may be relevant with immune activation-associated features.

To explore the relation between the IRG score and clinicopathological Characteristics, the interaction between IRG score and multitudinous clinical parameters (Age, Alcohol consumption, gender, TNM stage, Fetoprotein, Radiation therapy and survival status) was discussed ( Figures S2A–F ). We perceived that IRG scores increased along with the stage III-IV and higher level of fetoprotein. And Univariate and multivariate Cox regression analyses were further conducted to guesstimate the accuracy of the risk model and disclose whether IRG score could be considered as an independent prognostic factor for patients’ prognosis. Accordingly, Univariate Cox regression analysis revealed that both the IRG score and the stage were significantly correlated with OS of the patient ( Table S7 ). Furthermore, to uncover the prognostic significance of IRG score in virus-related HCC patients, the patients were assigned into different subgroups based on above clinical parameters. Totally, the high-risk patient’s survival was generally poorer compared to low-risk patients ( Figures S2A–F ).

As disclosed in Figure 6A , a nomogram was reciprocally constructed on the foundation of IRG scores, combined with clinical features. Followly, calibration curves defined the reliability and accuracy of nomogram to predict 3-, and 5-year prognosis ( Figures 6B–D ). As shown in Figures 6E–H , the AUC values were as expected, implying this nomogram had an excellent predictive ability for prognosis. Moreover, we also found that this prognostic model with diverse clinical factors presented more net benefits for predicting the prognosis.

We aimed to assess the relevance between IRG score and immune cells abundance with the CIBERSORT algorithm. As depicted in Figure 7A and Figure S3 , the IRG score was significantly associated with the infiltration of B cells memory, Eosinophils, Mast cells activated, Monocytes, Plasma cells, T cells CD4 memory activated, T cells CD4 memory resting, B cells naive, Macrophages M2, Dendritic cells activated, NK cells resting, T cells CD4 naive, T cells gamma delta, T cells CD8, and T cells follicular helper, while the negative performance appeared in relationship with IRG score and Dendritic cells resting, Macrophages M0, and Neutrophils. Then, the correlation between immune cell infiltration and expression status of six genes incorporated with the prognostic model construction was analyzed in Figure 7B . Also, high-risk patients experienced higher EstimateScore and StromalScore levels than those in low-risk group (p<0.05) ( Figure 7C ). Meanwhile, IRG score was positively associated with the expression of a series of immune checkpoints (such as CD200, CD70, and PDCD1) ( Figure 7E ) and the enrichment scores of immunotherapy response-related gene signatures ( Figure 7D ). Furthermore, we assessed the relationship between ICPs and risk group, demonstrating that ICPs (PD-1, LAIR1, and VTCN1, et al) were inconsistently distributed in two risk clusters ( Figure 7F ).

Present studies declared that ICP inhibitors were favorable to populations with increased TMB or higher MSI, uncovering that TMB and MSI were both ponderable indexes for predicting tumor immune response (38, 39). As Figure 8A demonstrated, that TMB in the low-risk cluster was higher than high-risk cluster suggested that immunotherapy provided more benefits for patients with high risk. And a negative correlation was drawn between IRG score and TMB in spite of the meaningless value (R=-0.15, p=0.08, Figure 8B ). To explore the impact of TMB status on prognosis in virus-related HCC patients, we also conducted survival analysis in various TMB classes. No significant difference of prognosis was revealed between High-TMB patients and low-TMB patients ( Figure 8C ). However, the survival analysis for combination of TMB and IRG score for virus-related HCC patients drew a conclusion that the prognostic benefit in the high-TMB group was eliminated by IRG score ( Figure 8D ). The measurement of RNA stemness score (RNAss) could represent cancer stemness, based on mRNA expression (40, 41). The relevance of IRG score and CSC score was presented in Figure 8E . Likewise, lower IRG score was connected with MSI-H pattern, while higher IRG score was linked with the microsatellite stable (MSS) pattern ( Figure 8F ), which also illustrated that low-risk patients may be more susceptive to immunotherapy. Meanwhile, genomic alterations in high and low groups were further analyzed. A rough similarity in the kinds of the top 30 genes with the highest mutation frequency between the low and high groups emerged ( Figures 8G, H ).

TIDE scores and IPS scores were conducted to make prediction for sufferers’ responsiveness for appraising the immune response of virus-related HCC patients. Analysis results in Figures 9A,B revealed that patients at low-risk generated a lower TIDE score and a higher IPS score, which demonstrated that they may suffer more sensitivity from immunotherapy (42, 43). In the following, aimed at analyzing the clinical application of IRG model, we calculated the alterations in terms of drug sensitivity between diverse risk clusters, reflecting that 5-fluorouracil, AZ628, AZD7762, Bortezomib, Camptothecin, Cisplatin, Cyclopamine, Dasatinib, Docetaxel, MG-132, Nilotinib, Obatoclax Mesylate, Paclitaxel, PHA-665752, Sunitinib, Vinblastine, Vorinostat, and VX-680 yielded advantageous effectiveness for low-risk patients ( Figure 9C ).

qRT-PCR was applied to verify the mRNA expression levels of six hub genes in plasma samples from 58 virus-related HCC patients and 50 normal people. The unpaired t-test was performed to compute the differences between the virus-related HCC plasma samples and normal plasma samples. And plasma samples validated that the significant differences existed in the expression levels of CCL20, and SLC7A2 between HCC and normal tissues, that is, CCL20 was highly expressed in most HCC plasma samples while the expression of SLC7A2 was significantly higher in normal plasma samples than in virus-related HCC plasma samples ( Figure 10A ). Subsequently, we extracted total RNA from different tumor cell lines (SMMC7721, Huh7, HepG2, and HCCLM3) and the normal liver cell lines (WLR68, and LO2) to measure the mRNA expression levels of CCL20, and SLC7A2. qRT-PCR assays were implemented and the results showed that the mRNA expression levels of CCL20 were significantly higher in liver cancer cells, like SMMC7721, Huh7, HepG2, and HCCLM3, than in normal cells, such as WLR68, and LO2. However, the consequence for SLC7A2 expression level was adverse ( Figures 10B, C ).

To observe the function of CCL20 and SLC7A2 during hepatocarcinogenesis, CCL20 was silenced by transfection with shCCL20 in SMMC7721 cells. In addition, SLC7A2 was further overexpressed based on the SMMC7721 cells with silence of CCL20. First, we successfully constructed two lentiviral vectors harboring shRNA-CCL20-1, and shRNA-CCL20-2, respectively, and established two stable knockdown cell lines in SMMC-7721. The two different shRNAs, especially shRNA-CCL20-2, effectively knocked down the expression of CCL20. Also, SLC7A2 overexpression was indicated in vitro ( Figures 10D, E ). Afterwards, we examined the effect of shCCL20, overSLC7A2, and shCCL20/overSLC7A2 on the morphology of SMMC7721 cells. Compared with control cells, shCCL20 cells, and overSLC7A2 cells, showing a spindle-like shape with scattered growth, cells with knockdown of CCL20 combined with SLC7A2 overexpression induced a cobblestone-like appearance with the significant dispersion change ( Figure 10F ). Furthermore, the study established nude mouse tumor xenograft models injected by SMMC7721 cells that were transfected with shCtrl or shCCL20/overSLC7A2. Tumor volume was measured every 3 days ( Figure 10G ). We found that CCL20 knockdown combined with SLC7A2 overexpression significantly lessened the tumor volume ( Figure 10H ). After 30 days, we measured the tumor weight and observed that tumor weight was distinctly lowered by CCL20 knockdown combined with SLC7A2 overexpression ( Figures 10I, J ).