Identification and analysis of a CD8+ T cell-related prognostic signature for colorectal cancer based on bulk RNA sequencing and scRNA sequencing data: A STROBE-compliant retrospective study.

Colorectal cancer (CRC) is one of the most common malignancies worldwide, leading to a large number of cancer-related mortalities. Aberrant CD8+ T cell infiltration plays a critical role in tumor progression and patient prognosis. This study aimed to identify a prognostic model for CRC based on CD8+ T cell-related genes. The infiltration levels of immune cells in CRC tissues were accessed by the ESTIMATE algorithm. Weighted gene co-expression network analysis (WGCNA) analysis was used to select CD8+ T cell-related genes. Prognostic genes were identified using Cox regression analysis and Kaplan-Meier curves. The least absolute shrinkage and selection operator (LASSO) algorithm was used to construct prognostic models. Gene set enrichment analysis (GSEA) was performed to annotate enriched gene sets. Single-cell RNA (scRNA) sequencing analysis was used to examine gene expression in different cell types. We found that the downregulated infiltration level of CD8+ T cells was an independent prognostic factor for CRC and selected a cluster of differentially expressed genes correlated with CD8+ T cell infiltration (CD8TDEGs). Subsequently, we identified 18 prognostic CD8TDEGs, according to which patients were reclassified into two clusters with distinct overall survival. Seven prognostic CD8TDEGs were selected to calculate the constructed prognostic model's risk scores. Interestingly, although CRC tissues with higher risk scores had higher infiltration levels of CD8+ T cells, the level of immune checkpoint genes was also high. Moreover, the scRNA-sequencing analysis showed that the expression levels of CD8TDEGs in the prognostic model varied among different types of cells. This study constructed a novel prognostic model for CRC and provided a foundation for targeting CD8+ T cell infiltration to improve the survival of CRC patients.

1. Introduction

Colorectal cancer (CRC) ranks the third and fifth leading cause of cancer-related mortality in the United States and China, respectively.[ Although therapeutic strategies advanced in recent decades, metastasis, recurrence, and drug resistance of CRC cells still threaten patient’s prognosis.[

Tumor cells live in a complex microenvironment consisting of multiple kinds of cells, such as stromal cells and immune cells.[ The infiltrated immune cells in the tumor microenvironment contain innate immune cells (macrophages, neutrophils, dendritic cells, mast cells, and NK cells) and adaptive immune cells (T cells and B cells).[ Tumor-infiltrating immune cells have been proved to participate in the multistep progression of cancer through breaking the balance between anti-and pro-tumorigenic signals.[

CD8+ T cells are a kind of lymphocytes that develop in the thymus and are committed to detecting antigens expressed by tumor cells.[ When the MHC class I molecules were cross presented to CD8+ T cells by dendritic cells, effector CD8+ T cells with cytotoxic capacity, namely CTLs, were generated.[ After being active, CD8+ T cells can bind and kill tumor cells via secreted death-inducing granules containing granzymes, perforin, cathepsin C, and granulysin.[ In addition, the interaction of Fas ligand expressed on CD8+ T cell membrane and Fas receptors on tumor cells leads to the fragmentation of tumorous DNA.[ The high infiltrating level of CD8+ T cells in most solid tumors have been proved to be a favorable prognosis biomarker for patients.[ However, cancer cells developed defending mechanisms. For example, CRC cells can express more immune checkpoint genes to inhibit CD8+ T cell function.[ Therefore, it is meaningful to construct a model for CRC patients to predict their prognosis, CD8+ T cell-infiltrated level, and expression of immune checkpoint genes.

In the present study, we identified a prognostic model for CRC based on CD8+ T cell-related genes and discovered that the poor prognosis of CRC patients with high-risk scores was mainly attributed to more abundant immune checkpoint genes in tumor tissues.

2. Material and Methods

2.1. Analysis of gene expression data

The RNA sequencing data and clinical characteristics of patients were downloaded from the GEO database (GSE44076)[ and the TCGA CRC cohort that was deposited in the University of California, Santa Cruz (UCSC) Xena browser (https://xenabrowser.net/datapages/). The differentially expressed genes (DEGs) were analyzed by R package “limma”. The cutoff was log fold change > 2 and P value < .05. The single-cell RNA (scRNA) sequencing data of CRC tissues were downloaded from the GEO dataset (GSE81861)[ and analyzed by R package “Seurat”.

2.2. Evaluating infiltrated immune cells in CRC tissues

To investigate the immune infiltration landscape in CRC tissues, the ESTIMATE algorithm was used to assess the immune infiltration level according to the expression levels of immune cell-specific marker genes. Marker genes of immune cells were obtained from the bulk transcriptome data of Bindea et al[ The prognosis values of infiltrated immune cells were analyzed by the TIMER online tool (https://cistrome.shinyapps.io/timer/).

2.3. Weighted gene co-expression network analysis (WGCNA)

Gene expression values in CRC tissue samples were used to construct the weighted co-expression network using the R package “WGCNA”. According to Pearson correlation value between genes, each transcript was first clarified into one similarity matrix which was subsequently transformed into an adjacency matrix. The adjacency matrix as calculated by amn = |cmn| β among which the cmn and amn mean Pearson correlation between paired genes and adjacency between paired genes, respectively. Parameter β was used to improve the correlation difference between genes. When the power of β = 4, a topological overlap matrix was generated from the adjacency matrix. To classify genes with similar expression patterns into distinct modules, a dynamic hybrid cutting method was applied by using a bottom-up algorithm whose module minimum size cutoff is 10.

2.4. Enrichment analysis and hub gene selection

GO and KEGG pathway enrichment analyses of genes were conducted in the DAVID database (https://david.ncifcrf.gov). The biological processes and KEGG pathway enrichment analysis of DEGs were visualized by using the Cytoscape and OmicShare tool (http://www.omicshare.com/tools).

2.5. Construction of the prognostic model

Cox regression analyses were performed to select prognostic DEGs correlated with CD8+ T cells. The prognostic model was constructed through least absolute shrinkage and selection operator (LASSO) regression based on the above DEGs using the R package “lasso”.

2.6. Gene set enrichment analysis (GSEA)

The GSEA analysis was used to search for gene sets that are correlated with the risk scores in CRC. The gene expression data were obtained from the TCGA CRC cohort. CRC tissue samples were clarified into a high and low group according to the risk score. The GSEA tool was implemented to explore the distribution of gene sets in the MSigDB database (http://www.gsea-msigdb.org/gsea/msigdb/index.jsp). The gene sets with a P value of < .05 and FDR value of < .25 were identified to be significantly correlated with the risk score.

2.7. Statistical analysis

Statistical analysis was performed by using GraphPad Prism 8.0 (GraphPad Software) using the Student t test. Survival analysis was performed using Kaplan–Meier curves with a log-rank test or Cox proportional hazard model. All data are presented as the mean ± standard deviation (SD). P < .05 was considered statistically significant.

3. Results

3.1. The infiltrated levels and prognostic values of immune cells in CRC

We first analyzed the infiltration levels of immune cells in CRC using the ESTIMATE algorithm, and the results showed that the infiltrated levels of immune cells were significantly lower in CRC tissues compared to normal tissues, including B cell, CD8+ T cell, CD4 + T cell, macrophage, and dendritic cell (Table 1). Besides, the multivariate Cox proportional hazard model exhibited that the infiltrated level of CD8+ T cell, rather than other immune cells, was an independent prognostic factor for CRC patients (Table 1).

Table 1

The relative infiltrated level and prognostic value of immune cells in CRC.

Cell	FC	P value	HR	95% CI	P value
B cell	0.76	<.001	2.02	0.022–184.35	.760
CD8+ T cell	0.91	<.001	0.010	0.00–0.47	.019
CD4+ T cell	0.91	.0014	0.397	0.0040–36.81	.689
Macrophage	0.79	<.001	36.861	0.45–2994.13	.108
Neutrophil	0.99	.87	0.083	0.00–110.24	.497
Dendritic cell	0.97	.048	5.634	0.30–104.46	.246

95%CI = 95% confidence interval, FC = fold change, HR = hazard ratio.

3.2. Identifying genes correlated with CD8+ T cell infiltration in CRC by using WGCNA

The genes deposited in the TCGA CRC cohort were applied to establish a co-expression network through the “WGCNA” package. To construct a scale-free network, β = 4 (scale-free R2 = 0.807) was identified as the soft-thresholding power (Fig. 1A and B). A hierarchical clustering tree was built by using dynamic hybrid cutting among which each branch represents genes with similar expression and each leaf was a single gene. Finally, eight modules were constructed (Fig. 1C). Among the eight modules, the turquoise module was highly correlated with CD8+ T cells (R2 = 0.74, P < .001) (Fig. 1D). Therefore, we identified the genes in the turquoise module as CD8+ T cell-correlated genes.

Figure 1.

Gene co-expression network of CRC. (A) The scale-free fit index of the soft threshold power (β) from 1 to 20. (B) The average connectivity of the soft threshold power (β) from 1 to 20. (C) Classify genes into different modules through hierarchical clustering. (D) Heatmap showing the correlation of module eigengenes with CD8+ T cell infiltration. CRC = colorectal cancer.

3.3. Functional analysis of DEGs correlated with CD8+ T cell in CRC

After analysis of the TCGA CRC cohort and GSE44076, we identified 1758 dysregulated genes in CRC tissues compared with normal tissues (Figs. 2A and B). Subsequently, we overlapped these DEGs with CD8+ T cell-correlated genes and further selected 834 CD8+ T cell-correlated DEGs differentially expressed genes correlated with CD8+ T cell (CD8TDEGs) in CRC (Fig. 2C). Consistently, GO analysis showed that most CD8TDEGs were associated with immune-related biological processes, such as T cell receptor signaling pathway, T cell activation, and immune response (Fig. 2D). In addition, KEGG pathway analysis showed that these CD8TDEGs mainly participated in immune-related signaling pathways, such as PI3K-AKT signaling pathway, Chemokine signaling pathway, and Cytokine-cytokine receptor interaction (Fig. 2E).

Figure 2.

Functional analysis of DEGs correlated with CD8+ T cell in CRC. (A) The differentially expressed genes in CRC tissues based on the TCGA database. (B) The differentially expressed genes in CRC tissues based on the GSE44076 dataset. (C) Overlapping differentially expressed genes and CD8+ T cell-related genes. (D) GO analysis of CD8TDEGs. (E) KEGG pathway analysis of CD8TDEGs. CD8TDEGs = differentially expressed genes correlated with CD8+ T cell, CRC = colorectal cancer, DEGs = differentially expressed genes.

3.4. Selecting prognostic DEGs correlated with CD8+ T cell in CRC

To select prognostic CD8TDEGs in CRC, we first carried out a univariate Cox regression analysis, and 32 prognostic CD8TDEGs were identified (Fig. 3A). Moreover, the Kaplan–Meier curve with a log-rank test was applied to further validate the prognostic values of CD8TDEGs. As shown in Figure 3B, CRC patients with low expression of GDF5, NAT2, and ACOT11 suffered from worse overall survival. On the contrary, CRC patients with low expression of CYRAB, CCBE1, RGS16, ITGBL1, FJX1, PCDH9, NOX4, TIMP1, SERPINH1, MEIS1, TNFAIP8L3, THBS2, CALB2, BGN, and BATF underwent favorable overall survival (Fig. 3C). Therefore, a total of 18 prognostic CD8TDEGs were identified in CRC.

Figure 3.

Selecting prognostic DEGs correlated with CD8+ T cell in CRC. (A) Univariate Cox regression analysis of CD8TDEGs in CRC. (B-C) The Kaplan–Meier curves of GDF5, NAT2, CYRAB, CCBE1, RGS16, ITGBL1, FJX1, PCDH9, NOX4, TIMP1, SERPINH1, MEIS1, TNFAIP8L3, THBS2, CALB2, BGN, and BATF in CRC. CD8TDEGs = differentially expressed genes correlated with CD8+ T cell, CRC = colorectal cancer.

3.5. Reclassifying CRC into 2 new clusters using the prognostic CD8TDEGs

Since above 18 CD8TDEGs were significantly associated with patients’ prognosis, we hypothesized that CRC patients may be reclassified into different subgroups with distinct overall survival. CRC samples were reclassified using the R package of “Consensus ClusterPlus”. The optimal clustering stability (k = 2–9) was determined by the proportion of ambiguous clustering measurements, and k = 2 was identified (Fig. 4A and B). Based on the unsupervised clustering, we eventually identified two distinct clusters (Fig. 4C). The Kaplan–Meier curve analysis showed that CRC patients in cluster 1 underwent worse overall survival than those in cluster 2 (HR = 2.57, P < .001) (Fig. 4D). Unexpected, the infiltrated levels of CD8+ T cells were significantly higher in CRC tissues in cluster 1 rather than in cluster 2 (Fig. 4E).

Figure 4.

Reclassify CRC according to the prognostic CD8TDEGs. (A-B) The optimal clustering stability and consensus index. (C) Consensus clustering analysis identified two subgroups. (D) Survival curves of cluster 1 and cluster 2 patients. (E) The infiltrated levels of CD8+ T cell in cluster 1 and cluster 2 tissues. CD8TDEGs = differentially expressed genes correlated with CD8+ T cell, CRC = colorectal cancer.

3.6. Constructing a prognostic model for CRC using the CD8TDEGs

Given that CRC patients can be reclassified into two subgroups with different prognoses based on CD8TDEGs, we subsequently constructed a prognostic model for CRC based on LASSO regression, and seven CD8TDEGs were selected (Fig. 5A and B). The risk score was calculated by the following formula: PCDH9 × (0.093) + GDF15 × (−0.23) + CALB2 × (0.0044) + ACOT11 × (−0.15) + FJX1 × (0.15) + TIMP1 × (0.12) + RGS16 × (0.033). According to the median value of risk scores, CRC patients were divided into high-risk and low-risk score groups. As shown in Figure 5C, CRC patients with higher risk scores trended to undergo worse overall survival. The time-dependent ROC analysis exhibited that the area under the curve was 0.70, 0.68, and 0.78 to predict patients’ prognosis at 1-year, 3-year, and 5-year, respectively (Fig. 5D). As expected, the risk scores in cluster 1 samples were significantly higher than that in cluster 2 (Fig. 5E). ROC analysis showed that the area under the curve was 0.87 to differentiate cluster 1 patients from cluster 2 patients (Fig. 5F). Moreover, we assessed the performance of our model in CRC patients with other clinicopathological characteristics by using multivariate Cox regression analysis and the results showed that the risk score and age were independent prognostic factors for CRC (Fig. 5G). Our results indicated that the prognostic model was robust in predicting the outcome of CRC patients.

Figure 5.

Constructing a prognostic model for CRC. (A) The coefficient profiles of the LASSO regression model. (B) Validating turning parameter selection of the LASSO regression model. (C) The Kaplan–Meier curve analysis of CRC patients in low- and high-risk score groups. (D) Time-dependent ROC analysis of the risk score to predict prognosis of CRC patients. (E) The risk scores in cluster 1 and cluster 2 tissues. (F) ROC analysis of the risk score to differentiate cluster 1 patients from cluster 2 patients. (G) The multivariate Cox regression analysis of risk score in CRC based on TCGA database. CRC = colorectal cancer, LASSO = least absolute shrinkage and selection operator.

3.7. Immune analysis of the prognostic model in CRC

We firstly compared the infiltered level of CD8+ T cell between the high-risk score group and low-risk score group. Intriguingly, CRC tissues in the high-risk score group have a higher infiltered level of CD8+ T cells (Fig. 6A). Besides, we found the expression of most CD8TDEGs in the prognostic model (PCDH9, CALB2, RGS16, TIMP1, and ACOT1) was either positively or negatively correlated with CD8+ T cell infiltration level (Fig. 6B–H). The GSEA analysis showed that the gene signatures of HALLMARK_ ALLOGRAFT REJECTION, HALLMARK_INFLAMMATORY_RESPONSE, HALLMARK_COMPLEMENT, HALLMARK_INTERFERON_GAMMA_RESPONSE, HALLMARK_INTERFERON_ALPHA_RESPONSE, HALLMARK_IL6_JAK_STAT3_SIGNALING, HALLMARK_TNFA_SIGNALING_VIA_NFKB, and HALLMARK_IL2_STAT5_SIGNALING were positively enriched in the low-risk score group (Fig. 6I). Moreover, we discovered that the immune checkpoint genes, including PD-1, PD-L1, and CTLA-4, were significantly higher in the high-risk score group compared to that in the low-risk score group (Fig. 6J–L). Therefore, we hypothesized that although more CD8+ T cells infiltrated CRC tissues in the high-risk score group, abundant immune checkpoint genes resulted in their exhaustion.

Figure 6.

Immune analysis of the prognostic model in CRC. (A) The infiltrated levels of CD8+ T cell in high- and low-risk score groups. (B–H) The correlations between infiltrated levels of CD8+ T cell and expression of PCDH9, CALB2, RGS16, TIMP1, ACOT1, GDF15, and FJX1 in CRC. (I) GSEA analysis of gene sets between low- and high-risk score groups. The expression levels of PD-1 (J), PD-L1 (K), and CTLA-4 (L) in high- and low-risk score groups. CRC = colorectal cancer.

3.8. scRNA-sequencing analysis of CD8TDEGs in prognostic model

To further explore the distribution of CD8TDEGs in CRC tissues, we conducted a scRNA-sequencing analysis of CRC tissues. As shown in Figure 7A, a total of nine clusters of cells were identified, including 5 CRC epithelial cell clusters, 1 B cell cluster, 1 T cell cluster, 1 fibroblast cluster, and 1 macrophage cluster. In addition, the expression values of CD8TDEGs in the prognostic model markedly varied among different cell clusters (Fig. 7B–H).

Figure 7.

scRNA-sequencing analysis of CD8TDEGs in the prognostic model. (A) The UMAP plot of identified cell clusters. (B-H) The expression values of PCDH9, CALB2, RGS16, TIMP1, ACOT1, GDF15, and FJX1 in various cell clusters. CD8TDEGs = differentially expressed genes correlated with CD8+ T cell, scRNA = single-cell RNA.

4. Discussion

CD8+ T cell is an important component of adaptive immunity.[ Emerging evidence suggests that enhancing the infiltration and response of CD8+ T cells can eliminate almost all cancer cells.[ In line with previous studies,[ our results exhibited that the infiltrating level of CD8+ T cells was significantly downregulated in CRC tissues was significantly related to the prognosis of patients.

Immunotherapy is expected to become a common treatment for solid cancer, alongside other treatments, such as chemoradiotherapy. Whereas CRC patients often occur tolerance to such therapy. For example, Xue et al discovered that intrinsic β-catenin signaling can suppress CD8+ T cell infiltration in CRC.[ In addition, cancer cells also result in T cell exhaustion by secreting immune escape factors such as PD‐L1 and IDO1 to inhibit the release of IFN‐γ by T cells.[ Consistently, our study revealed that due to more abundant immune checkpoint genes in tumor tissues, although part CRC patients own higher infiltrated levels of CD8+ T cells, their exhaustion leads to poor prognosis of patients. Recently, researchers discovered that the application of antibody‐based immunomodulation or vaccination can improve the priming ability of effector CD8+ T cells and establish a durable and efficient antitumor immunity.[ Monoclonal antibody‐based immunomodulation does not act on cancer cells directly but directs membrane ligands or receptors to regulate T cell responses.[ In this study, we reclassified CRC patients into two clusters with distinct immune-infiltrating cells and immunotherapy potential. Compared with those in the low-risk score, the CRC patients in the high-risk score group own higher CD8+ T cell-infiltrating levels and higher immune checkpoint gene expression (PD-1, PD-L1, and CTLA-4). Therefore, CRC patients with higher risk scores are more suitable to immune checkpoint blockade treatment.

We should take note that other immune cells in TME also regulate CD8+ T cell proliferation and infiltration. Tregs are a subset of the immunosuppressive CD4+ T cell family and can suppress the activity of CTLs by releasing TGF-β within the TME.[ Besides, to impair responses from CTLs, macrophages are often polarized to M2 to secrete STAT3 into the tumor microenvironment.[ Conversely, DCs are contributed to transfer stimulatory signals from CD4 + T cells to CD8+ T cells to trigger their tumor suppressor activity.[ Our scRNA-sequencing analysis exhibited that CRC cells contain 5 heterogenous subtypes and surrounded by various immune cells, such as macrophage. Therefore, we speculated that the effect of immunotherapy for CRC patients with high-risk scores may be improved by targeting different immune cells rather than one.

We must admit that there were several limitations in our study. First, we only analyzed the infiltrated levels and prognostic values of six types of immune cells. It will be better to enroll other immune cells, such as Tregs. Second, due to the limitation of cell number, scRNA sequencing of CRC tissues only screened nine cell clusters. Although the T cell cluster was identified, a detailed CD8+ T cell cluster was lacking. Third, the hypothesis that CRC patients with higher risk scores are more suitable for immune checkpoint blockade treatment should be further verified before its clinical application.

5. Conclusion

In summary, our study identified a novel prognostic model for CRC based on CD8+ T cell-related genes, according to which CRC patients can be classified into two clusters with distinct overall survival. In addition, our results revealed that immune checkpoint blockade treatments are more suitable for CRC patients with high-risk scores.

Acknowledgments

We acknowledge the TCGA and GEO databases for providing their platforms and contributors for uploading their meaningful datasets.

Author contributions

Conceptualization: Feng Pan.

Data curation: Zhenguo Pan.

Software: Qianjun Li, Yanling Feng.

Validation: Qianjun Li.

Visualization: Zhenguo Pan.

Writing – original draft: Feng Pan.

Writing – review & editing: Chengcheng Gao, Feng Pan.