Prevalence and chronology of colibactin-associated mutational processes and their microbiome spectra in Japanese colorectal cancer

Patients and samples

This study was approved by the Research Ethics Committees of the National Cancer Center, the University of Osaka and Institute of Science Tokyo as meeting the ethical guidelines for medical and health research involving human individuals (National Cancer Center Institutional Review Board, 2013-244; Research Ethics Committee, the University of Osaka, 20064-2: Human Subjects Research Ethics Review Committee, Institute of Science Tokyo, 2014018). Written informed consent was obtained from all participants before enrollment. All data were de-identified before analysis, and a double-anonymization procedure was applied to all participant identifiers.

This study enrolled a total of 600 individuals, comprising 454 patients with histologically confirmed CRC and 146 HCs, who were recruited at the National Cancer Center Hospital (Tokyo, Japan) and the University of Osaka Hospital (Osaka, Japan) between 2014 and 2025. Among these, cohort 1 consisted of cases prospectively collected specifically for the present study, whereas cohort 2 comprised cases previously described in our published work8 (Extended Data Fig. 1a).

These studies were designed to be observational and included participants for whom the biological specimens required for comprehensive analyses—such as stool samples, tumor tissues and matched normal DNA—were available between 2014 and 2025. Individuals with hereditary CRC syndromes, including familial adenomatous polyposis and hereditary non-polyposis CRC, or with inflammatory bowel disease, were excluded from the analysis.

WGS data were obtained from primary CRC tissues in 200 cases from cohort 1. WGMS data derived from stool samples were available for a total of 395 patients with CRC and 146 HCs across both cohorts. Detailed information on participants’ lifestyle and dietary habits was collected using a comprehensive self-administered questionnaire comprising 475 items across 25 pages. The questionnaire was developed based on the framework of the Japan Public Health Center-based Next Generation Study (JPHC-NEXT)46.

Stool sample collection and processing

Individuals classified as HC were those who presented with a positive fecal occult blood test during a hospital visit and subsequently underwent colonoscopy for secondary screening, with no clinically significant findings. Stool samples were obtained from HCs and CRC patients before the initiation of any treatment, including surgical resection or chemotherapy. The first stool passed at the hospital on the day of colonoscopy was collected for analysis47,48. In most cases, the stool samples collected on the day of colonoscopy were solid. Consistent with our previous observations47, the microbial taxonomic composition showed a high degree of concordance with that of standard frozen stool samples obtained without bowel preparation (Pearson’s correlation coefficient = 0.91, P < 0.01). All stool samples were immediately frozen on dry ice and subsequently stored at −80 °C until DNA extraction. Participants were instructed to consume a low-residue diet on the day preceding colonoscopy. On the day of examination, all participants received a bowel-cleansing agent followed by colonoscopic evaluation.

Tumor and matched normal sample collection

Tumor tissues were obtained from surgically resected specimens of patients with CRC. Matched non-tumor DNA was primarily extracted from peripheral blood lymphocytes; in selected cases, non-neoplastic colonic mucosa collected from regions anatomically distant from the primary tumor was also used as a normal control.

Clinical endpoint definitions

Overall survival was defined as the interval from the initiation of treatment to death from any cause, with censoring at the date of last confirmed follow-up for surviving patients (Supplementary Table 1). PFS was defined as the interval between the initiation of treatment and the first documented evidence of disease progression, as determined by serial radiological examinations (for example, computed tomography or magnetic resonance imaging) or by clinical confirmation of tumor recurrence or metastasis.

Microbial subtype classification and reproducibility assessment

For the combined test dataset (cohorts 1 + 2), comprising 199 CRC cases from our previous publication8 and 196 CRC cases newly collected from the National Cancer Center Hospital and the University of Osaka Hospital, we applied a previously published microbial subtype classification algorithm16. Cases were categorized into four microbial subtypes, and their bacterial community compositions were compared to assess reproducibility.

Baseline clinicopathological characteristics were summarized across subtypes, and overall survival was analyzed using the Kaplan–Meier method with comparisons performed with the log-rank test. These univariate comparisons assessed differences in survival that may reflect both subtype-related and background clinical variations. To determine whether microbial subtypes independently influence prognosis, we further applied multivariate Cox proportional hazards models adjusted for established prognostic covariates.

WGMS of fecal samples

Genomic DNA was extracted from frozen stool samples using a bead-beating protocol49 in conjunction with the GNOME DNA Isolation Kit (MP Biomedicals). DNA integrity and concentration were evaluated using the Agilent 4200 TapeStation system (Agilent Technologies). Following ethanol precipitation, the purified DNA was resuspended in TE buffer and stored at −80 °C until library preparation. Sequencing libraries were constructed from fecal DNA using the Nextera XT DNA Sample Prep Kit (Illumina) according to the manufacturer’s instructions. WGMS was performed on the NovaSeq 6000 platform (Illumina) with paired-end (PE) reads (2 × 150 bp), targeting an average sequencing depth of 5.0 Gb per sample.

Metagenomic sequencing quality control and taxonomic profiling

In total, 9,875,823,911 PE reads (51,843,293 on average) from cohort 1 and 10,156,584,408 PE reads (51,038,112 on average) from cohort 2 were generated from 150-bp sequencing and subjected to stringent quality control. Raw reads containing undetermined bases (‘N’) were removed. Reads containing bacteriophage phiX sequences were identified and filtered by alignment to the reference genome using Bowtie2 (v.2.2.9) with the ‘–fast-local’ preset. Adaptor and primer sequences were trimmed using Cutadapt (v.1.9.1) with the following parameters: for forward reads, -a CTGTCTCTTATACACATCTCCGAGCCCACGAGAC -O 33 -q 17; for reverse reads, -a CTGTCTCTTATACACATCTGACGCTGCCGACGA -O 32 -q 17. Within Cutadapt, reads with consecutive quality scores of ≤17 were trimmed at the 3′ end, and reads shorter than 50 bp after trimming were discarded. Reads with an average Phred quality score of ≤25 were also excluded. High-quality reads were aligned to the human reference genome (GRCh38; gi 568336000–568336023) using Bowtie2 (v.2.2.9), and all reads mapping to the human genome were removed. Unpaired reads were subsequently discarded. After these filtering steps, a total of 8,777,929,480 PE reads (46,899,434 on average) from cohort 1 and 8,896,915,982 PE reads (44,708,121 on average) from cohort 2 were retained for downstream analyses (hereafter referred to as high-quality reads). Quality-controlled reads were processed using MetaPhlAn3 with default parameters to generate species-level and genus-level taxonomic profiles50.

Identification of CRC subtypes based on metagenome-derived SHAP valuesConstruction of a CRC classifier using a random forest model

A random forest model was constructed to estimate the likelihood of CRC, hereafter referred to as the ‘microbial CRC trait score’, which was calculated for each sample. The model was trained on four independent metagenomic cohorts7,17,18,19 curated through the curatedMetagenomicData R package20, applying filtering thresholds of 1 × 10−5 for taxonomic abundance and 0.95 for prevalence. Model training and evaluation were performed using the scikit-learn library (v.1.2.2). Model performance was assessed by tenfold cross-validation, yielding an AUC of 0.84. The trained classifier was subsequently applied to MetaPhlAn3-derived species-level relative abundance profiles from 196 CRC samples (cohort 1) and 146 HCs previously reported8. This evaluation yielded an AUC of 0.74, confirming consistent discriminatory performance in an independent dataset (Extended Data Fig. 1b).

SHAP analysis and subtype identification

Interpretable artificial intelligence frameworks such as SHAP and LIME (local interpretable model-agnostic explanations) have been developed to elucidate the internal decision-making mechanisms of complex predictive models51. Among these, SHAP provides a mathematically rigorous and reproducible framework for feature attribution, quantifying the contribution of each input variable to the model’s output. SHAP analysis is grounded in cooperative game theory, in which Shapley values represent the fair distribution of contributions among all predictors. In this study, SHAP was used to interpret the internal decision process of the random forest model and to identify distinct microbial features associated with CRC. To further characterize CRC heterogeneity, unsupervised clustering of the SHAP values was performed using the k-means algorithm. The optimal number of CRC subtypes was determined by minimizing the within-cluster sum of squares (Extended Data Fig. 2a). Given that SHAP does not inherently account for potential confounding factors such as age, tumor stage or primary tumor site, we systematically examined these variables across SHAP-derived clusters to evaluate potential confounding influences. The reproducibility of the SHAP-based subtype classification was validated in an independent cohort of 199 CRC cases (cohort 2) (Extended Data Fig. 2b,c). All analyses were performed using the shapmat Python package (https://github.com/ryzary/shapmat).

WGS analysis of CRC

Genomic DNA was extracted from primary tumor tissues and matched peripheral blood samples. WGS libraries with an average insert size of 550 bp were prepared from 2 µg of genomic DNA using the TruSeq DNA PCR-Free Library Preparation Kit (Illumina). Sequencing was performed on the NovaSeq 6000 platform (Illumina) with 150 bp PE reads. The median sequencing depth was 52× (mean, 56×) for tumor samples and 32× (mean, 34×) for matched normal samples. The median estimated tumor purity was 0.52 (range, 0.11–0.97), and all cases met the inclusion criterion of ≥10% tumor purity.

Mutation calling

PE reads were aligned to the human reference genome (GRCh37) using BWA-MEM52. PCR duplicates were removed by eliminating paired reads that mapped to identical genomic coordinates, and pileup files were subsequently generated using SAMtools53. Somatic point mutations, including SNVs and short indels, were identified according to the following eight criteria: (1) mapping quality of ≥20 and (2) base quality of ≥10. Candidate somatic mutations were further filtered by applying the following conditions: (3) in each tumor sample, variants were required to be supported by at least four reads when the tumor variant allele frequency (TVAF) was ≥0.15, or by at least eight reads when 0.15 > TVAF ≥ 0.05, with at least one supporting read having a base quality ≥30; (4) the variant allele frequency (VAF) of the matched non-tumor sample had to be <0.03, with a minimum read depth of eight. To reduce context-dependent sequencing errors, reads from all non-tumor samples were aggregated to identify recurrent false-positive sites. For each genomic position with a non-tumor sequence depth of ≥10 and VAF < 0.2, the aggregate non-tumor VAF (NVAF) was calculated. Variants were retained only if (5) NVAF < 0.03 for sites with TVAF ≥ 0.15 or NVAF < 0.01 for sites with 0.15 > TVAF ≥ 0.05, and (6) the ratio of TVAF to NVAF was ≥20. To exclude germline single-nucleotide polymorphisms, (7) the proportion of non-tumor samples with a VAF ≥ 0.1 was required to be <0.002. Finally, (8) variants showing extreme strand bias (>95% of reads derived from a single strand) were excluded.

Copy number analysis

Copy number ratios, tumor purity and ploidy were estimated using Battenberg (v.2.2.10)54 with default parameters. The inferred copy number ratios were subsequently adjusted according to the estimated tumor purity for each sample to obtain purity-corrected copy number profiles.

Significantly mutated gene analysis

Significantly mutated genes were identified using three complementary approaches: the inactivation bias method22, the activation bias method and dNdScv55. For the inactivation bias test, the number of samples harboring inactivating mutations (nonsense, read-through, splice-site or frameshift variants) was compared with the number carrying alternative mutations using Fisher’s exact test. Conversely, for the activation bias test, the number of activating (hotspot missense) mutations was compared with the number of other mutations using Fisher’s exact test. Genomic positions exhibiting at least two identical mutations were defined as mutational hotspots. Mutations occurring within 5 bp of a hotspot were also classified as hotspot mutations. As hotspot missense mutations are generally activating but may occasionally be inactivating or functionally ambiguous, genes with P < 0.05 by the inactivation bias method were designated as putative tumor suppressor genes, and their activation bias P values were fixed at 1. Multiple hypothesis testing was corrected using the Benjamini–Hochberg false discovery rate procedure, and adjusted P values are reported as q values. Genes were considered significantly mutated if any of the three tests yielded a q value of <0.1.

Mutational signature analysis

De novo decomposition of SBS and small ID mutational signatures was performed using SigProfilerExtractor on WGS data from 1,153 CRC cases, comprising 200 Japanese CRC cases from the present study and additional non-Japanese CRC cases reported previously15. Given that mutational spectra differ substantially between hypermutated and non-hypermutated tumors, de novo decomposition was conducted separately for these two groups, following the approach described in a previous study15. The number of decomposed signatures was determined based on the optimal solution estimated by SigProfilerExtractor. For each decomposed signature, cosine similarity was calculated relative to the COSMIC reference signature set. Signatures with cosine similarity of ≥0.75 were assigned to known COSMIC signatures, whereas those with cosine similarity of <0.75 were interpreted as novel or undecomposed composite signatures and classified as undetermined.

For the 200 CRC cases analyzed in this study, reconstruction accuracy was evaluated by comparing the reconstructed and observed mutational spectra. SBS signatures demonstrated consistently high reconstruction fidelity (cosine similarity of ≥0.97 in all cases). By contrast, two cases exhibited ID signature cosine similarity values of <0.9, both of which contained a very small number of indels (44 and 68, respectively), indicating insufficient data for reliable signature decomposition. These two cases were therefore excluded from subsequent ID signature analyses.

Definition of colibactin-associated mutational signatures

The presence of a colibactin-associated mutational signature was defined as a relative contribution of SBS88 of ≥5% or ID18 of ≥5%. This threshold was adopted because ten out of 82 colibactin-exposed cases exhibited only one of the two signatures (SBS88 or ID18). Moreover, a 5% cutoff has been widely applied in prior investigations, including the first report characterizing the colibactin-associated mutational process12. Given that mutational signature decomposition algorithms may artifactually assign low-frequency components, contributions below 5% are generally regarded as unreliable and are typically truncated to zero, whereas contributions ≥5% are considered robust signals.

Expected probability of mutational signatures

The expected probability of each mutational signature for a given trinucleotide context in an individual patient was calculated using the following equations:

$${N}_{{ijk}}={S}_{{ik}}\times {C}_{{ij}}$$

$${P}_{{ijk}}={N}_{{ijk}}/\mathop{\sum }\limits_{i}{N}_{{ijk}}$$

where \({S}_{{ik}}\) denotes the number of mutations attributed to mutational signature i in patient k, \({C}_{{ij}}\) represents the proportion of trinucleotide context j within mutational signature i, \({N}_{{ijk}}\) is the estimated number of mutations corresponding to mutational signature i and trinucleotide context j in patient k, and \({P}_{{ijk}}\) is the expected probability of observing mutational signature i in trinucleotide context j for patient k (ref. 27).

Clonal analysis of mutational signatures

Mutations were categorized as clonal or subclonal using MutationTimeR (v.1.00.2)33. Mutations that occurred before copy number gains were annotated as clonal (early); those that arose following copy number gains were annotated as clonal (late); and mutations for which timing could not be determined were designated as clonal (NA). For each mutation group, mutational signature decomposition was performed using the deconstructSigs R package (v.1.9.0)56. The parameter signature ‘signature.cutoff’ was set to zero, and ‘signatures.ref’ was assigned to the SBS and ID signatures extracted by SigProfilerExtractor, as described above. Clonal (early), clonal (late) and clonal (NA) categories were collectively defined as clonal mutations. Cases with fewer than 50 mutations in any group were excluded from the analysis to ensure robust signature estimation.

Timing analysis of copy number gains

The relative timing of three classes of copy number gains was inferred from (1) CNN-LOH / loss + gain (N:0), (2) monoallelic gains (N:1) and (3) biallelic gains (N:M). Analyses were performed using MutationTimeR (v.1.00.2)33, which integrates somatic SNVs and copy number variation data derived from WGS. For each tumor, copy number timing estimates were generated using genomic segments containing at least 100 somatic mutations to ensure reliable inference. Weighted averages of mutation counts were computed across segments for each copy number gain category. Comparative analyses among the four CRC subtypes were performed after exclusion of hypermutated cases. To further reconstruct the temporal sequence of genomic events within individual tumors, PhylogicNDT SinglePatientTiming (v.1.0; https://doi.org/10.1101/508127) was applied. The inferred relative timing of events was scaled from zero, representing the zygotic (fertilization) stage, to one, corresponding to the time of clinical diagnosis.

Transcriptome sequencing (RNA-seq)

Total RNA was extracted from fresh–frozen tumor tissues using the miRNeasy Mini Kit (Qiagen) according to the manufacturer’s protocol. RNA integrity was assessed using the BioAnalyzer system (Agilent Technologies), and samples with an RNA integrity number of >6.0 were selected for sequencing. RNA-seq libraries were prepared using the SureSelect Strand-Specific RNA Library Preparation Kit (Agilent Technologies) with 300 ng of total RNA as input. Sequencing was performed on the HiSeq 2500 platform (Illumina) to generate 101 bp PE reads. Reads were aligned to the human reference genome (GRCh37/hg19) using STAR57, and gene-level read counts were quantified against the UCSC human transcriptome reference using HTSeq58. Expression levels were calculated as fragments per kilobase of exon per million mapped fragments, incorporating strand-specific information. Differential gene expression was assessed using the Wilcoxon rank-sum test.

Estimation of non-human bacterial reads in tumor tissues

We used the kneaddata tool (v.0.12.0)59 to eliminate human reads from human rRNA-depleted total RNA sequencing data of tumor tissues, encompassing human and bacterial transcripts. Relative abundance of bacterial reads was calculated as the ratio of non-human reads to all reads in each sample. The total RNA was extracted from fresh–frozen tumor tissues using the miRNAeasy kit (Qiagen). A total of 250 ng of total RNA was used to generate libraries for total RNA sequencing. The Universal Plus Total RNA-seq kit with human rRNA depletion (Tecan) was used for library preparation, followed by sequencing on a NovaSeq 6000 in PE 150 base mode.

Decomposition of tumor immune environments

Gene expression profiles, annotated using GENCODE v43lift37, were deconvoluted to estimate tumor immune cell composition. The proportions of 22 immune cell types were inferred using CIBERSORTx based on the LM22 reference signature matrix, which includes naive and memory B cells, plasma cells, CD8+ T cells, naive and memory CD4+ T cells (resting and activated) and follicular helper T cells. Additional immune cell subsets analyzed comprised Treg cells, γδ T cells, resting and activated natural killer cells, monocytes, M0, M1 and M2 macrophages, resting and activated dendritic cells, resting and activated mast cells, eosinophils and neutrophils. CIBERSORTx was executed in ‘B-mode’ with the ‘disable quantile normalization’ option enabled, as recommended by the developers for RNA-seq data. Only samples with a CIBERSORTx output P value of <0.05 were retained for downstream analyses, ensuring robust estimation of immune cell fractions.

T cell-inflamed GEP analysis

The T cell-inflamed GEP score60 was calculated as the weighted sum of the normalized expression values of 18 genes, each multiplied by its corresponding coefficient: CCL5 (0.008346), CD27 (0.072293), CD274 (0.042853), CD276 (−0.023900), CD8A (0.031021), CMKLR1 (0.151253), CXCL9 (0.074135), CXCR6 (0.004313), HLA-DQA1 (0.020091), HLA-DRB1 (0.058806), HLA-E (0.071750), IDO1 (0.060679), LAG3 (0.123895), NKG7 (0.075524), PDCD1LG2 (0.003734), PSMB10 (0.032999), STAT1 (0.250229) and TIGIT (0.084767). The resulting composite score reflects the degree of T cell-mediated inflammation within the tumor microenvironment.

Quantification of the clbP gene by qPCR

Quantification of the clbP gene (copies per ml) in genomic DNA extracted from tumor tissues or fecal samples was performed by Adenoprevent32. qPCR was conducted using the KAPA SYBR Fast qPCR Kit (Roche Molecular Systems) on the Thermal Cycler Dice Real-Time System (Takara Bio). The copy number of the clbP gene was determined from cycle threshold values using a standard calibration curve generated from clbP gene standards (Adenoprevent)30,31.

Immunohistochemistry

Formalin-fixed, paraffin-embedded tissue sections (4 μm thick) were deparaffinized and treated with 3% hydrogen peroxide for 15 min to quench endogenous peroxidase activity. Antigen retrieval was performed using an autoclave in 10 mM citrate buffer (pH 6.0) for 10 min. Slides were incubated with a mouse monoclonal anti-FOXP3 antibody (clone 236 A/E7; dilution 1:100; ab20034, Abcam) for 1 h at room temperature (20–25 °C), followed by detection using the EnVision system (Dako) with mouse linker (Dako), according to the manufacturer’s protocol. Signal visualization was achieved using the Chem-Mate EnVision detection method (Dako). Immunohistochemical staining for FOXP3 was performed in 28 CRC cases classified as subtype 3. For each case, five representative tumor regions exhibiting lymphocytic infiltration were randomly selected and imaged at ×200 magnification. FOXP3-positive lymphocytes were manually counted within each field.

Statistical analyses

When multiple comparison groups were present, overall differences were assessed using the Kruskal–Wallis test. Pairwise comparisons between individual groups were subsequently performed using the Wilcoxon rank-sum test or Fisher’s exact test, as appropriate. Correction for multiple hypothesis testing was performed using the Benjamini–Hochberg false discovery rate procedure, and adjusted P values are reported as q values.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.