Comprehensive analysis of scRNA-Seq and bulk RNA-Seq data reveals dynamic changes in tumor-associated neutrophils in the tumor microenvironment of hepatocellular carcinoma and leads to the establishment of a neutrophil-related prognostic model

Abstract

Background: Analysis of hepatocellular carcinoma (HCC) single-cell sequencing data was conducted to explore the role of tumor-associated neutrophils in the tumor microenvironment.

Methods: Analysis of single-cell sequencing data from 12 HCC tumor cores and five HCC paracancerous tissues identified cellular subpopulations and cellular marker genes. The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) databases were used to establish and validate prognostic models. xCELL, TIMER, QUANTISEQ, CIBERSORT, and CIBERSORT-abs analyses were performed to explore immune cell infiltration. Finally, the pattern of tumor-associated neutrophil roles in tumor microenvironmental components was explored.

Results: A total of 271 marker genes for tumor-associated neutrophils were identified based on single-cell sequencing data. Prognostic models incorporating eight genes were established based on TCGA data. Immune cell infiltration differed between the high- and low-risk groups. The low-risk group benefited more from immunotherapy. Single-cell analysis indicated that tumor-associated neutrophils were able to influence macrophage, NK cell, and T-cell functions through the IL16, IFN-II, and SPP1 signaling pathways.

Conclusion: Tumor-associated neutrophils regulate immune functions by influencing macrophages and NK cells. Models incorporating tumor-associated neutrophil-related genes can be used to predict patient prognosis and immunotherapy responses.

Keywords: CellChat; HCC; Immunotherapy; Single-cell; Tumor-associated neutrophils.

Fig. 1

Single-cell sequencing analysis. A UMAP plot colored by various cell clusters. B UMAP plot colored by subpopulation of cells after annotation. C The cell types identified by marker genes. D Heatmap showing the top ten marker genes in each cell cluster

Fig. 2

Tumor immune microenvironment analysis. A–D Immune cell and risk model correlation analysis. Tumor microenvironment score differences compared based on ESTIMATE algorithm, E stromal score, F immune score, G estimate score, and H tumorpurity score. I Results of the ssGSEA analysis

Fig. 3

Prediction of response to immunotherapy. A, B, C, D Results of the TIDE predictive analysis. E Results of correlation analysis of immune checkpoint inhibitor genes. F Results of immune checkpoint inhibition gene expression analysis

Fig. 4

A–J KEGG Pathway AUCell Score Results. K Results of differential analysis of KEGG pathway tumor tissues and paracancerous tissues

Fig. 5

Analysis of glycolysis-related genes. A, B Expression correlation analysis of model genes and glycolysis-related genes in normal and tumor tissues. C Differential expression of glycolysis-related genes in high- and low-risk groups

Fig. 6

Trajectory analysis of tumor-associated neutrophils. A–D Differentiation trajectory results for tumor-associated neutrophils. E Variation in the expression of model genes in tumor-associated neutrophil differentiation trajectories. F, G Glycolysis gluconeogenesis signaling pathway and B cell receptor signaling pathway differentially expressed gene along the pseudotime were hierarchically clustered into five subclusters

Fig. 7

A The number of interaction in cell–cell communication network. B The weight of interaction in cell–cell communication network. C, D, E Cell–cell communication interaction in SPP1, IL16, and INF-II signaling pathway. F, G, H Heatmap display of intercellular communication weights

Fig. 8

Analysis of gene expression. A–H Analysis of gene differential expression in TCGA database. I Verifying the expression of genes that constitute the risk model through RT-qPCR