Integrating Single-cell RNA-seq to construct a Neutrophil prognostic model for predicting immune responses in non-small cell lung cancer

Abstract

Non-small cell lung cancer (NSCLC) is the most widely distributed tumor in the world, and its immunotherapy is not practical. Neutrophil is one of a tumor's most abundant immune cell groups. This research aimed to investigate the complex communication network in the immune microenvironment (TIME) of NSCLC tumors to clarify the interaction between immune cells and tumors and establish a prognostic risk model that can predict immune response and prognosis of patients by analyzing the characteristics of Neutrophil differentiation. Integrated Single-cell RNA sequencing (scRNA-seq) data from NSCLC samples and Bulk RNA-seq were used for analysis. Twenty-eight main cell clusters were identified, and their interactions were clarified. Next, four subsets of Neutrophils with different differentiation states were found, closely related to immune regulation and metabolic pathways. Based on the ratio of four housekeeping genes (ACTB, GAPDH, TFRC, TUBB), six Neutrophil differentiation-related genes (NDRGs) prognostic risk models, including MS4A7, CXCR2, CSRNP1, RETN, CD177, and LUCAT1, were constructed by Elastic Net and Multivariate Cox regression, and patients' total survival time and immunotherapy response were successfully predicted and validated in three large cohorts. Finally, the causes of the unfavorable prognosis of NSCLC caused by six prognostic genes were explored, and the small molecular compounds targeted at the anti-tumor effect of prognostic genes were screened. This study clarifies the TIME regulation network in NSCLC and emphasizes the critical role of NDRGs in predicting the prognosis of patients with NSCLC and their potential response to immunotherapy, thus providing a promising therapeutic target for NSCLC.

Keywords: Immunotherapy response; NSCLC; Neutrophil; Prognosis; Tumor immune microenvironment; scRNA-seq.

Fig. 1

Single-cell analysis: the cell clusters and their Marker were obtained by reduced-dimensional and clustering. Twenty-eight cell clusters (A) were obtained after the first level classification, and ten cell types (B) were identified by marker gene annotation. Fourteen cell clusters (D) were obtained after the second-level classification of lymphoid immune cells, and seven cell types (E) were identified by marker gene annotation. Fifteen cell clusters (F) were obtained after the second-level classification of myeloid immune cells, and seven cell types (G) were identified by marker gene annotation. Eighteen cell clusters (H) were obtained from normal epithelial cells after secondary classification, and then nine cell types (I) were identified by marker gene annotation. (C) Heatmap of the expression level of Marker genes from twenty-eight cell types

Fig. 2

The abundance of 28 cell types in Bulk RNA-seq and cell interaction network in scRNA-seq. A Twenty-eight cell types were annotated to the TCGA queue by CIBERSORT. (*P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001). B The cellular communication network in which Cancer cells interact with other cells

Fig. 3

Pseudotime analysis of Neutrophils and mutational analysis of NDRGs. According to the pseudotime (A) of Neutrophils, the cell population was divided into four different differentiation states (B), and NDRGs (G) were obtained by difference analysis of differentiation states. GSVA (KEGG terms) analyzes four different differentiation states (C–F). Top 30% mutation frequency of NDRGs and mutation type (H) and mutation status of NDRGs in different states (I)

Fig. 4

Construction and verification of the prognostic risk model. The intersection (B) of the differential gene (A) and NDRGs. Eight NDRGs with prognostic characteristics were screened by the Elastic Net Regression algorithm (C, D), and six prognostic risk model genes were confirmed by Multivariate Cox (E). The risk score distribution, patient status, mRNA expression heatmap, ROC curve, and KM survival curve of the training sets (F), the internal validation set (G), and the external validation sets (H). (I) Nomogram of the prognostic risk model. (J) The nomogram calibration curves to predict the 1-, 3-, and 5-year survival

Fig. 5

Immune predictive performance and clinical predictive power of the prognostic model. A After grouping the risk score according to the median, check the abundance of 28 immune cells in the high-risk and low-risk groups. B Spearman correlation analysis between risk score and the abundance of 28 kinds of immune cells. C–H Age, Gender, M stage, N stage, T stage, and Stage distribution of the patients in the high-risk and low-risk groups. Univariate Cox Regression (I) and Multivariate Cox Regression (J) analysis of clinical information of TCGA cohorts. (*P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001)

Fig. 6

Expression levels, survival analysis and functional studies of six prognostic genes in the TCGA cohort. A Expression levels of six prognostic genes in the TCGA cohort. B–G KM survival curves of six prognostic genes in the TCGA cohort. After grouping MS4A7 (H), CXCR2 (I), LUCAT1 (J), CD177 (K), CSRNP1 (L) and RETN (M) at high and low levels, the enriched KEGG and GO pathways were scored for GSVA

Fig. 7

The docking results of proteins encoded by prognostic genes with small molecular compounds. The docking results of MS4A7 with Estradiol (A). The docking results of CXCR2 with Abrine (B). The docking results of RETN with Ionomycin (C). The docking results of CSRNP1 with Beclomethasone (D). The docking results of CD177 with XL147 (E)