StableMate: a statistical method to select stable predictors in omics data

BC gene expression data

We analysed the RNA-seq dataset from The Cancer Genome Atlas Program (TCGA-BRCA) to study the transcriptional regulation of ESR1 in ER+ BC, available from the R package TCGAbiolinks (18). The dataset includes the expression quantification of 60 660 genes on 113 normal samples and 1094 BC samples. We focused on the log-transcript-per-million (logTPM) of 19 937 protein-coding genes for analysis.

We used two other gene expression studies as validation: the microarray data of 2509 BC samples from the METABRIC cohort (19,20), available from cBioProtal (21) in the form of z-score relative to all samples (log), and RNA-seq data of 980 normal breast samples from GTEx (22) in the form of logTPM.

Colon cancer metagenomics data

We obtained nine colorectal cancer (CRC) case-control studies of faecal metagenome from the R package curatedMetagenomicData (23). We excluded two studies with a sequencing depth lower than the average 10 million reads per sample. The remaining studies included curated microbial species abundance and pathway abundance data from seven different countries and eight different cohorts: including 107 samples from Austria (24), 104 samples from the USA (25), 125 samples from Germany (26), 509 samples from Japan (27), 128 samples from China (28) and 114 samples from France (29), as well as two cohorts containing 53 and 60 samples from Italy (30). In total, all cohorts included 1429 samples. We filtered the species and pathway abundance data from each cohort down to 313 species and 431 pathways that were detected across all cohorts.

To normalize the abundance data, we applied rank transformation by calculating the within-sample ranking quantile of the abundance of each species (or pathway). A species ranked the qth most abundant in a sample is assigned the value (1 − q)/(p − 1), where p is the total number of species analysed. Therefore, the most abundant species of a sample has a rank transformed value of 1, and the least abundant species has a value of 0.

Glioblastoma single-cell RNA sequencing data

We analysed the glioblastoma (GBM) single-cell RNA sequencing (scRNA-seq) data from Darmanis et al. (31), who sequenced single cells sampled from four GBM patients at their tumour cores and surrounding peripheral tissues. The raw and curated read count data included 3589 cells measured on 23 368 genes available from http://gbmseq.org/. We retained 1874 cells of myeloid cell types, including 1329 cells sequenced from the core and 518 cells sequenced from the periphery for analysis. We used the R package Seurat to log-normalize the data and identify the most variable 2000 genes with the FindVariableFeatures function (32). We then imputed the log normalized data using Sincast imputation with default tuning (33) for StableMate variable selection and single-cell projection. Diffusion map (DM) and diffusion pseudotime (DPT) learning were performed on the original log-normalized data (without imputation).

The original SR

Briefly, SR examines all possible subsets of predictors in a brute force search, fits a regression function on each subset, and evaluates the subset’s stability across environments and its prediction ability. First, the stability of predictor subsets is constructed based on either a Chow test (testing for equal regression coefficients of the predictors between regression functions fitted in a specific environment) or a resampling approach. Subsequently, the prediction ability of stable subsets is evaluated based on negative mean squared prediction error combined with bootstrapping to define a cut-off for selecting the most predictive sets. The importance of each variable with respect to their stability, instability, and prediction ability is then assessed via frequency of selection. The final SR model is obtained as a weighted average of the regression functions fitted on the stable and predictive subsets [refer to (17) for more details].

We identified several limitations of SR in its current form.

• It is computationally infeasible to enumerate every possible subset of predictors in omics data where the number of predictors P is very large (i.e. >30). Pfister et al. (17) proposed the following solution: (i) pre-filter data to tens of predictors. Then from the pre-filtered predictor sets, (ii) randomly sample thousands of subsets to test for stability and subsequently prediction ability. However, we argue that this solution is inefficient, as thousands of subsets are not sufficient to represent the subset space of many predictors. A drastic pre-filtering is therefore required but can result in filtering out important predictors.

• Identifying first the stable predictor sets, and then assessing their prediction ability is not efficient. This is not only because the stable and predictive sets are included in the predictive sets, as we describe in Supplementary Section S3.1, but also because stability is more difficult to compute compared to prediction ability.

Because of these limitations, SR results lack both variable selection and prediction accuracy for large datasets, as we highlight in our simulation (Supplementary Figure S2).

Main steps of StableMate

We summarize the main steps of StableMate; a more detailed algorithm is presented in Algorithm 1 in Supplementary Section S3.1.

• Depending on the type of analysis, a base regressor for ST2* is first specified, for example, we used ordinary least square regression (OLS) for case studies in the ‘StableMate identifies genes associated with ESR1 expression in ER+ breast cancer using RNA-seq data’ and ‘StableMate characterizes cell identity transition of glioblastoma-associated microglia with scRNA-seq data’ sections and generalized linear models in the ‘StableMate discerns global microbial signatures for colon cancer in multi-cohort metagenomics data’ section.

• For each iteration k, k = 1, …, K

  • Apply random Lasso pre-screening, then add pseudo-predictor.

  • Run ST2* to select the most predictive predictor set denoted .

  • Run ST2*to select the stable predictors within such that .

• Define importance score for prediction and stability.

• Calculate significance cut-off scores to define stable, unstable and non-significant variables.

• Fit the final ensemble regression model as the weighted average of the fitted regressions on .

Pre-screening predictors based on random Lasso

We first pre-filter predictors based on a random Lasso procedure. For each ST2* run (described below), we randomly sample one-half of the samples to select the top p predictors with Lasso (34). As an example, we chose p = 100 for the ‘StableMate identifies genes associated with ESR1 expression in ER+ breast cancer using RNA-seq data’ and ‘StableMate characterizes cell identity transition of glioblastoma-associated microglia with scRNA-seq data’ sections. The advantages are of 2-fold. First, across the different resampling runs, the top p predictors are expected to differ, thus enabling us to cover a large and diverse range of predictors in our overall search. Second, we improve the stability of Lasso pre-screening when subjected to sample perturbation (35) (see more details in Supplementary Section S3.1).

ST2*: a new stochastic stepwise variable selection procedure

We improved the ST2 algorithm proposed by Xin and Zhu (36). ST2 is a stochastic version of the classic stepwise variable selection. ST2 searches for a set of predictors maximizing a particular objective function to quantify the predictive ability or stability of predictor sets. It uses a greedy approach with iterative forward and backward searching steps. ST2 starts with an initial predictor set (which can be empty). In the forward step, a collection of predictor sets is randomly sampled from predictors that are not included in the current model. The predictor set that yields the largest increase in the objective function is then added to the current model. In the backward step, a collection of predictor sets is randomly sampled from predictors that are included in the current model. The set that yields the largest increase of the objective function is then removed from the current model. The forward and backward steps alternate until the objective function does not improve further. However, the major drawback of ST2 is that it samples subsets of a fixed size that is randomized at each step. If a wrong size is sampled, ST2 may stop prematurely, leading to inaccurate variable selection.

In ST2*, we follow the ST2 algorithmic framework but we improved the procedure to sample different predictor subset sizes at each step (refer to Supplementary Section S3.1 for a detailed description of the ST2* algorithm). We also added objective functions that are well suited to assess prediction and stability; namely, we used the negative Bayesian information criterion (NBIC) and negative prediction sum of squares (NPSS) (Supplementary Section S3.1).

Finally, we run ST2* on K iterations (e.g. K = 2000 in our case studies) first to identify the predictive subsets of predictors, and second to identify the stable subsets within each predictive subset. As a result, we select an ensemble of stable and predictive predictor sets. These iterations address the stochastic nature of ST2*, which can yield potentially different predictor sets for each iteration. The sets of stable predictors are then used to build the final regression model (described below).

Cut-off prediction and stability scores

We calculate a prediction score for each predictor based on how often the predictor is selected as predictive across the ensembles. We do similarly for the stability score. The output can be represented in a variable selection plot such as Figure 2A, where the scores are represented on the x-axis (prediction) and y-axis (stability).

To define a significance cut-off of these scores, we create a pseudo-predictor as a negative control. A pseudo-predictor is represented as an artificial index P + 1 so that its inclusion in the regression model does not affect the model fitting nor the value of the objective function in ST2*, but it is still taken into account when calculating the scores of all predictor sets.

We applied a bootstrap procedure on the variable selections to compare the distributions of the prediction and stability scores of the predictors to that of the pseudo-predictor to assess their significance. A predictor with a prediction score larger than the pseudo-predictor’s in more than 97.5% times of the bootstrap iterations is considered as significantly predictive. We do similarly for the cut-off stability score. A predictor is considered environment-specific (unstable) if its stability score is lower than that of the pseudo-predictor. Finally, the rest of the predictors are assigned as ‘non-significant’, as shown in plot Figure 2A. Note that since these cut-off scores are based on the bootstrap of variable selections, the significance indicates the variability in the ST2* selections, and hence provide a reference on whether more ST2* runs need to be performed.

Final ensemble regression model generalizable to unseen environments

The final regression model is then built on the different sets of stable and predictive predictors. Each regression model is fitted by regressing the response variable on each stable and predictive subset in the ensemble. We then aggregate these models as the average of the fitted regression weighted by the ranking of objective functions NBIC and NPSS.

Gene modules

PCA (centred but unscaled) was used to identify metagenes in the form of PCs that represent gene modules from the TCGA BC RNA-seq data. The 23 most variable metagenes selected by the elbow method were then used as the predictors of ESR1 expression for the subsequent StableMate analysis (38). To avoid overfitting, ESR1 was removed from the data.

Aggregation of gene expression with similar expression patterns

We applied PCA (not centred nor scaled) on the set of genes, and extracted the loading coefficients of each gene on the first PC using a soft-thresholding approach to identify the top contributing genes with loading coefficients of the same sign. We then considered the absolute value of the loading coefficients of these top genes to obtain positive weights, which we then used for aggregating their expression by a linear combination.

Simulation study

The benchmark results are shown in Supplementary Figures S1 and S2.

We simulated systems of variables observed from four environments. A system is a model that describes the causal relationships between variables, and an environment is the probability distribution of variables that generates data. Therefore, a system of variables in different environments are generated by different probability distributions but with the same causal relationships. The simulations are described in Supplementary Section S3.2. For each simulation run, a variable in the system was randomly sampled as the response while the remaining variables were set as predictors. The regression models were trained on data generated in the first three environments to predict the response. The data of the fourth environment were used for testing.

Metagenomics case study (‘StableMate discerns global microbial signatures for colon cancer in multi-cohort metagenomics data’ section)

We trained each regression model to predict colon cancer disease status. We performed leave-one-dataset-out (LODO) cross-validation. We considered either all 313 species or 431 pathway abundances (no pre-filtering). In addition, we also performed LODO validation to evaluate the performance of the RF models trained using the different sets of predictors selected by either StableMate, Lasso, or RF.

Diffusion map

In case study 3 (see the ‘StableMate characterizes cell identity transition of glioblastoma-associated microglia with scRNA-seq data’ section) we visualized the scRNA-seq data using DM, which is a non-linear dimension reduction method highly suitable for single-cell data with potential cell state transitions. DM learns transition probabilities between cells and projects cells into a lower dimensional Euclidean space that approximates the ‘diffusion distances’ between cells accordingly. DM was run on the 2000 most variable genes of the log normalized data using the R function DiffusionMap with default parameters from the package destiny (43).

Diffusion pseudotime

DPT inference was then applied following DM learning (44). The cell that has the largest DC1 score in the tumour periphery was chosen as the root of the cell trajectory. The distance between the cumulative transition probabilities of any cell with the root cell is defined as its pseudotime.

Sincast imputation

The log-normalized scRNA-seq data was imputed using the sincastImp function with its default parameters, where the imputation of any cell is based on its nearest neighbouring cells.

Sincast projection

Only the most 2000 most variable genes of the query scRNA-seq data were considered. The query data were projected onto the reference PCA atlas after rank transformation. The projection is then reorganized and visualized via DM.

Quantification of query cell identity

To quantify the identity of each query cell based on the reference cell types, we used Sincast modified version of the Capybara cell score of Kong et al. (46). The approach is based on weighted restricted least square regression (33).

Data and StableMate setting.

We used the publicly available BC gene expression (BRCA) data from The Cancer Genome Atlas (TCGA) (48). We filtered the dataset to retain 113 normal and 778 ER+ tumour samples.

Since we were interested in the regulation of ESR1, we set ESR1 as the response and all other genes as predictors and fit simple linear regressions for StableMate analysis. We set the disease status (normal or ER+) as the environment variable so that we could identify stable genes, whose association with ESR1 did not change significantly between normal versus ER+ samples, and disease-specific genes, whose association with ESR1 significantly varied significantly between disease status. In addition to identifying individual genes, we also combined StableMate with PCA to identify stable and disease-specific gene modules. Namely, we still took ESR1 as the response and disease status as the environment variable, but we used meta genes (first a few PCs of all genes except ESR1) as predictors. Then, we defined the stable and disease-specific gene modules as the most important genes of stable and disease-specific metagenes.

StableMate selected genes proxy to ESR1 regulation

The StableMate variable selection results are summarized in Figure 2A. Among the most stable genes predictive of the ESR1 expression, CCDC170 and ARMT1 are the closest genes located to the upstream genomic region of ESR1 (Supplementary Figure S3A) and have been reported to fuse with ESR1 (49). Their proxy to ESR1 suggests that they might be subject to the same transcriptional regulation as ESR1, thus explaining their stability. On the other hand, the STC2 gene was identified as disease-specific. This might be explained by the fact that STC2 has been identified as a downstream target of ER signalling (50,51). In addition, the proximal promotor region of STC2 is not directly subject to ER binding but is dominated by other transcriptional activities such as hypoxia induced stress response (52,53). As a result, ER signalling is indirectly involved in the STC2 activation (51). This evidence supports our hypothesis that STC2 and ESR1 should be indirectly or distally related in transcriptional regulation as indicated by their environment-specific associations.

StableMate with PCA identified gene modules associated with ESR1

Feature selection from transcriptional data is often followed by gene set enrichment analysis. While the stability analysis on individual genes gave us some insights into ER regulation, we selected relatively few genes as either stable or disease-specific—this was insufficient for statistically meaningful enrichment analysis. To overcome this issue, we used the first 23 most significant PCs (selected by the elbow method) (38) of all genes (except ESR1) as predictors for ESR1 expression rather than individual genes. In this context, each PC is a linear combination of the expression of all genes except ESR1, and can be viewed as a metagene, which is useful for quantifying the activities of gene modules (54). Similar to our previous analysis, disease status (normal and ER+) was set as the environmental variable. Of note, we observed similar results with the metagene construction method of weighted gene co-expression network analysis (1), which identifies gene modules with hierarchical clustering first and then computes a metagene on each gene module by eigen decomposition, akin to PCA (see Supplementary Methods S1.1 and Supplementary Figures S4 and S5).

The StableMate variable selection results are summarized in Figure 2B. PC1 and PC2 were found to be highly stable and predictive, suggesting that the major source of variation they explain (15.56% and 10.83%, respectively) is closely related to ER regulation. All subsequent PCs up to PC6 were predictive but disease-specific. We considered the top 200 genes contributing to PC1 (most stable) and PC3 (disease-specific) and conducted an enrichment analysis. Genes from PC1 were mainly associated with biological processes related to hormone regulation (Figure 2C). The ESR1-mediated estrogen signalling is at the centre of hormone regulation, and hence the high prediction ability and stability of PC1 are manifest. Genes from PC3 were associated with basal cell-like transcriptional activities in epidermis development. The top genes contributing to PC3 (see details in Supplementary Figure S3B) included a high proportion of basal cytokeratins (BCKs), such as KRT5, KRT7, KRT14, and KRT17, suggesting that PC3 may reflect the ‘basalness’ of samples (Figure 2C). Interestingly, PC3 scores were positively correlated with ESR1 expression in normal samples but negatively correlated with ESR1 in ER+ samples (Figure 2C). This trend was also observed between the basal BC enriched genes [listed in Li et al. (55)] and ESR1 expression (Supplementary Figure S3C), confirming the PC3 characterization of basalness.

To validate the reproducibility of our findings, we queried the gene expression portals GTEx (22) for normal breast tissue and the METABRIC data from cBioPortal (21) for ER+ BC. Our analysis using these external datasets showed similar trends between ESR1 and PC3 (Figure 2E). The negative correlation between BCKs (contributing to PC3) and ESR1 expression may be explained by the fact that the BCK induction in ER+ BC requires low ER expression (55). However, to the best of our knowledge, no study so far has reported that this correlation may turn positive in normal breast tissues.

Data and StableMate setting

We retrieved eight CRC case-control faecal metagenomic datasets from the R package curatedMetagenomicData (23). The datasets were generated by eight different cohorts from seven countries (refer to Table 1 for the cohort used and for the number of CRC and controls in each cohort). Data were curated into abundance data using a standardized data processing pipeline by Pasolli et al. (23). In total, we collected 604 CRC and 596 control samples. Our analysis focused on the species abundance data measured on 313 microbial species. The analysis of pathway abundance data measured on 431 pathways is detailed in Supplementary Figure S6.

Since we were interested in the CRC diagnosis using metagenomics data, we set the disease status (CRC or normal) as the response and the microbial species as predictors and performed logistic regression in StableMate analysis. We implemented StableMate using the following two strategies. In the first analysis, we applied StableMate as in the toy example (see ‘Illustration of StableMate on toy example’ section) and our first case study (see the ‘StableMate identifies genes associated with ESR1 expression in ER+ breast cancer using RNA-seq data’ section), where all cohorts were pooled to select predictive species and assessed their stability by setting cohort as the environment variable. From this analysis, we found that the majority of the selected predictive species were stable and none of them was cohort-specific (as discussed later in the ‘StableMate discerns global microbial signatures for colon cancer in multi-cohort metagenomics data’ section). Therefore, in a second analysis, we applied StableMate on each cohort to identify cohort-specific predictive species and tested the stability of the species selected in the remaining seven cohorts combined. There we considered only two environments: the specific cohort and the remaining cohorts combined. These ‘cohort-specific analyses’ are useful for identifying species that are highly predictive in a specific cohort but their association with CRC in the specific cohort cannot be generalized to other cohorts.

StableMate identified stable microbial species predictive of colon cancer

From the pooled analysis, we identified 23 stable species to predict disease status (CRC or normal) (Figure 3A). To illustrate the strong cohort effect of the data and the ability of StableMate selection to identify predictors with a mild cohort effect, we used PCoA with all 313 species versus the 23 stable species. The PCoA results combined with a permutation ANOVA showed that the main source of variation was the cohort effect rather than disease status (Figure 3B1) when all species were used. This implies that a predictive model built using all species is likely to be affected by cohort (batch) effects. In contrast, the PCoA and ANOVA results of the 23 stable species selected by StableMate showed a decrease in cohort effects and an increase in the effects of disease status (Figure 3B2). In particular, the CRC and normal samples were better separated in the PCoA when using only the 23 stable species (left panel of Figure 3B2). A formal evaluation of the goodness of variable selection is presented in the ‘Benchmarking StableMate variable selection and prediction on metagenomics data’ section.

StableMate identified cohort-specific microbial species predictive of colon cancer

We conducted cohort-specific analyses for each of the eight cohorts to identify predictors with high cohort specificity. As an example, Figure 3C shows the results for the Austrian cohort. A number of species were found to be highly predictive in the Austrian cohort but with a very low stability score and, therefore, were identified as cohort-specific. Among them, P. copri was the most predictive and one of the most cohort-specific species, suggesting that P. copri might be a marker for CRC specific to the Austrian cohort only. It should be noted, however, that diet could be a confounder, as the Austria-specific species might be related to the low-fibre diet in that population (see Supplementary Section S1.2 for details).

Pooling data improves generalizability of prediction models

Most of the predictive species selected on the pooled data were stable (Figure 3A), whereas predictive species selected in the individual cohorts showed less stability (Supplementary Figure S7). This can be explained as the stability score of a predictor is correlated with the predictor’s influence on the generalizability of regression models. As the number of cohorts increases, regression models tend to put less weight on unstable predictors. As a result, unstable predictors have less influence on model fitting. Hence, we obtained high stability scores for these predictors. This analysis quantitatively confirms a common perception in classical meta-analysis that training regression models on pooled data rather than individual datasets can yield improved generalizability (26,30,58).

However, aside from pooling data, we were able to further improve generalizability of prediction models by taking into account the cohort effect through stability analysis. We conducted a benchmark study in the ‘Benchmarking StableMate variable selection and prediction on metagenomics data’ section, where we showed that the StableMate model built using the stable species outperformed several commonly used regression methods in the pooled analysis.

Data and StableMate setting

From the scRNA-seq dataset from Darmanis et al. (31) of four GBM tumours, we extracted and analysed 1847 myeloid cells from the tumour core (1329 cells) and from the tumour periphery (518 cells) with the cell annotation provided by the authors of the study.

We visualized the scRNA-seq data and observed a clear cell trajectory between the two locations, which may represent the cellular transition of microglia to TAM. We conducted a pseudotime analysis to quantify this trajectory and used StableMate to predict as a response the pseudotime based on expression of the genes as predictors. The base regressor used for StableMate analysis was simple linear regression. The cell location, core or periphery, was set as the environment variable. StableMate selected several cytokines as being predictive of the pseudotime. To further investigate the possible mechanism of these cytokines, we performed a second analysis to build a gene regulatory network for these cytokines. More specifically, we applied StableMate on each of the cytokines as a response and all other genes as predictors. We then summarized the selection results in the form of a network, where the cytokines are connected to their stable predictor genes.

DPT tumour periphery to core

We first visualized the myeloid cells projected onto a bulk RNA-seq reference atlas of myeloid cells (45) to assign the identity of the cells. We performed this using Sincast (33) (Figure 4A1). The projection showed that the cells from the periphery of the tumour closely matched foetal microglia in the reference, and the cells from the core of the tumour matched a wider range of monocytes and macrophages (Supplementary Figure S8A). The projection also showed a continuous state of transition, rather than discrete clusters. We also confirmed the transition by a separate DM analysis in Supplementary Figure S8B. This exploration suggested that the data were suitable for DPT analysis (Figure 4A2), where we set the cells at the tumour periphery as the root (start) of the trajectory (44). The inferred DPT was then used as the response for our first StableMate analysis described below.

StableMate analysis identifies cytokines that signify microglia pre-activation and polarization in tumour periphery

Among the genes selected by StableMate as predictive of DPT, we identified six cytokines whose expression were all negatively correlated with DPT (Supplementary Figure S8C). Among these cytokines, CCL3 and CCL4 were identified as stable (Figure 4B), while TNF, IL1B, CCL2, and CSF1 were identified as periphery-specific (Figure 4C). The selection of the six cytokines is interesting as they are important markers of microglia activation in response to disease (60).

In order to visualize the relationships of these cytokines, we ran StableMate on each cytokine as a response, where all the other genes were used as predictors, and built a gene regulatory network (Figure 4D—DPT was included as a ‘pseudo gene’ here to incorporate the result from the first analysis). This network showed that the two stable cytokines (CCL3 and CCL4) formed a community with DPT, whereas the four periphery-specific cytokines (TNF, IL1B, CCL2, and CSF1) formed another community.

Other stable genes predictive of DPT represented on this network include EGR2 and CD83, which were connected to both DPT and CCL3 (Supplementary Figure S8D). CCL3, CCL4, CD83, and EGR2 are all known to be associated with an immediate early inflammatory response by microglia in a pre-activated state, which is in between homeostasis to those fully activated under pathological conditions (61–66). These four genes showed consistent downregulation during the transition regardless of cell location. On the contrary, the periphery-specific cytokines, which are known markers for microglia polarization to either the pro-inflammatory M1 or anti-inflammatory M2 phenotype (60), exhibited stronger negative association with DPT in the periphery—resulting in low expression levels—and weak association with DPT in the core (Figure 4C).

Core-specific genes revealed reprogramming of tumour-infiltrating microglia into immunosuppressive TAM in GBM tumours

The core-specific genes identified by StableMate included two interesting cell surface markers: the macrophage marker MARCO and MHC class II antigen HLA-DOA (Supplementary Figure S8C). MACRO was lowly expressed in the tumour periphery but upregulated along the DPT trajectory towards the tumour core (Figure 4E, upper). HLA-DOA expression levels had a low–high–low pattern along the trajectory, with high expression levels at the boundary of tumour periphery and core (Figure 4E, lower; other MHC-II genes were also examined in Supplementary Figure S8E). The upregulation of MARCO and the downregulation of HLA-DOA towards the core may indicate the presence of MARCO^hi, MHC-II^lo macrophages, which are characteristic of the M2-like immunosuppressive TAM (67,68). In addition to these two cell surface markers, many pro-tumour markers were also identified by StableMate as core-specific and showed similar expression patterns as MACRO (Supplementary Figure S9A). One example is VCAN, which encodes a large extracellular protein contributing to the establishment of tumour microenvironments (Figure 4F). The expression patterns of these core-specific pro-tumour markers suggest that they responded specifically to the tumour microenvironment and hence are potentially good therapeutic targets.

In addition, we examined the immune activation state of the cells at the beginning of the core stage of the transition. We observed high expression of the stable cytokines CCL3 and CCL4 (Figure 4C), as well as the microglia marker TMEM119, which were all then gradually suppressed in the core (Supplementary Figure S8F). This may imply the reprogramming of activated microglia in the early stages of the core transition to TAMs.

StableMate based on logistic regression outperformed Lasso and performed comparably to RF in CRC classification

The LODO AUC values for all competing methods are shown in boxplots in Figure 5A. To illustrate the benefits of using stable predictors to build regression models, we considered two versions of the StableMate prediction model, one built using all selected predictive variables, the other using only the stable predictors. To further investigate if the differences in AUC values were statistically significant, we conducted a series of two-sided paired t-tests, and the P-values of these tests are shown in Figure 5A.

We first compared the performances of all methods except RF, since they all use variants of linear models to make predictions, whereas RF is a non-linear approach. From Figure 5A, we observed that the StableMate prediction using only stable predictors was significantly better than GLM, Lasso, and StableMate using all predictive variables. Among these methods, GLM had the worst generalizability. As GLM does not perform variable selection, the prediction model potentially included many noisy features. Lasso’s selected variables led to poorer prediction compared to the two variants of StableMate. The version of StableMate using only stable predictors led to slightly higher mean AUC compared to StableMate using all predictive variables, with a P-value indicating a significant difference. The superior performance of StableMate based on stable predictors highlighted the benefits of using such type of predictors to build prediction models.

Finally, we observed that the AUC performance of StableMate based on the stable predictors was indistinguishable from RF. This can be explained as StableMate only uses an ensemble of stringent logistic models to make predictions (see the ‘Materials and methods’ and ‘StableMate to identify stable and environment-specific statistical associations’ sections), whereas RF uses an ensemble of non-linear and highly flexible decision trees targeted for classification tasks.

One can also use StableMate as a variable selection method only and then separately build regression models based on the stable predictors. Filtering out irrelevant and unstable predictors in data can benefit model fitting by improving generalizability (71). Therefore, the notable classification performance of RF motivated the second benchmark study below, in which we evaluated whether StableMate can select predictors that lead to more generalizable RF models.

Species selected by StableMate lead to more generalizable prediction models

We applied StableMate, Lasso, and RF to select three species lists, each containing 23 species. Figure 5B shows a Venn diagram that compares these three lists (see also Supplementary Figure S10A for a more comprehensive comparison). The three lists included an overlap of nine species, all of which are well-known species associated with CRC (70). Among these three methods, StableMate and Lasso shared 20 species, many more than with RF species selection. Of note, StableMate selected three species that were neither selected by Lasso nor by RF (Figure 5B), including Bacteroides intestinalis whose connection with CRC has not been widely reported. Bacteroides intestinalis, along with its metabolic product, is considered as an important biomarker for human gastrointestinal health (72). Under carcinogenetic stress, B. intestinalis has been shown to exhibit increased activity in enhancing DNA integrity maintenance and suppressing central metabolic activities, suggesting B. intestinalis’ protective role against CRC (73,74). This hypothesis is concordant with our observation of a decreased abundance of B. intestinalis in CRC samples across cohorts (Supplementary Figure S11A).

For a quantitative evaluation, we built RF models using six different selections of species and computed their AUC using the LODO cross-validation approach. The results are summarized in Figure 5C, where the 23 stable species selected by StableMate led to the highest mean AUC (0.818). StableMate and Lasso selections had high AUC values since their selections were similar. The difference between StableMate and RF selections was not statistically significant, probably due to a lack of cohorts and statistical power. However, the StableMate selection led to less variable prediction performances (smaller interquartile range) compared to RF. Two species, Dialister pneumosintes and Parvimonas micra, were selected by Lasso and RF but not by StableMate (Figure 5B). In particular, P. micra is known to be associated with CRC as it promotes tumourigenesis (75,76). However, StableMate identified these two species as predictive but not significantly stable. The fact that there was no improvement in the prediction performance of RF trained on the StableMate selection with these two additional species (comparing ‘Lasso | SM’ and ‘SM’ in Figure 5C) justifies why StableMate did not select these two species.

Benchmark based on pathway abundance data

A similar benchmarking analysis based on the pathway abundance data showed that StableMate outperformed the other methods, including RF, in predicting CRC (highest mean AUC; see Supplementary Figure S6B). However, the pathway abundance data were less stable and less predictive of CRC compared to the species abundance data. All methods obtained lower AUC scores in LODO assessment. We observed strong differences in variable selections between the cohorts and the methods (Supplementary Figures S6A and S10).