来源:Prediction of gene expression levels in Saccharomyces cerevisiae based on chromatin accessibility using multiple machine learning models(Computational Biology and Chemistry)| 编译:DNA Lab Space | 原文许可:elsevier-subscription
导读
染色质可及区通常与转录因子结合和基因激活相关,但利用其序列特征预测基因表达的研究在酿酒酵母中仍属空白。本研究开发了Yeast-Gene——一个基于染色质可及区k-mer特征的监督机器学习模型,聚焦数百碱基对的局部序列,AUC达0.90。可解释性分析锁定AAGAA和CAAGA为高影响motif,二者可能与mRNA剪接相关。这些预测特征有望服务于合成生物学中高表达调控元件的理性设计。
研究背景
染色质状态在基因表达调控中扮演关键角色。染色质呈开放或闭合构象,通常分别对应基因表达的活跃或抑制。可及区更易被转录因子等调控蛋白结合,常与活跃表达共存;不可及区则限制因子结合、降低表达水平。酿酒酵母是重要真核模式生物和合成生物学底盘,但染色质可及性与基因表达水平关联的研究相对有限,这为通过可及区序列特征预测表达留下了空间。理解这一关系对解码酵母转录调控网络、优化生物合成中的基因组设计至关重要。当前基于DNA序列特征预测转录调控网络的研究多集中于启动子、增强子等非编码区,对染色质可及区的明确关注有限。例如DeepSEA等工具主要学习调控序列特征,但并非针对酵母可及区设计。
研究方法
2.1 数据来源
本研究使用的高通量测序数据(Hi-C、ATAC-seq和RNA-seq)来自本实验室,已存入美国国家生物技术信息中心(NCBI)序列读取存档(SRA),登录号为PRJNA1073072。所包含的六株酿酒酵母为yKJP020(H1)、yKJP058(H2)、J005(H3)、JCR27(H4)、BY4741(H5)和yKJP048(H6)。BY4741为野生型。yKJP020、yKJP048和yKJP058为经随机基因组重排的人工改造菌株,在多种环境胁迫条件下筛选获得。JCR27和J005为实验室筛选改造后用于生物合成的菌株。针对模型的下游应用,相关数据集已从公共数据库Gene Expression Omnibus(GEO)检索获得,登录号为GSE101290。该数据集涵盖酵母细胞周期不同阶段和多种营养条件的实验数据,包括相应的ATAC-seq和RNA-seq数据(GOWANS et al., 2018)。与野生型相比,工程菌株在染色质可及区和基因表达模式上表现出显著差异(GUO et al., 1997; PAN et al., 2011; TOCK and HENDERSON, 2018; AMES et al., 2010; BERGSTRÖM et al., 2014)。
【数据】菌株:yKJP020(H1)、yKJP058(H2)、J005(H3)、JCR27(H4)、BY4741(H5)、yKJP048(H6);数据库登录号:PRJNA1073072、GSE101290;数据来源:本实验室及GEO公共数据库
2.2 数据准备
使用MACS2软件对ATAC-seq数据进行peak calling。ATAC-seq中的双端片段堆积反映无核小体区域(NFRs),染色质可及区可通过NFRs处短片段积累来检测。Peak calling后,提取peak区间内的序列。此外,对RNA-seq数据进行基因表达定量分析。为标准化测序深度、基因长度和不同样本间的变异,使用每百万映射读取中每千碱基转录本的片段数(FPKM)对数据进行标准化,增强数据可比性和可解释性。为染色质可及区的基因表达设定阈值(FPKM > 100为高表达,FPKM ≤ 100为低表达)。FPKM值高于100的数据定义为正集,低于或等于100的定义为负集。预测模型的核心输入基于序列k-mer分析。在生物信息学中,k-mer指生物序列中k个连续核苷酸(或氨基酸)的子序列。长度为L的核酸序列,每次滑动一个碱基,产生(L-K+1)个k-mer,此外还可从核酸序列的反向互补链产生额外k-mer。该算法广泛应用于基因组学和蛋白质组学应用,作为分析的基本单元,具有多重优势(BREITWIESER et al., 2019; LEGGETT et al., 2013)。已有多种k-mer计数方法报道(MANEKAR and SATHE, 2018)。我们统计了从单核苷酸到六核苷酸的频率,覆盖k = 1至k = 6。为系统处理和分析这些数据,我们将所有k-mer频率统计存储在精心设计的数据集中,不仅包含不同长度k-mer的频率计数(MANEKAR and SATHE, 2019),还按类型分类数据,具体从开放和非开放区域序列中提取1–5mer信息以构建数据集。开放区域的k-mer频率信息定义为正集,非开放区域的定义为负集。
【数据】软件:MACS2;表达阈值:FPKM > 100(高表达/正集),FPKM ≤ 100(低表达/负集);k-mer范围:k = 1至k = 6;提取范围:1–5mer信息;标准化方法:FPKM
2.3 K-mer频率计算方法
在马尔可夫模型框架下,k-mer可理解为反映不同序列间符号的相对依赖性(YOON, 2009)。对于长度为n的序列S,记为S = (S1, S2, …, Sn),由字母表{A, T, C, G}生成,我们使用k-1阶马尔可夫链模型。此处,当前字符的概率依赖于前k-1个字符。概率公式如下:(1)P(s_i | s_{i-(k-1)}, ..., s_{i-1})。这表示第i个字符的概率取决于前k-1个字符。通过这种依赖关系,可计算整个序列的联合概率:(2)P(S) = P(s1) × P(s2 | s1) × P(s3 | s2) ... P(sn | s_{n-(k-1)}, ..., s_{n-1})。该表达式确定序列S的总体概率,从第一个字符的初始概率开始,通过后续字符的条件概率逐步计算。计算生物序列中所有不同k-mer的出现频率是许多生物信息学应用的关键步骤(MOECKEL et al., 2024)。
【数据】模型:k-1阶马尔可夫链;字母表:{A, T, C, G};公式:P(s_i | s_{i-(k-1)}, ..., s_{i-1});序列长度:n
2.4 数据集分割
为构建模型训练的样本集,我们设计了正负样本划分标准。具体而言,对于正样本,我们关注通过ATAC-seq分析鉴定的染色质可及区中FPKM值超过100的数据。鉴于不同位点染色质可及区的peak宽度各异,我们选取以开放区间peak点为中心、上下游各300 bp的序列,共600 bp DNA序列标记为正样本。该策略得到大量前期试验支持,表明此范围内的序列代表染色质可及性的局部或焦点区域(BUENROSTRO et al., 2013)。基于此标准,我们收集了1511个正样本,每个对应一个特定的开放区域基因。构建负样本集时,我们采用确保负因子真实性和多样性的策略。步骤包括:首先,系统搜索每条染色体中染色质可及区内FPKM值低于100的数据,以排除潜在表达干扰。随后,在这些非基因区域中,按相同位置规则(染色质可及区peak点上下游各300 bp)识别并提取连续序列片段,生成600 bp序列。这些序列被视为低表达基因区域,适合作为模型训练的负样本。共获得2834个负样本。为确保正负样本平衡并维持模型训练有效性,我们从上述候选序列中随机选取1511条,与正样本数量匹配,标记为负样本。在整个样本选择过程中,特别注意保持正负样本平衡,避免模型学习偏差。此外,所有序列均经过严格质量过滤,包括去除低复杂度序列、重复序列和潜在污染物,以确保最终数据集的纯度和可靠性。
【数据】正样本数:1511;负样本候选数:2834;最终负样本数:1511;序列长度:600 bp(peak点上下游各300 bp);正样本FPKM阈值:> 100;负样本FPKM阈值:< 100
2.5 模型构建
按照上述数据准备,从数据中提取k-mer特征。我们开发了Yeast-Gene,一个使用八种不同算法的二分类模型:随机森林(RF)、支持向量机(SVM)、梯度提升机(GBM)、逻辑回归(LR)、分类与回归树(CART)、K近邻(KNN)、朴素贝叶斯(NB)和AdaBoost。选择这八种分类器进行二分类是因为它们各有优势:逻辑回归提供计算效率、概率输出和高可解释性;朴素贝叶斯训练快速且在高维设置中稳健;KNN作为非参数方法,无需显式建模即可捕获局部决策结构;SVM通过核技巧有效处理非线性决策边界,在高维空间中泛化良好;CART通过单棵决策树提供模型简洁性和可解释性;基于树的集成方法……
【数据】模型名称:Yeast-Gene;分类器数量:8种(RF、SVM、GBM、LR、CART、KNN、NB、AdaBoost);任务类型:二分类;特征类型:k-mer特征

研究结果
3.1 Yeast-Gene的性能
#### 3.1.1 Yeast-Gene概述
为探索染色质可及区序列与基因表达的关联,我们设计了一种基于六个不同菌株基因表达FPKM水平的方法(图1和图S1)。我们根据FPKM值将基因分为五个表达水平类别(CUI et al., 2022; WU et al., 2022):未检测到(FPKM < 0.1)1.08%,极低(0.1 ≤ FPKM < 1)1.86%,低(1 ≤ FPKM < 10)13.92%,中(10 ≤ FPKM < 100)59.12%,高(FPKM ≥ 100)24.02%。FPKM = 100阈值超过所有基因的76.0%,位于表达分布的前24.0%。基于表达水平分布(图S2),我们设定FPKM ≥ 100为高表达阈值。K-mer特征已被应用于识别转录因子结合位点(ZHU et al., 2009)、发现序列特异性(HARBISON et al., 2004)、研究TATA框(BASEHOAR et al., 2004)以及识别物种(CSERHATI et al., 2019)。K-mer特征能有效捕获序列模式,从而更准确地定位调控区域并提取基因表达信息。随后我们将基因分为高表达和低表达组,分析各组对应染色质可及区内k-mer序列频率。目标是揭示不同基因表达水平与k-mer序列频率之间的潜在相关性。我们汇总了六个菌株的数据,采用严格验证策略,包括80%/20%数据划分和十折交叉验证,确保模型可靠性。我们进一步测试了k-mer长度从1到5,确认5-mer特征取得最佳预测结果(图S3和补充表1)。
【数据】表达分类:未检测到(FPKM<0.1)1.08%、极低(0.1≤FPKM<1)1.86%、低(1≤FPKM<10)13.92%、中(10≤FPKM<100)59.12%、高(FPKM≥100)24.02%;高表达阈值:FPKM≥100(超过76.0%基因,前24.0%);最佳k-mer长度:5-mer;数据划分:80%/20%;验证方法:十折交叉验证

3.2 预测任务中算法性能的综合评估
为全面评估不同算法在预测任务中的性能,我们选择了八种不同分类器进行模型训练和预测。这些分类器包括常见机器学习算法:支持向量机(SVM)、梯度提升机(GBM)、朴素贝叶斯(NB)、分类与回归树(CART)、K近邻(KNN)、逻辑回归(LR),以及集成方法如随机森林(RF)和AdaBoost(AB)。我们还测试了几种为哺乳动物开发的深度学习模型(如Enformer、DeepSEA),其预测性能表现不佳。我们推测这可能是因为酵母与哺乳动物在基因组大小、染色质结构、调控元件分布和非编码区特征方面存在显著差异,使得适用于哺乳动物的深度学习模型不适用于酵母。我们对各分类器预测基因表达水平的性能进行了详细比较。其中,RF分类器表现最佳,在测试集上达到AUC为0.90(图2A),表明预测性能强劲。SVM和AB分类器也表现良好,AUC分别为0.89和0.87。我们还比较了三个FPKM阈值(50、100和150;补充表2)的性能。100 FPKM值表现最佳,验证了我们参数设置的敏感性。为进一步评估模型的泛化能力和稳健性,我们计算了训练集和测试集上的准确率、灵敏度、特异度和F1分数。这些多维指标提供了模型性能的全面评估(图2)。结果展示了八种分类器在测试集上的主要性能指标。总体而言,RF和SVM模型在测试集上表现尤为出色。
【数据】最佳分类器:RF(AUC=0.90);SVM:AUC=0.89;AB:AUC=0.87;FPKM阈值比较:50、100、150(100最佳);评估指标:准确率、灵敏度、特异度、F1分数;深度学习模型:Enformer、DeepSEA(表现不佳)
3.3 Yeast-Gene各分类器的评估指标
结果表明,RF分类器在测试集上保持高且相对一致的性能,达到AUC为0.90、准确率为0.82、F1分数为0.84(图2和补充表3)。我们使用的k-mer特征与基因表达相关。这些特征是高维的,可能包含非线性模式,使RF适合对其建模。SVM表现略逊于RF,AUC为0.89、准确率为0.79、F1分数为0.80。其相对较强的性能可能归因于SVM能有效建模高维特征空间中的非线性决策边界。其他分类器如AdaBoost和GBM表现相对较弱,特别是在处理非线性和高维特征时。这可能是由于在当前样本量下,它们处理高维相关特征的能力有限。模型一致性分析旨在评估训练集和测试集之间性能是否保持稳定。根据测试集结果,RF保持高且相对一致的性能,与其强劲的训练表现匹配,反映低过拟合和可靠的内部一致性。相比之下,SVM在测试集上的性能略低于RF,特别是在准确率和F1分数方面。这可能是由于其捕获复杂特征交互的能力相对有限。尽管如此,其训练-测试性能差距较小,表明它也避免了严重过拟合。AdaBoost和GBM在测试集上表现相对较差,表明它们在高维特征空间中或样本量有限时可能不稳定。
【数据】RF:AUC=0.90、准确率=0.82、F1=0.84;SVM:AUC=0.89、准确率=0.79、F1=0.80;AdaBoost和GBM:表现相对较弱;模型一致性:RF训练-测试差距小,过拟合低
3.4 Yeast-Gene的稳健性
在机器学习中,向训练数据引入噪声是评估模型泛化和稳健性的常用策略。引入噪声后,几乎所有分类器的性能(以AUC表示)均出现一定下降(图3A)。为获得更全面的评估,我们还使用额外指标评估模型性能,包括准确率和F1分数(图3B)。这些结果提供了多维视角,允许评估各种数据场景并减少对特定数据集特征的依赖。基于这些评估,RF和SVM相比其他分类器表现出更好的稳定性和稳健性。这两种模型以较低的过拟合倾向和对不同数据分布的增强适应性而著称。
【数据】噪声类型:随机SNP诱导噪声;评估指标:AUC、准确率、F1分数;稳健性最佳:RF和SVM;性能变化:引入噪声后几乎所有分类器AUC均下降

3.5 短序列在Yeast-Gene上的性能
为研究是否可使用短于600 bp的序列稳健预测基因表达,我们使用不同长度序列(100 bp、200 bp、300 bp)进行建模和预测(图4)。预测结果显示,使用200 bp序列训练的模型(AUC = 0.89)优于使用100 bp(AUC = 0.86)和300 bp(AUC = 0.87)序列训练的模型。尽管使用200 bp序列的模型性能略逊于使用600 bp序列的模型(AUC = 0.90),这些结果表明即使短序列也能支持稳健预测。这表明高预测性序列特征已存在于染色质可及peak中心上下游100 bp范围内(总长度:200 bp)。尽管300 bp序列包含更多碱基对,其整体预测性能仍然较差,可能是因为300 bp序列引入了对预测信息量较少的区域,导致关键信号稀释或核心调控区域覆盖不理想。总体而言,增加序列长度可提高预测性能,但也带来过拟合风险。较长序列为基因表达模型提供更多预测信息,但短序列也有重要贡献。
【数据】序列长度比较:100 bp(AUC=0.86)、200 bp(AUC=0.89)、300 bp(AUC=0.87)、600 bp(AUC=0.90);最佳短序列:200 bp;核心预测区域:peak中心上下游100 bp(总长200 bp)

3.6 评估Yeast-Gene的泛化能力
机器学习模型的泛化能力取决于多个因素,包括模型复杂度、特征选择策略和样本量。为评估模型在不同数据集上的泛化能力,我们使用了来自酵母代谢周期不同阶段的独立数据集(GSE101290)进行外部验证(图7)。结果显示,RF和SVM在外部数据集上仍保持较好的预测性能,表明模型具有一定的跨数据集泛化能力。
【数据】外部验证数据集:GSE101290;外部验证最佳模型:RF和SVM;评估指标:AUC、准确率、F1分数

3.7 模型可解释性分析
为理解模型预测的序列特征基础,我们进行了特征重要性分析。图5A展示了预测分类结果最重要的前20个k-mer序列motif排名。其中AAGAA和CAAGA被识别为对基因表达预测高度有影响的motif,二者均可能与mRNA剪接相关。图5B展示了使用LIME对代表性样本进行的实例级模型可解释性分析。图5C展示了全局特征关系的部分依赖分析。
【数据】高影响motif:AAGAA、CAAGA;分析方法:特征重要性排名(top 20 k-mers)、LIME局部解释、部分依赖分析;潜在功能关联:mRNA剪接

3.8 特征统计分析
图6展示了正负样本间前五个判别性序列特征的统计比较。图6A通过箱线图展示特征值分布,图6B展示均值,图6C展示标准差,图6D展示每个特征的ANOVA分析结果。
【数据】分析特征数:前5个判别性序列特征;统计方法:箱线图分布、均值、标准差、ANOVA;比较组:正样本 vs 负样本

讨论与解读
随着生物信息学的快速发展,众多创新研究方法(MBOCK et al., 2025)和计算方法(ULLAH et al., 2023; OLUWAFEMI et al., 2024)不断涌现并广泛应用于各个领域。为探索染色质可及区序列与基因表达的关联,本研究提出了名为Yeast-Gene的模型,基于酿酒酵母染色质可及区的k-mer特征预测基因表达。与当前大多数聚焦哺乳动物的表达预测模型(AVSEC et al., 2021)不同,我们的模型靶向真核模式生物酿酒酵母。使用k-mer特征适合捕获酵母研究中与TF结合位点相关的序列模式(KATO et al., 2004),也有助于解决模型训练中数据量有限的问题(GHANDI et al., 2014)。为确保模型可靠性,我们进行了全面的机器学习基准测试:在统一数据预处理和划分下,系统比较了八种分类模型,并使用五个标准指标(AUC、准确率、灵敏度、特异度和F1分数)评估其性能。训练数据整合了多株酿酒酵母,包括野生型和人工工程菌株,实现了多菌株验证。结果表明,模型在不同遗传背景下保持稳定的预测性能,展现出稳健性。
编译者解读
该研究将染色质可及性序列与机器学习结合,为酵母基因表达预测提供了不依赖完整转录调控知识的实用工具。方法选择上,k-mer特征配合随机森林在小样本条件下表现稳健,避免了深度学习对大数据的需求;200 bp短序列即可达到接近600 bp的预测效果,降低了合成生物学中调控元件设计的序列合成成本。局限在于模型仅基于序列特征,未整合表观修饰、三维基因组互作等信息,且外部验证数据集规模有限。AAGAA和CAAGA motif与mRNA剪接的潜在关联值得进一步实验验证,若证实则可直接用于高表达元件的理性设计。
参考来源
期刊:Computational Biology and Chemistry, 2026
DOI: 10.1016/j.compbiolchem.2026.109015
DOI: 10.1016/j.compbiolchem.2026.109015