真核模式生物核小体分布参数模型分析

丰继华, 郭亚茹, 牟 锦, 黄月月， 刘 珂， 范力栋

(云南民族大学 电气信息工程学院, 昆明 650504)

在生物遗传和进化过程中，除了以DNA编码为代表的“硬”遗传物质外，还存在着另一种重要的“软”遗传物质，即基因组中的功能性蛋白质及其化学修饰。以染色体为例，其在微观上组织得井然有序，使得DNA序列在空间结构上非常紧凑，但又互不干扰，尤其是一些重要的编码片段会受到特殊“保护”。基因组实现这些功能主要依赖于一种染色质基本亚基——核小体(Nucleosome)[1]。

核小体在基因组中担负着两个重要作用[2]：1)压缩和折叠DNA，以适应细胞核空间的大小；2)限制DNA的易接近性。细胞广泛利用后一种功能来调控代谢过程，其中一个关键步骤是通过核小体定位对基因表达进行精细的调节[3-4]。核小体沿DNA的定位控制了遗传信息传递和获取，仅相隔几个核苷酸的替代位置就可能对基因表达产生深远的影响。因此，核小体作为染色体的基本单元，被认为是影响表观遗传状态的主要驱动因素。然而，单细胞生物和多细胞生物核小体在结构上的差异性仍然不明确[5-6]，影响了对基因组的认识进程。因此，如何定量获取核小体空间尺度特征尤为迫切。

研究发现，真核生物基因转录起始区域的核小体分布具有高度保守性[7-8]。无论是单细胞的酵母，还是属于高级哺乳动物的人类，其核小体在基因转录区和编码区的定位图谱总体上都呈现出一种周期性振荡衰减趋势[9-10]。但深入观察后会发现：不同物种的核小体占位图谱在细节上是有明显差别的，这种差异性可能代表了物种在染色质结构和功能上的进化印迹[11]。在此提出构建一个高精度数学模型对酵母和果蝇核小体分布进行空间尺度和频域分析。

1 核小体定位模型构建

1.1 酵母与果蝇核小体定位图谱

1.1.1 数据来源

本文使用的生物实验数据主要来源于两个研究团队。酵母数据来源于Lee等[8]于2007发布的高分辨率酵母核小体定位率实验数据；果蝇数据则来自于Mavrich等[12]获得的果蝇胚胎期核小体定位实验数据。

由于真核生物在基因转录起始区域周围有着共同的性质，本文选取的两个实验数据其实验精度已达到我们的实验要求。

1.1.2 转录起始位点周围核小体定位图谱

由于基因转录起始位点(TSS)周围核小体占位具有高度的保守性和代表性，我们选取了酵母和果蝇基因转录起始位点上、下游1000 bp范围，对齐平均后的核小体定位图谱(如图1)，并基于定位图谱分别建立拟合模型。

图1 中，纵坐标是以两个物种基因转录起始位点(TSS)为中心，对齐平均并归一化的核小体占位率曲线；横坐标是DNA的长度计量单位bp(核苷酸碱基对，Base Pair)。

1.2 拟合函数选取

根据实验核小体的定位特征，我们分别使用了多项式、傅里叶级数、高斯函数和正弦函数构建定位模型。表1为以上4种函数对酵母和果蝇核小体定位曲线的拟合性能指标，拟合结果如图2所示。

图1 两种模式生物核小体定位图谱

物种拟合函数和方差拟合优度标准差自由度校正决定系数酵母多项式34.80370.58660.13221991 0.5847傅里叶6.37800.92420.05671983 0.9236高斯7.24460.91390.0605 1977 0.9129正弦函数1.81380.99170.01881974 0.9916 果蝇多项式27.19420.70060.11691991 0.6992傅里叶6.61830.92710.05781983 0.9265 高斯2.10640.97680.03261977 0.9765正弦函数1.08940.98800.02351974 0.9878

综合比较4种函数的拟合性能和拟合效果后，发现由正弦复合函数构成的模型拟合精度最高，如图2-g、h所示，模型几乎能在整个区域跟踪核小体的定位特征；高斯函数次之；如图2-f所示只是在远离TSS区域出现了误差；而图2-e中除在远离TSS区域出现了误差外，TSS附近区域也存在较大误差。傅里叶函数能基本拟合占位图谱趋势，但误差较大(图2-c、d)，其精度不能满足后续分析需要；由多项式构建的模型则拟合度最低，无法捕获定位特征的细节特征(图2-a、b)。根据以上结论，最终确定以正弦函数构建核小体定位模型，并基于该模型进行后续的尺度和频域分析。

图2 4种函数的拟合曲线

1.3 定位模型及参数

本文构建的正弦复合函数如下：

(1)

(1)式由9个正弦函数线性组合而成，其中An为正弦函数的振幅，ωn为正弦函数的角速度，φn为正弦函数的初始相位。经过对酵母和果蝇核小体占位率曲线拟合后，得到模型参数如表2所示。

表2 核小体定位模型参数

2 定位模型尺度计算

利用经典信号处理方法挖掘基因组数据的内在机理，在生物信息学领域得到广泛应用[13]。为了研究核小体定位的动态与静态特性，在此分别对两种模式生物的定位模型(式1)进行微积分计算。

2.1 模型微分特性

对式(1)求导后，得：

(2)

f′(x)表示核小体分布的动态特性(占位变化)。令定位模型的导函数f′(x)=0，可得到核小体定位模型f(x)在TSS周围的各极值点(表3)。其中，极大值表示单个核小体的分布中心，极小值表示连接DNA(核小体之间的DNA)分布中心(图3)。

为便于分析，我们根据核小体分布中心位置与TSS距离的远近，对其依次进行了编号(图3和表3)。

图3 两种模式生物定位模型求导函数曲线

表3 核小体与连接DNA的分布中心位置(以TSS为坐标原点，单位为bp)

从表3可以观察到，尽管酵母和果蝇在TSS附近均形成一个最明显的核小体缺失区域，但区别在于，果蝇核小体缺失区域范围更大，宽度近400 bp，其特征更复杂，且存在两个较弱的分布峰。与之对应，酵母核小体缺失区域宽度仅300 bp左右，特征较为单一。

2.2 模型积分特性

为了确定两个物种在TSS附近单核小体的定位稳定性，我们以连接DNA分布中心作为积分区间，对式(1)计算定积分。

(3)

式(3)中，积分上、下限a、b取值为相邻连接DNA分布中心坐标，每个积分区间都包含有一个核小体分布中心。在此，我们计算了TSS上、下游共10个单核小体分布区间积分值(图4)，An、ωn、φn取值为模型参数(表2)。

图4是归一化后的单核小体分布区间定积分值，一定程度上代表了相应位置核小体的稳定性(定位强度)。从计算结果可看到，酵母和果蝇核小体在TSS下游的定位强度总体要高于上游区域。但有趣的是，在+1和+2核小体分布区域，酵母的+2核小体定位强度要略高于+1核小体，而果蝇在总体上呈现出从+1到+5核小体定位强度依次递减的规律。这一进化上的差异，再次证明+1核小体在多细胞生物中承担了更为关键的作用[14]。

图4 转录起始位点周围核小体的定位稳定性

3 定位模型频谱分析

3.1 模型傅里叶变换

根据定位模型式(1)，我们对其进行了周期信号傅里叶变换，用于分析核小体定位模型的频域特征。

(4)

3.2 模型的功率谱与相位谱

在式(4)基础上，可得到模型的功率谱表达式。

(5)

图5 两种模式生物功率谱和相位谱

图5横坐标为归一化频率，其中图5-a观察到两种模式生物的核小体定位频率(体现了占位变化速率)主要集中于0～0.007 Hz低频段，表明核小体在基因转录区的位置较稳定，进一步证实了这一区域核小体组织的保守性和稳健性，这一性质对于细胞的生存具有重要意义[15]。观察两个物种频谱后，发现单细胞酵母核小体占位变化只集中在频点0～0.002 Hz附近，而多细胞生物果蝇的核小体占位变化频点则在接近0.005～0.006 Hz附近出现一个明显峰值，这一特征可能代表了二者在进化上的差异，多细胞生物的核小体定位机制更为复杂。

图5-b、c为两种模式生物核小体定位模型的相位谱。总体上，相位谱随着频率升高呈现出周期性振荡趋势。在0～1.25×10-11Hz区间，酵母核小体的相位呈现出近似线性增长趋势，而果蝇则呈现出复杂的非线性变化。两个物种的相位谱在1.25×10-11Hz附近同时出现了一个拐点，但随着频率的增加，酵母核小体定位的相位向负方向急剧变化后，呈现出缓慢振荡上行的趋势；与之相反，果蝇核小体定位的相位却在向正方向急剧变化后，呈现出了缓慢振荡下行的趋势。这一现象与图3中的定位模型曲线是一致的，即以+1核小体位置为参照，在转录起始位点下游区域，果蝇核小体出现的位置要明显滞后于酵母核小体出现位置，其中的生物意义还有待进一步研究。

4 结论

本文研究了两种真核模式生物核小体定位模型的构建问题。基于定位模型的空间尺度和频谱特征，比较了两种生物在染色体结构上的差异性，探索了利用数学模型对核小体定位过程进行定量分析的方法，在此基础上得出以下结论：1)利用数学模型对酵母和果蝇核小体定位曲线进行拟合，为应用经典信号处理方法分析核小体定位特征提供了一种新思路；2)正弦作为一种单一频率的基本信号与作为染色质基本单元的核小体定位特征高度拟合，其背后可能隐藏着未知的自然规律;3)建立的核小体定位模型在一定程度上解释了两个物种在进化过程中存在的差异性和共性(物种间的保守性)，为后续研究其他生物核小体定位及染色质进化问题做了铺垫。

当然，本文使用核小体定位模型和时频分析方法所获得的结论，只是针对核小体定位的理论分析，其实际机理非常复杂，其间还有诸多未解决的问题需要进一步探索。