附件-代谢环境的大规模重建和系统发育分析（Large-scale reconstruction and phylogenetic analysis of metabolic environments）

作者：Elhanan Borenstein，Martin Kupiec， Marcus W. Feldman， and Eytan Ruppin
翻译：Wendy

种子集识别对丢失或错误数据的敏感性

  通过模拟检查了丢失或错误的代谢数据对获得的种子集的组成的影响。我们使用酿酒酵母的代谢网络（经过严格研究的人，其代谢代谢的种属）作为参考“完整”网络。使用大肠埃希氏菌的代谢网络，另一种代谢信息广泛的物种产生了几乎相同的结果。我们通过删除现有反应或以不同概率更改反应方向来扰乱该网络。 由于基于低至50％基因覆盖率的代谢图重建仍可检测到反应总数的70％，因此我们分析了高达30％的缺失概率。然后将为扰动网络获得的种子集与原始的不受扰动网络的种子集进行比较，测量真实阳性的百分比（仍被识别的原始种子集合的百分比）和假阳性（错误检测的种子的百分比）。从图S1A可以看出，正确鉴定出的种子化合物的百分比仅略低于不完全网络中包括的反应的百分比，这可以认为是鉴定出的种子的百分比的上限。

此外，对于所有丢失的数据级别，所获得的种子集的大小与原始网络的种子集的大小相比（可以计算为真阳性的百分比加上假阳性的百分比）几乎是恒定的。当错误地改变了许多反应的方向性（而不是完全删除反应）时，真实阳性的百分比更高（即，正确识别了原始种子集的较大部分），但错误检测到的种子化合物的百分比也增加了，导致种子集大小略有膨胀（图 S1B）。但是，总而言之，似乎种子的识别过程并没有放大代谢数据中的噪音或不完整程度。

  值得注意的是，观察到的鲁棒性在某种程度上可能与代谢网络的鲁棒性结构有关。具体地说，网络中存在替代途径-基因敲除时面对代谢稳健性的主要影响因素-可以确保在许多情况下，种子化合物参与一个以上的途径，从而使其对这些途径之一中的数据丢失的检测更加可靠。

紧密连接的组件统计信息

  根据先前的研究，对于大多数物种，我们的SCC分解会导致蝴蝶结结构，其中大部分化合物（在我们的分析中，平均为36.3％）构成巨大的强连接组分，而 其他化合物则排列成相对较小的组分（平均每组分1.4个化合物），且往返于巨大组分的路径较短（图S5）。

  包含不止一种化合物的源成分（即，没有进入边缘的成分）特别受关注。这些成分在种子预测中会产生歧义，因为无法仅从网络拓扑中确定要消除的化合物中包括哪种化合物。例如，此类歧义可能会妨碍对自动催化循环中必要输入和辅助因子的准确识别。在本文介绍的分析中，我们通过将这些组分中的所有候选化合物视为种子来解决此问题。但应注意的是，这样的源成分相对较少，每个物种中的大多数源成分（平均89％）都是单子，仅包含单一化合物，因此在种子集的组成上没有歧义（图S7A） 。

少于9％的化合物包括两种化合物（即两种可互换的化合物），而2.4％的化合物包括三种或更多种。 根据定义，由于较大的来源成分可能会为种子集贡献更多的化合物，因此我们还计算了较大来源成分中种子化合物的百分比。再次，我们发现每种物种中较高百分比的种子化合物形成单一成分（77％），只有三种化合物或更多化合物的源成分的一部分为8.5％（图S7B）。

环境属性

  分析中包括的环境特性是NCBI基因组计划提供的那些特性，它们是原核生物偏爱环境的描述（材料和方法）。利用这些特性的其他动机来自以下考虑因素：氧不仅显示出在塑造电子转移“市场”和跨膜蛋白合成中起着重要作用，而且还允许氧的主要转变、生化代谢网络和复杂性。同样，生化环境的可变性与我们分析中使用的栖息地特性（例如与宿主相关的，专门的，陆生的栖息地）密切相关，已显示出极大地影响细菌代谢网络的结构。因此，这两个环境属性（即氧气需求量和栖息地）预计将直接影响网络的拓扑结构，从而可能影响种子集的组成，使其成为分析的明显候选对象。另一方面，温度和盐度显然是主要的环境特征，引起了人们的极大兴趣，但它们对代谢网络拓扑的潜在影响可能不太直接。因此，这四个属性合在一起构成了一个适当且均衡的集合，用于研究种子集的组成与全球环境特征之间的关系，并研究了环境对代谢网络的直接和间接影响。

2-羟基戊二酸的系统形态

  已经发现柠檬酸（TCA）循环不完整，需要在几种专性细胞内物种（包括衣原体和蚜虫（Buchnera aphidicola））中外源获取2-氧戊二酸。其他专性寄生虫，例如Mollicutes，表现出极少的新陈代谢，并且完全没有TCA循环（16）。通过我们对2-氧戊二酸的分析获得的种系发生和种子模式与上述研究完全一致； 2-氧戊二酸酯在所有衣原体和布氏杆菌中均以种子的形式出现，在所有Mollicutes中完全不存在，并且在所有其他物种中均是存在的化合物（图S8）。数据表S1和数据集S2中有一张描述了所有化合物的系统谱的表，其中包括我们的分析，每种模式的几种计算方法的软相干性（SI材料和方法）。

成为种子的概率与拓扑特征之间的相关性

  某种化合物为种子的物种（网络）数量与该化合物出现的物种数量之比（Ns / No）提供了该化合物为种子的概率的度量，可用于检查 化合物的哪些特性使其更可能包含在种子集中。具体而言，将这种方法与该化合物在全球代谢网络中的拓扑特征进行比较[对应于生物圈的代谢组的集体潜力]，可以发现一些重要的相关性。此概率与化合物的作用范围（其具有路径的其他化合物的数量）及其中心（这些路径的平均长度）（0.59，P <10-300和0.52，P<10-300）分别相关；Spearman等级相关性）。此外，通过仅考虑在网络SCC分解过程中从未修剪过的化合物来滤除网络数据中的一些噪声，成为种子的概率与化合物的流出程度相关，而与其进入程度成反相关 （Spearman等级相关系数分别为0.36，P<10-5和0.47，P<10-9）。这些发现表明，种子化合物往往是位于网络外围（即远离核心代谢）的种子化合物，而是许多其他化合物的前体。

  可能是种子的化合物（Ns / No>0.5）倾向于与某些代谢途径（SI材料和方法）相关，包括脂肪酸生物合成和氨酰基-tRNA生物合成（多次测试校正后的P <0.05）。脂肪酸生物合成途径的富集与对代谢网络进化的最新研究一致，该研究发现许多新反应的获得主要发生在脂质代谢途径中。

种子集数据的主成分分析

  种子集数据的主成分分析（PCA）用于检查种子集空间中不同分类单元的分布。前两个主要成分对不同分类组的明确划分表明，结实成分是每个物种的良好特征（图S9A）。通过校正分析中包含的大量细菌类群，可以改善这种划分（图S9B）。

代谢网络的重构

  从KEGG LIGAND数据库的文件reaction_mapformula.lst中检索了数据库中主要反应的列表。 该文件还列出了每个反应的定义（即底物和产物化合物）以及其参与的每个途径的方向性（如果已知）。化合物仅限于主要反应物。对于每个物种，从PATHWAY数据库中的rn文件中检索每个路径中存在的反应列表。使用该反应-路径对列表，以及上面获得的反应的定义和方向，重构了每个物种的代谢网络。 该网络以有向图表示，其中节点表示化合物，边缘表示反应。从化合物a到化合物b的有向边表示化合物a是某种产生化合物b的反应的底物（即，对于每个给定的反应，代表其底物的所有结点都通过有向的边与代表其产物的所有结点相连）。图中省略了聚糖。 可逆反应或方向性未知的反应表示为两个方向上的有向边。我们还记录了每个物种网络中出现的反应和化合物的数量。对于每种化合物，还从LIGAND数据库中的化合物文件中检索了其参与的代谢途径。

  由于固有噪声和不完整的反应数据，重建的网络包含大量小的断开连接的组件（即未连接到网络主要部分的节点组），可能会明显干扰有意义的种子化合物的检测。在执行其余分析之前，将从网络中删除任何包含10种或更少化合物的此类成分。 我们将这些滴下成分中包含的化合物称为修剪化合物，在分析中将其视为种子状态未知的化合物。

紧密连接的组件分解

  给定网络G，由Kosaraju的算法执行强连接组件（SCC）分解，其工作原理如下：

1. 在G上运行深度优先搜索（DFS），以计算每个节点v的完成时间f [v]。
2. 计算转置的网络G’（每个边的方向都相反的网络G）。
3. 在G’上运行DFS，以f [v]的降序遍历节点。
   第二个DFS创建的DFS森林中的每棵树都形成一个单独的SCC。

系统发生模式和系统种子模式

  系统发育模式代表了所分析物种中存在和不存在模式特异的性状。例如，考虑到直系同源基因的集合，可以将某个基因的系统发育模式设想为布尔向量，表明可以在其中找到直系同源物的物种的集合。在本文的上下文中，每种化合物的系统发育模式都相关（图S2）：每个化合物的系统发生模式都是二元载体，表明该化合物在哪种物种中发生。同样，每种化合物的种系种子模式均为二元载体，表明该化合物为种子的物种（图S2）。考虑特定允许的化合物的物种分布格局，以研究其在现存物种中的状态（例如种子与非种子），并使用最大简约性或最大似然方法在树的内部节点（代表祖先物种）中追踪化合物的状态。

检测相关模态

  我们希望检测与主要环境变化或活生物体代谢变化相对应的化合物或种子模式。这些化合物的系统发育模式都应显示出与系统发育树拓扑的高度一致性，并诱导该物种有意义地划分为存在化合物/种子的物种和不存在该化合物/种子的物种。但是，考虑到数据的噪声性质（源于分类和注释原基因涉及的固有噪声），尤其是这种噪声对产生的事件和种子系统发育模式的潜在影响，简单的简约分析可能会产生误导。 为了检测这些模式，我们因此使用了一种基于 信息增益 的新颖方法。

  正式地，假设在给定的涵盖L种的系统发育模式中，该种的P具有1（“存在”）状态，而Q具有0（“不存在”）状态。使p=P/L，q=Q/L表示相对频率。然后通过H=-plog§-qlog(q)的模式给出熵的定义。对于树中的每个内部节点（假设树已植根），可以将这些物种划分为L1物种和L2物种，L1物种是内部节点的后代，而L2则不是。分别用P1，Q1和P2，Q2表示这两个组中每个组的存在/不存在种子的数量。此外，用p1=P1/L1，q1=Q1/L1，p2=P2/L2，q2=Q2/L2表示相关频率。每个组中的熵由下式给出H1=-p1log(p1)-q1log(q1)和H2=-p2log(p2)-q2log(q2)。因此，该节点的信息增益为IG = H - [H1(L1/L)+ H2(L2/L)]。我们遍历了树中的所有内部节点，发现其中一个具有最大的信息增益。将此信息增益值（及其划分的结果）分配给化合物及其系统发育模式。根据其信息增益对所有化合物进行分类，并检查那些信息增益值最高的化合物。 可以在数据集S2中找到描述每种化合物所有获得的量度的表格。

路径选取

  从KEGG中检索了每种化合物参与的代谢途径。给定一组化合物（例如，其Ns / No为0.5的那些），计算了该组化合物中参与每种途径的化合物的数量。我们对10,000个相同大小的随机化合物集重复了相同的过程，以计算在给定的集合中明显过量或不足的途径。通过错误发现率程序对所得的P值进行进一步校正以进行多次测试。

预测病原菌中外源氨基酸和辅因子

  从参考文献23中检索了有关埃希氏菌病病原体氨基酸和辅因子生物合成的数据，参见其中的表5。这些数据涵盖了三个新测序的病原体（吞噬嗜酸菌，沙眼埃里希氏菌和新里克立克体），以及来自立克次体的其他物种（Anaplasmaginginale，Ehrlichiaruminantium，Wolbachia pipientis wMel和立克次氏菌） Blochmannia floridanus和Wigglesworthialossinidia）。我们丢弃了布氏杆菌属种（以避免重复使用布氏菌相关数据进行验证），总共留下了9种。 在每个物种中，都有合成20种氨基酸和10种维生素/辅因子的能力。我们检查了通过分析获得的种子集，并检索了相应的数据集，描述了同一9个物种中这30种化合物中每种化合物的种子状态（表S3）。为了获得有关这些新测序物种的最大信息，我们使用了最新的KEGG汇编（发行45.0，2008年1月1日），并且没有修剪网络。 由于我们的重点是种子检测算法能够正确地区分来自合成（非种子）化合物的外部获得的种子，因此我们的分析仅限于网络中发现的化合物。我们比较了两个数据集，并检查了参考文献中是否报道了化合物。 我们的算法正确地将23种不会在特定物种中合成的种子识别为种子。为了评估种子预测的准确性，我们将其视为 二进制分类 问题，其中种子检测算法旨在对合成（非种子）与非合成（种子）化合物进行分类。因此，将分类精度定义为（TP+TN）/（TP+FP+FN+TN），其中TP表示真阳性的数量，TN表示真阴性的数量，FP表示假阳性的数量，FN表示 假阴性的数量。通过将物种和化合物的标签改组10,000,000次，并通过概率计算获得等于或更高的准确性的概率，可以计算出由此而来的准确性度量的统计显着性。 同样，着眼于正确的种子预测，精度以TP /（TP+FP）计算召回率以TP /（TP+FN）计算。

Tamura和Nei的替代率估算方法

  Wefollow Tamura和Nei（24）估计了整个系统发生树的替代数量和替代率（另见参考文献25）。他们的方法最初是为DNA序列中的核苷酸取代估计而开发的，但可以类似的方式应用于任何离散性状状态之间的转换。在研究的背景下，乙型肝炎的易感物种与状态序列相关，其中基因座l描述了该物种中化合物1的状态，作为物种DNA序列的类似物。 使用最大的简约原则从现存物种的状态中推断出特异种的始祖物种的状态（树的内部节点） 。通过将每个物种的性状状态与其直系祖先进行比较，计算出每种类型（从状态i到状态j）的替代数目。如果最简约的分配是模棱两可的（即有两个相等的简约分配），则两个状态（以及相关的替换）中的每一个都被认为概率为0.5。在网络重构过程中被修剪的化合物的状态（请参见“代谢网络重构”）也被认为是模棱两可的，可以采用种子状态也可以采用非种子状态。从我们的分析中删除了某些节点中最简约的赋值包括所有三个状态的化合物（参考文献24）。为了估计每种替换类型的相对频率，将替换次数除以分析中所有序列中原始（祖先）状态的频率。 然后表示相对频率，以使所有频率的总和为100％。

系统发育树的重建和评估

  我们再次将分析限制在可以与参考文献26中系统树出现的物种匹配的物种。 我们计算了每对物种的种子集之间的杰卡德距离矩阵，并将邻居连接算法和Fitch-Margoliash算法（分别由PHYLIP程序包的NEIGHBOR和FITCH实现）应用于该矩阵， 重建与物种相关的系统树。 同样，我们根据发生的化合物组与随机种子组（即具有相同大小的种子种子组的随机子集软混合化合物）之间的距离重建了系统树。这些树和参考树之间的距离都通过“分支得分距离”度量和“对称差异”度量进行了评估（由PHYLIP软件包中的TREEDIST程序实现）。 MEGA Tree Explorer（32）用于绘制系统发育树。

环境属性

  我们在这里提供环境属性数据中使用的每个类别的详细说明。 该描述改编自NCBI基因组项目帮助，可在网上找到（www.ncbi.nlm.nih.gov/genomes/static/gprj_help.html）。

  “盐度”描述了细菌的盐度要求（生长培养基中盐的氯化钠当量百分比）：非嗜盐性，0–2％NaCl； 一般嗜盐2–5％NaCl； 中度嗜盐菌5-20％NaCl ； 极端嗜盐菌20–30％NaCl 。

  “氧气”描述了生物体在各种氧气水平下生存的能力：零，未知的氧气需求； 有氧的，有机体可以在氧气的存在下生长，并可能使用氧气作为电子受体。 微需氧的，有机体可以忍受低水平的氧气，并且可能不使用氧气作为电子受体；兼性的，该有机体可以有氧或有氧地生长； 厌氧的生物体在没有氧气的情况下生长，并利用其他电子受体。

  “温度范围”描述了生物生长的温度范围的基本类别（摄氏度）。生长在多个类别上的生物根据其大部分温度范围重叠的类别进行分类：未知，未知该生物在什么温度下生长； 极度嗜冷的，生物体的生长速度为-30至-2；嗜冷的，生物体的生长速度为-1至+10；嗜温的，生物体的生长速度为11至45。 嗜热的，生物体以46至75的速度生长； 高嗜热性，有机体生长在75以上。

  “栖息地”描述了发现生物的基本环境：未知，未知。 与宿主相关，该生物通常或专心与宿主生物相关； 水生生物，这种生物经常或专心与淡水或海水环境相关； 在陆地上，这种生物通常或专心与陆地环境（例如土壤）联系在一起； 专门的，这种生物生活在专门的环境中，例如海洋热气口； 在上述多种环境中，可以发现多种生物。