Field cancerization profile-based prognosis signatures lead to more robust risk evaluation in hepatocellular carcinoma

Summary The development of reliable biomarkers has been an urgent issue as well as a hot spot of research on the diagnosis, treatment, and prognostic evaluation of hepatocellular carcinoma (HCC). Here, we established and validated two field cancerization profile-based prognostic signatures (gene expression score [GES] and immune score [IS]) for HCC. Our study confirmed that field cancerization profile-based models outperform conventional models on risk evaluation, offering insights for further studies on prognostic model construction. The nomogram constructed by combining GES, IS, and TNM stage was proved effective in improving the individualized prediction of the overall risk of patients. Distinct peritumoral characteristics were observed in several immune cells (e.g., CD8 T cells and dendritic cells), which might explain the diversified prognosis and clinical benefit of immunotherapy. Moreover, a series of drug targets, prognosis-associated genes, and pathways were identified, which may contribute to molecular mechanism studies as well as therapeutic target development of HCC.


INTRODUCTION
Hepatocellular carcinoma (HCC) is the fourth most common cause of cancer-related death worldwide, and its incidence is steadily increasing throughout the world (Bray et al., 2018;Yang et al., 2019). Unlike other solid malignancies, most HCC arises in a background of chronic liver disease, variably complicated with hepatitis and fibrosis, which makes HCC an occult, highly malignant disease with a high recurrence rate and poor prognosis (Zhu et al., 2013). Therefore, early screening and prevention, as well as outcome monitoring are of great importance to patient survival. The development of reliable biomarkers has been an urgent issue as well as a hot spot of research on the diagnosis, treatment, and prognostic evaluation of HCC.
With advances in global molecular profiling based on high-throughput methods, gene expression profiling came as a powerful diagnostic and prognostic tool in molecular medicine. Comparing to traditional staging system, the access to gene expression features enabled the development of specific biomarkers, thus providing new opportunities to disentangle the diversity and heterogeneity of cancers (Dakubo et al., 2007). Recently, gene expression profiling has been widely applied to the prediction of prognosis and response to therapy in a range of cancers, e.g., prostate cancer and breast cancer (Hofmann et al., 2002;Sorlie et al., 2003). In addition to gene expression profile, the immune microenvironment, which indicates the infiltration level of immune cells (e.g., T cells, B cells, and macrophages) was also found to be an indicator of clinical outcome and has been applied as biomarkers for multiple cancers (Bremnes et al., 2016;Song et al., 2019;Wang et al., 2019). However, despite that a multitude of prognostic signatures have been reported in HCC, none of these has become a tangible tool in clinical decision making owing to poor clinical reproducibility Zucman-Rossi et al., 2015). Therefore, it is crucial to identify new prognostic markers for HCC, thereby promoting individualized therapeutics for patients.
Owing to an increasing interest on tumor-host interactions, in-depth studies on the rewiring of biogenesis and metabolic and signaling pathways in an integrated level is arousing more and more interest in the past decade. The concept ''field cancerization'' was first proposed by Slaughter et al. by referring to the phenomenon that the normal tissue peripheral to malignant cells was histologically abnormal and therefore was part of the transformed cells in a particular tumor field, which may consequently lead to the recurrence of the local tumor (Dakubo et al., 2007;Slaughter et al., 1953). According to field cancerization theory, performance of the GES signature was also fully verified in two validation cohorts (LIHC cohort, 24-month survival prediction AUC = 0.70 (95% CI 0.54-0.87), 36-month survival prediction AUC = 0.79 (95% CI 0.64-0.93); LIRI cohort, 24-month survival prediction AUC = 0.65 (95% CI 0.54-0.76), 36-month survival prediction AUC = 0.59 (95% CI 0.46-0.71)) ( Figures 3C and 3D).
To validate the superiority of the GES model over models constructed from tumor-only or peritumor-only profiles, we constructed two types of control models. Type I control model was the tumoral or peritumoral part of GES (namely, TGS1 and PGS1), consisting of the tumor or peritumoral genes from GES with the same coefficient. Type II control model was built from the tumor or peritumoral expression profile (namely, TGS2 and PGS2), with marker genes selected using the same two-step procedure as GES. These two types of control models were necessary as they test the robustness of the field cancerization method in a complementary way. According to the result, although four models could significantly predict the survival of patients in the GSE14520 cohort (log rank test p < 0.001), none of them succeeded in splitting each validation cohort into groups with a significant difference in OS ( Figure S1). As for PGS1, its optimal cutpoint determined by maxstat was unable to separate the patients from the LIHC cohort into two datasets ( Figure S1).   (17) somatostatin receptor signaling (23) hormone mediated apoptotic signaling (21) opioid receptor signaling (3) serotonin receptor signaling (7) negati ve regulation of IL6 mediated signaling (19) positi ve regulation of b cell receptor signaling (38) positi ve regulation of toll li ke receptor 2 signaling (12) toll li ke receptor 2 signaling (1) regulation of toll li ke receptor 2 signaling (15) negati ve regulation of toll li ke receptor 2 signaling (11) negati ve regulation of toll li ke receptor 4 signaling (25) regulation of sphingolipid biosynthetic (9) dendritic cell apoptotic (31) negati ve regulation of amyloid precursor protein biosynthetic (32) negati ve regulation of phospholipid biosynthetic (37) adenylate cyclase inhibiting adrenergic receptor signaling (40) negati ve regulation of vascular associated smooth muscle cell apoptotic ( ++ + + +++ + + + + ++ + + ++ + + ++ + + + + + + + + + + + + + + + + + + + + + ++ + + + + + + + + ++ + + + + + + + + + + + + ++ + Our main objective was to explore the role of field cancerization profile in the risk evaluation of HCC. Prognostic models were established using a two-step procedure. Gene expression and immune cell fraction profile were applied to the model construction separately. The training cohort was derived from GSE14520. Two independent cohorts, the LIHC cohort from TCGA and the LIRI cohort from ICGA, were used to evaluate the external performance of each model. Only the patients with expression information profiles from both tumor and peritumor were included in the study. GES, gene expression score; IS, immune score; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; PEC, prediction error curve; DCA, decision curve analysis; GSVA, gene set variation analysis.  Article   CXCL10  CD52  IL2RB  IRF1  TBCCD1  ZSCAN31  SECTM1  KDM8  ZNF24  ZNF652  P2RX4  VEGFA  CYLD  DEF6  IDO1  LSM1  WHRN  NDUFA3  NECAB3  ACE2  PATJ  ZNF446  PRKRIP1  RBM28   IMPAD1  NMB  CDKN3  ICOS  ENTPD2  ZNF230  KLHL36  CATSPERG  CDH13  ANGPTL2  MFGE8  MSX1  DNMT3A  NCAM1  SPRED2  TGM1  FAIM  RPP30  ANAPC15  DAB2  SNRPF  PRIM2  ZNF415  MYH10  SETD4  DAPK3  TCTN2  ZFR  ZNF124  NEDD9  SMTN  SPATA5L1  CUEDC2  RNF111  PIBF1  GPATCH2  VPS45  CEBPD  GTF2F1  TNFAIP1  IQCK  RPLP2  BAX  RSU1  NUDT2  We further validate the superiority of field cancerization profile in model construction in an unbiased way. With 1,000 training sets subsetted from the GSE14520 cohort, we reconstructed the five models (GES, TGS1, TGS2, PGS1, and PGS2) and each was validated in two validation cohorts. During each validation, we reassessed the optimum cutpoint of each model and calculated the p value between the survival of two groups separated by the cutpoint. The p values of the GES model were compared with that of others using the pairwise Wilcoxon rank-sum test. According to the result in the LIRI cohort, the p values of GES models were significantly lower than that of other models (GES versus TGS1 p = 5.82e-19; GES versus TGS2 p = 5.78e-39; GES versus PGS1 p = 1.09e-05; GES versus PGS2 p = 2.67e-14, Figure 3E). Similar results were also obtained in the LIHC cohort (GES versus TGS1 p = 5e-24; GES versus TGS2 p = 8.58e-12; GES versus PGS1 p = 0.0182; GES versus PGS2 p = 0.361, Figure S2). These results proved that the GES model has better performance over models constructed from tumor-only or peritumor-only profiles.

OPEN ACCESS
iScience 25, 103747, February 18, 2022 5 iScience Article GES signature seems to be more promising in predicting the prognostic OS outcome of patients with HCC ( Figure 3F).
We obtained, in addition to gene expression information, a series of clinical indices for each patient, e.g., age, gender, TNM stage, hepatitis B virus (HBV), alanine aminotransferase (ALT), multinodular, cirrhosis, alpha-fetoprotein (AFP), and tumor size. Each clinical variable was applied for univariate Cox regression analysis, and we found that patients from the TNM stage II and TNM stage III subgroups got a significantly shorter OS than those from the TNM stage I subgroup (TNM stage II: HR = 2.04, 95% CI = 1.15-3.60, log-rank test p = 0.015; TNM stage III: HR = 5.89, 95% CI = 3.30-10.51, log-rank test p < 0.001) (Figure 4). Several clinical features (e.g., TNM stage and AFP level) were significantly different between the GES-high group and GES-low group (Table S2). Furthermore, GES levels were included in multivariate Cox regression to test the independence of the GES score over clinical indices. After correcting for TNM staging, age, gender, and other clinicopathological characteristics, the GES score remained a significant prognostic risk factor (GSE14520 cohort, HR = 0.07, 95% CI = 0.04-0.13, log-rank test p < 0.001; LIHC cohort, HR = 0.14, 95% CI = 0.

Immune cell infiltration spectrum across tumor and peritumor of HCC
With the CIBERSORTx algorithm, we systematically evaluated the fractions of 22 immune cells in both tumor and peritumor samples ( Figure 5). (Newman et al., 2019) Of the 22 immune cell components, 16 immune cells got a biased immune fraction (p < 0.05, Wilcoxon test) between the tumor group and peritumor group. Nine immune cells (plasma cells, naive CD4 T cells, activated memory CD4 T cells, regulatory T cells, macrophages [M0], activated dendritic cells, resting mast cells, neutrophils, and resting NK cells) got a significantly higher infiltration level in the tumor, and seven immune cells (CD8 T cells, T cells follicular helper, resting dendritic cells, macrophages (M1), activated NK cells, activated mast cells, and monocytes) got a significantly higher infiltration level in peritumor. According to the correlation analysis, the pairwise association of infiltration level between immune cells was generally consistent in tumor and peritumor, with several exceptions ( Figure 5A). For example, the infiltration level of macrophages M0 in tumor was negatively associated with tumoral activated dendritic cells, and it turned out to be a positive association in peritumor. Together these results revealed the distinct immuno-characteristics in peritumor tissue, which might contribute to a more precise prognosis prediction.
According to the analysis of association with clinical outcomes, 21 immune cells (13 in tumor and 8 in peritumor) appear to be significantly correlated with patient survival in univariate Cox regression ( Figure 5B).

For the 13 tumor prognosis-associated immune cells, the high infiltration level of 8 cells (B cells naive, T cells CD8, T cells follicular helper, T cells regulatory, NK cells activated, macrophages M1, and macrophages M2) is associated with better prognosis, as opposed to the performance of plasma cells, T cells CD4 naive, macrophages M0, dendritic cells resting, eosinophils, and neutrophils. For the eight peritumor prognosis-associated immune cells, five cells (B cells naive, T cells CD8, dendritic cells resting, dendritic cells activated, and mast cells resting) got a positive association between infiltration level and prognosis,
as opposed to the performance of macrophages M0, mast cells activated, and eosinophils. By comparison, several distinctive features were observed in peritumoral immune cells. For example, macrophages M1 and macrophages M2 in tumor were both significantly associated with high OS of patients (p < 0.05); however, this was not observed in peritumor immune cells. T cells CD8 were observed to have a beneficial effect on patient survival (p = 0.02), and a higher significance was observed in peritumor (p < 0.001), suggesting that T cells CD8 may play a greater role in peritumor ( Figure 5B). Moreover, the fraction of tumoral dendritic cells component was significantly (highly significant in dendritic cells resting (HR = 0.4 (0.23-0.69), p < 0.001) and marginally significant in dendritic cells activated (HR = 0.56 (0.30-1.04), p value = 0.067)) and positively associated with poor prognosis, suggesting a deleterious effect on clinical outcome, whereas a totally opposite beneficial effect was observed in peritumor samples (dendritic cells resting [HR = 1.97 (1.27-3.06), p value = 0.003]; dendritic cells activated [HR = 2.92 (1.07-7.98), p value = 0.037]) ( Figure 5B). Furthermore, we explored the interactions between gene expression and prognosis-associated immune cells ratios. By GO enrichment analysis ( Figure S3-S5), the gene sets that were significantly associated with the level of infiltration of each cell type were analyzed. For T cells CD8, most enriched GO terms were shared between tumor and peritumor (e.g., T cell activation, T cell receptor signaling pathway, and cell chemotaxis), which may explain their similar role in HCC. In contrast, several unique features (e.g., gene sets associated with neutrophil activation and neutrophil degranulation) were observed in peritumoral dendritic cells (resting). These results demonstrated that aberrant immune infiltration in not only the tumor itself, but also in the peritumor, might play important roles in the development of HCC.
By using a LASSO Cox regression, an immune score (IS) prognostic model was built from 22 survival-associated immune cells at the field cancerization level ( Figure 6A, Table S3). Of the 22 candidate immune cells, 8 tumoral immune cells and 4 peritumor immune cells were selected for model construction. Immune scores were then calculated for each patient using the infiltration level of each cell weighted by their regression coefficient (Table S3). The IS was associated with worse OS. We used R package maxstat to assess the optimum cutoff for immune scores in the GSE14520 cohort, which results in the highest log-rank statistic at 1.96 ( Figures 6B-6D). Patients with HCC in the GSE14520 cohort were then divided into risk low and risk high immune score groups. Of the 12 immune cells, the higher fraction of 4 was beneficial to clinical outcome and deleterious effect was observed in 7 immune cells. Kaplan-Meier survival analysis demonstrated that patients in the IS-low group had significantly longer OS than those in the IS-high group (HR = 0.28, 95% CI = 0.18-0.45, log-rank test p < 0.001; Figures 6E and 6F).      As with GES signature, we sought to validate the robustness of IS signature from several aspects. In the LIHC cohort, the patients could be stratified according to their IS score and OS time was markedly elongated in the low-risk group (HR = 0.28, 95% CI = 0.18-0.45, log-rank test p = 0.001, Figure 7A). The risk score was further validated in the LIRI cohort (HR = 0.33, 95% CI = 0.18-0.62, log-rank test p < 0.001, Figure 7B). We next carried out ROC analysis to assess the performance of IS signature in predicting 24-, and 36-month OS of patients from each cohort. The AUC of IS signature was high in all three cohorts (GSE14520 cohort, 24-month AUC = 0.75 ( Figure 6F, Figure 7C and 7D). These results proved the feasibility to determine patient prognosis based on immune cell infiltration. More importantly, our work demonstrated the diagnostic potential of IS signature for survival prediction in clinical application.

B cells naive B cells memory Plasma cells T cells CD8 T cells CD4 naive T cells CD4 memory resting) T cells CD4 memory (activated) T cells follicular helper T cells regulatory (Tregs) T cells gamma delta
We further compared the performance of IS model with four alternative models (TIS1, TIS2, PIS1, and PIS2). Four models with their corresponding cutpoints were able to stratify the GSE14520 cohort into two groups, and the patients from the low-risk group all got a significantly higher OS (p < 0.001). However, they failed to make proper stratification on the OS of patients from the LIHC and LIRI cohorts, with no significant difference in OS in groups divided by each model ( iScience Article We next assessed the prognostic association between IS signature and other known clinical risk factors. In the analysis of univariate Cox regression, the IS signature was found to be a significant indicator for OS, together with several already well-known risk factors (e.g., TNM stage, AFP, tumor size). Furthermore, we further analyzed the independence of IS by multivariate Cox regression. After correcting for TNM staging, age, gender, and other clinicopathological characteristics, the IS score remained a significant prognostic risk factor (GSE14520 cohort, HR = 0.33, 95% CI = 0.20-0.55, log-rank test p < 0.001; LIHC cohort, HR = 0.31, 95% CI = 0.12-0.82, p = 0.018; LIRI cohort, HR = 0.38, 95% CI = 0.20-0.75, p = 0.005, Figure 8), indicating the independence of IS signature in predicting patient survival. Taken together, these findings suggest that the risk score retains its prognostic relevance even after the classical clinicopathological prognostic features have been taken into account.23

Construction and evaluation of a prognostic nomogram
To provide a clinically relevant quantitative method to predict the probability of 24-and 36-month OS in patients with HCC, we constructed a nomogram that integrated the GES, IS, and TNM stages ( Figure 9A). iScience Article Calibration plots showed that the nomograms performed well compared with the performance of an ideal model ( Figure 9B). The C-index of the nomogram was 0.834 (SE = 0.02).
As expected, the result of decision curve analysis showed that the nomogram got a higher accuracy of prediction on 24-and 36-month OS, indicating that the nomogram was more clinically useful than using GES classifier, IS classifier, and TNM stage alone ( Figures 9C and 9D). ROC analysis to compare the sensitivity and specificity of the nomogram comprising the IS signature with clinicopathological risk factors showed that overall survival was more accurately predicted by the nomogram than by the risk factors in all three cohorts ( Figures 9E-9J).    (Table S5). Of these drugs, melatonin was selected as a full agonist of IDO1. Many studies have shown that melatonin's co-administration improves the sensitivity of cancers to inhibition by conventional drugs. For example, the combination of melatonin and an IDO inhibitor DL-1MT improves the efficacy of an immunotherapy (gDE7) targeting HPV-associated tumors (Moreno et al., 2018). Moreover, besides its effects on cancer cells, the protective effect of melatonin was observed in cells in the liver and the immune systems (Calvo et al., 2013;Carrillo-Vico et al., 2013;Molpeceres et al., 2007;Ren et al., 2017;Sardo et al., 2017), which further suggests the peritumoral effectiveness of melatonin on HCC.

Functional relevance and clinical potential of diagnostic markers for HCC
To identify the pathophysiological mechanisms underlying the difference in survival time in patients with HCC, we performed gene set variation analysis (GSVA) on both tumor and peritumor. From GSE14520 we selected 141 patients with consistently high or low risk levels in GES and IS classifiers, including 50 high-risk patients and 91 low-risk patients. From the tumor, we identified 36 gene sets with higher enrichment score (ES) in low-risk patients and four gene sets with higher ES in high-risk patients. These gene sets covered several proliferation-, angiogenesis-, and metastasis-associated pathways, e.g., WNT signaling pathway, MDA5 signaling pathway, blood coagulation, and fatty acid metabolism-associated pathways ( Figure S7). As for peritumor, we observed 33 gene sets with higher ES in high-risk patients and 7 gene sets with higher ES in low-risk patients ( Figure 10). Two significantly enriched gene sets, IL-6-mediated signaling pathway (log2 fold change = -0.50; padj = 2.91e-3) and G-protein-coupled receptor (GPCR) signaling pathway (log2 fold change = 0.46; padj = 1.34e-3, Figure 10), were identified in peritumor, whereas these two pathways were marginally significant in gene set variation analysis in tumor (padj = 5.8e-2 for IL-6-mediated signaling pathway and padj = 9.6e-2 for adenylate cyclase activating GPCR signaling pathway, respectively). This indicates some similarities between tumor and peritumor during tumorigenesis. Moreover, these two pathways got a higher log2 fold change than that in tumor (IL-6-mediated signaling pathway, -0.34 in tumor versus -0.50 in peritumor; adenylate cyclase activating GPCR signaling pathway, 0.30 in tumor versus 0.46 in peritumor), suggesting a stronger effect of these pathways in peritumor. The Toll-like receptor 2 (TLR2) signaling pathway, which plays an essential role in the innate immune response, got significantly higher expression in patients with high risk (log2 fold change = -0.38, padj = 1.7e-4, Figure 10). Notably, such biased expression TLR2 pathway gene sets were not observed in tumor (log2 fold change = -0.005, padj = 0.98). To some extent, the above-mentioned molecular characteristics may provide a reasonable interpretation of the prognostic value of the two signatures in our study and may further give us insights into the underlying mechanisms of tumor prognosis from the perspective of both tumor and peritumor.

DISCUSSION
HCC is a leading cause of cancer-related death in many parts of the world. It is characterized as an occult, highly malignant disease with a high recurrence rate and poor prognosis. Unlike other solid malignancies, the coexistence of inflammation and cirrhosis makes the early diagnosis and prognostic assessment of HCC much more difficult. Despite recent progress on the diagnosis of and therapy for HCC, its morbidity and mortality is still rising in many countries, which indicates an urgent need for valuable biomarkers for the diagnosis and treatment of HCC (Akoad and Pomfret, 2015;Zhu et al., 2013).
With the emergence of high-throughput technique as a powerful method for the profiling of gene expression, numerous genomic studies tried to provide gene expression-based prognostic markers for HCC (He et al., 2020;Liang et al., 2020;Zhang et al., 2020). However, recent prognostic markers do not inform the diagnosis, prognosis, or treatment decisions for HCC patients. Limitations of these markers include the lack of large-scale, independent validation and, more importantly, the neglect of peritumor features. In the current work, based on the theory of field cancerization, we selected candidate features from tumor and peritumor simultaneously. Through univariate Cox regression, we were able to uncover the detail of prognosis-associated gene and immune cell subsets in HCC. A series of prognosis-associated features iScience Article were captured, including a significant association between the fraction of CD8 T cells and patient survival in both tumor and peritumor tissue. More notably, a stronger association was observed between peritumor CD8 T cells and patient survival (tumor, HR = 1.73, p = 0.02; peritumor, HR = 3.27, p < 0.001, Figure 5B), suggesting a greater role of peripheral originated T cells CD8. Recent evidence showed that the clonotypic expansion of effector-like T cells within not only the tumor but also in normal adjacent tissue, and patients with gene signatures of such clonotypic expansion respond best to anti-PDL1 therapy, indicating that nonexhausted T cells and T cell clones supplied from the periphery may be key factors in explaining patient variability and clinical benefit from cancer immunotherapy (Wu et al., 2020). Put these together, a similar mechanism is likely shared in a range of solid malignancies, including HCC, which may hold promise for further therapies.
As an effort to overcome the underlying risk of overfitting, we applied a Cox regression model with a LASSO penalty for the selection of prognosis-associated features. Based on the cancerization field features from the GSE14520 cohort, the expression of 70 genes and infiltration levels of 22 immune cells were selected for the construction of GES and IS, respectively. Both signatures were able to perform robust prognosis-based classification in not only the training cohort but also two validation cohorts, indicating good reproducibility of the two positive regulation of adenylate cyclase activating GPCR signaling (6) vitamin d3 metabolic (10) cellular alcohol metabolic (17) somatostatin receptor signaling (23) hormone mediated apoptotic signaling (21) opioid receptor signaling (3) serotonin receptor signaling (7) negative regulation of interleukin 6 mediated signaling (19) positive regulation of b cell receptor signaling (38) positive regulation of toll like receptor 2 signaling (12) toll like receptor 2 signaling (1) regulation of toll like receptor 2 signaling (15) negative regulation of toll like receptor 2 signaling (11) negative regulation of toll like receptor 4 signaling (25) regulation of sphingolipid biosynthetic (9) dendritic cell apoptotic ( iScience Article signatures in HCC. Further analysis uncovered that the prognostic value of signature in our study was independent of the main prognostic factors in HCC. For example, AFP is most commonly employed in the clinic for HCC screening and as an important predictor for patient survival (Peng et al., 2004). The TNM staging system reports the pathological characteristics of patients (e.g., tumor grade, nodule, and metastasis), which had been widely applied in various cancers (Marrero et al., 2018). By performing multivariable Cox regression analysis, we identified that the prognostic values GES and IS are independent of AFP and TNM stage, which strongly demonstrated that each model could act as an independent prognostic factor for HCC. Furthermore, our study confirmed that field cancerization profile-based models could better predict a patient's prognosis. In this study, a series of tumor or peritumor profile-based models were trained and compared with GES and IS, which confirmed the advantage of field cancerization profile-based models. We also observed a lower prediction error in GES compared with other published HCC prognosis models ( Figures 3E and 3F). Considering that recent efforts on prognostic model development are mostly based on tumor features alone, the neglect of the peritumor profile might lead to their lack of robustness and limited clinical application. To provide a convenient and accurate tool for the evaluation of overall risk (including tumor recurrence and metastasis) at the individual level, we constructed a nomogram by integrating GES, IS, and TNM stage information. Although such methods generally use traditional prognostic factors, such as TNM stage and sex, we propose including our field cancerization profile-based model, which reflects the molecular heterogeneity of these tumors. The performance of the nomogram was verified in all validation cohorts. Thus, our nomogram might provide a simple and accurate method for predicting prognoses in HCC.
Cancer is a heterogeneous disease, and exploring the dysregulated genes involved in carcinogenesis and development might help to improve prognostic and therapeutic strategies. According to MalaCards (Rappaport et al., 2017), 28 genes from the GES signature are associated with human cancer, including 17 tumoral genes (CYLD, TGM1, IL2RB, NEDD9, VEGFA, IRF1, MFGE8, CDH13, DAB2, MSX1, DAPK3, CXCL10, CD52, LSM1, SMTN, NMB, and ZNF652) and 11 peritumoral genes (PIBF1, BAX, GPATCH2, CDKN3, ACE2, CUEDC2, DNMT3A, EGR2, IDO1, NCAM1, and ICOS). Among the 11 peritumoral genes, Bax (Bcl2-associated X protein) encoded by gene BAX was reported as a central apoptosis-inducing protein that works through the induction of permeabilization of the outer mitochondrial membrane (Benard et al., 2010). Bax interacts with Bcl-2 and their ratio determines survival or death following an apoptotic stimulus (Oltvai et al., 1993). Specifically, the overexpression of BAX had been reported to induce autophagy and enhance drug sensitivity of HCC in cell experiments. The CDKN3 gene encodes a cyclin-dependent kinase inhibitor, which dephosphorylates CDK2 kinase to prevent its activation. It has been reported to be deleted, mutated, or overexpressed in several kinds of cancers. It has been reported that the overexpression of CDKN3 in HCC could promote cell proliferation by stimulating G1-S transition. In this study, the high expression of peritumoral CDKN3 was associated with poor prognosis (coefficient = 0.466). Two CDKN3 antagonists, AT7519 and Kenpaullone, were identified through the search in GeneCards (Table S5). Our results suggest that these drugs may take effect in peritumor as well. In our research, a limited number of drugs were identified in peritumor genes, which is mainly due to the lack of study on the role of peritumor in cancerogenesis. Although the biological function of some of the marker genes has not been reported in HCC, their underlying effects as targets for HCC might be addressed in further biological and mechanistic investigations.
According to the GSVA, a series of pathways associated with patient outcomes were uncovered. Some cancer-associated pathways (e.g., the IL-6 signal pathway and the adenylate cyclase activating GPCR signaling pathway) were identified in both tumor and peritumor with varying degrees of enrichment. Moreover, several pathways with clinical potential were identified. For example, high expression of gene sets from the 3 0 -phosphoadenosine 5 0 -phosphosulfate (PAPS) biosynthetic process was observed in the peritumor of risk-high patients ( Figure 10). A previous family-based study in 240 families revealed that the 3 0 -phosphoadenosine 5 0 -phosphosulfate synthase 1 (PAPSS1) gene was a candidate HCC-susceptibility gene, which was associated with elevated serum a-fetoprotein and with poor survival in familial cases. Studies based on immunohistochemical analysis, western blot, and qRT-PCR revealed a significant change in the expression of PAPSS1 protein in patients under the treatment of IFN-a, and its relation to relapse and prognosis was also reported. Our study indicates that genes involved in this pathway (e.g., PAPSS1) may be applicable in clinical practice.
In this multicenter, retrospective cohort study, we developed and validated two weighted prognostic models (GES and IS) based on cancerization-field profile. Our study confirmed that field cancerization profile-based models could better predict the patient's prognosis, which offers insights for works on ll OPEN ACCESS iScience 25, 103747, February 18, 2022 15 iScience Article prognostic model construction. The nomogram constructed by combining GES, IS, and TNM stage has been proved effective in improving the individualized prediction of the overall risk (e.g., recurrence and metastasis) of patients. According to the peritumoral immune cell infiltration features, CD8 T cells in peritumor tissues showed a stronger correlation with the patient's prognosis than that in the tumor, suggesting that T cells provided by peritumor tissues may be a key factor explaining the diversified prognosis and clinical benefit of immunotherapy. Moreover, a series of prognosis-associated genes and pathways were identified in peritumor tissue, which may contribute to the study on the molecular mechanism of HCC as well as the development of therapeutic targets. With the prevailing of these prognostic signatures applied in the further classification of patients with HCC, they hold promise to reflect different biologic backgrounds with potential implications in patient selection for therapies and prediction of clinical outcomes.

Limitations of the study
We acknowledge several limitations that could be addressed in our study. First, as a retrospective study based on publicly available data, only a limited number of qualified studies were collected. As the role of peritumor has attracted a rising interest in carcinogenesis, further validation on a greater scale would be expected. In addition, the signatures were constructed based on a limited number of features, which include only 12,749 protein-coding genes, and the immune cell component involved did not represent all the HCC-associated immune features. As a growing number of researches have revealed the carcinogenesis-associated role of various features (e.g., peptides, long noncoding RNA, microRNA, and circle RNA), the inclusion of a wider range of features would not only enable the construction of a more robust model but also bring us closer to the underlying mechanisms of HCC. Finally, the mechanisms of the signature genes have not been clearly identified here, and experimental studies on these genes may provide important information to facilitate our understanding of their functional roles.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:  (Colaprico et al., 2016), from which 49 patients with paired expression profile were assigned to the LIHC cohort. Additionally, we obtained the gene expression data from the LIRI-JP project in ICGC (https://icgc.org/). Of the 232 patients included in the LIRI-JP project, paired expression profile was available in 199 patients, which were assigned to the LIRI cohort. For each cohort, the clinical data of patients were obtained for further analysis. As our analyses were based on deidentified data, institutional review board approval and informed consent were not required.
In silico decomposition of the immune microenvironment We adopted CIBERSORTx, a gene expression-based deconvolution approach to estimate the relative immune cell infiltration levels of tumor and peritumor of each patient. According to the CIBERSORTx algorithm, Monte Carlo sampling was applied to derive the deconvoluted p value for each sample (Newman et al., 2019). The proportions of 22 immune cells in both HCC tumor samples and peritumor samples were quantified using the LM22 signature (consisting of 547 marker genes) under 1,000 permutations. Significant cases (p < 0.05) which indicate accurate deconvolution were kept for further analysis. The fraction distribution of each immune cell in tumor and peritumor was compared using Wilcoxon rank-sum test. Spearman correlation test was applied to study the pairwise correlation in the 22 immune cells of tumor samples and peritumor samples, respectively.
We further evaluated the association between immune cell infiltration and patient prognosis. Immune cells significantly associated with patient survival (p < 0.05) were selected by univariate cox regression. For each cell, we select infiltration level-associated gene sets for GO enrichment analysis using R package cluster-Profiler (Yu et al., 2012).

Cancerization-field level GES prognosis signature
Unlike the previously published prognosis model constructed from tumor-only characteristics, we adopted the expression profile of both tumor and peritumor into model construction. A two-step procedure was used to select optimal marker genes for distinguishing risk-low patients from risk-high patients. In the first step, a univariate Cox proportional hazards model was applied on each of the profiled 12,749 genes in GSE14520 cohort to explore the association between gene expression and patient survival, significant genes (p < 0.05) were selected as candidates to allow more efficient feature selection. The number of significant genes was 1,876 and 1,940 in tumor and peritumor, respectively.
In the second step, candidate genes were subjected to further feature selection using the least absolute shrinkage and selection operator (LASSO) regression. LASSO uses L1 regularization to make the variables in the Cox regression collapsible. By changing the value of regularized parameter l, we are able to control the number of the final model (Goeman, 2010). By this method, we further screen out the most useful prognostic markers among the candidate genes. The R package glmnet was applied to perform the LASSO Cox regression model analysis (Friedman et al., 2010). Finally, a gene expression score (GES) model was constructed and normalized to predict the prognosis of HCC patients, the risk score for each patient was calculated based on the regression coefficient of each gene in the signature.
To make GES comparable between different profiling platforms, e.g., the signal intensity of microarray expression data and the FPKM value of RNAseq expression data, GES was calculated using the following formula: where n is the number of genes in the model, x i is the expression of the gene, and w i is the gene-specific weight, i.e., the coefficient for each gene.

Cancerization-field level IS prognosis signature
For the construction of a cancerization-field level immune score (IS) prognosis signature, a similar two-step procedure as GES was applied. Briefly, in the first step, univariate cox regression was performed on each of the 44 infiltrated immune cells ( iScience Article associated patient survival (p < 0.05) were selected as candidates. In the second step, the infiltration level of the candidate immune cells was subjected to LASSO Cox regression. Finally, an immune score (IS) model was constructed for the prediction of patient prognosis. IS was calculated using the following formula: where n is the number of immune cells in the model, x i is the infiltration level of the immune cell, and w i is the coefficient for each immune cell.

Assessment of the predictive value of the prognostic models
To figure out whether models constructed from cancerization-field profile (GES and IS) were more accurate on the prediction of patient survival, additional models were constructed, and their capability of risk stratification was compared.
1) For comparison with GES, four models were constructed, including tumor gene score 1 (TGS1), tumor gene score 2 (TGS2), peritumor gene score 1 (PGS1), and peritumor gene score 2 (PGS2). Model TGS1 is the tumoral part of GES, consisting of the tumor genes from GES with the same coefficient. Model TGS2 was built from tumor expression profile, with marker genes selected using the same two-step procedure as GES. For model TGS2, regularized parameter l was adjusted to keep the number of marker genes identical or similar to GES, ensuring that the prognosis prediction power did not result from a different gene number. Similar to that of TGS1 and TGS2, PGS1 and PGS2 were constructed using peritumor expression profile.
2) For comparison with IS, we constructed four models, tumor immune score 1 (TIS1), tumor immune score 1 (TIS1), tumor immune score 2 (TIS2), peritumor immune score 1 (PIS1), and peritumor immune score 2 (PIS2). Model TIS1 is the tumoral part of IS, consisting of the tumor immune cells from IS with the same coefficient. Model TIS2 was built by multivariate Cox regression on the profile of tumor immune cells. Similar to that of TIS1 and TIS2, PIS1 and PIS2 were constructed using peritumor immune cell infiltration profile.
To stratify patients of different risk scores into risk-high and risk-low groups, optimal cutpoint for each model was detected using the maximally selected rank statistics implemented in the R package maxstat (Hothorn and Lausen, 2003). The standardized log-rank statistic across different risk scores was calculated, thus providing a value of a cutpoint that corresponds to the most significant relation with outcome. HCC patients were then dichotomized into a risk-high and risk-low group according to the optimal cutpoint. Each model with its corresponding cutpoint was then validated in validation cohorts, with the survival of the high-risk and low-risk group estimated by the Kaplan-Meier method, and compared using the logrank test.
We further designed a random sampling approach to validate the superiority of cancerization-field profilederived models. For gene expression profile, 75% of patients in the GSE14520 cohort (156 individuals) were randomly sampled for 1,000 trials, which makes up 1,000 sets of training sets. For each training set, five prognostic models (GES, TGS1, TGS2, PGS1, and PGS2) were constructed, with corresponding optimal cutpoint determined based on the aforementioned methods. Each model was validated in LIHC and LIRI cohort, and log-rank p values were compared across five models using the pairwise Wilcoxon rank-sum test.
In an effort to compare GES and previously published models, the performance of prognosis prediction over time was compared. We used R package pec to perform the inverse probability of censoring weighting (IPCW) estimation of time-dependent Brier score based on ten-fold cross-validation, as well the computation of prediction error curves of each model (Mogensen et al., 2012).
We performed Cox multivariate regression to check whether GES and IS signature were independent of conventional clinicopathological factors (age, gender, TNM stage, HBV, ALT, cirrhosis, AFP, and tumor size). Furthermore, receiver operating characteristic (ROC) analysis was used in validating the GES and IS signature of their performance in predicting 24-and 36-month OS. Using R package forestplot and sur-vivalROC, the results were visualized by forest plots and ROC curves, respectively (Heagerty et al., 2000).

Development of nomogram and clinical use validation
We formulated a nomogram using the R package rms, which included GES classifier, IS classifier, and TNM stage in the nomogram (Nú ñ ez et al., 2011). We use calibration curves to assess the performance of the nomogram, which use the bootstrap method with 1000 resamples to get bias-corrected estimates of predicted versus Observed values, with a concordance index (C-index) calculated. Decision curve analysis was used to assess the clinical practicability of the nomogram, which applied the net benefits of threshold probabilities to quantitatively evaluate the clinical value of the nomogram (Vickers and Elkin, 2006). We further adopted the receiver operating characteristic (ROC) analysis to compare the predictive accuracy of TNM stage, GES classifier, IS classifier, and the nomogram on each of the three cohorts in our study.

Functional relevance exploration of diagnostic gene sets
To explore potential mechanisms underlying groups of distinct risk scores, patients with consistent high-or low-risk in both GES signature and IS signature were selected to analyze their enrichment for gene sets and pathways. The enrichment analysis was carried out by R package gsva, using gene sets from KEGG and GO (Ashburner et al., 2000;Hanzelmann et al., 2013;Kanehisa and Goto, 2000;Kanehisa et al., 2016). The gene sets were obtained from MSigDB database version 7.0 (Subramanian et al., 2005). Gene set enrichment scores (ES) of individual samples were calculated from the gene expression data using a Kolmogorov-Smirnov (KS) like random walk statistic with classical maximum deviation method (Subramanian et al., 2005). The GSVA ES of selected samples were then applied to fit the same linear model using R package limma and moderated t-statistics were estimated for each gene set (Ritchie et al., 2015). Significant gene sets were filtered with FDR adjusted p value < 0.05, and the top 40 gene sets with the highest fold change between groups were then visualized by heatmap and volcano plot, respectively.

QUANTIFICATION AND STATISTICAL ANALYSIS
Statistical analyses were performed using R (version 4.0.5). Kaplan-Meier analysis or multivariate cox analysis was performed for survival analysis. The spearman coefficient was applied to analyze correlations. Wilcoxon rank-sum test was used to compare the difference between the value from two groups. All reported p Values are two-tailed, and for all analyses, p < 0.05 was considered statistically significant. FDR correction was used for multiple test correction.