A catalog of genes and species of the human oral microbiota

Data sourcesThe oral gene catalog was built using three primary sources:Bacterial Genomes from the Human Oral Microbiome Database (HOMD).Fungal Genomes from the NCBI RefSeq database.Metagenomic Sequencing Data from multiple oral microbiome studies.The creation of the oral gene catalog was a multi-step process, combining and refining genes from each source. Bacterial GenesA total of 1,505 bacterial genomes were downloaded from HOMD (version 20170215, accessed in December 2017). Genes shorter than 60 nucleotides or containing ambiguous bases were filtered out. Redundancy was removed using CD-HIT-EST (v4.6; parameters: -aS 0.9 -c 0.95 -T 0 -M 0 -t 0 -d 0 -G 0). This process yielded 1,459,394 unique HOMD genes for the catalog. Fungal Genes1,017 fungal genomes were downloaded from NCBI RefSeq (May 2017). For the 492 genomes lacking existing annotations, gene calling was performed using Genemark-ES in fungi mode. After initial redundancy removal with CD-HIT-EST (v4.6; parameters: -aS 0.9 -c 0.95 -T 0 -M 0 -t 0 -d 0 -G 0), genes were selected for inclusion only if their corresponding genome was present in at least 20% of the samples in one of the metagenomic cohorts, determined by mapping reads with Bowtie2 (v2.2.3). This led to the selection of 2,440,644 fungal genes. Metagenomic Sequencing DataThe gene catalog was supplemented with data from 689 oral metagenomes, including newly sequenced samples, from the following studies:Human Microbiome Project (HMP): 382 samples (bioproject PRJNA255439).Chinese Cohort: 212 samples (bioproject PRJEB6997).TwinsUK Cohort: 48 newly sequenced samples (bioproject: PRJEB38483).Raw reads were subjected to quality control and trimmed using AlienTrimmer 0.4.0 (parameters: -k 10 -l 45 -m 5 -p 40 -q 20). Human sequences were removed by mapping against the human reference genome (GRCh38.p11) using Bowtie2 2.2.3. Metagenomic assembly was performed using SPAdes 3.9.0 (parameters: "-k 21,33,55 --only-assembler –meta" for Illumina paired-end data, or "--iontorrent -t 24 -m 300 -k 21,33,55 --only-assembler" for Ion Torrent single-end data). Contigs shorter than 500 bp or with coverage less than 2x were discarded. Gene calling was conducted with Prodigal (parameters: -m -p meta). Genes shorter than 60 bp were filtered out, and redundancy was removed with CD-HIT-EST (v4.6; parameters: -aS 0.9 -c 0.95 -T 0 -M 0 -t 0 -d 0 -G 0). Final Gene CatalogThe final gene catalog was assembled by sequentially adding non-redundant genes from each data source. Genes from HOMD and fungal genomes were combined first using cd-hit-est-2d. Then, non-redundant genes from the HMP, Chinese, and TwinsUK cohorts were sequentially added using cd-hit-est-2d (same parameters as cd-hit-est). A final redundancy removal step was performed. This process resulted in a catalogue of 8.4 million non-redundant genes MSPs RecoveryThe 689 metagenomic samples were aligned against the final gene catalog using the Meteor software suite to produce a gene abundance table. Then, co-abundant genes were binned into 853 Metagenomic Species Pan-genomes (MSPs) using MSPminer. MSPs Taxonomic AnnotationTaxonomic annotation for the MSPs was performed by aligning all core and accessory genes against representative genomes from the GTDB database (release r214) using blastn (task: megablast, word_size: 16).A species-level assignment was given if over 50% of the genes matched a representative genome with a mean nucleotide identity of at least 95% and a mean gene length coverage of at least 90%.The remaining MSPs were assigned to a higher taxonomic level (genus to superkingdom) if more than 50% of their genes shared the same annotation. Mapping rate distribution across public cohortsWe generated mapping rate distribution plots using Meteor2 (default parameters), comparing performance between: PRJEB6997, PRJNA255439, PRJNA48479 (cohort used in catalogue assembly) and PRJEB24090, PRJEB28422, PRJEB45799 (independent cohort not used in assembly).数据来源 本口腔基因目录构建采用三类主要数据源：人类口腔微生物组数据库（Human Oral Microbiome Database, HOMD）的细菌基因组、美国国家生物技术信息中心参考序列数据库（NCBI RefSeq）的真菌基因组，以及多项口腔微生物组研究的宏基因组测序数据。口腔基因目录的构建为多步骤流程，需整合并优化各数据源的基因序列。 细菌基因 于2017年12月获取HOMD 20170215版本数据，从中下载共计1505个细菌基因组。过滤掉长度小于60 nt或包含模糊碱基的基因，使用CD-HIT-EST（v4.6；参数：-aS 0.9 -c 0.95 -T 0 -M 0 -t 0 -d 0 -G 0）去除序列冗余，最终得到1,459,394个非冗余HOMD基因用于该目录。 真菌基因 于2017年5月从NCBI RefSeq下载1017个真菌基因组。针对其中492个无现有基因注释的基因组，采用Genemark-ES的真菌模式进行基因预测。先使用CD-HIT-EST（v4.6；参数同上）完成初始去冗余，随后仅保留对应基因组在至少一个宏基因组队列20%以上样本中存在的基因——该筛选通过Bowtie2（v2.2.3）比对测序读段（reads）实现。最终筛选得到2,440,644个真菌基因。 宏基因组测序数据 本基因目录补充了来自689份口腔宏基因组样本的数据，包含新测序样本，其来源如下：人类微生物组计划（Human Microbiome Project, HMP）：382份样本（生物项目编号PRJNA255439）；中国队列：212份样本（生物项目编号PRJEB6997）；TwinsUK队列：48份新测序样本（生物项目编号PRJEB38483）。 原始测序读段（reads）使用AlienTrimmer 0.4.0进行质控与修剪（参数：-k 10 -l 45 -m 5 -p 40 -q 20）。通过Bowtie2 2.2.3比对人类参考基因组GRCh38.p11，移除宿主人类序列。宏基因组组装采用SPAdes 3.9.0完成：Illumina双端数据使用参数"-k 21,33,55 --only-assembler –meta"，Ion Torrent单端数据使用参数"--iontorrent -t 24 -m 300 -k 21,33,55 --only-assembler"。丢弃长度小于500 bp或覆盖度低于2×的重叠群（contig）。采用Prodigal进行基因预测（参数：-m -p meta），过滤掉长度小于60 bp的基因，并使用CD-HIT-EST（v4.6；参数同上）去除序列冗余。 最终口腔基因目录 最终基因目录通过依次整合各数据源的非冗余基因构建得到：首先使用cd-hit-est-2d合并来自HOMD与真菌基因组的基因，随后依次使用cd-hit-est-2d（参数与cd-hit-est一致）整合来自HMP、中国队列及TwinsUK队列的非冗余基因，最后执行最终去冗余步骤。本流程最终得到包含840万个非冗余基因的口腔基因目录。 宏基因组物种泛基因组（MSPs）构建 使用Meteor软件套件将689份宏基因组样本比对至最终基因目录，生成基因丰度表。随后通过MSPminer将共丰度基因聚类为853个宏基因组物种泛基因组（Metagenomic Species Pan-genomes, MSPs）。 MSPs分类注释 MSPs的分类注释通过将所有核心基因与附属基因比对至GTDB数据库（版本r214）的代表基因组完成，比对工具为blastn（任务：megablast，字长：16）。若超过50%的基因与某一代表基因组匹配，且平均核苷酸一致性≥95%、平均基因长度覆盖度≥90%，则为该MSP赋予物种水平分类标签。对于剩余未完成物种注释的MSPs，若超过50%的基因具有相同分类注释，则将其归类至更高分类阶元（属至超界）。 公共队列比对率分布 使用Meteor2（默认参数）生成比对率分布图表，对比以下队列的比对表现：PRJEB6997、PRJNA255439、PRJNA48479（均为用于构建本基因目录的队列），以及PRJEB24090、PRJEB28422、PRJEB45799（未用于目录构建的独立验证队列）。