Analyzing aberrant DNA methylation in colorectal cancer uncovered intangible heterogeneity of gene effects in the survival time of patients

Colorectal cancer (CRC) involves epigenetic alterations. Irregular gene-methylation alteration causes and advances CRC tumor growth. Detecting differentially methylated genes (DMGs) in CRC and patient survival time paves the way to early cancer detection and prognosis. However, CRC data including survival times are heterogeneous. Almost all studies tend to ignore the heterogeneity of DMG effects on survival. To this end, we utilized a sparse estimation method in the finite mixture of accelerated failure time (AFT) regression models to capture such heterogeneity. We analyzed a dataset of CRC and normal colon tissues and identified 3406 DMGs. Analysis of overlapped DMGs with several Gene Expression Omnibus datasets led to 917 hypo- and 654 hyper-methylated DMGs. CRC pathways were revealed via gene ontology enrichment. Hub genes were selected based on Protein–Protein-Interaction network including SEMA7A, GATA4, LHX2, SOST, and CTLA4, regulating the Wnt signaling pathway. The relationship between identified DMGs/hub genes and patient survival time uncovered a two-component mixture of AFT regression model. The genes NMNAT2, ZFP42, NPAS2, MYLK3, NUDT13, KIRREL3, and FKBP6 and hub genes SOST, NFATC1, and TLE4 were associated with survival time in the most aggressive form of the disease that can serve as potential diagnostic targets for early CRC detection.

precision in each of the above steps leads to misleading conclusions.Among others, two issues affect precision: removing genomic positions with missing values or low read-depth and ignoring the heterogeneity of DMG effects on survival times.
To accurately predict the differentially methylated profiles in CRC, one must consider all biological and environmental factors such as dietary 11 , aging 12 , and hazardous behaviors 13 (e.g., smoking), among others.Such factors are often ignored by most studies when predicting methylation profiles.In addition, methylation data always suffer from heavy missing values that can affect subsequent analyses.For instance, 68% of CpG sites have missing values in at least one sample in our dataset (Section 2).Almost all DNA methylation pipelines, except a few such as the DMCHMM method 14 , filter out such positions from the analysis.We used DMCHMM to not only account for extra covariates but also efficiently impute the missing values.
Having identified the differentially methylated genes (DMG) associated with CRC and validating them, it is crucial to identify their underlying signaling pathways that regulate gene expression 15,16 .The main known CRC pathways are Wnt 17 , MAPK 18 , TGF-β 19 , and TP53 20 .Although significant progress has been made in understanding the biology of CRC, there are still many unknown pathways and mechanisms involved in this disease.Identification of hub genes, also known as driver genes is the next step in the analysis of biomarker detection.Hub genes play a critical role in regulating several genes in the biological network and have the potential to be regarded as therapeutic targets in CRC 21 .
In the next step, the relationship between identified DMGs and the survival time of CRC patients should be evaluated.Most studies employ a limited panel of biomarkers selected through conventional univariate Cox proportional hazard regression models and overlook the potential effects of the rest of the biomarkers [22][23][24] .In a recent study 25 , the Cox-LASSO survival model was used to account for a larger set of biomarkers but ignored the heterogeneity of covariate effects.To the best of our knowledge, none of the studies have taken into account the heterogeneity of DMG effects on survival time.To address this problem, one may use the sparse estimation method in the finite mixture of accelerated failure time (AFT) regression models 26 .Prior to this step, it is common to screen the number of genes to a manageable magnitude.This process can be done by selecting the top highly correlated genes with survival time of the patients using the correlation-adjusted scoring method 27 .
This study aimed to identify CRC-related DMGs to serve as potential biomarkers for early detection by including all the available information in the data and avoiding the exclusion of any genomic position.To this end, we acquired a high-throughput DNA methylation dataset which consists of patients with CRC and healthy individuals.Information on age, history of smoking, and drug abuse was also collected.A description of the data is provided in Section "Methods".Information on other datasets used for validation and survival analysis and all statistical and Bioinformatics methods are listed in this section.In Section "Results", a comprehensive analysis of data is conducted.Section "Discussion" gives a discussion and some concluding remarks.

Methods
In this section, we outline the data analysis process we followed to detect DMGs, hub genes, and their effects on the survival time and enriched pathways of CRC. Figure 1 depicts the flowchart of this process.

Phase I (pre-processing of discovery samples)
To identify methylation-based CRC biomarkers, information on 6 patients with adenocarcinoma of CRC and 6 normal males was obtained.Two groups were matched based on age, and family history of cancer 28 .
The methylation profiles in our dataset are derived from a three-step pre-processing phase conducted through SureSelectXT Human Methyl-Seq.Initially, the purity and quantity of 12 DNA tissue samples were assessed using specific criteria, including a minimum concentration of 50ng/µl , a purity ratio (A260/A280) ≥ 1.7 , a volume of at least 20ng/µl , and a total amount exceeding 3.0µg .Subsequently, global methylation profiles of CRC and normal samples were analyzed using SureSelectXT Human Methyl-Seq.www.nature.com/scientificreports/Pre-sequencing tasks, such as sample collection and DNA extraction, were consistently carried out by a single technician.Experimental conditions for all samples remained constant both before and after sequencing.During the sequencing process, all sample runs were executed simultaneously using the same device, employing Nextgeneration sequencing technology, a highly parallel sequencing method.This approach minimized the potential introduction of batch effects attributable to non-biological factors such as variations in laboratory conditions, personnel, and equipment used in the experiment.
In the second step, a quality control assessment of total reads using FastQ 29 was conducted.This step aimed to provide informative global and graphical representations of read quality in methylation sequencing, both pre and post-alignment.Notably, our data consistently exhibited high quality in raw sequencing reads across all samples.Subsequently, Trim Galore 30 was utilized to process the raw sequencing reads.This involved the removal of sequencing adapters, specifically the Illumina universal adapter, and discarding the low-quality bases (those with quality scores below 67, as per Illumina standards) located at the 3ʹ end of reads.Additionally, any ambiguous bases found in both reads were removed.
Finally, the raw bisulfite sequencing data were aligned to the human reference genome (GRCh37/19) using Bismark 31 .Several comparisons and visualization confirmed minimal to no presence of batch effects in our data.This discovery dataset includes methylation read counts and read-depth for each CpG site, generating 57 to 76 million Illumina sequencing reads per subject.Between 88.5% and 89.8% of sequenced reads were mapped to either strand of the human genome (GRCh37/19).On average, each CpG site was sequenced between 19 × and 24 × per sample.The sequencing details for the subjects are presented in Table 1.Approximately 68% of the 19,530,818 CpG sites have missing information in at least one sample.

Phase II (identification of differentially methylated genes)
We utilized the DMCHMM pipeline 32 to identify CpGs with differentially methylated patterns between CRC and normal discovery samples.We specifically did not remove any position with missing information or low readdepth.The missing information was imputed using DMCHMM via hidden Markov models 14 .Significant differentially methylated cytosines (DMCs) were selected based on the FDR threshold of 0.05.DMCs were aligned to the human reference genome (GRCh37/19) using the UCSC Genome Browser (https:// genome.ucsc.edu).A gene whose promoter was mainly hypo-or hyper-methylated was classified as hypo-or hyper DMG, respectively.

Phase III (cross-platform validation)
To validate our result, several methylation profiles (GSE53051 33 , GSE77718 34 , GSE101764 13 , GSE42752 35 , GSE48684 36 ) were extracted from the Gene Expression Omnibus (GEO, https:// www.ncbi.nlm.nih.gov/ geo/).Of these datasets, a total of 212 CRC and 242 normal mucosa tissue samples were selected based on setup conditions to minimize the confounding effect of other variables.These datasets have provided valuable insights into the molecular alterations that occur in CRC, and their findings have implications for the diagnosis and treatment of this disease.For the analysis of methyl array profiles of validation sets, the GEO2R (http:// www.ncbi.nlm.nih.gov/ geo/ geo2r/) web tool and the limma R-package 37 were used.To mitigate batch effects, we applied the 'removeBatchEffect' option from the package.A probe was considered differentially methylated if its adjusted p-value was less than 0.05, and the absolute of log 2 of methylation fold change was greater or equal to 1.The differentially methylated probes were aligned to the human reference genome (GRCh37/19) using the FDb.InfiniumMethylation.hg19 package 38 .In the last step, we compared the lists of DMGs based on the validation sets and our discovery samples to identify consistent hypo/hyper-methylated genes across different populations and platforms.

Phase IV (network construction and functional analysis)
In order to investigate the Protein-Protein Interaction (PPI) network and module analysis, we utilized the 'Search Tool for the Retrieval of Interacting Genes' (STRING) database.We set the interaction score threshold to 0.4 to screen for high-confidence interactions and visualized the resulting network using the Cytoscape www.nature.com/scientificreports/(Version 3.9.1).Next, we employed the Molecular Complex Detection (MCODE) algorithm to uncover densely connected substructures within the network.The MCODE score must be greater than 3 and the minimum number of nodes must be 4.In order to identify key hub genes within the network, we used the cytoHubba plugin and considered the degree of centrality as a parameter.
To gain insight into the biological mechanisms that are driving CRC and prioritize identified DMGs, we performed functional and pathway enrichment analysis using DAVID 40 (https:// david.ncifc rf.gov/).Gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG 41 ) pathways were considered significantly enriched if the p-values were less than 0.05 and the q-values were less than 0.1.The visualization of the identified GO terms and KEGG pathways were done with the clusterProfiler 42 , pathfindR 43 , and ShinyGO 44,45 (http:// bioin forma tics.sdsta te.edu/ go/) packages.

Phase V (Uncovering intangible heterogeneity of DMG effects on survival time)
To explore the relationship between identified DMGs and survival time, the DNA methylation profiles of 521 samples were obtained from The Cancer Genome Atlas (TCGA) network 46 .Complete information on clinical variables including days to follow-up and the status of the patient were analyzed.
We conducted several preliminary analyses on the overall survival time of patients with CRC.First, we estimated the density of the logarithm of the survival times using the Kaplan-Meier estimator.The density plot in Figure 2 shows a mixture distribution.Second, we applied mixture and non-mixture models of normal distributions using the mixtools package 47 .The BICs for the mixture of components K = 1, 2, 3, 4 , and 5 were estimated as 777.73, 709.02, 712.31, 722.51, and 721.68, respectively, with the lowest BIC observed for the mixture with 2 components.Finally, we employed mixture and non-mixture models of semiparametric scaled data using the stochastic EM algorithm 48 via mixtools.Please note that the package only outputs a mixture with two components, and the BICs for the non-mixture K = 1 and the mixture of 2 components were 1323.62 and 582.01, respectively.Similarly, the lowest BIC was observed for the mixture with 2 components.All of these preliminary analyses indicate heterogeneity in the overall survival time of patients with CRC.
Furthermore, a comprehensive literature review revealed numerous conflicting results.For example, researchers 49 found a significant association between the methylation level of RASSF1A and the overall survival of CRC patients, while other studies 50 did not observe such an association.
Moreover, several studies indicate stage-specific 51 and age-specific 52 effects of DNA methylation in certain genes on the survival outcomes of CRC patients.These results clearly suggest heterogeneity in the overall survival time of patients with CRC.
We hypothesize that the effect of identified DMGs and hub genes on the overall survival time varies in each subpopulation, but not all DMGs and hub genes have an effect in each subpopulation, implying that the underlying regression model is sparse.To capture such heterogeneity, we employed the sparse estimation method in the finite mixture of AFT regression models 26 .The details of the method are given in Section "Sparse finite mixture of AFT regression models to estimate the DMG effects on the survival times" below.The response variable is "Overall Survival Time", and the independent variables are the log-transformed average methylation of identified DMGs or hub genes that we discovered through Phase II-IV.The goal of such a regression model is to estimate the effects of each gene in different sub-populations of the response variable, providing insights into the effects of each gene on the survival time of patients with CRC.It is important to note that the response variable (survival time) is subject to right-censoring.The sparse estimation method requires tuning parameters, which are www.nature.com/scientificreports/estimated using a data-adaptive approach explained in Section "Sparse finite mixture of AFT regression models to estimate the DMG effects on the survival times".

Sparse finite mixture of AFT regression models to estimate the DMG effects on the survival times
As hypothesized above, the overall survival time of patients with CRC is heterogeneous; thus, we hypothesize that the relationship between overall survival time and DMGs and hub genes found in Phase II-IV is heterogeneous.Such heterogeneity cannot be detected using a regular AFT regression model for censored data.Therefore, we employ the finite mixture of the AFT regression model to capture intangible DMG and hub gene effects on survival time.To this end, we use the finite mixture of AFT regression model: where f k and S k are respectively the density of normal distribution and its survival function, y = log(t) , t is the overall survival time, δ is an indicator representing right-censored (i.e., δ = 0 if time is censored and 1 if it is not censored), X is the vector of all DMGs and hub genes discovered in Phase II-IV, β k is the vector of effects of these genes in Component k of the mixture model, σ 2 k is the variance, and π k is the proportion of the kth component.It is common to screen the number of genes prior to analysis in case of a large number of identified genes.To this end, we applied a correlation-adjusted score method using the carSurv package 27 to screen the genes.
Next, we used the fmrs package 53 to fit finite mixture and non-mixture of AFT regression models to the data.We employed the smoothly clipped absolute deviations (SCAD) penalty 26 .This sparse method requires K tuning parameters which are estimated via the data-adaptive component-wise BIC method proposed in Shokoohi et al. 26 .

Differentially methylated cytosine detection
We identified 2,691,019 DMCs between CRC and normal groups of the discovery dataset while adjusting for the potential confounding effect of smoking history or drug abuse.Of these identified DMCs, 1,985,557 positions were hypo-methylated and 705,462 CpGs were hyper-methylated in CRC vs normal samples.The heatmaps (see R package pheatmap 54 ) in Fig. 3a indicate a clear clustering pattern between the CRC and normal samples based on the predicted methylation levels of DMCs.
To explore the genomic location of the DMCs, we analyzed their distribution across different regions and summarized the results in Fig. 3b.Intergenic regions were found to harbor the majority of the detected DMCs both in the hypo and hyper categories.Notably, we observed that 32% of hyper-methylated DMCs were located in CpG islands, while only 9% of hypo-methylated DMCs were found in these regions.Additionally, the regions with the highest percentage of hyper-methylated DMCs were identified in introns, exons, and CGI shores.The Chord diagrams (see R package circlize 55 ) in Fig. 3c gives a comprehensive overview of how hyper and hypo-methylated DMCs were distributed across different genomic regions.Our findings suggest that many DMCs in intergenic regions were expanded to intronic regions in both hypo and hyper-methylated categories.
Given the potential significance of promoter methylation in cancer development and progression, we focused our subsequent analysis on DMCs located on gene promoters, which encompassed 268,978 CpGs.These CpGs resided on 3406 gene promoters, of which 1394 were hyper-methylated and 2012 were hypo-methylated.The list of DMGs is available as supplementary material.

Robust DMGs in CRC
To verify the robustness of identified DMGs, we performed a cross-platform procedure with DMGs identified in selected GEO datasets as depicted in Fig. 4a (see R package venn 56 ).The comparison revealed a total of 1571 overlapped DMGs that were consistently identified across multiple studies.As Fig. 4b (see R package karyoploteR 57 ) illustrated, the identified DMGs were spread almost evenly across different chromosomes, with chromosomes 1 and 7 having some dense regions of CRC-related DMGs.Within this set, 917 genes were hypo-methylated, and 654 genes were hyper-methylated.We focused our subsequent analysis on these identified DMGs to gain a deeper understanding of their role in CRC pathogenesis.
GO enrichment KEGG pathway analysis.The analysis of robust DMGs in CRC utilizing the DAVID tool yielded a variety of enriched biological processes, molecular functions, and cellular components.Specifically, the hyper-methylated DMGs were found to be principally involved in 'cell fate commitment' , 'regionalization' , 'embryonic organ morphogenesis' , 'embryonic organ development' , 'pattern specification process' , 'animal organ morphogenesis' , 'tube morphogenesis' , 'tube development' , and 'neurogenesis' in the context of biological processes (Fig. 5a, see R Shiny package ShinyGO 44 ).Enriched cellular components included 'basement membrane' , 'integral component of postsynaptic membrane' , and 'Collagen-containing extracellular matrix' (Fig. 5c).Additionally, KEGG pathway analysis indicated that hyper-methylated DMGs were significantly enriched in several pathways, including 'signaling pathways regulating pluripotency of stem cells' , 'axon guidance' , 'morphine addiction' , 'rap1 signaling pathway' , 'circadian entrainment' , and 'pathways in cancer' (Fig. 5e and Table 2).Regarding biological processes, the hypo-methylated DMGs were found to be associated with a number of processes including 'keratinization' , 'keratinocyte differentiation' , 'epidermal cell differentiation' , and 'epithelial cell differentiation' (Fig. 5b).Furthermore, analysis of the cellular component pathway revealed that the hypomethylated DMGs were most significantly enriched in the 'cornified envelope' , 'integral component of the synaptic membrane' , and 'integral component of the postsynaptic membrane' .Notably, these cellular components demonstrated the highest FDR and fold enrichment (Fig. 5d).Regarding molecular functions, the pathways  signaling in cardiomyocytes' , 'MAPK signaling pathway' , 'arrhythmogenic right ventricular cardiomyopathy' , and 'cell adhesion molecules' (Fig. 5f).These results offer valuable insights into the potential mechanisms of DMGs in CRC and identify possible therapeutic targets for this disease.A comprehensive summary of the KEGG pathways of hyper-methylated DMGs can be found in Table 2.
PPI network construction.We ran a PPI network to further investigate the complex interactions between DMGs and find important hub proteins.A total of 606 PPI nodes of the hyper-methylated DMGs were   constructed on the basis of the STRING database (Fig. 6, see R Shiny package ShinyGo 44 ).The 16 node proteins, including KIT, SEMA7A, BDNF, MEF2A, LDB2, GATA4, LHX2, SOST, CTLA4, NKX2-2, TLE4, BMP5, NFATC1, ZFPM1, DPYSL2, and ITGA2B that showed a close interaction with other node proteins were chosen as hub genes (Fig. 7a, see Cytoscape 39 ).The most important biological process and KEGG pathways of hub genes are shown in Fig. 7b and c.One important module was selected when the number of nodes is greater than 4. The key module demonstrated functions enriched in pathways such as Wnt signaling 58 (Table 2 and Fig. 8, see R Shiny package ShinyGo 44 )).
We performed a survival analysis using the TCGA-selected samples to investigate the association of selected hub genes with the survival time of CRC patients.Based on Fig. 9a-d (see GEPIA2021 59 ), those patients with gene SEMA7A ( p = 0.024 ), SOST ( p = 0.027 ), NFATC1 ( p = 0.017 ), and TLE4 ( p = 0.0061 ) being upregulated, had a significantly lower probability of survival.However, this conclusion is based on univariate analysis, and the effect of other genes and the potential heterogeneity of DMG effects were ignored.We reanalyzed these data by accounting for the heterogeneity of DMG effects and obtained different results as follows.
Intangible heterogeneity of DMG effects on survival time.We studied the relationship between the average promoter methylation of the identified DMGs and the survival time subject to right-censoring by accounting for the heterogeneity of gene effects using an independent set of 521 TCGA CRC samples.To this end, we screened all the 1571 candidate DMGs using the correlation-adjusted regression survival scores to obtain the list of top candidate covariates.This process led to the selection of 95 highly correlated DMGs.These genes were also dysregulated in the TCGA samples.In addition, 4 hub genes that were related to the survival time of CRC patients were added to the list of covariates.
Our analysis yielded a two-component mixture of AFT regression model.The estimated gene effects on the survival time are given in Table 3.The result showed that 46% of the subjects were classified into Component 1, which is the most aggressive form of the disease.Figure 10 (see R package fmrs 60 ) depicts the posterior probability of a subject belonging to Component 1. From this figure, we noticed that all living patients were classified into Component 2, which is the less aggressive form of the disease.A total of 83 and 18 DMGs were active in Components 1 and 2, respectively.Twelve genes including HLA-F, MMP2, MT1A, RFPL4B, SIX6, ZFAT, BCKDK, AMOTL1, ADCY10, KCNK10, STAU2, and NOC4L were not related to survival time in either of the components.These findings demonstrate the heterogeneity of DMG effects in CRC data and justify using a sparse mixture modeling rather than a univariate one.In addition, the DMGs with active promoters in Component 1 can be considered biomarkers for CRC prognosis.The bioinformatics and biological information of selected DMGs are given in Table 4.

Discussion
Colorectal cancer is one of the deadliest cancers in the world.Given that early stages of CRC do not display symptoms, proactive screening is the only viable approach to identify the disease 83 .As DNA methylation changes are closely associated with cancer, their role in CRC biomarker detection in the early stages of cancer is of great importance.Although many CRC biomarkers have been detected in the literature, only a few are used in practice.
Our findings resulted in identifying new biomarkers for CRC which can be used for diagnosis and prognosis.
We identified 1,571 DMGs most of which have been previously studied in the literature.Among them, SEPT9, NDRG4, VIM, APC, SFRP1, SFRP4, and SFRP5, 84 are the most important CRC-related ones.We also explored CRC-related hub genes.Fourteen functional modules that may play important roles in the early detection of CRC were highlighted and the sub-network of hub genes KIT, SEMA7A, BDNF, MEF2A, LDB2, GATA4, LHX2, SOST, CTLA4, NKX2-2, TLE4, BMP5, NFATC1, ZFPM1, DPYSL2, and ITGA2B was extracted.These hub genes were flagged as potential diagnostic and therapeutic targets for CRC in our analysis.
In addition to the diagnostic role of our identified hub genes such as NKX2-2, KIT, BNDF, and TLE4 in CRC and its sub-types 74,[85][86][87] , their roles in increasing CRC risk, tumor progression, and targeted therapy have been investigated.For instance, MEF2A 88 and BMP5 89 increase the CRC risk.Up-regulation of the expression of ITGB7 and ITGA2B has been found to be significantly associated with death by sodium butyrate-induced CRC organoids 90 .Moreover, some studies 91,92 have shown effective treatments by targeting CLT-4 and LDB2n.
There is a rich literature on the contribution of some of our identified hub genes in CRC and less evidence in support of some others such as LHX2, ZFPM1, and DPYSL2.For instance, the differences in tumor and corresponding adjacent benign tissues regarding LHX gene expressions have been investigated 93 .However, contrary to our findings, they did not find any statistical differences for LHX2 and LHX3 genes.Furthermore, the upregulation of ZFPM1 was revealed in molecular high-risk patients with cytogenetically normal acute myeloid leukemia 94 , yet its diagnostic value in CRC has not fully been confirmed 95 .SEMA7A is also one of our selected hub genes that play a key role in several cancers including pancreatic, breast, and lung cancers 69,[96][97][98] .However, there has been less attention on the role of SEMA7A in CRC.Further investigation is required on our flagged DMGs.
Although there are many mechanisms that drive CRC, only a handful of them have been discovered in past studies.As researchers continue to genotype large panels of CRC tumors, it can be expected that additional new pathways of CRC carcinogenesis will be revealed.SOST, an identified hub gene in our study, plays a vital role in inhibiting the Wnt signaling pathway by binding to the Wnt co-receptor, LRP5/6, and preventing its activation 99 .Therefore, decreased SOST expression could lead to an increase in Wnt signaling, promoting CRC cell proliferation, migration, and survival.Another identified hub gene is TLE4 which is involved in the negative regulation of the canonical Wnt signaling pathway.Only a few investigations provided evidence of TLE4 upregulation in CRC biopsies, partially through regulation of the JNK/c-Jun pathway 73 .Moreover, recent studies that focus on the NFAT signaling pathway showed a promising strategy for CRC treatment 64 .
Heterogeneity is one of the key features of genomic data.Specifically, there is evidence of the heterogeneity of DMG effects on the survival of CRC patients in the literature and in our dataset.The finite mixture of the AFT regression model is a plausible method to uncover such intangible heterogeneity.Our analysis suggested a mixture of two-component mixture of the AFT regression model in which patients were separated into two subgroups based on their vital status.In this model, almost all of the deceased patients were classified into the  www.nature.com/scientificreports/most aggressive form of the disease (Component 1).In Component 1, 83 DMGs including NMNAT2, ZFP42, NPAS2, MYLK3, NUDT13, KIRREL3, and FKBP6 had an effect on the survival time of the patients.The relation between some of these DMGs and survival time has been previously reported 100 .On the other hand, there are a few discoveries regarding other genes.For instance, significantly higher expression of NMNAT2 in CRC tissues compared to normal ones have been found, yet this gene was not a prognostic factor for overall survival 101 .Note that, while the hub genes SOST, NFATC1, and TLE4 were associated with survival in the univariate Cox model, they were only associated with survival time in the most aggressive form of the disease in our study.
Our study does not exclusively depend on bioinformatics analysis, as we have employed several statistical and machine learning analyses.These include modeling methylation profiles, identifying DMCs via DMCHMM, conducting statistical tests, performing multiple validation analyses, and applying statistical learning algorithms to survival times via fmrs.One of the advantages of the DMCHMM method is that it does not require a large number of samples or matched samples, as it is highly flexible and can accommodate various experimental designs.It demonstrates significant power, particularly when dealing with moderate to low sample sizes.

Figure 1 .
Figure 1.Study workflow for the analysis of CRC datasets.

Figure 3 .
Figure 3. Genomic location of identified differentially methylated CpGs and their predicted levels in CRC (T) and normal (N) samples using DMCHMM.The hierarchical clustering of CRC and normal samples in the heatmaps is based on complete linkage.

Figure 4 .
Figure 4. Summary of common identified DMG and their distribution.

39 softwareTable 1 .
Summary statistics of methylation sequencing reads of discovery samples.

Table 3 .
Estimated DMG effects in the two-component mixture of accelerated failure time regression model in the CRC data.Function Transcriptional corepressor that binds to various transcription factors.Inhibits transcriptional activation by PAX5, CTNNB1, and TCF family members in the Wnt signaling pathway.Cancer Overexpression may play a role in CRC development and progression, partly through the JNK/c-Jun pathway.It is a candidate for risk stratification of cancer recurrence after curative resection of early-stage CRC.Pathway Activation of anterior HOX genes in hindbrain development during early embryogenesis, Activation of HOX genes during differentiation Function Mutations linked to Charcot-Marie-Tooth disease type 1D, Charcot-Marie-Tooth disease type 4E, and Dejerine-Sottas syndrome.Cancer Targeting EGR2 may provide a therapeutic strategy to eliminate colon cancer stem cells and block nervous system-driven disease progression through differentiation.Function Crucial for trophoblast and placental development, possibly contributing to cytokinesis by cleavage furrow PIP2 hydrolysis.Controls neurite outgrowth by suppressing RhoA/Rho kinase signaling.and stability, maintaining neuronal polarity.Binds axonal microtubules and neural plasma membrane components, acting as a link between them.Its localization within the cell body helps define axonal polarity.Cancer Hyper-methylation in MAPT is associated with poor prognosis in stage II CRC patients.Function Among human counterparts of Drosophila melanogaster forkhead transcription factor.Expresses in the lung and placenta, activating transcription of several lung-specific genes.Cancer Regulates PRUNE2 transcription in CRC pathogenesis and is hyper-methylated in CRC samples.Reference Li et al. 80 , Hauptman et al. 81 Pathway Cytokine signaling in immune system, TNFR2 non-canonical NF-kB pathway Function Cytokine in the TNF ligand family, acting as a bidirectional signal transducer with TNFRSF9/4-1BB.Key role in antigen presentation, cytotoxic T cell generation, and T lymphocyte activation and proliferation.

Table 4 .
Bioinformatics and biological information of selected DMGs related to survival in colorectal cancer patients.