一种基于核小体脱氧核糖核酸模版的核小体预测方法

摘要

一种基于核小体脱氧核糖核酸模版的核小体预测方法，其特征在于，包括以下步骤：步骤1获取待预测DNA序列，长度为T，计算待预测DNA序列弯曲度信号Signal。步骤2建立核小体DNA模版信号P，P的长度为147bp，两端区域宽为50bp，高为0.07，中间区域宽为47bp，高为0.05，卷积P和Signal得到信号S_covn，从S_covn中部取出长为T-10的信号S_covn_keep。步骤3计算S_covn_keep的连续小波变换W(

说明书

技术领域

本发明涉及一种预测真核基因组核小体(nucleosome)的方法，该预测方法通过匹配核小体DNA的弯曲度模版实现预测。 

背景技术

真核生物DNA以染色质形式存在，其基本单元为核小体；核小体系由~147bp的DNA缠绕在组蛋白八联体(2倍的组蛋白H2A、H2B、H3和H4)上形成的结构，核小体之间以连接DNA相连。核小体定位是指DNA双螺旋相对于组蛋白核的位置，核小体DNA的蛋白结合位点被封闭，因此无法结合蛋白(转录因子)。核小体定位通过遮蔽或者暴露蛋白结合位点调节真核基因的复制、转录等基础生物学过程。因此，预测核小体对于分析基因的转录调节等至关重要，预测核小体就是预测基因组DNA哪些区域是缠绕在组蛋白核上的核小体DNA(nucleosome DNA)，哪些是处于连接区域的连接DNA(linker DNA)。 

利用湿实验手段直接检测核小体仍然需要一定的成本。近年来，研究发现核小体的位置部分地与DNA序列有关，即可以通过基因组DNA序列实现核小体的预测。目前，国际上有一种通过概率模型预测核小体的方法，这种方法首先需要计算四种核苷酸在核小体DNA序列的每个位置出现的频率，形成位点频率矩阵，然后计算待预测序列在该位点频率矩阵上的得分，实现核小体预测。这种方法的缺点在于：1，建立位点频率矩阵，首先需要知道一定数目的核小体DNA序列；2，对不同物种基因组核小体预测，需要建立不同的位点频率矩阵；3，计算方法较复杂，对于大规模的基因组预测，速度显得较慢。这些缺陷限制了这种方法的进一步应用。 

发明内容

本发明提供一种能够实现快速准确预测的基于核小体脱氧核糖核酸模版的核小体预测方法。 

本发明采用如下技术方案： 

一种基于核小体脱氧核糖核酸模版的核小体预测方法，包括以下步骤： 

步骤1获取待预测的DNA序列，核苷酸的长度为T，并利用 

$>C=v^0{(n_2-n_1)}^{-1}\underset{j=n_1}{\overset{n_2}{\Sigma }}({\rho }_j-{\mathrm{i\tau }}_j)\exp \left(\frac{2\mathrm{\pi ij}}{v^0}\right)$>

计算待预测的DNA序列弯曲度信号Signal， 

其中，C的模代表弯曲度，v⁰=10.4bp，是DNA双螺旋一个螺旋的碱基对平均数目，n₁和n₂代表加和范围的起始和终止位置，i为复数单位，ρ和τ是16种二联核苷酸在空间弯曲的幅度，ρ为在平面旋转的幅度，τ为扭曲的幅度，其值见表1。 

表1弯曲度计算公式中的ρ和τ的值 

详细的计算过程为：以一个宽度为10bp的滑动窗口从待预测序列的5’端选取序列，利用上述弯曲度公式，计算所选序列的C值，C的模为所选序列的弯曲度，作为Signal的第一个数值，然后向待预测序列的3’方向平移一个核苷酸，用窗口选却10bp的序列，计算所选序列的C值及其模，作为Signal的第二个数值。向3’方向滑动窗口，重复计算过程，记录C的模，直至待预测序列的末尾，便得到待预测序列的弯曲度信号Signal，其长度为T-10。 

步骤2建立核小体DNA模版信号P，所述的模板信号为一条147(bp)长的数 字信号，该模板信号两端区域宽度为50bp，高度为0.07，中间区域的宽度为47bp，高度为0.05，卷积弯曲度信号和模板信号得到卷积信号： 

给定模板信号P和弯曲度信号Signal，长度分别为147和T-10，二者的卷积为 

$>S\_\mathrm{covn}\left(k\right)=\underset{n=1}{\overset{k}{\Sigma }}P\left(n\right)·\mathrm{signal}(k-n)$>

其中，S_covn为卷积的结果，长度为147+T-10-1，从卷积结果S_covn中部取出长度为T-10的信号，标记为S_covn_keep，长度为T-10，k是一个变量，k∈[1,147+T-10-1]，n=1,2,3,…,k, 

步骤3检测卷积信号S_covn_keep的峰位置，具体过程为： 

步骤3.1计算S_covn_keep的连续小波变换，小波母函数为墨西哥帽函数；尺度范围为[2,8]，记连续小波变换结果为W(a,b)， 

$>W(a,b)=\frac{1}{\sqrt{a}}{\int }_{-\infty }^{+\infty }S\_\mathrm{conv}\_\mathrm{keep}\left(t\right){\psi }_{a,b}\left(t\right)\mathrm{dt}$>

其中，a和b分别为尺度因子和平移因子，Ψ为小波函数，W(a,b)为连续小波变换的结果，其中b∈[1,T-10]，S_conv_keep(t)为S_conv_keep第t个元素，W(a,b)为二维信号，表示尺度为a_i时，i∈[2,8]，当小波函数Ψ的平移量是b时，S_conv_keep为在小波函数空间的投影， 

步骤3.2在每个平移因子b处计算|W(a,b)|的最大值， 

M_W(b)=max(|W(a,b)|) 

其中，|W(a,b)|小波变换的绝对值，M_W(b)为在平移为b时，S_conv_keep在尺度范围[2,8]上的小波函数上投影的最大值，b∈[1,T-10]， 

步骤3.3M_W(b)中的峰即为核小体的二分点Dyad位置，这样便实现了核小体的预测。 

与现有技术相比，本发明具有如下优点 

1，本发明可以适用于多物种基因组核小体的预测，而不需要改变模版信号和其它参数。可以实现高等生物(如人)和简单真核生物(如酵母)基因组核小体预测的，适用范围广。 

2，本发明具有较好的预测准确性，在40bp误差范围内，预测结果优于基于概率的方法。 

3，本发明计算过程简单，因此计算速度较快，可以实现大规模基因组核小体的预测(如人类基因组)。 

附图说明

图1是核小体DNA的弯曲度模版信号P示意图； 

图2是弯曲度谱预测核小体的流程图； 

图3本发明方法对人类基因组20号染色体的一段区域(8000bp-28000bp)的核小体预测； 

图4本发明方法预测的人类20号染色体核小体的结果分析。在不同偏差范围内，本发明预测与实验测定的核小体和《自然》杂志方法预测的核小体比较，X坐标为偏差，Y为在该偏差下核小体的匹配率，匹配率指在此偏差下匹配的核小体的数目占总的实验检测的核小体的数目的百分比，覆盖率为某偏差下两个比较的核小体的重叠程度，当偏差为0bp时，覆盖率为100% 

具体实施方式

一种基于核小体脱氧核糖核酸模版的核小体预测方法，预测原理及流程见图4，包括以下步骤： 

步骤1获取待预测的DNA序列，长度为T，并利用 

$>C=v^0{(n_2-n_1)}^{-1}\underset{j=n_1}{\overset{n_2}{\Sigma }}({\rho }_j-{\mathrm{i\tau }}_j)\exp \left(\frac{2\mathrm{\pi ij}}{v^0}\right)$>

计算待预测的DNA序列弯曲度信号Signal，其中，C的模代表弯曲度，v⁰=10.4bp，是DNA双螺旋一周螺旋周期的碱基对平均数目，n₁和n₂代表加和范围的起始和终止位置，i为复数单位，ρ和τ是16种二联核苷酸在空间弯曲的幅度，ρ为在平面旋转的幅度，τ为扭曲的幅度，其值见表1。 

表1弯曲度计算公式中的ρ和τ的值 

详细的计算过程为：以一个宽度为10bp的滑动窗口从待预测序列的5’端选取序列，利用上述弯曲度公式，计算所选序列的C值，C的模为所选序列的弯曲度，作为Signal的第一个数值，然后向待预测序列的3’方向平移一个核苷酸，用窗口选却10bp的序列，计算所选序列的C值及其模，作为Signal的第二个数值。向3’方向滑动窗口，重复计算过程，记录C的模，直至待预测序列的末尾，便得到待预测序列的弯曲度信号Signal，其长度为T-10。 

步骤2建立核小体DNA模版信号P，所述的模板信号为一条147(bp)长的数字信号，该模板信号两端区域宽度为50bp，高度为0.07，中间区域的宽度为47bp，高度为0.05，卷积弯曲度信号和模板信号得到卷积信号S_covn： 

给定模板信号P和弯曲度信号Signal，长度分别为147和T-10，二者的卷积为S_covn： 

$>S\_\mathrm{covn}\left(k\right)=\underset{n=1}{\overset{k}{\Sigma }}P\left(n\right)·\mathrm{signal}(k-n)$>

其中，S_covn为卷积的结果，长度为147+T-10-1，从卷积结果S_covn中部取出长度为T-10的信号，标记为S_covn_keep，长度为T-10，k是一个变量，k∈[1,147+T-10-1]，n=1,2,3,…,k, 

步骤3检测卷积信号S_covn_keep的峰位置，具体过程为： 

步骤3.1计算S_covn_keep的连续小波变换，小波母函数为墨西哥帽函数；尺度范围为[2,8]，记连续小波变换结果为W(a,b)， 

$>W(a,b)=\frac{1}{\sqrt{a}}{\int }_{-\infty }^{+\infty }S\_\mathrm{conv}\_\mathrm{keep}\left(t\right){\psi }_{a,b}\left(t\right)\mathrm{dt}$>

其中，a和b分别为尺度因子和平移因子，Ψ为小波函数，W(a,b)为连续小 波变换的结果，其中b∈[1,T-10]，S_conv_keep(t)为S_conv_keep第t个元素，W(a,b)为二维信号，表示尺度为a_i时，i∈[2,8]，当小波函数Ψ的平移量是b时，S_conv_keep在小波函数空间的投影， 

步骤3.2在每个平移因子b处计算|W(a,b)|的最大值， 

M_W(b)=max(|W(a,b)|) 

其中，|W(a,b)|小波变换的绝对值，M_W(b)为在平移为b时，S_conv_keep在尺度范围[2,8]上的小波函数上投影的最大值，b∈[1,T-10]， 

步骤3.3M_W(b)中的峰即为核小体的二分点Dyad位置，这样便实现了核小体的预测。 

实例 

步骤1提取人类基因组20号染色体的DNA序列(NCBI，www.ncbi.nlm.nih.gov/)，表示为Seq，长度为T=62435964bp，以一个宽度为10bp，步长为1bp的窗口选取DNA序列，第t个窗口的序列为S_t。 

步骤2计算S_t的弯曲度。 

$>C=v^0{(n_2-n_1)}^{-1}\underset{j=n_1}{\overset{n_2}{\Sigma }}({\rho }_j-{\mathrm{i\tau }}_j)\exp \left(\frac{2\mathrm{\pi ij}}{v^0}\right)$>

其中，C的模代表弯曲度，v⁰=10.4bp，是DNA双螺旋一周螺旋周期的碱基对平均数目，n1和n2代表加和范围的起始和终止位置，i为复数单位，ρ和τ是16种二联核苷酸在空间弯曲的幅度，ρ为在平面旋转的幅度，τ为扭曲的幅度，其值见表1。 

步骤3重复1-2，直至窗口到序列Seq的末尾，便形成了Seq的弯曲度信号Signal，长度为T-10。卷积Signal和核小体周期模版信号P，S_covn为卷积的结果，长度为147+T-10-1，从卷积结果S_covn中部取出长度为T-10的信号，标记为S_covn_keep，长度为T-10。 

步骤4检测卷积信号S_covn_keep的峰位置，这样便实现了核小体的预测，图3显示了人类基因组20号染色体8000bp到28000bp的区域的预测的核小体。 

步骤5与实验检测的核小体做比较(D.E.Schones,K.R.Cui,S.Cuddapah,T.Y.Roh,A.Barski,Z.B.Wang,et al.,Dynamic Regulation of Nucleosome Positioning inthe Human Genome,Cell 132(2008)887-898.)。该实验数据是人类CD4+T细胞在激活(TCR信号)和休眠状态下测定的。同时，将本发明方法的预测也与Kaplan等的预测做了比较(Kaplan N,Moore IK,Mittendorf YF,Gossett AJ,Tillo D,Field Y，LeProust EM,Hughes TR,Lieb JD,Widom J,Segal E:The DNA-encodednucleosome organization of a eukaryotic genome.Nature 2009,458:362-366.)。 

从图4可见，在40bp的偏差范围内，本发明的方法可以预测超过53%的实验检测的核小体，而Kaplan等的预测只可以预测49%的核小体。而且，偏差小于40bp时，相比于Kaplan等的预测，本发明方的方法总是能准确地预测更多的核小体。这些结果说明本发明的结果优于Kaplan等的预测(图4)。