Tumor niche network-defined subtypes predict immunotherapy response of esophageal squamous cell cancer

Abstract

Despite the promising outcomes of immune checkpoint inhibitors (ICIs), resistance to ICI presents a new challenge. Therefore, selecting patients for specific ICI applications is crucial for maximizing therapeutic efficacy. Herein, we curated 69 human esophageal squamous cell cancer (ESCC) patients' tumor microenvironment (TME) single-cell transcriptomic datasets to subtype ESCC. Integrative analyses of the cellular network and transcriptional signatures of T cells and myeloid cells define distinct ESCC subtypes characterized by T cell exhaustion, and interleukin (IL) and interferon (IFN) signaling. Furthermore, this approach classifies ESCC patients into ICI responders and non-responders, as validated by whole tumor transcriptomes and liquid biopsy-based single-cell transcriptomes of anti-PD-1 ICI responders and non-responders. Our study stratifies ESCC patients based on TME transcriptional network, providing novel insights into tumor niche remodeling and potentially predicting ICI responses in ESCC patients.

Keywords: Cancer; Cancer systems biology; Immune response; Transcriptomics.

Graphical abstract

Figure 1

Schematic workflow for transcriptomic analysis of TME from ESCC patients

Figure 2

Immune cells analysis and classification (A) Uniform Manifold Approximation and Projection (UMAP) display of whole cells from 69 patients. Single-cell RNA-sequencing (scRNA-seq) results of the cells of TME were integrated and projected. (B and C) Non-tumor cells were isolated, and UMAP was redrawn with individual patient information (B) and five major cell types (C). (D) UMAP display of T cells subgroup with unique patients ID. T cells were isolated from immune cells and clustered again. (E) T cells were classified into four sub-groups by principal component analysis (PCA) and Pearson correlation. PCA result was clustered by the dendrogram, and Pearson correlation was displayed by color spectrum. (F) T cells were displayed in UMAP based on the sub-groups defined from PCA and Pearson correlation. (G) Each sub-groups of T cells were shown with subsets using stacked bar plots. (H) Myeloid cells of each patient were displayed with UMAP. Myeloid cells were isolated from immune cells and clustered independently. (I) Myeloid cells were categorized into four sub-groups by PCA and Pearson correlation. PCA results were displayed with a dendrogram, and Pearson correlation was shown by color spectrum. (J) Myeloid cells were displayed with sub-groups identified from PCA and Pearson correlation. (K) Each myeloid cell sub-group was displayed with subsets of myeloid cells.

Figure 3

Comparative analysis of patients by myeloid and T cell classifications (A) T_ex cell markers expression in each T cell sub-group. (B) The number and proportion of patients of T cell-sub-groups were displayed in each patient’s sub-groups categorized by myeloid cells. (C) The number of sub-grouped patients by myeloid cell was displayed in each patient’s sub-groups categorized by T cells. The proportion of myeloid-cell-based classified patients in each sub-group of T cells was shown with pie plots. (D) Individual patients were subjected to each sub-group of myeloid and T cells by Sankey plot. P009A patient was not included in the myeloid cell-based sub-group due to the lack of myeloid cells in the dataset. Each patient was classified into 12 groups (MT groups) by sub-groups of myeloid cells and T cells and categories. (E) T cells of each patient from 12 groups were displayed with UMAP. (F) T_ex cell markers expression in T cells in the MyT groups of patients. (G) GSEA analysis was performed in T cells of MT groups of patients. The results of GSEA from the m1t4 and m3t4 groups of patients were compared. GOBP and REACTOME databases were used, and the significant signaling pathways with positive values of NES were compared. Overlapped signaling pathways were displayed with a Venn diagram. (H and I) Cell-to-cell interactions were inferred using the CellChat analysis package. Tumor cells and TME cells of m1t4, m2t3, and m3t4 grouped patients were analyzed, and significant IL1, IL16, IL10, and IFN-II signaling pathways-related interactions were visualized with circle plot (H), m1t4 and m3t4 grouped patients were compared with m2t3 grouped patients to identify specific signaling pathways (I).

Figure 4

Validation of predicted responders by anti-PD-1 immunotherapy-treated patients (A–D) Peripheral blood immune cells transcriptomes of three responders (PBMC-R) and three non-responders (PBMC-NR) (to anti-PD-1 ICI) were integrated and presented with UMAP by cell types (A, B), patients (C), and response groups (PBMC-R vs. PBMC-NR) (D). (E) T_ex marker genes expression was compared by the anti-PD-1 response (R vs. NR). (F–H) GSEA analysis performed in T cells of PBMC-R vs. PBMC-NR groups using the REACTOME database. Significant results with positive NES and negative NES were compared with the results of m1t4 (F) and m3t4 (G). Enrichment plots of PD-1 signaling and Interferon signaling were displayed (H). (I) Bulk RNA-seq datasets of responders (Tumor-R) and non-responders (Tumor-NR) were collected and analyzed. (J) Pathway signatures of Tumor-R were analyzed using enrichr through comparison with the Tumor-NR group. (K–N) Transcriptomes of Tumor-R and Tumor-NR were compared with single-cell transcriptomes of three MyT groups using the Scissor package. Each MyT group dataset was prepared by integrating tumor cells and TME cells (K). Cells of Tumor-R phenotypes were displayed in red, and Tumor-NR phenotypes were in blue. Cells not assigned to Tumor-R or Tumor-NR were annotated as neutral cells (L). Results of the Scissor analysis in three groups were compared using barplots with absolute cell numbers (M) and proportion (N). ∗∗∗∗p values < 0.0001, as determined by Fisher’s exact test.

Figure 5

Specific similarity of m3t4 group with ICI-responders (A) scRNA-seq datasets, including tumor cells and TME cells of 69 ESCC patients, were prepared to compare with ICI responders and non-responders. (B) 69 patients’ datasets, including tumor cells and immune cells, were converted to a single matrix categorized by MT groups to integrate with transcriptomes with bulk RNA-seq datasets (Tumor-R and -NR). (C) Transcriptomes of MT groups and Tumor-R and -NR groups were compared using PCA analysis and displayed with circular dendrogram. (D) single-cell transcriptomes of PBMC-R and PBMC-NR were integrated with 69 ESCC patients’ TME transcriptomes. (E) Correlation of MT patient groups with the PBMC-R group was analyzed using PCA and shown by the dendrogram. (F) GSEA results of tumor cells of m1t4 and m3t4 were compared with that of the Tumor-R group using Venn diagrams. (G) T cells of m1t4 and m3t4 were analyzed with GSEA (database: GOBP) and compared with GSEA results of Tumor-R and T cells of PBMC-R groups. Overlapped pathways of m3t4-Tumor-R and m3t4-PBMC-R were also highlighted. (H) GSEA (database: REACTOME) results from T cells of m1t4 and m3t4, Tumor-R, and T cells of PBMC-R. m1t4 exclusive pathways were highlighted. Overlapped pathways of m3t4-Tumor-R and m3t4-PBMC-R were also highlighted. (I) Pathway scores of FGF, IFN-gamma, and HSF1 in m1t4, m3t4, Tumor-R, Tumor-NR, PBMC-R, and PBMC-NR groups were shown with dotplots. (J) Gene expression levels of OAS3, SCCPDH, and LIMK1 in m1t4, m3t4, Tumor-R, Tumor-NR, PBMC-R, and PBMC-NR groups were shown with heatmaps.