In this study, the microbiota profiles of samples from 10,362 cases (including 9,687 tumor tissues and 675 adjacent normal tissues across 32 cancer types) at the phylum, class, order, family, and genus levels were obtained from the BIC database (http://bic.jhlab.tw/ ); the corresponding genomic and clinical data for patients were obtained from the TCGA dataset (https://portal.gdc.cancer.gov/repository ). Figure 1 displays an overview of the study design, and the information on included tumor types is summarized in Table S1 .

First, we calculated the relative abundance of microbiota between the tumor and normal groups at the phylum, class, order, family, and genus taxonomic levels. Then, we applied Partial least squares discrimination analysis (PLS-DA) to visualize the microbiota profile landscapes between tumor and normal samples by using the package “mixOmics”. Next, we performed a statistical analysis of metagenomic profiles (STAMP) to investigate the overall differences in microbiota profiles between the tumor and normal groups. Moreover, Linear discriminant analysis effect size (LEfSe) analysis was used to identify specific enriched microbial biomarkers for each group. Linear discriminant analysis (LDA) was further applied to evaluate the microbial effects for different groups.

We first identified the cross-tumor abundant microbiota taxa at the phylum, class, order, family, and genus levels, respectively. Then, we calculated the relative abundance of microbiota across 31 tumor types (except for GBM, which only had normal tissues in the BIC program) at the phylum, class, order, family, and genus taxonomic levels.

Patients with survival data and follow time ≥ 30 days were chosen for survival analyses. First, we investigated the prognostic significance of microbiota by performing a univariate Cox proportional hazards model. Then, NMF was applied to identify distinct microbiome-based clusters based on the abundance of candidate prognostic microbes. The optimal number of clusters and their stability were determined by the consensus clustering algorithm. The R package “NMF” was used to perform the consensus clustering (32).

To quantify the microbiome features of individual patients, we explored a novel microbial-based scoring system (MS) to investigate the microbiota features of the individual patients. Specifically, candidate prognostic microbiomes were chosen from the Lasso regression algorithm to construct the microbial scoring system. The MS was calculated by the corresponding coefficients of selected microbiota signatures:

MS = Σi Coefficient (microbiota)* Abundance (microbiota)

Where i represents the selected microbial signatures.

GSVA enrichment analysis was performed to investigate the variation in biological processes between different MS subgroups by using “GSVA” R packages (33). The gene sets of “c2.all.v2022.1.Hs.symbols” and “h.all.v2022.1.Hs.symbols” were downloaded from the MSigDB database. Adjusted P value <0.05 was considered statistically significant. GSEA was used to explore the signaling enrichment between different MS subgroups by applying the “Clusterprofile” R package (34). The FDR q < 0.25 and P < 0.05 were considered statistically significant.

The ssGSEA algorithm was used to quantify the relative abundance of 29 immune cell types in the TME (35, 36). The relative abundance of each immune cell type in each sample was represented by the enrichment scores that were calculated by the ssGSEA analysis. The CIBERSORT algorithm was applied to analyze the compositions of 22 types of tumor-infiltrating immune cells among different MS subgroups (37).

Using the R package maftools (38), the overall mutation landscape was summarized in patients with high and low MS subgroups in the TCGA cohort. Then, TMB scores based on the TGCA somatic mutation data were calculated to evaluate the mutation status between different MS subgroups.

Student’s t-tests were applied to analyze normally distributed variables and the Wilcoxon rank-sum test was performed to evaluate non-normally distributed variables. One-way ANOVA and Kruskal-Wallis tests were used to conduct difference comparisons of more than two groups. Kaplan-Meier survival analysis and Cox proportional hazards model were chosen to investigate the prognostic significance of microbiota and microbiota-based subtypes by applying the survival and survminer packages. All statistical P values were two sided, with p < 0.05 being statistically significant. All data processing was done in R 4.0.5 software.

Overall, we collected and integrated the microbiota profile and clinical characteristics of 10,362 samples from 32 cancer types. A total of 47, 56, 127, 303, and 1,607 microbial taxa were obtained for each sample at the phylum, class, order, family, and genus levels, respectively. First, the average relative abundance of differential microbiota between tumor and normal tissues was calculated ( Figures 2A–E ), which shows the phylum, class, order, family, and genus levels, respectively. Then, the PLS-DA plot exhibited the microbiota profile landscapes between the tumor and normal samples ( Figure 2F ). We next explored the differential microbial compositions for each group. Overall, 8, 10, 17, 25, and 66 differential microbial components were found between tumor and normal tissues at the phylum, class, order, family, and genus levels, respectively ( Figures 2G ).

At the phylum level, Chordata, Firmicutes, Fusobacteria, Bacteroidetes, and Planctomycetes were the main bacteria groups in tumor samples, while Proteobacteria and Actinobacteria were the main groups in normal tissues. At the class level, the microbial composition of Actinobacteria, Alphaproteobacteria, and Betaproteobacteria were significantly abundant in normal tissues, whereas Fusobacteriia, Clostridia, Bacilli, and Bacteroidia were increased in tumor samples. Other differential microbial composition signatures between tumor and normal tissues are shown in Figure 2G and Table S2 . The top 30 differential microbial taxa between tumor and normal tissues are shown in Figure 2H .

The LDA score of specific microbial taxa in the tumor and normal groups showed the compositional abundances of p_Bacteroidetes, c_Bacteroidetes, o_Nostocalles, f_Prevotellaceae, and g_Ruminococcuss were higher in tumor tissues, while the compositional abundances of p_Actinobacteria, c_Actinobacteria, o_Actinomycetales, f_Moraxellaceae, and g_Acinetobacter were enriched in normal tissues ( Figure 3A and Table S3 ). Cladogram further showed the distinct tumor- and normal-enriched microbiota ( Figure 3B ). Next, we analyzed the cross-tissue abundant microbiome in 32 cancer types at the phylum, class, order, family, and genus levels, respectively. The average relative abundance of microbiota among 31 cancer types (except for GBM, which only had normal tissues in the BIC program) was further explored, which shown in Figures 3C–G . Moreover, the top 10 most differential microbial taxa across 31 cancer types were shown in Figures 3H–L at the phylum, class, order, family, and genus levels, respectively.

We first investigated specific candidate microbial taxa that correlated with survival outcomes by integrating the microbial abundance profile from BIC and survival information from TCGA. The results showed that a total of 182 microbiota were significantly associated with overall survival (OS) ( Table S4 ) and 112 microbiomes were significantly associated with disease-specific survival (DSS) ( Table S7 ). Then, Consensus Clustering analysis of the NMF algorithm was applied to classify patients with qualitatively different clusters based on the abundance of 182 prognostic microbiota ( Figure 4A ). Three distinct microbiome-based clusters were eventually identified, including 4,592 cases in cluster 1, 1,223 cases in cluster 2, and 3,460 cases in cluster 3. A tSNE plot further showed the performance of microbiome-based clusters to distinguish tumor samples in the TCGA cohort ( Figure 4B ). A Kaplan-Meier plot revealed a particularly prominent survival advantage in cluster 2, whereas the worst prognosis was found in cluster 3 ( Figures 4C, D ). The proportion of distinct clusters in 31 tumor types are shown in Figure 4J .

The Lasso regression algorithm based on the 182 prognostic microbiota was performed to find candidate microbial signatures ( Figures 4E–G ). A total of 64 selected microbial taxa were identified from the Lasso regression algorithm ( Table S5 ). Then, a novel microbiome-based scoring system was constructed to quantify the microbial profiles of individual patients, which we termed as microbial score (MS). Consistent with the NMF clustering analysis, two distinct MS subgroups were found and we named these two subgroups MS-low and -high. Further survival analysis indicated significant prognostic differences between the low- and high-MS subgroups ( Figures 4H, I ). The performance of 64 selected microbial signatures to classify different MS groups was further analyzed by using PLS-LA ( Figure 4K ). To better illustrate the association between the established MS with prognosis, an alluvial diagram was applied to visualize the attribute changes of individual patients ( Figure 5A ). LEfSe analysis was further performed to identify specific enriched microbial biomarkers for high- and low-MS groups ( Figures 5B, C , and Table S6 ).

The CIBERSORT algorithm was used to show the differences in the compositions of tumor microenvironment (TME) immune cell types between distinct MS subgroups. As shown in Figure 5D , remarkable differences in immune cell infiltration were observed between the high and low-MS groups, which suggested that intratumor microbiota plays an inevitable role in tumor microenvironment immune profiles. Furthermore, the GSVA algorithm showed significant differences in KEGG pathways, cancer hallmarks, and immune profiles among the distinct MS subgroups ( Figures 5E–G ). The GSEA algorithm showed nucleotide excision repair, cell cycle, DNA replication, homologous recombination, cytosolic DNA sensing pathway, pyrimidine metabolism, proteasome, spliceosome, and P53 signaling pathway were significantly enriched in the MS-high subgroup ( Figures 6A–I ).

We further analyzed the distribution differences of somatic mutation and TMB between low- and high-MS subgroups in the TCGA cohort. The top 20 genes of mutation frequency of the low- and high-MS subtypes are shown in Figure 7A, B , respectively. The top differential mutated genes between the low- and high-MS subtypes are shown Figure 7C . Moreover, a remarkable diversity of tumor mutation burden (TMB) was found between the distinct subgroups ( Figure 7D ). The abovementioned results indicated the potentially complex interaction between genomic variation and microbial components in cancers, which might be novel regulators for tumorigenesis and progression.