Blood-based metabolic markers for diagnosing colorectal cancer with supporting multi-omics evidence

A total of 715 serum samples were collected from the Zhuhai People's Hospital (the First Affiliated Hospital of Macau University of Science and Technology) between 2020 and 2023, including 248 patients diagnosed with CRC and 467 noncancer controls (NCCs) (Table 1). CRC diagnoses were confirmed histopathologically. Control participants underwent routine colonoscopy and were confirmed to be free of malignant or precancerous lesions. All participants provided informed consent, and the study protocol was approved by the institutional ethics committee.

Blood samples were collected by venipuncture, and the patient was required to fast for at least 8 hours but not more than 16 hours, preferably 12–14 hours. The serum was separated within 2h and centrifuged at 3000 rpm for 10min at room temperature. The supernatant was transferred to a centrifuge tube and centrifuged again at 14,000 rpm for 10min at 4 °C. Serum was obtained from the supernatant, and the serum samples were frozen and stored at −80 °C until sample processing.

All the serum samples stored at −80 °C were first thawed on ice before preparation. After the samples were vortexed for 30 s, a 10 μL aliquot of serum from each sample (including all patients with tumors and healthy controls) was mixed thoroughly for quality control (QC). Both the serum and QC samples were extracted at a volume ratio of 100 μL, and 400 μL of MeOH was added to each sample to initiate protein precipitation and metabolite extraction. After being vortexed for 30 s, the mixture was centrifuged at 14,000 rpm for 10min at 4°C. Two hundred microliters of the supernatant were transferred to new Eppendorf (EP) tubes, and the samples were subsequently dried on a speed vac for 150min at 37°C and stored at −80°C. Before UPLC-MS analysis, the dried samples were redissolved in 50 μL of ultrapure water. The samples were vortexed for 30 s and subsequently sonicated in a water bath for 30 s, followed by centrifugation at 14,000 rpm for 10min at 4°C. Finally, 20 μL of the supernatant was collected and analyzed immediately. The pooled QC sample was injected five times at the beginning of the analysis to ensure system equilibrium, after which it was injected every ten samples during serum sample detection to further monitor system stability.

Untargeted LC-MS profile analysis of polar metabolites was performed on a UPLC system (ACQUITY UPLC I-Class system, Waters Corp.) coupled with tandem ESI–QTOF mass spectrometry (Synapt G2-Si, Waters Corp.). A 2 μL sample was injected into the chromatograph and separated on an ACQUITY UPLC HSS T3 1.8 μm, 2.1× 100mm i.d. column (Waters Corp.). The column temperature was controlled at 30°C. Mobile phase A was H₂O containing 0.1% formic acid, and mobile phase B was 0.1% formic acid in CAN. During the entire analysis, the autosampler temperature was maintained at 4°C to avoid sample degradation. Two injections were performed for each sample to collect positive and negative data in full scan mode with a mass range of 50 to m/z at a resolution of 10,000. The electrospray ionization (EI) capillary voltages and cone voltages were set as 2.0 kV and 20V, respectively. The source temperature was 100°C, and the desolvation temperature was 200°C. The desolvation gas flow rate was 500 L/h.

The raw MS data were converted to mzXML format using MSConvert in the ProteoWizard software package (v3.0.23089) (15) and processed using the R-based XCMS package (16) for peak extraction, peak optimization, retention time alignment, feature formation, feature grouping, aggregation, spectrum extraction, spectrum aggregation, compound identification, and quantitative data generation. The following parameters were used: peak width = c (5, 20), noise = 1000, snthresh = 3, ppm = 20, binSize = 6, minFraction = 0.4, bw =20. The compound annotations of the metabolites were matched with the Human Metabolome Database (HMDB) (17) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (18) database by metID (19). The parameters used were ms1.matchh.ppm = 15, rt. match.tol = 30, threads = 30, column = rp (20). The ‘metid’ package was used for metabolite identification based on public databases.Candidates with the highest spectral similarity to fragmentation patterns from databases (HMDB and KEGG) were prioritezed. Candidates present in two databases were given higher priority, reducing reliance on single-database matches. Annotations consistent with expected adduct patterns (e.g., [M+H]+, [M+Na]+) were retained, while ambiguous matches were flagged for manual review. Redundant annotated metabolites are removed based on Level and score, retaining compounds with the smallest level and highest similarity score (SS). Before statistical analysis, the serum metabolomics peak intensity data were log2-transformed. QC-based robust LOESS signal correction (QC-RLSC) (21) was utilized for data normalization to correct for systematic bias, and features with relative standard deviations (RSDs) of more than 35% in the QC sample were filtered out.

After normalization, the data were analyzed using the R package ropls for multivariate statistical analysis and modeling, including principal component analysis (PCA) and orthogonal partial least squares discriminant analysis (OPLS-DA). The quality of the model was tested by 7-fold cross-validation, and 20 permutation tests further tested the validity of the model. R2 and Q2 represent the explanatory and predictive abilities of the model, respectively. The variable importance for projection (VIP) value denotes the contribution of the feature peaks to the classification. The Wilcoxon rank-sum test was used to compare the metabolite levels between patients with tumors and healthy controls, and the false discovery rate (FDR) was used for p-value correction. FDR < 0.05 and VIP > 1 were used to screen significantly changed metabolites. Spearman’s correlation analysis was performed to analyze the associations between metabolites. Enrichment pathway analysis was performed using the web-based MetaboAnalyst 6.0 software.

Classification was strictly based on annotations from two authoritative databases: HNDB and MedChemExpress. HMDB’s “Biological Role” and “Origin” fields served as key references, with metabolites classified as endogenous if annotated as “endogenous” (i.e., produced by human cells or resident microbiota) and detected in “blood” or “plasma” according to the database’s tissue localization data. For metabolites with limited annotation in HMDB, MCE’s “Compound Type” classification was utilized, excluding those explicitly labeled as “exogenous” (e.g., dietary phytochemicals, pharmaceutical drugs, or environmental pollutants). Additionally, literature evidence was consulted to identify the main sources of metabolites.

To explore the potential added value of integrating epigenetic information, cfDNA methylation profiling was conducted on a subset of 197 samples (68 CRC, 129 NCC). Targeted bisulfite sequencing was used to detect CRC-associated hypermethylated loci. cfDNA was isolated from the serum using a Magen cfDNA extraction kit following the manufacturer’s instructions and then ligated to a methylation adaptor using an NEBNext Ultra II DNA library Prep Kit for Illumina from NEB. Adaptor-ligated cfDNA was 12-to-1 mixed and hybridized with customized probes (Integrated DNA Technologies) using an xGen hybridization capture DNA libraries Kit (Integrated DNA Technologies). Hybridized mixture samples were eluted using the reagents and steps of the ‘‘Second Elution’’ part of the TruSeq Methyl Capture EPIC Library Prep Kit (FC-151-1003, Illumina) and then bisulfite converted using an EZ-96 DNA Methylation-Lightning Mag Prep Kit (D5047, ZYMO RESEARCH). Bisulfate-converted samples were amplified using the reagents and steps of the ‘‘Amplify Enriched Library’’ section of the TruSeq Methyl Capture EPIC Library Prep Kit (FC-151-1003, Illumina). The concentration of the prepared libraries was determined with a Qubit 2.0 fluorometer (Invitrogen, Life Technologies), and the library quality was assessed by capillary electrophoresis (Qsep100, Bioptic). The qualified libraries were sequenced on the Illumina Nova-seq platform (Illumina).

Features with more than 50% missing values and samples with more than 20% missing values were filtered out, and then the median of abundance was used to fill in the missing values. After data filtering and normalization, the limma function performed a difference analysis on the normalized data. Benjamini-Hochberg correction was performed to reduce the bias caused by multiple tests, and the probes with p-values < 0.05 were filtered. Probes with Δβ > 0.1 or Δβ < -0.1 were considered hypermethylated or hypomethylated, respectively. The genes corresponding to the differential probes were annotated using the org.Hs.eg.db R package with UCSC.hg19 as the reference genome file.

Feature-level integration concatenated filtered metabolomic and methylation features into a unified matrix. To ensure compatibility between the two omics layers, composite data were normalized using M-value transformation for better interpretability in linear models. The datasets were scaled to zero mean and unit variance to prevent bias toward higher-magnitude features. The final integrated feature set was then input into the downstream diagnostic model.

In our study, all participants were randomly stratified and sampled into training and test datasets at a ratio of 7:3. Feature selection and model construction were performed on the training set, and hyperparameters were optimized through cross-validation. Subsequently, validation was carried out on the test set.

Features were selected using the least absolute shrinkage and selection operator (LASSO) and random forest (RF) algorithms. We performed LASSO regression on the training dataset to select the features with nonzero coefficients as a small number of features capable of identifying patients with tumors based on the average misclassification error of 10 random cross-validations (according to default settings). Then, the RF algorithm was used to calculate the relation importance of individual differential features, and the differential features were rounded off by performing 5-fold cross-validation five times (according to default settings). Finally, a Venn diagram was used to identify the common features. LASSO regression and RF modeling were implemented using the R packages ‘glmnet’ and ‘randomForest’, respectively.

Following feature acquisition, different machine learning methods, including support vector machine (SVM), RF, eXtreme Gradient Boosting (XGBoost), and logistic regression (LR) algorithms, were employed to create a model that could differentiate between cancer and noncancer individuals. The R packages ‘e1071’, ‘randomForest’, ‘xgboost’, and ‘caTools’ were utilized to implement SVM, RF, XGBoost, and LR, respectively. The models were trained with hyper-parameters tuned by 5-fold cross-validations. The final models were validated on the test set. To address potential biases in the optimization process that might arise from the specific splits employed, we replicated the above-described procedure five times, each time using a distinct splitting seed. This approach yielded five alternative optimized models to evaluate the initial model to ensure robustness. Ultimately, the efficacy of the four models was evaluated in test sets by AUROC, area under the precision–recall curve (AUPRC), sensitivity, specificity, accuracy, precision, recall and F1 score.

The overall workflow of this study and participant recruitment information are illustrated in Figure 1. A total of 715 participants were included in this study, comprising 248 CRC patients and 467 NCCs. The baseline characteristics of the study participants were summarized in Table 1. No statistical differences in gender, age, and BMI among the two groups (p > 0.05). The patients consisted of 50 stage I, 71 stage II, 82 stage III, 33 stage IV, and 12 unidentified stage (Table 1). Serum samples were collected and subjected to untargeted metabolomic profiling using LC-MS. Based on rigorous quality control and statistical filtering, differential metabolites between CRC patients and noncancer controls (NCCs) were identified and functionally annotated through pathway enrichment analysis. To assess their diagnostic potential, we applied multiple machine learning algorithms to construct predictive models using a selected panel of metabolites. In a subset of 197 participants, cfDNA methylation data were further integrated to evaluate the added value of a multi-omics classification strategy.

To characterize metabolic reprogramming associated with colorectal cancer, untargeted serum metabolomic profiling was performed on samples from CRC patients and NCCs. Following raw data conversion to mzXML format, preprocessing was carried out using XCMS, including peak detection, retention time correction, outlier removal, and imputation of missing values. A total of 4,030 ion features were identified in serum positive electrospray ionization (ESI+), while 674 features were identified in serum negative electrospray ionization (ESI−). Additionally, putative features were identified based on the Human Metabolome Database (HMDB) and KEGG, resulting in 2,793 annotated metabolites in ESI+ and 407 in ESI− mode.

Univariate analysis using the Wilcoxon rank-sum test and fold change (FC) filtering was performed across all identified metabolites. In parallel, multivariate modeling with orthogonal partial least squares discriminant analysis (OPLS-DA) revealed clear separation between CRC and NCC samples (Figure 2A). The OPLS-DA model demonstrated strong performance with R² = 0.932 and Q² = 0.764, and permutation testing confirmed the absence of overfitting (Supplementary Figure 1). These results indicated robust global metabolic differences between the two groups.

A total of 75 metabolites were significantly altered in CRC (FDR < 0.05 and VIP > 1) (Supplementary Table 1), among which 65 were upregulated (log₂FC > 1) and 10 were downregulated (log₂FC < –1), as shown in the volcano plot (Figure 2B). Metabolite classification includes benzenoids: 8, homogeneous non-metal compounds: 2, lipids and lipid-like molecules: 20, nucleosides, nucleotides, and analogues: 4, organic acids and derivatives: 10, organic nitrogen compounds: 2, organic oxygen compounds: 5, Organic Polymers: 1, Organic salts: 1, organoheterocyclic compounds: 16, Phenylpropanoids and polyketides: 5. To identify cancer-specific endogenous markers, we excluded exogenous metabolites, such as dietary components and drug metabolites, resulting in a refined panel of 26 differentially abundant endogenous metabolites (Supplementary Table 2).

To explore their interrelationships, Spearman correlation analysis was performed (Figure 2C), revealing extensive co-regulation among several metabolites. Notably, lipoyl-AMP (HMDB0059635), a precursor involved in lysine lipoylation, a rare but evolutionarily conserved post-translational modification, was significantly upregulated. In mammals, only four core metabolic enzymes undergo lipoylation, and dysregulation of these enzymes has been implicated in various metabolic disorders. While the role of lipoyl-AMP in cancer remains unclear, its involvement in central energy pathways highlights its potential clinical relevance (22).

Pathway enrichment analysis based on KEGG annotations further revealed that several dysregulated metabolites were enriched in biologically significant pathways, including primary bile acid biosynthesis and taurine and hypotaurine metabolism (Figure 2D). These pathways are known to be modulated by the intestinal microbiota and play pivotal roles in CRC pathogenesis (23). In particular, taurocholic acid and its derivatives may alter microbial composition and promote the production of genotoxic substances such as hydrogen sulfide and deoxycholic acid (24). These findings underscore the close interplay between host metabolism and microbiota-derived metabolites in colorectal cancer development and support the translational value of these metabolic alterations as potential diagnostic biomarkers.

To optimize the diagnostic utility of serum metabolites, we applied feature selection methods to reduce redundancy among the 26 differentially abundant metabolites. Random Forest (RF) importance ranking identified 11 top-ranked metabolites (), while Least Absolute Shrinkage and Selection Operator (LASSO) regression identified 21 metabolites with minimal binomial deviance (). The intersection of both methods yielded 10 candidate biomarkers for model construction (), including sulfate, ubiquinone-1, deoxycholic acid glycine conjugate, demethylphylloquinone, 3-(3,5-diiodo-4-hydroxyphenyl)lactate, vitamin K1, O-decanoyl-L-carnitine, sedoheptulose 1,7-bisphosphate, CE(22:6(4Z,7Z,10Z,13Z,16Z,19Z)), and lipoyl-AMP (). 1 1 1 1

To evaluate the diagnostic performance of these metabolites, we constructed binary classification models using both linear (logistic regression, LR) and nonlinear algorithms, including support vector machine (SVM), RF, and eXtreme Gradient Boosting (XGBoost). All samples were randomly split into a training set (70%) and a test set (30%), and the hyperparameters were selected based on tenfold cross-validation on the training set. Evaluation metrics included area under the receiver operating characteristic curve (AUROC), sensitivity, specificity, and overall accuracy (Figure 3). Compared to SVM (Figures 3A–C) and LR (Figures 3J–L) models, the RF (Figures 3D–F) and XGBoost (Figures 3G–I) models achieved better predictive performance. In the test cohort (74 CRC, 140 NCC), the RF model achieved robust performance, with an AUROC of 0.97, sensitivity of 83.78%, and specificity of 94.29% (Figures 3E, F). Similarly, the XGBoost model achieved an AUROC of 0.96 in the test set, with 89.19% sensitivity and 95.71% specificity (Figures 3H, I). These results demonstrate that the selected serum metabolite panel enables accurate discrimination between CRC patients and healthy controls, supporting its potential for noninvasive clinical screening applications. The results of 5 models, obtained through setting different seeds, are consistent with our initial model, which we ensure the robustness of the model (Supplementary Table 4).

In addition, we conducted additional analysis to evaluate whether the model distinguishes Stage I CRC from controls. The AUROC of SVM, RF, XGBoost and LR models is 0.84 (95% CI: 0.73,0.95), 0.87 (95% CI: 0.80, 0.94), 0.95 (95% CI: 0.88, 1.00) and 0.79 (95% CI: 0.66, 0.93), respectively, with low sensitivity and AUPRC (). Because of the number of Stage I CRC cases in our current dataset is limited, which may introduce variability and limit the generalizability of these findings. These results should be provided as preliminary evidence of the model’s potential in early-stage detection. 1

To investigate whether DNA methylation features could complement metabolomic biomarkers and further improve diagnostic accuracy, we performed targeted cfDNA methylation profiling in a subset of participants (CRC, n = 68; NCC, n = 129) (Figure 4A, Table 2). Among the 3,816 high-quality CpG probes retained after preprocessing, 426 sites (11.16%) were found to be differentially methylated between CRC patients and controls (FDR < 0.05), including 110 hypermethylated (log₂FC > 0.1) and 316 hypomethylated sites (log₂FC < –0.1) (Figure 4B, Supplementary Table 5). Chromosomal distribution analysis revealed enrichment of differentially methylated CpG sites (DMCs) on chromosome 8, while relatively few were observed on chromosome 18 (Figure 4C). Given the regulatory importance of promoter methylation in gene expression, we focused our downstream analysis on DMCs located within promoter regions (Supplementary Table 6).

Unsupervised hierarchical clustering of these promoter DMCs effectively discriminated CRC patients from controls (Figure 4D), suggesting their potential diagnostic utility. To identify a robust diagnostic methylation panel, we applied RF and LASSO regression to the training cohort (n = 138), resulting in the selection of 11 high-priority DMCs (Supplementary Figures 2D-F, Supplementary Table 7). Using these 11 features, we built classification models with SVM, RF, XGBoost, and LR algorithms (Figure 5). In the training cohort (48CRC, 90NCC), the sensitivity of SVM and LR models is 83.33%, with AUROC values of 0.993 and 0.96, respectively (Figures 5A, C, J, L). However, both RF and XGBoost models achieved 100% sensitivity and excellent AUROC values (Figures 5D, G, F, I). In the independent test cohort (20 CRC and 39 non-cancer controls), the SVM model achieved the highest sensitivity, reaching 85% (Figure 5B). Although the sensitivity of the RF and XGBoost models decreased to 75% and 80%, respectively, both maintained high specificity (Figures 5E, H). “In contrast, the Logistic Regression (LR) model demonstrated both sensitivity and specificity below 80% (Figure 5K). These results confirm the diagnostic relevance of cfDNA methylation markers but also highlight limitations in sensitivity, likely due to the smaller sample size and biological heterogeneity.

To evaluate the potential benefit of combining modalities, we next integrated the selected 10 serum metabolites with the 11 methylation markers to construct a multi-omics classifier, denoted as META&METHY. For benchmarking, we also retained the metabolomics-only (META) and methylation-only (METHY) models. Using four machine learning algorithms, RF, SVM, XGBoost, and LR, we systematically compared the performance of single-omics and multi-omics classifiers. As shown in Figure 6A and Supplementary Table 7, the META&METHY model achieved the highest diagnostic performance, with AUROCs of 0.98, AUPRC of 0.97 and 100% sensitivity and precision in the test set when using either the RF or XGBoost classifiers. F1 score comparisons further demonstrated the superiority of the META&METHY model over both the META and METHY models (Figure 6B). Notably, the metabolomics-only models consistently outperformed the methylation-only models, with AUROCs ranging from 0.87 to 0.97 versus 0.84 to 0.88, respectively (Supplementary Table 8). These findings suggest that while methylation features provide valuable supplementary information, serum metabolites are the dominant contributors to model performance.

In summary, the integration of cfDNA methylation with serum metabolomics modestly enhances predictive power, particularly in terms of sensitivity and classification confidence. However, the high performance of metabolomics-based models alone underscores their potential as a standalone, noninvasive diagnostic tool for colorectal cancer. Future studies with expanded cohorts are warranted to further validate and optimize this multi-omics approach for clinical application.

The original contributions presented in the study are included in the article/, further inquiries can be directed to the corresponding author/s. 1