Nature | 原核生物基因的生物地理学研究

摘要

微生物基因编码了地球上生命的大部分功能。目前人们对基因在全球生物圈的分布以及对人类和地球健康的影响知之甚少。本研究构建了一个来自14种生境的公开可用的宏基因组的3.03亿个物种级基因非冗余基因集。在多个生境中发现的部分基因富含抗生素抗性基因和可移动遗传元件的标记。通过物种水平的基因进一步聚类为3200万个蛋白质家族，发现家族中包含了大部分基因，大多数物种水平的基因和蛋白质家族是罕见的。此外，物种水平的基因，尤其是稀有基因，显示出较低的适应性选择，也就是说大部分同源蛋白的编译基因都是无义突变（中性或近中性的变异）。

英文标题：Towards the biogeography of prokaryotic genes

中文标题：原核生物基因的生物地理学研究

期刊：Nature, December 2021

第一作者：Luis Pedro Coelho

通讯作者：Luis Pedro Coelho, Jaime Huerta-Cepas & Peer Bork

作者单位：

Institute of Science and Technology for Brain-Inspired Intelligence, Fudan University, Shanghai, China

内容

1.全球微生物基因集

本研究整合宏基因组和完整基因组，调查不同生境的原核基因，以了解基因的全球分布及其编码的分子功能。整理了14种生境无冗余的全球微生物基因集(GMGCv1)。GMGCv1来源于13174个公开的高质量宏基因组。每个宏基因组中组装contigs并预测开放阅读框(ORF)，得到2,007,736,046个ORF。并从proGenomics2数据库的84,029个高质量基因组中纳入了312,020,843个ORF。使用基于图的冗余去除算法，得到的2,319,756,889个序列，以95%的核苷酸一致性聚集302,655,267个聚类。从每个簇中保留了一个单一的序列，代表了所有核苷酸变异的95%核苷酸同源性——这相当于每个物种一个特定基因的副本，此后被称为“unigene”（单基因簇）。

2.大多数基因具有生境的特异性（Habitat-specific）

虽然MAGs通常是根据每个样本或每个生境建立的，但全球微生物基因集使我们能够识别生境之间共享的基因。由于物种水平的单基因簇代表多个序列(核苷酸同源性大于95%)，它们可能代表来自多个生境的基因(“多生境基因”)。这些可能包含在多种栖息地中茁壮成长的物种中，或者是可移动元件的一部分，也就是说，基因可以在基因组之间或跨栖息地边界水平转移。只有18,145,135个物种级单基因簇(占总数的5.8%，P < 10⁻³⁸)为多生境基因(图1b，附图5)。这与物种倾向于适应其环境的研究结果一致，在宿主相关微生物群落中，同种菌株包含宿主特异性基因。

抗生素耐药性基因(ARG)——被认为是频繁携带可移动元件，比其他单基因簇更有可能出现在多种环境中(3208187个ARG中有329857个（P < 10⁻³⁸)。为了量化栖息地间的物种重叠，考虑到许多物种尚未被了解，构建了每个生境的7443个宏基因组物种(MGSs)，只有1099个在不同生境间共享，这与观察到的单个unigene共享模式一致(Fig.1b)。

图1 全球微生物基因集 a,组装来自 14 个不同栖息地的宏基因组（标记大小代表短读长的总数）并提取 ORF,与proGenomes2 的 ORF 相结合，形成物种水平的单基因簇、蛋白质簇和蛋白质家族. b, 除了哺乳动物肠道微生物群之间的共享外，栖息地之间的单基因簇共享很少。每个色带的宽度代表左侧栖息地中共享基因的平均丰度。最宽的丝带将猫肠道与人类肠道连接起来，代表了猫肠道微生物组中 58.0% 的读长映射到人类肠道. c, 单基因簇积累曲线表明，某些栖息地的每个样本的累计基因多样性稳定，而其他栖息地（例如，海洋和土壤）仍处于采样不足. d, 最大的蛋白质家族包含 73,979 个单基因簇。所有单基因簇的一半只包含在 203,431 (0.6%) 个家族中（那些包含≥239 个物种水平的单基因簇），而 80% 的蛋白质家族只包含一个或两个基因，略低于总单基因簇的 8%

3.丰富结构也具有生境特异性

研究结果表明每个样本中每个物种的物种水平单基因簇平均数量不同(图2a, P < 10−^38,)。海洋和土壤环境表现出多种子模式的混合。在海洋样品中，这些子模式对应于不同的海洋深度，特别是比较浅水样品和那些在无法接触阳光的深水中收集的样品时，而土壤环境的差异伴随着酸度和湿度的差异。因此，在一个宏基因组中，每个物种的单基因簇数量是一个明确的生境的识别特征。

为了检验单基因簇丰富度主要是由包含多个同源单基因簇的群落驱动的，还是由多种功能群驱动的，通过计算蛋白质家族丰富度与物种水平的单基因簇丰富度的比率，发现不同生境之间的明显差异(图2b)。为控制测序深度差异，建立分类器 (taxonomic、系统发育、单基因簇和蛋白家族丰富度)预测每个样本的生境。通过交叉验证，该分类器在14种生境上的准确率为86.1%。

图2 同种基因的数量(每个物种的基因库)和每个宏基因组的功能冗余在不同生境之间的差异显著小于在不同生境之间的差异 a, 每个样本中同种基因数量的密度，表明每个样本泛基因组存在于环境样本中，而不是宿主中. b, 在每个栖息地的每个样本中检测到的每个蛋白质家族（功能冗余的代表）的单基因簇数量的密度显示了栖息地之间的明显差异。 

4.大多数基因都很罕见

研究单基因簇在宏基因组中的频率，发现不同生境的参数不同(图3)，结果表明大多数基因的频率很低。通过测量基因顺序和通路邻近的保守性跨越流行阶层，发现稀有的非基因确实比普遍的非基因表现出更少的功能交互作用(图4a)。 

分析序列变异分析数据是否与进化的中性模型兼容。中性意味着大多数观察到的基因差异对适应(几乎)没有影响，因此不是由于对特定生态位的适应(积极选择)，尽管纯化(消极选择)选择可能仍然是活跃的。蛋白质家族之间的选择方式是不同的，对每个蛋白质家族内测试了阳性(适应性)选择(方法)。结果表明，绝大多数单基因簇没有显示积极选择的证据(图4b)。罕见的单基因簇比普遍的单基因簇(高达10%)的适应性要低得多(图4b)。

研究使用了GMGCv1中5126个注释良好的大肠杆菌基因组的单基因簇，得到了阳性选择增加和基因流行率增加的非常相似的相关性(图4b)。结果表明，与更普遍的大肠杆菌相比，罕见大肠杆菌单基因簇中的位点处于更少的可检测的选择压力下(图4c)。

图3 大多数基因是罕见的 基因流行率的直方图在对数尺度上大致是线性的

图4 稀有单基因簇的选择压力较低 a, 操纵子结构更频繁地保留在流行基因中 b, 在可检测的正选择下的单基因簇的比例随着检测次数的增加而增加 c, 大肠杆菌泛基因组是唯一一个大小足以测试每个位点选择的基因组。大肠杆菌泛基因组中的高流行基因显示出更强的正（蓝色）和负（红色）选择。箱线图和点显示在总比对长度（每个类别的 n = 4,167）中每个单基因簇显著选择下的残基分数，灰线显示选择中至少有一个残基的基因比例

随着测序成本的持续降低，以后可能会捕获地球上所有丰富的原核生物物种。根据此研究的数据，这似乎对生物多样性非常高的生境也是可行的（例：土壤）。然而大量的稀有的、生境特异性的甚至是区域特异性的基因，以及从头基因产生、修饰和淘汰的可能的转换过程，全球基因库的相当一部分可能永远不会被捕获。

方法

1.宏基因组与基因组的筛选

本研究的宏基因组数据从均从European Nucleotide Archive (ENA)下载获得，仅2017年1月1日公开的样品被选入。宏基因组的识别遵循下述两个条件：(1) 样品的分类编号要么是408169，要么是在分类树上位于408169的子树上；（2）文库搜索Source Field必须设为“METAGENOMIC（宏基因组学）”。包含至少1M（百万）条序列、平均读长至少75bp且是在Illumina测序的样品被选作进一步分析的对象。样品按照ENA project项进行分组，至少包含100个样品的ENA project才被保留。之后对所有保留的样品进行人工筛查，其中5项研究的数据被剔除，因为它们包含了真核生物样品或是涉及了扩增子测序。

为了拓宽本研究的群落范围，手动加入了猫肠道和土壤宏基因组数据。这些样品符合上述样品要求，但包含少于100个样品的project同样被保留。

上述数据下载与筛选过程利用Python脚本完成，主要依赖requests模块。进一步的分析显示，有369个样品被误选、含有扩增子数据。这一部分数据仅在构建基因集过程中使用，没有参与本研究其余分析。

本研究的基因组数据选择同proGenomes2数据库一样，收集了一批NCBI数据库中部分高质量基因组。

图1a展示了所有样品的地理分布情况，该图是利用R中maptools包（1.1.0版本）绘制。

图1 全球微生物基因集 a,组装来自 14 个不同栖息地的宏基因组（标记大小代表短读的总数）并提取 ORF,与proGenomes2 的 ORF 相结合，形成物种水平的单基因簇、蛋白质簇和蛋白质家族。

2.重叠群组装与开放阅读框预测

收集到的序列使用NGLess进行处理，剪切质量值低于25的位点后，弃去短于60bp的序列。使用Megahit进行组装（使用宏基因组默认参数），开放阅读框（ORF）使用MetaGeneMark进行预测。这些步骤由NGLess脚本完成。

3.非冗余基因集的构建

非冗余基因集的构建分为以下4步：

1）利用滚动哈希（Rolling Hashes）算法严格去除能够完全作为其他基因序列的子序列的基因；

2）利用DIAMOND进行基因间的两两对比；

3）对上一步得到的匹配结果进行过滤，仅保留具有代表性的结果。简单来说，如果存在A′为B的子字符串，且A与A′的编辑距离（Edit Distance）小于等于A全长的5%时，那么认为A可作为B的代表。对于长度相等或相似的序列，这样的定义与物种水平95%核苷酸一致性标准是等效的。这一步的结果是以图（Graph）的形式保存的，其中节点为基因序列，有向边则对应代表性关系；

4）选择一个优势节点集（原图中所有节点均在该集合中或原图中所有节点均可由该集合中的节点代表），该步骤由贪婪算法完成。

4.GMGCV1的质量控制

尽管大量的去冗余后的单基因簇 (unigene)（共189,105,503个）仅在某一个样品中被组装，但其中74.9%的基因能够在多个样品中通过mapping检测到。大部分单基因簇是不完整的，但至少有91.7%能够并入蛋白家族，其中83.2%并入蛋白家族后至少包括了一个完整单基因簇，另8.5%并成了小的蛋白家族。

单基因簇通过DIAMOND向Uniprot数据库比对来检测可能存在的嵌合体。其中，比对达到70%以上的氨基酸一致性且序列重叠少于10个氨基酸的基因被认为是潜在的嵌合体。仅920,579个基因被检测为潜在嵌合体。

为进一步测试不完整的ORFs在基因集中的影响，对于基因末端重叠情况进行检查，因为如果多个不完整的ORFs共同源于同一片段同一端则很有可能出现这种情况（源于同一片段的两端则无法通过这种方式检查）。检查得到8%的单基因簇具有末端重叠。

另外，还测试了不完整ORFs的分布，为此比较了重叠群中相邻ORFs的比例。按理来说，不完整的ORFs分布不会太广泛，因为广泛分布的片段更有可能拼接识别成完整的ORFs。然而，重叠群上相邻ORFs的总体相关性呈现类似的规律（完整/完整=0.46，完整/片段=0.48，片段/片段=0.49）。

为评价人类宿主污染，基因集被切分为均含5万条序列的子文件进行对人类基因组参考数据库（GRCh38.p10）的blastn比对，比对区段包含由基因组、cdna与45S rRNA区段，e-value阈值设为0.00001。

使用AntiFAM进行了假ORFs的检测，37,428个单基因簇被识别为假阳性。

5.宏基因组组装基因组的构建

MAGs是利用Metabat2默认参数分箱得到的，共获得278,629个bins。MAGs质量由CheckM进行评估，并按照MIMAG阈值评价为高、中、低质量。

6.宏基因组物种构建

使用共丰度聚类方法（co-abundance clustering）对各个生物群中的MGS进行鉴定。只有在至少三个样本中都能观察到的完整单基因簇被聚类。以皮尔森相关系数＞0.9为阈值，计算样本层面的冠层概况（canopy profile）。根据大小、四分位GC范围、标记基因和分类信息，对共丰度基因簇进行筛选过滤。由此产生4773个簇，称为MGS，包含500多种基因。MGS中，如果超过80%的基因被注释为某一物种，且序列同一性大于95%，则认定该MGS属于这一物种。分类学信息不一致的MGS和四分位GC＞10%的MGS均被弃去。在界水平将MGS注释为细菌或古菌，删除标记基因数量小于六个的MGS。

7.MGCv1和参考基因组匹配率的评估

为了评估各生境微生物暗物质的数量，我们建立了一个非冗余基因集，该基因集完全基于全球基因集中使用的测序基因组中ORF的子集，包含44,098,640个非冗余的单基因簇。将宏基因组测序读长与该基因集比对，结果表明基因组测序已经获得了某些生境的大部分生物多样性。例如，人类肠道样本中平均80.3%的短读长可以比对到测序获得的基因组中，与先前研究一致。

8.蛋白质家族簇的计算

我们使用MMseq2的标准设置计算蛋白质家族簇，额外要求氨基酸识别阈值为50%、30%和20%，最小序列覆盖度为50%. 参数设置如下：--min-seq-id 0.2 -c 0.5 -cov-mode 2 -cluster-mode 0（将0.2替换为0.3和0.5分别表示30%和50%的同一性）。蛋白质簇也采用类似方法，最小识别阈值和最小序列覆盖度均为90%，参数设置如下：-min-seq-id 0.9 -c 0.9 -cov-mode 1 -cluster-mode 2。

9.分类信息预测

结合三种方法预测分类信息：

1）与单个物种序列聚类到一起的核苷酸同一性小于95%的单基因簇被归为该物种；

2）对于剩下的单基因簇，与Uniprot数据库的最佳比对结果被用于预测超界（细菌/古菌/真核生物/病毒）；

3）对于上一步预测为细菌或古菌的单基因簇，采用双重BLAST最小公共祖先方法获得最终预测结果。该方法将种水平的分配结果转化为属水平的分类信息。这套方法在从种到域的水平上对GMGCv1中78.4%的大基因簇进行了预测，其中细菌占94.6%，古菌占2.7%，真核生物占1.7%，病毒占0.9%。

10.种内和属内核苷酸同一性阈值的评估

使用Prokka对基因进行注释。对32个属的107个物种（物种簇）进行blastn检索。每个物种簇至少包含十个基因组，基因组数超过20的物种簇被重抽至仅含20个基因组。将某一基因组中的全部基因和同一物种簇下的其他基因组进行blastn，也和同一属下其他物种簇中的基因组进行blastn，总共进行了14686个物种簇内的基因组比较和51368个同属物种簇之间的基因组比较。核苷酸同一性结果见拓展数据Fig. 2a.

11.同源基因内部氨基酸同一性的评估

对于40个跨越细菌和古菌的通用标记基因，计算eggNOG中对应簇的平均氨基酸同一性。eggNOG5中预先计算好的比对信息被用于同一性计算，计算通过HMMER3软件包中的AliStat工具进行。

12.可移动遗传元件的注释

使用DDE重组酶的隐马尔可夫模型对数据集中的可移动遗传元件进行注释。如果有多种预测结果，只保留e值小于0.00001且得分最高的。

13.抗生抗性基因注释

基于综合抗生素抗性数据库（CARD）和ResFams数据库获得基因的抗性信息。使用CARD RGI工具将单基因簇分配到CARD模型中。如果单基因簇没有被这种方法注释，且ResFams隐马尔可夫模型的预测得分超过阈值，则使用ResFams对其进行注释。CARD的三个ARG模型中，我们只保留了CARD中同源的模型，弃去了基因缺失导致耐药性的模型和蛋白质变异模型，因为这些ARG无法可靠地识别。

14.基于k-mer的同源搜索

基因以7-mer作为索引，通过使用4位编码16种可能的氨基酸，每个7-mer被转换为0到2²⁸-1范围内的整数。每个序列以其包含的所有k-mer作为索引。对于所有7-mer，所包含的序列存储为一个递增的整数列表。检索待查询序列中所有7-mer的序列索引，并将其组合在一起，以检索数据库中与待查询序列共有7-mer最多的100个序列。然后，通过Smith-Waterman对这100个序列进行重新排序。这种方法由k-mer-find子目录中的代码实现。

15.宏基因组注释和丰度分析

利用eggNOG -mapper2 (version 2.0.1)对进行功能注释，将222,320,961个 (73.4%)物种水平的单基因簇分配到一个eggNOG直系同源组 (OG)。通过在未被选为单基因簇的冗余集合中随机选择的一组ORFs注释来验证这一方法。当它们被分配到OG时，这些基因的95.4%被注释到同一个OG。为了表征eggnog-mapper在部分ORFs上的性能，只考虑了单基因簇是一个完整ORF而冗余ORF是一个片段的情况。在这类情况下，93.7%的注释都指向同一OG。

在“Contig assembly and ORF prediction”中描述的序列裁剪和过滤之后，宏基因组数据使用minimap2映射到目录中。如果有明确的读长映射到，就被认为是样本中检测了单基因簇。使用NGLess和Jug脚本计算基因和功能丰度。简而言之，丰度被估计为短读长映射到一个给定序列的数量，多个映射(短读长映射到多个序列)被唯一的映射丰度分配。对于交叉样本比较，这些结果通过文库大小进行标准化。

此外，使用mOTUs2使用默认参数通过NGLess获得分类学信息。由于在低生物量样品中可能检测到污染，使用了一套阴性对照来获得可疑的mOTU簇，从结果中排除可能的跨生境物种。随后发现Janthinobacterium lividum [ref_mOTU_v2_1333]存在于多种生境中，这与之前在土壤和水生生境中广泛发现这种极端微生物结果相一致。

16.统计分析

在Python中使用NumPy、SciPy和Pandas进行统计分析。

为了检验多生境基因数量的显著性，将每个样本的生境打乱32次，统计打乱条件下的多生境基因数量。Wilks-Shapiro检验证实了结果为正态分布 (P = 0.98)。(理论上多生境基因总数是很多指示变量的和,对于每个单基因簇，每个编码各自代表性的单基因簇是一个多生境基因)。最终共有89,481,710±996,121 (mean±s.d.) 多生境基因。因此，观测值 (18,145,135)比随机预期值(P < 10⁻³⁰⁰)低71.6 s.d.。

箱式图显示了四分位数 (中位数为一条线)，胡须显示了数据的范围，不包括异常值。Tukey规则定义离群值，即Q1 - 1.5 × (Q3 - Q1)以下的数据点，其中Q1为第1个四分位数，Q3为第3个四分位数;或Q3 + 1.5 × (Q3 - Q1)以上的数据点。

17.单拷贝标记基因方法

使用fetchMG提取单拷贝标记基因。通过以下方法对每个样本内单拷贝OTUs进行估计：(1)计数，对于fetchMG识别的40个COGs中的每一个，需要至少有一个成对的读长被明确地分配，以获取特定COGs的物种估计；(2)平均特定COG的估计值，得到单拷贝OTUs的最终估计值。

采用COG 525 (valyl-tRNA synthetase)进行分类学丰富度的估算。之前的工作已经确定，特异性COG的物种分类阈值非常接近95%。这一选择优于COG 12 (一种GTP结合蛋白，也有接近95%的COG特异性阈值)，因为它平均长度更长 (COG 525和COG 12分别为2007和366个残基)。

为了进行验证，我们使用了上面描述的mOTUs2。在适宜使用mOTUs2估算物种多样性的生境中，两种方法估算的物种丰富度相关性较好 (人类肠道:r = 0.71, P < 10⁻³⁰⁰；人体阴道:r = 0.78, P = 1.1 × 10⁻¹⁰；人体皮肤:r = 0.86, P = 9.2 × 10⁻¹⁴⁰；人体口腔:r = 0.75, P = 3.3 × 10⁻²¹⁰；海洋:r = 0.63, P = 8.3 × 10⁻¹⁶;计算Spearman r，质量控制后样品序列≥100万)。对于其他栖息地的样本，相关性并不总是很高 (例如在猪肠道中，r =−0.08,P > 0.05)，因为这不是mOTUs2工具适当的使用环境。基于此，利用COG 525估算了所有样本的分类学丰富度。

18.多样性分析

通过随机抽样将基因计数表重抽到100万条序列。如果可用的序列数少于100万，那么在这组分析中就不会进一步考虑该样本——即使所有宏基因组在输入时包含≥100万个序列，经过质量过滤后，一些宏基因组包含的序列数少于100万个。该操作由profiles-all/gene中提供的diversity.py脚本执行。

蛋白质家族丰富度来作为功能丰富度的代表。仅使用eggnog-mapper推断的直系同源组的结果是相似的 (Spearman R = 0.83，比较了不同样本的蛋白质家族和直系同源组丰富度; 如果只使用来自研究充分的人类肠道样本R = 0.87)，确保这可以作为功能多样性的有效替代，即使某些蛋白质家族可能包含功能已经分化的非同源成员。

对于分类，使用随机森林分类器 scikit-learn建立100棵树 (使用默认参数)。采用10次分层交叉验证的方法评价分类精度。为了控制类别大小的不平衡，将较大的栖息地随机采样到最多200个样本(因此最大的栖息地最多代表11.8%的数据集)。这通过gmgc.analysis/profiles中的classifier-biome-from-div.py执行。

19.将基因频谱拟合到中性无限基因模型中

将基因频率ck定义为被检测k次的基因数量 (例如，c2是正好在两个宏基因组中检测到的基因数量)。在“无限基因模型 (infinite gene model)” 中，新的基因随机生成，现有基因随机丢失 (对适应度没有任何影响)，预测ck和1/k之间几乎是线性关系。

首先将单基因簇数量矩阵重抽到100万 (见“Diversity analyses”部分;这部分数据绘制在图2中)。排除了在质量控制后少于100万个序列的样本后，剩余的样本少于100个。对于与人类相关的栖息地，当来自同一个体的多个样本存在时，只使用一个 (因为来自同一个体的样本，即使在不同时间收集，也不是独立的样本)。

为了量化拟合优度(goodness of fit)，计算1/k与k = 1，…，100时的估计ck值之间的Pearson相关性。总体而言，相关性为0.989806 (P = 9.1 × 10⁻⁸⁵)，所有生境的相关性都很高(附表6)。

非常高的相关性使我们得出这样的结论：中性的“无限基因模型”很适合宏基因组的基因频谱，并且大多数基因不能处于强选择之下。在低端 (k = 1，…，10) 的匹配度特别高，我们称之为稀有基因 (附表6)。

这一结果与无限基因模型不是一个好的原核基因组模型的断言是一致的。正如在正文指出的那样，稀有基因代表了测序基因组的一小部分。

20.GMGC单基因簇和泛基因组 (pan-genome)簇的选择试验

使用ClustalOmega (version 1.2.4)对每个单基因簇下的所有ORFs的翻译产物进行了多重序列比对，共对一个包含198,208个GMGC 单基因簇的代表集进行了比对。将氨基酸比对反转录为密码子比对，并使用默认参数的FastTree2 (version 2.1)重建系统发育树。整个工作流使用ETE3 (version 3.1.1)执行，参数ete3 build -w standard_fasttree -nt-switch-threshold 0.0 -t 0.5 -launch-time 0.5 -noimg -clearall -nochecks。

我们还分析了大肠杆菌泛基因组中的127,618个单基因簇 (specI cluster 95)。利用Muscle v3.8.对单基因簇中的大肠杆菌蛋白序列进行比对，利用pal2nal将其转化为核苷酸比对。

对于GMGCv1 单基因簇和大肠杆菌基因簇，使用HyPhy version 2.5.1 (www.hyphy.org)进行选择试验。每个位点的选择试验用FUBAR模型(version 2.2)计算，该模型计算每个位点的dN/dS比以及每个密码子的阳性和阴性选择的后验概率。挑选后验概率≥0.95的阳性和阴性选择位点。每个基因选择的位点的比率是通过将选择的位点数除以所使用的比对的总长度来计算的。使用aBS-REL方法对蛋白质家族簇每个分支进行选择测试，该方法运行一个自适应分支位点模型，通过量化ω比 (dN/dS)，在系统发育中密码子位点和单个分支之间发生变化。对于GMGC家族内单基因簇的检测，对所有分支进行探索性分析， Holm-Bonferroni多重检验校正p值为0.05。在这项测试中分析限制在5912个蛋白质家族簇 (175395个单基因簇) 中，其中至少有一个完整的基因模型在比对结果中，并且已经被软件RNACode (version 0.3) 预测为 (P≤0.05) 代表基因比对。具有阳性选择证据的单基因簇的部分仅在完整ORFs代表的单基因簇内计算，以避免任何与不完整序列相关的混淆效应。对大肠杆菌簇使用相同的标准，除了只检测每个GMGC蛋白家族中的大肠杆菌分支，并假设所有簇都代表表达的基因。考虑到每个点位的选择测试可能会被序列采样 (即簇的大小) 和比对长度严重干扰，我们将这些测试的比对大小限制在109到361之间 (这些限制代表平均值±1 × s.d)，并重新平衡随机数据集，使每个稀有类别包含完全相同的簇大小分布。在更广泛的目录中，单基因簇的检测数量和它可用的序列数量之间有很强的联系，这是可以预期的。这种联系在分离株的基因中较弱，因为序列的数量既反映了其在宏基因组中的普遍性，也反映了其在分离株群体中的普遍性，但这并不能准确反映其在更广泛环境中的普遍性。在这里我们利用这种偏倚，对泛基因组进行了保存分析。

21.操纵子功能保守性

在单基因簇基因组环境中，KEGG通路的流行率被用作操纵子类似的功能保存的替代物。对于每个单基因簇，提取所有聚类的ORFs (即在95%核苷酸同源性聚类的ORFs) 的基因组上下文信息。

通过四个相邻基因(两个上游基因和两个下游基因)的窗口观察到的唯一KEGG通路与总KEGG通路的比值 (唯一的KEGG /总KEGG) 计算每个单基因簇的KEGG通路多样性。每个单基因簇的KEGG保存量通过1 - KEGG通路多样性来计算。对10组随机组合的10个稀有类别的GMGC 单基因簇进行KEGG保存评分，每个类别包含1万个单基因簇，其中至少3个或最多1000个ORFs。为了避免碎片序列产生的潜在偏差，从测试中排除了不完整的基因。

参考文献：

[1] Coelho, L. P., Alves, R., del Río, Á. R., Myers, P. N., Cantalapiedra, C. P., Giner-Lamia. (2021). Towards the biogeography of prokaryotic genes. Nature. doi:10.1038/s41586-021-04233-4

原文链接：

https://www.nature.com/articles/s41586-021-04233-4