Identification of a tissue resident memory CD8 T cell-related risk score signature for colorectal cancer, the association with TME landscapes and therapeutic responses

Abstract

Backgrounds: The tissue resident memory CD8 T cell (Trm) constitutes an important component of the local immunity. In the context of malignant tumors, mounting evidence also supports the potential anti-tumor property of this cell subset. Therefore, identification of Trm marker genes and exploration of the causative effect of Trm in shaping tumor microenvironment (TME) heterogeneity might provide novel insights for the comprehensive management of cancer patients. Methods: By dissecting a single T cell transcriptome dataset, we acquired marker genes for Trm, which were latter applied to bulk RNA sequencing profiles of two large colorectal cancer (CRC) patient cohorts downloaded from TCGA and GEO databases. First, colorectal cancer patients were divided into different Trm clusters using consensus clustering algorithm. Then, we established a Trm-related gene (TRMRG) risk score signature and tested its efficacy in predicting prognosis for colorectal cancer patients. Moreover, a sequence of rigorous and robust analyses were also carried out to investigate the potential role of Trm-related gene risk score in tumor microenvironment remodeling and therapeutic utility of it in colorectal cancer treatment. Results: A total of 49 Trm marker genes were identified by analyzing single cell RNA sequencing profiles. First, colorectal cancer patients were successfully classified into two Trm clusters with significant heterogeneity in functional enrichment patterns and tumor microenvironment landscapes. Then, we developed a Trm-related gene risk score signature and divided patients into different risk levels. High risk patients were characterized by attenuated immunogenicity, weakened sensitivity to immunotherapy, as well as adverse clinical outcomes. While low risk patients with advantages in survival exhibited increased immunogenicity, stronger metabolic activity and improved immunotherapeutic responses. Conclusion: Through combinatorial analysis of single cell and bulk RNA sequencing data, the present study identified Trm to play a non-negligible role in regulating the complexity and heterogeneity of tumor microenvironment for colorectal cancer. Moreover, the Trm-related gene risk score signature developed currently was corroborated to be tightly correlated with prognosis and therapeutic responses of colorectal cancer patients, thus exhibiting potential application value for clinical practice.

Keywords: colorectal cancer; immunotherapy; prognostic model; single cell RNA sequencing; tissue-resident memory T cells; tumor microenvironment.

FIGURE 1

Identification of Trm marker genes via scRNA-seq analysis (A,B) t-SNE plots colored by different (A) CRC samples and (B) cell clusters. (C) Identification of different T cell subsets. (D) Heatmap showing the top 5 marker genes of each T cell subtype. Tex: exhausted T cell; Tem: effector memory T cell; Temra: recently-activated effector memory T cell; MAIT: mucosal-associated invariant T cell; Trm: tissue-resident memory T cell.

FIGURE 2

Identification of Trm clusters in CRC patients. (A) Consensus matrix heatmap showing two Trm clusters. (B) CDF curve, (C) k-mean algorithm and (D) PCA analysis showing the robustness and stability of the clustering. (E) Heatmap showing the expression pattern of Trm marker genes in two clusters. (F) Bar plots showing the correlation of Trm clusters with clinicopathological factors including age, gender, tumor location, T, N, M status and tumor stage. CDF: Cumulative distribution function; PC: Principal component. Statistical Significance: *p < 0.05; ***p < 0.001; ns: not significant.

FIGURE 3

Correlation of Trm clusters with biological function and TME landscapes. (A) Heatmap showing the enrichment pattern of KEGG items for two Trm clusters. (B) Bar plot showing the enrichment pattern of Hallmark items for two Trm clusters. (C) Box plots comparing the infiltration differences of 28 immune cells between two Trm clusters. (D,E) Box plots comparing the differences in (D) TME score and (E) tumor purity between two Trm clusters. TME: Tumor microenvironment; KEGG: Kyoto Encyclopedia of Genes and Genomes. Statistical Significance: **p < 0.01; ***p < 0.001; ns, not significant.

FIGURE 4

Analysis of clinical correlation and biological function for TRMRG risk score signature. (A) Bar plots showing the frequencies of different clinical features including age (>65 vs. ≤65), gender (male vs. female), tumor location (left vs. right), T status (T1-T2 vs. T3-T4), N status (N0 vs. N1-N3), M status (M0 vs. M1) and tumor stage (I-II vs. III-IV) in high and low risk patients. (B) Box plots comparing risk score differences between multiple clinical subgroups. (C) GSEA analysis of TRMRG risk score signature. (D) Bar plot showing the enrichment pattern of Hallmark items for high and low risk patients. (E) Box plots comparing enrichment differences of gene sets developed by Mariathasan et al. Statistical Significance: *p < 0.05; **p < 0.01; ***p < 0.001; ns, not significant.

FIGURE 5

Analysis of TME landscape for TRMRG risk score signature. (A) Heatmap showing the immune cell infiltration level, TME score and tumor purity between high and low risk patients. (B–D) Box plot comparing differences in (B) immune cell infiltration level (C) TME scores and (D) tumor purity between high and low risk patients. (E,F) Pearson correlation analysis of the risk score with (E) immune cell infiltration levels, (F) TME scores and tumor purity. (G,H) Pearson correlation analysis of expressions of individual TRMRG with (G) immune cell infiltration levels, (H) TME scores and tumor purity. Statistical Significance: *p < 0.05; **p < 0.01; ***p < 0.001; ns: not significant.

FIGURE 6

Analysis of the activation status of seven-step anti-cancer immunity cycle for high and low risk patients. Statistical Significance: *p < 0.05; **p < 0.01; ***p < 0.001; ns, not significant.

FIGURE 7

Immunotherapeutic sensitivity analysis of TRMRG risk score signature. (A) Heatmap showing the expression patterns of HLA family genes, immunosuppressive molecules and ICGs between high and low risk patients. (B–D) Bar plots comparing expression differences of HLA family genes, immunosuppressive molecules and ICGs between high and low risk patients. (E) Differences in TCR richness and diversity between high and low risk patients. (F) Bar plots showing the frequencies of high, medium, low level IPS scores for anti-PD1 treatment, anti-CTLA4 treatment and the combinatorial therapy in high and low risk patients. (G) Box plot comparing differences in TIDE scores between high and low risk patients. (H) Frequencies of responders and non-responders to ICI therapies for patients with high and low risk. (I) Kaplan-Meier plot showing the OS for high and low risk patients after receiving anti-PD-L1 treatment in IMvigor210 cohort. (J) Box plot comparing the risk scores for responders and non-responsers in GSE78220. (K) Frequencies of responders and non-responders in high and low risk group in GSE78220. Statistical Significance: *p < 0.05; **p < 0.01; ***p < 0.001; ns: not significant.

FIGURE 8

Metabolism analysis of TRMRG risk score signature. Violin plots comparing the metabolic activities of (A) carbohydrate, (B) lipid, (C) amino acid and (D) drug between high and low risk patients. (E) Distinctions in metabolic activities between two metabolism subtypes. (F) Kaplan-Meier plot for patients with metabolism subtype A and B. (G) Risk score differences between patients with metabolism subtype A and B. (H) Frequencies of two metabolism subtypes in high and low risk patients. Statistical Significance: *p < 0.05; **p < 0.01; ***p < 0.001; ns, not significant.

FIGURE 9

Correlation of TRMRG risk score signature with sensitivity to anti-tumor drugs. (A) Differences in IC50 values for cisplatin, cytarabine, vinblastine and 5-fluorouracil between high and low risk patients. (B) Pearson correlation analysis of TRMRG expression levels with IC50 values for anti-tumor drugs in NCI-60 cell lines. Statistical Significance: *p < 0.05; ***p < 0.001.

FIGURE 10

Development of a nomogram based on TRMRG risk score signature. (A) Nomogram integrating TRMRG risk score, age, gender and tumor stage for predicting 5, 7, 10 years OS. (B–D) Time-dependent ROC curves evaluating the efficacy of the total nomogram score in predicting (B) 5, (C) 7, (D) 10 years OS, in comparison to age, gender and tumor stage. (E–G) DCA curves estimating the efficacy of the nomogram in predicting (E) 5, (F) 7, (G)10 years OS from the perspective of clinical benefit. (H–J) Calibration curves of the nomogram for predicting (H) 5, (I) 7, (J)10 years OS. AUC: area under curve.