Enlibio,Your biological research partner
Xinqidi Biotech Co.,Ltd,Wuhan,China 2008-2022
R&D 14th year

Predicting cancer prognosis and drug response from the tumor microbiome

Issuing time:2022-05-31 16:11

Abstract

Tumor gene expression is predictive of patient prognosis in some cancers. However, RNA-seq and whole genome sequencing data contain not only reads from host tumor and normal tissue, but also reads from the tumor microbiome, which can be used to infer the microbial abundances in each tumor. Here, we show that tumor microbial abundances, alone or in combination with tumor gene expression, can predict cancer prognosis and drug response to some extent—microbial abundances are significantly less predictive of prognosis than gene expression, although similarly as predictive of drug response, but in mostly different cancer-drug combinations. Thus, it appears possible to leverage existing sequencing technology, or develop new protocols, to obtain more non-redundant information about prognosis and drug response from RNA-seq and whole genome sequencing experiments than could be obtained from tumor gene expression or genomic data alone.

Introduction

The Cancer Genome Atlas (TCGA), available from the NCI Genomic Data Commons (GDC)1, provides RNA-seq and whole genomic sequencing (WGS) data for thousands of cases across dozens of cancer types. RNA-seq data is typically used to measure the expression of human genes, and there is a long history linking tumor gene expression to cancer outcomes2,3,4,5,6,7,8. Milanez-Almeida et al.9 recently showed that gene expression from TCGA RNA-seq data could predict overall survival (OS) or progression-free interval (PFI) better than classical clinical prognostic covariates—age at diagnosis, gender, and tumor stage. Importantly, Milanez-Almeida et al. used a data-driven machine learning (ML) based approach which selected features that were predictive of and correlated with prognosis, rather than approaches based on classical statistics or biological knowledge that chose features a priori.

Research into the human tumor microbiome has been rapidly expanding, and multiple laboratories have attempted to utilize existing technologies and data to identify microbes and quantify their abundance within human tumors compared to adjacent normal tissue. RNA-seq and WGS data not only contain human sequencing reads, but also reads from the local intratumor microbiome that are typically filtered out from the data when analyzing human gene expression or genomic alterations. Poore et al.10 recently developed a computational workflow, using two orthogonal microbial detection pipelines, to estimate, decontaminate, normalize, and batch effect correct microbial abundances from human high-throughput sequencing data. They applied this workflow to create a comprehensive dataset of pan-cancer tumor microbial abundances derived from WGS or RNA-seq data for the entire TCGA cohort.

Our central research questions then were, (1) does a data-driven ML approach reveal that tumor microbial abundances in TCGA data, quantified from these reads, are predictive of cancer prognosis or drug response, (2) what microbial genera are potentially predictive biomarkers of prognosis or drug response, (3) how do these models compare to equivalent models based on tumor gene expression data, and (4) does combining both microbial abundance and gene expression features produce models and select combinations of genes and microbial genera that are more predictive of prognosis or drug response than models from each individual data type? Here we use the processed microbial abundances directly from the Poore et al. dataset to build predictive models of prognosis and drug response for TCGA. We also use TCGA RNA-seq read counts to build equivalent predictive models for comparison.

We show that in four cancer types, adrenocortical carcinoma, cervical squamous cell carcinoma, brain lower grade glioma, and subcutaneous skin melanoma, tumor microbial abundances are better predictors of prognosis than clinical covariates alone. However, we find that tumor gene expression is a more powerful predictor of prognosis, across a wider range of cancer types, than microbial abundances. Moreover, we find five cancer-drug pairs where tumor microbial abundances are more predictive of patient drug response than clinical covariates alone. These five pairs include docetaxel treatment for breast invasive carcinoma and sarcoma, and several treatments for stomach adenocarcinoma. We find that tumor microbial abundances are similarly as predictive of drug response as gene expression, but in mostly different cancer-drug combinations.

Results

Tumor microbial abundances are substantially less predictive of prognosis than gene expression

An overview of the analytical workflow is presented in Fig. 1. It has four major parts, (1) data download and preprocessing, (2) prognosis and drug response ML modeling, (3) model evaluation and scoring, and (4) further feature analysis. A more detailed technical description of the analysis pipeline and computational methods is provided in Methods.

Fig. 1: Analysis pipeline overview.
figure 1

Download and data preprocessing (left) of Poore et al.10 TCGA primary tumor Kraken2 Voom-SNM microbial abundances with additional filters to reduce technical variation, NCI Genomic Data Commons (GDC) harmonized TCGA primary tumor RNA-seq counts and clinical data, TCGA curated overall survival (OS) and progression-free interval (PFI) outcome data45, and TCGA curated drug response clinical data14. Prognosis machine learning (ML) modeling (middle) of microbial abundance, gene expression, and combined data types with clinical covariates for each cancer using penalized Cox with elastic net penalties (Coxnet) against matched clinical covariate-only models using standard Cox regression. Drug response classification ML modeling of the same data types with clinical covariates for each cancer-drug combination using three ML approaches, (1) SVM-RFE, elastic net logistic regression (LGR), and limma-trend (microbial and combined data types) or edgeR (gene expression) differential analysis feature scoring and selection with L2 penalized LGR. Matched clinical covariate-only modeling performed with L2 penalized linear SVM or LGR. ML modeling generates 100 model instances for each model from 75/25 train/test randomly shuffled and stratified dataset splits. ML model instance scoring (right top) using concordance index (C-index) and time-dependent cumulative/dynamic AUC (C/D AUC(t)) for prognosis models and area under receiver-operating characteristic curve (AUROC) and area under the precision-recall curve (AUPRC) for drug response models. Significance of model performance improvement over matched clinical covariate-only model determined by signed rank test of C-index or AUROC scores between each matched model instance for prognosis and drug response models, respectively. Feature analysis (right bottom) performed using model instance coefficients and selection frequencies. Overall feature importance ranking and significance determined by signed rank test of model instance feature coefficients shifting from zero and filtering of top features for selection frequency ≥ 20%.

We built OS and PFI gene expression ML models of 32 TCGA tumor types (see Supplementary Table 1 for cohort information) using the Coxnet11 algorithm, which jointly selects the most predictive subset of features via cross-validation (CV) while simultaneously being able to control for prognostic clinical covariates. In our models, we included and controlled for the clinical covariates age at diagnosis, gender, and tumor stage. For comparison, we also built standard Cox regression models based on the clinical covariates alone. We evaluated the predictive performance of our models using Harrell’s concordance index (C-index), which is a metric of survival model predictive accuracy. Each model analysis generated 100 model instances and C-index scores from randomly shuffled train-test CV splits on the data. We found 33 OS and PFI models for 21 tumor types that had a mean C-index score ≥0.6 and significantly outperformed their corresponding clinical covariate-only models (Fig. 2a, c, Supplementary Figs. 1a, 2a). Our models were predictive of prognosis in 11 of the same 13 tumor types that were reported by Milanez-Almeida et al.9 (Supplementary Table 2). We did not analyze one tumor type that Milanez-Almeida did, acute myeloid leukemia (LAML), because Poore et al. excluded it from their analysis. Among the cancers and outcomes that Milanez-Almeida et al. analyzed, our methodology produced predictive models for four additional tumor types: breast cancer (BRCA), cervical squamous cell carcinoma (CESC), sarcoma (SARC), and uterine corpus endometrial carcinoma (UCEC), as well as quite a few predictive models for additional cancers and outcomes that were not analyzed in their study (Supplementary Table 2). We also evaluated prognosis model performance by calculating the time-dependent, cumulative/dynamic area under the curve (AUCC/D(t))12,13, which is an extension of the area under the receiver-operating characteristic curve (AUROC) for continuous outcomes and can provide a more detailed resolution picture of predictive performance throughout the test outcome time range compared to the C-index score. Although 33 of our OS and PFI gene expression models had a statistically significant C-index score improvement compared to clinical covariates alone, only 22 of these models showed an improvement in AUCC/D(t), where the improvement in mean AUCC/D(t) over the entire test time range after diagnosis was ≥ 0.025 (Supplementary Figs. 1b, 2b).

Fig. 2: Performance of gene expression and microbial abundance prognosis prediction models where features add predictive power to clinical covariates (a) gene expression with clinical covariate models (orange) and (b) microbial abundance with clinical covariate models (blue) vs clinical covariate-only models (grey).
figure 2

In both a and b data are presented as mean values +/− standard deviation of the mean (SDM) for \(n=100\) random training/test splits as described in Methods. Significance was computed by a paired two-sided Wilcoxon signed rank test, FDR adjusted for multiple comparisons: * \(p\le 0.01\), ** \(p\le 0.001\), ***\(p\le 0.001\). (c) C-index score violin density plots for \(n=100\)training/test splits for the six models where microbial abundance with clinical covariate features outperform clinical covariate-only models. Box plots within the violin plots show median as center, the lower and upper hinges that correspond to the 25th and the 75th percentile, and whiskers that extend to the smallest and largest value no more than 1.5 times the interquartile range from the median. Corresponding gene expression models shown for comparison. Lines connecting points (light grey) represent score pairs from same train-test split on the data. Mean C-index scores shown as red dots with red lines connecting the means. Significance for the prediction improvement over clinical covariate-only models was calculated using a two-sided Wilcoxon signed-rank test and adjusted for multiple testing using the Benjamini-Hochberg method with adjusted p-values shown at top. These are the same p-values indicated in panel a. Adjusted p-values colored in red signify difference where clinical covariate-only model is better. Source data and exact p values are provided as a Source Data file. The number of cases involved in each experiment are shown in Supplementary Table 1.

We applied Coxnet11 using the same methodology to build prognosis models using the microbial abundance estimates provided by Poore et al.10 We found six microbial abundance models that had a mean C-index score ≥0.6 and significantly outperformed their corresponding clinical covariate-only models (Fig. 2b, c, Supplementary Fig. 3a). We found that in only two of the six models, microbial abundances outperformed clinical covariates alone in terms of AUCC/D(t), where the improvement in mean AUCC/D(t) over the entire test time range after diagnosis was ≥0.025 (Supplementary Fig. 3b). In adrenocortical carcinoma (ACC), microbial features predicted OS significantly better than clinical prognostic covariates starting at approximately 6 years after diagnosis. In CESC, microbial abundances predicted OS better than clinical covariates from approximately 6 months to 10 years after diagnosis. Overall, we found that tumor microbial abundances from Poore et al. were only marginally predictive of prognosis across the TCGA cohort, and that gene expression was a significantly more powerful predictor of prognosis (Fig. 2, Supplementary Figs. 13).

Tumor microbial abundances are predictive of chemotherapy drug response in some cancers and in mostly different cancer-drug combinations than gene expression

We next asked whether tumor microbial abundances from pre-treatment biopsies could predict drug response better than the clinical covariates age at diagnosis, gender, and tumor stage alone. TCGA drug response clinical data were obtained from Ding et al.14 as described in Methods. Cases with complete response (CR) or partial response (PR) were labeled as responders and those with stable disease (SD) or progressive disease (PD) as non-responders. Thirty TCGA cancer-drug combinations met our minimum dataset size thresholds (see Supplementary Table 1 for cohort sizes and Supplementary Data 1 for a more detailed breakdown). We built drug response models using three different ML methods: (1) a variant of the linear support vector machine recursive feature elimination (SVM-RFE) algorithm15that we developed, (2) logistic regression (LGR) with elastic net16 (L1 + L2) penalties and embedded feature selection, and 3) logistic regression with an L2 penalty and limma17 (for microbial abundance and combined data type datasets) or edgeR18,19 (for RNA-seq count datasets) differential abundance/expression feature scoring and wrapper selection methods (see Methods for details). All three ML methods unconditionally included the clinical covariates—age at diagnosis, gender, and tumor stage – in the model (bypassing feature selection) while selecting the most predictive subset of microbial abundance or gene expression features. For comparison, we built standard linear SVM or LGR models using the clinical covariates alone. We evaluated the predictive performance of drug response models using AUROC. Each analysis generated 100 model instances, AUROC, and area under the precision-recall curve (AUPRC) scores from randomly shuffled train-test CV splits on the data.

We found five microbial abundance cancer-drug combinations that had a mean AUROC score ≥0.6 and performed better than clinical covariates alone in at least two out of three ML methods (Fig. 3). Three of these cancer-drug combinations involved stomach adenocarcinoma (STAD). We performed the same drug response modeling using TCGA gene expression data and here we found six cancer-drug response combinations that had a mean AUROC score ≥0.6 and significantly outperformed their corresponding clinical covariate-only models in at least two out of three ML methods (Fig. 4). Only one cancer-drug combination, SARC docetaxel, overlapped between the microbial abundance and gene expression drug response model results, suggesting that tumor microbial abundances have independent predictive power. Even though one of our thresholds for a significant drug response model was a mean AUROC score ≥0.6, the 11 total significant models that we found from both data types each had a mean AUROC > 0.7. We also found there was considerable overlap in the selected microbial abundance and gene expression features reported by each ML method (Fig. 5a, c) and frequently found a significant correlation between the feature importance rankings reported by each ML method when comparing the two most significant methods in each cancer-drug combination (Fig. 5b, d). These results suggest that our significant drug response models and their inferred important features are not the results a specific ML modeling methodology. Overall, our results support the notion that the tumor microbiome may contain information that is predictive of drug response in some cancers, consistent with recent reports20,21.

Fig. 3: Performance of microbial abundance drug response prediction models in the five cancer-drug combinations where models performed better than clinical covariates alone.
figure 3

a Mean AUROC scores for microbial abundance with clinical covariate models (blue) vs clinical covariate-only models (grey) Significance computed by a paired two-sided Wilcoxon signed-rank test, FDR adjusted for multiple comparisons: * \(p\le 0.01\), ** \(p\le 0.001\), ***\(p\le 0.001\). b mean AUROC scores for each ML method for pairs presented in a. In both a and b data are presented as mean values +/− SDM for \(n=100\) random training/test splits as described in Methods. c Violin density plots of AUROC scores for microbial abundance with clinical covariate models vs clinical covariate-only models for \(n=100\) training/test splits. Box plots within the violin plots show median as center, the lower and upper hinges that correspond to the 25th and the 75th percentile, and whiskers that extend to the smallest and largest value no more than 1.5 times the interquartile range from the median. Lines connecting points (light grey) represent score pairs from same train-test split on the data. Mean AUROC scores are shown as red dots connected by red lines. d Mean ROC (blue) and e precision-recall (PR) curves (purple) for microbial abundance with clinical covariate models vs clinical covariate-only models (grey). Mean AUROC and AUPRC scores shown in legends and shaded areas denote standard deviations. Significance for the prediction improvement over clinical covariate-only models was calculated using a paired two-sided Wilcoxon signed-rank test and adjusted for multiple testing using the Benjamini-Hochberg method with adjusted p-values shown at top of violin plots in c that are the same as the p-values indicated in panels a and b. In ce results for the modeling method that had the most significant Wilcoxon signed-rank test are shown. Source data are provided as a Source Data file. The number of cases involved in each experiment are shown in Supplementary Table 1.

Fig. 4: Performance of gene expression drug response prediction models in the six cancer-drug combinations where models performed better than clinical covariates alone.
figure 4

a Mean AUROC scores for gene expression with clinical covariate models (orange) vs clinical covariate-only models (grey) Significance was computed by a paired two-sided Wilcoxon signed-rank test, FDR adjusted for multiple comparisons: \(p\le\) 0.01, ** \(p\le\) 0.001, *** \(p\le\) 0.0001. b Mean AUROC scores for each ML method. In both a and b data are presented as mean values +/− SDM for \(n=100\) random training/test splits as described in Methods. c Violin density plots of AUROC scores for gene expression with clinical covariate models vs clinical covariate-only models for \(n=100\) training/test splits. Lines connecting points (light grey) represent score pairs from same train-test split on the data. Box plots within the violin plots show median as center, the lower and upper hinges that correspond to the 25th and the 75th percentile, and whiskers that extend to the smallest and largest value no more than 1.5 times the interquartile range from the median. Mean AUROC scores are shown as red dots connected by red lines. d Mean ROC (orange) and e precision-recall (PR) curves (green) for gene expression with clinical covariate models vs clinical covariate-only models (grey). Mean AUROC and AUPRC scores shown in legends and shaded areas denote standard deviations. Significance for the prediction improvement over clinical covariate-only models was calculated using a paired two-sided Wilcoxon signed-rank test and adjusted for multiple testing using the Benjamini-Hochberg method with adjusted p-values shown at top of violin plots in c that are the same as the p-values indicated in panel a. In ce results for the modeling method that had the most significant Wilcoxon signed-rank test are shown. Source data are provided as a Source Data file. The number of cases involved in each experiment are shown in Supplementary Table 1.

Fig. 5: Comparison of drug response model top-ranked selected features by each ML method. For each drug response model, we selected the two best ML methods by significance for the prediction improvement over their respective clinical covariate-only model.
figure 5

Venn diagrams for microbial abundance (a) or gene expression (c) models comparing the number of features individually selected by each ML method, and the intersection of the two ML methods. Spearman rank correlation plots for microbial abundance (b) or gene expression (d) models showing that the median rank of features (among the 100 model instances in which the feature was selected) often correlated between the two most significant ML methods; p-values are two-sided. The best method is shown on the x-axis, the second best on the y-axis. Source data are provided as a Source Data file. The number of cases involved in each experiment are shown in Supplementary Table 1.

Combining tumor microbial abundance and gene expression features adds a modest predictive improvement in some cancers

Finally, we investigated if models built from combining microbial abundance and gene expression features would result in an improvement in predictive power over their corresponding single data type models. Combining data types resulted in a modest predictive improvement in only three prognosis models: SARC OS, STAD PFI, and thymoma (THYM) OS (Supplementary Fig. 4a). Although this improvement was not statistically significant in terms of C-index score, the AUCC/D(t) metric showed a clear improvement in prognostic predictive power for these models, where the improvement in mean AUCC/D(t) over the entire time range after diagnosis was ≥0.025 compared to their respective single data type models. We also found five combined data type drug response models which performed significantly better than clinical covariates alone, although none of these models reached statistical significance when compared to their respective single data type models in terms of improvement in AUROC score, but one of these models, for BLCA cisplatin, did show an improvement in AUROC ≥0.025 compared to its corresponding single data type models (Supplementary Fig. 4b, c).

Evaluating the robustness of drug response models

Some of the TCGA drug response cohorts used in our study were of limited size and this could have an impact on the robustness of our analysis (see Supplementary Table 1 for cohort sizes and Supplementary Data 1 for a more detailed breakdown). To study this issue further, we evaluated the significance of model scores using a class label permutation test. We shuffled dataset class labels 1000 times and each time ran the outer CV procedure on the permuted dataset, where for each CV iteration we fit a model instance and calculated an AUROC score. We then calculated a p-value from the fraction of permuted scores that were greater than or equal to the true score. Three of the five microbial abundance drug response models that were reported above to have performed significantly better than clinical covariates alone had a permutation test p-value <0.05 and the remaining two, for stomach adenocarcinoma (STAD) cisplatin and oxaliplatin, had p-values <0.08 (Fig. 6a). Permutation test scores and significance for microbial abundance models were similar regardless of the modeling method used (Supplementary Fig. 5a). Five of the six gene expression drug response models that performed significantly better than clinical covariates alone had a permutation test p-value <0.05 (Fig. 6c). Again here, permutation test scores and significance were similar regardless of the modeling method used (Supplementary Fig. 6a). The testicular germ cell tumor (TGCT) bleomycin gene expression model did not quite reach significance, though it is worth mentioning that for the edgeR feature selection and L2 logistic regression modeling method it was close (p = 0.077).

Fig. 6: Evaluation of drug response model robustness. Model significance and robustness was further evaluated using a class label permutation test and examination of the effect feature selection had on model performance. Results for the modeling method which had the most significant Wilcoxon signed-rank test are shown.
figure 6

Permutation test result histograms and significance for microbial abundance (a) or gene expression (c) models showing the distribution of permutation mean AUROC scores. True mean AUROC score shown as dotted vertical grey line and kernel density estimate shown as a curve over the histogram. Curves showing the effect that model hyperparameters which control the number of selected features had on mean AUROC and average precision (AVPRE) scores during hyperparameter grid search across all 100 model instances for microbial abundance (b) or gene expression (d) models. Shaded areas denote standard deviations. Source data are provided as a Source Data file. The number of cases involved in each experiment are shown in Supplementary Table 1.

We further evaluated the robustness of our significant drug response models by examining the effect that the number of selected features had on model performance. During the hyperparameter grid search and tuning that occurred in the nested inner CV during each model instance fitting, scores for every combination of hyperparameter setting and inner CV train/validation fold were saved (see Methods for full details). We plotted how these scores were affected by the hyperparameters that controlled feature selection. Our decision to conservatively limit the feature selection search space in our drug response models to a maximum of 400 best scoring features, to reduce model complexity and the possibility of overfitting, appeared sufficient, as scores for our significant models reached a maximum or leveled off well within this search range (Fig. 6b & d). In the five microbial abundance models, predictive power was driven by a small number of features in three models, where selecting more features did not contribute to additional predictive power or it added noise (Fig. 6b). Even in the remaining two models, most of the predictive power was driven by the top 50 to 100 features. In the six gene expression models, this finding was even more stark, where all the predictive power was achieved by a small number of features in each model (Fig. 6d). In all the significant models from both data types, the variance in scores was not significantly affected by the number of selected features and feature-to-sample ratio within our chosen hyperparameter search range. As with the permutation test results, we found the effect that the number of selected features had on model performance was similar regardless of the feature selection or modeling method used (Supplementary Fig. 5b, 6b). In summary, these two comprehensive analyses suggest that the significant cancer-drug response combinations found in this study and the most important features inferred from their models represent a potentially real and robust biological signal.

Feature analysis reveals a wide range of predictive microbial genera

To learn more about the most predictive features, we determined the top microbial genera and top genes (Supplementary Data 2) selected by each of the significantly predictive microbial abundance and gene expression models, respectively, according to their selection frequency and model coefficients across the 100 model instances from each analysis. There were 428 distinct microbial genera appearing in at least one prognosis or drug response model. Of these 428 genera, 160 were individually significantly predictive of prognosis or drug response by a Wilcoxon test, indicating that the other genera were significantly predictive in combination. The median number of genera selected per model was 52, with a minimum of 3 (BRCA docetaxel) and a maximum of 78 (STAD cisplatin). Of the 428 genera, 95 were selected in more than one model and only 13 were selected in more than two models. This is consistent with the observation of Nejman et al.22 that the tumor microbiome is tumor type-specific. The predictive genera we found span all non-eukaryotic domains of life, in total encompassing 365 bacterial, 17 archaeal, and 46 viral genera (Supplementary Data 2).

Discussion

In summary, we find that the microbial abundance estimates generated by Poore et al.10 are predictive of cancer patient prognosis and response to chemotherapy in a subset of tumor types, survival outcomes, and treatments. Machine learning methods, such as those applied in this study, are not able to infer causality, but only inform on the positive or negative predictive associations covariates have with the response variable. The potential causal role that those covariates may play in determining patient prognosis or drug response can only be ascertained via dedicated mechanistic studies. Overall, in terms of the number of significant models, based on their cross-validated C-index or AUROC scores and improvement over clinical covariates alone, the tumor microbiome is considerably less predictive than the tumor human transcriptome at predicting patient prognosis, but notably, performs similarly to gene expression at predicting chemotherapy response and in mostly different cancer-drug combinations. Our investigation motivates future studies investigating the role of the tumor microbiome in predicting the response to targeted therapies and immunotherapies.

There are also some limitations to our current study. As we described previously, some TCGA drug response cohorts were of limited size or had relatively few responder or non-responder cases within these cohorts and this could have an impact on the interpretability of the results. Vabalas et al.23 conducted a literature review of ML algorithm validation of high-dimensional biological data models with limited sample size and performed their own independent simulation analyses evaluating different techniques. They found that, consistent with previous literature, nested CV was the optimal validation method and gives unbiased performance estimates regardless of sample size. They also found that performing feature selection and other model development steps (e.g., normalization, outlier removal) fully within the inner nested CV is essential to avoid overfitting and to produce unbiased results, and that hyperparameter tuning should ideally also be performed in nested fashion. Finally, they found that performing an adequate number of CV folds was important to reduce bias. Our analyses have followed their observations and recommendations, employing them at every level of model development and evaluation, including additional techniques not reviewed in their work (see Methods for full details).

There are further limitations to this study inherited from limitations in the data, originally raised by Poore et al.10 First, the study was retrospective, using existing data from the TCGA. As such, it did not involve any specific protocols to capture microbial reads or to control for contamination. Second, decontamination of such retrospective data is a highly involved and dataset-specific process, which they made great effort to validate. Poore et al. conclude from this validation that the retrospective study of TCGA was successful, and that similar retrospective studies would be valuable. A third point, which they touch on briefly, is that the protocols that were used have limitations with respect to capturing microbial reads and cannot distinguish if the source of microbial reads is intracellular or extracellular, or alive or dead when the sample was taken. Poore et al. suggest, correctly we believe, that additional protocols need to be developed for prospective studies.

Accepting the limitations of the study, we observed certain trends. Proteobacteria and Firmicutes were the most frequent phyla identified as predictive features (Supplementary Data 2), followed by Actinobacteria and Bacteroidetes. Among viruses, Herpesvirales were the most frequent. More microbial genera were negatively predictive of drug response or prognosis than those positively predictive (negative for 306/537 features; two-sided binomial test p-value = 0.0014). Firmicutes reversed this trend, being more often positively predictive (positive for 49/82 features, two-sided Fisher’s exact test p-value = 0. 0.0036; Supplementary Table 3).

Further examining the predictive features of our significant models and their cancer types, we found that several genera of Firmicutes were predictive of OS in CESC, including genera of Lactobacillaleswere found to be negatively predictive of survival. We also found that the genus Chlamydia had an even stronger negatively predictive association with OS in CESC. Notably, though CESC is known to often arise from HPV infection, the presence of other microbial species, in particular the genera Chlamydia and Lactobacillales, have been reported in the literature to be associated with the risk of developing CESC24,25.

Our prognosis analysis results were different than two recent reports26,27 which found some intratumor microbes that were potentially correlated with prognosis in three TCGA cancers. We did not find that the Poore et al. tumor microbial abundances estimated from TCGA were predictive of OS or PFI in these three cancers using our data-driven, regularized ML computational approach. A few important possible reasons for this difference in results are that different source data and methods were used to perform prognosis analysis compared to these studies. Gnanasekar et al.26 analyzed the THCA cohort by tumor subtype, they used harmonized and normalized GDC TCGA data instead of legacy TCGA followed by normalization and batch effect correction as in Poore et al., they only used RNA-seq data instead of WGS and RNA-seq data, they applied different methods for extraction of microbial reads and decontamination, and finally they did not perform any direct analysis of correlation of their derived microbial abundances with survival outcomes. Dohlman et al.27 analyzed colorectal cancers (colon (COAD) and rectum (READ) adenocarcinomas) also using harmonized and normalized GDC TCGA data, they used WGS and whole exome sequencing (WXS) data instead of WGS and RNA-seq, they also used different methods for extraction and decontamination of microbial reads, and finally they also applied classical univariate statistics on their entire data to infer correlation with overall survival (OS). While we believe the use of harmonized GDC TCGA data is superior to legacy TCGA, Poore et al. applied robust computational methods to remove technical variation from legacy TCGA data and validated that their approach was effective. We also applied additional filters of TCGA samples to further remove technical variation. We also believe that, in general, applying classical univariate statistics on the entire data to find correlations has the potential to overfit the specific dataset and it does not consider the multivariate nature of high-dimensional biological data like intratumor microbial abundances. A data-centric, multivariate, and regularized ML approach focused on fitting models on training data and evaluating on unseen test data has potential to generalize better and discover whether features are potentially predictive of and correlated with the response variable, such as survival outcomes or drug response.

Looking at our drug response model results, in STAD, tumor microbial abundances were predictive of response to three different drugs: cisplatin, leucovorin, and oxaliplatin. The genus Helicobacter was a quantified microbial abundance feature in the Poore et al. dataset although notably, even though it is well established that patients infected with H. pylori have an increased risk of developing gastric cancer28, Helicobacter was not identified as a predictive feature of drug response in our STAD models. This finding is in line with recent research indicating reduced microbial diversity, decreased abundance of H. pylori, and enrichment of other mostly commensal bacterial genera in gastric carcinoma29. Instead, in STAD we found that known opportunistic bacteria Cedecea and Sphingobacteriumwere both strongly negatively predictive of leucovorin response, Sphingobacterium was strongly negatively predictive of cisplatin response, and the opportunistic bacteria Rouxiella was strongly negative predictive of oxaliplatin response. Cedecea and Sphingobacterium have been implicated in bacteremia in immunocompromised individuals in rare cases, including cancer30,31,32,33. As dysbiosis is frequent in stomach cancer34,35, and considering the mechanism of action of leucovorin, it may be of interest to study whether organisms from these two genera may sequester or prevent the bacterial production of folinic acid36.

We found three microbial genera whose abundances were strongly associated with breast cancer response to docetaxel. Indeed, the involvement of the tumor microbiome in breast cancer (BRCA)22,37 has recently received considerable attention. In BRCA, we found that the genus containing Epstein-Barr virus (EBV) was negatively associated with response to docetaxel, which is concordant with previous findings that EBV is associated with chemoresistance to docetaxel in gastric cancer38. Interestingly, Cyanobacteria were predictive features in several cancers in our study and we identified a genus of Cyanobacteria as predictive of response to docetaxel in BRCA. Notably, the presence of Cyanobacteria in BRCA was recently confirmed by Nejman et al.22 by 16S-rRNA sequencing. While the genus we identified, Raphidiopsis, a planktonic Cyanobacteria that produces toxins harmful to human health and found in freshwater, is possibly a taxonomic identification error in the original microbial abundance estimates, our findings may point to a related genus under the recently discovered clade Melainabacteria of Cyanobacteria39, which is present in humans. Though Melainabacteria are difficult to culture, we believe that confirmation of the relationship between BRCA to response to docetaxel and Melainabacteria should be tested, and a first step would be to confirm our computationally derived findings in a dedicated 16S-rRNA analysis.

Interestingly, in sarcoma (SARC), among the most predictive microbial features, we found the genus Lactococcus to be positively associated with response to docetaxel. Lactococcus contains species that can sometimes cause opportunistic infections in humans, as Lactococcus are similar to Streptococcus and formerly belonged to that genus. The result that this genus was positively associated with response in our model initially appeared counterintuitive, although while the use of therapeutic bacteria as antitumor agents has not been an extensively studied field, there have been some limited findings in the literature that suggest the use of bacteriotherapy as anticancer agents40. Historically, the intentional use of the toxins of various Streptococcus species showing significant antitumor activity in SARC has been documented41,42,43. One possible testable explanation for some microbes being strongly positive predictive of docetaxel response in our model is that they might produce some extracellular products or toxins that could work as an adjuvant to the chemotherapy.

In summary, while these findings and others reported in this study are computationally derived associations, we believe that they can serve as leads for further experimental studies of the role of microbial species in modulating patient survival and drug response, potentially by metabolizing drug levels in the tumor microenvironment as suggested above, or by altering the immune response, either by changing the levels of specific immunometabolites or by having the tumors present specific bacterial antigens44.


Article classification: Biological abstract
Share to:
Tel:+86-027-87610298
Tel:+86-027-87610297
Add:Room A11-329, 1st Floor, No.1, SBI Venture Street, Optics Valley, East Lake
New Technology Development Zone, Wuhan, China.
Certificate NO.:U18Q28010569R0S