基于希尔伯特-黄变换的核小体定位特征提取

刘 珂，丰继华，黄月月，郭亚茹，牟 锦

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

遗传信息的传递和获取是通过核小体的定位来实现的，仅相隔几个bp就可能对基因表达产生重大的影响。因此，染色体的基本单元——核小体[1]，是影响表观遗传状态的主要因素。核小体的定位对基因的表达调控有重要的影响，它的定位变化总是伴随着基因从抑制到转录状态的转变。大量实验结果表明，核小体的形成和在染色质的精准定位是真核生物基因表达所必需的。有人提出核小体的形成及其在染色质上的精准定位有以下两方面的作用：(1)提供一个支架结构，使转录因子之间的信息传递更有效；(2)染色质结构的不均一性，即某些区域不形成核小体，保证了转录因子易于接近染色质模板。经研究表明真核生物基因转录起始区域的核小体定位分布具有高度保守性[2-3]，任何真核生物其核小体在基因转录区和编码区的定位图谱上总体都呈现一种周期性振荡衰减趋势[4-5]。但进一步观察后会发现：在细节上，不同物种的核小体占位图谱是有明显差别的，这种差异性可能代表了物种在染色质结构和功能上的进化印迹[6]。但是目前，对于核小体定位预测和定位性质的分析大多数为定性分析，单细胞生物和多细胞生物核小体在结构上的差异性仍然不明确[7-8]，影响了我们对基因组的认识进程。因此，如何定量获取核小体定位分布的空间尺度特征和频谱特征尤为迫切。

利用传统的频谱分析方法——快速傅里叶变换(FFT)获得了分布模型的功率谱和相位谱，并在此基础上对酵母和果蝇两种模式生物的频谱特征进行了比较和分析[9]。但傅里叶变换存在的问题是：(1)没有时间-频率的定位功能；(2)不能获得瞬时频率；(3)傅里叶变换在时间和频率分辨率上具有局限性。尽管通过快速傅里叶变换得到了初步的核小体定位信号的频谱特征，但想要进一步研究核小体组织结构，必须要有一种具有自适应性的高分辨率频谱分析方法的介入。

希尔伯特-黄变换(HHT)的独特之处在于分析信号过程中不需要基函数，其由两部分组成：经验模式分解(EMD)[10-11]和Hilbert谱分析。首先，采用EMD方法将定位信号序列分解为一组本征模态函数，再通过希尔伯特变换可得到具有明显物理意义的参量：瞬时频率谱和边际谱。其次，本文通过加入类型强度不同的白噪声改善希尔伯特-黄变换存在的模态混叠。通过改进的希尔伯特-黄变换定量的提取核小体定位信号的特征，从而明确单细胞生物和多细胞生物核小体在结构上的差异性，加快我们对基因组的认识进程。

1 希尔伯特-黄变换

1.1 数据来源

研究的主要目的是比较单细胞生物和多细胞生物在核小体组织结构上的差异性，同时又要考虑其可比较性。由于哺乳动物细胞已分化，在选择分化时期的标准上比较困难，所以选择果蝇胚胎期细胞，既具有单细胞的典型特征又具有多细胞分化过渡的特征，满足需要的可比较性。

选取的两个实验数据其实验精度和分辨率均已达到实验要求。

综上考虑，选择来源于William Lee[12]等人于2007年发布的高分辨率酵母核小体定位率实验数据和来自于2008年Travis N. Mavrich[13]等人获得的果蝇胚胎期核小体定位实验数据(见图1)。

图1 两种模式生物核小体定位图谱Fig.1 Nucleosome occupancy map of two model organisms

1.2 希尔伯特-黄变换原理

Wu和Huang提出的希尔伯特-黄变换(HHT)主要分为两步：经验模态分解与希尔伯特变换。经验模式分解(EMD)主要用于将非线性和非平稳信号分解为一系列本征模态函数(IMF)。将通过EMD分解的IMF进行Hilbert变换可以得到具有实际物理意义的瞬时频率。

(1)

(2)

其中：

(3)

其中ai(t)是ci(t)的瞬时幅度，它可以反映ci(t)的能量随时间变化，θi(t)是ci(t)的瞬时相位。很容易获得相位，每个IMF的瞬时频率可以通过相位的导数来定义，如公式(4)所示：

(4)

我们可以用以下形式表示数据x(t)，其不包含残留物rn(t)

(5)

H(ω,t)定义为希尔伯特振幅谱：

(6)

通过定义希尔伯特振幅谱，经验模态分解得到的本征模态函数进行希尔伯特变换即可获得时频分析中两个重要的参量：瞬时频率与瞬时幅度谱。黄锷博士在此基础上引进两个新的物理参量：边际谱h(ω)与瞬时能量密度级IE：

(7)

(8)

1.3 希尔伯特-黄变换的改进

图2 加入四种白噪声后的相关系数曲线(酵母)Fig.2 Correlation coefficients of four white noise distributions(Yeast)

由图2可知，酵母核小体定位谱在加入幅值系数为0.08的高斯分布白噪声后，分解得到的IMF整体与原始信号相关性最大，分解效果最好；加入指数分布白噪声后，分解分解得到的IMF整体与原始信号相关性最小，分解效果最差。

对比图2与图3，果蝇核小体定位分布谱同样是在添加指数分布白噪声时分解效果最差，但当加入幅值系数为0.06的瑞丽分布白噪声时，IMF整体与原始信号的相关性最高，效果最好。

图3 四种白噪声分布下相关性对比(果蝇)Fig.3 Correlation comparison of four white noise distributions (Drosophila)

为了便于描述，将添加特定白噪声后的经验模态分解称为NEEMD。分别以酵母和果蝇核小体定位分布谱作为输入信号，比较EMD、EEMD与NEEMD分解后IMF的方差与贡献率，以此作为模态混叠现象消除程度的评判标准。

图4(a)为EMD分解法的处理结果，通过EMD分解将信号分为4层，图中r4为分解4层之后余留下来的残余信号。从图中我们可以看到IMF1到IMF5都包含原始信号，我们无法区分哪一层具有原始信号信息最多。表1与表2显示IMF3的方差与贡献率是最大的，但是从图像上看IMF1具有原始信号信息最多，因此唯一的解释是使用EMD方法出现了模式混合现象。

图4(b)为EEMD分解法的处理结果，添加的白噪声遵循很多参考文献中的建议值：添加白噪声次数为100，白噪声遵循高斯分布，幅值系数ε为0.2。由图像可以看出，IMF6和IMF7波形具有很强的相似性。

图4(c)为NEEMD分解法处理的结果，根据核小体定位分布特点调整参数，选取添加高斯分布白噪声，噪声的次数为100，幅值系数ε为0.08。从实验结果看，该方法可以极大改善模态混叠现象。表1与表2证实了这一点：EEMD方法显示IMF1方差最大、贡献率最高，应该具有原始信号信息最多， NEEMD方法显示IMF1同样为方差最大、贡献率最高的一项，并且NEEMD 中IMF1的贡献率比在EEMD更高。所以NEEMD相对于EEMD来说分解效果更好。

与酵母分解结果相对应，图5(a)为果蝇核小体定位分布谱的EMD分解结果， 图5(b)为EEMD分解结果。在EEMD实验中，添加的白噪声参考了其它文献的建议值：添加白噪声次数为100，白噪声遵循高斯分布，幅值系数ε为0.2。图5(c)为NEEMD分解结果。我们在实验中根据核小体定位分布特点调整了参数，选取添加噪声的次数为100，白噪声的幅值系数ε为0.08，噪声类型为瑞丽分布。尽管从分解图谱中不能直观看出改善程度，但是从表3与表4中可以得出以下结论：EMD方法显示IMF3贡献率最高，而图像直观显示IMF1与原始信号最为接近，所以具有模态混叠现象；EEMD方法显示IMF1方差最大、贡献率最高，在NEEMD方法中 IMF1同样为方差最大、贡献率最高的一项，并且在NEEMD中IMF1贡献率比EEMD方法产生的IMF1更高，因此在应用于果蝇核小体定位分布时，NEEMD相较于EEMD就具有更好的分解效果。

图4 核小体定位分布分解谱(酵母)Fig 4 Decomposition spectrum of nucleosome positioning distribution (Yeast)

以上实验表明：不同模型中加入白噪声的类型和幅度与所分解的数据相关，只有加入合适的白噪声，才能最大程度上改善模态混叠现象，提高分解的准确性。

图5 核小体定位分布分解谱(果蝇)Fig.5 Decomposition spectrum of nucleosome positioning distribution (Drosophila)

表1 三种分解方法中IMF的方差(酵母)Table 1 IMFs’ variance of three decomposition methods (Yeast)

表2 三种分解方法中IMF的贡献率(酵母)Table 2 IMFs’ normalized variance of three decomposition methods (Yeast)

表3 三种分解方法中IMF的方差(果蝇)Table 3 IMFs’ variance of three decomposition methods (Drosophila)

表4 三种分解方法中IMF的贡献率(果蝇)Table 4 IMFs’ normalized variance of three decomposition methods (Drosophila)

2 核小体定位特征提取

希尔伯特-黄变换(HHT)与傅里叶变换的相似之处在于，二者都是将时域中的信号进行解析后，便于从频域进行分析。在此，我们用希尔伯特-黄变换对核小体定位分布进行分析的目的是有二个：一是弥补傅里叶变换在频谱分析中的不足，二是通过HHT更深入研究两个物种核小体分布的频谱差异性。

2.1 酵母与果蝇核小体定位信号的时域分解特征

图6(a)、(b)分别表示酵母与果蝇核小体定位信号的时域本征模态分解三维谱。从图中可以看出二个物种核小体定位信号在空间尺度上的细微差异，与果蝇相比较，在相同IMF分层下，酵母信号总体要平坦一些，信号曲面在复杂程度上要比果蝇低。以第4层分解结果为例，在同一层中，果蝇核小体定位信号的突起比酵母多，说明多细胞生物核小体组织的复杂程度要高于单细胞生物。

图6 两种模式生物三维采样点-IMF层数-幅值谱Fig.6 The 3D sample point-IMF layer number-amplitude spectrums of two model organisms

2.2 酵母与果蝇核小体定位信号的频域分解特征

图7(a)、(b)分别表示酵母与果蝇核小体定位信号的频域本征模态分解三维谱。从分层信号频谱分布可以看到，果蝇核小体定位信号成份中的高频部分要多于酵母。同样是在第五层信号中，酵母还保留了多个频率分量，而果蝇则少得多。这与时域分解特征得到的结论相一致，即果蝇的核小体排列相对酵母有更多的变化，也更为复杂。

图7 两种模式生物三维采样点-IMF层数-频率谱Fig.7 The 3D sample point-IMF layer number-frequency spectrums of two model organisms

2.3 酵母与果蝇核小体定位信号的二维幅频曲线

图8(a)、(b)分别为酵母与果蝇核小体定位信号的二维幅频曲线。从图中每个IMF分量及剩余分量的幅频曲线可以看出两者在大体趋势上是相近的，在同等情况下都是分为9个IMF分量，且每个分量的幅频曲线趋势大致相同，这从某种程度上反应了生物在核小体组织上进化的保守性。

图8 两种模式生物IMF幅频曲线Fig.8 IMF amplitude frequency curve of two model organisms

2.4 酵母与果蝇核小体定位信号的边际谱

边际谱的定义是对Hilbert谱(二元函数)进行时间积分，积分结果是自变量w的一元函数，即幅值谱，它也描述信号的幅值在整个频率段上随频率的变化情况。边际谱的意义是信号中瞬时频率的总幅值的大小，将所有时刻某一频率的幅值加起来就是该频率的总幅值，即边际谱线的高度。

图9(a)、(b)分别为酵母和果蝇核小体定位信号的边际谱，其走势进一步验证了两种模式生物进化过程具有的保守性，但是通过进一步研究发现虽然两者走势相同但是在0.002～0.003 Hz位置酵母只出现了一次最高峰，而果蝇在此处具有两个几乎相同的波峰，并且在0.005 Hz附近酵母的核小体定位模型只有一个小幅度的波峰，而果蝇的核小体定位模型出现了相对而言比较大幅度的上升，并且两个波峰也是几乎相同的。在信号幅值上酵母的骤减程度明显大于果蝇，相较于酵母而言，果蝇核小体定位模型信号幅值变化较为平缓。

在某种程度上而言，边际谱的差异反映了生物从单细胞进化到多细胞过程中核小体组织形式发生的微妙变化。

图9 两种模式生物边际谱Fig.9 Marginal spectrum of two model organisms

2.5 酵母与果蝇核小体定位信号的瞬时频率

图10(a)、(b)分别为酵母和果蝇核小体定位信号的瞬时频率。从图中看是大致相似的，在中心位置即转录起始位点TSS附近略有不同，酵母的核小体定位模型信号在转录起始位点急剧减小，而果蝇在TSS附近相对来说频率幅值缩降程度较小，或许体现了较为高等的生物在转录起始位点的自我保护机制，使得较为高等的真核模式生物在遭遇某种变故时能延缓自身的突变，以达到保护自身的作用。

图10 两种模式生物瞬时频率Fig10 Instantaneous frequency of two model organisms

2.6 酵母与果蝇核小体定位信号的三维HHT时频谱

图11(a)、(b)分别为酵母和果蝇核小体定位信号的HHT三维时频谱。其中X轴为归一化后的频率，Y轴为采样点数，Z轴为信号幅值。从图中可直观看出，果蝇的核小体定位分布比酵母更为复杂，但是它们的分布趋势是大致相似的，这从一定程度上再次印证了生物在核小体组织上进化的保守性。

图11 两种模式生物HHT时频谱三维显示Fig.11 The 3D HHT spectrums of two model organisms

3 结 语

在生物进化过程中，由单细胞到多细胞核小体组织的演变既存在变异性也存在某种保守性，使得生物可以在保留优良特性的基础上加入某种抗干扰能力，适应生存需要。本文通过希尔伯特-黄变换(HHT)方法证实了生物进化结构与其复杂程度是一致的，通过HHT定量提取单细胞生物和多细胞生物核小体定位信号特征，进一步验证两物种之间在核小体组织结构上显著的差异性。