Data from 188 patients with stage II CRC and tissue samples were collected from the institutional database of the Fudan University Shanghai Cancer Center between 2013 and 2019. The American Joint Committee on Cancer staging system was used to determine each patient’s stage [15]. Tested by next-generation sequencing, 24 patients were classified as MSI-H. None of the patients had radiation therapy, chemotherapy, or immunotherapy before tumor resection. Clinical and pathological data were obtained from patient records and postoperative pathology reports.

Sections (4 mm thick) were cut from formalin-fixed, paraffin-embedded CRC tissue and control tonsil tissue for multiplex immunohistochemistry (mIHC). The slides were dewaxed in xylene, rehydrated, and rinsed in graded ethanol solutions and tap water. Antibody diluent/block (72424205; PerkinElmer) was applied to block endogenous peroxidase. The slides were boiled in a Tris-EDTA buffer (pH: 9; 643901; Klinipath) and underwent microwave treatment (MWT) for antigen retrieval. Information on the primary antibodies and the corresponding fluorophores is provided in Table S1 in Multimedia Appendix 1, including 2 panels (Figure 1). One antigen required 1 round of labeling, including primary antibody incubation, secondary antibody incubation, and tyramide signal amplification (TSA) visualization, followed by labeling of the subsequent antibody. After incubation with the primary antibody for 1 hour at room temperature, the slides were incubated with Opal Polymer HRP Ms+Rb (2414515; PerkinElmer) at 37 ℃ for 10 minutes. TSA visualization was performed with the Opal 7-Color IHC Kit (NEL797B001KT; PerkinElmer) containing the fluorophores 4,6-diamidino-2-phenylindole (DAPI; Thermo Scientific) and the TSA Coumarin system (NEL703001KT; PerkinElmer). MWT was performed to remove the antibody-TSA complex with the Tris-EDTA buffer (pH: 9). TSA single-stained slides were finished with MWT, counterstained with DAPI for 5 minutes, and enclosed in Antifade mounting medium (I0052; NobleRyder).

Multiplexed and single-color control slides were scanned at an absolute magnification of 200× by the PerkinElmer Vectra automated multispectral microscope. Representative fields from the single-color slides were imaged, and a spectral library for unmixing was generated by inForm image analysis software (version 2.1; PerkinElmer). Index cases were stained using the multiplex method and then imaged. Channels were unmixed using the spectral library. All settings were saved within an algorithm to allow for batch analysis of multiple original multispectral images of the same tissue [16].

The nuclear morphological features were based on DAPI staining. The numbers of immune cells in each image were scored as percent cellularity (number of positive cells/number of nucleated cells). Five representative fields at 200× magnification of tissue area were selected. The densities of immune cells were segmented independently by 2 pathologists. Immune variables were classified based on the patterns of fluorochrome intensity.

Patients were monitored every 3‐6 months for 3 years, then every 6‐12 months up to 5 years. Follow-ups included rectal exams, carcinoembryonic antigen (CEA) tests, annual radiological studies, and colonoscopies as needed.

The ColonCore panel (Burning Rock) is designed for simultaneous detection of MSI status and mutations in 37 CRC-related genes (Table S2 in Multimedia Appendix 1). The MSI detection method was a read-count-distribution–based approach, using the coverage ratio of a specific set of repeat lengths as the main characteristic of each microsatellite locus. The MSI status of a sample was determined by the percentage of unstable loci in the given sample [17].

Statistical analysis was performed and visualized by R (version 3.4.3; R Foundation for Statistical Computing), SPSS software (version 25.0; IBM Corp), and GraphPad Prism 7 software (GraphPad Software Inc). All group-wise comparisons were conducted by the 2-sided unpaired Mann-Whitney U test, followed by the Bonferroni procedure. The Cox proportional hazards regression model was used to assess the hazard ratios, 95% CIs, and P values for univariate and multivariate analysis. Variables with P<.10 after adjusting for common clinicopathological parameters were included in the multivariate analysis. Survival times were compared using the log-rank test. A P value of <.05 was considered statistically significant, and all P values corresponded to 2-sided statistical tests.

The process from feature engineering to model evaluation is depicted in Multimedia Appendix 2. Categorical features were one-hot encoded as dummy variables. The mutation landscape was also one-hot encoded based on gene classes, with 2 classification stringencies, using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) gene functional classification tool [18]. To further engineer the mutation landscape, we calculated the joint posterior mutation probability P(MSI|M) with the following equation:

where P(MSI|M) is the probability that a patient has MSI-high status given their mutation landscape and is based on previous probabilities and the frequencies of mutated genes in MSI-H and MSS populations, gi is a mutated gene in a patient’s sample, M=gi is a set of all detected mutated genes in the sample, PMSI is the probability of MSI-H in a CRC sample (0.83), PMSS is the probability of MSS in a CRC sample (0.17), Pgi|MSI is the frequency of a mutated gene gi in a CRC MSI-H population, and Pgi|MSS is the frequency of a mutated gene gi in a CRC MSS population. Pgi|MSI and Pgi|MSS were based on the findings of Serebriiskii et al [19] and 2 datasets [20,21] on cBioPortal [22,23]. Our Bayesian-based metric can explicitly incorporate previous biological knowledge, including MSI-H/MSS prevalence in populations with CRC and microsatellite status–specific mutation frequencies, and probabilistic reasoning into the modeling process. Leveraging these priors potentially enhances the model’s ability to distinguish MSI-H from MSS cases. This metric was included in the dataset along with the other features and used for model training.

Though no missing data were presented in the dataset, our model can handle missing input from users because we trained several models with varied complexities, as elaborated in the following section.

Sample microsatellite status was one-hot labeled (MSI-H as 1 and MSS as 0). Multiple extreme gradient boosting (XGBoost) models [24] were trained to identify MSI-H–like tumors with different combinations of features, that is, patient metainformation, mutational landscape–derived features, and mIHC results, such as PD-L1, CD163, and CD8 mIHC staining results. The predicted likeliness of MSI-H by the models was defined as MSI-H score. Specifically, 44 models, with different combinations of features shown in Table S3 in Multimedia Appendix 1, were trained. To address potential bias caused by class imbalance, the scale_pos_weight parameter was set to the class ratio during model training.

The models were then deployed as a public web interface. As users’ data privacy is prioritized, users’ input is never stored on our server. Each model ensemble consists of 10 submodels, each trained with a distinct random seed. The final prediction for any user input is the average of the middle 6 submodel outputs, excluding the extremes. One of the XGBoost tree models in pseudocode is shown in Multimedia Appendix 3.

To understand the feature importance in an unbiased and holistic way, a massive exploratory XGBoost model with all features was trained. These features include the features patient metainformation, mutational landscape–derived features, and all mIHC results. The massive model can capture the full spectrum of the MSI-H score variation and avoid the potential bias or noise introduced by the feature selection process. The feature importance was then computed using both the XGBoost built-in function and the Shapley additive explanations package [25].

Following classification by this model (threshold=0.3 defining MSI-H, chosen so that predicted MSI-H proportion approximates the epidemiologically documented prevalence of MSI-H CRC), we visualized all mIHC features of all samples by grouped box plots. For each cell type, we compared 4 groups (all MSS vs MSI-H, other MSS vs MSI-H, other MSS vs MSI-H–like MSS, and MSI-H–like MSS vs MSI-H) by 2-sided unpaired Mann-Whitney U test and Benjamini-Hochberg adjustment [26]. To cluster cells based on 4 comparisons, we projected each cell type into a 4D latent space using a formula measuring similarity between cell percentages of former and latter populations:

where r_former and r_latter represent the percentages of a specific cell type (eg, CD8+ cells) measured in individual samples belonging to the “former” and “latter” groups, respectively. p_adj is the adjusted P value of a comparison. We then performed hierarchical clustering based on Euclidean distance in latent space with complete linkage [27].

Model generalizability was assessed by training models with identical hyperparameters through stratified 5-fold cross-validation. Cohen κ coefficients were computed on each hold-out fold, and the mean Cohen κ was computed based on the 5 κ’s, with greater κ’s indicating better model performance. This training and validation process was repeated 10 times with different random stratified splits and model initializations.

Ethics approval was obtained from the Ethics Committee of Fudan University Shanghai Cancer Center, and informed consent was obtained from all participants (1808190‐12). Neither the patients nor the public were involved in this study (ie, only database tissue samples and data from patient records and postoperative pathology reports were used). All patient data collected for this study were deidentified prior to analysis.

TILs were analyzed using the mIHC method. Significant differences were found between the stromal region and the tumor region (Figure 2). The stromal region showed a higher prevalence of CD3+ T cells (P<.001), CD8+ T cells (P<.001), memory T cells (P<.001), CD8+ memory T cells (P<.001), CD3+ PD-1+ T cells (P<.001), CD4+ T cells (P=.048), regulatory T cells (Tregs; P<.001), macrophages (P=.001), M1 macrophages (P=.003), M2 macrophages (P=.001), and PD-L1+ cells (P=.007) than the tumor region. However, no significant difference was observed for CD8+ PD-1+ T cells and PD-L1+ macrophages between the stromal region and tumor region.

Variables (Table 1) commonly collected in clinics and related to prognosis or with a P value <.10 in univariate analysis (Tables S4 and S5 in Multimedia Appendix 1) were analyzed using the Cox proportional hazards regression model (Table S6 in Multimedia Appendix 1). Microsatellite status was not linked to overall survival or disease-free survival in multivariate analysis. Significant overall survival predictors included age, CEA, M1 macrophage (CD68+ CD163–) infiltration in stromal region, and PD-1+ T cell (CD3+ PD-1+) infiltration in tumor region. Significant disease-free survival predictors included age, CEA, tumor differentiation, and CD8+ T cell (CD8+) infiltration in tumor region (Table S6 in Multimedia Appendix 1).

The mIHC experiment was performed to examine the TILs in all CRCs (Figure 3). MSI-H CRCs exhibited significantly higher infiltration of TILs compared to MSS CRCs in both the tumor region and stromal region. Specifically, MSI-H CRCs had a more abundant presence of PD-L1+ M2 macrophages (P=.001), CD163+ cells (P=.001), PD-L1+ macrophages (P=.01), M2 macrophages (P=.001), and macrophages (P=.03) in the stromal region, as well as M2 macrophages in the tumor region (P=.02).

MSI-H–like MSS CRCs exhibited TIL infiltration patterns akin to MSI-H CRCs but distinct from other MSS CRCs. No significant difference was observed between MSI-H–like MSS and genuine MSI-H CRCs (the fourth column in Figure 3). Compared to other MSS CRCs, MSI-H–like MSS CRCs showed higher infiltration of CD163+ cells in the stromal region (P=.001; the box plot and the third column in Figure 3) and potentially increased levels of PD-L1+ M2 macrophages (P=.13), FOXP3+ cells (P=.09), Tregs (P=.09), PD-L1+ macrophages (P=.16), M2 macrophages (P=.09), and macrophages (P=.16) in the stromal region, as well as M2 macrophages (P=.13) in the tumor region. The distinct infiltration patterns of TILs indicate that heightened presence of macrophages and Tregs are key factors in distinguishing MSI-H–like MSS CRCs from MSS CRCs.

Macrophages were also found to be significantly more abundant in genuine MSI-H CRCs than in other MSS CRCs (the second column in Figure 3). Specifically, in the stromal region, PD-L1+ M2 macrophages (P<.001), CD163+ cells (P<.001), PD-L1+ macrophages (P=.004), M2 macrophages (P<.001), M1 macrophages (P=.045), CD4+ T cells (P=.045), and macrophages (P=.01) were significantly more abundant in MSI-H CRCs than in other MSS CRCs. In the tumor region, M2 macrophages (P=.007), M1 macrophages (P=.045), and macrophages (P=.045) were found to be significantly increased in MSI-H CRCs compared with other MSS CRCs.

The TIL distribution shows that the model performed well. The scorer, which was trained and validated on only 3 types of lymphocytes, classified MSI-H–like MSS CRC samples with similar TIL distributions as MSI-H CRC samples (the fourth column in Figure 3) rather than MSS CRC samples (the third column in Figure 3), despite most features of TILs (15 other lymphocytes) being unknown to the model. In addition, as anticipated, the heat map in Figure 3 (second and third columns) revealed that other MSS CRC samples exhibited a slightly closer TIL distribution to MSI-H–like MSS CRC samples than to MSI-H CRC samples.

The feature importance in a large predictive model is described in Multimedia Appendix 4. TIL features from the whole or stromal region were more predictive than the tumor region alone for the MSI-H status. Top features for MSI-H score predictor included macrophage subtypes, mutational landscape, and immune cell distributions.

Increasing the number of features generally enhanced κ¯, indicating better generalization performance (Figure 4). However, models incorporating PD-L1 mIHC staining tended to exhibit lower κ¯ compared to those without, likely due to noise in PD-L1 measurements, as evidenced by the high SD of κ¯ for models 2 to 4. This noise effect was mitigated by increasing model complexity; for example, model 44 had a greater κ¯ than model 41, despite including PD-L1. Feature importance analysis (Figure 5) revealed that while PD-L1 remained relevant, its importance diminished as models became more complex, suggesting that sophisticated models learned to filter out noise and extract useful information from PD-L1. The variable Spearman correlation matrix heat map is shown in Multimedia Appendix 5. In addition, an MSI-H scorer web interface is freely accessible [28].

Immunotherapy has been successful for treating MSI-H CRCs but is not as effective in MSS CRCs, which comprise the majority of CRCs. Thus, we developed a machine learning–based MSI-H predictor to generate a robust and reliable score that can capture the complexity and heterogeneity of CRC and better target patients with MSS CRC who may benefit from immunotherapy. Our study also provides insights into the immune landscape of CRC and the role of immune cell distributions, clinical features, and gene mutations in influencing MSI status. This CRC prognostic study mostly agrees with our previous research [29] and with findings from other authors [30]. For example, according to our results, TIL infiltration, primarily by macrophages or CD163+ cells, was significantly higher in MSI-H CRCs than in MSS CRCs (Figure 3), consistent with previous studies [31].

We observed a higher abundance of TIL subsets in the stromal region than in the tumor region, indicating a more active immune response in the stromal region (Figure 2). Comparing the MSI predictive performance of models 42, 43, and 44 (Figure 4) also highlights the importance of stromal TILs. Regional disparities underscore the importance of analyzing the complete tumor region for comprehensive insights. Our scorer successfully identified MSI-H–like MSS samples with TIL distributions similar to genuine MSI-H CRCs (Figure 3). In addition, the balance of proinflammatory and anti-inflammatory scale is an important feature for immunological characters. Macrophages can be classified into 2 main subtypes: M1 macrophages with proinflammatory and antitumor functions and M2 macrophages with anti-inflammatory and protumor functions. The ratio of M1/M2 macrophages may influence immunotherapy outcomes, reflecting the balance between proinflammatory and anti-inflammatory signals in the tumor microenvironment [32]. Tregs are frequently known to be immunosuppressive and can predict both the host immune response and chemotherapeutic response [33]. Both macrophages and Tregs are important in the regulation of immunoactivity. As is shown in our results, the distribution of macrophages and Tregs appears to be important in differentiating MSI-H–like MSS CRCs from other MSS CRCs based on TIL infiltration patterns (Figure 3). By comparing κ̄ variations within model sets and between specific model set pairs (Figure 4)—including model 5/6/7 versus 13/16/19, 27/28/29 versus 35/38/41, 2/3/4 versus 12/15/18, 24/25/26 versus 34/27/40, 11/14/17 versus 20/21/22, and 33/36/39 versus 42/43/44—we observed that CD163 mIHC result increased the predictive value for whole-tumor MSI scores but reduced it for tumor region scores. To better understand the differences in M2 macrophages and Tregs between MSS CRCs and MSI-H CRCs, further research on their function in CRCs is necessary.

Our analysis revealed that PD-L1+ M2 macrophages in the total region, mutational landscape, CD163+ cells in the stromal region, PD-L1+ M2 macrophages in the stromal region, and tumor site were the most important features for predicting MSI-H status (Multimedia Appendix 4), aligning with other research results [2,3]. Macrophages can express PD-L1 and interact with PD-1+ T cells, which may affect the response to immunotherapy [34]. However, PD-L1+ macrophages potentially indicate M1-like polarization profiles [35]. The stroma is important because it can influence the extracellular matrix formation, angiogenesis, immune response, and therapeutic resistance of tumors [3]. The importance of the mutational landscape for prediction is widely known [36]. We did not find an obvious immunological explanation for why tumor site would impact the similarity between MSI-H and MSS. Further study is needed to clarify the underlying mechanisms. Moreover, we observed that feature number and type influenced the generalization ability of the MSI-H score prediction models (Figure 4 and Figure 5). This suggests that the omission of diverse variables requires specific computational models, and our machine learning scorer is adept at incorporating all such considerations, thereby highlighting our advantage.

On the basis of our results, we proposed a hypothesis regarding the changes that occur in MSI-H–like MSS CRCs compared to other MSS CRCs. MSI-H–like MSS CRCs foster an immunosuppressive microenvironment with M2 macrophages, Tregs, and PD-L1 that inhibits T cell responses [37]. However, there are enough T cells present that can be reactivated upon PD-1/PD-L1 blockade, leading to the sensitivity of MSI-H–like MSS CRCs to ICIs. The abundance of macrophages suggests that there may be some M1-like populations that, when disinhibited, promote antitumor immunity. Detailed differences in immune cell populations and their functions in MSI-H–like MSS CRCs and other MSS CRCs should be further investigated to understand the mechanisms underlying the differential response to immunotherapy. Furthermore, future clinical trials could be conducted to evaluate ICI treatment between patients with MSI-H–like MSS CRC and other patients with MSS CRC with low MSI-H scores.

Limitations of our study include the lack of internal or external validation of the MSI-H score in patients with MSS CRC receiving immunotherapy and the absence of investigation into the underlying molecular mechanisms. Further research and clinical trials are needed to validate our MSI-H score and elucidate the associated mechanisms.

In conclusion, our study revealed significant variations in TIL distribution across tumor regions and MSI status. Integrating clinical, TIL, and mutational data, we developed a robust MSI-H scorer that captures CRC’s complexity and heterogeneity. Macrophages, gene mutations, and tumor site emerged as key predictors. MSI-H–like MSS CRCs exhibited TIL infiltration patterns with high immunoactivity similar to MSI-H CRCs, distinctly different from other MSS CRCs. Our privacy-protected MSI-H score predictor is freely available on the web, enabling clinical and research applications.