The colon cancer microbiome risk score for diagnosis of colon cancer

Abstract

Methods for predicting survival in patients suffering from colon cancer are described.

Description

THE COLON CANCER MICROBIOME RISK SCORE FOR DIAGNOSIS OF COLON CANCERCROSS REFERENCE TO RELATED APPLICATIONS

[0001] This application claims the benefit of, and priority to, U.S. Provisional Application No. 63 / 559,612, filed on February 29, 2024, which is incorporated herein by reference in its entirety.SEQUENCE LISTING

[0002] This application contains a sequence listing having the filename 5600461- 00024_SL.xml, which is 5 KB in size, and was created on February 21 , 2025. The entire content of this sequence listing is incorporated herein by reference.FIELD

[0003] This disclosure is related to testing for colon cancer.BACKGROUND

[0004] Identifying accurate biomarkers that correlate with clinical outcomes in colon cancer can be difficult because diagnostic tests available detect presence or absence of cancer but do not provide information about prognosis for a patient. The lack of multi-omics cancer datasets with extensive follow-up information hinders the identification of accurate biomarkers of clinical outcomes in colorectal cancer. There is a need for better understanding of colon cancer biology which can facilitate the discovery of personalized therapeutic approaches and prediction of clinical outcomes.SUMMARY

[0005] Described herein are diagnostic and prognostic methods comprising combining a patient’s microbiome signature, e.g., presence of Ruminococcus bromii, with an Immunologic Constant of Rejection, to arrive at a composite score (mICRoScore), which identifies a group of patients as having higher probability of survival compared to patients who lack the microbiome signature (e.g., presence of Ruminococcus bromii).

[0006] In one aspect, provided is a method for predicting survival in patients suffering from colon cancer, the method comprising: a) sequencing 16S rRNA genes using DNA extracted from tumor samples of patients suffering from colon cancer;b) detecting the presence of Ruminococcus bromii DNA sequences in the tumor samples of patients suffering from colon cancer; and c) identifying patients suffering from colon cancer having higher levels of Ruminococcus bromii DNA sequences in the tumor samples as having a higher probability of survival.

[0007] In another aspect, provided is a method for predicting survival in patients suffering from colon cancer, the method comprising: a) measuring expression levels of the genes CCL5, CD274, CD8A, CD8B, CTLA4, CXCL9, CXCL10, FOXP3, GNLY, GZMA, GZMB, GZMH, IDO1 , IFNG, IL12B, IRF1 , PDCD1 , PRF1 , STAT 1 , and TBX21 in tumor samples of patients suffering from colon cancer and assigning an immunologic constant of rejection (ICR); b) detecting the abundance of specific microbes in the tumor samples of patients suffering from colon cancer and assigning a microbiome signature risk (MBR) score; and

[0008] c) identifying patients suffering from colon cancer who have a high ICR and low MBR score as having a higher probability of survival. In a further aspect, provided is a method of identifying a patient as suffering from colon cancer, the method comprising detecting the presence of one or more genes selected from SEC31A, BAX, ZMYM2, BCR, ATF7IP, BCL1 B, ERCC5, or PARK2 in a tumor sample of the patient.

[0009] Further provided is the use of diagnostic methods described herein for identifying and / or treating a patient suffering from colon cancer.BRIEF DESCRIPTION OF THE DRAWINGS

[0010] FIG. 1A describes the AC-ICAM study. Samples from a total of 348 patients with colon cancer were included in AC-ICAM. The number of profiled samples and resulting analytes are indicated for each platform, including RNA-seq, WES, TCR sequencing (immunoSEQ TCRp assay), 16S rRNA gene sequencing and metagenomic analysis from whole-genome sequencing (WGS) to profile microbiome composition. An additional 42 tumor samples were profiled with 16S rRNA gene sequencing that did not have any matched normal tissue available (ICAM42).

[0011] FIG. 1 B-1 Hazard Ratios (HRs) of deconvoluted abundancies of distinct infiltrating cell populations by Single Sample gene set enrichment (ssGSEA) with Overall Survival (OS) and Progression Free Survival (PFS). Corresponding 95% confidence intervals (error bars) ascalculated by Cox proportional hazard regression are displayed as a forest plot (n = 346 independent samples from 346 patients).

[0012] FIG. 1 B-2 provides P values for the associated HRs of FIG 1 B-1 as indicated in the bar chart (— log 10 P value).

[0013] FIG. 2A depicts correlation between immune gene signatures and TCR metrics from immunoSEQ DNA sequencing. With ICR having the highest correlation with productive clonality compared to ssGSEA scores for deconvoluted cell populations.

[0014] FIG. 2B depicts correlation of proportion of tumor-enriched T cell clones in the tumor (in percent) versus normal paired tissue, with ICR score. Pearson’s r and P value of the correlation are indicated in the plot. All P values are two-sided.

[0015] FIG. 3 depicts ESTIMATE scores in AC-ICAM (tumor tissue and matched healthy colon tissue, AC-ICAM cohort, see Methods section) and in The Cancer Genome Atlas (TCGA) colon adenocarcinoma (COAD) (TCGA-COAD) cohorts. Unpaired two-sided Student’s t-test.

[0016] FIG. 4A depicts Kaplan-Meier OS curve for the combination of ICR cluster and mutational load category. Mutational load high is defined as nonsynonymous mutation frequency of >12 per Mb. Overall P value is calculated by log-rank test.

[0017] FIG. 4B depicts a scatter-plot of ICR score by genetic immunoediting (GIE) value for ICR-high and ICR-low samples. Numberof samples in each quadrant is indicated in the graph. Gray area delineates ICR scores from 5-9.

[0018] FIG. 4C depicts Kaplan-Meier plots for OS by IES. Censor points are indicated by vertical lines and corresponding table of number of patients at risk in each group is included below the Kaplan-Meier plot. Overall P value is calculated by log-rank test.

[0019] FIG. 5A depicts coefficients of the 41 taxa in the MBR classifier as selected by the OS elastic-net Cox regression model. Family is indicated between parentheses. *means taxonomical order is indicated between brackets, as family was unassigned (uncultured).

[0020] FIG. 5B depicts Forest plot showing the HR, 95% confidence intervals (error bars) and corresponding P value calculated by Cox proportional hazard regression analysis for OS of the 16S or WGS MBR classifier scores in training and test sets.

[0021] FIG. 5C depicts MBR score in tumorversus healthy paired colon tissue. The gray band reflects the 95% confidence interval for predictions of the linear regression model betweenthe plotted variables. P value for Spearman correlation for relative abundance and P value for Pearson correlation for MBR scores are indicated. OS. All P values are two-sided; n reflects the independent number of samples.

[0022] FIG. 6A depicts Kaplan-Meier curves of OS by mICRoScore in AC-ICAM. Overall P value is calculated by log-rank test. Vertical lines indicate censor points. HRs and 95% confidence intervals are calculated by Cox proportional hazard regression. All P values are two-sided; n reflects the independent number of samples.

[0023] FIG. 6B depicts Kaplan-Meier curves of OS by mICRoScore in TCGA-COAD. Overall P value is calculated by log-rank test. Vertical lines indicate censor points. HRs and 95% confidence intervals are calculated by Cox proportional hazard regression. All P values are two-sided; n reflects the independent number of samples.

[0024] FIG. 6C depicts Kaplan-Meier curve of OS in ICR-high samples by mICRoScore in AC-ICAM. Overall P value is calculated by log-rank test. Vertical lines indicate censor points. HRs and 95% confidence intervals are calculated by Cox proportional hazard regression. All P values are two-sided; n reflects the independent number of samples.

[0025] FIG. 6D depicts Kaplan-Meier curve of OS in ICR-high samples by mICRoScore in TCGA-COAD. Overall P value is calculated by log-rank test. Vertical lines indicate censor points. HRs and 95% confidence intervals are calculated by Cox proportional hazard regression. All P values are two-sided; n reflects the independent number of samples.DETAILED DESCRIPTION

[0026] Current guidelines for cancer treatment rely on the tumor-node-metastasis staging and the detection of DNA mismatch repair (MMR) deficiency or microsatellite instability (MSI), in addition to standard clinicopathological variables, to determine treatment recommendations. MSI is caused by somatic or germline defects of MMR genes and leads to the accumulation of somatic mutations, neoantigens resulting in immune recognition, and high density of tumor infiltrating lymphocytes. The Immunologic Constant of Rejection signature (ICR) is a 20-gene expression signature of cytotoxic immune response with prognostic value in some solid cancers. Another variable used is the evaluation of the density and spatial distribution of T cells (Immunoscore), which is associated with a reduced risk of relapse and death independently of other clinicopathological variables, including MSI status. The Immunoscore captures the strength of the in situ adaptive immune reaction. Though such quantitative phenotypic features of primary colon cancer, including those that are cancer cell intrinsic,immunological, stromal or microbial in nature, have been reported to be associated with clinical outcomes individually, it is not clear how their interactions impact patient outcomes.

[0027] Described herein are multi-omic analyses from orthogonal genomic platforms to rigorously profile a large collection of primary colon cancer specimens (unselected for tumor cell purity) and matched healthy colon tissue, complemented with clinical and pathological data follow-up.

[0028] “Immunologic Constant of Rejection” or “ICR” is a classifier based on consensus clustering (CC) analysis of the expression levels of 20 immune genes (namely, CCL5, CD274, CD8A, CD8B, CTLA4, CXCL9, CXCL10, FOXP3, GNLY, GZMA, GZMB, GZMH, IDO1 , IFNG, IL12B, IRF1 , PDCD1 , PRF1 , STAT1 , and TBX21). ICR typically correlates with overall survival (OS) and progression free survival (PFS). A high expression of these genes indicates an active immune engagement, i.e., at least a partial rejection of the cancer tissue.

[0029] “Microbiome signature risk (MBR) score” refers to a calculation using the coefficients shown in FIG. 5A, which are weights. The MBR is a weighted average based on abundance of microbes (such as the microbes described in FIG. 5A). A higher positive coefficient in FIG. 5A means high hazard risk of death, whereas a negative coefficient corresponds to lower risk of death for that microbe. Thus, by way of example only, R Bromii has the largest negative weight on the MBR score calculation because R Bromii has the largest negative coefficient.

[0030] “Genetic immune editing” is calculated as the ratio of the observed versus the expected number of neoantigens in a tumor sample. When fewer neoantigens than what would be expected based on the tumor’s mutational burden are found, a tumor has been immune edited.

[0031] “Stromal material” in a tumor sample refers to supporting connective tissue which may be present in a tumor environment at the time of obtaining a tumor sample.

[0032] “Immune content” in a tumor sample refers to proteins such as, and not limited to, antibodies, immunoglobulins, or cells such as, and not limited to, infiltrated lymphocytes.

[0033] “Tumor sample purity” can be determined by ESTIMATE (Estimation of STromal and Immune cells in MAIignant Tumor tissues using Expression data) which is a tool for predicting tumor purity, and the presence of infiltrating stromal / immune cells in tumor tissues using gene expression data. ESTIMATE algorithm is based on ssGSEA and generates three scores: a stromal score (that captures the presence of stroma in tumor tissue), an immune score (thatrepresents the infiltration of immune cells in tumor tissue), and an estimate score (that infers tumor purity).

[0034] An Immune Constant of Rejection (ICR) signature is generally associated with metastasis-free survival and / or pathological response to chemotherapy. “High ICR” in patients suffering from colon cancer refers to a high expression of a panel of genes (CCL5, CD274, CD8A, CD8B, CTLA4, CXCL9, CXCL10, FOXP3, GNLY, GZMA, GZMB, GZMH, IDO1 , IFNG, IL12B, IRF1 , PDCD1 , PRF1 , STAT1 , and TBX21) that are used as a genetic signature to assign an immune constant of rejection (ICR). A high ICR is indicative of an active immune engagement, and at least a partial rejection of the cancer tissue. In some embodiments, the high expression of genes means that the expression levels of the genes are at least 10%, at least 15%, at least 20%, at least 25%, at least 30%, at least 35%, at least 40%, at least 45%, at least 50%, at least 55%, at least 60%, at least 65%, at least 70%, at least 75%, or at least 80% higher than the levels of expression of the same genes in healthy colon tissue. ICR can be used as a constant variable by summing up the gene expression or using ssGSEA, or it can be used to cluster a population using consensus clustering.

[0035] Provided is a method for predicting survival in patients suffering from colon cancer, the method comprising sequencing 16S rRNA genes using DNA extracted from tumor samples of patients suffering from colon cancer. In some embodiments, the method comprises identifying the microbiome profiles of the patients suffering from colon cancer. The method further comprises detecting the presence of Ruminococcus bromii DNA sequences in the tumor samples of patients suffering from colon cancer. The method further comprises comparing levels of Ruminococcus bromii DNA sequences in the tumor samples of patients suffering from colon cancer with levels of Ruminococcus bromii DNA sequences in healthy colon tissue samples. Patients suffering from colon cancer who have increased levels of Ruminococcus bromii DNA sequences in tumor samples are identified as having a higher probability of survival. In some embodiments, the levels of Ruminococcus bromii may be compared with levels of Ruminococcus bromii in healthy colon tissue samples to identify patients having increased levels of Ruminococcus bromii DNA sequences in tumor samples.

[0036] Provided is a method for predicting survival in patients suffering from colon cancer, the method comprising a) measuring expression levels of the genes CCL5, CD274, CD8A, CD8B, CTLA4, CXCL9, CXCL10, FOXP3, GNLY, GZMA, GZMB, GZMH, IDO1 , IFNG, IL12B, IRF1 , PDCD1 , PRF1 , STAT1 , and TBX21 in tumor samples of patients suffering from colon cancer and assigning an immunologic constant of rejection (ICR). The method furthercomprises b) detecting the abundance of specific microbes in the tumor samples of patients suffering from colon cancer and assigning a microbiome signature risk (MBR) score. The method further comprises c) identifying patients suffering from colon cancer who have a high ICR, and who also have a low MBR score as having a higher probability of survival. In some embodiments, the abundance of specific microbes is measured by detecting the presence of genes associated with the microbes.

[0037] In some embodiments of the method, in step b), the abundance of one or more microbes from genera UBA1819 (Ruminococcaceae), Lachnospiraceae NK4A136 group (Lachnospiraceae), Erysipelotrichaceae UCG-003 (Erysipelotrichaceae), Sellimonas (Lachnospiraceae), Coprococcus 1 (Lachnospiraceae), Terrisporobacter (Peptostreptococcaceae), Olsenella (Atopobiaceae), Oscillibacter (Ruminococcaceae), Dialister (Veillonellaceae), Parabacteroides (Tannerellaceae), Parvimonas (Family XI), [Eubacterium] eligens group (Lachnospiraceae), Bifidobacterium (Bifidobacteriaceae), Hungatella (Lachnospiraceae), Mogibacterium (Family XIII), Caproiciproducens (Ruminococcaceae), Erysipelatoclostridium (Erysipelotrichaceae), Halomonas (Halomonadaceae), and Odoribacter (Marinifilaceae) is detected.

[0038] In some embodiments of the method, in step b), the presence of one or more microbes from genera Ruminococcaceae NK4A214 group (Ruminococcaceae), Solobacterium (Erysipelotrichaceae), Ruminococcus 2 (Ruminococcaceae), Treponema 2 (Spirochaetaceae), Streptococcus (Streptococcaceae), Phascolarctobacterium (Acidaminococcaceae), Moryella (Lachnospiraceae), Coprococcus 2 (Lachnospiraceae), Ruminococcus 1 (Ruminococcaceae), [Eubacterium] ventriosum group (Lachnospiraceae), Actinomyces (Actinomycetaceae), CAG-56 (Lachnospiraceae), Ruminiclostridium 6 (Ruminococcaceae), Barnesiella (Barnesiellaceae), Prevotella 9 (Prevotellaceae), Family XIII AD3011 group (Family XIII), Prevotellaceae, Peptostreptococcus (Peptostreptococcaceae), Candidatus Soleaferrea (Ruminococcaceae), Rhodospirillales, Leptotrichia (Leptotrichiaceae), Veillonella (Veillonellaceae), UBA1819 (Ruminococcaceae), Lachnospiraceae NK4A136 group (Lachnospiraceae), Erysipelotrichaceae UCG-003 (Erysipelotrichaceae), Sellimonas (Lachnospiraceae), Coprococcus 1 (Lachnospiraceae), Terrisporobacter (Peptostreptococcaceae), Olsenella (Atopobiaceae), Oscillibacter (Ruminococcaceae), Dialister (Veillonellaceae), Parabacteroides (Tannerellaceae), Parvimonas (Family XI), [Eubacterium] eligens group (Lachnospiraceae), Bifidobacterium (Bifidobacteriaceae), Hungatella (Lachnospiraceae), Mogibacterium (Family XIII),Caproiciproducens (Ruminococcaceae), Erysipelatoclostridium (Erysipelotrichaceae), Halomonas (Halomonadaceae), and Odoribacter (Marinifilaceae) is detected.

[0039] In some embodiments, the method further comprises detecting the presence of genetic immunoediting (GIE), where the presence of GIE identifies a patient suffering from colon cancer as having a higher probability of survival. In some embodiments, the method further comprises detecting the presence of TH1 cell / cytotoxic immune activation and / or detecting concurrent expansion of TCR clonotypes.

[0040] In some embodiments the method does not depend on tumor sample purity, e.g., the method can be performed on crude tumor samples without any further processing even if the tumor samples include infiltrated cells such as lymphocytes or if the tumor samples include collagen or other connective tissue. Thus, tumor samples do not need to be processed prior to using the diagnostic methods described herein.

[0041] In some embodiments, the tumor sample may include stromal material. In some embodiments, the tumor sample may include immune content. In some embodiments, the tumor sample may include a combination of stromal material and immune content.

[0042] In some embodiments, the MBR score stratifies patients within a group of patients having high ICR. In some embodiments, a low MBR score identifies patients within a group of patients having high ICR as having a higher probability of survival.

[0043] Provided is a method for treating colon cancer in a patient, the method comprising administering a microbiome composition comprising Ruminococcus bromii, or other MBR species with a negative weight in an MBR calculation, to the patient. In some embodiments, the composition is administered orally. In some embodiments, the composition is administered intravenously. In some embodiments the composition is administered via other parenteral routes such as intratumoral, peritumoral, intradermal, subcutaneous, intramuscular, or intraperitoneal injection.

[0044] In some embodiments, the microbiome composition further comprises one or more microbes selected from one or more genera depicted as low risk (e.g., negative coefficients, low weight in MBR calculation) in FIG. 5A. In some embodiments, the microbiome composition further comprises one or more microbes selected from genera UBA1819 (Ruminococcaceae), Lachnospiraceae NK4A136 group (Lachnospiraceae), Erysipelotrichaceae UCG-003 (Erysipelotrichaceae), Sellimonas (Lachnospiraceae), Coprococcus 1 (Lachnospiraceae), Terrisporobacter (Peptostreptococcaceae), Olsenella(Atopobiaceae), Oscillibacter (Ruminococcaceae), Dialister (Veillonellaceae), Parabacteroides (Tannerellaceae), Parvimonas (Family XI), [Eubacterium] eligens group (Lachnospiraceae), Bifidobacterium (Bifidobacteriaceae), Hungatella (Lachnospiraceae), Mogibacterium (Family XIII), Caproiciproducens (Ruminococcaceae), Erysipelatoclostridium (Erysipelotrichaceae), Halomonas (Halomonadaceae), or Odoribacter (Marinifilaceae).

[0045] Also provided is a method of identifying a patient as suffering from colon cancer, the method comprising detecting the presence of one or more genes selected from SEC31A, BAX, ZMYM2, BCR, ATF7IP, BCL1 B, ERCC5, or PARK2 in a tumor sample of the patient.

[0046] Further provided is a method of identifying a patient suffering from colon cancer and having a high ICR as having at least 95% probability of overall survival over a period of at least 5 years, the method comprising: b) detecting the abundance of specific microbes in the tumor samples of patients suffering from colon cancer and assigning a microbiome signature risk (MBR) score; and c) identifying patients suffering from colon cancer and having a high ICR and having a low MBR score as having at least 95% probability of overall survival over a period of at least 5 years.

[0047] In further embodiments, the disclosure provides a system for identifying patients suffering from colon cancer and having a high ICR as having at least 95% probability of overall survival over a period of at least 5 years, the method comprising: a) using AC-ICAM246, n=246, as a testing set to develop MBR scores; b) using AC-ICAM42 and TCGA-COAD, as the validation sets to validate the MBR scores of step a); and c) identifying patients suffering from colon cancer and having a high ICR and having a low MBR score as having at least 95% probability of overall survival over a period of at least 5 years.

[0048] In various embodiments, to detect clinically relevant associations between the microbial repertoire and clinical outcome, at the methods include identifying a microbiome signature predictive of survival using genus level data from 16S rRNA gene sequencing. Two datasets for microbiome analysis may be used as described herein in the examples section. 1) AC-ICAM246, n=246, as the testing (or training) set and 2) AC-ICAM42 and TCGA-COAD, as the validation set. In some embodiments, the methods include running a multivariable elastic-net OS Cox regression model in AC-ICAM246, and selecting 41 features (taxa) with acoefficient different to zero (i.e., associated with differential risk of death,). In some embodiments, this list of taxa and associated coefficients is deemed an MBR classifier. In some embodiments, the methods include assigning a score to each sample (the MBR score) by applying the MBR classifier (e.g., as shown in FIG. 5A).

[0049] In some embodiments, the methods include diagnostic or prognostic predictions as follows. A low MBR score (MBR<0, MBR Low), in our training cohort was associated with a considerable (85%) reduction of risk of death. The associate between MBR Low (risk) and prolonged OS was confirmed in two independent testing sets (ICAM42 and TCGA-COAD cohorts), individually and combined.

[0050] The methods described herein provide a quantitative diagnostic test and / or accurate biomarkers for clinical outcomes. In this cohort study, comprehensive genomic analyses on fresh-frozen samples from 348 patients affected by primary colon cancer, encompassing RNA, whole-exome, deep T cell receptor and 16S bacterial rRNA gene sequencing on tumor and matched healthy colon tissue, complemented with tumor whole-genome sequencing for further microbiome characterization was performed. A type 1 helper T cell, cytotoxic, gene expression signature, called Immunologic Constant of Rejection, captured the presence of clonally expanded, tumor-enriched T cell clones and outperformed conventional prognostic molecular biomarkers, such as the consensus molecular subtype and the microsatellite instability classifications. Quantification of genetic immunoediting, defined as a lower number of neoantigens than expected, further refined its prognostic value. A microbiome signature, driven by Ruminococcus bromii, associated with a favorable outcome was identifed. By combining microbiome signature and Immunologic Constant of Rejection, a composite score (mICRoScore), which identifies a group of patients with high survival probability was developed and validated.

[0051] In some embodiments, provided herein in a composite score (mICRoScore) obtained from ICR and MBR score. A ICR high and MBR low score is deemed mICRoScore high. While ICR high has been reported as being associated with improved survival, a composite score of ICR high and MBR low identified a subgroup of patients within ICR high patients who had a 97% 5-year OS, with only three deaths detected at a later follow-up that were not related to colon cancer. No deaths were observed during the entire follow-up period of 5 years in patients with mICRoScore high in the TCGA-COAD cohort (n = 107, testing set). In both the training (AC-ICAM) and the testing (TCGA-COAD) sets, the mICRoScore segregated patientswithin ICR high group who had a very high (e.g., 97%) probability of survival, highlighting the value of using the MBR score in conjunction with the ICR.EXAMPLES

[0052] Materials: Fresh-frozen tumor samples and matched neighboring healthy colon tissues (tumor-normal pairs) from systemic treatment-naive, patients with histological diagnosis of colon carcinoma were profiled with orthogonal genomic platforms. After crossplatform quality control (based on whole-exome sequencing ( WES) and RNA-seq data) and inclusion criteria checking, genomic data from 348 patients were retained and used for downstream analyses (FIG. 1A). The median follow-up time was 4.6 years. This resource is the Sidra-LUMC AC-ICAM: an Atlas and Compass of Immune-Cancer-Microbiome interactions.

[0053] Methods

[0054] Samples used in this observational cohort study (tumor tissue and matched healthy colon tissue, AC-ICAM cohort) are from patients with colon cancer diagnosed at Leiden University Medical Center, the Netherlands, from 2001 to 2015 that did not object for future use of human tissues for scientific research and that were consented on biospecimen protocol ‘Immunology and Genetic of colon Cancer’ approved by the Committee on Medical Ethics of Leiden University Medical Center (study protocol no. POO.193 (06 / 2001)). Snap-frozen tumor and healthy colon tissue were stored at -80°C until processing for DNA and RNA extraction. DNA and RNA from those samples were extracted at Leiden University Medical Center and then transferred to Sidra Medicine for sequencing together with de-identified clinicopathological data of the corresponding patients (Sidra Medicine IRB study protocols no. 1768087-1 (04 / 2016)71602002725 (06 / 2022)). All genomic assays ( WES, WGS, 16S RNA gene sequencing, RNA-seq, TCR sequencing and PCR) were performed at Sidra Medicine.

[0055] Patient information was de-identified and patient samples were anonymized and handled according to the medical guidelines described in the Code of Conduct for Proper Secondary Use of Human Tissue of The Federation of Dutch Medical Scientific Societies. This research was performed according to the recommendations outlined in the Helsinki Declaration.

[0056] Each assay included all samples that had sufficient material (for example, DNA or RNA) available at the time of processing considering the need to preserve aliquots for additional / future assays.

[0057] Collection of biological samples

[0058] Snap-frozen tumor and healthy colon tissue were collected from patients with colon cancer who under went surgical resection of the primary tumor between 2001 and 2015 at Leiden University Medical Center. Patients who received radiotherapy and / or chemotherapy before resection and patients with a primary tumor of non-epithelial origin were excluded. Based on tissue availability, successful nucleic acid extraction and subsequent sequencing quality control (QC), data from 348 patients were retained in the final AC-ICAM cohort. Clinicopathological and follow-up data were retrospectively collected from hospital records. Extensive clinicopathological and survival data of the cohort were available.

[0059] Statistical analysis

[0060] Details of the statistical analysis are described in each method section. All P values were two-sided. Multiple testing corrections were performed by calculating the false discovery rates (FDRs) using the Benjamini-Hochberg method, as appropriate. For missing data, no data imputations were used.

[0061] Survival analysis

[0062] Kaplan-Meier curves were generated using ggsurvplot from R package survminer (v.0.4.9). Hazard Ratios (HRs) between any two groups of interest and corresponding P values based on a Cox proportional hazard regression analysis and 95% confidence intervals (95% Cl), were calculated using R package survival (v.2.41-3). Cox proportional hazard analysis was only computed when both groups of comparison consisted of at least ten patients. Overall P value for comparison of survival between two or more groups was also calculated by log-rank test.

[0063] Multivariate Cox regression was performed using conventional clinical and biological variables. Separate multivariate Cox regression analyses were run including age (continuous), pathological stage (ordinal), microsatellite instability (MSI) status (binary) and consensus molecular subgroups (CMS) (categorical). Additional variables that were found significant in univariate Cox proportional hazard regression analysis were added to these models. These variables included, ICR score (continuous) or ICR cluster (ordinal), genetic immune editing (GIE) (binary) and MBR group (binary). Forest plots were generated using ‘forestplot’ (v.1.7.2).

[0064] Tissue processing

[0065] T umor and healthy tissue samples (unselected for tumor cell purity) were sectioned in a cryostat until the surface area was sufficient to assess tissue morphology by H&E staining. Non-target tissue was removed by macrodissection, including necrotic or adipose tissue and for tumor tissue samples, healthy colon tissue. When macrodissection was required, an H&E- stained slide was examined after this to confirm removal of unwanted tissue types. Frozen tissue was then sectioned at 20 pm until approximately ~10-15 mg was collected per sample. A final section post-sample processing was made for H&E staining. The collected tissue was stored at -80 °C for a few months until DNA and RNA extraction.

[0066] QC metrics of RNA and DNA data were superimposable between samples collected over the years.

[0067] DNA and RNA extraction

[0068] Nucleic acid extraction from fresh-frozen tissue sections was performed using the QIAGEN AllPrep DNA / RNA Mini kit following the manufacturer’s protocol. This process was fully automated on a QIAGEN QIAcube. p-mercaptoethanol (P-ME) was added to the lysis buffer on the day of use. Lysis was performed by completely submerging the sections in 350 pL lysis buffer. Tubes were rotated for at least 1 h at room temperature to allow complete homogenization. QIAcube AllPrep DNA / RNA Mini kit Standard (v.2) program was run, after which DNA and RNA samples were stored at -80 °C. The same DNA was used for human and microbiome sequencing. Samples were shipped from Leiden University Medical Centre (LUMC), The Netherlands to Sidra Medicine, Qatar under a temperature-controlled environment at -80 °C (for 4 d). Samples from 361 patients were sequenced by WES and RNA-seq. Samples from 13 patients were excluded as they did not pass QC, including concordance between healthy and tumor samples.

[0069] The final cohort included 348 patients, for which RNA-seq for tumor samples was possible and passed QC. A subset of samples from these patients were processed with additional assays including WGS, TCR sequencing and 16S RNA gene sequencing, based on the availability of samples for these assays, as described in the following sections.

[0070] RNA sequencing

[0071] The integrity and concentration of the extracted RNA was assessed on the LabChip GXII Touch HT using the RNA Assay and the DNA 5K / RNA / Charge Variant Assay LabChip (PerkinElmer). Sequencing mRNA libraries were constructed from 500 ng of total RNA using the Illumina TruSeqStranded mRNA kit (Illumina). cDNA was synthesized using SuperscriptIV Reverse Transcriptase (Thermo Fisher) and amplified for 15 cycles after ligating with TruSeq RNA Combinatorial Dual-Index adapters. Clonal amplification and cluster generation was performed using Illumina’s cBot 2 System. Sequencing libraries were run on Illumina HiSeq platforms using 75 bp (93% of samples) or 150 bp (7% of samples) paired-end reads at the Clinical Genomics Laboratory, Sidra Medicine. A coverage of 20 M reads per sample was targeted. Obtained coverage was 18.4 M (s.d. 4.7 M).

[0072] Transcriptomic data processing

[0073] Data conversion and demultiplexing was performed using bcl2fastq2 conversion software (v.2.20). FastQC was run to perform QC checks on the raw sequence data (Python v.2.7.1 , FastQC v.0.11.2). Trimming of adaptor sequences was performed using flexbar (v.3.0.3) using Illumina primers FASTA file. Subsequently, reads were aligned to reference genome GRCh38.93 by Hisat2 (v.2.1.0) using SAMtools (v.1.3). After alignment, QC was performed to verify quality of the alignment and paired-end mapping overlap (Bowtie2, v.2.3.4.2). Finally, the feature-Counts function of subreads (v.1 .5.1) was used to count paired reads per genes. Gene expression normalization was performed within lanes, to correct for gene-specific effects (including GC content) and between lanes, to correct for sample-related differences (including sequencing depth) using R package EDASeq (Exploratory Data Analysis and Normalization for RNA-seq) (v.2.12.0). The resulting expression values were quantile normalized using R package preprocessCore (v.1.36.0). All downstream analysis of the expression data was performed using R (v.3.5.1 , or later).

[0074] Whole-exome sequencing

[0075] DNA concentrations were quantified using Quant-iT broad range dsDNA Assay (Thermo Fisher) on the FlexStation 3 Microplate reader (Molecular Devices). DNA of both tumor and matched normal samples was available for 294 patients. Whole-exome libraries were constructed with the Agilent SureSelect XT Target enrichment kit and the exonic DNA was captured using the Agilent SureSelect XT Human All Exon V6r2 capture library for 60-Mb exonic regions. Libraries were constructed using 250 ng of DNA and were sequenced on Illumina’s HiSeq 4000 platform using 150 bp paired-end reads (150PE) at the Genomics Core, Sidra Medicine. Reads were mapped to reference genome hs37d5 (1000 Genomes Phase2 Reference Genome Sequence) based on GRCh37 / hg19 using BWA (v.0.7.12)65.

[0076] WES (200* for tumor and 100* for normal) had an on-target sequencing rate of 65- 70%. The median (across samples) of the average target coverage (per sample) was 129* (interquartile range (IQR) 18) for tumor samples and 69* (IQR 10) for normal samples.

[0077] In tumors, sequencing achieved >20-fold coverage of at least 99% of targeted exons and >70-fold in at least 81 % targeted exons. In healthy samples, sequencing achieved >20- fold coverage of at least 94% of targeted exons and >30-fold in at least 84% targeted exons. Adaptor trimming was performed using the tool trimadap (v.0.1.3). ConPair was run to evaluate concordance and estimate contamination between matched tumor-normal pairs. In eight of the pairs a mismatch was detected and for five pairs, a potential contamination was indicated. HLA typing data were used to validate these results. All potential mismatches and contaminations were excluded, retaining 281 patients for data analysis.

[0078] TCGA data

[0079] RNA sequencing. RNA-seq data (raw counts) from TCGA were downloaded and processed using R package TCGAbiolinks (v.2.18.0). Gene symbols were converted to official HGNC gene symbols and genes without symbol or gene information were excluded. Normalization was performed within lanes, to correct for gene-specific effects (including GC content) and between lanes, to correct for sample-related differences (including sequencing depth) using R package EDASeq (v.2.12.0) and quantile normalized using preprocessCore (v.1.36.0). After normalization, samples were extracted to obtain a single primary tumor tissue (TP) sample per patient. Clinical data were sourced from the TCGA Pan-Cancer Clinical Data Resourcel 1 and survival events OS and progression-free interval (relabeled here as PFS) were used. ICR clustering and calculation of ICR score was performed exactly as described for the AC-ICAM cohort. For the TCGA-COAD cohort, the optimal number of clusters for best segregation based on the Calinski-Harabasz criterion was three. CMS classification of TCGA- COAD samples was performed as described for the AC-ICAM cohort. The Single Sample Predictor by ‘CMSclassifier’ (v.1.0) was used for comparison of CMS classification between AC-ICAM and TCGA-COAD.

[0080] A renormalized matrix of both TCGA-COAD and AC-ICAM datasets was generated by merging the raw counts matrices and performing the EDASeq normalization, as described above, on this combined matrix. These data were used to calculate ssGSEA scores for deconvoluted immune cell subpopulations, immune signatures and oncogenic pathways, to compare between cohorts.

[0081] Somatic mutation data

[0082] Somatic mutation calls from the TCGA MC3 Project were downloaded using R package TCGAmutations (v.0.3.0) using the function tcga_load() with parameters ‘COAD’ for study and ‘MC3’ for source. The downloaded Mutation Annotation Format (MAF) file contained 406 distinct TCGA tumor sample barcodes and 18,183 genes (Hugo Symbol). This file was filtered to only include nonsynonymous mutations (‘Frame_Shift_Del’, ‘Frame_Shift_lns’, ‘ln_Frame_Del’, ‘ln_Frame_lns’, ‘Missense_Mutation’, ‘Nonsense_Mutation’, ‘Splice_Site’, ‘T ranslation_Start_Site’, ‘Nonstop_Mutation’), analogous to the variant filter applied to the AC- ICAM somatic mutation calls.

[0083] Microbiome

[0084] Microbiome genus relative abundance matrix for TCGA-COAD cohort (125 tumor samples and 221 genera, WGS data) was downloaded from The Cancer Microbiome Atlas website. TCGA-COAD relative abundance matrix was filtered to exclude duplicated samples (samples from vial B, eight samples). Overall, 81 genera were present with a nonzero abundance in at least one of the 117 samples (main matrix). The same filter as the one used for AC-ICAM 16S RNA gene-sequencing data (presence in at least 10% of the samples with at least 1 % relative abundance in one sample) was applied, 27 taxa at the genus level were retained.

[0085] NHS and the HPFS study data

[0086] Somatic mutation data. Somatic mutations in NHS and HPFS Colorectal Cancers were downloaded from the supplementary data of the Giannakis et al. study (Giannakis, M. et al. Cell Rep. 15, 857-865 (2016), Supplementary Table 3). The downloaded file contained 619 distinct tumor sample barcodes and 19,208 genes (Hugo Symbol). Samples with tumor anatomic site specified as rectum (anatomic site is available in Giannakis Supplementary Table 1) were excluded, and 482 colon cancer samples were retained. Only nonsynonymous mutations were included at the variant filter (‘Frame_Shift_Del’, ‘Frame_Shift_lns’, ‘ln_Frame_Del’, ‘ln_Frame_lns’, ‘Missense_Mutation’, ‘Nonsense_Mutation’, ‘Splice_Site’, ‘Translation_Start_Site’, ‘Nonstop_ Mutation’), analogous to the variant filter applied to the AC-ICAM and TCGA-COAD somatic mutation files.

[0087] Cancer-related gene annotation

[0088] A cancer-related gene list was constructed from using different sources, as previously described (Saad, M. et al. Lancet Oncol. 23, 341-352 (2022)). (1) genes used by two consortiato define germline genetic variations in pediatric cancers (n = 159;34 n = 565 (Zhang, J. et al. N. Engl. J. Med. 373, 2336-2346 (2015)); (2) genes with at least one pathogenic or likely pathogenic germline variants in the TCGA cohort (n = 99)( Huang, K. et al. Cell 173, 355-370 (2018)); (3) genes classified as driver genes according to the most updated TCGA analysis (n = 299) (Bailey, M. H. et al. Cell 173, 371-385 (2018)); (4) genes included in the MSK- IMPACT (n = 505), MSK-IMPACT HEME (n = 575), Foundation One CDx (n = 324) and Foundation One Heme (n = 593) panels; (5) cancer genes cataloged as tier 1 by the Sanger Cancer Gene Census (n = 576); and (6) cancer genes defined as such by Vogelstein, B. et al. Science 339,1546-1558 (2013). Sources 4-6 were downloaded from OncoKB (Chakravarty, D. et al. JCO Precis. Oncol 1 , 1-16 (2017)). Original sources’ gene names were converted into Ensemble GRCh37 gene symbols. The final list included 1 ,219 unique cancer genes.

[0089] Transcriptome analysis

[0090] ICR score and clustering. Consensus clustering based on 20 a priori selected ICR genes (IFNG, IRF1 , STAT1 , IL12B, TBX21 , CD8A, CD8B, CXCL9, CXCL10, CCL5, GZMB, GNLY, PRF1 , GZMH, GZMA, CD274 / PDL1 , PDCD1 , CTLA4, FOXP3 and IDO1)21 , was applied to the normalized log2-transformed expression matrix using R package ConsensusClusterPlus (v.1.42.0) using 5,000 repeats, agglomerative hierarchical clustering with Ward criterion inner and complete outer linkage. The optimal number of clusters allowing for the best segregation of samples was based on the Calinski-Harabasz criterion. Optimal number of clusters used for segregation was three. Colon cancer samples in the cluster with the highest expression of ICR genes were designated as ‘ICR high’, the intermediate cluster as ‘ICR medium’ and the cluster with the lowest expression was designated ‘ICR low’. The mean Iog2-transformed expression value of the 20 ICR genes is referred to as the ICR score.

[0091] CMS classification

[0092] Samples were classified according to CMS by R package ‘CMSclassifier’ (v.1.0) using random forest method. The obtained CMS labels (from the column ‘RF.predictedCMS’ in output dataframe) were used for all downstream analyses with the exception of the comparison of CMS subtypes between AC-ICAM and TCGA cohort.

[0093] To allow between-cohort comparison, the CMSclassifier was run using the ‘singlesample predictor’ method. This method makes it possible to predict unique samples, with a constant output whether the sample is predicted alone or within a series of samples (Guinney,J. et al. Nat. Med. 21 , 1350-1356 (2015)) and can therefore be used for comparison across cohorts.

[0094] Dimension-reduction of the complete expression matrix was performed using t-SNE by ‘Rtsne’ (v.0.15) and visualized using ggplot2 (v.3.3.2). The t-SNE plot was annotated with distinct colors to visualize the distribution of samples of different CMS (using random forest method) in high-dimensional space. The same t-SNE plot was annotated by ICR cluster in a separate panel. A circos plot to visualize the relation between CMS and ICR classifications was generated using the chord- Diagram function from R package ‘circlize’ (v.0.4.8).

[0095] Immune cell deconvolution and ESTIMATE

[0096] Consensus tumor microenvironment cell estimation (ConsensusTME) was performed to estimate relative abundancies of specific immune cell subsets from bulk transcriptome data. This method relies on integrated gene sets from multiple sources that have been curated and validated on a per-cancer-type basis, using benchmark datasets and seems to outperform previously published methods (Jimenez-Sanchez, A., et al., Methods Cancer Res. 79, 6238- 6246 (2019)). The ConsensusTME was applied using R package ConsensusTME (v.0.0.1 .9) using parameters ‘COAD’ to specify cancer type and ‘ssgsea’ as statistical method.

[0097] The median of each ConsensusTME score was calculated per CMS stratified by ICR cluster and was displayed in a dotted heat map using R package ComplexHeatmap (v.2.1 .2). The association of each ConsensusTME score with OS and PFS was calculated by Cox proportional hazard regression. HR and corresponding 95% Cis as are displayed as forest plots (forestplot v.2.0.1).

[0098] To infer estimated levels of overall stromal and immune cell infiltration to the tumor, the ESTIMATE algorithm (v.1.0.13) was applied to the expression data in R. ESTIMATE was run for both TCGA-COAD dataset and the AC-ICAM cohort. The combined ESTIMATE score for both the stromal and immune signature was compared between cohorts and a box-plot was generated using ggplot2 (v.3.3.2).

[0099] Analysis of tumor-related signatures and immune traits.

[0100] Single-sample gene set enrichment analysis (ssGSEA) was applied to the log2- transformed, normalized gene expression matrix71 (GSVA, v.1.38.2). Gene sets that reflect specific tumor-related pathways were selected from multiple sources as described in detail in Roelands et al. and Supplementary Source Data Table 6a. Enrichment scores of each of these 48 pathways by CMS were visualized using Complex-Heatmap (v.2.1 .2). To better understandthe interactions between tumor-intrinsic signaling and the immune microenvironment, the Pearson correlation between the ICR score and the scores of the 48 tumor-related pathways was calculated. This analysis was performed in the total cohort as well as across CMS subtypes.

[0101] Immune traits considered for analysis were based on a collection of well-characterized immune traits. This collection includes 68 gene signatures related to immunomodulatory signaling, including IFN signaling, TGF-p, wound healing (core serum response) and T cell / B cell response. Gene expression values were median centered and gene symbols were mapped to EntrezIDs (org.Hs.eg.db_3.6.0). Signatures scores were then mean centered and their s.d. values were scaled to one. For all other immune traits, ssGSEA was applied. These included signatures for antigen-presenting machinery (APM1 and APM2) and angiogenesis and nine TCGA-based coexpression signatures (metagene attractors). This collection was supplemented with the tumor inflammation signature and two non-overlapping signatures of IFN-stimulated genes (ISGs), including IFNG hallmark gene set IFNG.GS and ISG resistance signature (ISG. RS), calculated using ssGSEA. Finally, the deconvoluted immune cell abundancies by ConsensusTME70 and ICR score were included among the immune traits. In total 103 immune traits (including ConsensusTME) were used.

[0102] The pairwise Pearson correlation between all immune traits was calculated and the resulting correlation matrix was plotted using ComplexHeatmap (v.2.1.2) with hierarchical clustering. Co-clustering immune traits that formed distinct modules were visualized and labeled according to the immune traits’ enrichment. The clustering was compared to previously defined immune trait modules within a pan-cancer setting, by annotation of the correlation matrix with the previously defined clusters in Sayaman, R. W. et al. Immunity 54, 367-386 (2021).

[0103] Survival analysis on AC-ICAM subsampling.

[0104] AC-ICAM was subsampled hundreds of times in two ways, one was random, the other was on a subgroup of samples with an ESTIMATE distribution that approximates that of the TCGA-COAD. The function ‘approxfun’ in R was used to generate a function to approximate the density of ESTIMATE scores in TCGA-COAD. Cases were sampled from AC-ICAM using the ‘sample’ function in R with prob argument set to sample points with probability distribution of the TCGA-COAD. Each subsampled cohort consisted of 200 samples. The number of subsets in which the Cox proportional regression for ICR score was significant was comparedbetween the two ways of subsampling, statistical significance was determined using a chi- squared test.

[0105] TCR targeted sequencing by immunoSEQ assay

[0106] This sensitive and specific dedicated assay requires high quantity of genomic DNA(>2 pg) and sample selection was exclusively based on DNA availability. TCR sequencing was performed using extracted DNA of 114 primary tissue samples and ten matched healthy colon tissues with sufficient DNA available. DNA samples were normalized to a concentration of 125 ng pL-1 using 3.840 pg of DNA as input per sample. The immunoSEQ assay from Adaptive Biotechnologies was used to amplify all possible variable, diversity and joining (VDJ) gene rearrangements of the TCRp locus (TRB) using a multiplex PCR method. PCR and magnetic bead cleanup were performed according to manufacturer’s instructions. Recommended QC was performed after the first PCR and second PCR amplification steps by running the PCR product on an agarose gel. Purified second PCR amplification products were pooled and the library pool was quantified using Agilent Bioanalyzer 2100. Subsequently, pools were diluted to a concentration of 1 pM and sequenced on Illumina NextSeq 500 / 550 system with Mid Output kit (150 cycles) and Custom NextSeq Sequencing Primer (P / N, M150) (read 1 , 156 cycles and read 2, 9 cycles). Sequencing was performed using survey resolution (two replicates per sample). A sample manifest was created in immunoSEQ Analyzer and the raw sequencing data were uploaded to the Adaptive Biotechnologies cloud following the manufacturer’s instructions. Data were processed using the company’s proprietary pipeline. Number of total templates analyzed per sample ranged 1 ,906-95,834 (median 21 ,258). The average read coverage per sample ranged 11.4-80.6 (median 36.2).

[0107] TCR analysis

[0108] TCR immunoSEQ data analysis. ImmunoSEQ sample-based output variables, as made available by the immunoSEQ Analyzer, include the total number of templates analyzed, number of productive templates, fraction productive templates, number of total rearrangements, number of productive rearrangements, productive clonality and the maximum productive frequency. Herein, the total number of templates reflects the total number of T cells analyzed, of which only the rearrangement in the sample are inframe and do not contain a stop codon).

[0109] The total number of productive rearrangements is the total number of unique T cell clones and clonality is calculated by normalizing the productive entropy using the total numberof productive rearrangements and subtracting the result from 1. Values for (productive) clonality range from 0 to 1 , with values near 0 reflecting more polyclonal samples and values near 1 representing samples with just few predominant rearrangements dominating the observed T cell repertoire (TRB gene).

[0110] A high T cell clonality implies presence of expanded T cell clones. Relationships between ICR score, immune traits, number of productive templates and productive clonality were tested using Pearson’s correlation and visualized by scatter-plots using ggplot2 (v.3.3.2). Similarly, Pearson’s correlation coefficient was calculated between productive clonality and each of the 18,270 genes in the expression matrix. A volcano plot was used to visualize significant results (ggplot2). The top 50 genes with the highest correlation with TCR productive clonality were mapped to the Global Molecular Network and core network analysis was performed using Ingenuity Pathway Analysis software. Data on all productive rearrangements per sample were exported from the immunoSEQ Analyzer Rearrangement Details View. This file includes the exact nucleotide sequence generated through V(D) J recombination, corresponding amino acid sequence, number of templates and productive frequency. Overlapping TCR sequences between tumor samples and matched healthy colon tissues (n = 9) were evaluated and visualized by scatter-plots (ggplot2). Sequences with a productive frequency at least 32-fold higher in the tumor compared to the healthy colon tissue and a tumor productive frequency >0.1 % were defined as tumor-enriched sequences, as previously implemented by Beausang, J. F. et al. Proc. Natl Acad. Sci. USA 114, E10409-E10417 (2017). The fraction of tumor-enriched TCR sequences in the tumor was calculated by dividing the number of productive templates of tumor-enriched sequences by the total number of productive templates per tumor sample. Pearson’s correlation coefficient between the fraction tumor-enriched TCR sequences and ICR score was calculated.

[0111] MiXCR for TCR repertoire derived from bulk RNA-seq.

[0112] The software MiXCR (v.3.0.13)30 was used to retrieve the VDJ repertoire from bulk RNA-seq data aligned to reference genome GRCh37. MiXCR was run through docker and with the single command analyze shotgun. The R package ‘immunarch’ was used to analyze the MiXCR output into the R environment. For the TCRp locus (TRB), the TCR clonality was calculated as 1 - normalized Shannon entropy for all samples, except seven cases for which MiXCR failed to identify clones.

[0113] Whole-exome-sequencing data analysis

[0114] Somatic mutation calling and small insertions and deletions. SNVs were called using mutect (v.1.1.7) and somatic small insertions and deletions (indels) using strelka2 (bcbio- nextgen v.1 .1 .1). An optimized variant filtering pipeline was applied. To filter out false-positive single-nucleotide polymorphism calls, fpfilter was used, the applied filtering parameters are specified in the fpfiler. pl script shared on GitHub. Subsequently, MAF files were generated using VCFtoMAF tool (v.1.6.16), which also appended the SIFT (sorting intolerant from tolerant), PolyPhen and Exome Aggregation Consortium annotations. MAF files were loaded into R where indels with low complexity regions were excluded. For both SNVs and indels, a cutoff for minimum allele fraction of 5% and tumor depth of more than three reads was applied. The Exome Aggregation Consortium data were then used to filter out common variants that are encountered in >1 % in the general population. After these technical exclusion criteria, biological filters were applied, including selection of nonsynonymous mutations (frame shift deletions, frame shift insertions, inframe deletions, inframe insertions, missense mutations, nonsense mutations, nonstop mutations, splice site and translation start site mutations).

[0115] The resulting number of variants / mutations per Mb (capture size is 40 Mb) per sample is referred to as the nonsynonymous tumor mutational burden (TMB). Next, to identify most frequently mutated genes in our cohort that might play a role in cancer, variants that are predicted to be tolerated according to SIFT annotation or benign according to PolyPhen (polymorphism phenotyping) were excluded. Finally, all artifact genes, which are typically encountered as bystander mutations in cancer that are mutated for example as a consequence of a high homology of sequences in the gene, were excluded. The OncoPlot function from ComplexHeatmap (v.2.1.2) was used to visualize the most frequent somatic mutations.

[0116] Comparison of TMB with TCGA datasets.

[0117] To compare the TMB in the AC-ICAM with all 33 TCGA cohorts derived from the MC3 project, The tcgaCompare function from maftools (v.2.6.05, R) was used. For AC-ICAM, the filtered MAF for nonsynonymous mutations was used as input with specified capture size of 40.

[0118] Comparison of somatic mutations with other cohorts.

[0119] To define mutated genes in the AC-ICAM that were not previously described in colon cancer, a comparison of the most frequently mutated genes in AC-ICAM (>5% of the tumor samples) with frequencies detected in previously published datasets containing colon cancersamples (TCGA-COAD and NHS-HPFS) as well as reported cancer driver genes (Bailey, M. H. et al. Cell 173, 371-385 (2018)) or colon oncogenic mediators (Colaprico, A. et al. Nat. Commun. 11 , 69 (2020)) was performed. First, extracted genes with a nonsynonymous mutation frequency >5% in the AC-ICAM cohort. Subsequently, only genes that are likely involved in cancer development, were retained. All artifact genes (mutations typically encountered as bystander mutations in cancer that are mutated for example as a consequence of a high homology of sequences in the gene), were excluded.

[0120] Genes that have previously been reported as colon cancer oncogenic mediator or cancer driver gene for colorectal cancer (COADREAD) were also excluded. Finally, only genes with a mutation frequency <5% in the NHS-HPFS colon cancer cohort (Giannakis, M. et al. Cell Rep. 15, 857-865 (2016)) and <5% in TCGA-COAD were maintained. As a final filter, only genes that had a nonsynonymous mutation frequency of at least twofold in AC- ICAM compared to TCGA-COAD were labeled as potentially new in colon cancer.

[0121] Estimation of MSI from whole-exome sequencing data.

[0122] MANTIS (v.1 .0.4), a tool for rapid detection of microsatellite instability was applied on our WES data (Bonneville, R. et al. JCO Precis. Oncol. 1 , 1-15 (2017)). Briefly, a bed file suitable for use by MANTIS was created using RepeatFinder function of the MANTIS tool, to find microsatellites regions within the reference genome (GRCh37). MANTIS was then run for each tumor and matched normal BAM file pair using these detected microsatellite loci. The instability score between the two samples within the pair was used to classify samples either as MSI-H (MANTIS score > 0.4) or MSS (MANTIS score < 0.4).

[0123] Somatic mutations associated with ICR. The association of specific somatic alterations was investigated, including SNVs and small insertions or deletions (indels) and ICR immune phenotype. Binomial linear regression models were fitted to define which specific mutations associate with ICR score using the glm function with family ‘binomial’ (R). This analysis was performed in the total cohort (n = 281) as well as within hypermutated (n = 69) and nonhypermutated (n = 212) subgroups separately. The estimate and P value were extracted for each gene and FDR was calculated using the Benjamini-Hochberg method. Significant genes with an FDR < 0.1 and that were mutated in at least five patients in the analysis subgroup were plotted as OncoPrints (ComplexHeatmap, v.2.1.2).

[0124] Mutations in homologous recombination genes, mucinous histology, and ICR.

[0125] Genes with an inverse association with ICR score within hypermutated colon cancer included genes involved in homologous recombination repair. The frequency of mutations in either of the identified genes (BRCA2, BRCA1 and FANCA) genes were compared between hypermutated cases of mucinous histology with hypermutated cases with other histological classifications. An unpaired Student’s t-test was used to compare ICR score between hypermutated cases of mucinous histology with hypermutated cases with other histological classifications.

[0126] Somatic copy-number alteration segmentation

[0127] A segmentation file was generated for each sample and later a merged file for all samples was uploaded to IGV (v.2.11 .0). A pipeline using GATK (gatk-package-4. beta.6) to generate each tumor sample’s segmentation file was used. The following steps were performed: 1. Calculated the coverage of tumor and normal BAM files for each interval using GATK CalculateTargetCoverage. 2. Generated the panel of ‘normal’ using normal samples by GATK CreatePanelOfNormals options. 3. Normalizing the tumor data using GATK NormalizeSomaticReadCounts methods using PON generated during the step above. 4. Performed the segmentation of tumor data using input files from the above steps using GATK PerformSegmentation. 5. The merged segmentation file of all the samples was uploaded to IGV and snapshots were generated.

[0128] Overview of SCNAs.

[0129] The prevalence of Somatic Copy Number Alterations (SCNAs) among ICR clusters and hypermutated and non-hypermutated subgroups was examined by exploration of the segmentation file in IGV. Briefly, the log2-transformed segmentation file was loaded in IGV with reference genome GRCh37, including an annotation text file including mutational load category (hypermutated, non-hypermutated), POLE mutation status, ICR cluster, CMS and MSI status. The samples were ordered consecutively by MSI status, CMS, ICR, POLE and mutational load category. Prevalence of amplification and deletions was visually inspected and compared between groups.

[0130] Genetic immunoediting and immunoediting score HLA typing, neoantigen prediction and GIE.

[0131] HLA typing was performed on both WES and RNA-seq data using OptiType (bcbio- nextgen v.1.1.5 in Python v.2.7.0)78. Neoantigen prediction tool pVACseq from pVACtools was run using the following predictors: MHCnuggetsI, NNalign, NetMHC, SMM, SMMPMBECand SMMalign. The obtained vcfs from our somatic mutation calling pipeline were used as input for pVACseq, along with the predicted HLA type from WES data. Gene expression data aligned to GRCh37 in transcripts per million was annotated to the vcfs using vcf-expression- annotator. Mutant-specific binders, relevant to the restricted HLA-I allele, are referred to as neoantigens, as described in detail by Zhang et al.79. Mutated epitopes with a median IC50 binding affinity across all prediction algorithms used <500 nM, with a corresponding wild-type epitope with a median IC50 binding affinity > 500 nM, were used as criteria to infer neoantigens.

[0132] Predicted neoantigens were used to calculate the GIE value. The GIE value was calculated by taking the ratio between the number of observed versus the number of expected neoantigens. The expected number of neoantigens was based on the assumption of a linearity between TMB and the number of neoantigens. It was therefore assumed that samples that have a lower frequency of neoantigens than expected (lower GIE values), display evidence of immunoediting. A higher frequency of neoantigens than expected indicates a lack of immunoediting.

[0133] IES classification and analysis.

[0134] The IES is a composite score based on both ICR and GIE. Tumors of IES4 are those predicted to be the most immune active, as they are ICR high and display GIE. Tumors of IES1 are expected to be most immune silent, classified both as ICR low and an absence of GIE. Tumors of the intermediate groups IES2 and IES3 reflect ICR-low and GIE and ICR-high and non-GIE tumors, respectively.

[0135] Mutational load category, MSI status and pathological stage distribution was compared between IES groups using a chi-squared test. The OS was compared between patients with different IES and between GIE and non-GIE tumors in the ICR-medium group using Cox proportional hazard regression analysis. A Cox proportional hazard’s multivariate model was fitted with IES (ordinal) and pathological stage (ordinal).

[0136] Association between IES and TCR clonality.

[0137] The Spearman correlation between IES as ordinal variable and TCR clonality from immunoSEQ as well as MiXCR-based clonality was calculated. Several additional analyses were performed to assess whetherthe relationship between TCR productive clonality and IES was driven by ICR. Multiple regression analysis was performed with ICR score and immunoSEQ TCR clonality as continuous variables to predict productive TCR clonality(immunoSEQ). Second, the data were modeled through local polynomial regression fitting of the productive TCR clonality (immunoSEQ) by IES category (ordinal variable).

[0138] Microbiome: bacterial 16S rRNA PCR sequencing

[0139] The 16S rRNA gene sequencing was performed at the Host-microbe Interaction Laboratory, Sidra Medicine. Hypervariable regions V3-V4 of 16S rRNA gene were amplified with PCR using the amplicon primers with Illumina adaptors (underlined):

[0140] Forward:

[0141] 5'TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG'3 (SEQ ID NO: 1)

[0142] Reverse:

[0143] 5'GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAA TCC'3 (SEQ ID NO: 2)

[0144] In brief, PCR was performed in a 25-pL reaction mixture containing 5 pL each forward and reverse primer (1 pM), 2.5 pL template DNA for the samples and 12.5 pL 1 * Hot Master Mix (Phusion Hot Start Master Mix). No human DNA depletion was used. The amplifications were performed on a Veriti 96- well Thermal Cycler (Thermo Scientific) with the following program: initial denaturation at 95 °C for 2 min, followed by 30 cycles of denaturation at 95 °C for 30 s, primer annealing at 60 °C for 30 s and extension at 72 °C for 30 s, with a final elongation at 72 °C for 5 min. The presence of PCR products was confirmed by electrophoresis in a 1.5% agarose gel conducted at 80 V / cm in Tris-borate-EDTA (TBE) buffer. Amplicons were then purified using AgenCourt AMPure XP magnetic beads (Beckman Coulter) according to the Illumina MiSeq

[0145] 16S Metagenomic Sequencing Library Preparation protocol.

[0146] As positive controls, included DNA from stool samples (extracted with QIAGEN QIAmp Fast DNA Stool Mini kit), using the same input of DNA as the one used for the AC-ICAM samples. Similar 16S rRNA amplicon PCR products across the tissue samples and the positive controls were obtained, indicating that the DNA extraction protocol used resulted in enough recovery of the microbial DNA from our specimens.

[0147] Samples were multiplexed using a dual-index approach with the Nextera XT Index kit (Illumina) according to the manufacturer’s instructions. The concentration of amplicons was determined using the Qubit HS dsDNA assay kit (Life Technologies,) followed by pooling toachieve an equimolar library concentration. The final pooled product was paired-end sequenced at 2 x 300 bp using a MiSeq Reagent kit v3 on Illumina MiSeq platform (Illumina) at the Sidra Medicine research facility. Sequencing was also performed on 27 empty wells across plates to exclude the occurrence of large-scale cross-contamination among samples during sequencing procedures: the minimum and maximum read counts were 2 and 234, respectively and the average and median reads counts were 37 and 18, respectively. No negative controls for sampling or DNA extractions were included. Samples were aliquoted randomly in the plate.

[0148] Microbiome: 16S rRNA gene sequencing and data processing

[0149] Sequenced data were demultiplexed using MiSeq Control Software. The overall quality of sequencing quality was evaluated using FastQC and the demultiplexed sequencing data were imported into Quantitative Insights into Microbial Ecology (QIIME2; v.2019.4.0) software package. The data were denoised with DADA2, which includes a multi-step process, including read filtering, dereplication and chimera removal. Paired 250-bp reads were trimmed of the initial five low-quality bases and further processed to generate the amplicon sequence variant, interchangeably called operational taxonomic units (OTUs). The data were subsampled at a depth of 22,704 and then normalized using the rarefaction on OTUs count at even depth. Taxonomic classification was performed utilizing 16S rRNA gene database from Silva classifier (silva-132-99-515-806-nb-classifier). The data were imported into R in a Biological Observation Matrix (biom) format, before further evaluation with Phyloseq (v.1 .34.0). The 16S rRNA gene sequencing was performed on all samples with sufficient DNA available: 246 tumor samples and 246 matched healthy colon tissues from the same patients (AC-ICAM246) and on additional 42 tumor samples (ICAM42) for which there was no sufficient DNA available from the healthy colon counterpart.

[0150] The minimum and maximum read counts were 25,868 and 351 ,069, respectively. The average and median reads counts were 82,506 and 75,668, respectively. No samples were excluded from the analysis. Alpha diversity (within sample community) was assessed by observed OTUs (sum of unique OTUs per sample), Chaol (Chao 1987) an abundance-based richness estimator that is sensitive to rare OTUs, Shannon (Shannon 1948) and inverse Simpson (InvSimpson) (Simpson 1949), the last one being more dependent on highly abundant OTUs and less sensitive to rare OTUs. Indices were read into R using R package vegan (v.2.5-6).

[0151] Relative abundance of distinct microbiome elements was determined using the transform_sample_counts function from Phyloseq, such that sum of all abundance values per sample is equal to one (Microbiome_Relative = transform_sample_counts (pyloseq_object, function(x) x I sum(x))). OTU tables were aggregated by taxonomic ranks including phylum (26 unique phyla), class (48 classes), order (97 orders), family (207 families), genus (562 genera) and species (846 species). As the confidence for annotation of reads decreases with decreasing rank, some reads were only annotated with higher ranks.

[0152] Microbiome: WGS and data processing

[0153] Library construction and sequencing was performed at the Sidra Clinical Genomics Laboratory Sequencing Facility. DNA was quantified using the Quant-iT dsDNA Assay (Invitrogen) on the FlexStation 3 (Molecular Devices). The library was constructed from 250 ng of DNA with the Illumina TruSeq DNA Nano kit. Library quality and concentration was assessed using the DNA 1 k assay on a PerkinElmer GX2 and qPCR using the KAPA Library quantification kit on a Roche LightCycler 480 II. Genomic libraries were sequenced with paired-end 150 bp on HiSeq X (32% of samples) and Novaseq 6000 (68% of samples) systems (Illumina) following the manufacturer’s recommended protocol to achieve a minimum average coverage 60* for tumor samples. Quality passed reads were aligned to the human reference genome GRCh38 using BWA. Human sequencing reads were removed and unaligned nonhost reads were extracted using SAMtools. Low-quality unaligned reads were trimmed and samples were processed for taxonomic profiling using MetaPhlAn2 (ref. 80). MetaPhlAn2 uses a library of unique clade-specific marker genes to estimate bacterial relative abundance at the species level. The program was run with default parameters except analysis type set to relative abundance and restricted to bacterial organisms only. WGS was targeted to achieve >60* coverage per sample. The median (across samples) of the average target coverage (per sample) was 76x (range of 50-92). Of 3.2 x 1011 total reads (median 1.9 x 109 reads per sample; IQR 2.1 X 108), 1.5 x 108 (median 1 x105 reads per sample; IQR 3.4 x 105) were aligned to bacteria. A total of 132 taxa, at genus level were detected, of which 3 were excluded as possible contaminants (Deinococcus, Ralstonia and Enhydrobacter \2 (main matrix). When the same filter as the one used for 16S RNA gene-sequencing data (presence in at least 10% of the samples with at least 1% relative abundance in one sample) was applied, 54 taxa at the genus level were retained. WGS was performed in all samples with sufficient DNA available (n = 167).

[0154] Ruminococcus bromii PCR

[0155] PCR was performed based on (Wang, R.-F.et al., Mol. Cell. Probes 11 , 259-265 (1997)) using R. bromii 16S rDNA forward primer (GAAGTAGAGATACATTAGGTG (SEQ ID NO: 3)) and R. bromii 16S rDNA reverse primer (ACGAGGTTGGACTACTGA (SEQ ID NO: 4)). PCR was performed using AmpliTaq Gold 360 Master Mix (Thermo Fisher, 4398881), 20 ng of sample DNA and 5 nM of each primer. The amplification conditions were one cycle of 95 °C for 10 min, then 35 cycles of 95 °C for 30 s, 50 °C for 30 s and 72 °C for 30 s and finally one cycle of 72 °C for 7 min before storing at 4 °C. PCR products (10 pL each) were separated by electrophoresis in 2% agarose gels (Sigma, A4718) containing ethidium bromide (1 pg mL-1) (Sigma, E1510) using a 100-bp DNA ladder (New England Biolabs, N0551 G) for size verification. PCR band intensity was defined as negative when intensity was absent or extremely faint. PCR was considered positive if band was gradually more intense (graded from 2 to 4). PCR was performed in all samples from the AC-ICAM246 cohort with sufficient amounts of DNA available (n = 126).

[0156] Microbiome data analysis

[0157] Genus-level filtering. On tumor samples, microbiome genera were filtered to include genera which are present in at least 10% of the samples with at least 1% relative abundance in one sample; 138 out of 562 were retained. These included 137 genera and the genus labeled ‘unknown’ that reflects all reads for taxa with insufficient confidence at the genus level. The same filtering was applied to normal samples; 129 genera were retained. A total of 120 genera overlapped between normal and tumor samples, 9 genera were unique in normal samples and 18 genera were unique in tumor samples.

[0158] This set of filtered genera were used for all downstream analysis except for the comparison between tumor and normal pairs. For this analysis included any genera that passed the filtering approach described above for either normal or tumor groups (if taxa passed the filtered in tumor samples they were retained in normal samples and vice versa; total 147 genera).

[0159] Contaminant assessment.

[0160] To remove putative contaminants from the 16S rRNA gene-sequencing data, used a list of microbial taxa that are typically found in negative blank reagents, as described by Salter, S. J. et al. BMC Biol. 12, 87 (2014). This list has previously been curated and annotated by Poore, G. D. et al. Nature 579, 567-574 (2020). by manual review of the literature. Thiscuration allowed the discrimination of taxa that are ‘likely contaminants’, ‘potentially pathological or commensal genera’ and ‘mixed evidence’ genera that have been described both as pathogens as well as contaminants. Those taxa that were ‘likely contaminants’ as well as ‘mixed evidence’ were flagged for potential exclusion from our 16S rRNA gene-sequencing microbiome abundance matrix.

[0161] In total, detected 25 taxa that were ‘likely contaminants’ and 10 taxa with ‘mixed evidence’ in at least one out of the 492 samples. To evaluate the extent of potential contamination by these 35 taxa, calculated the sum of these taxa for each sample. On average, only 0.04732% of the total microbial abundance per sample consisted of ‘flagged’ taxa (min, 0%; first quartile, 0%; median, 0%; third quartile, 0.03485%; and max, 4.46%). Furthermore, most of these putative contaminant taxa (n = 33) were detected in only fewer than 20 (out of 492) samples. Potential contaminating bacteria that were detected in the highest numbers of samples were Oxalobacter in 39 samples and Micrococcus in 28 samples. Detected putative contaminants and taxa with mixed evidence from the 16S rRNA-sequencing data were removed and applied the minimal abundance filter (presence in at least 10% of the samples with at least 1% relative abundance in one sample).

[0162] Microbiome comparison between tumor and healthy colon tissue.

[0163] At the phylum level, the overall distribution of microbiome composition was visualized using stacked bar charts. The order of samples was determined by descending relative abundance of the phylum Fusobacteria in tumor samples and the matching healthy colon samples from corresponding patients were ranked in the same order as the tumor stacked bar chart.

[0164] A paired Mann-Whitney U-test (two-sided) was used to determine microbial phyla / genera with significantly different relative abundance between tumor and paired normal samples. FDR was calculated using the Benjamini-Hochberg method. Results were visualized in volcano plots.

[0165] Microbiome comparison between ICR groups.

[0166] An unpaired Mann- Whitney U-test (two-sided) was used to calculate which filtered genera (n = 138) were differentially abundant between ICR-high and ICR-low samples. FDR was calculated using the Benjamini-Hochberg method. Results were visualized in volcano plots.

[0167] Co-abundance network inference.

[0168] Co-abundance analysis was performed in tumor samples from the AC-ICAM246 cohort. Co-abundance analysis, which involves studying the presence of multiple components within a composition, can be difficult to perform accurately when using relative abundance. This is because the relative abundance of the different components is constrained to sum to 1 , which can lead to the appearance of false correlations. To address this issue, techniques such as co-abundance network inference can be used to more accurately infer relationships between the components. Before co-occurrence analysis, the genus labeled ‘unknown’ was excluded. SparCC was used to calculate the co-occurrences between the 137 remaining taxa using centered log-ratio (clr)-transformed OUT counts in tumor samples (Python, SparCC3). A total of 500 inference and 10 exclusion iterations were used to estimate the median correlation of each pairwise. The statistical significance of the correlations was calculated using a bootstrapping procedure to generate 500 simulated data. For each component pair, pseudo P values (two-sided) were assigned as the proportion of simulated bootstrapped data with a correlation at least as extreme as the one computed for the original data.

[0169] Benjamini-Hochberg FDR was used for multiple testing correction. All the correlations were then sorted using a statistically significant cutoff (FDR < 0.05) and SparCC correlation coefficient > ±0.3. Clusters among the networks (groups of at least three correlated genera using the cutoffs specified above) were defined via a fast greedy clustering algorithm. All cooccurrence networks were made using the R package ‘NetCoMI (v.1.1.0) - Network Construction and Comparison for Microbiome Data’84 and visualized using Cytoscape (v.3.9.1). Within each cluster, the total relative abundance was calculated by summing up the relative abundance values for genera that positively correlated with each other. For each of the identified clusters, survival analysis was performed by binarizing each sample into high and low abundance based on the median total relative abundance of each cluster.

[0170] MBR model development, training set.

[0171] First normalized the genus abundance matrix by converting each genus column into a z score using mean and s.d. and treating the normalized abundance matrix as the training set. Built a relaxed multivariable elastic-net OS Cox regression model using the glmnet R package (v.4.1.4) on the training set. The optimal hyper parameters (y and A) for the best model were identified through fivefold cross-validation via a grid-search technique using the ‘cv.glmnet’ function. Used the concordance index as a performance metric. The parameters for which the mean cross-validation concordance index was the highest were selected asoptimal hyper parameters. Next, the final model was built using these hyper parameters on the complete training set. To calculate risk scores in the training dataset (MBR scores), passed the training set and best model to the ‘predict’ function. A total of 41 features (genera) were present in the best model with nonzero coefficients; these features are referred to as the ‘MBR classifier’, which represents the final model.

[0172] A positive or negative coefficient of each genus of the MBR classifier can be binarized into ‘high-risk’ and ‘low-risk’ groups using the cutoff threshold of 0 and attributed to the strength of association with survival. A higher positive coefficient means high hazard risk of death, whereas a negative coefficient corresponds to lower risk of death.

[0173] MBR model validation, testing sets.

[0174] Validated the final model on two datasets. Both datasets consist of samples that were not used for model training (unseen data). One is an independent internal (ICAM42) dataset, referred to as testing cohort 1 and the other is an external cohort (TCGA-COAD), referred to as testing cohort 2. The ICAM42 consists of 42 samples and TCGA-COAD consists of 117 samples. Processed the two datasets to convert the abundance values for each genus into z scores using the mean and s.d. derived from the training set. These abundance matrices were passed to the ‘predict’ function along with the best model to estimate corresponding risk scores. The risk score (MBR score) of any tested sample is only dependent on the relative abundance of the list of genera that overlap with the ones included in the MBR classifier (the risk score for each sample is not dependent on one of the other samples). Finally, the MBR scores are binarized using the cutoff threshold 0 to categorize the test sample into ‘high-risk’ (>0) and ‘low-risk’ (<0) groups as performed in the training set. Therefore, no cutoff optimization occurred in the validation phase.

[0175] MBR model performance assessment.

[0176] Tested the concordance index (1) in the training set using the final MBR model; (2) in the training set using the cross-validation of the best MBR model (five permutations, 80% training and 20% validation partition); and (3) in each test set cohort separately (ICAM42 and TCGA-COAD) and in the full test set (ICAM42 and TCGA-COAD combined) using the final MBR model.

[0177] Taxa used for the MBR score calculation in other cohorts.

[0178] To calculate the MBR score in each additional dataset, used taxa that overlapped with the 41 genera of the MBR classifier, which was developed using 16S rRNA gene sequencingon tumor samples. There were 16 of 41 taxa in the TCGA-COAD (WGS data) and 18 of 41 taxa in the AC-ICAM WGS data (tumor sample) main matrices. All the 41 taxa were available in the ICAM42 cohort (tumor samples) and the MBR score for AC-ICAM healthy colon tissue samples was based on 36 genera that passed the applied genus-level filtering for healthy tissue (the list of the taxa used for each platform is available in Supplementary Table 11). The Silva classifier used for genus attribution in the 16S rRNA gene-sequencing data includes ‘Ruminococcus 1’ and ‘Ruminococcus 2’, whereas WGS- WES TCGA data only include ‘Ruminococcus’ as genus-level taxa. Therefore, for matching purposes, when calculating the risk score, replaced the labeling of ‘Ruminococcus 1’ and ‘Ruminococcus 2’ with ‘Ruminococcus’. In WGS AC-ICAM ‘R. bromii’ was used instead.

[0179] R. bromii validation analysis.

[0180] The specific species underlying the reads supporting the Ruminococcus 2 taxon from 16S sequencing data was characterized. Previously, a high degree of sequence similarity was reported between the Ruminococcus 2 taxa from the Silva classifier and the species R. bromii. The subset of samples that had both 16S sequencing and WGS data available was used to calculate the Spearman correlation between each Ruminococcus species (from WGS data) and the 16S Ruminococcus 2 (16S) relative abundance. In addition, the proportion of WGS reads that mapped to each specific Ruminococcus species was calculated as fraction of all WGS reads that were assigned to the Ruminococcus genus.

[0181] To confirm the presence of R. bromii, performed a PCR specific to R. bromii on the 126 AC-ICAM tumor samples for which sufficient DNA was still available (see section R. bromii PCR for technical details on PCR). The concordance between detection of R. bromii in PCR and 16S Ruminococcus 2 was defined as the percentage of samples for which both methods had identical results. The discordant cases were further examined by evaluation of WGS results. Furthermore, the concordance between detection of R. bromii in PCR and R. bromii in WGS was assessed in the 86 samples for which data from both methods were available.

[0182] mICRoScore development.

[0183] In view of the individual contribution of analytes extrapolated by individual platforms such as the ICR (RNA-seq data), the GIE ( WES data) and the MBR scores (16S data) and TCR clonality (immunoSEQ and MiXCR), development of a multi-omics parameter was sought that could capture a subgroup of patients with exceptional survival. Each parameter that was significant in the univariate Cox regression analysis (ICR, as ordinal variable, low,medium, high; GIE as binary variable, non-GIE versus GIE; and MBR score, as binary variable, low versus high), was assessed within a multivariable Cox regression model adjusted for age (as continuous variable), CMS subtypes (as categorical variable, CMS1-CMS3, versus CMS4), stage (as ordinal variable, I, II, III and IV) and MSI status (as binary variable, MSS versus MSI-H). The parameters that were retained by the multivariable Cox models were combined into an integrated score. For univariate analysis, used RNA-seq, WES and TCR clonality data from the entire AC-ICAM cohort and MBR score derived from 16S rRNA genesequencing data of the AC-ICAM246 cohort. The mICRoScore reflects the co-presence of ICR high and MBR low risk, defined as mICRoScore high. On the AC-ICAM246 (training set), all samples with MBR-high risk and / or in ICR-medium or ICR-low group are defined as mICRoScore low. The survival between patients with mICRoScore high and mICRoScore low was compared using Cox proportional hazard regression analysis and a log-rank test. mICRoScore validation. Used data from TCGA-COAD as external validation cohort to test the mICRoScore (testing set). The TCGA-COAD cohort includes 107 patients with both tumor microbiome data (downloaded from Dohlman et al. (Cell Host Microbe 29, 281-298 (2021)) and RNA-seq data available (used for ICR estimation). ICR assignments from this cohort (see section TCGA data) were combined with the MBR classification to classify patients as mICRoScore high and mICRoScore low. The survival between patients with mICRoScore high and mICRoScore low was compared using a log-rank test.

[0184] Calculations

[0185] TCR clonality calculation by immunoSEQ assay data (targeted DNA). Entropy (H) is calculated by a standard Shannon entropy calculation with log base 2. Clonality is the inverse of the normalized entropy calculation. The equations are below:

[0186] Shannon entropy : H (x) = -ZP (x) Iog2 [P (x)]

[0187] Specifically, for our data:

[0188] For a productive (inframe) sequence x,

[0189] P (x) = sequence count / total productive count

[0190] Entropy =

[0191] -1 x the sum over all unique productive (inframe) sequences of ( (sequence count / total productive count)

[0192] x|og2(sequence count / total productive count) )

[0193] Normalized entropy =

[0194] entropy / log2(productive unique inframe sequences)

[0195] Clonality = 1 - normalized entropy

[0196] TCR clonality calculation on bulk RNA-seq data (MiXCR).

[0197] Entropy (H) is calculated by a standard Shannon entropy calculation with log

[0198] base 2. The equations are below:

[0199] Shannon entropy H (x) = -ZP (x) Iog2 (P (x))

[0200] For a sequence x,

[0201] P (x) = sequence count / total count

[0202] The Shannon entropy was normalized so that it can assume a value between 0 and 1 . The normalized Shannon entropy is referred to as Pielou’s evenness and is calculated as below:

[0203] Pielou’s evenness : J = H / log (S)

[0204] where S is the number of unique TCR / CDR3 sequences.

[0205] Clonality is calculated as the inverse of the normalized entropy

[0206] (J) calculation:

[0207] Clonality = 1 -J

[0208] Genetic immunoediting value. The GIE value is calculated by taking the ratio between the observed (O) versus the expected (E) number of neoantigens:

[0209] GIE value = O / E

[0210] in which E is a function of the number of nonsynonymous mutations in that specific sample (x):

[0211] E (x) = -2.38770 + 0.09171 xx

[0212] Results

[0213] The ICR outperforms conventional molecular classifications

[0214] As a first objective, conducted a validation of the ICR signature on the AC-ICAM cohort. This objective was predefined before data were generated (prospective validation ofretrospectively collected samples; Methods provides detail). A consensus-clustering approach based on the ICR genes segregated the cohort in three clusters / immune subtypes: ICR high (hot tumors), ICR medium and ICR low (cold tumors). Systematic transcriptomic analysis using 103 previously defined immune traits (Methods) revealed co-clustering of these traits into seven different modules (M1-M7), with ICR belonging to M2 (lymphocyte infiltration signature), together with other immune signatures, including the tumor inflammation signature. Then characterized the immune disposition in relation to Consensus Molecular Subtypes (CMS), a well-defined transcriptomic-based classification of colon cancers. CMS categories include CMS1 / immune, CMS2 / canonical, CMS3 / metabolic and CMS4 / mesenchymal. Overall, t-distributed stochastic neighbor embedding (t-SNE) plotting of the whole expression data segregated CMS1-CMS3 samples, but a high heterogeneity was observed for CMS4. Within CMS subtypes, ICR varied considerably. While most of the CMS1 samples were ICR high, implying immune activation, CMS4 samples were spread across the three ICR immune subtypes. According to the anatomic location, a progressive right-to-left colon enrichment (for CMS2) and depletion (for CMS1) was evident. ICR score (average of the 20 ICR genes) and leukocyte subsets enrichment scores, sho wed only a modest decrease from right-to-left colon, with ICR high being more prevalent in cecum versus rectosigmoid tumors. The enrichment scores of cancer-cell-related pathways were clearly distinct across CMS subtypes. ICR score correlated negatively with certain cancer-cell pathways in all CMS subtypes (for example, WNT-p catenin and NOTCH signaling), whereas a positive correlation with immunosuppressive and stromal-related pathways (for example, transforming growth factor (TGF)-p, epithelial to mesenchymal transition and vascular endothelial growth factor signaling) was only observed in CMS4 tumors.

[0215] The abundance of natural killer (NK) cell and T cell subsets was the highest in the ICR-high immune subtype across all CMS, whereas other leukocyte subsets were more variable. Conversely, the abundance of fibroblast and endothelial cells was increased in CMS4, irrespective of ICR assignment, confirming the increased stromal content in these tumors. Based on statistical significance, the association between ICR score and progression- free survival (PFS) was stronger than what observed for any stromal cell or leukocyte subsets; similar results were obtained for the association with overall survival (OS) (FIG. 1 B-1 , forest plot, FIG. 1 B-2 provides P values for forest plot). ICR immune subtypes had distinct OS and PFS, which gradually increased from ICR low to high. As expected, CMS4 was associated with poor survival; ho wever, ICR reverted this negative trend in survival, with ICR high being associated with longer OS even within the CMS4 group. Conversely, CMS did not stratify theICR-high cluster. ICR remained significantly associated with improved OS in the Cox multivariate analysis (together with pathological stage and age), whereas microsatellite instability (MSI) status and CMS did not. The relationships between ICR and CMS were confirmed in the TCGA colon cancer cohort (TCGA-COAD). Overall, in TCGA, the survival differences were attenuated (in the PFS analysis) or absent (in the OS analysis) for ICR, immune infiltrates and CMS. Nevertheless, ICR still stratified survival in patients with CMS4 cancers. Overall, validated the prognostic role of ICR in colon cancer.

[0216] ICR captures tumor-enriched, clonally expanded T cells

[0217] A dedicated deep sequencing of the TRB gene by immunoSEQ was performed on all samples (114 tumors and 9 healthy colon tissues) with sufficient DNA for this assay. TRB gene sequence information was also extracted from bulk RNA-seq using MiXCR (n = 341)30. Among stromal cell and leukocyte subsets (measured by RNA-seq), the strongest correlation with the number of conventional (ap) T cells with a productive TCR (immunoSEQ TCR productive DNA templates), was observed for estimates of T cell subsets (FIG. 2A), implying robustness of DNA and RNA-based measurements; ho wever, the strongest correlation with immunoSEQ TCR productive clonality was observed for ICR score (r = 0.61), substantiating the ability of ICR to capture additional features beyond T cell abundance. Despite the inherent limitation in terms of sensitivity and specificity of TCR repertoire analysis using bulk RNA-seq, MiXCR TCR clonality correlated well with immunoSEQ TCR clonality (r = 0.64) as well as with ICR (r = 0.40). Consistently, among ICR clusters (overall and within CMS categories), the immunoSEQ TCR clonality was the highest in the ICR-high group and in the CMS1 / immune group among CMS subtypes, which has the highest proportion of ICR-high tumors. Using the whole transcriptome (18,270 genes), six out of the top ten genes positively correlating with TCR immunoSEQ clonality were represented by ICR genes (IFNG, STAT1 , IRF1 , CCL5, GZMA and CXCL10). Furthermore, the network of the top 50 genes correlating with immunoSEQ TCR clonality were centered on the ICR master regulators IRF1 and STAT1. The correlation of immunoSEQ TCR clonality with most of the ICR genes was stronger compared to the one observed with markers of tumor-reactive CD8+ T cells defined by singlecell sequencing approaches.

[0218] For nine patients, immunoSEQ TCR profiles were available on both the tumor and matched healthy colon tissue. This allowed the definition of overlap between T cell clones observed in the tumor and healthy colon sample for each of these patients. The proportion of tumor-enriched T cell clones correlated with ICR score (r = 0.75, P = 0.019; FIG. 2B). Thisimplies that the T cell clones infiltrating ICR-high tumors are highly divergent from those infiltrating healthy tissue, whereas T cells in ICR-low tumors are also present in healthy tissue.

[0219] In conclusion, our analyses demonstrated that the ICR signature captures the presence of tumor-enriched, clonally expanded T cells, possibly explaining its prognostic connotation.

[0220] Somatic alterations associated with weak immune response

[0221] Identification of potential drivers of immune responsiveness related to cancer cell somatic alterations, such as mutations and copy-number variations was sought by performing WES on 281 tumor samples and corresponding healthy tissue.

[0222] In terms of somatic mutations, the tumor mutational burden (TMB) of the AC-ICAM dataset was highly comparable to the TCGA-COAD cohort, as were the clinicopathological parameters. Unlike the TCGA-COAD cohort, ho wever, inclusion of samples in our study did not depend on tumor purity. In fact, stromal and immune content (ESTIMATE score) and the infiltration of individual lymphocyte subpopulations (FIG. 3) was significantly increased in the AC-ICAM compared to the TCGA-COAD datasets, whereas the opposite was observed for cancer-cell-intrinsic signatures. This was paralleled by a lower proportion of CMS1 and a higher proportion of ICR low in the TCGA-COAD compared to AC-ICAM. While the same proportion of MSI-high (MSI-H) cases was observed in the two cohorts, MSI-H TCGA-COAD samples displayed lower levels of CD8+ T cells, which is consistent with a positive selection of less-immune-infiltrated specimens. Then subsampled the cohort 100 times using two methodologies: one was random and the other was on a subgroup of samples with an ESTIMATE distribution that approximates that of the TCGA-COAD. The random subsampling resulted in tripling the number of subsets in which the Cox proportional regression sho wed a statistically significant survival benefit of the ICR score compared to the sampling method approximating the TCGA-COAD ESTIMATE distribution (P < 0.0001 , chi-squared test) (Supplementary Figs. 7 and 8). These findings suggest that a lower immune-stroma infiltration could have an impact on survival analysis, contributing to the lack of correlation between immune traits and OS observed in TCGA-COAD.

[0223] An overview of the somatic alterations landscape of the AC-ICAM cohort is as follows. Identified eight cancer-related genes with a mutation frequency of <5% in TCGA-COAD36 and Nurses’ Health Study (NHS)-Health Professionals Follow-up Study (HPFS) cohorts that were enriched in AC-ICAM and that had not been previously reported as colon canceroncogenic mediators or cancer driver genes for colorectal cancer. These genes were SEC31A, BAX, ZMYM2, BCR, ATF7IP, BCL1 B, ERCC5, and PARK2.

[0224] Overall, observed somatic mutations in 42 genes associated positively (P < 0.05) with ICR score, whereas no mutations were enriched in samples with a lower ICR score. When stratified the analysis according to the hypermutation status, identified gene mutation frequencies that were associated with both a higher or lower ICR score. Mutations of MAP3K1 , which were previously associated with low ICR in breast and pan-cancerTCGA analysis'! 0,21 , were the only ones with a negative correlation with ICR score in both hypermutated and nonhypermutated cancers in AC-ICAM. In hypermutated tumors, mutations in the homologous recombination repair genes BRCA1 , BRCA2 and FANCA and the mucinous histology were associated with a lower ICR score, consistently with the previously reported enrichment of BRCA1 and BRCA2 somatic mutations in mucinous colorectal tumors.

[0225] With respect to somatic copy-number genomic aberrations (SCNAs), no clear association was observed with ICR immune classification as they were dependent primarily on the mutational load / MSI status and secondarily on the CMS status.

[0226] Altogether, this analysis identified a relationship between specific cancer-related genes and / or histological characteristics and a lower level of intratumoral immune activation.

[0227] Genetic immune editing refines the prognostic value of ICR

[0228] Proceeded by integrating ICR and TMB data. While hypermutated samples frequently displayed an ICR-high phenotype, a considerable proportion of ICR-high samples (46%) had a low TMB, which did not impact the OS within or across ICR classes (FIG. 4A), conforming with what had been previously observed for Immunoscore.

[0229] While no difference observed in OS between high versus low TMB tumors, the presence of genetic immunoediting (GIE; calculated as the ratio of the observed versus the expected number of neoantigens) was nevertheless associated with improved OS. Then explored a composite score, called the immunoediting score (IES), based on both ICR cluster assignment and presence or absence of GIE (IES1 = ICR low and no GIE; IES2 = ICR low and GIE; IES3 = ICR high and no GIE; IES4 = ICR high and GIE) (FIG. 4B), similar to what was proposed in metastatic colon cancer by combining the Immunoscore and GIE41 .

[0230] It is proposes that the combination of the two parameters may more accurately reflect the presence of an active, antitumor immune response. Consistently with this hypothesis, a progressive increase of OS was observed from IES1 to IES4 (FIG. 4C). The additive value ofcombining ICR with GIE was confirmed in ICR-medium samples, which served here as an internal validation. While the TMB was higher in GIE versus non-GIE samples, GIE was observed in a significant proportion of both hypermutated and non-hypermutated tumors (55.1 versus 38.7%). Patients with IES4 tumors, of which ~50% were hypermutated or MSI-H, indeed demonstrated improved survival, with similar survival across stage I— III. No conclusion could be made in the IES4 stage IV subgroup as it only included two patients. No statistically significant difference was observed in terms of stage distributions and IES (chi-squared test, P = 0.46). IES remained significantly associated with OS in a multivariable Cox model corrected by stage (P = 0.045). IES categories also differed in term of TCR clonality, with increasing clonality from IES1 to IES4. The same trend was observed within the ICR-medium subgroup, in which the TCR clonality was increased (although not significantly) in the GIE samples compared to the non-GIE samples. The positive correlation between IES and TCR clonality was statistically significant when corrected for ICR score using multiple regression analysis and was confirmed by local polynomial regression analysis. Overall, these results suggest that the level of immune editing (IES) accurately reflects the level of a protective antitumor immune response driven by clonally expanded T cells.

[0231] Microbiome composition in healthy and colon cancer tissue

[0232] Sequenced the 16S rRNA gene using DNA extracted from matched tumor and healthy colon tissues from 246 patients (AC-ICAM246 cohort). This dataset was used for the microbiome landmark analysis.

[0233] Whole-genome sequencing (WGS, median coverage 76x) was performed in a subgroup of these samples (n = 167) for technical validation. For validation purposes, once the landmark analysis was completed, analyzed 16S rRNA gene-sequencing data from 42 additional tumor samples for which no matched normal DNA was available for this assay (referred to herein as ICAM42 cohort). After applying the same abundance filter to AC- ICAM246 and TCGA-COAD datasets, AC-ICAM captured all the genera detected in TCGA- COAD13, which displayed almost identical co-correlation patterns in the two cohorts, in additional to several other genera.

[0234] First, compared the relative abundance of taxa between matched tumor and healthy colon tissues. At the phylum level, observed a significant increase of Fusobacteria in tumor compared to healthy samples with a high concordance between the two methods. At the genus level, the strongest changes were observed for Fusobacterium, which was mostly represented by F. nucleatum. Our analysis captured several additional taxa highly enrichedin either tumor or healthy tissues (false discovery rate (FDR) < 0.05 and fold change > 2). No major difference in a diversity (the variety and abundance of species within an individual sample) was observed between tumor and healthy samples and only a modestly reduced microbial diversity was observed in ICR-high versus ICR-low tumors. Selenomonas and Selenomonas 3 were the taxa most significantly increased in ICR-high versus -low tumors. In terms of survival analysis, the highest number of nominally significant associations was obtained using tumor data (rather than healthy colon data) and OS as the end point. Fusobacterium and F. nucleatum abundances were associated with advanced stage, presence of BRAF mutations, MSI-H status, and a trend toward worse PFS survival, as previously observed. Instead of a negative correlation with T cells, Fusobacterium or F. nucleatum abundances were associated with cytotoxic T cells and NK cells paralleled by an increase of myeloid markers and signaling (for example, CD68, TREM1 and IL8 signature). The lack of association with a favorable outcome might be explained by the ability of F. nucleatum to inhibit T and NK cell killing of tumor cells by binding and activating the inhibitory receptors TIGIT45 and CEACAM1 or by induction of IL-8-mediated myeloid activation.

[0235] A microbiome signature (MBR score) predictive of survival

[0236] To detect clinically relevant associations between the microbial repertoire and clinical outcome, aim was identifying a microbiome signature predictive of survival using genus-level data from 16S rRNA gene sequencing, as part of our landmark microbiome analysis (AC- ICAM246, n = 246, testing set). On the AC-ICAM246, ran a multivariable elastic-net OS Cox regression model that selected 41 features (taxa) with a coefficient different to zero (associated with differential risk of death; Methods). This list of taxa and associated coefficients was termed the MBR classifier (FIG. 5A). A score was assigned to each sample (MBR score) by applying the MBR classifier. The MBR score displayed stability across different anatomic locations (in both tumor and healthy samples, despite the variable abundances of some taxa with respect to anatomic location.

[0237] Co-abundance network inference using SparCC48 correlation coefficients revealed five distinct clusters of taxa. Taxa enriched in ICR-high versus ICR-low samples or in tumor versus healthy colon samples displayed high co-abundance (enriched in C3) and the same was observed for taxa enriched in healthy colon or in ICR-low samples. Low and high-risk taxa (according to MBR classifier) were spread across the different clusters. Only marginal differences in survival were observed using estimates based on the cumulative abundance of genera belonging to each cluster identified by the network analysis. The only survivalassociation with an FDR <0.1 was detected for C5 (OS analysis, P = 0.017, hazard ratio (HR) 1.6, high versus low abundance, FDR = 0.085). C5 was constituted by three taxa, including one MBR-high-risk genera and no MBR- low- risk genera. Overall, these results suggest that clinical outcome is influenced by microbiome diversity, which is captured by the MBR classifier.

[0238] Consistently, a high a diversity was associated with a prolonged OS FDR < 0.05 for all the a diversity estimates. Because of the strong contribution of Ruminococcus 2 to the MBR classifier, sought to identify the actual Ruminococcus species. In WGS data, the Ruminococcus genus mostly consisted of Ruminococcus bromii, which also had the strongest correlation with Ruminococcus 2. R. bromii presence was confirmed by PCR, which had strong correlation with sequencing data (for example, 91 % concordance between WGS and PCR).

[0239] Validation of the MBR score

[0240] A low MBR score (MBR < 0, MBR low), in our training cohort (ICAM246, training set) was associated with a considerable (85%) reduction of risk of death (FIG. 5B). Confirmed the association between MBR low (risk) and prolonged OS in two independent testing sets (ICAM42 and TCGA-COAD cohorts), individually and combined (Fig. 5B, testing sets). The performance of the final MBR model was lower on the test sets than on the training set, which is typical for machine-learning models; ho wever, the concordance index of the final MBR model in both the test sets were superimposable to the ones obtained via cross-validation of the best MBR model on the training set, substantiating that the model can generalize well to new (unseen) data.

[0241] A similar, but less-pronounced trend in terms of reduction of the risk of death was detected by simply using intratumoral Ruminococcus 2 (based on 16S data) or R. bromii presence (based on either PCR or WGS data). Intratumoral Ruminococcus 2 and MBR score, which strongly correlated with each other, were similar in tumor and healthy colon tissues (FIG. 5C).

[0242] The relationship between the microbiome and clinical outcome pointed to an interaction between the microbiome and biological processes occurring in the tumor. When correlating immune trait values with the MBR score, the strongest (inverse) correlation with the MBR score was observed for signatures capturing the prevalence of CD103+ dendritic cells (DCs) with unique antigen processing and presentation capabilities for efficient antigencross-presentation to CD8+ T cells (CD103+, mean signature (P = 0.003) and CD103+ signature to CD103- signature ratio (P = 0.001)). Consistently, correlation analyses between individual taxa included in the MBR classifier and immune traits demonstrated, with few exceptions, a positive correlation with myeloid signatures and a negative correlation with the CD103+ / - ratio for taxa with positive MBR coefficient (higher risk of death), while the reverse was observed for taxa with a negative MBR coefficient.

[0243] Development and validation of the mICRoScore

[0244] Development of a multi-omics parameter was sought that could capture a subgroup of patients with exceptional survival. Among single-omics parameters that were significant in the univariate Cox regression OS analysis (ICR, MBR and GIE categories), only ICR and MBR were retained by the multivariable Cox models (P < 0.05; Supplementary Table 9) adjusted for age, CMS subtypes, stage and MSI status. MBR and ICR were therefore combined into an integrated score (mICRoScore).

[0245] Indeed, in the training cohort (AC-ICAM246), the co-presence of ICR high and MBR low (mICRoScore high) identified a subgroup of patients with a 97% 5-year OS, with only three deaths detected at a later follow-up (FIG. 6A) that were not related to colon cancer. No deaths were observed during the entire follow-up in patients with mICRoScore high in the TCGA- COAD cohort (n = 107, testing set; FIG. 6B). In both the training (AC-ICAM) and the testing (TCGA-COAD) sets, the mICRoScore-high group consisted of patients at different stages. The additive effect of the two parameters was due to the ability of MBR to segregate ICR high into two distinct risk categories (FIG. 6C and FIG. 6D).

[0246] Conclusion

[0247] A multi-omics approach allowed thorough examination of the molecular characteristics of immune responsiveness in colon cancer and uncovered interactions between the microbiome and the immune system. It was found that a TH1 cell / cytotoxic immune activation, as captured by the ICR, immunoediting, concurrent expansion of TCR clonotypes and specific intratumoral microbiome composition, were associated with a favorable clinical outcome. ICR was associated with OS independently of MSI and CMS, which both lost statistical significance in the multivariate analysis.

[0248] Its prognostic impact increased when combined with a metric capturing the genetic immunoediting (GIE / IES). Using deep TCR sequencing in tumor and healthy tissues, thedata showed that the prognostic effect of ICR could be due to its ability to capture the presence of tumor-enriched and possibly tumor-antigen specific, T cell clones.

[0249] The AC-ICAM addressed the limitations of the TCGA colon cancer cohort noted by the scientific community and corroborated by our comparative analyses. While several studies have described associations between response to immunotherapy and the gut microbiome and identified cancer-specific microbiome compositions, comprehensive microbiome analyses focused on patients with primary colon cancer are lacking. By analyzing the tumor microbiome composition using 16S rRNA gene sequencing in AC-ICAM samples, identified a microbiome signature (MBR risk score) with strong prognostic value. This signature was derived from tumor samples, but there was a strong correlation between the healthy colon and tumor MBR risk scores, suggesting that this signature may capture the patient’s gut microbiome composition.

[0250] By combining the ICR and MBR scores, identified and validated a multi-omics biomarker (mICRoScore) that could predict exceptionally long survival in patients with colon cancer.

[0251] Both the mICRoScore and IES could be tested in the context of cancer immunotherapy as predictive biomarkers. Data from the NIBIT-M4 trial and publicly available datasets suggest that the combination of the genetic immunoediting and ICR (IES) has predictive value in melanoma patients treated with immune checkpoint inhibitors.

[0252] By analyzing the tumor microbiome composition using 16S rRNA gene sequencing in AC-ICAM samples, identified a microbiome signature (MBR risk score) with strong prognostic value. This signature was derived from tumor samples, but there was a strong correlation between the healthy colon and tumor MBR risk scores, suggesting that this signature may capture the patient’s gut microbiome composition.

[0253] Additional analysis and technical validation using orthogonal platforms such as WGS and PCR indicated that the detected signal was driven by R. bromii. Correlation analyses between the MBR risk score and immune traits suggest a specific positive modulation of CD103+dendritic cells, which are critical for antitumor immune responses. It is believed that the identified consortium of bacteria favors optimal T cell priming mediated by CD103+dendritic cell activation and suppression of the myeloid compartment, leading to the induction of a partially protective antitumor immunity.

[0254] Unless otherwise indicated, all numbers expressing quantities of ingredients, properties such as molecular weight, reaction conditions, and so forth used in the specification and claims are to be understood as being modified in all instances by the term “about.” As used herein the terms "about" and “approximately” means within 10 to 15%, preferably within 5 to 10%. Accordingly, unless indicated to the contrary, the numerical parameters set forth in the specification and attached claims are approximations that may vary depending upon the desired properties sought to be obtained by the present invention. At the very least, and not as an attempt to limit the application of the doctrine of equivalents to the scope of the claims, each numerical parameter should at least be construed in light of the number of reported significant digits and by applying ordinary rounding techniques. Notwithstanding that the numerical ranges and parameters setting forth the broad scope of the invention are approximations, the numerical values set forth in the specific examples are reported as precisely as possible. Any numerical value, however, inherently contains certain errors necessarily resulting from the standard deviation found in their respective testing measurements.

[0255] The terms “a,” “an,” “the” and similar referents used in the context of describing the invention (especially in the context of the following claims) are to be construed to cover both the singular and the plural, unless otherwise indicated herein or clearly contradicted by context. Recitation of ranges of values herein is merely intended to serve as a shorthand method of referring individually to each separate value falling within the range. Unless otherwise indicated herein, each individual value is incorporated into the specification as if it were individually recited herein. All methods described herein can be performed in any suitable order unless otherwise indicated herein or otherwise clearly contradicted by context. The use of any and all examples, or exemplary language (e.g., “such as”) provided herein is intended merely to better illuminate the invention and does not pose a limitation on the scope of the invention otherwise claimed. No language in the specification should be construed as indicating any non-claimed element essential to the practice of the invention.

[0256] Groupings of alternative elements or embodiments of the invention disclosed herein are not to be construed as limitations. Each group member may be referred to and claimed individually or in any combination with other members of the group or other elements found herein. It is anticipated that one or more members of a group may be included in, or deleted from, a group for reasons of convenience and / or patentability. When any such inclusion ordeletion occurs, the specification is deemed to contain the group as modified thus fulfilling the written description of all Markush groups used in the appended claims.

[0257] Specific embodiments disclosed herein may be further limited in the claims using consisting of or consisting essentially of language. When used in the claims, whether as filed or added per amendment, the transition term “consisting of’ excludes any element, step, or ingredient not specified in the claims. The transition term “consisting essentially of’ limits the scope of a claim to the specified materials or steps and those that do not materially affect the basic and novel characteristic(s). Embodiments of the invention so claimed are inherently or expressly described and enabled herein.

[0258] Furthermore, numerous references have been made to patents and printed publications throughout this specification. Each of the above-cited references and printed publications are individually incorporated herein by reference in their entirety.

[0259] In closing, it is to be understood that the embodiments of the invention disclosed herein are illustrative of the principles of the present invention. Other modifications that may be employed are within the scope of the invention. Thus, by way of example, but not of limitation, alternative configurations of the present invention may be utilized in accordance with the teachings herein. Accordingly, the present invention is not limited to that precisely as shown and described.

Claims

CLAIMSWe claim:1 . A method for predicting survival in patients suffering from colon cancer, the method comprising: a) sequencing 16S rRNA genes using DNA extracted from tumor samples of patients suffering from colon cancer; b) detecting the presence of Ruminococcus bromii DNA sequences in the tumor samples of patients suffering from colon cancer; and c) identifying patients suffering from colon cancer having higher levels of Ruminococcus bromii DNA sequences in the tumor samples as having a higher probability of survival.

2. A method for predicting survival in patients suffering from colon cancer, the method comprising: a) measuring expression levels of the genes CCL5, CD274, CD8A, CD8B, CTLA4, CXCL9, CXCL10, FOXP3, GNLY, GZMA, GZMB, GZMH, IDO1 , IFNG, IL12B, IRF1 , PDCD1 , PRF1 , STAT 1 , and TBX21 in tumor samples of patients suffering from colon cancer and assigning an immunologic constant of rejection (ICR); b) detecting the abundance of specific microbes in the tumor samples of patients suffering from colon cancer and assigning a microbiome signature risk (MBR) score; and c) identifying patients suffering from colon cancer who have a high ICR and low MBR score as having a higher probability of survival.

3. The method of claim 2, wherein, in step b), the presence of one or more microbes from genera UBA1819 (Ruminococcaceae), Lachnospiraceae NK4A136 group (Lachnospiraceae), Erysipelotrichaceae UCG-003 (Erysipelotrichaceae), Sellimonas (Lachnospiraceae), Coprococcus 1 (Lachnospiraceae), Terrisporobacter (Peptostreptococcaceae), Olsenella (Atopobiaceae), Oscillibacter (Ruminococcaceae), Dialister (Veillonellaceae), Parabacteroides (Tannerellaceae), Parvimonas (Family XI), [Eubacterium] eligens group (Lachnospiraceae), Bifidobacterium (Bifidobacteriaceae), Hungatella (Lachnospiraceae), Mogibacterium (Family XIII), Caproiciproducens (Ruminococcaceae), Erysipelatoclostridium (Erysipelotrichaceae), Halomonas (Halomonadaceae), and Odoribacter (Marinifilaceae) is detected.

4. The method of claim 2, wherein, in step b), the presence of one or more microbes from genera Ruminococcaceae NK4A214 group (Ruminococcaceae), Solobacterium (Erysipelotrichaceae), Ruminococcus 2 (Ruminococcaceae), Treponema 2 (Spirochaetaceae), Streptococcus (Streptococcaceae), Phascolarctobacterium (Acidaminococcaceae), Moryella (Lachnospiraceae), Coprococcus 2 (Lachnospiraceae), Ruminococcus 1 (Ruminococcaceae), [Eubacterium] ventriosum group (Lachnospiraceae), Actinomyces (Actinomycetaceae), CAG-56 (Lachnospiraceae), Ruminiclostridium 6 (Ruminococcaceae), Barnesiella (Barnesiellaceae), Prevotella 9 (Prevotellaceae), Family XIII AD3011 group (Family XIII), Prevotellaceae, Peptostreptococcus (Peptostreptococcaceae), Candidatus Soleaferrea (Ruminococcaceae), Rhodospirillales, Leptotrichia (Leptotrichiaceae), Veillonella (Veillonellaceae), UBA1819 (Ruminococcaceae), Lachnospiraceae NK4A136 group (Lachnospiraceae), Erysipelotrichaceae UCG-003 (Erysipelotrichaceae), Sellimonas (Lachnospiraceae), Coprococcus 1 (Lachnospiraceae), Terrisporobacter (Peptostreptococcaceae), Olsenella (Atopobiaceae), Oscillibacter (Ruminococcaceae), Dialister (Veillonellaceae), Parabacteroides (Tannerellaceae), Parvimonas (Family XI), [Eubacterium] eligens group (Lachnospiraceae), Bifidobacterium (Bifidobacteriaceae), Hungatella (Lachnospiraceae), Mogibacterium (Family XIII), Caproiciproducens (Ruminococcaceae), Erysipelatoclostridium (Erysipelotrichaceae), Halomonas (Halomonadaceae), and Odoribacter (Marinifilaceae) is detected.

5. The method of claim 2, further comprising detecting the presence of genetic immunoediting (GIE), where the presence of GIE identifies a patient suffering from colon cancer as having a higher probability of survival.

6. The method of claim 2, wherein the low MBR score identifies patients within a group of patients having high ICR who have higher probability of survival.

7. A method for treating colon cancer in a patient, the method comprising administering a microbiome composition comprising Ruminococcus bromii, or other microbe species with a negative weight in an MBR calculation, to the patient.

8. The method of claim 7, wherein the microbiome composition comprises one or more microbes selected from genera UBA1819 (Ruminococcaceae), LachnospiraceaeNK4A136 group (Lachnospiraceae), Erysipelotrichaceae UCG-003 (Erysipelotrichaceae), Sellimonas (Lachnospiraceae), Coprococcus 1 (Lachnospiraceae), Terrisporobacter (Peptostreptococcaceae), Olsenella (Atopobiaceae), Oscillibacter (Ruminococcaceae), Dialister (Veillonellaceae), Parabacteroides (Tannerellaceae), Parvimonas (Family XI), [Eubacterium] eligens group (Lachnospiraceae), Bifidobacterium (Bifidobacteriaceae), Hungatella (Lachnospiraceae), Mogibacterium (Family XIII), Caproiciproducens (Ruminococcaceae), Erysipelatoclostridium (Erysipelotrichaceae), Halomonas (Halomonadaceae), or Odoribacter (Marinifilaceae).

9. A method of identifying a patient as suffering from colon cancer, the method comprising detecting the presence of one or more genes selected from SEC31 A, BAX, ZMYM2, BCR, ATF7IP, BCL1 B, ERCC5, or PARK2 in a tumor sample of the patient.

10. A method of identifying a patient suffering from colon cancer and having a high ICR as having at least 95% probability of overall survival over a period of at least 5 years, the method comprising: b) detecting the abundance of specific microbes in the tumor samples of patients suffering from colon cancer and assigning a microbiome signature risk (MBR) score; and c) identifying patients suffering from colon cancer and having a high ICR and having a low MBR score as having at least 95% probability of overall survival over a period of at least 5 years.

11. A system for identifying patients suffering from colon cancer and having a high ICR as having at least 95% probability of overall survival over a period of at least 5 years, the method comprising: a) using AC-ICAM246, n=246, as a testing set to develop MBR scores; b) using AC-ICAM42 and TCGA-COAD, as the validation sets to validate the MBR scores of step a); and c) identifying patients suffering from colon cancer and having a high ICR and having a low MBR score as having at least 95% probability of overall survival over a period of at least 5 years.

Images