We examined patients who underwent surgical resection or liver transplantation between 2013 and 2019 at Beijing Hospital. Patients included in this study had to meet the following criteria: (a) be at least 18 years old; (b) have a pathological diagnosis of HCC; (c) not receive any preoperative treatment; (d) have no history of prior malignancy, autoimmune disease, or immune deficiency disease; and (e) provide well-preserved formalin-fixed paraffin-embedded (FFPE) slides with hematoxylin-eosin (H&E) staining. Patients with incomplete clinical information were excluded. Ultimately, 64 patients were included in the study.

To analyze the TME of HCC and validate the prognostic model, pathological images, clinical information, and RNA-sequencing data were retrieved from the TCGA-LIHC database via Genomic Data Commons (http://gdc.cancer.gov/). Data preprocessing was conducted to enhance the quality of data and ensure the reliability of further analysis. The gene expression RNAseq data were normalized, duplicated values and missing values were eliminated. Patients without complete survival data or pathological images were excluded. Finally, 198 patients were included.

Our prior study proposed a reliable and effective pathological images processing pipeline (14). Each image was digitally captured at 40× magnification and labeled as a region of interest (ROI), defined as the major malignant region, using the ImageScope annotation tool. Subsequently, we randomly sampled 20 patches within each ROI and calculated the number of tumor cells, lymphocytes and stromal cells within these patches.

To categorize tumors into distinct immune phenotypes, we initially determined the percentage of lymphocytes and the ratio of lymphocytes to tumor cells in each image. Subsequently, we conducted a test for normality to identify the parameter with the least dispersion, selected based on the interquartile range (IQR), for further analysis (15). Patients were then stratified into two subtypes based on immune infiltration levels: high and low. This categorization was determined using the median lymphocyte percentage as the threshold. Finally, we compared clinical features and prognosis between these two subtypes.

We identified differentially expressed genes (DEGs) among various subtypes using the DESeq2 R package, employing criteria of a base mean > 10, log2 Fold Change > 1, and adjusted P value < 0.05 (17). Subsequently, we conducted Gene Ontology (GO) functional pathway enrichment analysis using the clusterProfiler R package, with significance determined at a P value < 0.05 (18). Furthermore, we obtained HALLMARK- and KEGG-related gene datasets from the Gene Set Enrichment Analysis (GSEA) official website. We then performed GSEA utilizing the GSEA algorithm (19) and Gene Set Variation Analysis (GSVA) employing the GSVA R package (20).

The ESTIMATE (Estimation of STromal and Immune cells in Malignant Tumor tissues using expression data) analysis was performed to assess the level of immune infiltration, utilizing the estimate R package (21). Additionally, Cell type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) analysis was employed to determine the relative abundance of 22 different immune cell types within the tumor tissue (22). Furthermore, Tumor Immune Dysfunction and Exclusion (TIDE) analysis was carried out to evaluate the potential for tumor immune escape, utilizing the TIDE website (http://tide.dfci.harvard.edu/) (23).

To further investigate the prognostic value of immune classification, we categorized patients into two groups: a good prognosis group and a poor prognosis group, defined as having a DFS > 1 year (16). From the Beijing Hospital cohort, we collected 55 variables comprising clinical and pathological features. We then developed Receiver Operating Characteristic (ROC) curves for each variable using the pROC R package and extracted the Area Under the Curve (AUC) for evaluation (24). Variables with an AUC exceeding 0.6 were selected, and Least Absolute Shrinkage and Selection Operator (LASSO) analysis was employed to reduce the number of variables in the risk model using the Glmnet R package (25). A 20-fold cross-validation was conducted to identify the optimal lambda value. “Lambda.1se” was utilized to determine the minimum number of independent variables required for a well-performing model. Subsequently, we employed the Beijing Hospital cohort as the training set and 58 patients from the TCGA-LIHC database, who provided complete clinical information, as the test set. Logistic regression was then applied to establish a prognostic model, and ROC curves were generated for both the training and test sets. Finally, a nomogram was constructed based on the prognostic model.

OS was defined as the period between the day of pathological diagnosis and the day of death, while DFS was defined as the duration between the day of pathological diagnosis and the occurrence of tumor recurrence, metastasis, or death. Patients who remained free of recurrence were censored at the final follow-up. Survival analysis was conducted using Kaplan–Meier (K-M) analysis employing the Survival and Survminer R packages. Categorical and non-normally distributed measurement variables were compared using the Wilcoxon test, whereas normally distributed measurement variables were compared using the t-test. All statistical analyses were performed using R software (version 4.1.3).

Our prior research had developed an adapted CNN model to recognize cells autonomously and accurately in pathological images of patients with HCC, with classification accuracies of 95.7%, 92.3%, and 77.6% for tumor cells, lymphocytes, and stromal cells, respectively (14).

The adapted CNN model was utilized in this study, applied to both the Beijing Hospital and the TCGA cohort. Analysis of cell type percentages revealed no significant disparities between the two cohorts (Figure 1A). Lymphocyts, stromal cells and tumor cells account for 6.26%, 38.76% and 70.52% in the Beijing Hospital cohort, and 5.56%, 37.70% and 66.44% in the TCGA cohort respectively.

Lymphocyte percentage and lymphocyte/tumor cell ratio were computed as potential parameters for further analyses, the parameter with a lower degree of dispersion serves as the basis for subsequent grouping. Both parameters constituted non-normally distributed continuous variables, with the dispersion of lymphocyte percentage being less pronounced (IQR 0.07 vs. 0.14 in the Beijing cohort, 0.04 vs. 0.09 in the TCGA cohort, as depicted in Figure 1B), so that the lymphocyte percentage was selected as the parameter for stratification (26). Images with lymphocyte percentages above or below the median were categorized as having high or low immune infiltration (median = 0.039 in the Beijing cohort, median = 0.011 in the TCGA cohort). Lymphocyte percentages falling 1.5 times below Q1 or exceeding 1.5 times above Q3 were identified as outliers, 2 outliers in Beijing cohort and 14 outliers in TCGA cohort were excluded.

We categorized the patients into high- and low-immune cell infiltration subtypes based on the threshold described above. For patients in the Beijing Hospital cohort, we collected data on 17 parameters, including epidemiological factors, indicators of liver function, medical history, tumor stage, and pathohistological features of the tumor, and compared the two subtypes (refer to Table 1). The low-infiltration subtype exhibited a higher incidence of hypertension (51.4% vs. 17.2%, p = 0.010) and fatty liver (25.7% vs. 3.4%, p = 0.017), and displayed a larger tumor diameter (median 5.0 vs 4.5, p = 0.048). Additionally, the low-infiltration subtype demonstrated a higher incidence of satellite nodules, elevated AFP levels, and a more advanced TNM stage, but showed no statistical significance. These findings suggest that lower immune infiltration may be associated with a history of metabolic syndrome and may promote tumor progression.

We further performed K-M analysis on both the Beijing Hospital and TCGA cohorts to assess the prognostic value of immune classification. In the Beijing Hospital cohort, patients with the high-infiltration subtype exhibited a favorable prognosis in terms of DFS (p=0.013), but no significant difference was observed in OS (Figures 2A, B). In the TCGA cohort, the high infiltration subtype demonstrated a favorable prognosis in both OS and DFS (Figures 2C, D, p=0.012 for OS, p=0.026 for DFS).

RNA sequencing data were gathered from the TCGA-LIHC database. We identified the DEGs between the two subtypes using the DESeq2 R package and annotated genes associated with immune pathways according to the KEGG database (Figures 3A, B). The analysis revealed that most of the immune-related genes were down-regulated in the low infiltration subtype. To delve deeper into the discrepancies in cellular function between the subtypes, we performed functional enrichment analyses utilizing the GO, GSEA, and GSVA methodologies (Figures 3C, D, E). The top 10 pathways enriched in the GO analysis (sorted by qvaule, increased) were all linked to immune function. Meanwhile, the top two pathways enriched in the GSEA analysis (ranked by absolute NES, decreased) were the chemokine and cytokine signaling pathways. The extent of immune infiltration was quantified using ESTIMATE analysis, and the estimated immune and stromal scores were compared between the two subtypes using the Wilcoxon test (Figure 3F). The high-infiltration subtype exhibited higher scores, indicating a greater degree of immune cell infiltration in the TME. Subsequently, CIBERSORT analysis was performed to assess immune cell abundance in the two subtypes (Figure 3G). The findings revealed that the low-infiltration subtype manifested a higher level of type 2 macrophages (M2), monocytes, and resting natural killer (NK) cells, suggesting a propensity towards immune suppression. Finally, TIDE analysis was employed to evaluate the potential for tumor escape, indicating no significant difference between the two subtypes (Figure 3H).

Patients with a DFS shorter or longer than 1 year were classified into poor or good prognosis groups. A total of 55 variables, including clinical and pathological features, were collected from the patients in the Beijing Hospital cohort. ROC curves were developed for each variable to evaluate their predictive value, and variables with an AUC greater than 0.6 were listed (Figure 4A). The immune classification ranked 8 (AUC = 0.67), while other variables with high rankings were mostly tumor features, such as TNM stage, microvascular invasion (MVI), and serum AFP level, which were consistent with previous findings (7, 8). To develop a prognostic model for patient outcomes, we used all patients from the Beijing Hospital cohort as the training set and 58 patients from the TCGA cohort with complete clinical information as the test set. We performed LASSO analysis and cross-validation to reduce the number of variables and determine the minimum number of variables needed for a model with favorable performance (Figure 4B). Finally, five variables were included in the logistic regression analysis: immune classification, age, AFP level, vascular invasion, and TNM stage, and a nomogram was developed (Figure 4C). ROC curves were generated for both the training and test sets, with AUCs of 0.918 and 0.814, respectively (Figure 4D).