Analysis of multi-omics differences in left-side and right-side colon cancer

Data collection and preprocessing

First, we downloaded colon cancer patients’ transcriptome data, clinical information, and somatic mutation data from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov). The transcriptome data were comprised of 390 colon cancer patients (172 LCC samples and 218 RCC samples), and the somatic mutation data were comprised of 329 colon cancer patients (142 LCC samples and 187 RCC samples). According to Dwertmann et al., the LCC consists of the descending colon, sigmoid colon, and splenic flexure of colon and the RCC consists of the ascending colon, cecum, and hepatic flexure of colon (Hsu et al., 2019). We used genecode.v22.annotation (https://www.gencodegenes.org/) to comment on the transcriptional data downloaded from TCGA database.

Clinical analysis in LCC and RCC

We used R to classify the data used to analyze the differences between LCC and RCC in terms of age, gender, pT, pN, pM, pStage and survival. Pearson’s chi-square (χ2) test was used to calculate differences in clinical characteristics between LCC and RCC. We compared the overall survival rates of LCC and RCC in different clinical subtypes using the survival package in R.

Immune microenvironment in LCC and RCC

We obtained immune-related gene sets with 29 immune cell types and immune-related functions from previous studies (Bedognetti et al., 2015; Ren et al., 2020; Liu et al., 2018; Maibach et al., 2020; Aran et al., 2016) to explore the differences between LCC and RCC in the immune microenvironment. We used the single sample gene set enrichment analysis (ssGSEA) algorithm to obtain the scores of 29 immune cell types and immune-related functions with the ‘GSVA’ package in R (Hänzelmann, Castelo & Guinney, 2013). We visualized the results using the pheatmap package in R (Galili et al., 2018). We used the estimate package in R to analyze the differences between LCC and RCC in the immune microenvironment to calculate the immune score, stromal score, ESTIMATE score, and tumor purity. Then we compared the differences between the two groups using the Mann–Whitney U test. We also compared the expression levels of the human leukocyte antigen (HLA) gene family and immune check-point genes, and the abundance of immune cell infiltration in LCC and RCC. We obtained the immune cell infiltration data using CIBERSORT (Chen et al., 2018).

Screening differential genes in LCC and RCC

We used the |log2 fold change(LogFC)|>1 and the false discovery rate (FDR) <0.05 with the Wilcox test to identify the differences in expression of mRNAs in LCC and RCC. The results were visualized using heatmap and volcano diagrams. Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) were used to determine the enrichment of the differential genes through the Database for Annotation, Visualization and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/) (Yang et al., 2019). The top 20 biological processes (BP) of GO enrichment analysis (Chen et al., 2017) were depicted in a circle diagram, and the top 15 KEGG pathways (Altermann & Klaenhammer, 2005) were depicted in a bubble diagram.

Screening prognostic mRNAs in LCC and RCC by univariate COX analysis

We included different genes in our study. We used the survival package in R with P < 0.005 to identify the prognostic mRNAs in LCC and RCC, respectively, using univariate COX (Tibshirani, 2009). According to univariate COX analysis, there were 22 genes associated with the prognosis of RCC, with a potential collinear relationship among them. We used the LASSO regression algorithm with a penalty term to delete genes with multicollinearity for additional analysis.

Construction and verification of the prognostic signature and validation of prognostic models

Prognostic genes related to LCC and RCC were included in our study. We set up a random number seed in order to divided LCC patients from TCGA into a training set and an internal testing set with a 1:1 ratio and established a 4-mRNA LCC prognosis model using multivariate COX regression analysis with a noose penalty (Grant, Hickey & Head, 2019; Zhou et al., 2019). We used the same method to establish a 6-mRNA RCC prognosis signature (Marisa et al., 2013; Zhang et al., 2016; Zhang et al., 2020). The samples was divided into two groups using the median risk score. We judged the efficacy of the model by plotting the Kaplan–Meier (KM) curve and receiver operating characteristic (ROC) curve (Obuchowski & Bullen, 2018). The GSE39582 data set was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39582) (Barrett et al., 2013) and was used as an external validation set. We validated the model by plotting the Kaplan–Meier (KM) curve. The GSE39582 data set included 566 colon cancer samples (342 LCC samples and 224 RCC samples) and their corresponding survival information in accordance with the GPL570 (Affymetrix Human Genome U133 Plus 2.0 Array) (Table S1). We performed an independent prognostic analysis of the risk score in the total TCGA set to further verify the model’s efficacy. The risk score was calculated as (Cho et al., 2019):

(1) $$Riskcore=\sum _{i=1}^N(Expi\times Coef)$$

with N representing the number of signature genes, Expi representing the gene expression levels, and Coef representing the estimated regression coefficient value from the Cox proportional hazards analysis.

Single gene mutation analysis in LCC and RCC

On the JAVA8 platform, we analyzed the number of variants and the length of exons for each sample using Perl scripts to calculate mutation frequency (Tabibzadeh et al., 2020). Samples were divided into two groups according to the location of colon cancer and the Mann–Whitney test (McGee, 2018) was used to compare the tumor mutation burden (TMB) difference between two groups. We used the maftools package (Mayakonda et al., 2018) for visualization and performed Fisher’s exact test in pairs between the top 25 mutated genes to analyze the mutational exclusion and co-occurrence. We also used oncoplot in R to visualize the top 30 mutated genes from the 142 LCC samples and 187 RCC samples to produce waterfall plots. Then we used the ggplot2 and boxplot packages to visualize the classification and frequency of mutation types, frequency of variant types, frequency of SNV classes, the tumor mutation burden in specific samples, and the top 10 mutated genes in LCC and RCC. The top 10 mutated genes in LCC were: APC, TP53, TTN, KRAS, MUC16, SYNE1, FAT4, RYR2, PIK3CA, and OBSCN. The top 10 mutation genes in RCC were: TTN, APC, MUC16, SYNE1, TP53, KRAS, FAT4, PIK3CA, PCLO, and ZFHX4.

Differences in clinical characteristics between LCC and RCC

The LCC and RCC data in the TCGA database and the results of the chi-square test on clinical characteristics are shown in Table 1. We classified the data by stage, T, N, M, and age after separating the data by LCC and RCC. We used the Kaplan–Meier (KM) curve of over survival (OS) to compare the survival differences of different clinical characteristics between the two groups. The results indicated that RCC had a worse prognosis than LCC, which was also seen in stages III-IV, T3-4, and N1-2 (Figs. 1A–1D). The survival rate of RCC was worse than that of LCC (Figs. 1E–1F) although there was no statistical difference between the M1 and age > 65 subgroups.

Immune microenvironment landscape between LCC and RCC

The ssGSEA algorithm showed that 29 types of immune cells and their functions were enriched in each sample. We then obtained the immune score, stromal score, ESTIMATE score, and tumor purity. The heatmap indicated that the RCC had a higher immune invasion than the LCC (Fig. 2A). Comparing the two groups’ scores, we found that only the immune scores were significantly different (Fig. 2B). We confirmed that RCC had a higher immune infiltration than LCC by further comparing the expression levels of the HLA gene family and immune checkpoint-related genes and the abundance of immune cell infiltration (Figs. 2C–2D). Previous studies have shown that the changes in HLA class I genes in colon cancer are closely related to RCC, suggesting microsatellite instability (MSI). In addition, the high expression of PD-L1 also occurs more frequently in RCC, indicating MSI (Kikuchi et al., 2019). Our results support the conclusion that RCC has more immune infiltration and is highly correlated with MSI. Therefore, this result suggests that right-side colon cancer was significantly more reactive than left-side colon cancer in immune response, which might provide new treatments for colon cancer.

Differential gene analysis between LCC and RCC

The Wilcox test was used to extract differential mRNAs to obtain 360 differential genes, which included 218 up-regulated genes in RCC and 142 up-regulated genes in LCC (Figs. 3A–3B). All of the differential genes are shown in Table S2. All of the differentially expressed genes were enriched by the biological processes of GO and KEGG pathways in the DAVID database (Tables S3 and S4). The top 20 biological processes of GO enrichment analysis were graphed in a circle diagram, while the top 15 KEGG pathways were displayed in a bubble diagram (Figs. 3C–3D). The top three biological pathways were ‘associative learning’, ‘arachidonic acid secretion’, and ‘anterior/posterior pattern specification’. The differentially expressed genes were significantly enriched in the steroid hormone biosynthesis pathway.

Univariate COX screening of prognostic genes in LCC and RCC

We screened the genes related to the prognosis of LCC and RCC using univariate Cox analysis in the LCC and RCC patients with P < 0.005. We obtained 9 genes related to prognosis in LCC and 22 genes related to prognosis in RCC (Tables 2–3). In order to avoid model overfitting, we performed LASSO regression analysis with the penalty term on RCC to solve the multicollinearity problem again by dimension reduction, and finally obtained 12 genes related to prognosis in RCC (Fig. 5B).

Construction of prognosis signature in LCC and RCC

TCGA LCC patients were divided into a training set and an internal testing set at a 1:1 ratio. Multivariate COX regression analysis with noose penalty was then used to establish a 4-mRNA LCC prognosis signature and a 6-mRNA RCC prognosis signature.

The 4-mRNA LCC prognosis signature and risk score were calculated as: C1orf105*0.458+FAM132B*1.703+TNNT1*0.130+RSPO4*0.268 (Table 4). The median risk score (0.622) in the training set was used to assign patients to the high risk or low risk group. Patients with a high risk score had significantly worse survival rates than those with low-risk scores (P = 0.046, Fig. 4A). Furthermore, the AUC of the risk score for 1-year, 2-year, 3-year, and 5-year OS were 0.751, 0.810, 0.860, and 0.904, respectively (Fig. 4B). The survival status, risk scores, and gene expression data of LCC patients in the training group are shown in Figs. 4C–4E. RSPO4, FAM132B, and TNNT1 were highly expressed in the high-risk group, while C1orf105 was not well-expressed in the high-risk group.

The risk score of the 6-mRNA RCC prognosis signature was calculated as: OFCC1*4.834+KLRG2*0.195+PAX5*0.461+SYNGR3*0.096+SLC22A31*1.232+CCDC160*0.368 (Table 5). The median risk score (0.689) in the training set was used to assign patients to the high risk or low risk group. Patients with a high-risk score was had a significantly worse survival rate than those with a low-risk score (0.012, Fig. 5A). Furthermore, the AUC of the risk score for 1-year, 2-year, 3-year, and 5-year OS were 0.776, 0.714, 0.670, and 0.792, respectively (Fig. 5C). The survival status, risk scores, and gene expression data of RCC patients in the training group are shown in Figs. 5D–5F. All six genes were highly-expressed in the high-risk group.

Validation of the prognosis signature in LCC and RCC

The prognostic accuracy of the prognosis signature was validated in three independent cohorts, including the testing set, the total TCGA data set, and the GSE39582 data set.

The OS in the high-risk group was significantly worse than that of the low-risk group in the testing set in the 4-mRNA LCC prognosis signature (P = 0.016, Fig. 6A). The predicted 1-year, 2-year, 3-year, and 5-year OS was 0.731, 0.760, 0.779, and 0.700, respectively (Fig. 6B). The total TCGA set also validated the prognostic accuracy of the signature (P = 0.001, Fig. 6C), with respective AUCs of 0.732, 0.776, 0.820, and 0.793 for 1-year, 2-year, 3-year, and 5-year OS (Fig. 6D).

The OS of the high-risk group was significantly worse than that of the low-risk group in the testing set for the 6-mRNA RCC prognosis signature (P = 0.042, Fig. 7A). The predicted 1-year, 2-year, 3-year, and 5-year OS was 0.770, 0.754, 0.689, and 0.646, respectively (Fig. 7B). The total TCGA set validated the prognostic accuracy of the signature (P = 0.002, Fig. 7C), with respective AUCs of 0.760, 0.718, 0.663, and 0.718 for 1-year, 2-year, 3-year, and 5-year OS (Fig. 7D).

The GSE39582 data set showed the same conclusion in the 4-mRNA LCC prognosis signature (P = 0.185) and the 6-mRNA RCC prognosis signature (P = 0.25) (P = 0.018, Fig. 6E; P = 0.025, Fig. 7E). The survival status, risk scores and gene expression data of LCC and RCC patients in the testing set and total TCGA set are shown in Figs. S2 and S3.

The prognosis signature confers additional prognostic power for LCC and RCC patients

Clinical characteristics, including the pStage (P < 0.001), pN (P < 0.001), pM (P = 0.004), and the risk score (P < 0.001), were closely associated with patient survival in LCC (Fig. 8A). The pStage (P < 0.001), pT (P < 0.001), pN (P < 0.001), pM (P < 0.001), age (P = 0.013), and risk score were closely associated with patient survival in RCC (Fig. 8B). Multivariate Cox regression analysis further showed that the our signature is an independent prognostic indicator for OS in LCC and RCC (Figs. 8C–8D, Tables S5–S6).

Single gene mutation landscape in LCC and RCC

The most obvious mutations, including missense mutations, were: deletion, nonsense mutation, splice site, insertion, translation start site, and nonstop mutation. The missense mutation was the most obvious. We also found that single nucleotide polymorphisms (SNP) were more frequent than insertions or deletion and the most common single nucleotide variant (SNV) was C > T (Li et al., 2011). The number of altered bases in each sample was counted and the mutation types were plotted in a box plot. The 10 most prevalent mutated genes in LCC and RCC were shown with ranked percentages (Figs. S4 and S5). The Mann–Whitney test was used to compare the TMB of LCC and RCC, and the results showed that the RCC had a higher TMB (Fig. 9A). The mutation information of each sample in LCC and RCC was graphed in a waterfall plot (Figs. 9B–9C), which showed that the mutation frequency of RCC was generally higher than that of LCC. APC, TP53, TTN, and KRAS mutations were present in both LCC and RCC (Cappell, 2008). We found that BRAF mutations were more pronounced in RCC, and APC mutations were significantly higher in LCC. The higher immune infiltration and the higher BRAF mutation in RCC suggested that RCC is closely related to MSI. The study of Lochhead P et al. showed that BRAF mutations in colorectal cancer were linked to MSI through the methylation of CIMP and MLH1 promoter methylation (Lochhead et al., 2013). This is consistent with results from previous research (Popescu et al., 2021). The high APC mutation in LCC suggests that it may be related to the inactivation of the Wnt pathway (Faux et al., 2021). LINC02418 has been shown to be a tumor driver in colon cancer (Tian et al., 2020), and whether there is an inherent relationship between LINC02418 and the mutated gene has not been investigated. The association between the mutated gene and LINC02418 may become the direction of future research.

The correlation analysis of top 25 mutated genes in LCC and RCC was conducted using the maftools package. Molecular interactions were more frequent in RCC than in LCC (Figs. 9D –9E ). In RCC, the co-occurrence of APC and KRAS and the mutually exclusive relationship of BRAF with APC and KRAS further indicated a potential relationship with CIN and CIMP (Issa, 2008).The different molecular mechanisms of RCC and LCC suggest that they may require different therapeutic approaches and prognoses.