基于ChIP-seq數(shù)據(jù)HMM方法識(shí)別全基因組的差異組蛋白修飾位點(diǎn)

1、基于ChIP-seq數(shù)據(jù)HMM方法識(shí)別全基因組的差異組蛋白修飾位點(diǎn)摘要目的：表觀遺傳修飾是調(diào)控基因表達(dá)和基因組功能的一個(gè)主要因素。在不同的表觀遺傳修飾中，差異組蛋白修飾位點(diǎn)（DHMSs）是不同細(xì)胞類型、時(shí)期和環(huán)境影響時(shí)，表觀遺傳動(dòng)態(tài)性質(zhì)和基因表達(dá)調(diào)控的一個(gè)研究熱點(diǎn)。為了測(cè)定全基因組的組蛋白修飾，ChIP-seq技術(shù)是一種有效的方法。因此，通過(guò)比較兩個(gè)ChIP-seq文庫(kù)可以識(shí)別潛在的DHMSs。結(jié)果：我們的目的是識(shí)別DHMSs，提出一種稱為ChIPDiff的方法來(lái)通過(guò)ChIP-seq測(cè)定的數(shù)據(jù)全基因組比對(duì)組蛋白修飾位點(diǎn)?；谟^察的ChIP片段數(shù)，提出了一個(gè)隱馬模型的方法推斷每個(gè)基因組位置的組蛋

2、白修飾變化狀態(tài)。我們通過(guò)比對(duì)小鼠ESC和NPC細(xì)胞的H3K27me3修飾位點(diǎn)來(lái)評(píng)估ChIPDiff的效果。我們證明了此方法確定H3K27me3 的DHMSs具有高靈敏度，特異性和重復(fù)性。進(jìn)一步應(yīng)用ChIPDiff揭示不同細(xì)胞時(shí)期的差異H3K4me3和H3K36me3位點(diǎn)。我們研究中的比對(duì)有很多有趣的生物學(xué)發(fā)現(xiàn)。1.介紹真核DNA是被打包到一個(gè)由周圍環(huán)繞組蛋白的DNA的重復(fù)核小體組成的染色質(zhì)結(jié)構(gòu)。組蛋白可以發(fā)生大量的翻譯后修飾如，甲基化，乙?；?，磷酸化和泛素化。組蛋白修飾影響基因表達(dá)和基因組功能。大量實(shí)驗(yàn)證明一些組蛋白甲基化類型在生物學(xué)過(guò)程中起主要作用。一個(gè)典型的例子是在哺乳動(dòng)物胚胎干細(xì)胞通過(guò)H

3、3K27me3抑制發(fā)育調(diào)控維持干細(xì)胞多能性。在癌癥中也特異的發(fā)現(xiàn)一些表觀遺傳K27干細(xì)胞標(biāo)記。此外，H3K9me3、H3K9me2和癌細(xì)胞中沉默腫瘤抑制基因相關(guān)。因此，特異基因組位置的差異組蛋白修飾密度，文中稱為差異組蛋白修飾位點(diǎn)“DHMS”，在不同細(xì)胞類型，時(shí)期和環(huán)境影響是比較研究的重點(diǎn)。我們可以用染色質(zhì)免疫共沉淀（ChIP）來(lái)測(cè)定組蛋白修飾信號(hào)，抗體用于富集修飾位點(diǎn)的DNA片段。在過(guò)去的幾年開(kāi)發(fā)了幾種基于ChIP的技術(shù)，包括ChIP-chip, ChIP-PET and ChIP-SAGE，用于大規(guī)?；蚪M區(qū)域的組蛋白修飾和轉(zhuǎn)錄因子結(jié)合位點(diǎn)研究。隨著最近超高通量測(cè)序技術(shù)如Illumina/

4、 Solexa GA 測(cè)序的產(chǎn)生，ChIP-seq成為一個(gè)主要的高覆蓋、高分辨率和低成本的方法。ChIP-seq的基本思想是讀取ChIP富集的序列的一端，接著映射這些短讀稱為tag到基因組上以找到這些片段的基因組位置。一個(gè)ChIP文庫(kù)中有百萬(wàn)個(gè)tag標(biāo)簽測(cè)序，形成一個(gè)代表全基因組與組蛋白修飾位點(diǎn)和轉(zhuǎn)錄因子結(jié)合位點(diǎn)的ChIP片段數(shù)的譜。受到ChIP-seq在單個(gè)文庫(kù)識(shí)別組蛋白修飾位點(diǎn)的鼓舞，我們想是否可以通過(guò)計(jì)算的比較不同細(xì)胞類型和實(shí)驗(yàn)條件的兩條ChIP-seq文庫(kù)來(lái)識(shí)別DHMS。Mikkelsen等人測(cè)定了小鼠ESC、NPC和MEF細(xì)胞的H3K4me3 (K4) 和 K27位點(diǎn)，比較三種類型啟

5、動(dòng)子區(qū)域修飾位點(diǎn)的發(fā)生。他們研究的局限在于修飾位點(diǎn)是定性的比較而非定量。一個(gè)例子說(shuō)明了這種局限，K4調(diào)控K1f4,已知其和基因表達(dá)正相關(guān)。K1f4在ESC和NPC啟動(dòng)子定性分析中都標(biāo)記K4，因此不能解釋在ESC的K1f4上調(diào)。另一方面，定量比較表明ESC的K1f4啟動(dòng)子的K4密度比NPC多5倍，這和表達(dá)變化是一致的。據(jù)我們所知，幾乎沒(méi)有全基因組定量比較兩個(gè)ChIP-seq文庫(kù)的文獻(xiàn)。受芯片分析的啟發(fā)，一個(gè)簡(jiǎn)單的解決這個(gè)問(wèn)題的方法是將基因組分為箱bins，計(jì)算每個(gè)binChIP片段數(shù)的倍數(shù)變化。然而，fold-change方法對(duì)由ChIP片段隨機(jī)樣本的技術(shù)變化時(shí)敏感的。本文中，我們提出的方法稱為

6、ChIPDiff通過(guò)考慮連續(xù)bin之間的相關(guān)性改進(jìn)了fold-change方法。我們用隱馬模型建立相關(guān)性，轉(zhuǎn)移概率用一種無(wú)監(jiān)督方式自動(dòng)訓(xùn)練。接下來(lái)通過(guò)訓(xùn)練HMM參數(shù)來(lái)推斷組蛋白修飾狀態(tài)的變化。為了評(píng)估ChIPDiff的性能，我們首先比較Mikkelsen數(shù)據(jù)ESC和NPC的K27文庫(kù)。在全基因組識(shí)別了4277個(gè)k27的DHMS區(qū)域。三個(gè)標(biāo)準(zhǔn)顯示效果是令人滿意的：（a）敏感性：2006年在高度保守的非編碼元件中，80%的從基因表達(dá)推斷的DHMSs被ChIPDiff確定。（b）特異性：基于非細(xì)胞特異性控制比對(duì)，我們估計(jì)識(shí)別的DHMS區(qū)域的假陽(yáng)性率是0.19%。（c）重復(fù)度：檢查兩個(gè)獨(dú)立的子集的結(jié)果

7、的交集，顯示3-4百萬(wàn)個(gè)tags測(cè)序的57.4的DHMSs在技術(shù)上重現(xiàn)，評(píng)價(jià)結(jié)果還表明，在所有三個(gè)方面的定性分析，該方法優(yōu)于fold-change的方法。我們進(jìn)一步應(yīng)用ChIPDiff到H3K4me3（K4）和H3K36me3（K36），發(fā)現(xiàn)這兩種類型組蛋白修飾的DHMSs和研究了他們?cè)诟杉?xì)胞分化潛在的生物的作用。研究中有幾個(gè)有趣的生物學(xué)發(fā)現(xiàn)。2.方法2.1確定組蛋白修飾位點(diǎn)給定來(lái)個(gè)ChIP-seq文庫(kù)，L1和L2，識(shí)別DHMSs的第一步是確定L1和L2的組蛋白修飾假定的位點(diǎn)。這部分詳述這一步。ChIP-seq實(shí)驗(yàn)產(chǎn)生的原始數(shù)據(jù)的tags被映射到基因組，獲得它們的位置和方向。由于ChIP-se

8、q實(shí)驗(yàn)的PCR過(guò)程，大量的tags可能源于一個(gè)單一的ChIP片段。為了移除這一重復(fù)性，映射到相同位置和相同方向的tags被作為一個(gè)單一的copy。注意到在ChIP-seq協(xié)議一個(gè)單一的tag是通過(guò)測(cè)序一個(gè)ChIP片段的末端得到的，平均長(zhǎng)度是200bp。因此我們通過(guò)其方向的100bp轉(zhuǎn)移tag的位置近似估計(jì)響應(yīng)ChIP片段的中心。全基因組被分成1k-bp的bin，計(jì)算每個(gè)bin的ChIP片段中心數(shù)。預(yù)處理過(guò)程之后，產(chǎn)生ChIP片段數(shù)譜?？紤]到基因組有m個(gè)bin，譜L1和L2分別表示為X1 = x1.1，x1.2，.x1.m 和X2 = x2.1，x2.2，.x2.m。其中xij是在Li中第j個(gè)b

9、in的片段數(shù)。為了描述每個(gè)bin中片段的結(jié)合富集，我們定義F值標(biāo)準(zhǔn)化測(cè)序的深度：其中n1和n2是L1和L2測(cè)序片段的總數(shù)，如圖。 Mikkelson et al.(2007)和Robertson et al.(2007)指出有與重復(fù)序列區(qū)域的存在，并不是所有的bin都能在tag映射程序中檢測(cè)到。讓記為基因組“有效”的bin，分值F的期望在有效bin時(shí)是F(i)/(m×)，等于2/(m×)。Mikkelson et al.(2007)估計(jì)小鼠基因組的等于0.7。如果一個(gè)bin的F值大于2/(m×)，我們標(biāo)記其為一個(gè)推測(cè)的組蛋白修飾位點(diǎn)。1k bp內(nèi)的連續(xù)修飾位點(diǎn)彼此

10、分開(kāi)被合并為組蛋白修飾區(qū)域。2.2用Fold-change方法定量的比較修飾強(qiáng)度為了便于定義和描述，文章其他部分將介紹的基于推定的組蛋白修飾區(qū)域在2.1介紹，假設(shè)一個(gè)區(qū)域包含k個(gè)bin，我們定義L1和L2的ChIP片段數(shù)分別為x1.i，x2.i，在區(qū)域的第i個(gè)bin（i=1，1，k）。組蛋白修飾表現(xiàn)出對(duì)各種動(dòng)力性和化學(xué)計(jì)量性。對(duì)一個(gè)ChIP實(shí)驗(yàn)，我們定義文庫(kù)Lj的第i個(gè)bin的修飾強(qiáng)度是任意ChIP片段來(lái)自ChIP過(guò)程第i個(gè)bin的概率，定義為pj，i。由于提取和測(cè)序ChIP片段是一個(gè)隨機(jī)抽樣過(guò)程，文庫(kù)Lj的第i個(gè)bin的觀察片段xj，i的后驗(yàn)概率，強(qiáng)度的條件概率pj，i，近似服從二項(xiàng)分布：（

11、1）我們接下來(lái)估計(jì)先驗(yàn)概率pj，i服從beta分布：（2）B(,)是beta函數(shù)。注意到beta分布先于二項(xiàng)是共軛的，所以條件概率也服從beta分布，期望等于。在我們的應(yīng)用中，參數(shù)和設(shè)為1和m，m是基因組中bin的總數(shù)（詳見(jiàn)補(bǔ)充方法）。我們定義一個(gè)DHMS，當(dāng)一個(gè)bin內(nèi)L1和L2的強(qiáng)度比值大于（L1富集DHMS）或者小于1/(L2富集DHMS)。是一個(gè)預(yù)先確定的閾值，值1。一個(gè)簡(jiǎn)單識(shí)別DHMSs的方法是估計(jì)ChIP片段數(shù)的期望強(qiáng)度（更好的是對(duì)數(shù)比）的倍數(shù)變化，如下：（3）基于方程（3）的對(duì)數(shù)比估計(jì)顯示圖1（a）。fold-change法的一個(gè)缺陷是由于隨機(jī)抽樣引起技術(shù)差異。圖1（b）顯示一個(gè)

12、RI-plot描述了依據(jù)強(qiáng)度的log比值變化。當(dāng)強(qiáng)度相對(duì)較小，log值的變化太高，這可能引起大量的假陽(yáng)性。2.3一個(gè)基于隱馬模型的方法識(shí)別DHMSs組蛋白修飾通常發(fā)生在連續(xù)區(qū)域范圍是幾百甚至上千個(gè)核苷酸。因此可以期望連續(xù)的bin測(cè)量的強(qiáng)度變化可能強(qiáng)相關(guān)。通過(guò)觀察ChIP-seq譜支持這一觀點(diǎn)。例如，圖1（a）的log比值譜的自相關(guān)是0.84。在ChIP-chip數(shù)據(jù)分析中，Li et al.(2005)年設(shè)計(jì)的HMM模型構(gòu)建連續(xù)探針之間的信號(hào)相關(guān)成功的應(yīng)用于識(shí)別p53結(jié)合位點(diǎn)，表示HMM在我們研究中應(yīng)用的潛在可能性。在此我們提出一個(gè)基于HMM的方法，ChiPDiff來(lái)解決這一問(wèn)題。我們定義Si

13、為第i個(gè)bin的組蛋白修飾變化狀態(tài)（i=1到k），基于2.2對(duì)于DHMS的定義，狀態(tài)Si為以下三個(gè)值之一：0：無(wú)差別位點(diǎn)，if 1/p1，i/p2，i ；1：L1富集DHMS，if p1，i/p2，i；2：L2富集DHMS，if p1，i/p2，i1/。我們建模bin間的相關(guān)性作為一個(gè)一階馬爾可夫鏈Pr( Si|S0,S1，., Si-1)= Pr(Si|Si-1)，S0是區(qū)域內(nèi)第一個(gè)bin前的起始狀態(tài)。一個(gè)HMM實(shí)施是通過(guò)觀察片段數(shù)推斷狀態(tài)的后驗(yàn)概率分布。HMM的三個(gè)特征:起始狀態(tài)S0的先驗(yàn)概率，emission發(fā)射概率，和狀態(tài)轉(zhuǎn)移概率。初始狀態(tài)S0采用固定值0，因?yàn)槲覀兗俣▋蓚€(gè)文庫(kù)中區(qū)域起

14、始位置是組蛋白修飾缺乏的基因組位置。我們通過(guò)整合所有可能的Si值的p1，i和p2，i得到emission發(fā)射概率讀者可以參考補(bǔ)充方法的詳細(xì)推導(dǎo)。在等式（4）中，服從二項(xiàng)分布（1），服從beta分布（2）。轉(zhuǎn)移概率列表由Baum-Welch算法訓(xùn)練得到，采用期望最大化（EM）步驟以無(wú)監(jiān)督的方式從隱藏狀態(tài)迭代估計(jì)HMM的參數(shù)。訓(xùn)練過(guò)程中，傳輸參數(shù)初始化是統(tǒng)一的，初始狀態(tài)S0和狀態(tài)傳輸概率如以上描述確定。因?yàn)檗D(zhuǎn)移概率表在整個(gè)基因組是相同的，是通過(guò)所有推定的組蛋白修飾區(qū)域轉(zhuǎn)移頻率累加訓(xùn)練的（train）。在ChiPDiff的最后一步，每個(gè)bin中的概率分布狀態(tài)由forward-backward算法推斷

15、。如果bin的后驗(yàn)概率大于置信閾值（0<<1)當(dāng)Si=1或Si=2定為一個(gè)DHMS區(qū)。連續(xù)的沒(méi)有縫隙的DHMS被合并為一個(gè)DHMS。ChiPDiff最大計(jì)算量的一步是訓(xùn)練轉(zhuǎn)移概率表。兩個(gè)策略可以減少計(jì)算量(a) 訓(xùn)練HMM之前，發(fā)射概率的積分被數(shù)值計(jì)算的而且被編寫(xiě)成一張查詢列表。(b)我們?cè)试S轉(zhuǎn)移概率列表基于從推定組蛋白修飾區(qū)域隨機(jī)選擇子集訓(xùn)練。3.結(jié)果我們應(yīng)用ChIPDiff處理Mikkelson實(shí)驗(yàn)數(shù)據(jù)，ChIPDiff的的性能通過(guò)比較小鼠ESC和NPC的H3K27me3文庫(kù)評(píng)估。我們又應(yīng)用ChIPDiff處理H3K4me3和H3K36me3數(shù)據(jù)發(fā)現(xiàn)了DHMSs而且研究它們?cè)诟?/p>

16、細(xì)胞分化中潛在的生物學(xué)作用。3.1 H3K27me3數(shù)據(jù)評(píng)估選用H3K27me3評(píng)估的原因是因?yàn)樗腄HMSs在高度保守的非編碼元件（HCNEs）已經(jīng)有人研究。而且，K27優(yōu)先標(biāo)記基因區(qū)域功能作為抑制子，這有利于我們利用表達(dá)數(shù)據(jù)間接的驗(yàn)證。我們用ChIPDiff比較ESC和NPC的K27文庫(kù)，fold-change閾值設(shè)為3.0置信閾值為0.95.HMM隨機(jī)訓(xùn)練10000次選定組蛋白修飾區(qū)域，26230bins認(rèn)定為DHMS是，對(duì)應(yīng)于4722連續(xù)區(qū)域。它們中3,833 (81.2%)區(qū)域ESC富集，889 (18.8%)NPC富集，這意味著細(xì)胞分化時(shí)期K27消耗的整體趨勢(shì)。我們首次評(píng)估了ChI

17、PDiff的性能通過(guò)確定其生物學(xué)意義，如敏感性。Bernstein發(fā)現(xiàn)K27在ESC中富集在高度保守的非編碼元件（HCNEs），抑制發(fā)育調(diào)控子的數(shù)量來(lái)維持細(xì)胞的stemness。這些組蛋白標(biāo)記在不同分化細(xì)胞中消失。HCNEs中，我們選擇了223個(gè)基因，Mikkelson研究了它們的表達(dá)。因?yàn)镵27作為功能抑制子，這些中的一些被K27標(biāo)記的HCNEs在NPC中上調(diào)，我們認(rèn)為在這些基因DHMSs被確定。與預(yù)期相同，一個(gè)包含30個(gè)上調(diào)基因的子集被確定，標(biāo)準(zhǔn)化超過(guò)4倍。它們中80%被標(biāo)記的由ChIPDiff識(shí)別的DHMSs在啟動(dòng)子區(qū)域轉(zhuǎn)錄起始位點(diǎn)±1kb。相反，193個(gè)基因中只有19.2%在

18、NPC中被DHMSs標(biāo)記的不上調(diào)。為了檢驗(yàn)ChIPDiff的特異性，我們需要評(píng)估錯(cuò)誤識(shí)別的不是細(xì)胞特異的DHMS區(qū)域的片段。針對(duì)這一目的，我們將這一文庫(kù)分為兩個(gè)技術(shù)復(fù)制本：Lesc，rep1和Lesc，rep2，Lnpc，rep1和Lnpc，rep2。復(fù)制樣本的tag組成取自ChIP-seq實(shí)驗(yàn)的不同通道，有相似的測(cè)序depth（見(jiàn)補(bǔ)充表四復(fù)制本的通道）。通過(guò)合并tags產(chǎn)生兩個(gè)新的文庫(kù)分別是Lesc，rep1和Lnpc，rep1，Lesc，rep2和Lnpc，rep2。因?yàn)閺?fù)制本有相似的測(cè)序depth，兩個(gè)文庫(kù)的差異不是細(xì)胞特異的可能只是實(shí)驗(yàn)技術(shù)變化的影響。比較這些非細(xì)胞特異的控制集，Chi

19、PDiff識(shí)別出9個(gè)差異的區(qū)域，因此我們估計(jì)在識(shí)別細(xì)胞特異比較時(shí)DHMS區(qū)時(shí)假陽(yáng)性率為0.19（9/4722）。我們通過(guò)構(gòu)建兩個(gè)獨(dú)立的細(xì)胞特異比對(duì)途徑passes檢測(cè)重復(fù)性：Lesc，rep1vs.Lnpc，rep1，Lesc，rep2vs.Lnpc，rep2。為了測(cè)量重復(fù)性，我們定義一個(gè)分?jǐn)?shù)作為兩個(gè)passes均識(shí)別的DHMSs數(shù)目與單個(gè)pass識(shí)別DHMSs平均值的比值。結(jié)果得ChIPDiff重復(fù)分?jǐn)?shù)為57.4%。注意到重復(fù)性的條件是重復(fù)本的測(cè)序depth在3到4百萬(wàn)個(gè)tags（補(bǔ)充表4）。為了比較不同方法的效果，我們重復(fù)比較fold-change法和定量方法的敏感性、特異性和重復(fù)度。在定

20、量方法中，ESC和NPC的K27修飾位點(diǎn)單個(gè)識(shí)別用Mikkelson的方法，K27位點(diǎn)只在單個(gè)細(xì)胞類型標(biāo)記識(shí)別為DHMSs。連續(xù)的DHMSs合并為DHMS區(qū)。為了公正的比對(duì)，閾值調(diào)整使所有三種方法DHMSs區(qū)域相似數(shù)目被確定（因?yàn)殚撝颠x的離散值，所以這個(gè)數(shù)目不同）。評(píng)估結(jié)果總結(jié)見(jiàn)表1。ChIPDiff所有三方面均優(yōu)于其他兩種方法。Fold-change和定量的方法都有高的假陽(yáng)性率，表示這些方法對(duì)技術(shù)變化和實(shí)驗(yàn)偏差是敏感的。3.2應(yīng)用H3K4me3和H3K36me3數(shù)據(jù)我們擴(kuò)展我們的數(shù)據(jù)研究到H3K4me3和H3K36me3。這兩種修飾類型以不同的方式正向調(diào)控基因表達(dá)。Guenther2007發(fā)

21、現(xiàn)K4在基因轉(zhuǎn)錄起始標(biāo)記活性啟動(dòng)子，而K36發(fā)生在基因區(qū)作為延伸的標(biāo)志。我們之前的研究（Zhao等人）也顯示K4和K27在人類ESC中建立不同的基因組區(qū)域活性和非活性染色質(zhì)結(jié)構(gòu)。這引起我們研究NPC和ESC之間這些組蛋白修飾的DHMSs的興趣。此外，K4位點(diǎn)通常在ChIP-seq譜轉(zhuǎn)錄起始位點(diǎn)周圍急劇以點(diǎn)狀的模式出現(xiàn)，而K36位點(diǎn)出現(xiàn)以更寬的模式，提供一種全面的檢驗(yàn)-bed評(píng)估我們方法對(duì)于不同組蛋白修飾的適用性。我們用3.1提到的ChiPDiff運(yùn)行相同的文庫(kù)。結(jié)果見(jiàn)表2。連續(xù)的DHMSs合并在一起。值得注意的是，K4在ESC富集的DHMSs遠(yuǎn)遠(yuǎn)大于在K4在NPC富集的DHMSs?？紤]到這種不

22、平衡也出現(xiàn)在K27，我們假設(shè)被K4和K27標(biāo)記的二價(jià)染色質(zhì)結(jié)構(gòu)可能相關(guān)（Bernstein2006）。接下來(lái)的分析中，我們發(fā)現(xiàn)3833個(gè)ESC的K27 DHMSs中的1961（51.2%）個(gè)與K4的DHMSs重疊。相反，K36和K27傾向于互相排斥：只有8個(gè)（0.21%）個(gè)DHMS重疊。為了研究DHMSs和基因表達(dá)的相關(guān)性，我們注釋DHMSs區(qū)的Refseq基因和基因表達(dá)數(shù)據(jù)。為了去除冗余，基因注釋最長(zhǎng)ORF選擇是如果相同的轉(zhuǎn)錄本注釋到相同的基因，結(jié)果篩選一共18795個(gè)唯一的基因。如圖2所示，K4和K36共調(diào)控基因表達(dá)強(qiáng)顯著相關(guān)。這與之前的研究相一致。1085個(gè)基因在ESC中上調(diào)，791（72.9%）和ESC富集K4和K36相關(guān)，表示基因表達(dá)可能潛在的由DHMS預(yù)測(cè)。顯然，ESC兩個(gè)關(guān)鍵的轉(zhuǎn)錄因子，Nanog和Oct4，由K4和K36的DHMSs標(biāo)記，也暗示這些ESC的組蛋白修飾標(biāo)記在影響轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)其主