《Molecular Horticulture》：Integrative machine learning and gene regulatory network analysis identifies novel genes for leafy head formation in Brassica rapa and B. oleracea
编辑推荐：
叶球形成是芸薹属（Brassica）作物中的关键发育过程。研究人员采用一种结合机器学习（ML）与基因调控网络（GRN）分析的整合策略，以鉴定参与大白菜（Brassica rapa）和甘蓝（Brassica oleracea）叶球形成的新基因。以47个已知叶球相
叶球形成是芸薹属（Brassica）作物中的关键发育过程。研究人员采用一种结合机器学习（ML）与基因调控网络（GRN）分析的整合策略，以鉴定参与大白菜（Brassica rapa）和甘蓝（Brassica oleracea）叶球形成的新基因。以47个已知叶球相关基因为训练集的随机森林（Random Forest）模型表现出稳健性能，在B. rapa和B. oleracea中的平均曲线下面积（AUC）分别达到0.87和0.85。通过对模型预测结果进行筛选，研究鉴定出11个在两个物种间共享且具有高置信度的候选基因。为进一步揭示其调控机制，研究人员分别构建了两个物种的基因调控网络。将ML预测结果与这些网络整合后，研究鉴定出与叶球形成特异相关的关键调控簇。网络中心性分析显示，在关键簇中存在多个核心基因，其中包括ANT、GRF2、GRF3和TCX3等重要转录因子，提示这些基因在两个物种的叶球形成过程中具有关键作用。在两个物种中平行检测到相同基因及相似网络结构，支持了研究结果的可靠性。该整合策略为叶球形成的遗传调控提供了新的认识，并为后续芸薹属物种的功能研究奠定了基础。
该研究发表于《Molecular Horticulture》，聚焦于芸薹属蔬菜叶球形成这一重要农艺性状的遗传调控机制。叶球形成决定了大白菜和甘蓝等作物的商品性、产量潜力与消费者偏好，其形成过程涉及叶片生长、叶片极性建立、分生组织维持、激素调控以及细胞扩张等多个层面的协同作用。既往研究已通过全基因组重测序、转录组分析和miRNA表达分析等手段鉴定出部分参与叶球形成的基因，并指出生长素、叶片近轴—远轴极性以及分生组织活性是该过程的重要调控因素。然而，叶球形成并非由单一基因或单一路径控制，而是依赖复杂的多基因调控网络。当前的主要问题在于：尽管已有若干关键基因和通路被报道，但它们之间如何在系统层面被整合，如何共同塑造叶球发育这一复杂性状，仍不清楚。与此同时，芸薹属物种已积累大量基因组和转录组数据，传统分析方法难以充分挖掘其中潜在的高维信息，因此亟需引入能够处理复杂多组学特征并解析调控关系的计算方法。
基于这一背景，研究人员构建了一种整合机器学习与基因调控网络分析的研究框架，分别在B. rapa和B. oleracea中平行开展分析，以识别叶球形成新基因并解析其潜在调控模块。之所以选择两个物种进行并行研究，是因为这种设计既可以比较不同物种中的共同调控机制，又能够通过结果一致性增强预测可信度，还可用于观察叶球类群形成过程中可能存在的保守性与分化特征。研究最终表明，随机森林模型能够有效预测新的叶球形成候选基因，而与基因调控网络结合后，可进一步定位与叶球形成密切相关的关键调控簇和核心转录因子。研究还发现，两个物种在关键网络模块、富集生物过程和部分核心调控因子上表现出明显一致性，说明叶球形成在不同芸薹属物种中具有一定保守的分子基础。该研究的重要意义在于，不仅扩展了叶球形成相关候选基因库，而且建立了一个适用于复杂园艺性状解析的计算分析范式，为后续功能验证与分子育种提供了高价值靶标。
在技术方法上，研究主要采用了以下几类关键策略。首先，研究人员整合来源于公共数据库的RNA-seq原始数据，分别获得B. rapa的634个标准化表达特征和B. oleracea的762个标准化表达特征。其次，以47个文献报道的叶球形成相关基因为正类样本，并依据表达中位数匹配策略筛选负类样本，构建随机森林分类模型，采用五折交叉验证评估模型性能。再次，利用PlantTFDB转录因子集及标准化表达矩阵，通过GENIE3推断转录因子—靶基因调控关系，并借助iGraph进行网络聚类。最后，结合SynOrths共线性基因配对、Jaccard指数、GO富集分析和CytoNCA中心性分析，系统比较两个物种中叶球形成相关网络的保守性与核心节点。
在结果部分，研究首先提出了整体研究流程，并通过四个步骤完成叶球形成相关基因与调控网络的识别。该流程包括知识库构建、随机森林预测、基因调控网络构建以及两个物种之间的比较分析，为后续结果提供了系统框架。
在“Model construction for predicting genes involved in leafy head formation in B. rapa and B. oleracea”部分，研究人员基于两个物种的大规模转录组特征和47个已知叶球相关基因建立模型。为避免正负样本表达水平差异带来的偏倚，研究以表达分布相似为原则构建负类基因集合。模型评估显示，B. rapa模型的平均AUC为0.87 ± 0.08，B. oleracea模型的平均AUC为0.85 ± 0.04，说明两套模型均具有较强的分类性能。进一步的特征重要性分析采用平均精度下降值（MDA）评估不同组织表达特征的贡献，结果表明花和角果组织在两物种间存在显著差异，而叶、根、茎端分生组织（SAM）和幼苗组织差异不显著，提示这些组织在叶球相关基因预测中可能具有相似的重要性。偏依赖图（PDP）分析进一步显示，根组织表达大多对预测呈负向影响，而叶组织的影响方向则更为复杂。
在“Identification and functional analysis of predicted genes involved in leafy head formation”部分，模型被应用于全基因组范围的基因预测。研究在B. rapa中鉴定出5,201个预测候选基因（preGenes），在B. oleracea中鉴定出8,101个。参数敏感性分析表明，不同树数量设置下预测结果具有高度一致性，说明模型结果稳定。基于物种间共线性关系，研究发现有2,614对共线性基因在两个物种中均被预测为候选基因，且二者预测概率呈显著正相关。进一步筛选显示，共有11个高概率候选基因（概率>0.9）值得重点关注。对所有preGenes进行GO富集分析后，两个物种共同富集于单向细胞生长、近轴/远轴模式指定、生长素运输和分生组织发育调控等过程，这些过程均与叶球形成密切相关。研究还利用Guo等人提供的大白菜结球期24个空间分割组织转录组数据，对11个高概率基因的表达模式进行了验证，发现这些基因可分为内叶高表达型、SAM偏好表达型和普遍低表达型三类，从空间表达层面支持其参与叶球形成的可能性。留一交叉验证进一步证明，该方法能够较稳定地重新识别已知叶球相关基因。
在“Construction of gene regulatory networks for transcription factors in B. rapa and B. oleracea”部分，研究人员分别从PlantTFDB获得2,840个B. rapa转录因子和2,918个B. oleracea转录因子，并基于转录组数据构建全局基因调控网络。经过阈值筛选后，B. rapa网络包含332,429条互作关系，B. oleracea网络包含382,103条互作关系。进一步聚类后，前者被划分为10个子簇，后者被划分为9个子簇，表明叶球形成相关调控关系具有明显模块化结构。
在“Integration of machine learning predictions and gene regulatory networks”部分，研究将ML预测得到的preGenes映射到上述转录因子调控网络中。结果显示，B. rapa中有4,856个preGenes进入GRN，其中473个为转录因子；B. oleracea中有7,597个preGenes进入GRN，其中724个为转录因子。提取preGenes特异调控关系后，B. rapa子网络保留18,792条互作，B. oleracea子网络保留44,139条互作。对已知基因与高概率候选基因的分布进行考察发现，B. oleracea的cluster 7和B. rapa的cluster 8均富集最多的已知叶球相关基因，并同时容纳全部高概率候选基因，因此被界定为叶球形成的关键调控模块。依据这些簇中基因的已知功能，研究进一步指出其主要涉及叶片极性建立、分生组织调控及生长素介导的生长三个核心过程。
在“Comparative analysis of heading network in B. rapa and B. oleracea”部分，研究利用Jaccard指数比较两个物种GRN中各簇的共线性候选基因组成。结果显示，B. oleracea cluster 7与B. rapa cluster 8之间的Jaccard指数最高，为0.20，提示这两个关键簇在基因组成和功能上具有较高保守性。两簇均包含YAB1.1、YAB1.2、AS1.2、PHB.1、PHB.2、KNAT6.1和KNAT6.2等已知叶球相关基因，进一步支持其关键地位。GO富集分析显示，这两个关键簇在近轴/远轴模式指定、轴向指定、表皮发育和生长素运输等功能上高度一致，而部分其他簇则在光合作用光反应、光强和碳水化合物代谢等方面表现出相似性。为识别关键簇中的核心调控因子，研究进一步开展网络中心性分析，综合介数中心性、接近中心性、度中心性和特征向量中心性四项指标筛选核心基因。结果表明，B. oleracea cluster 7中鉴定出81个核心基因，其中39个为转录因子；B. rapa cluster 8中鉴定出50个核心基因，其中20个为转录因子。若干核心基因在两个物种中具有共线性对应关系，提示其可能具有保守调控作用。尤其是ANT、GRF2、GRF3和TCX3等转录因子在两个物种关键簇中均被识别，说明这些因子可能是叶球形成共享调控网络中的核心节点。
讨论部分围绕模型构建策略、组织特异性特征、候选基因功能以及网络保守性与分化展开。研究指出，采用与正类基因表达水平相匹配的负类样本筛选方法，有助于减少训练偏倚，使模型更可能捕捉真实生物学差异。组织特征分析提示，叶、幼苗和根组织在两个物种中对预测的重要性较为相似，而根组织普遍对叶球形成候选基因预测呈负向贡献，符合叶球形成主要发生于地上部组织的生物学认知。11个高概率候选基因在结球期空间表达上的差异模式，也与其潜在功能相一致。网络分析方面，研究强调两个物种中均存在保守的关键调控模块，同时也存在某些物种特异的功能富集现象，反映出叶球形成机制兼具保守性与分化性。核心基因分析进一步支持GRF2、GRF3、ANT和TCX3等转录因子在叶球形成中的重要作用。总体而言，研究建立的整合分析框架扩展了叶球形成的候选基因谱系，并为复杂园艺性状的遗传解析提供了可迁移的方法学基础。
研究结论部分可概括为：该研究通过整合机器学习与基因调控网络分析，在B. rapa和B. oleracea中成功鉴定出大量新的叶球形成候选基因，并解析出与叶球形成相关的关键调控模块和核心转录因子。两个物种中相似候选基因与网络结构的平行出现，验证了该策略的稳健性，并揭示了叶球形成背后保守的遗传调控机制。这一研究为芸薹属叶球性状的功能基因发掘、分子机制研究和定向育种提供了重要理论基础与候选资源。
