Comparative analyses of the prognosis, tumor immune microenvironment, and drug treatment response between left-sided and right-sided colon cancer by integrating scRNA-seq and bulk RNA-seq data

Data download

The scRNA-seq dataset GSE200997, comprising 7 normal and 15 tumor samples, was downloaded from the US National Center for Biotechnology Information’s Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The original data included 23,828 genes and 49,859 cells, of which 18,273 cells were from normal samples, and 31,586 were from tumor samples. And the cells from tumor samples were used for further analysis, of which 16,448 cells were LCC cells and 15,138 were RCC cells.

And another scRNA-seq dataset, GSE144735, was also downloaded from GEO database for verification. The dataset included 27,414 cells from 6 patients in the core and border regions, as well as in matched normal mucosa. Based on a previously reported study [7], the cecum, ascending, and hepatic belong to the right side, and the splenic, descending, sigmoid, rectum, and junction belong to the left side. Finally, 5,780 LCC cells and 2,474 RCC cells from the core region of tumor samples were selected for subsequent validation analysis.

The Cancer Genome Atlas Colon Adenocarcinoma (TCGA-COAD) bulk RNA-seq data were downloaded from the University of California-Santa Cruz Xena platform (https://xenabrowser.net/datapages/). Only samples with survival and location information were included. Similarly, based on the published literature [7], we selected 394 samples, including 149 LCC samples, 224 RCC samples, and 21 transverse sectional samples. Their detailed information is provided in Table 1. Furthermore, we also downloaded their corresponding somatic mutation profiling data.

The independent dataset GSE103479 was downloaded from GEO database (https://www.ncbi.nlm.nih.gov/geo/), which included 77 LCC and 59 RCC samples, detailed information could be seen in Supplementary Table 1.

Identification of the location-related differentially expressed genes (DEGs) based on the dataset GSE200997

The various functions in Seurat package (version=4.1.1) for the R statistical software (v.4.0.2; http://www.R-project.org) was used to analyze the scRNA-seq data [12–14]. Genes expressed in <3 cells and cells expressing <50 genes or >6,000 genes were excluded. In addition, the mitochondrial content was <20%.

The filtered data were normalized using the NormalizeData function, and the first 3000 highly variable genes were screened using the FindVariableFeatures function. Next, the features were scaled using the ScaleData function, and their dimensionality was reduced using the RunPCA function. Then, we set dim = 20 and clustered the cells using the FindNeighbors and FindClusters functions. Importantly, we conducted the t-distributed stochastic neighbor embedding and uniform manifold approximation and projection (UMAP) algorithms to reduce dimensionality further and visualize the cluster classification based on the selected top 20 principal components. In addition, the cell types of clusters were annotated using R’s SingleR package, with HumanPrimaryCellAtlasData as the reference [15]. Moreover, we identified DEGs between LCC and RCC cells in each cell type using the FindMarkers function with logfc.threshold = 0.585. Finally, we identified location-related DEGs for CC by merging the DEGs identified in all annotated cell types. Additionally, location-related DEGs were subjected to gene ontology (GO) functional and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses using R’s clusterProfiler package [16].

Prognostic LCC and RCC signature identification and validation

The tumor samples containing expression data and survival information in the TCGA-COAD cohort were used to establish the prognostic model. Finally, 123 LCC samples and 185 RCC samples were included. Kaplan–Meier survival and univariate Cox regression analyses were performed on the location-related DEGs in the LCC and RCC expression profile data with a cutoff criterion of P < 0.05. The genes with maximum prognostic value were selected for least absolute shrinkage and selection operator (LASSO) regression analysis. Next, we established separate risk score models predicting LCC and RCC patient prognosis using a multivariate Cox regression analysis. Then, each patient’s risk score was calculated according to the following formula:

$$\begin{array}{l}risk score\\ ={\sum }_i^{n}coefficient of \chi i\\ \times scaled expression value of \chi i\end{array}$$

where “χi” represents the current signature, “i” represents the location of the current signature, and “n” represents the number of the whole signatures. Moreover, the prognostic signatures for LCC and RCC patients were validated using an independent cohort GSE103479.

Additionally, we investigated the expression of these prognostic genes in different LCC and RCC cell types of scRNA-seq datasets (GSE200997 and GSE14473), and compared the expression differences in LCC and RCC samples of bulk RNA-seq datasets (TCGA-COAD cohort and GSE103479).

The CMS analysis

We predicted the CMS group of each tumor sample of TCGA-COAD cohort using R’s CMScaller package [17], the samples with an FDR > 0.05 were filtered out. Then, we calculated the proportion of LCC and RCC samples in each CMS group, and compared the relationship between each CMS group and OS. Moreover, the expression differences of the prognostic markers among CMS groups were compared using Kruskal–Wallis test.

Statistical analysis

We compared the overall survival (OS) status between LCC and RCC patients in the TCGA-COAD cohort and GSE103479 cohort by using the Kruskal–Wallis test. Additionally, receiver operating characteristic (ROC) and Kaplan–Meier curves were plotted to assess the predictive capabilities of the established risk score prognostic models. The most significant difference between true and false positive points on the ROC curve was selected as the best critical value for grouping patients. Moreover, we evaluated the performance of prognostic models in predicting tumor staging and microsatellite status for LCC and RCC patients. All results with P < 0.05 were considered statistically significant.

Exploring TIM characteristics in LCC and RCC patients

Consistent with the samples used in constructing LCC and RCC prognosis models, the CIBERSORT algorithm was used to estimate the proportion of 22 tumor infiltrating immune cell types in each sample based on expression profiling dataset [18]. Differences in the immune landscape between the LCC and RCC patients were assessed using an unpaired t-test. The cor.test function in R was used to assess the correlations between the estimated proportions of immune cell types and the risk score of each sample in LCC and RCC patients. Moreover, we compared the mRNA levels of immune checkpoint proteins and their ligands between LCC and RCC patients and high-risk and low-risk patients.

Additionally, tumor samples containing somatic mutation data and survival information in the TCGA-COAD cohort were used to explore the mutation characteristics of LCC and RCC samples. Finally, 102 LCC samples and 166 RCC samples were included. First, we analyzed and visualized LCC and RCC patients’ somatic mutation profiles using the maftools R package [19]. Genes with significant mutation differences between LCC and RCC patients were assessed for their association with OS. Moreover, the fraction of affected oncogenic pathways was predicted using the OncogenicPathways function in the maftools R package. Next, we calculated the tumor mutation burden (TMB) value and visualized the LCC and RCC patient mutation profiling dataset. Then, TMB values and mutation profiles were compared between LCC and RCC patients. The optimal cutoff value of TMB was determined by the surv_cutpoint algorithm of survivial R package and then all samples were categorized into TMB-high group and TMB-low group. Subsequently, Kaplan–Meier curves were plotted to visualize the association between TMB values and OS in LCC and RCC patients. Finally, we calculated the median value of the sample’s TMB, and the samples with TMB values higher than this median value were classified as TMB-high group, while the samples with TMB values lower than this median value were classified as TMB-low group. Subsequently, the relationship between TMB and OS were further analyzed. All results with P < 0.05 were considered statistically significant.

Prediction and comparison of drug sensitivity in LCC and RCC patients

Consistent with the samples used in constructing LCC and RCC prognosis models, and based on the expression profile dataset of these samples, we predicted the half-maximal inhibitory concentration (IC₅₀) of cancer drugs in LCC and RCC patients using the oncoPredict R package [20] and selected drugs with an average IC₅₀ <5 and related them to CC treatment. Next, we compared responses to these drugs in LCC and RCC patients and high-risk and low-risk LCC and RCC patients.

Data availability statement

The TCGA-COAD cohort was available in the UCSC Xena (https://xenabrowser.net/datapages/). The scRNA-seq datasets GSE200997, GSE14473 and bulk RNA-seq data GSE103479 were download from NCBI-GEO database (https://www.ncbi.nlm.nih.gov/geo/).

Location-related DEGs identification in different cell types

Figure 1 highlights the whole analysis workflow. After filtering according to the selection criteria, there were 31,586 tumor cells in the scRNA-seq dataset (GSE200997), of which 16,448 were from LCC patients and 15,138 were from RCC patients. The quality control data of scRNA-seq data, such as the range of RNA features, counts, and mitochondrial gene expression percentages for each cell, were shown in Figure 2A. The cluster tree with a resolution range of 0 to 1.6 showed that when RNA_ Snn_ was equal to 0.6, the number of branches was minimized, making it the optimal choice for dimensionality reduction (Figure 2B). Then, all cells were classified into 21 clusters using the UMAP algorithm and further automatically annotated into eight main cell types using the singleR R package (Figure 2C, 2D). LCC and RCC cell distributions in each annotated cell type were shown in Figure 2E. LCC and RCC DEGs statistics in each annotated cell type were provided in Table 2. More detailed DEGs information was provided in Supplementary Table 2. Finally, 690 location-related DEGs were obtained by merging the DEGs identified from each cell type. Functional and pathway enrichment analyses were performed on the location-related DEGs using the clusterProfiler R package. KEGG analysis indicated that the location-related DEGs were significantly enriched for terms associated with the interleukin-17 and tumor necrosis factor signaling pathways and apoptosis (Figure 2F). GO analysis indicated they had functions related to cytoplasmic translation, positive cell activation regulation, and cell-cell adhesion regulation (Figure 2G).

Construction and validation of prognostic models for LCC or RCC patients

Based on TCGA-COAD cohort, we observed that the survival time was significantly longer in LCC patients than in RCC patients (P = 0.038; Supplementary Figure 1A). However, survival time did not differ significantly between patients with CC in the transverse section and LCC (P = 0.830) or RCC (P = 0.530) patients. In the GSE103479 cohort, the prognosis of LCC patients was also significantly better than that of RCC patients (P =0.043, Supplementary Figure 1B). In TCGA-COAD cohort, the relationship between CMS and OS was not significant (P >0.05, Supplementary Figure 1C). In GSE103479 cohort, CMS3 had the best OS, while CMS1 had the worst OS (Supplementary Figure 1D).

The Kaplan–Meier survival analysis identified 17 genes significantly associated with LCC patient prognosis and 20 genes significantly associated with RCC patient prognosis (P < 0.05; Supplementary Table 3). Univariate Cox regression analysis indicated that five were independent indicators for LCC patients and eight for RCC patients (Table 3). After further shrinkage through LASSO regression, we selected five genes for LCC patients and four for RCC patients. The coefficient of each prognostic signature was shown in Supplementary Figure 1E, 1F. The Kaplan-Meier curves of those genes were shown in Supplementary Figure 2.

In the prognostic model for LCC patients, two genes (FosB proto-oncogene AP-1 transcription factor subunit [FOSB] and chromosome 11 open reading frame 96 [C11orf96]) were risk factors (hazard ratio [HR] >1), while three (ribosomal protein L35 [RPL35], regenerating family member 1 alpha [REG1A], and tescalcin [TESC]) were protective factors (HR <1; Figure 3A). The heatmap showed that the risk factors were downregulated in the low-risk group and upregulated in the high-risk group (Figure 3B). In contrast, the protective factors showed the opposite pattern. The scatter diagram in Figure 3B indicated that the OS of patients was better in the low-risk group than in the high-risk group, consistent with the results in Figure 3C (P = 9.7×10^-4). The areas under the ROC curve (AUC) for the prognostic model at 3-, 4-, and 5-year OS are 0.721, 0.844, and 0.926, respectively (Figure 3D).

There were four genes in the prognostic model for RCC patients. One (heat shock protein family A member 1A [HSPA1A]) was a risk factor (HR >1), and three (cluster of differentiation 69 [CD69], growth differentiation factor 15 [GDF15], and galectin 2 [LGALS2]) were protective factors (HR <1; Figure 3E). Consistent with the prognostic model for LCC patients, we found that the expression of risk factor was low in low-risk patients and high in high-risk patients. In contrast, the protective factors showed the opposite pattern (Figure 3F). In addition, patients in the low-risk group had better prognoses than those in the high-risk group (P < 1.0×10^-4; Figure 3G). The model’s AUCs were 0.806 for 3-year OS, 0.816 for 4-year OS, and 0.836 for 5-year OS (Figure 3H).

We further validated the effectiveness of the identified prognostic markers based on GSE103479 cohort. Consistent with the results obtained in the prognosis model based on TCGA-COAD cohort, in the LCC prognosis model, the prognosis of patients with low-risk patients was significantly better than that of patients with low-risk score (P = 0.014, Supplementary Figure 3A), and the AUC for the prognostic model at 3-, 4-, and 5-year OS are 0.680, 0.688, and 0.615, respectively (Supplementary Figure 3B). In the RCC prognosis model, the low-risk patients’ OS was significantly longer than high-risk patients’ OS (P =0.012, Supplementary Figure 3C), and the model’s AUCs were 0.615 for 3-year OS, 0.710 for 4-year OS, and 0.639 for 5-year OS (Supplementary Figure 3D).

We also evaluated the performance of prognostic model in predicting tumor staging and microsatellite status based on TCGA-COAD cohort. In the RCC prognostic model, higher risk scores were significantly associated with greater metastasis (P = 2.6×10^-5; Figure 4A), N stage (P = 0.012; Figure 4B), advanced pathological stage (P = 1.4×10^-4; Figure 4C), and more stable microsatellite status (P = 0.007; Figure 4D) but not T stage (P = 0.200; Figure 4E). However, in the LCC prognostic model, risk scores were not significantly associated with tumor stage and microsatellite status (P >0.05; Supplementary Figure 3E–3I), potentially reflecting the small sample size at a certain level (<30; Table 1), leading to no statistical significance.

The expression of the candidate prognostic genes in scRNA-seq and bulk RNAseq datasets

Based on scRNA-seq dataset GSE200997, the differential expression of the marker genes between LCC cells and RCC cells in each annotated cell type was shown in Table 4. We observed that HSPA1A was expressed in all RCC cell types, CD69 was mainly expressed in immune-related cell types, including T, B, and natural killer cells; while GDF15 was mainly expressed in epithelial cells, and LGALS2 was expressed in epithelial cells and macrophages (Figure 4F). FOSB and RPL35 were expressed in all LCC cell types; while REG1A was expressed in epithelial cells; TESC was mainly expressed in T, B, and epithelial cells; and C11orf96 was mainly expressed in fibroblasts and tissue stem cells (Figure 4G).

Based on another scRNA-seq dataset GSE14473, we validated the expression differences of the candidate markers between LCC and RCC cells, and we found that all candidate markers were significantly differentially expressed except for C11orf96 (Table 5). Moreover, we also observed the expression patterns of the candidate markers in the cells types was similar to those of the dataset GSE200997 (Supplementary Figure 4A, 4B).

Additionally, we investigated the expression differences of the candidate biomarkers in two bulk RNA-seq dataset. Based on the TCGA-COAD cohort, we observed that only REG1A and CD69 were significantly differential expressed between LCC and RCC patients (P <0.05, Figure 4H). In the dataset GSE103479, all the candidate markers were not significantly differential expressed between LCC and RCC patients (P > 0.05, Supplementary Figure 4C).

Based on the two bulk RNA-seq datasets, we exhibited the proportions of LCC and RCC patients in each CMS group (Figure 4I, 4J). We observed that FOSB and GDF15 were significantly expressed among CMS groups in both the two bulk RNA-seq cohorts (P <0.05, Figure 4K and Supplementary Figure 4D).

TME differences between LCC and RCC patients

Based on the TCGA-COAD cohort, the proportions of 22 immune cell types were estimated in LCC and RCC patients using the CIBERSORT algorithm. The immune-infiltrating profiles of LCC and RCC patients are shown in Supplementary Figure 3A, 3B. When we compared the composition of immune cell types between LCC and RCC patients, we found significant differences in CD8 and regulatory (Tregs) T cells and M0, M1, and M2 macrophages (Figure 5A).

In the LCC prognostic model, the immune-infiltrating degree of M2 macrophages was positively correlated with risk score (cor = 0.181, P = 0.045; Figure 5B), while the immune-infiltrating degree of activated CD4 memory T cells was negatively correlated (cor = −0.185, P = 0.041; Figure 5C). In the RCC prognostic model, the immune-infiltrating degree of Tregs was positively correlated with risk score (cor = 0.229, P = 0.002; Figure 5D), while the immune-infiltrating degree of resting CD4 memory T cells was negatively correlated (cor = −0.332, P = 4.0×10^-6; Figure 5E).

Additionally, we plotted and compared the mutation profiles of LCC and RCC patients using the maftools R package. While the top 10 mutated genes in LCC and RCC patients were similar, their ranking differed. The median number of mutations in LCC patients was higher than in RCC patients. In LCC and RCC patients, the most common variant type was single nucleotide polymorphisms (SNP), the most common SNP classification was missense, and the most common SNP class was C>T (Figure 5F, 5G). Among the top 10 mutated genes, adenomatous polyposis coli regulator of WNT signaling pathway (APC) and tumor protein p53 (TP53) mutation rates were significantly higher in LCC patients (P < 0.05). In contrast, Kirsten rat sarcoma virus proto-oncogene GTPase (KRAS), spectrin repeat containing nuclear envelope protein 1 (SYNE1), and mucin 16 cell surface associated (MUC16) the mutation rates were significantly higher in RCC patients (P < 0.05; Figure 6A). Moreover, the mutation status of TP53 was significantly associated with OS in RCC patients (HR = 1.9, P = 0.036; Figure 6B) but not in LCC patients (HR = 1.19, P = 0.695; Figure 6C). The mutation statuses of APC, KRAS, SYNE1, and MUC16 were not significantly associated with OS in LCC or RCC patients (P >0.05; Supplementary Figure 5C–5J). Figure 6D, 6E shows that the fractions of pathways and samples affected by the oncogenic pathways in RCC patients were higher than in LCC patients.

We also calculated and visualized each sample’s TMB value. The median TMB value was 1.81/Mb in LCC patients (Supplementary Figure 5K) and 2.47/Mb in RCC patients (Supplementary Figure 5L). TMB values were significantly higher in RCC patients than in LCC patients (P = 5.9×10^-8; Figure 6F). Moreover, LCC patients with lower TMB values had better OS (P = 0.019; Figure 6G). However, TMB values were not significantly associated with OS in RCC patients (P = 0.110; Figure 6H).

Prediction and comparison of the drug response in LCC and RCC patients

Figure 7A exhibited that the expression level of immune checkpoint targets in TMB-high group was significantly higher than that in TMB-low group (P < 0.05). And the expression of immune checkpoint targets in RCC patients was significantly higher than in LCC patients (P < 0.05; Figure 7B). Among LCC patients, the expression of immune checkpoint targets was significantly higher in the high-risk group than in the low-risk group, indicating that the high-risk group would be more sensitive to immunotherapy (P < 0.05, Figure 7C). Among RCC patients, cluster of differentiation 274 (CD274), cytotoxic T-lymphocyte associated protein 4 (CTLA4), and T cell immunoreceptor with Ig and ITIM domains (TIGIT) expression were significantly higher in the low-risk group than in the high-risk group. However, lymphocyte activating 3 (LAG3), programmed cell death 1 (PDCD1), and hepatitis A virus cellular receptor 2 (HAVCR2) expression did not differ significantly between the high- and low-risk groups (Figure 7D). Therefore, low-risk RCC patients are more likely to benefit from immunotherapy drugs targeting CD274, CTLA4, or TIGIT.

Additionally, we found that CC patients were sensitive to 44 drugs (average IC₅₀ <5; Supplementary Table 4), of which eight were related to CC treatment. Figure 7E shows that sensitivity to these eight drugs does not differ significantly between LCC and RCC patients (P > 0.05). In addition, drug response sensitivity did not differ significantly between high-risk and low-risk LCC patients (P > 0.05, Figure 7F). However, the average IC₅₀ values for camptothecin, teniposide, vinorelbine, and mitoxantrone were significantly lower in low-risk than in high-risk RCC patients, indicating that low-risk patients will be more sensitive to them (Figure 7G–7J). Sensitivity to the other four drugs did not differ significantly between high-risk and low-risk RCC patients (Figure 7K–7N).