Discovering biomarkers associated with infiltration of CD8+ T cells and tumor-associated fibrosis in colon adenocarcinoma using single-cell RNA sequencing and gene co-expression network

Abstract

Background: Colorectal adenocarcinoma (COAD) is a prevalent malignant tumor associated with a high mortality rate. Within the tumor microenvironment, CD8⁺ T cells play a pivotal role in the anti-tumor immune response within the human body. Fibrosis directly and indirectly affects the therapeutic response of tumor immunotherapy. However, the significance of regulatory genes associated with tumor-associated fibrosis and CD8⁺ T cell infiltration remains uncertain. Therefore, it is imperative to identify biomarkers with prognostic value and elucidate the precise role of CD8⁺ T cells and tumor-associated fibrosis.

Methods: We performed a single-cell transcriptome analysis of COAD samples from the GEO database. To evaluate immune infiltration in COAD samples, we utilized CIBERSORT and ESTIMATE. Furthermore, we analyzed the correlation between CD8⁺ T cells and immune infiltration. To analyze COAD expression's quantitative immune cell composition data, we conducted a Weighted Gene Correlation Network Analysis and utilized a deconvolution algorithm. The data for these analyses were obtained from the GEO database. We utilized univariate Cox regression and LASSO analysis to create a prognostic model. The predictive model was assessed through Kaplan-Meier analysis, and a survival prediction nomogram was created. Additionally, we analyzed the correlation between the prognostic model and chemotherapy drug sensitivity. To estimate the expression of hub genes, we employed immunohistochemistry, real-time PCR, and western blot techniques.

Results: Single-cell transcriptome analysis has indicated a higher prevalence of CD8⁺ T cells in COAD tumor samples. The connection between COAD and CD8⁺ T cells was further confirmed by WGCNA and deconvolution analysis using the GEO database. The Protein-Protein Interaction network analysis revealed three hub genes: LARS2, SEZ6L2, and SOX7. A predictive model was subsequently created using LASSO and univariate COX regression, which included these three genes. Two of these hub genes (LARS2 and SEZ6L2) were found to be upregulated in COAD cell lines and tissues, while SOX7 was observed to be downregulated. The prognostic model demonstrated a significant association with CD8⁺ T cells, suggesting that these genes could serve as potential biomarkers and targets for gene therapy in treating COAD.

Conclusion: This study has identified three key genes associated with CD8⁺ T cells and the prognosis of COAD, providing new prognostic biomarkers for diagnosing and treating COAD.

Keywords: CD8 + T cells; colon adenocarcinoma; fibrosis; prognostic biomarkers; single-cell RNA sequencing.

Figure 1

Resolution of the immune microenvironment in COAD and matched normal samples by scRNA-seq. (A) Cells were clustered into nine types via the tSNE dimensionality reduction algorithm; each color represented the annotated phenotype of each cluster. (B) Proportions of CD8⁺ T cells, CD4⁺ T cells, and Treg cells in each sample or tumor type. (C, D) Numbers and strength of interactions among cell types in COAD and healthy single cells.

Figure 2

Weighted gene co-expression network analysis (WGCNA) for hub module construction and validation. (A, B) PCA plots showing batch effects before and after correction using the SVA algorithm. (C) Soft-threshold power selection for scale-free network construction, depicting the scale-free topology model fit and mean connectivity. (D) Hierarchical clustering dendrogram of genes based on topological overlap, with distinct color-coded modules. (E) Module-trait relationships, with heatmaps demonstrating strong correlation between black/turquoise modules and CD8⁺ T cell infiltration (p < 0.05).

Figure 3

Enrichment analysis for black and turquoise modules. (A, B) GO and KEGG enrichment analysis for black module genes, revealing significant pathways in membrane trafficking, protein localization, and metabolism. (C, D) Enrichment analysis for turquoise module genes, highlighting pathways involved in ECM remodeling, vascular development, and growth factor signaling.

Figure 4

Identification of significant genes related to prognosis. (A, B) PPI networks of genes in the black and turquoise modules. (C) Random survival forest analysis using TCGA data, identifying genes with a relative importance > 0.4 as prognostic markers. (D) GeneMANIA interaction network illustrating functional relationships among the four hub genes. (E-H) Kaplan-Meier survival curves for LARS2, SEZ6L2, and SOX7, demonstrating significant associations between high/low expression levels and patient survival outcomes (log-rank test, p < 0.05).

Figure 5

Immune and clinical characteristics of these hub genes. (A) Expression profiles of LARS2, SEZ6L2, and SOX7 across different immune cell populations from TCGA data. (B-D) Immune infiltration analysis from the TIMER database, showing the correlation between hub gene expression levels and CD8⁺ T cell infiltration (p < 0.05). * represents a p-value < 0.05, ** represents a p-value < 0.01, and *** represents a p-value < 0.001. The smaller the p-value, the more statistically significant the result.

Figure 6

Gene Set Enrichment Analysis for hub genes. (A-C) GSEA enrichment plots for top three signaling pathways related to LARS2, SEZ6L2, and SOX7. (D-F) Bubble plots highlighting key genes involved in the enriched pathways, where larger circles represent higher enrichment scores.

Figure 7

Drug sensitivity data is based on the GDSC database. (A-C) Correlation of hub gene expression with drug response in the GDSC database. Results indicate that SEZ6L2 and SOX7 are associated with sensitivity to cisplatin, dasatinib, and gefitinib, while LARS2 expression correlates with response to cisplatin and gemcitabine.

Figure 8

Nomogram for predicting overall survival in COAD patients. (A) Nomogram-based prognostic scoring system, integrating clinical factors and hub gene expression for 3-year and 5-year survival prediction. (B) Calibration curves demonstrating the high predictive accuracy of the nomogram (AUC = 0.72), compared to TNM staging (AUC range: 0.62–0.68).

Figure 9

Experimental verification of signature gene expression in COAD. (A-C) RT-qPCR analysis of LARS2, SEZ6L2, and SOX7 mRNA expression in COAD and normal colonic epithelial cell lines. (D-F) Western blot analysis confirming higher LARS2 and SEZ6L2 protein levels and lower SOX7 protein expression in COAD cells compared to normal cells. (G) IHC analysis of LARS2, SEZ6L2, and SOX7 in COAD patient tissues, showing significant overexpression of LARS2 and SEZ6L2, and downregulation of SOX7 compared to matched normal tissues (p < 0.001). * represents a p-value < 0.05, ** represents a p-value < 0.01, and *** represents a p-value < 0.001. The smaller the p-value, the more statistically significant the result.