All human tumor samples and matched adjacent normal samples were collected from patients diagnosed with CC at the Harbin Medical University Cancer Hospital in October 2020. Patients were diagnosed with adenocarcinoma by endoscopy and had not undergone any preoperative chemotherapy or radiotherapy. Tissue samples were immediately stored in liquid nitrogen after resection. This study was approved by the Ethical Committees of Harbin Medical University Cancer Hospital, and patients signed an informed consent form.

In this study, we downloaded 629 CC RNA-sequencing data from two high-throughput sequencing platforms; 473 cases from TCGA and 156 cases from GEO (GSE103479) (http://www.ncbi.nlm.nih.gov/geo/). The above data included information about each patient’s somatic mutation, primary tumor site, clinical stage, and survival status. The tumor was characterized as RCC when it was located in the cecum, ascending colon, and hepatic flexure of the colon. The primary location of the tumor was the splenic flexure of the colon, descending colon, sigmoid and recto-sigmoid junction and were defined as LCC. After excluding patients with incomplete survival data, 444 samples were included in the study; TCGA 322 and GEO 122. We obtained the normalized matrix files from GEO for the microarray data. For data from TCGA, we downloaded the RNA-sequencing data (FPKM value) of gene expression, which was then converted into transcripts per kilobase million (TPM) values for combined analysis. To adjust the batch effect caused by non-biotechnology deviation, the “Combat” function of the R package “SVA” was used.

iTRAQ was used to examine six pairs of left- and right- colon cancer and matched normal samples, and 243 DEPs (|log2foldchange| > 0.5, p-adj < 0.05) were identified (Detail iTRAQ method and protein information in Supplementary Data). Then, DEPs related mRNA expression levels were extracted from the integrated data set. The DEP mRNA was used to distinguish between tumor types. The hierarchical agglomerative cluster was used to cluster all tumor samples in the R package “ConsensusClusterPlus” (Wilkerson and Hayes, 2010). Cluster count and membership were identified through unsupervised analysis using stability evidence. To ensure the cluster’s stability, the analysis was repeated 1,000 times.

To compare the biological processes of DEP mRNA in distinct tumor types, we utilized the R package “GSVA” to perform enrichment difference analysis. GSVA uses a nonparametric and unsupervised method to estimate the variation pathway and enrichment of biological processes across several expression datasets. Following that, we download Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)-related databases from the MSigDB database (http://softwar.broadinstiture.org/gsea/msigdb/) to finish the GSVA analysis. The heat map displays the significantly distinct pathway in the analysis results (p-adj < 0.05).

To estimate immunocyte infiltration of samples, the R package “GSVA” was used to analyze single-sample gene-set enrichment analysis (ssGSEA). We acquired information on immune cell marker gene expression from Charoentong’s Charoentong et al., 2017 research and calculated the enrichment coefficients, which showed the relative abundance of immune cells in samples. Finally, we compared immune infiltration patterns across various DEP mRNA clusters.

R package “Limma” was used to identify DEGs (p-adj < 0.01) in different DEP mRNA clusters (Ritchie et al., 2015). Following that, numerous sets of DEGs were obtained, and the intersection was used to extract their expression in all samples for subsequence processing. The intersecting DEGs were analyzed using the R package “clusterProfiler” (Yu et al., 2012), which included GO and KEGG functional annotations analysis.

We developed an algorithm and defined it as DEPRS, for quantification of DEP mRNA-related types in LCC and RCC. The entire process was as follows: Firstly, the intersecting DEGs were subjected to univariate Cox proportional hazard regression analysis (COX) using the R package “glmnet” (Friedman et al., 2010). The gene that had a significant effect on prognosis was standardized. Based on these DEP mRNA-associated DEGs, the unsupervised clustering method was used to cluster all patients for further analysis. Meanwhile, we performed principal component analysis (PCA) to identify the key components of these genes to construct a DEP mRNA-related gene signature. Principal component 1 and principal component 2 were chosen as the signature scores. Finally, we used an approach comparable to the gene expression level index to calculate each patients’ DEPRS: DEPRS= ∑PCA1i+∑PCA2i (i is the expression of DEP mRNA-related DEGs). The best cut-off value for dividing high and low DEPRS groups for prognosis was obtained using the R package “maxstat” (Laska et al., 2012).

Immunotherapy sensitivities were estimated in the high and low DEPRS groups from two perspectives: immune checkpoint-related genes immunophenoscores (IPS) and open-label immunotherapy cohort. Various immune-related genes, including effector cells, immune suppressor cells, MHC molecules, and immunoregulatory cytokines were used to determine immunogenicity. Machine-learning algorithms can accurately estimate and quantify immunogenicity. The IPS of CC in TCGA was downloaded from the TCIA database (https://tcia.at/). Then, we compared the immunophenoscore differences between high and low DEPRS groups in immunotherapy, to predict immunotherapy sensitivities. Simultaneously, a thorough search for gene expression profiles in publicly available immunotherapy cohort of metastatic urothelial tumors (IMvigor210: http://research-pub.gene.com/IMvigor210CoreBiologies) Mariathasan et al., 2018 was performed. Data were pretreated using the R package “IMvigor210CireBiologies.” The RNA-SEQ data was filtered and normalized using the R package “edgeR,” and transformed using voom in the R package “limma.” In addition, we downloaded and organized prognostic status and therapeutic effect data. Based on the above computation, the DEPRS for each sample in this cohort was calculated and divided into high and low score groups to compare the difference in therapeutic response between the two groups.

To screen and identify DEP-related prognostic key genes in LCC and RCC, two methods were used: weighted gene co-expression network analysis (WGCNA), random forest, and artificial neural network. First, we used the combination of DEP mRNA and DEPRS to run the WGCNA. To transform the adjacency matrix (AM) to a topological overlap matrix with DEGs, the appropriate power index was selected. The higher the value of the mean connectivity, the more the network conforms to the scale-free characteristics. According to the relationship between soft threshold (power) and mean connectivity, the minimum index when the R ² of scale-free network reaches 0.8 was taken as the appropriate index. Then, a correlation between gene consensus modules and DEPRS was established, and gene significance (GS) was defined as the mediated p-value of each gene (GS = lgP) in a linear regression between gene expression and the scores. Subsequently, GS >0.6 genes in the module with the highest positive correlation coefficient with DEPRS were screened. Following that, the random Forest software package was used to classify all DEP mRNA. The parameter mtry (Optimal variable number of binary trees in the nodes) was set to 6, and the optimal number of trees included in the random forest was 49. After constructing the random forest module, the dimensional importance value (IV) was determined using the module’s decreasing accuracy method (Gini coefficient method). Genes with a length of more than two were treated as special genes to construct subsequent modules. Unsupervised hierarchical cluster analysis was performed on these specific genes and a heat map was generated to demonstrate their classification effect. Following that, we classified these genes into high and low expression groups based on their median expression and categorized them as gene scores (low expression as 0, high expression as 1). On significant variables, the R package “Neuralnet” was used to construct an artificial neural network model. The model’s hidden layer parameter was set to 5, and the outcomes of weight score multiplied by gene score were used to construct a classification model for HDEPRSG and LDEPRSG. The R package “pROC” was used to create the receiver operating characteristic curve (ROC), as well as calculate the area under the curve to confirm classification performance. Finally, the specific gene was obtained from an intersection of the particular gene and the screened gene. To repeatedly validate the important gene at the protein level, we downloaded the proteomics cohort in the TCGA COAD sample (including 29 normal samples and 64 tumor samples) from The Clinical Proteomic Tumor Analysis Consortium (CPTAC) (https://proteomics.cancer.gov/programs/cptac). Additionally, we investigated the prognostic significance and expression of key genes.

All the statistical analyses were performed using R-version 4.0.5. The Wilcoxon test and Kruskal-Wallis test were used to compare two groups, and more than two groups respectively. The prognosis curve was constructed using the Kaplan-Meier plotter, and the log-rank test was used to determine whether there was a statistically significant difference in prognosis. Spearman correlation coefficient was used to measure the relationship between variables. The mutation gene status in different groups was demonstrated using the R package “maftool” (Mayakonda et al., 2018). p < 0.05 was considered statistically significant.

The CC samples and matched normal samples were subjected to iTRAQ analysis. Then, 243 DEPs (|log2foldchange| > 0.5, p-adj < 0.05) were identified and the corresponding mRNA data were extracted from the integrated dataset. The network was used to depict the comprehensive DEP mRNA landscape. Univariate COX analysis was used to investigate all of the DEP mRNAs, and the interaction between them was demonstrated (Figure 2A). As shown in the figure, there was an explicit negative association between the favorable factor and risk factor. Next, Following that, all CC samples were clustered using the R package “ConsensusClusterPlus” based on 243 DEP mRNA, and 3 clusters were identified (Figure 2B).

A Kaplan-Meier (K-M) curve was used to compare the prognosis of various DEP mRNA clusters. We found that DEP mRNA cluster B had a significantly worse prognosis than the other two clusters (p = 0.031) (Figure 2C). Then, DEP mRNA expression in the three clusters was shown using a heat map, and the expression of the 3 clusters was significantly different (Figure 2D). The KEGG-related GSVA was used to investigate the biological function of these DEP mRNA clusters. When cluster A was compared to cluster B, we observed that DEP mRNA cluster B was significantly more abundant in carcinogenic activation pathways, such as small cell lung cancer, renal cell carcinoma, adherens junction, and TGF-β pathway (Figure 2E). Then, GSVA analysis between clusters A and C revealed that cluster A enriched for a variety of pathways not previously associated with it, including ECM receptor interaction, chemokine pathway, and calcium signaling pathway (Figure 2F). In comparison to clustering C, DEP mRNA cluster B enriched for a variety of carcinogenic activation pathways including the MAPK signaling pathway, focal adhesion, etc. (Figure 2G). As a result of the aforementioned finding, DEP mRNA cluster B was found to have the worst prognosis among three clusters and an enrichment analysis with numerous carcinogenic pathways.

The PAC clustering method was used to confirm the DEP mRNA clustering result previously obtained. As a result of the PCA grouping, three distinct DEP mRNA groups were identified, indicating the consistency and accuracy of the test (Figure 3A). The DEP mRNA cluster B showed a poor prognosis and was associated with the cancer signaling pathway. The ssGSEA was used to evaluate each sample to determine the link between the three clusters and immune infiltration. Surprisingly, an investigation of immunological infiltration revealed that DEP mRNA cluster B was significantly abundant in immune cells such as MDSC, activated B cell, activated CD4 T cell, and natural killer cell (Figure 3B). However, when compared to clusters A and C, cluster B showed no advantage in terms of survival time. To further investigate the differences in the immunological microenvironment of DEP mRNA, the R package “Limma” was used to identify DEG between different DEP mRNA clusters. There were 1,577 overlapping DEGs (Figure 3C). Additionally, using the R package “clusterProfiler” the DEGs were subjected to GO and KEGG analysis. In GO enrichment, the overlapped DEGs were enriched in immune response (Figure 3D), while in KEGG they were enriched in inflammation-related pathways and carcinogenic pathways (Figure 3E).

To conduct in-depth research on the various expression types of DEP that contribute to the differences in the immune microenvironment, the TCGA-COAD sample was re-divided into three gene clusters based on overlapping DEGs (Figure 4A). To determine the survival difference between the three gene clusters, a prognostic analysis was performed. Although the results presented on the survival curve weren’t perfect, but the novel gene clusters were shown the significant difference of prognosis among different clusters (p = 0.04) (Figure 4B). Following that, the DEGs expression in distinct gene clusters and DEP mRNA clusters was visualized using the heat map (Figure 4C). Taking into account the complexity and variability of individual differential protein expression patterns and the subsequent identification of important genes, we developed a novel algorithm called DEPRS to quantify the DEP mRNA expression in individual patients. The optimal cut-off value was determined using the R package “maxstat,” and patients were divided into high DEPRS and low DEPRS groups (HDEPRSG and LDEPRSG). As shown in the prognostic analysis, LDEPRSG had a better prognosis than HDEPRSG (p < 0.001) (Figure 4D). Following that, a Sankey diagram illustrating the distribution of patients with tumor sites, DEP mRNA cluster, gene clusters, and DEPRS was displayed (Figure 4E). As indicated in Figure 4E, CC was divided into LCC and RCC, which were further divided into three DEP mRNA clusters. Following that, the DEP mRNA clusters were stratified into three gene clusters. Surprisingly, samples in the DEP mRNA cluster B and gene cluster B were classified as having a high DEPRS score, indicating a poor prognosis. A portion of DEP mRNA cluster B occupied a section of gene cluster B, but partial samples in DEP mRNA cluster A were retained in gene cluster A distributed in low DEPRS, which was associated with a better prognosis. It was indirectly demonstrated that multiple clustering modes produced consistent results. Following the above result, the difference between LCC and RCC wasn’t distinct that means the DEPRS model own the special advantage to ignore tumor site.

To further illustrate the relationship between DEPRS and the preceding two cluster modes, we examined the correlation between the obtained clusters and DEPRS. There was a significant difference in DEPRS between various DEP mRNA clusters (Figure 5A). The DEP mRNA cluster B had had a much higher median DEPRES value than the other two clusters. Then, when compared to other clusters, gene cluster C had significantly increased DEPRS (Figure 5B). Because the DEG was shown to be enriched in immune-related pathways, we examined the correlation between TMB and DEPES. There was a significant correlation between TMB and DEPRS using association analysis (Coefficient: R = 0.018, p = 0.0025) (Figure 5C). With an increase in DEPRS, the distribution of gene clusters followed a significant rule. TMB showed an increasing tendency as DEPRS increased, indicating that patients with high TMB have a poor prognosis. According to the preceding data, gene cluster A had a low DEPRS and a low TMB, indicating a better prognosis. Additionally, we examined the effect of TMB and DEPRS integration on prognosis. All samples were stratified into high- and low- TMB subgroups. K-M curves were used to analyze the combined effect of TMB and DEPRS on prognosis. Within the same DEPRS status, the group with a high TMB had a worse outcome than the group with a low TMB. Nonetheless, TMB status did not affect the DEPRS predictive ability of prognosis, such that patients with high DEPRS invariably had a poor prognosis (Figure 5D). We also created a correlation map to visualize the DEPRS and immune cell interaction in TME, owing to the relationship between TMB and DEPRS. It is demonstrated unequivocally that the characteristics of DEPRS were highly correlated with high levels of immunocyte infiltration (Figure 5E). We aimed to determine the association between MSI status and DEPRS because it was critical for immunotherapy sensitivity. MSI-H differed significantly from MSS and MSI-L in this study, indicating that MSI-H was associated with a high DEPRS (Figure 6A). Consistent with this finding, 23% of MSI-H patients had high DEPRS compared with 6% with low DEPRS (Figure 6B). In addition, we compared the differences in somatic variation driver genes between individuals with high and low DEPRS. The top 20 driver genes exhibiting the highest mutation frequency were selected (Figures 6C,D). The mutation rate of the majority of driver genes was higher in the HDEPRSG than in the LDEPRSG.

The purpose of this study was to compare immunotherapy sensitivity parameters in two DEPRS groups. For immune checkpoints, PD-1, PDCD1LG2, CD274, CTLA4, HAVCR2, and LAG3 were significantly expressed in HDEPRSG than LDEPRSG (p < 0.05) (Figures 7A–F). Therefore, the parameters indicated that different groups responded differently to immunotherapy. Additionally, we demonstrated the stability of DEPRS as a predictor of immunotherapeutic efficacy in a public immunotherapy cohort. The relationship between DEPRS and treatment outcome was shown in Figure 7G. The results showed that the stable disease (SD)/progression disease (PD) proportion was significantly higher in the high DEPRS than in the low DEPRS group (79 vs 67%) and that SD/PD patients had a higher DEPRS than the complete response (CR)/partial response (PR) group (p =0.071) (Figure 7H). On the other hand, the relationship between IPS and DEPRS was used to estimate the predictive potential of DEPRS. The IPS, which measures immunogenicity, was a significant difference between DEPRS groups, with the IPS in the high DEPRS group being lower than the IPS in the low DEPRS (p = 0.0044) (Figure 7I). Based on these findings, we hypothesized that DEPRS may possess the ability to predict prognosis and immunotherapy efficacy.

The random forest module was constructed to identify key genes based on DEPRS, and Figure 8A depicts the relationship between reference model error and the number of decision trees (Figure 8A). A total of 49 decision trees were selected since the error rate was the lowest and the most relatively stable. The IV > 2 genes (HTRA3, S100A8, NNMT, and FBLN2) were chosen for further analysis using the Gini coefficient method (Figure 8B). Unsupervised hierarchical cluster analysis was used to study these specific genes, and the resulting heat map was sued to illustrate the relationship between their expression and DEPRS status (Figure 8C). The expression of these four genes was significantly increased in the high DEPRS group. Next, these four genes were used to construct a neural network prediction model as illustrated in Figure 8D. The ROC curve was used to determine the sensitivity of the module, and the AUC value was 0.886 (Figure 8E). The results suggested that the neural network prediction model based on four specific key genes produced satisfactory results and that HDEPRSG was more accurately classified than LDEPRSG. On the other hand, a gene co-expression network based on DEPRS was constructed to identify the key DEP mRNA. By selecting number 5 as the appropriate soft threshold (Figure 9A), a scale-free co-expression network was constructed (Figure 9B), yielding six modules. The brown module had the highest correlation with DEPRA (coefficient = 0.83, p < 0.001) (Figure 9C). After interacting with the four special genes identified by the neural network, two key genes HTRA3 and FBLN2 were obtained. The K-M curves were used to demonstrate that the high expression of two genes was associated with a poor prognosis (Figures 9E,F). At the protein level, high FBLN2 expression was associated with a poor prognosis (p < 0.001) (Figure 9G), and FBLN2 expression was significantly higher in tumor samples than in normal samples (p < 0.001) (Figure 9H).