高通量 RNA 甲基化測(cè)序數(shù)據(jù)處理與分析研究進(jìn)展

1、高通量RNA 甲基化測(cè)序數(shù)據(jù)處理與分析研究進(jìn)展*劉戀1張紹武1*孟佳2陳潤生1,3(1西北工業(yè)大學(xué)自動(dòng)化學(xué)院,信息融合技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,西安710072;2西交利物浦大學(xué)生物科學(xué)系,吳江太湖新城研究院,蘇州215123;3中國科學(xué)院生物物理研究所,北京100101摘要隨著高通量測(cè)序技術(shù)快速發(fā)展,MeRIP-seq (methylated RNA immunoprecipitation sequencing測(cè)序技術(shù)開啟了RNA 表觀遺傳學(xué)研究新局面,能夠在全基因組范圍內(nèi)描述RNA 甲基化.從MeRIP-seq 高通量數(shù)據(jù)中挖掘RNA 甲基化模式,有助于揭示mRNA 甲基化在調(diào)控基因表達(dá)、剪切

2、等方面所發(fā)揮的潛在功能,有效指導(dǎo)癌癥的干預(yù)治療.本文從MeRIP-seq 測(cè)序原理出發(fā),較全面地綜述MeRIP-seq 數(shù)據(jù)處理和分析方法研究現(xiàn)狀,并對(duì)其所面臨的計(jì)算問題進(jìn)行討論和展望.關(guān)鍵詞MeRIP-seq 測(cè)序,數(shù)據(jù)處理與分析,RNA 甲基化,表觀遺傳學(xué)科分類號(hào)Q5,Q6,Q7DOI :10.16476/j.pibb.2015.0078生物化學(xué)與生物物理進(jìn)展Progress in Biochemistry and Biophysics 2015,42(10:891899 *國家自然科學(xué)基金資助項(xiàng)目(91430111,61473232,61401370,61170134.*通訊聯(lián)系人.Te

3、lE-mail:zhangsw 收稿日期:2015-03-23,接受日期:2015-07-01表觀遺傳學(xué),包括組蛋白共價(jià)修飾(covalent histone modification、DNA 甲基化修飾(DNA methylation、RNA 甲基化修飾(RNA methylation、基因組印記(genomic imprinting、基因沉默(gene silencing、RNA 編輯(RNA editing及非編碼RNA (noncoding RNA等,是指在核苷酸序列不發(fā)生改變的情況下,生物表型或基因表達(dá)發(fā)生了穩(wěn)定的可遺傳變化1.RNA 甲基化作為表觀遺傳學(xué)

4、研究的重要內(nèi)容之一,是指發(fā)生在RNA 分子上不同位置的甲基化修飾現(xiàn)象,6-甲基腺嘌呤(N 6-methyladenosine ,m 6A和5-甲基胞嘧啶(C 5-methylcytidine ,m 5C是真核生物中最常見的兩種RNA 轉(zhuǎn)錄后修飾.RNA 甲基化在調(diào)控基因表達(dá)、剪接、RNA 編輯、RNA 穩(wěn)定性、控制mRNA 壽命和降解等方面可能扮演重要角色.相對(duì)于DNA 甲基化,RNA 甲基化更加復(fù)雜、種類繁多、普遍存在于各種高級(jí)生物中2-4.由于缺乏有效檢測(cè)手段,相關(guān)研究多局限于非編碼tRNA 和rRNA ,或小部分編碼轉(zhuǎn)錄片段1,且多數(shù)RNA 甲基化功能未知.隨著高通量測(cè)序技術(shù)的發(fā)展5及一

5、些RNA 甲基化功能的發(fā)現(xiàn)6-11,人們開始關(guān)注RNA 甲基化研究.尤其MeRIP-seq (methylated RNAimmunoprecipitation sequencing高通量測(cè)序技術(shù)的出現(xiàn),能夠高效精確檢測(cè)全轉(zhuǎn)錄組不同的RNA 甲基化,奠定了RNA 甲基化研究基礎(chǔ).如何有效處理和分析MeRIP-seq 技術(shù)生成的海量數(shù)據(jù),是成功發(fā)現(xiàn)RNA 甲基化機(jī)理及功能的關(guān)鍵.本文較全面介紹MeRIP-seq 測(cè)序原理、數(shù)據(jù)處理及分析基本流程、關(guān)鍵方法、現(xiàn)有算法軟件,重點(diǎn)討論MeRIP-seq 數(shù)據(jù)處理和分析過程中所面臨的挑戰(zhàn).1MeRIPseq 技術(shù)測(cè)序原理MeRIP-seq 技術(shù)將甲基化D

6、NA 免疫共沉淀(methylated DNA immunoprecipitation ,MeDIP技術(shù)12、RNA 結(jié)合蛋白免疫共沉淀(RNA immunoprecipitation ,RIP技術(shù)和RNA 測(cè)序(RNA sequencing ,RNA-seq技術(shù)13組合起來,高精度地檢測(cè)全基因組(或全轉(zhuǎn)錄組范圍內(nèi)的RNA 甲基生物化學(xué)與生物物理進(jìn)展Prog.Biochem.Biophys.2015;42(102MeRIPseq 測(cè)序文庫制備和測(cè)序平臺(tái)數(shù)據(jù)輸出本小節(jié)將針對(duì)Illumina/Solexa 測(cè)序平臺(tái),介紹MeRIP-seq 測(cè)序文庫制備及測(cè)序平臺(tái)數(shù)據(jù)輸出.2.1MeRIPseq 測(cè)

7、序文庫制備MeRIP-seq 測(cè)序文庫制備過程如下:首先從樣本細(xì)胞組織中分離出RNA ,考慮到總RNA 中含有大量的rRNA 序列,因此需要結(jié)合不同的方法去除其中的rRNA .對(duì)于真核生物而言,常采用Poly(T寡核苷酸提取出帶Poly(A的RNA 去除rRNA ;而對(duì)不含Poly(A尾的轉(zhuǎn)錄本序列以及存在部分降解的總RNA 樣本而言,需要試劑盒去除rRNA ,從而得到除rRNA 外的全部RNA ,然后將提取出的RNA 隨機(jī)打斷.MeRIP-seq 技術(shù)對(duì)帶有甲基化修飾的片段(IP 樣本進(jìn)行測(cè)序時(shí),需要平行對(duì)一個(gè)對(duì)化.MeRIP-seq 技術(shù)采用免疫共沉淀方法,即甲基化RNA 特異性抗體與被隨

8、機(jī)打斷的RNA 片段進(jìn)行孵育,抓取有甲基化修飾的片段進(jìn)行測(cè)序;同時(shí)需要平行測(cè)序一個(gè)對(duì)照(control樣本,對(duì)照樣本用于消除抓取帶有甲基化片段過程中的背景.然后將免疫共沉淀(IP樣本和對(duì)照樣本中的序列片段對(duì)比(或定位到參考基因組/轉(zhuǎn)錄組上,檢測(cè)RNA 甲基化位點(diǎn).對(duì)照樣本測(cè)量對(duì)應(yīng)RNA 的表達(dá)量,本質(zhì)上是RNA-seq 數(shù)據(jù).圖1為MeRIP-seq 技術(shù)檢測(cè)m 6A RNA 甲基化過程示意圖.MeDIP-seq 和ChIP-seq 測(cè)序技術(shù)均是將免疫共沉淀與測(cè)序相結(jié)合.MeRIP-seq 技術(shù)主要應(yīng)用于RNA 甲基化研究,而ChIP-seq 、MeDIP-seq 主要應(yīng)用于DNA 甲基化研究

9、.MeRIP-seq 技術(shù)要求必須有對(duì)照樣本,而MeDIP-seq 和ChIP-seq 技術(shù)對(duì)于對(duì)照樣本沒有要求.表1為MeRIP-seq 、MeDIP-seq 和ChIP-seq 三種測(cè)序技術(shù)對(duì)比.Fig.1The work flow of detecting m 6A RNA methylationusing MeRIPseq technology圖1MeRIPseq 技術(shù)檢測(cè)m 6A RNA 甲基化過程輸入的Control 樣本Control 樣本IP 樣本凈化&測(cè)序片段(100ntRNA 樣本使用m 6A 抗體得到的IP 樣本AAAAAAAAAAAAAAAAAAGGAAAAGG

10、CGCGC AAAGGCGCGCGTTATCTTGGAAAA TCAGAAAAGGCGC TTCTAACTTGGATable 1Comparison of MeRIPseq,MeDIPseq,ChIPseq sequencing technologies 表1MeRIPseq 、MeDIPseq 、ChIPseq 三種測(cè)序技術(shù)對(duì)比ChIP-seqMeDIP-seq MeRIP-seq 研究對(duì)象化學(xué)修飾化學(xué)修飾化學(xué)修飾分子DNA DNA RNA 比對(duì)器非拼接非拼接拼接特征蛋白質(zhì)綁定位點(diǎn)或峰CpG 島甲基化位點(diǎn)或峰量化相對(duì)量(與絕對(duì)量線性相關(guān)相對(duì)量(與絕對(duì)量線性相關(guān)相對(duì)量(與絕對(duì)量不相關(guān)處理流程一

11、般不做峰檢測(cè)代表軟件MACS 14,CisGenome 15Batman 16,MeQA 17exomePeak 18,MeRIP-PF 19差異分析僅需要免疫沉淀樣本僅需要免疫沉淀樣本需要免疫沉淀樣本和對(duì)照樣本模體雙鏈雙鏈鏈特異性顯著性分析Peak 富集分析讀段定位標(biāo)簽平移序列延伸Peaks 的顯著性分析在外顯子上進(jìn)行讀段定位比較2個(gè)泊松分布的均值,確定Peaks每個(gè)滑窗中對(duì)IP 和control 樣本的reads 數(shù)泊松建模892··劉戀,等:高通量RNA 甲基化測(cè)序數(shù)據(jù)處理與分析研究進(jìn)展2015;42(10照樣本(Control 樣本進(jìn)行測(cè)序,其IP 樣本和Contr

12、ol 樣本的片段選擇方法主要有以下2種:a .將打斷的RNA 片段分成兩份,一份直接用于制備Control 樣本的cDNA 文庫,另一份采用抗m 6A 抗體與被打斷的RNA 進(jìn)行孵育,抓取帶有m 6A 修飾的片段,用于制備IP 樣本的cDNA 文庫.由于測(cè)序得到的結(jié)果不以所有RNA 片段為背景,稱這樣得到的IP 樣本和Control 樣本是非成對(duì)的(unpair,在進(jìn)行數(shù)據(jù)處理時(shí)需先對(duì)Control 樣本進(jìn)行處理.b .取兩份相同的RNA 進(jìn)行打斷,其中一份所有的RNA 片段都進(jìn)行測(cè)序,作為Control 樣本,另一份采用抗m 6A 抗體抓取帶有m 6A 修飾的片段進(jìn)行測(cè)序作為IP 樣本.由于

13、測(cè)序得到的IP 樣本背景為當(dāng)前測(cè)序得到的Control 樣本,稱這樣得到的IP 樣本和Control 樣本是成對(duì)的(pair,可直接用于數(shù)據(jù)處理.獲取測(cè)序片段后(包括IP 樣本測(cè)序片段和Control 樣本測(cè)序片段,用隨機(jī)引物和反轉(zhuǎn)錄酶從RNA 片段合成雙鏈cDNA .然后,對(duì)合成的cDNA 進(jìn)行末端修復(fù)并在3端加“A ”,使用特定測(cè)序接頭(adapter連接cDNA 片段兩端,從而得到用于測(cè)序的cDNA .通常情況下,為了得到更高的測(cè)序效率,一般采用電泳切膠法獲取一定長度的cDNA ,再對(duì)其進(jìn)行PCR 擴(kuò)增,得到所需的cDNA 文庫20.2.2測(cè)序平臺(tái)數(shù)據(jù)輸出將制備好的測(cè)序文庫放入測(cè)序平臺(tái)的

14、各通道(lane,通過橋式擴(kuò)增,形成數(shù)以億計(jì)的簇,開始測(cè)序.測(cè)序時(shí),將4種聚合酶加入到單分子陣列中,每個(gè)被加入熒光標(biāo)記的核苷釋放出相對(duì)應(yīng)的熒光.測(cè)序儀通過捕獲熒光標(biāo)記核苷酸所釋放的熒光信號(hào),利用計(jì)算機(jī)軟件確定測(cè)得的堿基及順序,根據(jù)測(cè)序順序連成讀段(read/fragment,輸出以FASTQ 格式記錄讀段序列及測(cè)序質(zhì)量分?jǐn)?shù).在FASTQ 文件中,每4行為一個(gè)讀段,其中第1行以“”開頭,后面是reads 的ID 以及其他信息,第2行為測(cè)序得到的read 的堿基序列,第3行以“+”開頭,跟隨著該read 的名稱(一般與后面的內(nèi)容相同,但有時(shí)可以省略,而“+”一定不能省,第4行代表reads 的質(zhì)量

15、21.3MeRIPseq 測(cè)序數(shù)據(jù)處理MeRIP-seq 技術(shù)主要用于mRNA 甲基化檢測(cè),其測(cè)序數(shù)據(jù)處理主要包括讀段定位、峰檢測(cè)(peak calling、差異甲基化檢測(cè)及剪接異構(gòu)體層次的相關(guān)處理.圖2為MeRIP-seq 測(cè)序數(shù)據(jù)處理流程. Fig.2The process of treating the MeRIPseq data 圖2MeRIPseq 測(cè)序數(shù)據(jù)處理流程(a單樣本MeRIP-seq 測(cè)序數(shù)據(jù)處理.(b雙樣本MeRIP-seq 數(shù)據(jù)比較分析流程.(a(b基因注釋可視化及讀段注釋讀段定位基因?qū)哟紊系姆鍣z測(cè)剪接異構(gòu)體層次上的峰檢測(cè)后續(xù)處理MeRIP-seq 測(cè)序數(shù)據(jù)(Contr

16、ol 樣本MeRIP-seq 測(cè)序數(shù)據(jù)(IP 樣本MeRIP-seq 測(cè)序數(shù)據(jù)(樣本AMeRIP-seq 測(cè)序數(shù)據(jù)(樣本BMeRIP-seq 數(shù)據(jù)處理基本流程MeRIP-seq 數(shù)據(jù)處理基本流程其他高層數(shù)據(jù)分析基因?qū)哟蔚牟町惣谆治黾艚赢悩?gòu)體層次的差異甲基化分析893··生物化學(xué)與生物物理進(jìn)展Prog.Biochem.Biophys.2015;42(103.1讀段定位獲得Control及IP兩樣本測(cè)序數(shù)據(jù)后,首先對(duì)讀段數(shù)據(jù)進(jìn)行預(yù)處理(如將測(cè)序質(zhì)量較差的讀段過濾,然后將兩個(gè)樣本的所有讀段序列映射(mapping定位到參考基因組上,這是后續(xù)數(shù)據(jù)處理和分析的基礎(chǔ).目前,RNA數(shù)

17、據(jù)的讀段定位算法主要采用以下三種技術(shù)22:空位種子索引(spaced-seed indexing、Burrows-Wheeler轉(zhuǎn)換(Burrows-Wheeler transform,BWT、Smith-Waterman動(dòng)態(tài)規(guī)劃23.空位種子索引算法基本原理:將讀段切成片段,形成種子片段,從中選取一部分作為種子建立索引,然后利用查找、延伸等方法來定位讀段.其代表軟件包括MAQ24、ZOOM25、RMAP26. BWT算法基本原理:通過B-W轉(zhuǎn)換對(duì)參考基因組進(jìn)行一次有規(guī)律的重新排序并建立索引,然后利用查找和回溯定位等方法進(jìn)行讀段定位.在查找過程中,可以利用堿基替代來實(shí)現(xiàn)允許的錯(cuò)配.其代表軟件包

18、括Bowtie27、BWA28、SOAP229. Smith-Waterman動(dòng)態(tài)規(guī)劃算法基本原理:利用初始條件和迭代關(guān)系計(jì)算兩條序列所有可能的比對(duì)分值,對(duì)相同位點(diǎn)加分,不同位點(diǎn)減分.采用空隙懲罰機(jī)制處理片段中存在的間隙,并將結(jié)果存放于一個(gè)矩陣中,利用動(dòng)態(tài)規(guī)劃方法回溯尋找最優(yōu)比對(duì)結(jié)果.其代表軟件有BFAST30、SHRiMR31.MeRIP-seq測(cè)序數(shù)據(jù)實(shí)際上是一種RNA讀段數(shù)據(jù),讀段定位時(shí)需要進(jìn)行拼接定位,且讀段定位中會(huì)面臨跨越兩個(gè)外顯子結(jié)合區(qū)域的定位問題.為解決此問題,人們采用以下三種方法進(jìn)行RNA讀段定位:a.基于已知剪接點(diǎn)的比對(duì)定位.該方法在已知基因注釋信息基礎(chǔ)上實(shí)現(xiàn),剪接點(diǎn)在已知接

19、合區(qū)域數(shù)據(jù)庫中可檢測(cè)到.此類方法不能確定新的剪接點(diǎn).代表性軟件工具包括SpliceSeq32、SAMMate33.b.從頭拼接比對(duì)定位.此方法不需已知的注釋信息,且允許新剪接點(diǎn)的檢測(cè).代表性軟件工具包括MapSplice34、SpliceMap35.c.使用注釋信息進(jìn)行從頭拼接的比對(duì)定位.代表性軟件工具包括TopHat36、STAR37.TopHat軟件首先采用Bowtie比對(duì)非拼接的讀段,然后采用Maq組裝已比對(duì)的讀段形成序列的島;在島嶼序列中,TopHat 根據(jù)之前未映射的讀段、可能的標(biāo)準(zhǔn)供體以及接受位點(diǎn)來確定剪接點(diǎn).讀段定位后,通常采用SAM38或BAM文件存儲(chǔ).BAM格式是對(duì)SAM文件

20、的壓縮,可以將SAM格式壓縮到接近原來的20%.SAMTools38、BEDTools39、IGV40為SAM和BAM文件常用處理軟件.3.2峰檢測(cè)算法IP樣本中甲基化位點(diǎn)抓取的讀段較多,將其映射到參考基因組上,會(huì)在甲基化位點(diǎn)附近形成一個(gè)讀段富集區(qū)(enrichment region或者一個(gè)“峰(peak,因而甲基化富集點(diǎn)檢測(cè)算法稱之為峰檢測(cè)(peak calling算法41.峰檢測(cè)過程中,經(jīng)常遇到兩種比較特別的讀段:一種為同一個(gè)讀段可映射到基因組的多個(gè)位置上,稱之為“多映射讀段(multimapping reads”;另一種為一些完全相同的讀段,稱之為“復(fù)制讀段(duplicated rea

21、ds”,該類讀段可能是由PCR擴(kuò)增引起的.對(duì)于“多映射讀段”,常采用下面2種方法處理:a.在不同位置根據(jù)周圍區(qū)域情況按比例分配;b.完全刪除這種讀段,這是最簡單并且最有效的方法42.對(duì)于“復(fù)制讀段”,采用SAMTools38處理.MeRIP-PF19和exomePeak18是目前檢測(cè)MeRIP-seq數(shù)據(jù)讀段富集區(qū)的兩個(gè)主要工具. MeRIP-PF首先將IP樣本數(shù)據(jù)及對(duì)照樣本數(shù)據(jù)映射到參考基因組上,并把參考基因組分割成25bp的固定窗口,通過比對(duì)該窗口上IP樣本和Control樣本的讀段(read數(shù)目,確定m6A甲基化區(qū).但該MeRIP-PF19以固定窗口分割參考基因組,對(duì)于跨窗口的“峰”及跨

22、外顯子的“峰”不能有效地處理,假陽性較高.exomePeak18采用Przyborowski43和Wilenski43-44方法比較兩個(gè)泊松分布的均值(或C-test,對(duì)特定基因外顯子集合進(jìn)行峰檢測(cè),可檢測(cè)跨越外顯子連接區(qū)域的峰,該方法可有效解決轉(zhuǎn)錄豐度問題.由于基因剪接異構(gòu)體的多樣性, exomePeak算法沒有考慮轉(zhuǎn)錄復(fù)雜性,所涉及的諸如平移(shifting、延伸(extension、平滑(smoothing、檢測(cè)等計(jì)算操作相對(duì)直接簡單.盡管exomePeak 算法目前存在這些不足,但該算法仍可以較好地檢測(cè)RNA甲基化位點(diǎn),并對(duì)其進(jìn)行注釋.3.3差異甲基化檢測(cè)基于MeRIP-seq數(shù)據(jù)進(jìn)

23、行差異甲基化檢測(cè),有助于確定2種實(shí)驗(yàn)/顯性條件(如正常和癌癥下的mRNA表觀遺傳調(diào)控差異.ChIP-Seq數(shù)據(jù)與MeRIP-seq數(shù)據(jù)的差異甲基化檢測(cè)有其本質(zhì)區(qū)別.在ChIP-Seq數(shù)據(jù)中,由于DNA總數(shù)在兩種情況下(加刺激、未加刺激是相同的,那么修飾DNA分子的百分比與其數(shù)量保持相同的變化趨勢(shì),因此無論使用相對(duì)量(百分比還是絕對(duì)量,其差異是一致894··劉戀,等:高通量RNA甲基化測(cè)序數(shù)據(jù)處理與分析研究進(jìn)展2015;42(104MeRIPseq數(shù)據(jù)處理面臨的生物信息學(xué)挑戰(zhàn)MeRIP-seq技術(shù)為RNA表觀遺傳學(xué)開啟了新的研究領(lǐng)域,但數(shù)據(jù)分析及處理方法的發(fā)展滯后于實(shí)驗(yàn)技術(shù)的

24、進(jìn)步,現(xiàn)有DNA甲基化數(shù)據(jù)分析及處理方法不能直接用來分析RNA甲基化數(shù)據(jù),急需在以下幾方面發(fā)展有效的計(jì)算方法,分析MeRIP-seq高通量RNA甲基化數(shù)據(jù).4.1甲基化位點(diǎn)預(yù)測(cè)與ChIP-Seq數(shù)據(jù)類似,基于MeRIP-seq數(shù)據(jù)的RNA甲基化位點(diǎn)預(yù)測(cè)需要消除背景讀段分布噪聲,如GC含量、映射能力、抗體非特異性結(jié)合、局部拷貝數(shù)變異等因素引起的實(shí)驗(yàn)誤差和測(cè)序誤差.ChIP-Seq數(shù)據(jù)的背景偏差相對(duì)較小,其轉(zhuǎn)錄因子或DNA甲基位點(diǎn)預(yù)測(cè)不需要對(duì)照樣本,僅通過估計(jì)鄰居基因組區(qū)域的背景就可實(shí)現(xiàn)DNA甲基化位點(diǎn)預(yù)測(cè)45-46.與此相反,由于mRNA片段轉(zhuǎn)錄豐度變化較大及其在3和5端的衰減,MeRIP-se

25、q 數(shù)據(jù)的背景讀段分布變化非常大,必須通過對(duì)照樣本測(cè)量背景轉(zhuǎn)錄豐度.因此MeRIP-seq數(shù)據(jù)甲基化位點(diǎn)預(yù)測(cè)需要檢測(cè)相對(duì)于對(duì)照樣本轉(zhuǎn)錄豐度的IP 樣本“富集峰(peak enrichment”.因而,mRNA 甲基化位點(diǎn)檢測(cè)與常用的DNA甲基化位點(diǎn)檢測(cè)有本質(zhì)上區(qū)別.另外,當(dāng)RNA甲基化位點(diǎn)處于外顯子連接區(qū)的.但在MeRIP-seq數(shù)據(jù)中,由于mRNA差異表達(dá)影響,MeRIP-seq數(shù)據(jù)的背景(如mRNA轉(zhuǎn)錄豐度差異較大.有可能同時(shí)出現(xiàn)“過甲基化(hypermethylation”和“甲基化RNA總量下降”情況,如圖3所示.在DNA甲基化中,未加刺激的情況下,3個(gè)DNA分子中有2個(gè)被修飾,而加刺

26、激情況下,3個(gè)DNA分子中只有1個(gè)被修飾.在修飾的DNA分子質(zhì)量下降的同時(shí),其百分比也是下降的.但在RNA甲基化中,未加刺激的情況下,4個(gè)RNA分子中有2個(gè)被甲基化,而在加刺激情況下,僅有1個(gè)RNA分子,且被修飾.即相對(duì)于未加刺激情況下的RNA甲基化,加刺激情況下的RNA甲基化數(shù)量雖然減少了,但其RNA甲基化百分比卻增加.圖3表明:由于DNA總量保持不變,甲基化DNA總量和其在總DNA中的相對(duì)量保持相同的變化趨勢(shì);由于RNA總量可能變化,甲基化RNA總量和相對(duì)量的變化可能完全不同.另外,圖3中所示帶有甲基化的RNA在加刺激中的總量雖然下降,可是其相對(duì)量卻上升,表明了一種過甲基化現(xiàn)象,同時(shí)RNA

27、表達(dá)量下調(diào)了.exomePeak18工具包含差異甲基化區(qū)域檢測(cè)功能,其檢測(cè)原理基于超幾何測(cè)試計(jì)算兩種情況下的峰值富集顯著性差異,且與一般情況下的ChIP-Seq和RNA-seq計(jì)算的絕對(duì)峰值差異不同.Fig.3The difference of DNA andRNA differential methylation圖3DNA與RNA 差異甲基化區(qū)別DNA RNA 修飾DNA RNA加刺激后未加刺激895··· 896 · 生物化學(xué)與生物物理進(jìn)展 Prog. Biochem. Biophys. 2015; 42 (10 附近時(shí)，“峰”將跨越外顯子連接區(qū) 因

28、此 RNA 甲基化位點(diǎn)預(yù)測(cè)算法需要確定跨越 2 個(gè)或多個(gè)外顯 ，否則，當(dāng)采用現(xiàn)有諸如用于 ChIP-Seq 子的“峰” 數(shù)據(jù)的 MACS14峰檢測(cè)算法時(shí)，會(huì)錯(cuò)誤檢測(cè)出多個(gè) 孤立 “峰” 雖然 exomePeak 18能夠?qū)崿F(xiàn)跨越外顯子連接的 RNA 甲基化位點(diǎn)檢測(cè)，但 exomePeak 并沒有完全 解決上述 MeRIP-seq 所存問題 由于 exomePeak 通過泊松模型計(jì)算讀段數(shù)目，沒有考慮生物學(xué)差 異，會(huì)遺漏過離散的讀段 因此，需要發(fā)展新的 RNA 甲基化位點(diǎn)檢測(cè)算法， 以更加準(zhǔn)確地 進(jìn)行 RNA 甲基化位點(diǎn)檢測(cè) 4.2 基因剪接異構(gòu)體層次上 mRNA 甲基化預(yù)測(cè) 眾所周知，高等真核

29、生物中，通過可變性剪接 相同的基因會(huì)被轉(zhuǎn)錄成不同的異構(gòu)體 (isoform 47， 在基因剪接異構(gòu)體上也會(huì)發(fā)生 RNA 甲基化 當(dāng) IP 樣本中的一個(gè)峰處于異構(gòu)體共享外顯子時(shí)，進(jìn)行 RNA 甲基化位點(diǎn)檢測(cè)前，需對(duì)峰讀段進(jìn)行去卷積 運(yùn)算，確定每個(gè)異構(gòu)體的相對(duì)貢獻(xiàn) 另外，還需要 確定對(duì)照樣本中異構(gòu)體表達(dá)的數(shù)量及它們的相應(yīng)豐 度 總之，如何應(yīng)用 RNA-seq 對(duì)照數(shù)據(jù)預(yù)測(cè)不同 異構(gòu)體甲基化位點(diǎn)，是 MeRIP-seq 數(shù)據(jù)分析中一個(gè) 迫切需要解決的挑戰(zhàn)性問題 4.3 基因及其剪接異構(gòu)體層次上的 mRNA 差異甲 基化預(yù)測(cè) 不同實(shí)驗(yàn)條件下， ChIP-Seq 數(shù)據(jù)的背景 ( 基因 組 DNA 通常非

30、常相似，而由于 mRNA 的差異表 達(dá)， MeRIP-seq 數(shù)據(jù)的背景 (mRNA 轉(zhuǎn)錄豐度 差異 較大 因而，現(xiàn)有適合于 ChIP-Seq 數(shù)據(jù)的差異分 析算法 48，不能直接用來比較兩個(gè) IP 樣本中的讀 段數(shù) 需要研究包括相應(yīng) RNA-seq 對(duì)照樣本的新 計(jì)算框架，比較富集峰的相對(duì)數(shù)量 另外，針對(duì)某 一轉(zhuǎn)錄異構(gòu)體，需要研究有效的算法檢測(cè)其差異甲 基化 4.4 基于分子網(wǎng)絡(luò)的 RNA 甲基化功能注釋 RNA 甲基化可通過調(diào)控基因表達(dá)而實(shí)施重要 生物學(xué)功能，但 RNA 甲基化如何調(diào)控基因、究竟 有哪些生物學(xué)功能，目前缺乏深入研究 我們可通 過整合其他組學(xué)數(shù)據(jù)、構(gòu)建與 RNA 甲基化相關(guān)的

31、 分子網(wǎng)絡(luò)，采用相關(guān)的分子動(dòng)態(tài)信息及網(wǎng)絡(luò)分析方 法 49-53，研究 RNA 甲基化的基因調(diào)控機(jī)制及其所 發(fā)揮的生物學(xué)功能 但如何與其他組學(xué)數(shù)據(jù)整合、 如何構(gòu)建 RNA 甲基化分子網(wǎng)絡(luò)及如何挖掘分析也 是目前急需解決的挑戰(zhàn)性問題 除上述迫切需要解決的問題之外，RNA-seq 分 析中諸如多種讀段、轉(zhuǎn)錄水平上的測(cè)序變化及比對(duì) 偏差等因素，對(duì)于 MeRIP-seq 甲基化峰的檢測(cè)同樣 重要，而 ChIP-Seq 數(shù)據(jù)分析方法中不需要考慮這 些因素 另外，發(fā)展一些有效的方法將 RNA 甲基 化數(shù)據(jù)與其他組學(xué)數(shù)據(jù)進(jìn)行整合，深入研究 RNA 甲基化機(jī)理及其生物學(xué)功能也是生物信息學(xué)今后的 一個(gè)重要研究方向

32、 因而迫切需要發(fā)展新的針對(duì) MeRIP-seq 數(shù)據(jù)的分析方法和計(jì)算工具解決上述問 題，促進(jìn)表觀轉(zhuǎn)錄組學(xué)這一新興領(lǐng)域的快速發(fā)展 5 總結(jié)與展望 RNA 甲基化在調(diào)控基因表達(dá)、剪接、 RNA 編 輯、RNA 穩(wěn)定性、控制 mRNA 的壽命和降解等方 面可能扮演重要角色，其甲基化機(jī)理、位點(diǎn)預(yù)測(cè)和 差異表達(dá)研究，有助于進(jìn)一步揭示細(xì)胞發(fā)育、疾病 等生物學(xué)現(xiàn)象，幫助藥物研發(fā)者設(shè)計(jì)出能夠調(diào)節(jié)基 因表達(dá)、殺死或控制疾病細(xì)胞的小分子 本文從 MeRIP-seq 高通量測(cè)序技術(shù)出發(fā)，首先介紹此技術(shù) 測(cè)序原 理， 在技 術(shù) 特 征 和 數(shù) 據(jù) 處 理 流 程 方 面 與 MeDIP-seq、 ChIP-seq 2

33、 種高通量測(cè)序技術(shù)進(jìn)行了 對(duì)比，然后對(duì) MeRIP-seq 高通量測(cè)序數(shù)據(jù)的讀段定 位、峰檢測(cè)、差異甲基化檢測(cè)及剪接異構(gòu)體等相關(guān) 處理方法進(jìn)行歸納總結(jié)，最后，對(duì) RNA 甲基化位 點(diǎn)檢測(cè)、剪接異構(gòu)體層次上的甲基化位點(diǎn)檢測(cè)、 RNA 差異甲基化分析及基于分子網(wǎng)絡(luò)的 RNA 甲基 化功能注釋所面臨的生物信息學(xué)挑戰(zhàn)問題進(jìn)行了展 望 希望本文能夠?qū)φ诨蚣磳⒉捎?MeRIP-seq 實(shí) 驗(yàn)進(jìn)行科學(xué)研究的學(xué)者和 MeRIP-seq 高通量數(shù)據(jù)處 理研究者提供參考 參 考 文 獻(xiàn) 1 Fu Y, He C. Nucleic acid modifications with epigenetic signif

34、icance. Current Opinion in Chemical Biology, 2012, 16(5: 516-524 2 Desrosiers R, Friderici K, Rottman F. Identification of methylated nucleosides in messenger RNA from Novikoff hepatoma cells. Proc Natl Acad Sci USA, 1974, 71(10: 3971-3975 3 Harris R A, Wang T, Coarfa C, et al. Comparison of sequenc

35、ingbased methods to profile DNA methylation and identification of monoallelic epigenetic modifications. Nature Biotechnology, 2010, 28(10: 1097-1105 4 Dubin D T, Taylor R H. The methylation state of poly A-containing-messenger RNA from cultured hamster cells. Nucleic Acids Research, 1975, 2(10: 1653

36、-1668 5 Meyer K D, Saletore Y, Zumbo P, et al. Comprehensive analysis of mRNA methylation reveals enrichment in 3 UTRs and near stop 2015; 42 (10 劉戀, 等：高通量 RNA 甲基化測(cè)序數(shù)據(jù)處理與分析研究進(jìn)展 22 王 · 897 · codons. Cell, 2012, 149(7: 1635-1646 6 Dominissini D, Moshitch-Moshkovitz S, Salmon-Divon M, et al.

37、Transcriptome-wide mapping of N6-methyladenosine by m6A-seq based on immunocapturing and massively parallel sequencing. Nature Protocols, 2013, 8(1: 176-189 7 Meyer K D, Jaffrey S R. The dynamic epitranscriptome: N6-methyladenosine and gene expression control. Nature Reviews Molecular Cell Biology,

38、2014, 15(5: 313-326 8 Liu J, Yue Y, Han D, et al. A METTL3-METTL14 complex mediates mammalian nuclear RNA N6-adenosine methylation. Nature Chemical Biology, 2013, 10(2: 93-95 9 Ping X L, Sun B F, Wang L, et al. Mammalian WTAP is a regulatory subunit of the RNA N6-methyladenosine methyltransferase. C

39、ell Research, 2014, 24: 177-189 10 Jia G, Fu Y, Zhao X, et al. N6-methyladenosine in nuclear RNA is a major substrate of the obesity-associated FTO. Nature Chemical Biology, 2011, 7(12: 885-887 11 宋述慧, 李語麗, 于 軍. RNA 中 6- 甲基腺嘌呤的研究進(jìn)展. 遺 傳, 2013, 35(12: 1340-1351 Song S H, Li Y L, Yu J. Hereditas, 2013

40、, 35(12: 1340-1351 12 Weber M, Davies J J, Wittig D, et al. Chromosome-wide and promoter-specific analyses identify sites of differential DNA methylation in normal and transformed human cells. Nature Genetics, 2005, 37(8: 853-862 13 Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for tra

41、nscriptomics. Nature Reviews Genetics, 2009, 10(1: 57-63 14 Feng J, Liu T, Qin B, et al. Identifying ChIP-seq enrichment using MACS. Nature Protocols, 2012, 7(9: 1728-1740 15 Ji H, Jiang H, Ma W, et al. An integrated software system for analyzing ChIP-chip and ChIP-seq data. Nature Biotechnology, 20

42、08, 26(11: 1293-1300 16 Down T A, Rakyan V K, Turner D J, et al. A Bayesian deconvolution strategy for immunoprecipitation-based DNA methylome analysis. Nature Biotechnology, 2008, 26(7: 779-785 17 Huang J, Renault V, Sengenes J, et al. MeQA: a pipeline for MeDIP-seq data quality assessment and anal

43、ysis. Bioinformatics, 2012, 28(4: 587-588 18 Meng J, Cui X, Rao M K, et al. Exome-based analysis for RNA epigenome sequencing data. Bioinformatics, 2013, 29 (12: 15651567 19 Li Y, Song S, Li C, et al. MeRIP-PF: An easy-to-use pipeline for high-resolution peak-finding in MeRIP-Seq data. Genomics, Pro

44、teomics & Bioinformatics, 2013, 11(1: 72-75 20 孫 磊, 張 林, 劉 輝. 基于 RNA-Seq 的長非編碼 RNA 預(yù)測(cè). 生物化學(xué)與生物物理進(jìn)展, 2012, 39(12: 1156-1166 Sun L, Zhang L, LIU H. Prog Biochem Biophys, 2012, 39 (12: 1156-1166 21 Cock P J, Fields C J, Goto N, et al. The Sanger FASTQ file format for sequences with quality scores,

45、and the Solexa/Illumina FASTQ variants. Nucleic Acids Research, 2010, 38(6: 1767-1771 曦, 汪小我, 王立坤, 等. 新一代高通量 RNA 測(cè)序數(shù)據(jù)的處 理與分析. 生物化學(xué)與生物物理進(jìn)展, 2010, 37(8: 834-846 Wang X, Wang X W, Wang L K, et al. Prog Biochem Biophys, 2010, 37(8: 834-846 23 楊 燁, 劉 娟. 第二代測(cè)序序列比對(duì)方法綜述. 武漢大學(xué)學(xué)報(bào): 理學(xué)版, 2012, 58(5: 463-470 Yan

46、g Y, Liu J. J Wuhan Univ (Nat.Sci.Ed., 2012, 58(5: 463-470 24 Li H, Ruan J, Durbin R. Mapping short DNA sequencing reads and calling variants using mapping quality scores. Genome Research, 2008, 18(11: 1851-1858 25 Lin H, Zhang Z, Zhang M Q, et al. ZOOM! Zillions of oligos mapped. Bioinformatics, 20

47、08, 24(21: 2431-2437 26 Smith A D, Xuan Z, Zhang M Q. Using quality scores and longer reads improves accuracy of Solexa read mapping. BMC Bioinformatics, 2008, 9(1: 128-136 27 Langmead B, Trapnell C, Pop M, et al. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. G

48、enome Biol, 2009, 10(3: R25.1-R25.10 28 Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 2009, 25 (14: 17541760 29 Li R, Yu C, Li Y, et al. SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics, 2009, 25(15: 1966-1967 30 Home

49、r N, Merriman B, Nelson S F. BFAST: an alignment tool for large scale genome resequencing. PloS One, 2009, 4(11: e7767 31 Rumble S M, Lacroute P, Dalca A V, et al. SHRiMP: accurate mapping of short color-space reads. PLoS Computational Biology, 2009, 5(5: e1000386 32 Ryan M C, Cleland J, Kim R, et a

50、l. SpliceSeq: a resource for analysis and visualization of RNA-Seq data on alternative splicing and its functional impacts. Bioinformatics, 2012, 28 (18: 23852387 33 Xu G, Deng N, Zhao Z, et al. SAMMate: a GUI tool for processing short read alignments in SAM/BAM format. Source Code for Biology and M

51、edicine, 2011, 6(1: 2-13 34 Wang K, Singh D, Zeng Z, et al. MapSplice: accurate mapping of RNA-seq reads for splice junction discovery. Nucleic Acids Research, 2010, 38(18: e178-e178 35 Au K F, Jiang H, Lin L, et al. Detection of splice junctions from paired-end RNA-seq data by SpliceMap. Nucleic Ac

52、ids Research, 2010, 38(14: 4570-4578 36 Kim D, Pertea G, Trapnell C, et al. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol, 2013, 14(4: R36 37 Dobin A, Davis C A, Schlesinger F, et al. STAR: ultrafast universal RNA-seq aligner. Bi

53、oinformatics, 2013, 29(1: 15-21 38 Li H, Handsaker B, Wysoker A, et al. The sequence alignment/map format and SAMtools. Bioinformatics, 2009, 25(16: 2078-2079 39 Quinlan A R, Hall I M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics, 2010, 26(6: 841-842 ·

54、898 · 生物化學(xué)與生物物理進(jìn)展 Prog. Biochem. Biophys. 2015; 42 (10 40 Robinson J T, Thorvaldsdóttir H, Winckler W, et al. Integrative genomics viewer. Nature Biotechnology, 2011, 29(1: 24-26 41 Valouev A, Johnson D S, Sundquist A, et al. Genome-wide analysis of transcription factor binding sites based

55、 on ChIP-Seq data. Nature Methods, 2008, 5(9: 829-834 42 Meng J, Cui X, Liu H, et al. Unveiling the dynamics in RNA epigenetic regulations/BIBM. 2013 IEEE International Conference on Bioinformatics and Biomedicine. Shanghai: BIBM, 2013: 139-144 43 Krishnamoorthy K, Thomson J. A more powerful test fo

56、r comparing two Poisson means. Journal of Statistical Planning and Inference, 2004, 119(1: 23-35 44 Przyborowski J, Wilenski H. Homogeneity of results in testing samples from Poisson series with an application to testing clover seed for dodder. Biometrika, 1940, 31(3-4: 313-323 45 Zhang Y, Liu T, Me

57、yer C A, et al. Model-based analysis of ChIP-Seq (MACS. Genome Biol, 2008, 9(9: R137 46 Kharchenko P V, Tolstorukov M Y, Park P J. Design and analysis of ChIP-seq experiments for DNA-binding proteins. Nature Biotechnology, 2008, 26(12: 1351-1359 47 Pan Q, Shai O, Lee L J, et al. Deep surveying of al

58、ternative splicing complexity in the human transcriptome by high-throughput sequencing. Nature Genetics, 2008, 40(12: 1413-1415 48 Xu H, Wei C L, Lin F, et al. An HMM approach to genome-wide identification of differential histone modification sites from ChIP-seq data. Bioinformatics, 2008, 24(20: 2344-2349 49 Zhang X, Zhao J, Hao J K, et al. Conditional mutual inclusive information enables accurate quantification of associations in gene regulatory networks. Nucleic Acids Researc