Pathway-Specific Insights into Colorectal Cancer through Comprehensive Multi-Omics Data Integration

1. Introduction

Development in the next generation sequencing and bioinformatics have revolutionized our ability to understand the complex biological systems and their relationship with disease occurrence [1,2,3]. Multi-omics data integration is a computational approach that combines multiple layers of omics datasets. It has emerged as a significant tool to unravel biological mechanisms, and this provides comprehensive understanding of the molecular changes in disease occurrence [4].

Omics studies like metagenomics, genomics, transcriptomics, proteomics, metabolomics etc. provides remarkable insights to comprehend the disease states [5,6]. However, only one omics analysis is not enough to evaluate the diseases [7,8]. Especially cancer like complex diseases needs the combination of the omics study outputs to provide the evaluation of disease from all perspectives [9]. PCA (Principal Component Analysis), CCA (Canonical Correlation Analysis), CoIA (Co-Inertia Analysis), MFA (Multivariate Factor Analysis), PMA (Penalized Multivariate Analysis), PLS (Partial Least Square), Recursive Feature Elimination (RFE), Graph-Fused methods, Lasso Penalization methods, Bayesian classification methods, Regression methods are most preferred multi-omics integration techniques [10,11,12,13]. They generally focus on collecting the most explanatory features from different layers of omics datasets and merging all this information to obtain interactions. Despite the number of integration studies increasing dramatically, it is still difficult to obtain a golden standard. Because, handling large number of datasets containing huge number of variables is still challenging [14,15]. However, the results obtained from developed multi-omics analyzing packages provide promising results especially to explain complex cancer diseases [16,17,18]. In research conducted by Vartika Bisht and his colleagues, Recursive Feature Elimination (RFE), BART method and Bayesian classification were used for integration of omics datasets, and they discovered Kiaa1199, Cdh3, Guca2b, Lgals4, Ca7, Nrc3c2, Abcg2 and Aqp8 genes associated with cancer pathophysiology. Furthermore, they also concluded that Staphylococcus, Blautia and Roseburia were interacting with cholesteryl ester metabolite which is prominent in colorectal cancer patients [19]. All have shown remarkable cancer hallmarks such as resistance to apoptosis, angiogenesis, cell evasion, cancer pathogenicity etc. for enlightened understanding of cancer mechanism of action.

Colorectal cancer (CRC) is the third most common cancer type across the world and it demonstrates quite dynamic changes caused by environmental effects such as geography, climate, lifestyle habits (good/bad), diet, genetic information etc. [20,21]. Scientists have been trying to find potential biomarkers based of stages information and classify the CRC patients into four groups named Consensus Molecular Subtypes (CMS1, CMS2, CMS3, and CMS4) [22]. In a study published in 2017, transcriptome dataset of the 377 patients (data obtained from The Cancer Genome Atlas) were analyzed and expression level of Nek4, Rnf34, Hist3h2bb, Nudt6, Lrch4, Glb1l, Hist2h4a, Tmem79, Amıgo2, C20orf135 and Spsb3 genes were significantly differed between cancer stages (p= 0.0001) [6]. CMS groups are constructed based on immune infiltration and activation, signaling pathway information, metabolic dysregulation and stromal infiltration. Another study revealed the association of microbial changes between CMS groups. According to the results, Fusobacteria and Bacteroidetes were found to be enriched, whereas Firmicutes and Proteobacteria decreased in CMS1 subgroups. In addition to that Fusobacterium hwasookii and Porphyromonas gingavilis demonstrated the highest enrichment in CMS1 subgroup. Selenomonas and Prevotella species were discovered to indicate significant increase in CMS2 subgroup [24]. This grouping approaches are quite important to better diagnose of the patients and selecting the best treatment methodology. Studies that focused on the assessment of gut microbiome and host transcriptome provided in depth knowledge about colorectal cancer occurrence. However, it still needs the integration of omics results to provide pathway specific interactions.

In this study, we aimed to integrate microbiome, transcriptome and microbial pathway datasets obtained from CRC patients to unveil the specific interactions of host genome and microbial community. The results of the integration were consolidated by the host pathway dataset obtained from transcriptome data and the databases like genecard, metacyc, ncbi, KEGG (Kyoto Encyclopedia of Genes and Genomes) and Uniprot (Universal Protein Knowledgebase) [25,26,27,28,29]. Our analysis consisted of two approaches, one focusing on the integration of all information of the omics analysis results to provide a more comprehensive perspective not to neglect any potential interaction. The second one was more specific, it only focused on the integration of statistically significant outputs of the omics analysis. Both approaches presented non-negligible outcomes by providing latent interactions between datasets.

2. Materials and Methods

2.1. Patient Cohort and Experimental Design

Transcriptome and microbiome 16SrRNA raw datasets were obtained from recent study “Distinct gut microbiome patterns associate with consensus molecular subtypes of colorectal cancer” published by Rachel V. Purcell et. al. in 2017 [24]. Tumor samples were obtained from 34 colorectal cancer patients. 20 of them were collected from female and remaining 14 samples from male patients with the age ranged from 44-88 years (mean 74, S.D. ± 8.4). Beside age and gender, metadata also contains “Site” referring to the site of the tumor samples (one sample was obtained from rectum and others were collected from colon), “Side” referring to the position of the tumor samples (21 from right-side and 12 from left-side), “Stage” referring to the stage of the cancer patient (5 tumor tissue were stage 1, 14 were stage 2, 13 were stage 3 and 1 sample was stage 4) and “CMS” referring to which consensus molecular subtype that the tumor sample belongs (6 of them belonged to CMS1 group, 13 of them belonged to CMS2 group, 9 samples belonged to CMS3 group and 5 weren’t classified properly and were reported as unclassified (UC)). Both transcriptome and microbiome raw datasets were analyzed again including the pathway analysis. Microbial pathway abundance data was constructed because of pathway analysis of 16SrRNA analysis. Microbial pathway abundance data, transcriptome data and microbiome datasets were integrated by using two approaches. The first approach includes normalization and sparse partial least square discriminant analysis, and the second one includes normalization, feature selection and N-integration method. Dataset obtained from transcriptome analysis, genecard database and metacyc databases were used for the validation of the information obtained from integration of three datasets.

Figure 1. The chart demonstrates the experimental flow and adventages of both aproaches at the end.

Figure 1. The chart demonstrates the experimental flow and adventages of both aproaches at the end.

2.2. Microbiome and Pathway Analysis

16SrRNA analysis were performed for determining the microbial community on tumour samples. Raw datasets that were generated by Illumina Miseq platform contained 250bp long reads that cover V3-V4 region of 16SrRNA gene. FLASh (ver. 1.2.11) was used for merging forward and reverse reads. Pre-processing analysis was done by using QIIME2 (ver. 2023.5) [30]. Quality assessment and removal of chimeric sequences were performed by DADA2 (ver. 3.19) [31]. The taxonomic classification was performed using the q2-feature-classifier plugin within the QIIME2 framework, leveraging the Greengenes database as a reference. Operational Taxonomic Unit (OTU) tables were subsequently generated. Differential abundance analysis was conducted using ANCOM-BC (version 3.19).

Furthermore, the biom files derived from the QIIME2 pipeline were utilized for microbial functional pathway prediction analysis that was carried out by using PICRUST2 (ver. 2.4.1.) plugin of QIIME2 analysis platform [32]. KEGG orthology database was used for estimation of molecular functions. 16SrRNA dataset and microbial pathway abundance dataset were saved for further integration analysis.

2.3. Transcriptome and Pathway Analysis

Tumor samples were sequenced using the Illumina HiSeq platform, generating 125 bp paired-end reads. Adapter sequences were removed using Fastq-mcf (v1.1.2.537), followed by quality-based trimming of long reads, removal of reads shorter than 50 bp, and exclusion of low-quality reads using SolexaQA++ (v3.1.6). Filtered reads were aligned to the human genome reference GRCh38 using the STAR mapping tool (v2.5.2b). The resulting alignments were processed with Samtools (v1.3.1) for read merging and quantification via htseq-count.

Raw read counts were normalized to transcripts per million (TPM) using the DESeq2 package (v1.10.1) in R. RNA-Seq analysis culminated in pathway enrichment analysis using the pathfindR package (v2.1.0) in R [33]. Pathway annotations were retrieved from the KEGG database, with significant pathways filtered based on a false discovery rate (FDR) threshold of 0.05. Pathways meeting the FDR-adjusted p-value criteria were selected for downstream analyses. Transcriptome datasets were reserved for further integration analysis, whereas host (human) pathway datasets were kept for validation of the results obtained from multi-omics data integration.

2.4. Data Standardization and Pre-Processing

16SrRNA microbiome data (at each taxonomic level), microbial pathway abundance data and the transcriptome data were normalized by using centered log ratio (CLR) normalization method. Spls-DA was applied to data to fulfil our first approach in the multi-omics integration. Preprocessing was applied to each dataset for the improvement of the machine learning model accuracy and correlation coefficient obtained at the end of the first approach integration. Samples that belong to the UC group were excluded from the dataset directly by using dplyr package (ver. 1.1.4.) and tidyverse package (ver. 2.0.0.). Then, feature selection was applied based on the significance of the variables in the datasets. First, Shapiro-Wilk test was applied to check the normality of the data and Kruskal-Wallis’s test was performed . Variables that were found to be non-significant were filtered out from the datasets. We constructed new microbiome (phylum level, n=5 where n refers to number of variables), transcriptome (n=851) and microbial pathway abundance (n=26) datasets containing only significant variables at the end of the preprocessing step. CLR normalization was applied to new datasets to be used in N-Integration method.

2.5. Data Integration by Using sPLS Analysis

Sparse Partial Least Square-Discriminant Analysis (sPLS-DA) was performed by mixomics package (version 6.10) in R, for 16SrRNA microbiome, transcriptome and the microbial pathway abundance datasets as pairs to understand the association between them [34]. PLS models were tuned based on 10-fold cross validation, the parameters were used as default applied for two matrices X and Y which corresponds with omics datasets in our case. Selected explanatory variables were used for visualization of CMS groups via Principal Component Analysis (PCA) and correlation between variables that were selected from different datasets were displayed as heatmaps constructed based on hierarchical clustering with Euclidean distance method.

2.6. Data Integration by Using N-Integration

Data integration Analysis for Biomarker discovery (DIABLO) method which is also named as N-Integration method was used for specific detection of biomarkers. Preprocessed datasets were used for increasing the accuracy of the model compared to previous approach. Supervised model was tuned based on 10-fold cross validation for 10 times. Most explanatory variables obtained from constructed components were separated and underwent correlation. Heatmaps based on hierarchical clustering with Euclidean distance method and network visualization based on correlations between variables were used to illustrate the association of different omics datasets.

2.7. Statistical Analysis

All statistical analyses were conducted using R (version 4.4.1) [35]. Data distribution was assessed using the Shapiro-Wilk test. For subgroup comparisons, the Kruskal-Wallis test and pairwise Wilcoxon tests were applied, with multiple testing correction performed using the Benjamini-Hochberg method. A false discovery rate (FDR)-adjusted p-value threshold of <0.05 was set for statistical significance. Mixomics package was used for discovering the correlation of variables by using Machine learning methodologies were implemented utilizing mixOmics (v6.10) and caret (v6.0.94) packages, facilitating advanced data integration and predictive modeling workflows.

3. Results

3.1. This Microbiome Analysis and Functional Pathway Prediction

The results of 16SrRNA microbiome analysis indicated that CRC patients that belong to CMS1 group had enriched level of Fusobacteria (p < 0.05) and Bacteroidetes (p > 0.05) and decreased level of Firmicutes (p > 0.05) and Proteobacteria (p < 0.05) at phylum level (p < 0.05) (Figure 2). Results at genus level demonstrate that the abundance of Prevotella significantly increased in CMS1 group. Moreover, Bacteroides was significantly enriched in CMS2 and CMS3 groups. The most abundant genera were Bacteroides in each CMS group in CMS1 is 48%, 65% in CMS2 and 29% in CMS3. It was followed by Faecalibacterium as the second highest abundant genera in CMS1 and CMS2, with 15% and 4% percentage respectively. Clostridium and Faecalibacterium had almost the same relative abundance in CMS3 patient group with 5.5%. The data normality was checked by Shapiro test and the non-linear distribution of variables were confirmed. So, pairwise Wilcoxon test was applied to each CMS group to observe a significant differences between each group. The results showed no significant difference between CMS subgroups, the statistical results obtained from pairwise Wilcoxon test were indicated in Table 1. Picrust2 plugin was used for functional characterization of OTU tables. Pathway information was driven from Kyoto Encyclopedia of Genes and Genomes (KEGG) database. 478 pathways (shown in Supplementary_File_1) were found by pathway prediction analysis, univariate feature selection was performed to determine significant variables (pathways) that were decreased to 26 after selection. These pathways were kept for the multi-omics integration analysis.

3.2. Transcriptome and Pathway Enrichment Analysis

Results generated from RNA-Seq analysis were demonstrated in a volcano plot. Rn7sk, Rpph1, Rn7sl1, Rn7sl2, Znf460, Snord17, Snora73b, H1-4, H4c12, H3c7 genes were discovered to be highly upregulated in tumor tissue collected from colorectal cancer patients (Figure 3A). PCA was also carried out for RNA-Seq results. Variables that are included in transcriptome data did not cluster each CMS group properly (Figure 3B). RNA-Seq analysis was followed by pathway enrichment analysis to observe which gene is upregulated or downregulated in their related pathways. 179 pathways (from KEGG database) were found to be associated with over 5000 genes that are either downregulated or upregulated indicated in Supplementary_File_2. The pathway enrichment results were not used for the integration of omics datasets, it was driven for the validation of the selected genes attained from sPLS and DIABLO multi-omic integration analysis.

3.3. Integration of Microbiome and Transcriptome Datasets Using sPLS Method

Pairwise multi-omics integration was carried out by sPLS model for dataset generated as a result of 16SrRNA and transcriptome analysis. This first approach of our integration analysis aimed to find association by including all variables to be able to observe the relation in a comprehensive way. Variables that were selected by regression were used for creating a heatmap to illustrate the correlation between each dataset.

The heatmap analysis exhibits intricate correlation between phylum level microorganisms and gene expression profiles derived from transcriptomic data (Figure 4A). Excluding the Spirochaetes and Verrucomicrobia, the heatmap directly discriminates patients into two groups. The first cluster includes Klf3, Cldn7, Plekhg3, Mall, Rna5-8n2, Snord10, Sytl2, Dnm2, Cxadr, Tff3, Serpina1, Fam120a, Rtn4, Mıb1, Oip5-as1, S100a6, and H1-5 genes, along with the phyla Actinobacteria, Cyanobacteria, Tenericutes, TM7, Bacteroidetes, Firmicutes, and Proteobacteria. The second group encompasses the genes Adar, Csnk2a1, Rpl37, Ahctf1, Odc1, Nup98, Ccdc14, Znf292, Ccnd2, Igkc, Prdx5, Gpsm2, Gpx2, Pan3, Tkt, Rpl22, Ahcy, Nme2, Rbm12, Cct6a, Rplp2, and Slc7a1, together with the bacterial phyla Synergistetes and Fusobacteria. Notably, a strong positive correlation was observed between the phylum Spirochaetes and the expression of Oip5-as1, S100a6, H1-5, and Adar genes (r = 0.71), indicating a potential regulatory relationship. Likewise, Fusobacteria exhibited a significant positive association with Klf3, Cldn7, Sytl2, Dnm2, and Tff3 genes (r = 0.70), underscoring its potential involvement in related molecular pathways. Correlations with a coefficient below 0.5 were excluded from further investigation. Negative correlations were also identified: Spirochaetes indicates a strong negative relationship with Tff3, Serpina1, and Prdx5 (r = -0.71). In addition, Fusobacteria is inversely correlated with Odc1, Nup98, Ccdc14, Ahcy, and Cct6a genes (r = -0.69). Furthermore, Synergistetes shows negative associations with Rpl22, Ahcy, and Gpsm2.

These findings highlight the complex interplay between microbial taxa and host gene expression, emphasizing key correlations that may provide insights into functional and pathological mechanisms. The PCA that was conducted on selected variables obtained from microbiome and transcriptome datasets did not manage to discriminate each CMS group (Figure 4B).

3.4. DIABLO Integration for preprocessed microbiome, transcriptome and microbial pathway datasets

Unclassified samples were extracted from microbiome, transcriptome and microbial pathway abundance datasets and non-significant variables from each data were excluded via pairwise Wilcoxon test based on CMS subgroups. This second approach with cleanup preprocessing was applied for both increasing the model accuracy and for obtaining more precise correlations by minimizing data heterogeneity. DIABLO was carried out to integrate omics data. Pre-processed data underwent sparse partial least square regression discriminant analysis (sPLS-DA) first and the components are constructed that have variables from each omics block. The majority of explanatory variables were identified using a supervised learning model, and these variables were consolidated into components. For downstream analysis, the most informative components (the first and second ones in our case) were selected. To evaluate the discriminative power of the variables comprising these components, sample plot was employed, revealing their ability to delineate and characterize group differences effectively. Notably, variables derived from the transcriptome block demonstrated a robust capacity to capture the segregation and intricate relationships among the molecular subtypes (Figure 5).

Significant features are selected as explanatory variables for the classification of CMS subgroups as a result of the sPLS-DA regression model. Rothia from microbiome data and L.lysing and CRNFORCAT pathways were selected from microbial pathway data to explain the patient status of CMS1 subgroup. The model did not identify any transcriptomic features as discriminative for the CMS1 subgroup. However, for the CMS2 subgroup, the transcriptomic analysis highlighted key genes, including Top2a, Hsp90ab1, Ncl, Ahcy, Hspd1, Psma7, Eif2s2, Cbx3, Cse1l, Lbr, and Tpx2 as critical for distinguishing this molecular subtype. Additionally, the microbial pathway analysis identified the p-Cymene degradation pathway and the glutaryl-CoA degradation pathway, while Dialister from the microbiome data showed the highest explanatory value in classifying the CMS2 subgroup.

For the CMS3 subgroup, the model pinpointed Atp8a1, Akap13, and Reg4 from the transcriptomic dataset as highly informative. Furthermore, microbial taxa such as Desulfovibrionaceae, Gamellaceae, and Paraprevotella were identified as significant contributors from the microbiome dataset. From the microbial pathway dataset, the myo-Inositol degradation pathway was selected for its relevance in explaining the CMS3 subgroup. These multi-omics features collectively provide a comprehensive framework for understanding the molecular and microbial underpinnings of the CMS subgroups (Figure 6).

Circle plot was used to observe the relation of the variables that are obtained from different omics datasets (Figure 6). Variable coordinates are determined by correlation of component 1 and component 2 obtained from features extracted by 10-fold cross validation Features that construct the component 1 and component 2 are plotted in coordinate system. Proximity of the variables at the poles demonstrates the strong biological relation.

Further investigation was made into understanding the relation between variables. Association was demonstrated by network visualization and circos plot. Glutaryl Coa degradation and p-cymene degradation pathways were positively associated with Ahcy, Eif2s2, Hsp90ab1, Psma7, Lbr, Rpl7l1, Cse1l, Cbx3, Ncl, Hspd1, Tpx2 and Top2a genes. Especially, association of these pathways (Glutaryl Coa degradation and p-cymene degradation pathways) with Cse1l and Tpx2 genes were quite remarkable with highest correlation coefficient value (r = +0.94) (Figure 8).

Circos plot was generated to observe the correlation between variables attained from each omics data, in a closer way (Figure 9). It indicates positively and negatively correlated features by connecting them with the lines. In our study, the myo-Inositol degradation pathway exhibited a significant negative association with the genes Top2a, Tpx2, Hspd1, and Cse1l. Similarly, the glutaryl-CoA degradation pathway and p-cymene degradation pathway were negatively correlated with Paraprevotella (from microbiome dataset) and the transcriptomic genes Akap13, Reg4, and Atp8a1.

Conversely, the methionine acetate pathway demonstrated a positive correlation with the Sema5a gene. Positive associations were also identified between the glutaryl-CoA degradation pathway and p-cymene degradation pathway and the Dialister genus. Additionally, Ruminococcaceae and Lachnospiraceae were positively linked to Sema5a, while Dialister showed a positive correlation with Cse1l. In contrast, Paraprevotella displayed a negative association with the genes Tpx2 and Cse1l, underscoring a complex interplay between these microbial and transcriptomic elements within the dataset. These findings provide nuanced insights into the interconnected pathways and their relationships across multiple data modalities.

4. Discussion

Huge number of variables, content and structure of different omics datasets have been an issue to deal with. Nevertheless, the number of multi-omics integration studies increases because it provides great insights to understand complex mechanisms. In our study we aimed to produce interactions of transcriptome, microbiome, and microbial pathway datasets obtained from colorectal cancer patients. Pathway enrichment analysis was done for transcriptome data and human pathways associated with up-regulated and down-regulated genes were recorded. Both microbial pathway dataset obtained from 16SrRNA analysis and pathway enrichment table results obtained from transcriptome analysis gave promising outputs in mechanisms of colorectal cancer progression. Metacyc, Biocyc, Genecard and KEGG databases were used to evaluate the functional interactions came by the multi-omics integration analysis.

4.1. sPLS for Multi-Omics Data Integration

Partial Least Squares (PLS) was selected among numerous available approaches due to its proven utility in multi-omics integration studies and its robustness for continuous development. The sparse Partial Least Square (sPLS) provided a comprehensive overview of the data, despite the known limitation of reduced explanatory power for individual components. This limitation arises from the method’s emphasis on retaining as many variables as possible to maximize covariance between datasets. Nevertheless, the results derived from the subsequent correlation analysis following sPLS regression were highly interpretable and biologically relevant.

Notably, the analysis revealed a remarkable association between the Cldn7 gene and Fusobacteria. Cldn7, a critical member of the claudin family, plays a fundamental role in the structural integrity of tight junctions. The reduced expression level of this gene has been strongly linked to the progression of cancer [36]. Fusobacteria species, frequently detected in colorectal cancer, are known to disrupt the intestinal barrier. Previous studies have demonstrated that Fusobacterium nucleatum compromises the oral mucosal epithelial barrier by inhibiting Cldn4, another claudin family member [37,38]. Our findings suggest that the observed strong correlation between Fusobacteria and Cldn7 may indicate a similar mechanism leading to tight junction disruption and subsequent barrier collapse.

In addition, Klf3 was identified as positively correlated with Fusobacteria. Klf3 is a transcription factor that plays a pivotal role in activating WNT1 and the WNT/β-catenin signaling pathway, which are key drivers of colorectal cancer progression [39]. Supporting this observation, previous research has demonstrated that Fusobacterium nucleatum promotes colorectal cancer progression through Cdk5-mediated activation of the WNT/β-catenin pathway [40].

RNA-Seq enrichment analysis further supported these associations, revealing downregulation of Wnt2, Wnt3, Notum, Lrp5, Lrp6, Csnk2a2, Gsk3b, Chd8, Crebbp, Ep300, Tbl1xr1, Invs, Plcb3 and upregulation of Wnt2b, Wnt5b, Serpinf1, Dkk1, Sfrp1, Sfrp2, Sfrp4, Rspo2, Rspo3, Fzd1, Fzd4, Frat1, Tcf7l1, Cby1, Sox17, Ctnnd2, Ccnd3, Prkaca, Skp1, Rbx1, Ror1, Prickle1, Daam2, Rhoa, Rac2, Mapk10, Ppp3ca, Ppp3cb all of which are involved in the WNT signaling pathway.

Without the integrative capability of sPLS, the identification of the interplay between Klf3 and Fusobacteria, or Cldn7 and Fusobacteria, would have been overlooked. Specifically, the association of Fusobacteria with Cldn7 suggests a potential role in tight junction destabilization, while its correlation with Klf3 points to its involvement in WNT/β-catenin pathway dysregulation in colorectal cancer. These findings underscore the value of sPLS in revealing complex biological relationships, despite its lower explanatory power for individual components. By offering a comprehensive integrative perspective, sPLS provides essential insights for pathway-driven diagnostics, which are crucial for understanding multifaceted diseases such as cancer.

4.2. Impact of the Pre-Processing Datasets

First approach aimed to include all the variables not to miss any detailed information. Despite the fact that, even potential outcomes were received by sPLS analysis, we still needed to get rid of the noise of the datasets to produce more significant interactions. Unclassified samples were extracted from the datasets and Wilcoxon test was performed to keep only significant features. This preprocessing provided the reduction of noisiness in data and performing DIABLO for this preprocessed data demonstrated increased accuracy of the model. As well as it also increased the correlation coefficients values obtained from correlation of the variables constructed the components of the model. Especially, explanatory variables obtained from transcriptome dataset provided best discrimination of CMS subgroups.

4.3. Multi-Block sPLS for Multi-Omics Data Integration for Pre-Processed Datasets

As we interpret the results from microbial pathway point, we discovered that Glutaryl Coa degradation and p-cymene degradation pathways were positively associated with Ahcy, Eif2s2, Hsp90ab1, Psma7, Lbr, Rpl7l1, Cse1l, Cbx3, Ncl, Hspd1, Tpx2 and Top2a genes (r > 0.65). Particularly, correlation of these pathways with Cse1l and Tpx2 genes were significantly high (r = 0.94). P.cymene degradation pathway demonstrates antitumor influence on cancer patients. The mechanism behind this is thought of as a reduction of inflammatory factors like IL-1 and LEP, whereas upregulation of IL-6 that promotes growth of Bifidobacteria and Clostridium IV [41]. Another prominently identified microbial pathway was the Glutaryl-CoA degradation pathway, which plays a pivotal role in the metabolism of lysine, tryptophan, and hydroxylysine, and has been recognized for its anti-tumor properties. Prior research demonstrated that inhibition of this pathway induces NRF2 glutarylation, thereby enhancing its stability. This process subsequently activates the ATF4-ATF3 signaling cascade, ultimately leading to cancer cell apoptosis [42]. In the present study, the observed increase in p-cymene levels within tumor tissue suggests that this compound may represent a host immune response mechanism, potentially mediated by commensal microbiota to counteract tumor development and progression. Evidence of p.cymene mediated inhibition of MMP-9 (matrix metalloproteinases) expression shows its anti-tumor activity. Besides that, MMP-2 and MMP-9 inhibition play an important role in degradation of extracellular matrix. Moreover, ERK1/2 and p38 MAPK signal pathways were also inhibited by p-cymene [43]. Direct or the indirect impact of the p-cymene on the inhibition of these genes emerges promising therapeutic strategy. In our study we have found that MMP-2 genes were downregulated in diabetic cardiopathy I, preteoglycyns in cancer, estrojen signaling, fluid shear stress and atherosclerosis, relaxing signaling, bladder cancer, and GnRH signaling pathways (Table S1). So that, our results that are strengthened by literature strongly propose the potential interaction of p.cymene on MMP-2 inhibition.

Hspd1 gene is a heat shock protein which is related to transcription and peroxisomal lipid metabolism. The Hspd1 gene has been identified as essential for maintaining metabolic fitness, with its activity being critical for sustaining cellular energy homeostasis. A study published in 2021 demonstrated that loss-of-function mutations in Hspd1 lead to a profound energetic collapse, which directly impairs the ability of cancer cells to proliferate and expand, thereby disrupting tumor progression. [44]. Cbx3 is taking role in transcription, RNA polymerase and promoter opening pathways. Also, studies indicate that this gene is overexpressed in a variety of cancers but still a poor prognosis biomarker in cancer patients [45,46]. Cse1l gene is chromosome segregation 1 like gene that play role in p53 effector and Golgi to endoplasmic reticulum traffic pathways. Researchers discovered that Cse1l directly promotes the development of tumors in various malignancies [47]. Lbr gene is Lamin B receptor that takes role in super pathway of cholesterol biosynthesis and transcription pathways. Its role in cancer was revealed as Lbr protects the genome from tumorigenesis and chromosomal instability [48]. Tpx2 gene is a microtubule nucleation factor found in microtubule organization pathways. It depicts upregulated expression level profiles in CRC patients [49]. In our study, significantly discovered microbial pathways were having a notable protective role against cancer, whereas significantly found genes were generally oncogenes. This indicates the impact of the microbial community on tumor tissue to suppress tumor growth. Further research is needed to unravel the exact mechanism of action in these microbial pathways against cancer and their relation through these oncogenes.

Circos plot was generated to observe the correlation between variables in a closer way. It indicates positive and negatively correlated features by connecting them with the lines with the correlation cut off = 0.7. In our study, MYO-inostol degradation pathway was found to be negatively associated with Top2a, Tpx2, Hspd1 and Csel1 genes. In addition to that, The Glutaryl-CoA degradation and p-cymene degradation microbial pathways demonstrated a negative association with Paraprevotella in the microbiome dataset and with the Akap13, Reg4, and Atp8a1 genes in the transcriptome dataset. Conversely, the meth-acetate pathway exhibited a positive correlation with the Sema5a gene. Additionally, Glutaryl-CoA degradation and p-cymene degradation pathways were positively linked to the Dialister genus.

Notably, both Ruminococcaceae and Lachnospiraceae showed a positive correlation with the Sema5a gene, while Dialister was positively associated with the Csel1 gene. In contrast, Paraprevotella exhibited a negative association with Tpx2 and Csel1 genes, suggesting potential functional interactions between microbial communities and host gene expression.

Components of the last integration model with implementation of preprocessing datasets was constructed by Hspd1, Cbx3, Csel1, Lbr and Tpx2genes from transcriptome data, Actinobacteria, Spirochaetes, Fusobacteria, Bacteroidetes and Proteobacteria from microbiome data and cholorthocleavage pathway, crnforcat pathway, MYOinotsitol degradation pathway, p.cymene degradation pathway and glutaryl.CoA degradation pathway from microbial pathway abundance dataset. Contribution of the datasets to explain CMS groups were increased and weights for transcriptome, microbiome and pathway are reached to 0.80, 0.73 and 0.60 respectively. Model that was designed to predict CMS groups achieve over 95%. Correlation analysis made for the explanatory components and Hspd1, Cbx3, Csel1, Lbr and Tpx2 genes were discovered to be highly associated with p.cymene degradation and glutaryl-CoA degradation pathways.

Hspd1, Cbx3, Csel1, Ncl, Top2a and Tpx2 genes and Proteobacteria indicated negative association. Moreover, Proteobacteria was also found to be negatively correlated with p.cymene degradation and glutaryl CoA degradation pathways. P.cymene promotes the growth of genus Bifidobacterium and genus Clostridium which they belong to Actinomycetota and Firmicutes respectively[41]. Including this negative association, all output shows reasonable significant outcomes which are potential candidates raking part in colorectal cancer progression. However experimental validation is still necessary to comprehend the exact mechanism of action at molecular level.

5. Conclusions

Cancer is one of the most studied diseases worldwide, however there is still no exact treatment approach to get over it. Rather than discovering new biomarkers, it is quite significant to understand a person’s specific outputs to unravel the mystery behind complex biological mechanisms in cancer like complex diseases. Strong interactions of Klf3 and Cldn genes with Fusobacteria were quite promising output to protect the tight junction structure on intestinal barrier. Also, the role of p-cymene degradation and Glutaryl-Coa degradation pathways was discovered as very important defensive mechanisms against tumor. All these findings are quite promising to be validated in further experimental research. Increasing number of multi-omics integration approaches in research will provide to evaluate the details of the molecular mechanisms and eventually we will be able to serve person specific, pathway dependent solutions to cure the cancer.