dataesr/scientific-paragraphs-categorization huggingface dataset:
5× the IQR (n = 4). c, Importance of each listed feature (belonging to the four-kingdom model) by the cross-validation of predictive performance for each population dataset as estimated using the internal random forest 'Gini importance' method. The 'aggregation' column shows the integrated ranks (using a rank aggregation algorithm) of listed markers within each cohort along with changes in abundance (differentials), with red indicating a species increase and blue indicating a species decrease in patients with CRC compared to controls. d, AUROC matrix of models built with the panel of 16 multi-kingdom features for CRC detection. Values on the diagonal refer to the average AUROC of 20× repeated fivefold stratified cross-validations. Values off the diagonal refer to the AUROCs obtained by training the model on the population of the corresponding row and applying it to the population of the corresponding column.
The LOCO row refers to the performances obtained by training the model on the 16 microbial features using all but the population dataset of the corresponding column and applying it to the dataset of the corresponding column. The asterisk represents the significance of models assessed with 1,000 permutations (two-sided test). *P = 0. 001. e, AUROC matrix of models built with the panel of 16 multi-kingdom features for CRC early detection. Values on the diagonal refer to the average AUROC of 20× repeated fivefold stratified cross-validations. Values off the diagonal refer to the AUROCs obtained by training the model on the population of the corresponding row and applying it to the population of the corresponding column. The LOCO row refers to the performances obtained by training the model on the 16 microbial features using all but the population dataset of the corresponding column and applying it to the dataset of the corresponding column.
Validation of the 16-marker multi-kingdom panel in independent cohorts. To externally validate the predictive performance of our 16-marker multi-kingdom panel and avoid overoptimistic reporting of model accuracy, we analysed 3 independent datasets (Supplementary Data 1 and 2) from China (CHN_SH, 86 controls and 80 patients with CRC), Italy (ITA, 52 controls and 61 patients with CRC) and the United States (USA, 52 controls and 52 patients with CRC) (Extended Data Fig. 7a ). The average AUROC of the cross-validation models was 0. 88 for the CHN_SH cohort and 0. 81 for the ITA cohort, respectively, while the AUROC was relatively lower for the USA (0. 68 on average). The latter may be related to long-time frozen storage of samples (over 25 years). The average AUROC for the cohort-to-cohort analysis was relatively decreased (CHN_SH:0.
75, ITA:0. 71 and USA:0. 66; Extended Data Fig. 7b ), while the AUROCs for the LOCO analysis were slightly improved, ranging from 0. 70 to 0. 76 (Extended Data Fig. 7c ). Altogether, this additional cohort analysis validated the robustness of our multi-kingdom marker panel across a total of eight cohorts from seven countries. Specificity of the CRC predictive models based on the multi-kingdom marker panel. In light of shared microbiota alterations across various diseases  28  , it is necessary to verify the disease specificity for the identified microbial biomarkers panel, thereby ensuring a low false positive rate for CRC diagnosis. For this purpose, several non-CRC disease datasets were assessed, including those from gastrointestinal disease (inflammatory bowel disease (IBD)) and non-gastrointestinal diseases (type 2 diabetes (T2D) and Parkinson's disease (PD)) (Extended Data Fig.
Data are presented as the mean ± s. d. of three biological replicates. P values were calculated using a two-sided Wilcoxon signed-rank test and were Bonferroni-adjusted. The box plots show the IQRs as boxes, with the median as a black horizontal line and the whiskers extending up to the most extreme points within the 1. 5× the IQR. e, AUROC matrix of models built with the 175 important EggNOG genes. Values on the diagonal refer to the average AUROC of 20× repeated fivefold stratified cross-validations. Values off the diagonal refer to the AUROCs obtained by training the model on the population of the corresponding row and applying it to the population of the corresponding column. The LOCO row refers to the performances obtained by training the model using all but the cohort dataset of the corresponding column and applying it to the dataset of the corresponding column.
Discussion Most studies have primarily focused on the bacterial microbiota and its effects on human health and disease  12, 13, 30, 31. Recently, investigations have revealed the critical roles of non-bacterial microorganisms in human diseases  32, 33  (Supplementary Discussion). In this study, we performed a comprehensive analysis on the multi-kingdom microbiome using CRC metagenomic datasets across eight different cohorts. We discovered a series of both bacterial and non-bacterial markers and evaluated their performance in detecting patients with CRC across cohorts. We showed that fungal, archaeal and viral species could separate patients with CRC and healthy controls across multiple geographical cohorts (Fig. 2  and Extended Data Figs. 2 3 4 ). However, the predictive value of different kingdom models varied and the bacteria-and fungi-based models, respectively, showed superior accuracy over the archaea-and virus-based models generally (Fig.
2b ). Notably, these models showed some preferences at the population scale (Fig. 2b, c ), which may be due to differences in geography and lifestyle (Supplementary Discussion). Nevertheless, our findings emphasize the need for integrated analysis to identify universal cross-cohort microbial features for accurate CRC diagnosis. Previous studies proposed paradigms to identify reproducible microbial biomarkers across multiple datasets and populations by developing machine learning models, followed by cross-study and leave-one-out likely validations  12, 13. Similarly, we developed diagnostic models with multi-kingdom species that significantly improved predictive accuracy (Supplementary Discussion, Fig. 3  and Extended Data Fig. 5 ). Moreover, models based on the 16-feature panel achieved very high predictive values for CRC diagnosis (average AUROC = 0.
TIS individuals without physician-diagnosed dementia, stroke, cancer, cardiovascular disease, psychiatric disorders or other serious illnesses were recruited. Written informed consent was obtained from all individuals before data and biospecimen collection. Use of TIS individuals was approved by the Ethics Committee of the School of Life Sciences, Fudan University (institutional review board approval no. 496). Stool samples were collected in faecal collection tubes and were immediately transferred to a -80 °C freezer until time for use. The gDNA of the faecal specimens was extracted with a Stool Genomic DNA kit (catalogue no. CW2092S; CWBIO) according to the manufacturer's instructions except for the modification of step 4 with bead-beating for 10 min (Glass beads, acid-washed; catalogue no.
G8772; Sigma-Aldrich) to better extract fungal DNA. The details of the DNA extraction method are given in the Supplementary Methods. Sequencing libraries were generated with the NEBNext Ultra DNA Library Prep Kit for Illumina (New England Biolabs) and library quality was confirmed with an Agilent 2100 Bioanalyzer and quantified using real-time PCR. Whole-genome sequencing was carried out on the NovaSeq 6000 system (Illumina). All samples were paired-end sequenced with a 150-base pair (bp) read length to a targeted data set size of 12 Gb. No statistical methods were used to predetermine sample size but our sample sizes are similar to those reported in previous publications  [11] [12] [13]. Public populations of patients with CRC and controls. Raw sequencing data of eight populations from seven countries were downloaded from the Sequence Read Archive (SRA) (details shown in Supplementary Data 1), mainly from two recently published CRC papers  12, 13  and the Japanese cohort  9.
Metadata were manually curated from the published papers. Study design. We included a total of 1,368 samples from 9 geographical populations of faecal shotgun metagenomic sequencing data, including publicly available and in-house (CHN_SH) sequencing data. To obtain universal microbial features across different countries, we divided these samples into discovery and validation datasets, with broad regional origin (Supplementary Data 1). Sequencing data of the populations from Austria (AUS, PRJEB7774), France (FRA, PRJEB6070), Germany (GER, PRJEB27928), China (CHN, PRJEB10878) and Japan (JPN, PRJDB4176) were included in our discovery dataset. In total, there were 494 controls and 491 patients with CRC, which included 318 patients with early-stage (stages I and II) and 173 patients with advanced-stage CRC (stages III and IV).
The validation dataset consisted of populations from the USA (USA, PRJEB12449), Italy (ITA, SRP136711) and China (CHN_SH, in-house). Sequencing data preprocessing. The KneadData (http://huttenhower. sph. harvard. edu/kneaddata) v. 0. 6 tool was used to ensure data consisting of high-quality microbial reads free from contaminants. Low-quality reads were removed using Trimmomatic (v0. 39) (SLIDINGWINDOW:4:20 MINLEN:50 LEADING:3 TRAILING:3). The remaining reads were mapped to the mammalian genome (hg38, felCat8, canFam3, mm10, rn6, susScr3, galGal4 and bosTau8; UCSC Genome Browser) and 21,288 bacterial plasmids (National Center for Biotechnology Information (NCBI) RefSeq database accessed in January 2020), 3,890 complete plastomes (NCBI RefSeq database accessed in January 2020) and 6,093 UniVec sequences (NCBI RefSeq database accessed in January 2020) by bowtie2 v.
2. 3. 5 (ref. 46 ); matching reads that were potentially host-associated and laboratory-associated sequences were removed as contaminant reads. Microbial taxonomic and functional profiles. Taxonomic profiling. Taxonomic classification of bacteria, archaea, fungi and viruses was assigned to metagenomic reads using Kraken2, an improved metagenomic taxonomy classifier that utilizes k-mer-based algorithms  47. A custom database consisting of 18,756 bacterial, 359 archaeal and 9,346 viral reference genomes from the NCBI RefSeq database (accessed in January 2020) and 1,094 fungal reference genomes from the NCBI RefSeq database (accessed in January 2020), FungiDB (46) (http://fungidb. org) and Ensemble (accessed in January 2020) (http://fungi. ensembl. org) (accessed in January 2020) was built using Jellyfish (v2.
2. 6. 3 was used to predict genes via the metagenome mode (-p meta). A non-redundant microbial gene reference was constructed with CD-HIT using a sequence identity cut-off of 0. 95 and a minimum coverage cut-off of 0. 9 for the shorter sequences. The reference was annotated with EggNOG mapper v. 2. 0. 1 based on EggNOG orthology data. Moreover, gene abundance was estimated with CoverM v. 0. 4. 0 (https://github. com/ wwood/CoverM) by mapping high-quality reads to reference sequences. An index was created against contigs from the non-redundant genes that originated via the Burrows-Wheeler Aligner (BWA). Clean reads were then mapped to the contig index (BWA MEM) and SAM files were converted into BAM files via SAMtools. Then, CoverM was used to calculate the coverage of genes in the original contigs (coverm contig).
Nature Microbiology Generalization of microbial markers. To further test the generalization of CRC microbial markers across technical and geographical differences in multiple populations, we extensively validated the diagnostic models with cohort-to-cohort transfer validation and LOCO validation as described previously  12, 13. Briefly, in cohort-to-cohort transfer validation, the models were trained on a single cohort and their performances were assessed. In LOCO validation, four out of five cohorts in the discovery dataset were pooled as a training set and the remaining cohort was used as an external validation set. Independent validation with external datasets. Furthermore, we used three additional datasets from Italy (ITA), the USA and China (CHN_SH) to perform independent validation analysis and test the robustness of features as CRC diagnostic markers.
Like model construction in the discovery cohorts, fivefold cross-validation models were constructed with the identified best panel of multi-kingdom microbial markers and evaluated with the average AUROC. Additionally, we performed cohort-to-cohort and LOCO analyses to further test the robustness of the identified markers. In cohort-to-cohort analysis, models were trained with each cohort in the discovery dataset and tested with each cohort in the validation dataset; in LOCO analysis, models were trained with the combined five cohorts from the discovery dataset and tested with each cohort in the validation dataset.
Associations with an FDR < 0. 00001 were included in the downstream analysis. Network was visualized with Gephi v0. 9. 2. Associations between species and function. Spearman associations between microbial species and their functions were performed using the Hierarchical All-against-All method v. 0. 8. 17 (http://huttenhower. sph. harvard. edu/ halla), a computational method used to find multi-resolution associations in high-dimensional, heterogeneous datasets. Associations with an FDR < 0. 01 were included in the downstream analysis. qPCR of potential CRC-associated genes. To quantify the abundance of the oraS, oraE, bdhA and bdhB genes, qPCR as outlined by Wirbel et al. 13  was performed on a subset of gDNA prepared from randomly selected samples of the CHN_SH cohort (24 controls and 24 patients with CRC).
The primers used for validation are listed in Supplementary Data 12 and the patient characteristics are summarized in Supplementary Data 13. Total microbial DNA was extracted using a Stool Genomic DNA kit according to the manufacturer's instructions; DNA concentration was determined using NanoDrop. The PCR reactions were prepared with the TB Green Premix Ex Taq II (Tli RNaseH Plus) (catalogue no. RR820A; Takara Bio) containing 0. 6 μM of primer and 5 ng of gDNA in a 25 μl final reaction volume. Reactions were performed on a CFX96 Real-Time PCR Detection System (Bio-Rad Laboratories). The cycling programme was set as indicated: initial denaturation at 95 °C for 30 s; 40 cycles of 95 °C for 5 s; 55 °C for 30 s; and 72 °C for 30 s, followed by melting curve analysis. Gene expression levels were evaluated using the Ct method described previously  13.
Differential species of archaea and virus among populations. a, UpSet plot showing the number of differential archaea species identified via MaAsLin2 in each population and shared by combinations of datasets. The number on top of each column represents the size of differential species. The set size on the right represents the number of differential species in each cohort and connected dots represent the common differential species across connected cohorts. b, UpSet plot showing the number of differential viral species identified via MaAsLin2 in each population and shared by combinations of datasets. The number on top of each column represents the size of differential species. The set size on the right represents the number of differential species in each cohort and connected dots represent the common differential species across connected cohorts.
Acknowledgements We thank all the enrolled participants who contributed specimens and data to this study. We thank  G. Zhao  for helpful comments and discussions. We acknowledge funding from the  MOST Key R&D Program of China  (no. 2020YFA0907200  to N. -N. L. ; no. 2018YFC2000700  to H. W. ; no. 2017YFC0907500  to X. C. ),  National Natural Science Foundation  (nos. 82030099  and  81630086  to H. W. ; no. 31900129  to N. -N. L. ; no. 82073637  to X. C. ; no. 81774152  to R. Z. ; no. 81770571  to L. Z. ; no. 82000536  to N. J. ),  Shanghai Public Health System Construction Three-Year Action Plan  (no. GWV-10. 1-XK15 ), the  Program for Young Eastern Scholar  at the  Shanghai Institutions of Higher Learning program  (no. QD2018016  to N. -N. L. ), Innovative research team of high-level local universities in Shanghai, Medicine and  Engineering Interdisciplinary Research Fund of Shanghai Jiao Tong University  (no.
Data availability The metagenomic sequencing data of the China-SH validation cohort are deposited in both the NCBI SRA under accession no. PRJNA731589 and the National Omics Data Encyclopedia under accession no. OEP001340. The other raw metagenomic data are available in the SRA (https://www. ncbi. nlm. nih. gov/sra) and European Nucleotide Archive (https://www. ebi. ac. uk/ena/) under accession nos. PRJEB7774, PRJEB10878, PRJEB6070, PRJEB27928, PRJDB4176, PRJEB12449 and PRJNA447983. Source data are provided with this paper.
2 | 2 Differential species of archaea and virus among populations. a, UpSet plot showing the number of differential archaea species identified via MaAsLin2 in each population and shared by combinations of datasets. The number on top of each column represents the size of differential species. The set size on the right represents the number of differential species in each cohort and connected dots represent the common differential species across connected cohorts. b, UpSet plot showing the number of differential viral species identified via MaAsLin2 in each population and shared by combinations of datasets. The number on top of each column represents the size of differential species. The set size on the right represents the number of differential species in each cohort and connected dots represent the common differential species across connected cohorts.
RESULTS Data set and site location. The data set consisted of 21 publicly available environmental metagenomes encompassing marine and continental serpentinite-hosted systems, hot springs, and basalt-hosted hydrothermal vents spread worldwide (see Table  S1  and Fig. S1  at https://doi. org/10. 5281/zenodo. 6597409). They originated from 10 well-studied sites comprising PBHF, LCHF, CROMO, Santa Elena ophiolite, Voltri Massif, Cabeço de Vide aquifer (CVA), Tattapani and Yellowstone hot springs, Axial Seamount (Juan de Fuca Ridge), and the Piccard hydrothermal field (mid-Cayman rise). For each metagenome, the type of habitat, the depth, the temperature, and the pH are reported in Table  1. All selected serpentinizing systems emit highly alkaline fluids, with pHs ranging between 9 and 12. 3, except for sample QV1.
2 in CROMO (pH 7. 9). Fluids from other ecosystems had relatively circumneutral pH values. De novo assembly of each metagenome generated between 10,084 and 401,248 contigs, depending on the sample (see Table  S2  at https://doi. org/10. 5281/zenodo. 6597409). This wide variability of contigs resulted from variations in sequencing depth variability between metagenomes. Consequently, the number of predicted genes ranged from 16,780 to 579,028. Their rate of annotation to known functional categories reached, on average, 55% and 38% based on the Cluster of Orthologous Groups (COG) and KEGG databases, respectively. Microbial community composition. The overall microbial taxonomic diversity identified in the set of metagenomes was mostly affiliated with 22 phyla (19 bacterial and 3 archaeal phyla [Fig.
Additionally, the hierarchical clustering of all genes annotated with KEGG Orthology (KO) did not regroup the serpentinizing sites (Fig. S3  at https://doi. org/10. 5281/zenodo. 6597409) but pointed out functional similarities between LCHF and the basalt-hosted hydrothermal vents. Indeed, their metagenomes were enriched for genes encoding transporters (e. g. , for sugar, phosphate, and Mg 21 ) and ion channels as well as for genes often associated with coping strategies under fluctuating environmental conditions, including the comFC genes, related to bacterial competence  (39) , or the dppA gene, linked to chemotaxis  (40). Likewise, two genes encoding alginate biosynthesis and transport (alg8 and eexD, respectively) were recovered exclusively in metagenomes from deep-sea hydrothermal fields and could contribute to the formation of biofilms  (41).
The same overrepresentation of phn genes was obtained using a COG-based annotation in the random forest clustering (Fig. S4  at https://doi. org/10. 5281/zenodo. 6597409). The carbon-phosphorus lyase operon. The seven identified genes (i. e. , phnGHIJKLM) belong to the operon coding for the C-P lyase and, therefore, are involved in the breakdown of various phosphonates  (43). These genes constitute the minimal catalytic core unit, essential for C-P bond cleavage  (44). Remarkably, the seven genes (phnGHIJKLM) of the C-P lyase were almost exclusively present in serpentinite-hosted ecosystems (Fig. 5  and Table  S3  at https://doi. org/10. 5281/zenodo. 6597409). The even abundance of phn genes within each metagenome most likely reflected their organization in an operon. In the serpentinizing sites studied, we estimated that 5% to 52% of microbial genomes (according to the normalization against the abundances of 10 single-copy housekeeping genes) could contain the C-P lyase essential genes, a rate comparable to the 30% detected in the metagenomes of the Sargasso Sea, where phosphonate degradation is a prominent process (45) (see Table  S4  at https://doi.
org/10. 5281/zenodo. 6597409). Taxonomic annotation of the phn genes identified in the metagenomes distributed them in several microbial phyla (Fig. 5 ). They were predominantly affiliated with Alphaproteobacteria and Firmicutes at LCHF and PBHF, respectively. In the continental serpentinizing ecosystems, most of phn genes were affiliated with Betaproteobacteria members. Phylogenetic inferences of the phnJ gene, encoding the catalytic subunit involved in the cleavage of the C-P bond, were used to ascertain these assignments (Fig. S6  at https://doi. org/10. 5281/zenodo. 6597409). However, a high percentage of phn genes could not be assigned to known taxa ("not classified" in Fig. 5 ) with the lowest-common-ancestor (LCA) algorithm, especially in submarine sites. In our metagenomic data set, the genomic organization of the phn operon differed from the canonical form (phnGHIJKLM) by gene order rearrangements and frequent insertion of additional genes encoding acetyltransferases, transporters, or phosphatases (Fig.
S5  at https://doi. org/10. 5281/zenodo. 6597409). Such rearrangements are relatively common, as well as is the integration of other genes related to phosphorus cycling  (44). In particular, phn gene operons reconstructed from the metagenomes of our study included genes coding for phosphonate transporters (phnCDE). Although the phosphonate transporter-encoding genes seemed particularly enriched in serpentinizing ecosystems than for other metagenomes (except for TAT. 3 [Table  S3  at https://doi. org/10. 5281/zenodo. 6597409]), we observed a notably low detection of genes involved in different phosphonate biosynthesis pathways. In some metagenomic phn operons, the phnF gene (encoding a transcriptional regulator) is replaced by other genes (dasR, yvoA, or yycF) encoding transcription regulators of the GntR family normally implicated in the regulation of genes from the N-acetylglucosamine-degrading pathway  (46).
Finally, an atypical operon identified in the metagenomes of PBHF showed an insertion of a NADphosphite oxidoreductase gene (ptxD, encoding phosphite degradation) between the genes phnE and phnG (Fig. S5  at https://doi. org/10. 5281/zenodo. 6597409). This novel genomic organization could degrade phosphonate or phosphite since the Pn transporter encoded by phnCDE genes can also transport the phosphite ion inside the cells  (29). Compared to the case for the C-P lyase pathway, the genes involved in alternative phosphonate degradation pathways  (47)  were more sparse in serpentinizing systems, in particular genes involved in 2-aminoethylphosphonate degradation (reported in Table  S3  at https://doi. org/10. 5281/zenodo. 6597409). Although we noted that genes involved in phosphate transport (pstSABC) were widely distributed in all metagenomes, the stress marker gene (phoB) for the inorganic phosphorus (Pi) depletion regulatory system  (48)  was enriched in serpentinite-hosted environments (in addition to TAT.