基于希爾伯特-黃變換的核小體定位特征提取

劉 珂，豐繼華，黃月月，郭亞茹，牟 錦

(云南民族大學 電氣信息工程學院, 昆明 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定量提取單細胞生物和多細胞生物核小體定位信號特征，進一步驗證兩物種之間在核小體組織結構上顯著的差異性。