酵母基因组核小体定位序列预测

Prediction of Nucleosome Positioning Sequence for Yeast Genome

摘要: 核小体是染色质结构的基本单位，其在整条DNA序列上的定位分布情况，对于真核生物的基因表达调控起关键作用。用机器学习方法预测核小体定位成为近年来的研究热点。以DNA序列6-mer组分为参数，采用我们提出的多样性增量特征选择技术，筛选出8个6-mer作为分类特征。进一步，采用支持向量机算法，10折交叉检验的总精度达到98.2%。结果表明，核小体定位序列和连接序列核苷k-mer组分的特异化分布，是影响酵母核小体定位的主要因素。

Abstract: Nucleosome is a basic unit of chromatin structure. Its location and distribution on the entire DNA sequence play a key role in the regulation of gene expression in eukaryotes. The prediction of nucleosome positioning with machine learning method has become a hot topic in recent years. Taken the 6-mer component of DNA sequence as the parameter, we used the increment of diversity feature selection technique proposed by us to select eight 6-mers as the classification characteristics. Furthermore, the total accuracy of the 10 fold cross validation is 98.2% using the support vector machine algorithm. The results show that the specific distribution of the k-mer component in the nucleosomal and linker sequences is the main factor that affected nucleosome positioning in yeast.

1. 引言

核小体是真核生物特有的DNA包装方式，是染色质结构的基本组成单位。约147bp的DNA序列紧密缠绕在组蛋白8聚体上约1.65圈，形成核小体核心颗粒(nucleosome core particle)，这个核心颗粒阻断了大部分蛋白质分子与DNA的接触 [1]。核小体核心颗粒之间由长度不等的连接DNA(linker DNA)相连 [2]。核小体和连接DNA在基因组上的精确定位，对调控基因的表达过程起到关键作用 [3] [4]。

随着对表观遗传学研究热度的提升，基于DNA序列信息或序列依赖的各种性质的核小体定位预测模型不断被提出 [5] - [10]。核小体在基因组中的分布要受到多种因素的影响，Segal组的研究表明 [5] ，DNA的序列组成是基因组中核小体组织的主要决定因素。Chen等通过计算序列依赖的二核苷变形能 [8] ，发现对于酵母基因组，150 bp长的核小体定位序列与连接序列的变形能存在显著差异，基于此，他们以98.1%的精度分类了酵母核小体定位序列与连接序列。如果我们仔细分析酵母基因组核小体定位序列与连接序列的序列组成，就不难发现变形能的这个差异，来源于二者序列组成的差异。最近，Awazu提出一个多元回归模型 [9] ，以3核苷(或其互补)频次为参数，在Chen等使用的数据集上 [8] ，声称分类精度达到100%，但Awazu的方法参数有33个之多 [9] ，且这些参数的取值强依赖性于数据集。2016年，Teif综述了现有的核小体数据资源和在线预测工具 [10]。该文总结了近些年分别从生物信息学，物理学和生物学角度研究核小体定位的近百篇论文，为今后的核小体定位研究提供了一个手册。

采用实验手段测定基因组核小体定位情况耗时耗力，如果能预先由生物信息学方法给出一个基因组核小体定位预测的图谱，可能为实验的实施提供预设的靶标，进而节约实验成本。此外，基于生物信息学方法的核小体定位预测，还能发现隐藏在碱基序列中的规律性信息，为解开生命奥秘提供新的线索。因此，在分子生物学研究中，生物信息学手段是必要的。本文，在Chen等使用的酵母核小体定位数据集上 [8] ，我们基于DNA序列的6核苷组成，采用多样性增量特征选择技术，选出仅8个6-mer为分类参数，使用支持向量机算法进行分类，获得较高的分类精度。本文的创新点在于，采用了我们研究组新近发展的特征选择方法，在保证预测精度的同时，极大地压缩了模型的输入参数个数，使得该模型具有简单易用，参数少，精度高，鲁棒性好的特点。

2. 材料与方法

2.1. 酵母核小体定位序列和连接序列数据集

实验确定的酵母核小体定位数据取自Lee等的实验结果 [11]。这个实验结果中给出4 bp分辨率的1,206,683个DNA片段，采用“套索模型(lasso model)”对每个片段分配一个核小体形成能力得分，得分的高低反映了核小体形成的倾向性高低。Chen等从其中选出5000个长度为150 bp得分最高的片段作为核小体定位序列 [8] ，和5000个长度为150 bp得分最低的片段为连接序列。然后他们使用CD-HIT软件过滤相似序列 [12] ，最后得到序列相似性低于80%的3620个样本，其中1880个正样本，即核小体定位序列，1740个为负样本，即连接序列。用具有过高序列相似性的数据集检验预测模型，可能给出虚高的性能评价结果，本文，我们基于Chen等的结果，并进一步采用BLASTClust (https://toolkit.tuebingen.mpg.de/blastclust)程序，使数据集中的序列相似性低于30%。最后我们得到2377个样本，其中1334个核小体定位序列和1043个连接序列，序列长度均为150 bp。

2.2. 方法

2.2.1. 以序列K-Mer组分为特征

由Segal组的研究结论可知 [5] ，序列组成是基因组中核小体组织的主要决定因素。那么，基于序列的碱基分布来预测核小体定位序列成为当然的想法。当k取值较大时，基于序列k-mer组分分布可能重构整个基因组序列，因此，一条序列的k-mer分布可以作为该序列的充分表示。我们以150 bp长序列中所含的6-mer频次作为序列特征。k = 6对于长度150 bp的序列而言，应该是一个充分表示了。

显然，当k = 6时，提取的特征总数将达到4096个。如果这些特征均输入分类器用于分类，严重的过拟合现象将不可避免。因此，应用特征选择技术实现降维是必要的步骤。在应对小样本或高维数据的模式识别问题时，特征选择是许多机器学习方法的重要组成部分。特征选择技术的主要目标是，从整体特征集中发现一个能够有效描述数据的特征子集。因此，对于从序列提取到的4096个6-mer频次特征，我们要应用特征选择技术对其进行筛选，选择有效分类信息，使模型具有更好的鲁棒性。

2.2.2. 多样性增量特征选择技术

最近，我们研究组发展了一种新的特征选择技术，称为多样性增量特征选择(feature selection based on incremental of diversity, FSID) [13]。在此之前，我们主要将多样性增量方法结合二次判别分析直接用于模式分类 [14] [15]。近些年，众多新的特征选择技术得以发展。2015年，Drotár等比较了10个目前最好的特征选择技术 [16] ，结果表明，信息增益(information gain, IG)方法导致最好的稳定性 [17] ，而最小冗余最大相关(minimum Redundancy Maximum Relevance, mRMR)方法导致最高的预测精度 [18]。FSID方法类似于IG方法，但不同之处在于，FSID方法无需用特征在序列中出现的相对频次去近似特征出现的概率，而直接使用特征的绝对频次。这样的做法对于小样本问题优势明显，不存在理论困难 [13]。因为我们知道，当对一个随机事件仅作有限几次观测后，就以观测的频率代替概率会存在很大风险。下面我们对FSID方法做一简要描述。

给定一个两分类问题C_i(i = 1,2)，特征X出现在类别C₁的频次记为n，除特征X以外的其它特征出现在类别C₁的频次记为 $\bar{n}$。则特征X在C₁类别中的多样性量定义为

$D\left(X,{\text{C}}_1\right)=\left(n+\bar{n}\right){\log }_2\left(n+\bar{n}\right)-n{\log }_2\left(n\right)-\bar{n}{\log }_2\left(\bar{n}\right)$ (1)

类似地，特征X出现在类别C₂的频次记为m，其它特征出现在类别C₂的频次记为 $\bar{m}$ 。按照与(1)式同样的方式，可以分别定义特征X在C₂类别中的多样性量D(X,C₂)，以及在混合系统C₁ + C₂中的多样性量D(X,C₁ + C₂)为

$D\left(X,{\text{C}}_2\right)=\left(m+\bar{m}\right){\log }_2\left(m+\bar{m}\right)-m{\log }_2\left(m\right)-\bar{m}{\log }_2\left(\bar{m}\right)$ (2)

和

$\begin{array}{l}D\left(X,{\text{C}}_1+{\text{C}}_2\right)=\left(n+m+\bar{n}+\bar{m}\right){\log }_2\left(n+m+\bar{n}+\bar{m}\right)\\ -\left(n+m\right){\log }_2\left(n+m\right)\\ -\left(\bar{n}+\bar{m}\right){\log }_2\left(\bar{n}+\bar{m}\right)\end{array}$ (3)

特征X在C₁和C₂类别之间的多样性增量(increment of diversity, ID)定义为

$ID\left(X\right)=D\left(X,{\text{C}}_1+{\text{C}}_2\right)-D\left(X,{\text{C}}_1\right)-D\left(X,{\text{C}}_2\right)$ (4)

由以上多样性增量的定义可以看出，当给定样本集时，某个特征X在C₁和C₂两个类别中出现的频次差异越大，ID(X)的值越大，而频次差异越小，ID(X)的值越小。如果特征X是类别无关的，那么一般地，特征X在两个类别中出现的频次应几乎无差别。因此，ID(X)就可作为特征X是否与类别相关的度量。也就是说，如果ID(X) > ID(Y)，表明特征X与类别相关性要强于特征Y。当这种类别相关性的强度达到我们的预期(ID₀)时，即当ID(X) > ID₀，特征X被选择，ID₀为特征选择阈值。阈值ID₀是以选出的特征使得预测结果总精度最大化来确定的。这个特征选择方案被我们称为多样性增量特征选择技术(FSID)。

2.2.3. 预测算法

预测算法采用支持向量机，算法实现采用R语言“e1071”包中的svm函数完成 [19]。核函数采用径向核函数，参数c和γ分别采用默认值。

2.2.4. 分类性能评价指标

分类性能采用如下4个指标度量，分别是敏感性(Sensitivity, Sn)，特异性(Specificity, Sp)，总精度(Accuracy, ACC)和马氏相关系数(Matthews correlation coefficient, MCC)。这些量定义如下

$\begin{array}{l}\text{Sn}=\frac{\text{TP}}{\text{TP}+\text{FN}},\text{Sp}=\frac{\text{TN}}{\text{TN}+\text{FP}},\text{ACC}=\frac{\text{TP}+\text{TN}}{\text{TP}+\text{FN}+\text{TN}+\text{FP}},\\ \text{MCC}=\frac{\text{TP}\times \text{TN}-\text{FN}\times \text{FP}}{\sqrt{\left(\text{TP}+\text{FP}\right)\times \left(\text{TP}+\text{FN}\right)\times \left(\text{TN}+\text{FP}\right)\times \left(\text{TN}+\text{FN}\right)}}\end{array}$ (5)

这里，TP是被正确预测的核小体定位序列数，FN是被错误预测的核小体定位序列数，TN是被正确预测的连接序列数，FP是被错误预测的连接序列数。

上面4个分类性能评价指标分别表明了一个预报器的四个不同方面的性能。Sn是在全体正样本中能够被正确预测为正样本的频率，它用来衡量一个预报器识别正样本的能力。类似地，Sp是用来衡量一个预报器识别负样本的能力。ACC测度正确识别全部样本的能力。MCC是预测性能的一个最佳平衡测度。MCC的取值范围是[−1,+1]。MCC = 0表明预报器实际执行了一个随机猜测，也即它的预测结果与样本的真实分类标签不相关。MCC = ±1表明预报器是完美的。同时给出一个预报器的4个性能指标，可以较全面地反映出预报器的输出性能。

3. 结果与讨论

3.1. FSID特征选择结果

基于核小体序列中6-mer频数信息，采用FSID方法对酵母核小体定位序列和连接序列进行分类预测特征选择。将全部数据样本随机分割为10份，保证每份中正负集样本数之比大致与正负总样本数之比相当。合并其中的9份样本作为训练样本。在训练样本中，统计所有特征出现的频次，采用FSID方法选择ID值大于阈值ID₀的特征，然后轮换训练样本。如果10次轮换某一特征均被选出，则特征被最终选择。将最终选择出的特征送到SVM中进行核小体定位序列预测。本文中当阈值ID₀ = 615时，预测精度达到

Table 1. Prediction results in 10 fold cross-validation

表1. 10折交叉检验预测结果

最大，此时选取出8个特征参数分别是AAAAAA、AAAAAT、ATATAT、ATTTTT、TAAAAA、TATATA、TTTTTA和TTTTTT。

可以看到，选出的8个特征均为poly(dA:dT)。进一步的分析表明，这些特征均来自连接序列。或者说在酵母基因组中，核小体之间的连接序列中普遍地存在着poly(dA:dT)片段，而核小体定位序列中则罕有。Poly(dA:dT)序列是刚性的，不利于核小体的形成 [20]。这可能是Liu等和Chen等可以成功地基于序列依赖的二核苷弯曲能方法 [7] [8] ，来预测酵母核小体定位的原因所在。

3.2. SVM预测结果

采用10-fold交叉检验，酵母核小体定位序列预测结果列于表1。由表1结果可见，我们的模型对酵母核小体定位序列的预测敏感性(Sn)是99.1%，特异性(Sp)为97.7%，总精度(ACC)为98.2%，马氏相关系数(MCC)值达到0.963。我们的结果与Chen等基于序列依赖的二核苷变形能模型的结果，具有相当的精度 [8]。但Chen等数学模型更为复杂 [8] ，且其使用的数据集的序列相似性更高。在2012年，Chen等还曾提出一个基于氨基酸物化性质的酵母核小体预测模型 [21] ，该模型经特征选择后保留了884个特征作为输入，在与文献同样的数据集上进行10折交叉检验 [8] ，其预测结果的4个性能指标均低于我们的模型，如表1所示。

Awazu提出一个多元回归模型 [9] ，以3核苷(或其互补)频次为参数，在Chen等使用的数据集上 [8] ，声称分类精度达到100%，但Awazu的方法需拟合的参数有33个之多 [9] ，且这些参数的取值强依赖性于数据集，预测精度极不稳定，作者将该方法用于其他物种的核小体预测，发现获得的精度极低 [9]。

从以上的分析看出，我们给出的模型具有参数少且精度高的特点，更少的参数将使得模型具有更高的鲁棒性和泛化能力。这个良好的性能得益于我们采用了高效的特征选择技术FSID。

4. 结论

本文，我们采用多样性增量特征选择技术FSID，对以核苷六联体(6-mer)为参数的特征集进行筛选，筛选出8个poly(dA:dT)特征对酵母核小体定位序列进行分类预测，在序列相似性低于30%的数据集上，10折交叉检验获得98.2%的高精度结果。模型具有数学方法简单，使用参数少，预测精度高的优点。模型给出的结果还表明，酵母核小体之间的连接序列的主要序列组成特点是，存在普遍的poly(dA:dT)片段。这些片段具有很强的刚性，不易弯曲，难于形成核小体结构。

尽管由FSID方法所选择的特征与类别之间显著相关，但FSID方法中并未考虑特征与特征之间的相关性。如何将减少特征之间相关性的算法也融合在现有的FSID模型中，是今后研究中需要解决的问题。

基金项目

本项目由内蒙古自治区自然科学基金项目(2015MS0331和2016MS0306)资助。

NOTES

^*通讯作者。