Constructing a neural network model based on tumor-infiltrating lymphocytes (TILs) to predict the survival of hepatocellular carcinoma patients

Clinical cohorts

We retrospectively evaluated independent cases from the following two cohorts: (a) patients who underwent liver cancer resection surgery at Qingdao University Affiliated Hospital from 2013 to 2016, and (b) The Cancer Genome Atlas (TCGA). We presented pathological images from both cohorts in whole-slide image (WSI) format. The Qingdao University Affiliated Hospital study cohort included 331 hepatocellular carcinoma (HCC) patients, corresponding to 331 digital pathological slides (WSIs). The study has been approved by the Ethics Committee of The Affiliated Hospital of Qingdao University (approval number QYFY WZLL 28523). After excluding poor-quality images, we ultimately retained 285 high-quality digital pathological slides for subsequent in-depth analysis and research. We selected diagnostic slides for download from the TCGA cohort, which included 378 HCC patients (https://portal.gdc.cancer.gov/analysis_page?app=Downloads, accessed on November 27, 2022).

In the cohort from (a) Qingdao University Affiliated Hospital, patients’ risk characteristics (including test results) were obtained from admission to surgery, and pathological report data came from professional pathology tests.

Ethics statement

The studies involving human participants were reviewed and approved by The Affiliated Hospital of Qingdao University Research Ethics Committee approved the study protocol. The patients/participants provided their written informed consent to participate in this study.

Digital scanning

For each patient’s diagnostic slide from the Qingdao University Affiliated Hospital study cohort, we selected three diagnostic slides per patient from a nine-point sample based on site. To ensure the quality and representativeness of the slides, we leveraged the expertise of two experienced pathologists. During the selection process, our criteria were as follows:

1. The sections are clear and without blurriness.

2. The tissue structure is clear, with good staining of cell nuclei and cytoplasm.

3. Abundant and well-defined cancer nests with clear layers.

After the selection was completed, we utilised the NanoZoomer-XR C12000 Hamamatsu scanner at a 20× magnification (0.5 µm/pixel) to ensure the clarity and resolution of the images.

Training cell classifier for cell annotation

In our study, we utilised the open-source software QuPath, which integrates machine learning algorithms, to conduct digital analysis on whole-slide images (WSI). The selection of regions of interest (ROIs) was based on the positional characteristics of the tumour and was performed by two independent pathologists to ensure consistency in selection. Ultimately, we took the intersection of the ROIs selected by both experts, ensuring that the chosen ROIs covered at least 50% of the tumour area while excluding the fibrous capsule surrounding the tumor. Before conducting the research, we performed image normalisation and intensity correction. Subsequently, we used QuPath for cell identification within the ROIs (parameters are detailed in the Supplemental Material), followed by detailed annotation of the tumour area. When classifying all cells within the ROI, we categorised cells into four main groups: red for tumour cells, yellow for lymphocytes, blue for stromal cells, and black for other cells (Fig. 1B).

We trained an efficient cell classifier (with eight hidden layers (maximum number of iterations: 1,000)) by manually annotating cells, enabling it to automatically classify all cells within the ROI. To achieve more accurate cell classification, we conducted multiple rounds of training (Fig. 1B). Finally, after independent evaluation by two pathologists, the QuPath-trained cell classifier was considered sufficiently accurate for further analysis and research when it reached an accuracy rate of 95%. However, there were variations in staining effects and tumour differentiation levels across different cancer nests. In cases where there were significant differences in staining effects and tumour differentiation levels among different cancer nests, and if the classifier’s accuracy did not meet expectations during training, model performance could be optimised by manually correcting some misclassified areas. This process can be considered as one cycle of model iteration.

Quantification of TILs indices

We adopted four currently recognized indices to quantify the presence of tumor-infiltrating lymphocytes (TILs):

eTILs% (Percentage of TILs in Tumor Cells):

• Definition: Calculate the proportion of TILs in tumor cells.

• Formula: eTILs% = (TILs/(TILs + Tumor Cells)) × 100

esTILs% (Percentage of TILs in Stromal Cells):

• Definition: Calculate the proportion of TILs in stromal cells.

• Formula: esTILs% = (TILs/(TILs + Stromal Cells)) × 100

etTILs% (Percentage of TILs in Total Cells):

• Definition: Calculate the proportion of TILs in total cell count.

• Formula: etTILs% = (TILs/Total Cell Count) × 100

eaTILs (Density of TILs in Tumor Area):

• Definition: Calculate the infiltration density of TILs in the tumor area.

• Formula: eaTILs = TIL/Analyzed Tumor Area (mm²)

These indices, calculated with precise methods, help us comprehensively quantify the presence of TILs and reveal their distribution and density in the tumor microenvironment (Zheng et al., 2022).

Statistical analysis

We utilised the Neuralnet and NeuralNetTools packages in R, version 4.3.2 (R Core Team, 2023), to incorporate the quantified pathological slide features for constructing a neural network model. Through the model, we obtained the high- and low-risk groupings for the patients. Subsequently, we further utilised the IOBR and survival packages in R for subsequent statistical analyses, including immune infiltration, immune checkpoint prediction, immune scoring, as well as univariate and multivariate Cox regression analyses.

Immunohistochemistry (IHC)

We performed IHC staining on 51 paraffin-embedded liver cancer tissues collected from Qingdao University Affiliated Hospital between January 2015 and June 2015 to determine the expression level of TNFSF4 protein. The staining intensity (0 = negative, 1 = weak, 2 = moderate, or 3 = strong) and the percentage of strongly stained cells in the observed field (0 = 0%, 1 = 1–25%, 2 = 26–50%, 3 ≥ 51%) were used to calculate the IHC staining score by multiplying the intensity and percentage scores. A final score <4 was defined as low expression, and a score of 4–9 was defined as high expression.

Patient characteristics

This study included a total of 465 HCC patients from two cohorts (median [IQR] age range of 57 [49–65] years for Cohort A and 61 [51–70] years for Cohort B; gender ratio ranged from 241 males [84.6%] out of 285 patients in Cohort A to 117 males [65%] out of 180 patients in Cohort B). The majority of diagnoses were for stages II and III of the disease. The datasets from both cohorts had similar distributions of clinical and pathological characteristics, and Table 1 summarises the baseline characteristics of our study cohorts (Table 1).

Risk factor analysis

In the risk factor analysis, we first focused on the patient cohort from Qingdao University Affiliated Hospital. By tracking follow-up data until January 1, 2022, we obtained the survival time and survival status of the patients and conducted comprehensive statistics on numerous clinical data.

In this study, we used eTILs% (percentage of tumor-infiltrating lymphocytes in tumor cells) as a feature of TILs, and employed the Cox proportional hazards model to identify significant risk factors in survival analysis (P < 0.05) (Fig. 1C). These factors included eTILs%, satellite lesions (surgical record/pathology), tumor diameter, AST, microvascular invasion, lymphocytes, number of lesions, liver capsule, etc. (Table 2).

Through multivariate Cox regression analysis, we further revealed that eTILs%, satellite lesions, and ALB remained significant independent prognostic factors (Fig. 1D) (Table 3).

The protective effect of TILs

In our study, through multivariate Cox analysis, we observed that the features of TILs exhibited the most significant differences in prognosis, with the highest hazard ratio (HR) and the smallest p-value. This highlights the importance of TILs in the survival analysis of liver cancer patients. To further explore the clinical significance of TILs and verify their independence in prognosis, we conducted subgroup analyses with two other features (ALB (Figs. 2A–2C) and satellite lesions (surgical record/pathology) (Figs. 2D, 2E)), confirming the significant prognostic value of eTILs% in different subgroups. This further substantiates that high levels of TILs infiltration can significantly extend patient survival time.

Subsequently, we performed Kaplan-Meier (K-M) curve analyses based on the four characteristics of TILs, all of which showed p-values less than 0.01. This indicates that, whether considering the overall level of TILs or their different characteristics, high levels of TILs infiltration are significantly associated with better patient survival outcomes.

Neural network model

Based on the patient slide characteristics from the Qingdao University Affiliated Hospital study cohort, we utilized seven quantified features related to TILs, including four TILs features, the nuclear-cytoplasmic ratio of TILs, the presence of TILs-aggregated lymph nodes within cancer nests, and the infiltration of TILs within cancer nests. The model was constructed using the neuralnet package in R to build a multilayer perceptron (MLP) neural network, which includes one hidden layer with five neurones (hidden = 5), and the output layer corresponds to the risk stratification of patients (high risk and low risk) to predict their survival status (Fig. 3A).

The survival time endpoints were chosen at 1, 3, and 5 years according to mainstream standards (Li & Ma, 2021; Albano, Bilfinger & Nemesure, 2018). Following model validation, the ROC for the 3-year model reached 0.975, indicating the best performance. Therefore, we chose 3 years as the survival time endpoint for patients (Fig. 3B).

We then selected 378 patients from the TCGA training cohort, excluding those with follow-up times less than 3 years and without transcriptome data, resulting in 180 patients for the external validation set. We performed the same cell segmentation and annotation processes on the validation set, resulting in a final ROC value of 0.836, indicating the successful construction and performance of our model in external validation (Fig. 3C). This professional analysis emphasizes the model’s reliability and applicability, providing strong support for future applications of TILs in liver cancer prognosis research.

After completing the neural network model validation, we classified the patients in the TCGA validation cohort into high-risk and low-risk categories based on the model results, and subsequently performed Kaplan-Meier (K-M) survival curve analysis. The K-M curve showed significant results, with a p-value less than 0.01, further proving the model’s accuracy in risk stratification (Fig. 3D).

Immune analysis and checkpoint screening

We conducted an in-depth study on the gene expression differences between two groups of patients predicted by the model in the TCGA cohort, filtering out 44 differentially expressed genes that are highly expressed in the high-risk group according to the criteria of p < 0.05 and LogFC < −0.5 (Fig. 4A). After obtaining the set of differential genes, we performed gene functional enrichment analysis using tools such as the clusterProfiler package in R, which indicated that these genes play a key role in biological processes involving cell proliferation and cell division (Fig. 4B). We speculate that these genes may promote further division and proliferation of tumour cells, leading to adverse outcomes. Subsequently, we conducted immune scoring for both high- and low-risk groups, and the analysis showed significant differences in the low-risk group in aspects such as stromal signature, immune signature, immune cell subsets, immune signalling molecules, immune CYT, IPRES.1, HLA, TIL, cytotoxic cells, T cells, CD8 T cells, X13 T cell signature, CD8 T effector, nd inflamed T cells. GEP, IFNG signature, cytotoxic cells. Danaher, and exhausted. CD8.Danaher (Fig. 4C). Immune infiltration analysis revealed a higher presence of macrophages M2 in the low-risk group, suggesting that the immune response may be associated with the polarisation of macrophages M2 (Fig. 4D) (Cao et al., 2023).

Finally, we conducted an immune checkpoint analysis based on the differences in TILs using the aforementioned differential genes, screening for 46 immune checkpoints. Among these checkpoints, we identified TNFSF4 as a potential target that may benefit, offering possible assistance for immune therapy in liver cancer patients (Fig. 4E).

TNFSF4 in IHC and expression analysis

To further validate our bioinformatics findings that identified TNFSF4 as a potential therapeutic target, we selected pathological slides from 50 consecutive patients within 1 year from the Qingdao University Affiliated Hospital study cohort (Table 4). We re-sectioned these patients’ paraffin blocks and performed immunohistochemistry (IHC) staining using a TNFSF4 antibody (Fig. 5A).

Through IHC staining, we found that the expression levels of TNFSF4 in HCC patients were significantly associated with prognosis. There was a significant difference in prognosis between the high expression and low expression groups (p < 0.05) (Fig. 5B). This result suggests that TNFSF4 could serve as a potential biomarker for the prognosis of HCC patients.

To further confirm the expression pattern of TNFSF4 in HCC, we analyzed data from The Cancer Genome Atlas (TCGA) database. The mRNA expression of TNFSF4 was much higher in HCC tissues compared to liver tissues that did not have a tumor (Fig. 5C). Paired sample analysis further validated this finding, revealing significantly higher TNFSF4 expression in tumor tissues compared to adjacent non-tumor tissues (Fig. 5D).

These findings indicate that TNFSF4 is not only highly expressed in HCC tissues, but that its high expression is also associated with a poor patient prognosis.

Conclusion

This study looked back at liver cancer patients from two large research groups, using machine learning and neural network models to look at clinical data and digitally scanned pathological slides. We deeply investigated the relationship between the immune status of liver cancer patients and their prognosis. The study quantified TILs and other key factors in a precise manner, significantly improving the analysis’s efficiency and accuracy. Our successfully established neural network model not only reliably predicted patient prognosis, but also performed excellently in independent validation cohorts, achieving high ROC results and demonstrating the model’s robustness and generalizability. The study comprehensively examined the distribution of TILs across different cell types, covering various dimensions such as cell proportions and infiltration densities. In addition, the study looked closely at immune scores, immune infiltration, and immune checkpoints to find out how TILs are connected to the immune microenvironment in liver cancer and what that means for patients. The study used differential gene screening to find that TNFSF4 was highly expressed in the high-risk group. This suggests that TNFSF4 could be a personalized therapeutic target for immunotherapy and points the way for new treatment strategies in the future.