The role of cancer stem cells in the progression of colorectal cancer and the tumor microenvironment: new perspectives from single-cell RNA sequencing and spatial transcriptomics analyses

Introduction

Globally, colorectal cancer (CRC) is the third most commonly diagnosed cancer, representing 9.6% of all cancer cases, with approximately 1.93 million new cases in 2022. It ranks as the second leading cause of cancer-related deaths, accounting for 9.3% of all cancer fatalities (1). CRC is influenced by a multitude of factors, encompassing genetic predispositions, environmental and lifestyle elements—such as smoking, heavy alcohol consumption, and obesity—as well as the composition of the colonic microbiota (2). Over the past fifty years, studies on various types of cancer have shown that there are many heterogeneous cells in tumors, and certain types of cells play a very important role in the occurrence and development of tumors (3,4).

Cancer stem cells (CSCs) have long been considered a key component of tumor development (5). It is reported that CSCs are tumorigenic, can self-renew and differentiate, thereby generating tumor heterogeneity, and are resistant to conventional treatments, including chemotherapy and radiotherapy (6,7). The latest research shows that CSCs no longer represent just one cell type, but cells with stem cell-like characteristics. The CSC state is dynamic and can exist in different phenotypic states (8,9). In addition, a recent study has shown that CRC cells acquire stemness through embryonic reprogramming (10). Its plasticity increases the ability to invade and metastasize, accelerating the progression of cancer (11).

Over the past few decades, with the continuous advancement of sequencing technology, people have developed from conventional transcriptome sequencing to single-cell transcriptome sequencing, achieving high-throughput gene expression data of single cells in tissues (12,13). Although single-cell sequencing has solved many problems for us, single-cell RNA sequencing (scRNA-seq) alone cannot provide spatial information. In addition, due to the lack of multi-region sampling, the spatial heterogeneity within tumors at single-cell resolution is still poorly understood. Spatial transcriptomics (ST) is a new biotechnology that can visualize and quantitatively analyze the transcriptome at spatial resolution in tissue sections, making up for the lack of spatial information in scRNA-seq (14). However, the main limitation of ST is the lack of cellular resolution, whereas scRNA-seq can compensate for this shortcoming. However, systematic studies on the stemness of CRC utilizing single-cell transcriptomics and ST sequencing are still limited.

In this study, we performed scRNA-seq and ST analysis on primary CRC. We identified the most stem-like tumor epithelial cell population in CRC and analyzed its interactions with other cell types and its spatial distribution characteristics. We found that fibroblasts co-localize with and regulate the development of stem cells. By correlating CSCs’ plasticity states with spatial niches, we are seeking therapeutic targets within CSCs/fibroblasts-rich regions, and hope that the results of this study would provide guidance for finding new therapeutic targets for CRC.

Methods

Data source

The publicly available scRNA-seq data (GSE132465) and bulk RNA-seq data (GSE33113, GSE39582) can be acquired from the Gene Expression Omnibus database (GEO, https://www.ncbi.nlm.nih.gov/geo/). Transcriptomic data and clinical information of The Cancer Genome Atlas Colorectal Adenocarcinoma (TCGA-COAD) cohort were downloaded from the UCSC Xena data portal (https://xenabrowser.net). Spatial transcriptome data from patients were sourced from the 10x Genomics official website, with specific datasets available at the following links: (https://www.10xgenomics.com/datasets/human-colorectal-cancer-11-mm-capture-area-ffpe-2-standard) and (https://www.10xgenomics.com/datasets/human-colorectal-cancer-whole-transcriptome-analysis-1-standard-1-2-0).

Gene Ontology (GO) gene sets were retrieved from the MsigDB website (http://software.broadinstitute.org/gsea/msigdb).

Functional enrichment analysis

To investigate the differences in gene expression between the identified subtypes, we employed the DESeq2 R package for differentially expressed genes (DEGs) analysis, applying rigorous thresholds of an adjusted P value below 0.01 and an absolute log₂ fold change (FC) exceeding 2 (15). Upon identifying the DEGs, we proceeded with a comprehensive functional analysis, encompassing GO, Kyoto Encyclopedia of Genes and Genomes (KEGG), Gene Set Enrichment Analysis (GSEA). This analysis, facilitated by the ClusterProfiler R package (16). The analysis utilized GMT files sourced from reputable public databases such as MSigDB and consensus PathDB, which are based on pathway and gene set annotations.

Analysis of the CSC-associated prognostic genes

Initially, we refined the prognostic characteristics of the TCGA-COAD cohort by univariate Cox regression analysis, using CSC-related biomarkers and overall survival (OS) data. We selected 10 genes with differences in prognosis and analyzed the expression location and prognosis analysis of genes with hazard ratios (HRs) >1 to clarify the function of related genes.

TIDE analysis

In order to explore the sensitivity of clusters to immunotherapy, TIDE analysis was performed. TIDE score files were obtained from TIDE web resource (http://tide.dfci.harvard.edu/) (17).

Copy number variation (CNV) analysis via scRNA-seq

We applied the InferCNV R package (version 1.6.0, https://github.com/broadinstitute/inferCNV) to estimate CNVs from scRNA-seq data of CRC cells. This method sorts genes based on chromosomal position and applies a moving average to the relative expression, enabling the inference of CNVs.

Drug sensitivity prediction

To assess patients’ predisposition towards molecular agents, we employed the Genomics of Cancer Drug Sensitivity (GDSC) database for comprehensive evaluation (18). Leveraging the R package pRRophetic, we estimated the half-maximal inhibitory concentration (IC₅₀) and conducted rigorous cross-validation to validate the accuracy of our predictions.

Robust cell type decomposition (RCTD) analysis

We used RCTD to map cell types identified from reference scRNA-seq datasets onto ST data (19). We utilized the Seurat function FindAllMarkers to isolate marker genes for each cell type, focusing only on those with positive FCs post-log₂ transformation. Once the marker genes were identified, we standardized the processing of both the reference datasets and the Visium ST data to full duplex mode, adhering closely to the established RCTD analytical protocol.

ST data analyses

Utilizing the Seurat package (version 4.4.0), we delved into the analysis of spatial RNA-seq data. Commencing with normalization, we applied the sctransform function to standardize gene expression matrices for each sample. Next, we performed principal component analysis (PCA) and uniform manifold approximation and projection (UMAP) for dimensionality reduction (20). For clustering, we implemented the Louvain algorithm, optimizing with parameters dims =20 and res =0.2. To identify DEGs across clusters, we rigorously applied the Wilcoxon test, selecting genes with an adjusted P value <0.05 and a logFC >0.25 as markers. To assess the proportions of distinct cell types within clusters, we conducted a Fisher’s exact test, using all genes as a backdrop, to statistically evaluate the overlap between ST marker genes and those identified through scRNA-seq. This approach enabled us to compute odds ratios and associated P values, thereby uncovering patterns of cell type enrichment.

Multimodal intersection analysis (MIA)

We combined publicly accessible scRNA-seq data with ST data using the previously reported MIA method (21). The hypergeometric cumulative distribution method was used to determine the significance of overlap between cell type marker genes and ST genes, and the P value was calculated using all genes as the background.

Spatial map of cell dependencies

We employed the MISTy implementation in mistyR (v1.8.1) to assess the importance of individual cell type abundances. This assessment aimed to explain the abundance patterns of other major cell types. We formulated a multi-view model to estimate cell type positions across all slides using three different spatial views: (I) an intrinsic view that analyzes the relationship between deconvolution estimates within a point; (II) a juxtaposition view that sums deconvolution estimates of immediate neighbors (with distance thresholds up to 2); and (III) a peer view that weights deconvolution estimates of more distant neighbors of each cell type (within an effective radius of 5 points).

Pseudotime analyses

Pseudotime analyses, trajectory construction, and calculation of pseudotime-dependent gene expression dynamics along these trajectories, were executed utilizing Monocle 2 (version 2.14.0) (22). The Seurat-object comprising the integrated CRC dataset was transformed into a format compatible with Monocle 2, specifically a CellDataSet. The analysis adhered to the standard procedural guidelines.

Cellchat analysis

We utilized CellChat (version 1.4.0) to deduce intercellular interactions among various cell types within CRC microenvironment (23). Leveraging the database integrated within CellChat, which encompasses categories such as “Secreted Signaling” and “Cell-to-Cell Contact”, we conducted comparative analyses between CRC patient samples and control samples in accordance with the guidelines provided in the CellChat tutorial (https://htmlpreview.github.io/?https://github.com/sqjin/CellChat/blob/master/tutorial/Comparison_analysis_of_multiple_datasets.html). This approach facilitated the identification of dysregulated signaling pathways, characterized by an upregulation or downregulation of specific ligand-receptor interactions between the disease and control states.

Statistical analyses

This study used R software version 4.3.0 for statistical analysis. When dealing with quantitative data, we used the Student’s t-test to evaluate the statistical significance of variables that were normally distributed; for variables that were not normally distributed, the Wilcoxon rank sum test was used for analysis. When it comes to comparisons between two or more groups, we selected the non-parametric Kruskal-Wallis test and the parametric one-way analysis of variance (ANOVA) as statistical methods based on the distribution of the data. For the analysis of contingency table data, we selected the chi-squared test and Fisher’s exact test as appropriate statistical tools based on the specific grouping conditions and sample size. Furthermore, to gain insight into the relationship between metabolic transcriptomic patterns and patient prognosis, we applied Kaplan-Meier (K-M) survival analysis and Cox proportional hazards models with the help of the “Survminer” package in R (version 0.4.6). All statistical comparisons were performed as two-sided tests, and the significance level (alpha) was set at 0.05.

Ethical statement

This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Results

Mapping cellular heterogeneity in CRC with scRNA-seq analysis

By leveraging the Seurat algorithm to integrate the GSE132465 datasets, we identified a comprehensive set of 23,970 genes across 57,068 cells. These datasets encompass a total of 33 patient samples, comprising 23 primary CRC tissues and 10 normal colorectal tissues. After meticulously integrating diverse samples to mitigate batch effects, we proceeded with standardization and centralization to prepare the data for further analysis. To reduce dimensionality, we applied PCA, which allowed us to capture the major sources of variability in our dataset. Subsequently, we employed the UMAP method for clustering and visualization. This approach effectively segmented the 57,068 cells into 16 distinct clusters (Figure 1A). Additionally, we visualized the cellular distribution patterns between tumor and normal tissues (Figure 1B), revealing a predominant increase in cell abundance within tumor tissues compared to the relatively lower cell density observed in normal tissues. Employing specific highly expressed genes for cell annotation, we delineated nine predominant cell types: T cells, B cells, plasma cells, myeloid cells, epithelial cells, mast cells, glial cells, endothelial cells, and fibroblasts (Figure 1C-1E). The results indicated a higher proportion of epithelial and myeloid cells in tumor tissues, contrasted with a relatively diminished presence of fibroblasts, B cells, and plasma cells (Figure 1F). The expression levels of the top 5 DEGs in 9 major cell types showed specificity, mainly concentrated in the corresponding cell types with high expression, while in other cell types with low expression (Figure 1G).

Figure 1 Mapping cellular heterogeneity in CRC with scRNA-seq analysis. (A) UMAP algorithm was applied to the top 20 principal components for dimensionality reduction, and 16 cell clusters were successfully classified. (B) UMAP plots of 57,068 cells from normal and tumor tissues. (C) Bubble plots of the marker genes expressed in each cell type. Size of each bubble represents the percent of cells expressing marker and color represents the level of expression. (D) UMAP of 57,068 cells colored by major cell types. (E) Feature plots for the canonical marker genes of epithelial cell (EPCAM), endothelial cell (PECAM1), fibroblast (COL3A1), myeloid cell (CD68), T cells (CD3D), B cells (CD79A), mast cell (TPSB2), plasma B cell (TNFRSF17), glial cell (PLP1). (F) Bar plot visually represents a comparison of the proportions of the primary cell types present in both normal and tumor tissues. (G) Heatmap of the expression levels of the top 5 DEGs among 9 major cell types in CRC. CRC, colorectal cancer; DEGs, differentially expressed genes; scRNA-seq, single-cell RNA sequencing; UMAP, uniform manifold approximation and projection.

Identification of malignant and non-malignant epithelial cells in CRC

As depicted in Figure 1G, there is a significant difference in the expression patterns of epithelial cells between normal and tumor tissues. Consequently, we extracted all epithelial cells for re-clustering analysis (Figure 2A). Clusters 7, 13, and 14 predominantly consist of epithelial cells from normal tissues (Figure 2B,2C). To explore the differences between epithelial cells in tumor and normal tissues, we conducted differential expression analysis on epithelial cells from tumor and normal samples to identify DEGs associated with tumor samples. We performed enrichment analysis on the highly expressed DEGs. KEGG enrichment analysis revealed that the DEGs are primarily involved in the ribosome and IL-17 signaling pathways (Figure 2D). GO enrichment results indicated significant enrichment in processes such as cytoplasmic translation, ribonucleoprotein complex biogenesis, focal adhesion, cytoplasmic ribosome, cell-substrate junction, cadherin binding, structural constituent of ribosome, and unfolded protein binding (Figure 2E). To further determine whether the epithelial cells within the tumor tissue are malignant, we extracted and re-clustered the epithelial cells for re-annotation. Post-clustering, the epithelial cells from the tumor samples were categorized into 14 distinct clusters (Figure 2F). We employed InferCNV to perform CNV scoring, which served to gauge the degree of malignancy within the cell populations (Figure 2G). Subsequently, the CNV scores were subjected to K-means machine learning clustering, revealing that a cluster count of 9 yielded optimal results (Figure 2H). Notably, Cluster 9, with the lowest score (Figure 2I), was preliminarily identified as non-malignant cells (Figure 2J). We proceeded to analyze the metabolic profiles between non-malignant and malignant cells (Figure 2K, Figure S1). In malignant cells, the functionality of numerous metabolic pathways was typically diminished, encompassing amino acid metabolism, carbohydrate metabolism, and certain aspects of fatty acid metabolism. In contrast, purine metabolism was found to be relatively enhanced.

Figure 2 Profiling cellular complexity in CRC: an scRNA-seq analysis of epithelial cells and malignant identification. (A,B) UMAP of epithelial cells, color-coded by clusters (A) and tissue group (B). (C) The distribution of cellular proportions within each cluster across both normal and tumor tissues. (D,E) Enrichment analysis of upregulated DEGs in epithelial cells between tumor and normal tissues using KEGG pathways (D) and GO terms (E). (F) UMAP of epithelial cells in tumor tissue. (G) Large-scale CNVs of epithelial cells from tumor tissue were displayed by hierarchical heatmaps to identify malignant cells. (H) Elbow method was used to determine the optimal clustering number. (I) Violin plot visualizing the CNV scores of K-means clusters. (J) UMAP plot of epithelial cells, color-coded by malignant condition. (K) Analysis of metabolism between malignant and non-malignant epithelial cells. ns, not significant; *, P<0.05; **, P<0.01; ***, P<0.001. BP, biological process; CC, cellular component; CNV, copy number variation; CRC, colorectal cancer; DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; scRNA-seq, single-cell RNA sequencing; UMAP, uniform manifold approximation and projection.

Quantitative analysis of malignant epithelial cell stemness

The malignant epithelial cells were clustered and annotated (Figure 3A, Figure S2). Employing CytoTRACE, we evaluated the stemness and differentiation potential across various cell clusters, identifying clusters 0, 1, 2, 3, and 4 as those with notably higher stemness, with cluster 4 standing out with the highest stemness score (Figure 3B,3C). Malignant epithelial cells were then stratified into high and low stemness groups based on their respective scores. Using the Single Sample Gene Set Enrichment Analysis (ssGSEA) algorithm, we assessed these cells in relation to key pathways. Notably, malignant epithelial cells with elevated stemness scores exhibited significant activation of pathways like AKT, EMT, MAPK, NF-κB, and TGF-β. Interestingly, the high-stemness group showed reduced immune checkpoint scores, suggesting a potential vulnerability to immune evasion, which could be instrumental in disease progression (Figure 3D). Subsequently, we performed a correlation analysis between the stemness scores and the scores of the pathways. This analysis revealed that the stemness score had the strongest correlation with the oxidative stress score among all the pathway scores analyzed (Figure 3E). To further explore the differences between epithelial cells in the high and low CytoTRACE value groups, differential analysis was conducted, followed by KEGG enrichment analysis on the DEGs that were upregulated in the high CytoTRACE group. The results indicated that the TNF, ErbB, and NF-κB signaling pathways were significantly activated (Figure 3F).

Figure 3 Quantitative analysis of malignant epithelial cell stemness. (A) UMAP of the malignant epithelial cells, color-coded by clusters. (B,C) The tSNE plot illustrates the stemness states of malignant epithelial cells as determined by CytoTRACE analysis. (D) Violin plot visualizing the ssGSEA scores of tumor epithelial cells between high- and low-CytoTRACE values groups. (E) The correlation analysis between the ssGSEA scores and CytoTRACE values. (F) KEGG pathway enrichment analysis of upregulated DEGs in high-CytoTRACE values tumor epithelial cells. DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; ssGSEA, Single Sample Gene Set Enrichment Analysis; tSNE, t-distributed Stochastic Neighbor Embedding; UMAP, uniform manifold approximation and projection.

Deciphering the role and regulatory mechanisms of CRC subtypes in tumor progression

Elevated stemness in tumor cells is linked to unfavorable prognoses, indicating that an enriched presence of these cells could be a determinant of treatment resistance and diminished survival rates (24). Therefore, we further analyzed cluster 4 with the highest stemness and named it colorectal cancer stem cell subtype (Co-CSS). First, we enriched and analyzed the marker genes of Co-CSS, revealing enrichment in DNA-binding transcription factor binding, RNA polymerase II transcription regulator complex, granulocyte migration and significant activation of the TNF, IL-17, and TGF-β pathway (Figure 4A,4B). Next, we identified 33 stemness-related genes by intersecting the 45 marker genes of Co-CSS with the 6,995 genes of TCGA-COAD (Figure 4C). Subsequently, we conducted univariate Cox regression analysis on these 33 stemness-related genes, identifying 10 genes with varying prognoses. TRIP6 is a prognostic risk factor for CRC (Figure 4D). TRIP6 is highly expressed in tumors, and patients with high TRIP6 expression have a poor prognosis (Figure 4E,4F). At the scRNA-seq level, TRIP6 was found to be highly expressed in tumor epithelial cells exhibiting high stemness and CytoTRACE scores (Figure 4G,4H). In addition, TRIP6 has been shown to promote the metastasis and tumor progression of CRC (25). In summary, our findings underscore the significant role of Co-CSS in tumor progression. To investigate the transcriptional regulators that exert substantial control over Co-CSS, we conducted SCENIC analysis, pinpointing that MYB, NUAK2, and RUNX3 have the most pronounced transcriptional activity within Co-CSS (Figure 4I,4J). Inhibiting these transcription factors, which are highly active in Co-CSS, could be a strategic approach to impede the malignant progression of these cells, presenting a potential therapeutic target for intervention.

Figure 4 Deciphering the role and regulatory mechanisms of CRC stem subtypes in tumor progression. (A) Enrichment analysis of GO terms for the marker genes of Co-CSS. (B) KEGG pathway enrichment analysis of the marker genes of Co-CSS. (C) Intersection of the gene sets obtained from the differential genes of TCGA-COAD and marker genes of Co-CSS, resulting in a final set of 33 genes. (D) Univariate Cox regression analysis of OS. (E) Comparative analysis of TRIP6 expression between paired tumor and normal tissues from the TCGA-COAD. (F) Kaplan-Meier curve result. (G-H) Dot plots of TRIP6 expressed in the malignant epithelial cells. Dot color reflects expression level and dot size represents the percent of cells expressing marker genes in different cell clusters (G) and CytoTRACE groups (H). (I) Violin plot visualizing the top 5 TF of regulon specificity score in Co-CSS. (J) UMAP plot illustrates the transcription factors of MYB, NUAK2, and RUNX3 in single-cell profiles of Co-CSS. ***, P<0.001. AUC, area under the curve; BP, biological process; CC, cellular component; Co-CSS, colorectal cancer stem subtype; COAD, colon adenocarcinoma; CRC, colorectal cancer; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; OS, overall survival; TCGA, The Cancer Genome Atlas; TF, transcription factor; UMAP, uniform manifold approximation and projection.

Interplays between Co-CSS and tumor microenvironment (TME) in CRC

To identify potential interactions between Co-CSS and other cells in the TME, we engaged CellChat analysis, uncovering a rich tapestry of intercellular signaling links. Co-CSS showed strong interactions with endothelial cells and fibroblasts, with a significant number of interactions, particularly with fibroblasts, where the number and intensity of interactions were high (Figure 5A,5B). Delving deeper, we scrutinized the specific ligand-receptor dynamics at play between Co-CSS and these interacting cells. Notably, we observed significant upregulation in the ligand-receptor pairs associated with endothelial cells (CD99-CD99, COL4A1-SDC4, COL4A2-SDC4) and fibroblasts (COL1A1-SDC4, COL1A2-SDC4, COL6A1-SDC4, COL6A2-SDC4), underscoring the pivotal role these pathways serve in the modulation of Co-CSS by their microenvironmental neighbors (Figure 5C). Expanding our inquiry, we investigated the broader signaling pathways at work and discovered that fibroblasts exert regulatory control over Co-CSS via the WNT and THBS1 pathways (Figure 5D, Figure S3A-S3C). Intriguingly, the ligand WNT2 was found in the exclusive purview of fibroblasts, while its receptors, FZD5, FZD7, and LRP5, were distinctly present within Co-CSS (Figure 5E,5F). The IGF pathway emerged as a potential conduit through which endothelial cells might influence Co-CSS, with the IGF2 ligand expressed in endothelial cells and its receptors ITGA6 and ITGB4 identified in Co-CSS (Figure 5G-5I). Furthermore, myeloid cells were uncovered as regulators of Co-CSS through the TNF pathway, and mast cells appeared poised to engage Co-CSS via the LIGHT pathway (Figure S3D-S3I). These pathways, each with established ties to tumorigenesis and progression, highlight the WNT signaling pathway as a critical catalyst for the proliferation of tumor stem cells and the preservation of their stemness (26). In aggregate, these findings suggest that diverse cellular actors within the TME wield regulatory influence over the genesis and evolution of Co-CSS through a complex network of signaling pathways.

Figure 5 Cross-cellular regulatory networks within the tumor microenvironment dictate the role of Co-CSS in CRC progression. (A,B) Heatmaps depict the differential number (A) and strength (B) intercellular interactions between lymph nodes and tumors, accompanied by an overview of cell-cell interactions. The direction of the interactions is signified by the color of the arrows and edges, with the size of the circles representing the proportion of cells within each cell group. The thickness of the edges corresponds to the frequency and intensity of interactions between the cell populations. (C) Bubble plots of the significant differentially expressed ligand-receptor pairs. Dot color reflects communication probabilities, and dot size represents computed P values. Each cell group is shown in different colors. Empty space means the communication probability is zero. P values are computed from a two-sided permutation test. (D) Contribution heatmap of cell communication in WNT signaling pathways. (E) Dot plots of the significant differentially expressed ligand-receptor pairs in WNT signaling pathways. (F) Violin plot visualizing ligand-receptor pairs WNT signaling pathways in all cell clusters. (G) Contribution heatmap of cell communication in IGF signaling pathways. (H) Dot plots of the significant differentially expressed ligand-receptor pairs in IGF signaling pathways. (I) Violin plot visualizing ligand-receptor pairs of IGF signaling pathways in all cell clusters. Co-CSS, colorectal cancer stem cell subtype; CRC, colorectal cancer.

Spatial heterogeneity of Co-CSS unveiled by ST mapping

To explore the spatial distribution characteristics of CRC stem subtypes, we downloaded ST data from a patient’s sample obtained through the 10x Genomics official database. Slide 1 contained 9,080 spots with an average count of 18,085 transcripts per spot (nCounts). Following stringent quality control measures, genes expressed in fewer than 10 spots were filtered out. This process resulted in a refined dataset with 18,043 nCounts for slide 1 and 16,855 for slide 2. The hematoxylin and eosin (H&E) stained images of the tissue slides were meticulously assessed. Utilizing dimensionality reduction and clustering techniques on the spots of each ST array, we discerned 10 predominant cell clusters within slide 1, distinguished by their unique gene expression signatures (Figure 6A). To delineate the cellular composition across various regions within the tumor, we employed the MIA algorithm for cell annotation. We identified that clusters 1, 7, 9, and 0 predominantly contained epithelial cells, while cluster 6 was characterized by myeloid cells. Additionally, clusters 2, 5, and 1 were found to harbor endothelial cells, fibroblasts, B cells, T cells, and other cell types, primarily distributed in the interstitial areas (Figure 6B,6C). The primary cellular constituents of clusters 3, 4, and 8 remained ambiguous, prompting further investigation with cell-specific markers. Notably, epithelial cell markers were observed to be highly expressed in clusters 0, 3, 4, 7, 8, and 9, and macrophage markers were notably elevated in cluster 6. Consistent patterns were also observed for other cell types (Figure S4). Subsequently, we applied the ssGSEA algorithm to score the top 50 marker genes of the Co-CSS subtypes within individual cells of the ST data from slide 1. This analysis revealed that the area of cluster 9 exhibited the highest score, suggesting a rich Co-CSS presence (Figure 6D,6E). To bolster the validity of our findings, we re-engaged the MIA algorithm to integrate single-cell stem cells and other cell subtypes with the spatial transcriptome data. Interestingly, this approach predominantly identified stem cells within the area of cluster 9 (Figure 6D). Finally, to provide auxiliary validation, Co-CSS marker genes ASCL2 and PTPRO were assessed, corroborating that Co-CSS-associated genes were highly expressed in cluster 9 (Figure 6F). Collectively, our results indicate a significant concentration of Co-CSS in cluster 9 (Figure 6G).

Figure 6 Spatial heterogeneity of Co-CSS unveiled by ST mapping. (A) Hematoxylin and eosin staining of tissue sections and UMAP plot showing the spot clusters identified using unbiased clustering of the ST spots in slide 1. (B) The MIA map of all scRNA-seq identified cell types and ST-defined regions in slide 1. (C) Expression of the epithelial markers in each cluster identified in slide 1. (D) Violin plots and dot plots of CRC scores for each cluster. (E) The MIA map of scRNA-seq identified cell types (Co-CSS and other cell types) and ST-defined regions in slide 1. (F) Violin plot visualizing the marker genes of Co-CSS. (G) Spatial feature plots of Co-CSS. The spatial transcriptomics data were obtained from 10x Genomics public repository (https://www.10xgenomics.com/datasets/human-colorectal-cancer-11-mm-capture-area-ffpe-2-standard and https://www.10xgenomics.com/datasets/human-colorectal-cancer-whole-transcriptome-analysis-1-standard-1-2-0) and processed following their standard FFPE protocol (CG000518, CG000520, CG000495). Scale bars are not shown due to a lack of scale information in the raw dataset. CRC, colorectal cancer; Co-CSS, colorectal cancer stem cell subtype; CSS, cancer stem-like signature cells; MIA, multimodal intersection analysis; scRNA-seq, single-cell RNA sequencing; ST, spatial transcriptomics; UMAP, uniform manifold approximation and projection.

Spatial co-localization of fibroblasts and their regulatory impact on Co-CSS

Based on the expression patterns of malignant epithelial cells and Co-CSS, we observed that cluster 9, where Co-CSS reside, was surrounded by cluster 2 (Figure 7A). Consequently, we hypothesize that Co-CSS are spatially co-localized with these diverse cell types. To elucidate the cellular distribution within cluster 9, we harnessed the RCTD method for the classification of spatially resolved transcriptomic profiles obtained through scRNA-seq. This sophisticated analytical technique enabled us to dissect the heterogeneous cell populations present at each spot, culminating in the identification of five unique cell types (Figure 7B). Moreover, a cell proportion graph revealed that cluster 2 has a relatively high proportion of fibroblasts, endothelial cells, myeloid cells, and other cell types, while cluster 9 contains Co-CSS, fibroblasts, myeloid cells, and endothelial cells, with fibroblasts being the most abundant in cluster 2 (Figure 7C).

Figure 7 Spatial co-localization of fibroblasts and their regulatory impact on Co-CSS. (A) Hematoxylin and eosin staining of tissue sections and UMAP plot showing the cluster 2 and 9 in slide 1. (B) Cell types of each spatial spot revealed by RCTD. (C) The percentage of all cell types in each cluster. (D) Mean importance of the abundance of major cell types within a spot. (E) Mean importance of the abundance of major cell types within a two-spot radius. (F) Mean importance of the abundance of major cell types within a 5-spot radius. (G) KEGG enrichment analysis of cluster 9. (H) Correlation between CSS score and WNT signaling pathway score. (I) Spatial feature plots display the ssGSEA scores for CSS and WNT signaling pathway in slide and corresponding line graphs show the ssGSEA scores along the trajectory, providing a detailed view of CSS and WNT signaling pathway activity changes. The spatial transcriptomics data were obtained from 10x Genomics public repository (https://www.10xgenomics.com/datasets/human-colorectal-cancer-11-mm-capture-area-ffpe-2-standard and https://www.10xgenomics.com/datasets/human-colorectal-cancer-whole-transcriptome-analysis-1-standard-1-2-0) and processed following their standard FFPE protocol (CG000518, CG000520, CG000495). Scale bars are not shown due to lack of scale information in the raw dataset. Co-CSS, colorectal cancer stem cell subtype; CSS, cancer stem-like signature cells; KEGG, Kyoto Encyclopedia of Genes and Genomes; RCTD, Robust Cell Type Decomposition; ssGSEA, Single Sample Gene Set Enrichment Analysis; UMAP, uniform manifold approximation and projection.

Further, to more accurately determine the co-localization relationship between Co-CSS and other cells, we utilized the MISTy algorithm to assess three different neighborhood sizes: (I) the number of cells within a single spot; (II) the number of cells within a two-spot radius; and (III) the number of cells within a five-spot radius. This multi-tiered analysis not only underscored the importance of cell type distribution and their co-localization (Figure 7D) but also provided insights into the microenvironment’s influence on Co-CSS at different scales—from the local context of a two-spot radius (Figure 7E) to the wider interactions within a five-spot radius (Figure 7F). We observed that fibroblasts can best predict the abundance of Co-CSS in spots and local areas, indicating that the distribution of fibroblasts and Co-CSS is spatially co-localized and likely engaged in intercellular communication. To substantiate the reliability of our findings, we validated them in slide 2. The MISTy algorithm (Figure S5A-S5G) evaluated the distribution of Co-CSS in cluster 8 of slide 2, also identifying the co-localization of fibroblasts and Co-CSS in specific spots. Interestingly, we enriched marker genes in cluster 9 and found that the WNT signaling pathway is significantly activated (Figure 7G). Moreover, in the Co-CSS, we found that the CSS score was positively correlated with the WNT signaling pathway score (Figure 7H), and it also showed a positive correlation trend in space. We used SPATA spatial trajectory analysis to define an arrow in a direction from the low CSS score region to the high score region, and found that the WNT signaling pathway score table showed a trend similar trend to the CSS score (Figure 7I). The same result was also presented in slide 2 (Figure S5H,S5I). To confirm the reliability of our findings again, we repeated the test in slide 3 and got the same results as slide 1 and 2 (Figure S6). In summary, given that the single-cell results suggested that fibroblasts may regulate Co-CSS through the WNT pathway, and the present study found a spatial correlation between fibroblasts and Co-CSS, and a spatial co-localization between WNT signaling pathway and Co-CSS. At the spatial transcriptome level, it was also speculated that Co-CSS might be regulated by WNT pathway in fibroblasts.

Clinical potential of NPDC1 and NSMF in treatment of CRC

The cluster 9 in slide 1 was enriched with Co-CSS and fibroblasts that facilitate the progression of Co-CSS. Consequently, we aimed to identify clinically relevant therapeutic targets from the marker genes in cluster 9, enabling more precise treatment strategies. Initially, we intersected the top 100 highly expressed genes in cluster 9 with those highly expressed in tumor tissues from TCGA data, yielding 57 overlapping genes (Figure 8A). Univariate Cox analysis was employed to sift through these genes for those with prognostic significance. We identified NPDC1 and NSMF as prognostic risk factors for patients (Figure 8B), and K-M survival analysis revealed that patients with high expression of either NPDC1 or NSMF had a poorer prognosis compared to those with low expression (Figure 8C). Additionally, in slide 1, both NPDC1 and NSMF were highly expressed in the Co-CSS-enriched regions (Figure 8D). In single-cell data, NPDC1 and NSMF were found to be highly expressed in the epithelial cells of tumor tissues (Figure 8E), a trend consistent with TCGA-COAD data (Figure 8F). Notably, NPDC1 showed particularly robust expression in tumor cells with high stemness or high CNV characteristics (Figure 8G). Pathway expression scoring indicated that high expression of NPDC1 is also associated with high scores in the MAPK, TGF, oxidative stress, and EMT pathways (Figure 8H). Pseudotime series analysis of tumor cells demonstrated an initial increase followed by a decrease in NPDC1 expression during tumor cell development (Figure 8I), suggesting that NPDC1 may play a role in driving tumor cell differentiation. In summary, our findings indicate that in CRC, tumor cells with high NPDC1 expression are more malignant and may be associated with tumor cell metastasis.

Figure 8 Targeting cancer stemness: the clinical potential of NPDC1 and NSMF in treatment of CRC. (A) A Venn diagram illustrates the intersection of gene sets derived from the DEGs of TCGA-COAD and the marker genes of cluster 9, culminating in a final set of 57 genes. (B) Univariate Cox regression analysis based on OS in TCGA-COAD cohort. (C) K-M curve of NPDC1 (left) and NSMF (right). (D) Bubble plots of NPDC1 and NSMF expressed in the spot cluster. (E) Bubble plots of NPDC1 and NSMF expressed in the malignant epithelial cells. Dot color reflects expression level and dot size represents the percent of cells expressing marker genes in normal and tumor tissues. (F) Expression of NPDC1 and NSMF between normal and tumor tissues. (G,H) Dot plots of NPDC1 and NSMF expressed in the malignant epithelial cells. Dot color reflects expression level and dot size represents the percent of cells expressing marker genes in different groups. (I) Potential trajectory of all malignant epithelial cells identified two distinct cell fates colored by state. The arrow shows the potential evolutionary direction in the trajectory. Dot plots of dynamic expression of NPDC1. States 1-7 represent different developmental stages. **, P<0.01; ***, P<0.001. CI, confidence interval; COAD, colon adenocarcinoma; CRC, colorectal cancer; DEGs, differentially expressed genes; HR, hazard ratio; OS, overall survival; TCGA, The Cancer Genome Atlas.

The classifier based on NPDC1 and NSMF predict the prognosis of CRC patients

To delineate the roles of NPDC1 and NSMF in CRC, we stratified CRC patients based on the expression of these genes, seeking to uncover potential differences among the identified clusters. Through consensus clustering, we discerned 2 distinct clusters within the TCGA-COAD, each characterized by a unique gene expression signature and labeled as A and B (Figure 9A). Cluster A represented patients with high expressions of NPDC1 and NSMF, while cluster B represented patients with low expressions of both NPDC1 and NSMF (Figure 9B). Survival analysis showed that cluster B had a great prognosis compared with cluster A (Figure 9C). Enrichment analysis found that in cluster A, glycosphingolipid biosynthesis, P53 pathway, IL-17 pathway and purine metabolism were significantly activated (Figure 9D). To assess the potential responsiveness of clusters A and B to immunotherapy, a Tumor Immune Dysfunction and Exclusion (TIDE) analysis was conducted. Figure 9E illustrates that there is no significant disparity in the TIDE scores between clusters A and B. However, cluster A exhibits a higher dysfunction score, suggesting an overabundance of immunosuppressive elements such as PD-L1 and CTLA-4 within the TME of patients. This overexpression may facilitate the tumor cells’ evasion of immune detection. Furthermore, cluster A has a lower microsatellite instability (MSI) score, which may foreshadow a suboptimal response to immunotherapy. Through immune checkpoint analysis, it was found that LAG3, CD276, TNFSF9, TNFRSF14, LGALS9, TNFRSF4, TNFRSF18, and TNFRSF25 were highly expressed in cluster A, indicating possible effective immunotherapy targets (Figure 9F). Finally, the oncoPredict package was used for drug sensitivity analysis. First, we found that cluster A showed higher sensitivity to the commonly used chemotherapy drug fluorouracil, while cluster B showed higher sensitivity to oxaliplatin and regorafenib (Figure 9G,9H). In clinical practice, patients who are resistant to fluorouracil alone can be tested for NPDC1 and NSMF. In addition, we screened drugs with different sensitivities in clusters A and B and drugs that are very sensitive in both groups for clinical reference.

Figure 9 Exploring the clinical significance of NPDC1 and NSMF isoforms in CRC. (A) Consensus clustering based on NPDC1 and NSMF. (B) Box plots show the expression of several immune checkpoints between A and B groups (Wilcoxon rank sum test) (C) Survival analysis of the patients grouped by cluster A and B in TCGA-COAD. (D) KEGG enrichment analysis of upregulated DEGs in cluster A. (E) TIDE algorithm analysis for cluster A and B. (F) Immune checkpoint genes expression. (G,H) Estimated IC₅₀ values of the therapy drugs, chemotherapy drugs and targeted molecular drugs were compared between 2 clusters. ns, not significant; *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. Cluster A represented patients with high expressions of NPDC1 and NSMF, while cluster B represented patients with low expressions of both NPDC1 and NSMF. COAD, colon adenocarcinoma; CRC, colorectal cancer; DEGs, differentially expressed genes; IC₅₀, half-maximal inhibitory concentration; KEGG, Kyoto Encyclopedia of Genes and Genomes; MSI, microsatellite instability; TCGA, The Cancer Genome Atlas; TIDE, Tumor Immune Dysfunction and Exclusion.

Discussion

Similar to other malignancies, CRC is a heterogeneous disease containing cell subtypes with different phenotypes and functions (26,27). CSCs refer to the stem-cell-like phenotype of cancer cells, which have the ability to self-renew, differentiate and initiate tumors, and have been widely recognized as an important player in tumor occurrence and development (24,28,29). An increasing number of studies have shown that cancer stemness increases with tumor progression. Such self-renewing tumor cells drive tumor occurrence, metastasis and drug resistance, and these phenotypes lead to poor clinical outcomes (30-32). With further exploration, we found that the stem-cell-like phenotype of cancer cells is a plastic phenotype, and the degree of cancer stemness can be regulated by many factors, such as hypoxia, metabolism, redox reaction, PI3K-AKT-mTOR pathway, TGF-β pathway, etc. (33,34). However, the exact mechanism by which cancer stemness drives tumor development is still not fully understood. In recent years, with the advancement of various omics technologies such as single-cell transcriptomics and ST, our understanding of the roles that various phenotypes play in the progression of diseases has greatly expanded (35,36). Therefore, to explore the relationship between cancer stemness and the occurrence and development of CRC, we conducted an integrated analysis using scRNA-seq and ST. Our aim is to identify novel therapeutic strategies for CRC by harnessing the insights gained from these cutting-edge approaches.

This study employed public scRNA-seq data to delve deeply into the cellular heterogeneity between tumor and normal tissues in CRC. We identified key CSC subtypes and explored their stemness characteristics as well as their interactions with other cells within the TME. Although several markers, such as CD133, CD44, CD166, and CD24 have been shown to belong to CSCs in a variety of solid tumors, but the expression of these markers is not uniform among tumor types (37,38). Therefore, in order to find the CSC group in CRC, we used CytoTRACE to find the cell group with the highest degree of stemness in the defined malignant tumor cells and defined it as Co-CSS, instead of using CRC markers for enrichment. To verify the role of Co-CSS in CRC, we intersected the maker genes of Co-CSS with the genes of TCGA-COAD, and then performed Cox regression analysis on the stemness-related genes to obtain 10 genes related to the prognosis of CRC. We found that TRIP6 was highly expressed in CRC and was associated with poor prognosis, indicating that Co-CSS plays a very important role in tumor progression, which is consistent with previous studies (25,39).

In recent years, the TME has received increasing attention because it plays a key role in tumor immunosuppression, progression and metastasis, local drug resistance and targeted therapy response (40,41). TME includes tumor cells, rich and diverse immune cells, fibroblasts, endothelial cells and other cell types that vary depending on the tissue, as well as a variety of signaling molecules (42,43). In the TME, cells do not work alone. They can communicate between cells through cell-to-cell contact or by secreting signal ligands to achieve ligand-receptor binding (44,45). In tumor tissue, Co-CSS can maintain a certain degree of stemness, which is related to the TME in which it is located (46,47). Therefore, we focused on the interaction between Co-CSS and other cells in CRC. We found that Co-CSS strongly interacts with endothelial cells and fibroblasts, and there are numerous interactions. Regarding the signaling interactions between Co-CSS and other cells, fibroblasts regulate Co-CSS through the WNT signaling pathway, while endothelial cells may affect Co-CSS through the IGF pathway. Previous studies have shown that fibroblasts and endothelial cells can promote cancer stemness by secreting relevant ligands that act on Co-CSS receptors (34,48,49). These findings suggest that fibroblasts and endothelial cells play a significant role in supporting the maintenance and development of Co-CSS stemness.

Although scRNA-seq can analyze cell heterogeneity and reveal the relationship between cell clusters in the TME, ST is still needed to understand the spatial heterogeneity within the tumor (50). In this study, we found that Co-CSS and fibroblasts co-localized, and fibroblasts may regulate the progression of stem cells through the WNT signaling pathway. Through MIA, we integrated ST data and scRNA-seq data to determine Co-CSS. Compared with the main distribution areas of other cells, Co-CSS was mainly distributed in cluster 9, while fibroblasts were mainly distributed in cluster 9. We found that Co-CSS and fibroblasts co-localized. For stem cells, cell-to-cell contact can position stem cells in an ecological niche that is conducive to maintaining their stemness, maintaining their phenotype and exerting their functions (51,52). Obviously, fibroblasts play a very important role in the niche anchoring of Co-CSS. The WNT signaling pathway plays an important role in tumor development and has been shown to regulate CSCs and induce stemness in colon cancer and other cancers (53,54). We found that WNT was activated in the cluster 9 region and the WNT signaling pathway showed high activity scores in fibroblasts, indicating that fibroblasts may regulate Co-CSS through WNT, which is consistent with our results in scRNA-seq.

Although the current National Comprehensive Cancer Network (NCCN) treatment guidelines for CRC have included targets such as MSI, BRAF, KRAS, NRAS, RAS, HER2, and NTRK in genetic testing, we still need to find more new targets to provide more strategies for the treatment of CRC (55,56). Because the cluster 9 region is rich in tumor stem cells and fibroblasts that promote the progression of Co-CSS, we hope to find meaningful clinical treatment targets from marker genes in the cluster 9 region. Among them, we found that NPDC1 and NSMF (NELF) were associated with poor prognosis. Previous studies have identified NPDC1 as a prognostic marker in CRC patients, with its expression indicating poor clinical outcomes. Meanwhile, NSMF plays a critical role in promoting cancer stemness (57,58). In addition, we found that NPDC1 was highly expressed in tumor cells with high stemness and high CNV, and its expression showed a trend of first increasing and then decreasing during the development of tumor cells, which indicates that NPDC1 may play a driving role in the differentiation of tumor cells. Therefore, we speculate that targeting NPDC1 can inhibit cancer stemness and block the progression of tumor cells, and may be a potential therapeutic target for CRC, which requires further research to test its possibility.

Limitations

Although this study offers valuable analysis and insights, it has not been experimentally verified. Future research could further validate its conclusions by incorporating experimental data, increasing the sample size, and updating the data sources, thereby enhancing both the depth and breadth of the study.

Conclusions

This study offers insights into the role of Co-CSS in CRC through ST and gene expression profiling, uncovering novel markers and pathways. It redefines the interplay between Co-CSS and the TME, offering a foundation for future mechanistic studies and therapeutic interventions aimed at targeting Co-CSS to improve CRC treatment.

Acknowledgments

None.

Footnote

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-586/prf

Funding: This study was supported by Natural Science Foundation of Shanghai Municipality (16ZR1400800).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-586/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Bray F, Laversanne M, Sung H, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2024;74:229-63. [Crossref] [PubMed]
2. Dekker E, Tanis PJ, Vleugels JLA, et al. Colorectal cancer. Lancet 2019;394:1467-80. [Crossref] [PubMed]
3. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell 2011;144:646-74. [Crossref] [PubMed]
4. Marusyk A, Almendro V, Polyak K. Intra-tumour heterogeneity: a looking glass for cancer? Nat Rev Cancer 2012;12:323-34. [Crossref] [PubMed]
5. Batlle E, Clevers H. Cancer stem cells revisited. Nat Med 2017;23:1124-34. [Crossref] [PubMed]
6. Reya T, Clevers H. Wnt signalling in stem cells and cancer. Nature 2005;434:843-50. [Crossref] [PubMed]
7. Dean M, Fojo T, Bates S. Tumour stem cells and drug resistance. Nat Rev Cancer 2005;5:275-84. [Crossref] [PubMed]
8. Meacham CE, Morrison SJ. Tumour heterogeneity and cancer cell plasticity. Nature 2013;501:328-37. [Crossref] [PubMed]
9. Marjanovic ND, Weinberg RA, Chaffer CL. Cell plasticity and heterogeneity in cancer. Clin Chem 2013;59:168-79. [Crossref] [PubMed]
10. Mzoughi S, Schwarz M, Wang X, et al. Oncofetal reprogramming drives phenotypic plasticity in WNT-dependent colorectal cancer. Nat Genet 2025;57:402-12. [Crossref] [PubMed]
11. Saygin C, Matei D, Majeti R, et al. Targeting Cancer Stemness in the Clinic: From Hype to Hope. Cell Stem Cell 2019;24:25-40. [Crossref] [PubMed]
12. Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet 2009;10:57-63. [Crossref] [PubMed]
13. Macosko EZ, Basu A, Satija R, et al. Highly Parallel Genome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets. Cell 2015;161:1202-14. [Crossref] [PubMed]
14. Ståhl PL, Salmén F, Vickovic S, et al. Visualization and analysis of gene expression in tissue sections by spatial transcriptomics. Science 2016;353:78-82. [Crossref] [PubMed]
15. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 2014;15:550. [Crossref] [PubMed]
16. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284-7. [Crossref] [PubMed]
17. Brinkman EK, Chen T, Amendola M, et al. Easy quantitative assessment of genome editing by sequence trace decomposition. Nucleic Acids Res 2014;42:e168. [Crossref] [PubMed]
18. Cokelaer T, Chen E, Iorio F, et al. GDSCTools for mining pharmacogenomic interactions in cancer. Bioinformatics 2018;34:1226-8. [Crossref] [PubMed]
19. Cable DM, Murray E, Zou LS, et al. Robust decomposition of cell type mixtures in spatial transcriptomics. Nat Biotechnol 2022;40:517-26. [Crossref] [PubMed]
20. Becht E, McInnes L, Healy J, et al. Dimensionality reduction for visualizing single-cell data using UMAP. Nat Biotechnol 2018; Epub ahead of print. [Crossref]
21. Moncada R, Barkley D, Wagner F, et al. Integrating microarray-based spatial transcriptomics and single-cell RNA-seq reveals tissue architecture in pancreatic ductal adenocarcinomas. Nat Biotechnol 2020;38:333-42. [Crossref] [PubMed]
22. Qiu X, Mao Q, Tang Y, et al. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods 2017;14:979-82. [Crossref] [PubMed]
23. Jin S, Guerrero-Juarez CF, Zhang L, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun 2021;12:1088. [Crossref] [PubMed]
24. Zhang D, Tang DG, Rycaj K. Cancer stem cells: Regulation programs, immunological properties and immunotherapy. Semin Cancer Biol 2018;52:94-106. [Crossref] [PubMed]
25. Gou H, Wong CC, Chen H, et al. TRIP6 disrupts tight junctions to promote metastasis and drug resistance and is a therapeutic target in colorectal cancer. Cancer Lett 2023;578:216438. [Crossref] [PubMed]
26. Wang W, Kandimalla R, Huang H, et al. Molecular subtyping of colorectal cancer: Recent progress, new challenges and emerging opportunities. Semin Cancer Biol 2019;55:37-52. [Crossref] [PubMed]
27. Dagogo-Jack I, Shaw AT. Tumour heterogeneity and resistance to cancer therapies. Nat Rev Clin Oncol 2018;15:81-94. [Crossref] [PubMed]
28. Weissman IL, Anderson DJ, Gage F. Stem and progenitor cells: origins, phenotypes, lineage commitments, and transdifferentiations. Annu Rev Cell Dev Biol 2001;17:387-403. [Crossref] [PubMed]
29. Pattabiraman DR, Weinberg RA. Tackling the cancer stem cells - what challenges do they pose? Nat Rev Drug Discov 2014;13:497-512. [Crossref] [PubMed]
30. Singh A, Settleman J. EMT, cancer stem cells and drug resistance: an emerging axis of evil in the war on cancer. Oncogene 2010;29:4741-51. [Crossref] [PubMed]
31. Malanchi I, Santamaria-Martínez A, Susanto E, et al. Interactions between cancer stem cells and their niche govern metastatic colonization. Nature 2011;481:85-9. [Crossref] [PubMed]
32. Li F, Tiede B, Massagué J, et al. Beyond tumorigenesis: cancer stem cells in metastasis. Cell Res 2007;17:3-14. [Crossref] [PubMed]
33. Ito K, Suda T. Metabolic requirements for the maintenance of self-renewing stem cells. Nat Rev Mol Cell Biol 2014;15:243-56. [Crossref] [PubMed]
34. Paul R, Dorsey JF, Fan Y. Cell plasticity, senescence, and quiescence in cancer stem cells: Biological and therapeutic implications. Pharmacol Ther 2022;231:107985. [Crossref] [PubMed]
35. Ding J, Adiconis X, Simmons SK, et al. Systematic comparison of single-cell and single-nucleus RNA-sequencing methods. Nat Biotechnol 2020;38:737-46. [Crossref] [PubMed]
36. Marx V. Method of the Year: spatially resolved transcriptomics. Nat Methods 2021;18:9-14. [Crossref] [PubMed]
37. Prasetyanti PR, Medema JP. Intra-tumor heterogeneity from a cancer stem cell perspective. Mol Cancer 2017;16:41. [Crossref] [PubMed]
38. Visvader JE, Lindeman GJ. Cancer stem cells in solid tumours: accumulating evidence and unresolved questions. Nat Rev Cancer 2008;8:755-68. [Crossref] [PubMed]
39. Gou H, Liang JQ, Zhang L, et al. TTPAL Promotes Colorectal Tumorigenesis by Stabilizing TRIP6 to Activate Wnt/β-Catenin Signaling. Cancer Res 2019;79:3332-46. [Crossref] [PubMed]
40. Quail DF, Joyce JA. Microenvironmental regulation of tumor progression and metastasis. Nat Med 2013;19:1423-37. [Crossref] [PubMed]
41. Chen F, Zhuang X, Lin L, et al. New horizons in tumor microenvironment biology: challenges and opportunities. BMC Med 2015;13:45. [Crossref] [PubMed]
42. de Visser KE, Joyce JA. The evolving tumor microenvironment: From cancer initiation to metastatic outgrowth. Cancer Cell 2023;41:374-403. [Crossref] [PubMed]
43. Cadamuro M, Fabris L, Zhang X, et al. Tumor microenvironment and immunology of cholangiocarcinoma. Hepatoma Res 2022;8:11. [Crossref] [PubMed]
44. Armingol E, Officer A, Harismendy O, et al. Deciphering cell-cell interactions and communication from gene expression. Nat Rev Genet 2021;22:71-88. [Crossref] [PubMed]
45. Hanahan D, Coussens LM. Accessories to the crime: functions of cells recruited to the tumor microenvironment. Cancer Cell 2012;21:309-22. [Crossref] [PubMed]
46. Kise K, Kinugasa-Katayama Y, Takakura N. Tumor microenvironment for cancer stem cells. Adv Drug Deliv Rev 2016;99:197-205. [Crossref] [PubMed]
47. Su S, Chen J, Yao H, et al. CD10(+)GPR77(+) Cancer-Associated Fibroblasts Promote Cancer Formation and Chemoresistance by Sustaining Cancer Stemness. Cell 2018;172:841-856.e16. [Crossref] [PubMed]
48. Ma Z, Li X, Mao Y, et al. Interferon-dependent SLC14A1(+) cancer-associated fibroblasts promote cancer stemness via WNT5A in bladder cancer. Cancer Cell 2022;40:1550-1565.e7. [Crossref] [PubMed]
49. Zhang Z, Dong Z, Lauxen IS, et al. Endothelial cell-secreted EGF induces epithelial to mesenchymal transition and endows head and neck cancer cells with stem-like phenotype. Cancer Res 2014;74:2869-81. [Crossref] [PubMed]
50. Zormpas E, Queen R, Comber A, et al. Mapping the transcriptome: Realizing the full potential of spatial data analysis. Cell 2023;186:5677-89. [Crossref] [PubMed]
51. Sneddon JB, Werb Z. Location, location, location: the cancer stem cell niche. Cell Stem Cell 2007;1:607-11. [Crossref] [PubMed]
52. Borovski T, De Sousa E, Melo F, Vermeulen L, et al. Cancer stem cell niche: the place to be. Cancer Res 2011;71:634-9. [Crossref] [PubMed]
53. Vermeulen L, De Sousa E, Melo F, van der Heijden M, et al. Wnt activity defines colon cancer stem cells and is regulated by the microenvironment. Nat Cell Biol 2010;12:468-76. [Crossref] [PubMed]
54. Nusse R, Clevers H. Wnt/β-Catenin Signaling, Disease, and Emerging Therapeutic Modalities. Cell 2017;169:985-99. [Crossref] [PubMed]
55. Benson AB, Venook AP, Al-Hawary MM, et al. Colon Cancer, Version 2.2021, NCCN Clinical Practice Guidelines in Oncology. J Natl Compr Canc Netw 2021;19:329-59. [Crossref] [PubMed]
56. Benson AB, Venook AP, Al-Hawary MM, et al. Rectal Cancer, Version 2.2022, NCCN Clinical Practice Guidelines in Oncology. J Natl Compr Canc Netw 2022;20:1139-67. [Crossref] [PubMed]
57. Li J, Sun Y, Cao L, et al. Correlation of NPDC1 Expression and Perineural Invasion Status with Clinicopathological Features in Patients with Colon Cancer. Int J Gen Med 2023;16:4549-63. [Crossref] [PubMed]
58. Aoki K, Nitta A, Igarashi A. NELF and PAF1C complexes are core transcriptional machineries controlling colon cancer stemness. Oncogene 2024;43:566-77. [Crossref] [PubMed]