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

(云南民族大學(xué) 電氣信息工程學(xué)院, 昆明 650504)

遺傳信息的傳遞和獲取是通過(guò)核小體的定位來(lái)實(shí)現(xiàn)的，僅相隔幾個(gè)bp就可能對(duì)基因表達(dá)產(chǎn)生重大的影響。因此，染色體的基本單元——核小體[1]，是影響表觀遺傳狀態(tài)的主要因素。核小體的定位對(duì)基因的表達(dá)調(diào)控有重要的影響，它的定位變化總是伴隨著基因從抑制到轉(zhuǎn)錄狀態(tài)的轉(zhuǎn)變。大量實(shí)驗(yàn)結(jié)果表明，核小體的形成和在染色質(zhì)的精準(zhǔn)定位是真核生物基因表達(dá)所必需的。有人提出核小體的形成及其在染色質(zhì)上的精準(zhǔn)定位有以下兩方面的作用：(1)提供一個(gè)支架結(jié)構(gòu)，使轉(zhuǎn)錄因子之間的信息傳遞更有效；(2)染色質(zhì)結(jié)構(gòu)的不均一性，即某些區(qū)域不形成核小體，保證了轉(zhuǎn)錄因子易于接近染色質(zhì)模板。經(jīng)研究表明真核生物基因轉(zhuǎn)錄起始區(qū)域的核小體定位分布具有高度保守性[2-3]，任何真核生物其核小體在基因轉(zhuǎn)錄區(qū)和編碼區(qū)的定位圖譜上總體都呈現(xiàn)一種周期性振蕩衰減趨勢(shì)[4-5]。但進(jìn)一步觀察后會(huì)發(fā)現(xiàn)：在細(xì)節(jié)上，不同物種的核小體占位圖譜是有明顯差別的，這種差異性可能代表了物種在染色質(zhì)結(jié)構(gòu)和功能上的進(jìn)化印跡[6]。但是目前，對(duì)于核小體定位預(yù)測(cè)和定位性質(zhì)的分析大多數(shù)為定性分析，單細(xì)胞生物和多細(xì)胞生物核小體在結(jié)構(gòu)上的差異性仍然不明確[7-8]，影響了我們對(duì)基因組的認(rèn)識(shí)進(jìn)程。因此，如何定量獲取核小體定位分布的空間尺度特征和頻譜特征尤為迫切。

利用傳統(tǒng)的頻譜分析方法——快速傅里葉變換(FFT)獲得了分布模型的功率譜和相位譜，并在此基礎(chǔ)上對(duì)酵母和果蠅兩種模式生物的頻譜特征進(jìn)行了比較和分析[9]。但傅里葉變換存在的問(wèn)題是：(1)沒(méi)有時(shí)間-頻率的定位功能；(2)不能獲得瞬時(shí)頻率；(3)傅里葉變換在時(shí)間和頻率分辨率上具有局限性。盡管通過(guò)快速傅里葉變換得到了初步的核小體定位信號(hào)的頻譜特征，但想要進(jìn)一步研究核小體組織結(jié)構(gòu)，必須要有一種具有自適應(yīng)性的高分辨率頻譜分析方法的介入。

希爾伯特-黃變換(HHT)的獨(dú)特之處在于分析信號(hào)過(guò)程中不需要基函數(shù)，其由兩部分組成：經(jīng)驗(yàn)?zāi)Ｊ椒纸?EMD)[10-11]和Hilbert譜分析。首先，采用EMD方法將定位信號(hào)序列分解為一組本征模態(tài)函數(shù)，再通過(guò)希爾伯特變換可得到具有明顯物理意義的參量：瞬時(shí)頻率譜和邊際譜。其次，本文通過(guò)加入類(lèi)型強(qiáng)度不同的白噪聲改善希爾伯特-黃變換存在的模態(tài)混疊。通過(guò)改進(jìn)的希爾伯特-黃變換定量的提取核小體定位信號(hào)的特征，從而明確單細(xì)胞生物和多細(xì)胞生物核小體在結(jié)構(gòu)上的差異性，加快我們對(duì)基因組的認(rèn)識(shí)進(jìn)程。

1 希爾伯特-黃變換

1.1 數(shù)據(jù)來(lái)源

研究的主要目的是比較單細(xì)胞生物和多細(xì)胞生物在核小體組織結(jié)構(gòu)上的差異性，同時(shí)又要考慮其可比較性。由于哺乳動(dòng)物細(xì)胞已分化，在選擇分化時(shí)期的標(biāo)準(zhǔn)上比較困難，所以選擇果蠅胚胎期細(xì)胞，既具有單細(xì)胞的典型特征又具有多細(xì)胞分化過(guò)渡的特征，滿足需要的可比較性。

選取的兩個(gè)實(shí)驗(yàn)數(shù)據(jù)其實(shí)驗(yàn)精度和分辨率均已達(dá)到實(shí)驗(yàn)要求。

綜上考慮，選擇來(lái)源于William Lee[12]等人于2007年發(fā)布的高分辨率酵母核小體定位率實(shí)驗(yàn)數(shù)據(jù)和來(lái)自于2008年Travis N. Mavrich[13]等人獲得的果蠅胚胎期核小體定位實(shí)驗(yàn)數(shù)據(jù)(見(jiàn)圖1)。

圖1 兩種模式生物核小體定位圖譜Fig.1 Nucleosome occupancy map of two model organisms

1.2 希爾伯特-黃變換原理

Wu和Huang提出的希爾伯特-黃變換(HHT)主要分為兩步：經(jīng)驗(yàn)?zāi)B(tài)分解與希爾伯特變換。經(jīng)驗(yàn)?zāi)Ｊ椒纸?EMD)主要用于將非線性和非平穩(wěn)信號(hào)分解為一系列本征模態(tài)函數(shù)(IMF)。將通過(guò)EMD分解的IMF進(jìn)行Hilbert變換可以得到具有實(shí)際物理意義的瞬時(shí)頻率。

(1)

(2)

其中：

(3)

其中ai(t)是ci(t)的瞬時(shí)幅度，它可以反映ci(t)的能量隨時(shí)間變化，θi(t)是ci(t)的瞬時(shí)相位。很容易獲得相位，每個(gè)IMF的瞬時(shí)頻率可以通過(guò)相位的導(dǎo)數(shù)來(lái)定義，如公式(4)所示：

(4)

我們可以用以下形式表示數(shù)據(jù)x(t)，其不包含殘留物rn(t)

(5)

H(ω,t)定義為希爾伯特振幅譜：

(6)

通過(guò)定義希爾伯特振幅譜，經(jīng)驗(yàn)?zāi)B(tài)分解得到的本征模態(tài)函數(shù)進(jìn)行希爾伯特變換即可獲得時(shí)頻分析中兩個(gè)重要的參量：瞬時(shí)頻率與瞬時(shí)幅度譜。黃鍔博士在此基礎(chǔ)上引進(jìn)兩個(gè)新的物理參量：邊際譜h(ω)與瞬時(shí)能量密度級(jí)IE：

(7)

(8)

1.3 希爾伯特-黃變換的改進(jìn)

圖2 加入四種白噪聲后的相關(guān)系數(shù)曲線(酵母)Fig.2 Correlation coefficients of four white noise distributions(Yeast)

由圖2可知，酵母核小體定位譜在加入幅值系數(shù)為0.08的高斯分布白噪聲后，分解得到的IMF整體與原始信號(hào)相關(guān)性最大，分解效果最好；加入指數(shù)分布白噪聲后，分解分解得到的IMF整體與原始信號(hào)相關(guān)性最小，分解效果最差。

對(duì)比圖2與圖3，果蠅核小體定位分布譜同樣是在添加指數(shù)分布白噪聲時(shí)分解效果最差，但當(dāng)加入幅值系數(shù)為0.06的瑞麗分布白噪聲時(shí)，IMF整體與原始信號(hào)的相關(guān)性最高，效果最好。

圖3 四種白噪聲分布下相關(guān)性對(duì)比(果蠅)Fig.3 Correlation comparison of four white noise distributions (Drosophila)

為了便于描述，將添加特定白噪聲后的經(jīng)驗(yàn)?zāi)B(tài)分解稱(chēng)為NEEMD。分別以酵母和果蠅核小體定位分布譜作為輸入信號(hào)，比較EMD、EEMD與NEEMD分解后IMF的方差與貢獻(xiàn)率，以此作為模態(tài)混疊現(xiàn)象消除程度的評(píng)判標(biāo)準(zhǔn)。

圖4(a)為EMD分解法的處理結(jié)果，通過(guò)EMD分解將信號(hào)分為4層，圖中r4為分解4層之后余留下來(lái)的殘余信號(hào)。從圖中我們可以看到IMF1到IMF5都包含原始信號(hào)，我們無(wú)法區(qū)分哪一層具有原始信號(hào)信息最多。表1與表2顯示IMF3的方差與貢獻(xiàn)率是最大的，但是從圖像上看IMF1具有原始信號(hào)信息最多，因此唯一的解釋是使用EMD方法出現(xiàn)了模式混合現(xiàn)象。

圖4(b)為EEMD分解法的處理結(jié)果，添加的白噪聲遵循很多參考文獻(xiàn)中的建議值：添加白噪聲次數(shù)為100，白噪聲遵循高斯分布，幅值系數(shù)ε為0.2。由圖像可以看出，IMF6和IMF7波形具有很強(qiáng)的相似性。

圖4(c)為NEEMD分解法處理的結(jié)果，根據(jù)核小體定位分布特點(diǎn)調(diào)整參數(shù)，選取添加高斯分布白噪聲，噪聲的次數(shù)為100，幅值系數(shù)ε為0.08。從實(shí)驗(yàn)結(jié)果看，該方法可以極大改善模態(tài)混疊現(xiàn)象。表1與表2證實(shí)了這一點(diǎn)：EEMD方法顯示IMF1方差最大、貢獻(xiàn)率最高，應(yīng)該具有原始信號(hào)信息最多， NEEMD方法顯示IMF1同樣為方差最大、貢獻(xiàn)率最高的一項(xiàng)，并且NEEMD 中IMF1的貢獻(xiàn)率比在EEMD更高。所以NEEMD相對(duì)于EEMD來(lái)說(shuō)分解效果更好。

與酵母分解結(jié)果相對(duì)應(yīng)，圖5(a)為果蠅核小體定位分布譜的EMD分解結(jié)果， 圖5(b)為EEMD分解結(jié)果。在EEMD實(shí)驗(yàn)中，添加的白噪聲參考了其它文獻(xiàn)的建議值：添加白噪聲次數(shù)為100，白噪聲遵循高斯分布，幅值系數(shù)ε為0.2。圖5(c)為NEEMD分解結(jié)果。我們?cè)趯?shí)驗(yàn)中根據(jù)核小體定位分布特點(diǎn)調(diào)整了參數(shù)，選取添加噪聲的次數(shù)為100，白噪聲的幅值系數(shù)ε為0.08，噪聲類(lèi)型為瑞麗分布。盡管從分解圖譜中不能直觀看出改善程度，但是從表3與表4中可以得出以下結(jié)論：EMD方法顯示IMF3貢獻(xiàn)率最高，而圖像直觀顯示IMF1與原始信號(hào)最為接近，所以具有模態(tài)混疊現(xiàn)象；EEMD方法顯示IMF1方差最大、貢獻(xiàn)率最高，在NEEMD方法中 IMF1同樣為方差最大、貢獻(xiàn)率最高的一項(xiàng)，并且在NEEMD中IMF1貢獻(xiàn)率比EEMD方法產(chǎn)生的IMF1更高，因此在應(yīng)用于果蠅核小體定位分布時(shí)，NEEMD相較于EEMD就具有更好的分解效果。

圖4 核小體定位分布分解譜(酵母)Fig 4 Decomposition spectrum of nucleosome positioning distribution (Yeast)

以上實(shí)驗(yàn)表明：不同模型中加入白噪聲的類(lèi)型和幅度與所分解的數(shù)據(jù)相關(guān)，只有加入合適的白噪聲，才能最大程度上改善模態(tài)混疊現(xiàn)象，提高分解的準(zhǔn)確性。

圖5 核小體定位分布分解譜(果蠅)Fig.5 Decomposition spectrum of nucleosome positioning distribution (Drosophila)

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

表2 三種分解方法中IMF的貢獻(xiàn)率(酵母)Table 2 IMFs’ normalized variance of three decomposition methods (Yeast)

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

表4 三種分解方法中IMF的貢獻(xiàn)率(果蠅)Table 4 IMFs’ normalized variance of three decomposition methods (Drosophila)

2 核小體定位特征提取

希爾伯特-黃變換(HHT)與傅里葉變換的相似之處在于，二者都是將時(shí)域中的信號(hào)進(jìn)行解析后，便于從頻域進(jìn)行分析。在此，我們用希爾伯特-黃變換對(duì)核小體定位分布進(jìn)行分析的目的是有二個(gè)：一是彌補(bǔ)傅里葉變換在頻譜分析中的不足，二是通過(guò)HHT更深入研究?jī)蓚€(gè)物種核小體分布的頻譜差異性。

2.1 酵母與果蠅核小體定位信號(hào)的時(shí)域分解特征

圖6(a)、(b)分別表示酵母與果蠅核小體定位信號(hào)的時(shí)域本征模態(tài)分解三維譜。從圖中可以看出二個(gè)物種核小體定位信號(hào)在空間尺度上的細(xì)微差異，與果蠅相比較，在相同IMF分層下，酵母信號(hào)總體要平坦一些，信號(hào)曲面在復(fù)雜程度上要比果蠅低。以第4層分解結(jié)果為例，在同一層中，果蠅核小體定位信號(hào)的突起比酵母多，說(shuō)明多細(xì)胞生物核小體組織的復(fù)雜程度要高于單細(xì)胞生物。

圖6 兩種模式生物三維采樣點(diǎn)-IMF層數(shù)-幅值譜Fig.6 The 3D sample point-IMF layer number-amplitude spectrums of two model organisms

2.2 酵母與果蠅核小體定位信號(hào)的頻域分解特征

圖7(a)、(b)分別表示酵母與果蠅核小體定位信號(hào)的頻域本征模態(tài)分解三維譜。從分層信號(hào)頻譜分布可以看到，果蠅核小體定位信號(hào)成份中的高頻部分要多于酵母。同樣是在第五層信號(hào)中，酵母還保留了多個(gè)頻率分量，而果蠅則少得多。這與時(shí)域分解特征得到的結(jié)論相一致，即果蠅的核小體排列相對(duì)酵母有更多的變化，也更為復(fù)雜。

圖7 兩種模式生物三維采樣點(diǎn)-IMF層數(shù)-頻率譜Fig.7 The 3D sample point-IMF layer number-frequency spectrums of two model organisms

2.3 酵母與果蠅核小體定位信號(hào)的二維幅頻曲線

圖8(a)、(b)分別為酵母與果蠅核小體定位信號(hào)的二維幅頻曲線。從圖中每個(gè)IMF分量及剩余分量的幅頻曲線可以看出兩者在大體趨勢(shì)上是相近的，在同等情況下都是分為9個(gè)IMF分量，且每個(gè)分量的幅頻曲線趨勢(shì)大致相同，這從某種程度上反應(yīng)了生物在核小體組織上進(jìn)化的保守性。

圖8 兩種模式生物IMF幅頻曲線Fig.8 IMF amplitude frequency curve of two model organisms

2.4 酵母與果蠅核小體定位信號(hào)的邊際譜

邊際譜的定義是對(duì)Hilbert譜(二元函數(shù))進(jìn)行時(shí)間積分，積分結(jié)果是自變量w的一元函數(shù)，即幅值譜，它也描述信號(hào)的幅值在整個(gè)頻率段上隨頻率的變化情況。邊際譜的意義是信號(hào)中瞬時(shí)頻率的總幅值的大小，將所有時(shí)刻某一頻率的幅值加起來(lái)就是該頻率的總幅值，即邊際譜線的高度。

圖9(a)、(b)分別為酵母和果蠅核小體定位信號(hào)的邊際譜，其走勢(shì)進(jìn)一步驗(yàn)證了兩種模式生物進(jìn)化過(guò)程具有的保守性，但是通過(guò)進(jìn)一步研究發(fā)現(xiàn)雖然兩者走勢(shì)相同但是在0.002～0.003 Hz位置酵母只出現(xiàn)了一次最高峰，而果蠅在此處具有兩個(gè)幾乎相同的波峰，并且在0.005 Hz附近酵母的核小體定位模型只有一個(gè)小幅度的波峰，而果蠅的核小體定位模型出現(xiàn)了相對(duì)而言比較大幅度的上升，并且兩個(gè)波峰也是幾乎相同的。在信號(hào)幅值上酵母的驟減程度明顯大于果蠅，相較于酵母而言，果蠅核小體定位模型信號(hào)幅值變化較為平緩。

在某種程度上而言，邊際譜的差異反映了生物從單細(xì)胞進(jìn)化到多細(xì)胞過(guò)程中核小體組織形式發(fā)生的微妙變化。

圖9 兩種模式生物邊際譜Fig.9 Marginal spectrum of two model organisms

2.5 酵母與果蠅核小體定位信號(hào)的瞬時(shí)頻率

圖10(a)、(b)分別為酵母和果蠅核小體定位信號(hào)的瞬時(shí)頻率。從圖中看是大致相似的，在中心位置即轉(zhuǎn)錄起始位點(diǎn)TSS附近略有不同，酵母的核小體定位模型信號(hào)在轉(zhuǎn)錄起始位點(diǎn)急劇減小，而果蠅在TSS附近相對(duì)來(lái)說(shuō)頻率幅值縮降程度較小，或許體現(xiàn)了較為高等的生物在轉(zhuǎn)錄起始位點(diǎn)的自我保護(hù)機(jī)制，使得較為高等的真核模式生物在遭遇某種變故時(shí)能延緩自身的突變，以達(dá)到保護(hù)自身的作用。

圖10 兩種模式生物瞬時(shí)頻率Fig10 Instantaneous frequency of two model organisms

2.6 酵母與果蠅核小體定位信號(hào)的三維HHT時(shí)頻譜

圖11(a)、(b)分別為酵母和果蠅核小體定位信號(hào)的HHT三維時(shí)頻譜。其中X軸為歸一化后的頻率，Y軸為采樣點(diǎn)數(shù)，Z軸為信號(hào)幅值。從圖中可直觀看出，果蠅的核小體定位分布比酵母更為復(fù)雜，但是它們的分布趨勢(shì)是大致相似的，這從一定程度上再次印證了生物在核小體組織上進(jìn)化的保守性。

圖11 兩種模式生物HHT時(shí)頻譜三維顯示Fig.11 The 3D HHT spectrums of two model organisms

3 結(jié) 語(yǔ)

在生物進(jìn)化過(guò)程中，由單細(xì)胞到多細(xì)胞核小體組織的演變既存在變異性也存在某種保守性，使得生物可以在保留優(yōu)良特性的基礎(chǔ)上加入某種抗干擾能力，適應(yīng)生存需要。本文通過(guò)希爾伯特-黃變換(HHT)方法證實(shí)了生物進(jìn)化結(jié)構(gòu)與其復(fù)雜程度是一致的，通過(guò)HHT定量提取單細(xì)胞生物和多細(xì)胞生物核小體定位信號(hào)特征，進(jìn)一步驗(yàn)證兩物種之間在核小體組織結(jié)構(gòu)上顯著的差異性。