CoExpPhylo：整合共表达与系统发育分析的新一代植物特殊代谢通路基因发现工具

本研究针对植物特殊代谢通路基因鉴定困难的问题，开发了一种整合多物种共表达分析与系统发育聚类的新流程CoExpPhylo。该工具通过建立直系同源共表达簇（OCGs），成功在花青素、原花青素和黄酮醇等代谢通路中识别出已知基因并发现新候选基因，为植物特殊代谢研究提供了高效的计算生物学解决方案。

　　

随着测序技术的飞速发展，植物基因组和转录组数据呈现爆炸式增长，研究重点已从数据生成转向功能解析。然而，特殊代谢通路中的基因鉴定仍然面临巨大挑战。虽然共表达分析是识别同一通路或生物学过程中相关基因的常用方法，但其存在明显局限性——难以区分因共享调控触发而共表达的基因与真正参与同一通路的基因。

为了提升功能预测的可靠性，整合系统发育分析通过考虑进化保守性提供了额外验证维度。研究人员开发了CoExpPhylo这一计算流程，系统性地结合共表达分析和系统发育方法，基于一个或多个诱饵基因候选，跨多个物种识别参与特殊生物合成通路的候选基因。

CoExpPhylo流程包含多个计算步骤：物种特异性共表达分析、局部序列比对识别直系同源基因、候选基因聚类为直系同源共表达簇（OCGs）、功能注释、全局序列比对、系统发育树生成，以及可选的可视化功能。该工作流程高度可定制，用户可调整相关性阈值、过滤参数和注释来源。

在花青素、原花青素、黄酮醇以及叶黄素和玉米黄质生物合成等多个通路上的基准测试表明，CoExpPhylo能够有效识别已知基因，同时发现新的候选基因。

研究人员主要采用了以下关键技术方法：基于240个物种的转录组数据集，使用kallisto进行表达量定量；通过DIAMOND进行局部序列比对；利用NetworkX构建共表达网络；采用MAFFT或MUSCLE进行多序列比对；运用FastTree、RAxML-NG或IQ-TREE构建系统发育树；使用iTOL进行树可视化。

实施过程

CoExpPhylo采用Python3开发，依赖SciPy v1.13.0、NumPy v1.26.4、Networkx v3.3和Plotly v5.22.0等库，并通过bash脚本调用GNU parallel进行并行处理。流程使用多种外部工具：DIAMOND v2.0.14.152进行局部比对，MAFFT v7.490或MUSCLE v3.8.1551进行全局比对，FastTree v2.1.11、RAxML-NG v1.2.2或IQ-TREE v2.0.7进行系统发育树推断。

输入数据收集

分析基于来自不同目240个物种的转录组数据集。编码序列从Phytozome、NCBI GenBank和物种特定网站收集，RNA-seq数据从NCBI SRA获取并使用fastq-dump处理。使用kallisto v0.44处理计数表，合并为每个物种的单一计数表，所有值为TPM（每百万转录本数）。过滤低质量样本后，为每个数据集生成注释文件。

共表达分析

加载输入数据后，对每个物种单独进行共表达分析。使用SciPy的spearmanr函数计算斯皮尔曼等级相关系数r和相应调整后p值。默认相关系数截止值为0.7，p值阈值默认为0.05。还可调整最大检索基因数和最小累积表达量参数。

局部比对共表达序列

通过DIAMOND识别候选序列，使用默认参数进行blastp运行。候选序列需满足：e值不超过10^-5，比特分数大于100，比对长度超过100，相似性大于80%。这些值可通过参数调整。

生成直系同源共表达簇

使用DIAMOND blastp将序列集合作为查询和数据库进行自比对，结果过滤后通过NetworkX构建图。将具有高序列相似性的序列连接在一个图中，最终将每个OCG的序列写入单独FASTA文件。OCG按共表达分析获得的序列比例排序，包含少于10个序列或来自少于三个物种共表达序列的OCG被排除。

OCG功能注释

可选地，使用参考肽段FASTA文件对OCG进行功能注释。每个OCG随机选择50%序列用于注释，至少使用5个序列。通过DIAMOND将每个OCG的序列与参考数据库比对，根据比特分数检索最佳匹配，选择最常见注释作为该OCG的功能注释，并显示可靠性评分。

全局比对

使用MAFFT或MUSCLE对每个OCG进行全局比对，默认使用MAFFT。比对后修剪序列，仅保留足够占有度的位置，默认占有度截止值为10%。

系统发育树生成

使用FastTree、RAxML-NG或IQ-TREE从修剪后的比对文件构建系统发育树。默认使用FastTree，选项为-wag和-nosupport。输出文件为Newick格式，支持上传至iTOL等可视化程序。

批量上传至iTOL

分析完成后，可将树文件自动上传至iTOL，需要Perl脚本和有效API。同时上传注释文件，标记通过共表达分析收集的序列。

输出文件

CoExpPhylo生成多个最终输出文件：记录脚本版本、参数设置、外部工具版本和输入文件路径的文档文件；显示每个物种共表达序列分布的物种计数直方图；提供OCG概述的排序表；每个OCG的Newick格式树文件；以及用于iTOL可视化的注释文件。

概念验证：应用类黄酮生物合成

在花青素、原花青素和黄酮醇三个类黄酮生物合成分支中应用CoExpPhylo，使用相应基因作为输入诱饵，成功检测到多个注释为已知通路成员的OCG以及可能的候选基因。

花青素

聚焦于花青素分支，选择后期酶作为诱饵候选基因，评估CoExpPhylo识别特殊代谢通路基因的能力。输出显示与花青素生物合成中多个关键酶相关的OCG：一个OCG注释为PAL，另一个为C4H，两者都是苯丙烷通路的酶。CHS、CHI、F3H和F3'H被分配到不同的序列簇。初始诱饵基因（DFR、ANS/LDOX和arGSTs）也被检索为OCG。

原花青素

使用后期PA生物合成基因（包括LAR和ANR）作为诱饵，成功检索到PA直接生物合成中所有酶编码基因的注释OCG。与密切相关的通路（如arGST和FLS）相关的基因未被检索，表明CoExpPhylo在区分代谢分支方面具有一定特异性。

黄酮醇

以FLS为诱饵序列研究黄酮醇分支，检索到所有黄酮醇生物合成通路核心基因的注释OCG。还发现了含有通常与花青素生物合成相关基因注释序列的额外OCG，如DFR和ANS/LDOX，表明黄酮醇和花青素分支在共表达分析中不能完全区分。

概念验证：应用类胡萝卜素生物合成

将CoExpPhylo应用于类胡萝卜素生物合成通路，聚焦叶黄素和玉米黄质分支，使用下游生物合成基因作为输入诱饵，流程成功检测到多个与已知酶对应的OCG以及可能参与通路的候选基因。

叶黄素

聚焦类胡萝卜素生物合成的叶黄素分支，使用LUT1/CYP97C1作为诱饵，所有已知叶黄素通路基因（PSY、PDS、ZDS、CRTISO、β-LCY、ε-LCY、LUT5/CYP97A3和LUT1/CYP97C1）都被恢复。

玉米黄质

聚焦玉米黄质分支，恢复了所有已知生物合成基因（除CRTISO和BCH2外）。CoExpPhylo还识别了LUT5/CYP97A3、LUT1/CYP97C1和ε-LCY，这些主要与叶黄素生物合成相关。

工具基准测试

CoExpPhylo允许用户选择MAFFT或MUSCLE进行步骤5中的多序列比对。评估比对方法对系统发育分析的影响发现，两种工具在整体进化结构上一致，但在详细系统发育距离上存在差异。不同建树工具的比较显示更大变异性，尽管整体树拓扑保持一致，但拓扑细节因不同建树方法而异。

建树软件的选择显著影响运行时间性能。使用MUSCLE代替MAFFT进行全局比对步骤导致运行时间增加9.5或12.5分钟。用另一种系统发育树推断工具替代FastTree对计算时间的影响更大：使用IQ-TREE分别增加超过三和四小时的计算时间。MAFFT与RAxML-NG组合执行需要超过140小时，接近六天。

并行处理减轻了大OCG对运行时间的影响。由于所有OCG同时处理，最大簇决定了总计算时间。对于IQ-TREE和RAxML-NG，这种影响特别明显，最大OCG（包含212个序列）分别占总运行时间的77.7%和99%。考虑到数据集中有20个OCG各自包含超过1,000个序列，较大数据集的计算需求将大幅增加。

MAFFT用于全局比对和FastTree用于系统发育树构建的组合代表了最有效的方法。鉴于不同方法获得的生物学结果没有实质性差异，但替代工具的计算时间随OCG大小急剧增加，使用MAFFT和FastTree确保了准确性和性能之间的最佳平衡。

研究表明，CoExpPhylo能够基于每个通路中的单个诱饵基因发现多个通路，跨多个物种识别有前景的特殊生物合成候选基因，为进一步功能和比较分析提供了坚实基础。通过整合共表达分析与系统发育关系，该流程能够从更广泛的、进化 informed的视角理解通路组织和基因功能。

CoExpPhylo的一个关键优势是即使在共表达模式分散的情况下也能识别功能相关基因，这使其成为扩展超出充分表征模式物种的通路注释的强大工具。生成的OCG为突出可能参与生物合成通路的候选基因提供了数据驱动的基础，独立于先前的通路知识。这些候选基因应作为进一步实验验证的起点。

虽然CoExpPhylo有效解析了许多通路组件，但仍有进一步改进的机会。例如，大型基因家族（如MYB转录因子）可能表现出高序列相似性，导致在某些情况下聚类较宽。此外，在某些情况下，来自密切相关的代谢通路的基因被分组在同一OCG中。这突显了需要额外的功能验证来精确确定通路特异性，特别是对于多功能酶或调控蛋白。

注释质量取决于使用的参考数据库，这可能限制阿拉伯芥中缺失基因的分类。将参考数据扩展到包括多个物种，并改进高度同源序列的聚类方法，可以进一步增强功能基因组的解析度。此外，整合蛋白质-蛋白质相互作用预测器将增强程序对于形成metabolon的通路的功能性。准确的蛋白质-蛋白质相互作用预测可以强调物理相互作用的可能性，从而支持识别的候选基因的功能相关性。

尽管存在这些考虑因素，CoExpPhylo代表了生物合成通路探索的宝贵工具，能够识别保守和谱系特异性基因。通过促进跨多样植物物种的候选基因发现，它为采用系统性和可扩展方法研究特殊代谢提供了新机会。

该研究发表在《BMC Genomics》期刊，为植物特殊代谢研究提供了重要的计算生物学工具，填补了共表达分析与系统发育整合分析领域的空白，对未来植物代谢工程和合成生物学研究具有重要推动意义。