8種表觀(guān)轉(zhuǎn)錄組m6A MeRIP-seq的差異甲基化區(qū)域（DMR）分析軟件比較 | 生信專(zhuān)區(qū)
大家好，這里是專(zhuān)注表觀(guān)組學(xué)十余年，領(lǐng)跑多組學(xué)科研服務(wù)的易基因。
RNA甲基化是近年來(lái)研究基因表達(dá)調(diào)控轉(zhuǎn)錄后變化的重要研究領(lǐng)域，包括N6-甲基腺苷（m6A）在內(nèi)的各種類(lèi)型RNA甲基化參與人類(lèi)疾病發(fā)展。MeRIP-seq作為一種新興的在轉(zhuǎn)錄組范圍內(nèi)定量檢測(cè)m6A水平的測(cè)序技術(shù)，拓展了RNA表觀(guān)遺傳學(xué)研究的基礎(chǔ)和臨床應(yīng)用，且呈上升趨勢(shì)。RNA甲基化數(shù)據(jù)分析的基本問(wèn)題之一是通過(guò)對(duì)比病例和對(duì)照來(lái)鑒定差異甲基化區(qū)域（DMR）?，F(xiàn)有開(kāi)發(fā)了多種用于DMR檢測(cè)的分析方法，但缺乏對(duì)這些分析方法的綜合評(píng)估。
本文利用模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)，全面評(píng)估了DMR calling的所有8種現(xiàn)有方法：FET-HMM、exomePeak（2.16.0版）、MeTDiff（1.1.0版）、DRME、QNB（1.0版）、exomePeak2（1.9.1版）、RADAR（0.2.4版）和TRESS（1.4.0版）。
模擬分析采用Gamma–Poisson模型和logit線(xiàn)性框架，并調(diào)試適應(yīng)各種樣本量和DMR比例進(jìn)行基準(zhǔn)檢測(cè)。所有8種方法在低input水平區(qū)域中觀(guān)察到低靈敏度，但樣本量增加會(huì)大大提高靈敏度。TRESS和exomePeak2在檢測(cè)精確度、FDR（False DiscoveryRate）、I型錯(cuò)誤（type-I error）調(diào)控和運(yùn)行時(shí)間等指標(biāo)上表現(xiàn)最好，但受限于低靈敏度。DRME和exomePeak以FDR（假發(fā)現(xiàn)率）和I型錯(cuò)誤（type-I error）膨脹為代價(jià)獲得高靈敏度。三個(gè)真實(shí)數(shù)據(jù)集分析表明，這些方法在鑒定DMR長(zhǎng)度和唯一發(fā)現(xiàn)區(qū)域有不同偏好性。
圖1：MeRIP-seq實(shí)驗(yàn)和DMR檢測(cè)示意圖
MeRIP-seq從RNA樣本生成配對(duì)的IP數(shù)據(jù)和input對(duì)照數(shù)據(jù)。將測(cè)序reads比對(duì)至參考基因組，然后通過(guò)最近開(kāi)發(fā)的統(tǒng)計(jì)方法鑒定差異甲基化區(qū)域（DMR）。其核心統(tǒng)計(jì)模型和特征列于餅圖內(nèi)圓，下游對(duì)DMR基因進(jìn)行peak注釋、biomarker發(fā)現(xiàn)、通路（pathway）鑒定和基因本體（GO）分析。
表1：現(xiàn)有差異RNA甲基化分析方法，方法按時(shí)間順序排列。
TDR：真發(fā)現(xiàn)率(True Discovery Rate)，即在某個(gè)截止點(diǎn)前被鑒定區(qū)域中所占的真陽(yáng)性。
FDR：假發(fā)現(xiàn)率（False Discovery Rate）。
隨著MeRIP-seq在過(guò)去幾年的廣泛使用，已經(jīng)開(kāi)發(fā)了幾種檢測(cè)DMR的計(jì)算方法，所有已知8種方法的詳細(xì)信息見(jiàn)表1（按開(kāi)發(fā)時(shí)間順序排列）。作為第一個(gè)發(fā)布的工具，exomePeak在兩種實(shí)驗(yàn)條件下對(duì)input對(duì)照和IP樣本的歸一化reads數(shù)應(yīng)用Fisher精確檢驗(yàn)(Fisher’s exact test，F(xiàn)ET)，因?yàn)槭褂盟兄貜?fù)的總reads，它忽略了生物學(xué)重復(fù)之間的異質(zhì)性。隨后改進(jìn)exomePeak開(kāi)發(fā)了MeTDiff和FET-HMM方法，MeTDiff假設(shè)β二項(xiàng)分布，并通過(guò)似然比檢驗(yàn)(Likelihood Ratio Test，LRT)比較不同條件下的甲基化水平，然而MeTDiff沒(méi)有很好解決測(cè)序深度的技術(shù)差異；FET-HMM采用了FET的改進(jìn)版本，并使用FET的二值決策作為差異甲基化狀態(tài)的觀(guān)察，隨后在檢測(cè)到的甲基化區(qū)域內(nèi)的小箱子上擬合隱馬爾可夫模型(Hidden Markov Model，HMM)，以合并沿基因組的依賴(lài)性，然而FET-HMM在合并每組中的重復(fù)檢測(cè)，忽略了生物學(xué)重復(fù)之間的組內(nèi)變化。2016年開(kāi)發(fā)的DRME解決了這個(gè)問(wèn)題，尤其是在小樣本量的情況下，DRME假設(shè)IP和input對(duì)照計(jì)數(shù)數(shù)據(jù)均為負(fù)二項(xiàng)式模型，且僅使用input對(duì)照數(shù)據(jù)來(lái)預(yù)測(cè)背景基因表達(dá)，通過(guò)計(jì)算基于IP數(shù)據(jù)統(tǒng)計(jì)顯著性來(lái)檢測(cè)DMR。DRME的作者后來(lái)改進(jìn)了他們的模型并開(kāi)發(fā)了QNB，QNB也使用負(fù)二項(xiàng)式模型，與DRME不同的是，QNB在背景表達(dá)式的預(yù)測(cè)和檢測(cè)統(tǒng)計(jì)的計(jì)算中結(jié)合了input對(duì)照和IP數(shù)據(jù)，DRME和QNB的共同限制是，兩者都將ip內(nèi)和input內(nèi)的變化作為目標(biāo)變化，但在MeRIP-seq中信號(hào)為IP/input比，應(yīng)該嚴(yán)格建模該比方差。后來(lái)exomePeak的作者又提出了exomePeak2，與exomePeak相比，exomePeak2解釋了IP效率和GC含量偏差變化，當(dāng)存在多個(gè)重復(fù)時(shí)，exomePeak2 calling DESeq2通過(guò)將IP和input視為成對(duì)樣本來(lái)鑒定DMR。
上述6種DMR分析方法僅適用于兩組間的比較。在真實(shí)的生物學(xué)實(shí)驗(yàn)中，特別是在大型研究中，經(jīng)常出現(xiàn)混雜協(xié)變量（如年齡或性別），但在上述方法中無(wú)法正確解釋。為了解決這個(gè)問(wèn)題，最近提出了兩種方法：RADAR和TRESS，兩種方法都使用線(xiàn)性框架將甲基化水平與實(shí)驗(yàn)因子相關(guān)聯(lián)。RADAR采用泊松隨機(jī)效應(yīng)模型（Poisson random effect model），而TRESS采用伽馬-泊松分布（Gamma–Poisson distribution）。TRESS與RADAR在兩個(gè)方面不同，第一個(gè)區(qū)別是TRESS假設(shè)原始reads數(shù)遵循負(fù)二項(xiàng)式分布，通常用于建模各種測(cè)序數(shù)據(jù)類(lèi)型。而RADAR假設(shè)預(yù)處理（從文庫(kù)大小歸一化開(kāi)始，然后進(jìn)行input對(duì)照調(diào)整）的計(jì)數(shù)數(shù)據(jù)遵循泊松分布，預(yù)處理后的數(shù)據(jù)不再是計(jì)數(shù)格式，因此泊松假設(shè)模棱兩可。另一個(gè)區(qū)別在于一旦模型擬合，TRESS可以檢測(cè)所有包含因子或其中任何線(xiàn)性組合的影響；而使用RADAR檢測(cè)不同因子，需要重新提供設(shè)計(jì)矩陣并再次擬合模型，在計(jì)算上比較低效?？偟膩?lái)說(shuō)，上面描述的方法列在表1中，顯示了input數(shù)據(jù)類(lèi)型、算法簡(jiǎn)要描述和使用中的各自?xún)?yōu)缺點(diǎn)。
數(shù)據(jù)生成模型和模擬（DATA GENERATIVE MODEL AND SIMULATION）
模擬框架的核心是伽瑪-泊松分布，并為適應(yīng)MeRIP-seq數(shù)據(jù)進(jìn)行了適當(dāng)修改。假設(shè)總共有10000個(gè)候選DMR，其中10%在處理和未處理的條件下差異甲基化。每種條件下的重復(fù)次數(shù)從2到10不等，以評(píng)估樣本量對(duì)DMR calling的影響。
整體差異peaks比較（OVERALL DIFFERENTIAL PEAK COMPARISON）
對(duì)8種m6A DMR檢測(cè)方法進(jìn)行基準(zhǔn)檢測(cè)，每種場(chǎng)景下進(jìn)行20次模擬。在一個(gè)共同實(shí)驗(yàn)設(shè)計(jì)下評(píng)估所有方法，作為整體比較基線(xiàn)，分別研究每個(gè)影響因子。
在每次模擬中，10000個(gè)候選DMR中的1000個(gè)被設(shè)計(jì)為真DMR，在病例組和對(duì)照組中都有三個(gè)重復(fù)。使用幾個(gè)評(píng)估指標(biāo)來(lái)評(píng)估八種DMR檢測(cè)算法的性能，如使用真發(fā)現(xiàn)率（TDR）來(lái)分析生物標(biāo)志物發(fā)現(xiàn)的精確度。還研究了ROC曲線(xiàn)（Receiver Operating Characteristic curve）、靈敏度和假發(fā)現(xiàn)率等經(jīng)典指標(biāo)。值得注意的是，F(xiàn)ETHMM共有三種策略：“FHB”、“FHC”和“FastFHC”。本文使用FHC而非默認(rèn)設(shè)置（FastFHC），因?yàn)槟J(rèn)設(shè)置中的編碼有對(duì)P值的異常操作。
圖2顯示了在上述基線(xiàn)模擬場(chǎng)景中所有方法之間的DMR calling性能比較。如圖2A所示，所有方法的TDR中，每個(gè)截止點(diǎn)，TRESS和exomePeak2的TDR值是最高的，且?guī)缀跸嗤?/strong>。exomePeak、FETHMM、DRME和QNB生成的TDR相當(dāng)，但略低于TRESS和exomePeak2。圖2B顯示TRESS、exomePeak2、exomePeak、FETHMM和DRME的AUC最高，而MeTDiff的AUC最低。值得注意的是，盡管TDR和ROC在方法比較中都是有效的指標(biāo)，但TDR更具信息性，因?yàn)榕琶壳暗闹笜?biāo)與生物標(biāo)志物發(fā)現(xiàn)更相關(guān)。圖2C中TRESS、exomePeak2、exomePeak、FETHMM和DRME的P值高度相關(guān)（Spearman相關(guān)性≥0.93），而MeTDiff得到的P值與其他方法差異較大。分析了最高和最低相關(guān)性的四對(duì)方法，并證明了FETHMM、exomePeak和exomePeak2之間的相似性。在圖2D、E中，Benjamini–Hochberg調(diào)整后的P值計(jì)算靈敏度和假發(fā)現(xiàn)率，以0.05為截止值。盡管DRME、FETHMM和exomePeak具有高靈敏度，但其FDR值也很高。結(jié)合這兩個(gè)指標(biāo)，表明DRME、FETHMM和exomePeak的I型錯(cuò)誤（type-I error）膨脹。MeTDiff和QNB在模擬中表現(xiàn)出不穩(wěn)定性。MeTDiff很難鑒定出真陽(yáng)性，因此并不是在所有的比較中都表現(xiàn)良好。TRESS和exomePeak2實(shí)現(xiàn)了幾乎相同的最佳整體性能，在保持低FDR的同時(shí)發(fā)現(xiàn)了相當(dāng)一部分真陽(yáng)性。圖2F中還總結(jié)了使用平均靈敏度和FDR的聯(lián)合分布的總體性能。理想情況下，好的方法應(yīng)該具有高靈敏度，同時(shí)保持低FDR，因此位于左上角區(qū)域的方法是首選方法。
圖2：m6A-DMR檢測(cè)方法的性能比較。
每種方法鑒定的排名靠前區(qū)域的真發(fā)現(xiàn)率(TDR)。TDR定義為調(diào)整后p值排名前靠前區(qū)域中真DMR占比。
DMR檢測(cè)方法的受試者工作特征(ROC)曲線(xiàn)。
8種方法的p值相關(guān)性熱圖。
D-E. 每種方法的靈敏度和FDR分布的小提琴圖，用BH調(diào)整后的p值計(jì)算。
F. 每種方法檢測(cè)DMR的平均靈敏度與FDR。模擬在三個(gè)病例組、三個(gè)對(duì)照組、10%真DMR的情況下進(jìn)行。N=20次模擬。
樣本量（SAMPLE SIZE）
接下來(lái)研究樣本量對(duì)DMR calling精確度的影響，因?yàn)闃颖玖客ǔＪ菍?shí)驗(yàn)設(shè)計(jì)中的主要參數(shù)。本研究?jī)山M的模擬樣本量分別為2、3、5、7和10，每種條件下2、3、7、10個(gè)樣品的TDR分別如圖3A-D所示。幾乎所有方法在靠前排名 (如前100或前200)calling區(qū)域獲得高TDR(>0.8)，且在排名靠后時(shí)顯示出精確性下降。具體來(lái)說(shuō)，TRESS和exomePeak2在所有截止點(diǎn)上都保持最高的精確度，而MeTDiff表現(xiàn)最差，沿秩遞減的精確度最低。隨著樣本量增加，所有方法的精確度都有所提高。當(dāng)N=7和10時(shí)，這種趨勢(shì)尤其明顯，其中方法報(bào)告的TDR值相似。在圖3E中，TDR以熱圖的形式呈現(xiàn)，包括所有模擬場(chǎng)景下的結(jié)果(N=2、3、5、7、10)，按排名前400、700、1000和1500區(qū)域進(jìn)行分層?？傮w而言，所有方法中TDR值隨區(qū)域排名提高和樣本量增加而增加。大樣本量可以大大提高檢測(cè)精確度，即使是排名中等區(qū)域(如前1000名)。RADAR和MeTDiff在小樣本量中 (N=2和3)的檢測(cè)精確度較低，但隨著樣本量增加，其性能幾乎相同。即使在極小的樣本量下(N= 2)， TRESS和exomePeak2的TDR也大于0.8。在經(jīng)驗(yàn)貝葉斯框架（empirical Bayes framework）下，TRESS和exomePeak2在全基因組中實(shí)現(xiàn)了信息借用，因此其在小樣本量中的表現(xiàn)優(yōu)于其他方法。在其他基因組學(xué)研究中，這種建模技術(shù)已被證明是有效的統(tǒng)計(jì)框架，特別是對(duì)于小樣本量。總之，對(duì)于小樣本量的項(xiàng)目，TRESS和exomePeak2是首選。
圖3：不同樣本量的DMR檢測(cè)精確度比較。
A-D. 每組中進(jìn)行2、3、7和10次重復(fù)的樣本量下，每種方法鑒定的排名靠前區(qū)域的真發(fā)現(xiàn)率（TDR）。
E. 不同樣本量和TOP區(qū)域截止值組合下的TDR值熱圖。樣本量標(biāo)注在右側(cè)，每組2個(gè)、3個(gè)、5個(gè)、7個(gè)和10個(gè)。排名靠前區(qū)域截止線(xiàn)標(biāo)注在左側(cè)，范圍從前400名、前700名、前1000名到前1500名。方法在熱圖中按列排序。在10%真DMR下進(jìn)行了N=20次模擬，取平均TDR值。
分層評(píng)估（STRATIFIED ASSESSMENT）
高通量測(cè)序數(shù)據(jù)（如批量RNA-seq）的差異表達(dá)分析準(zhǔn)確性高度依賴(lài)于表達(dá)水平，因此本研究按input范圍分層檢測(cè)DMR準(zhǔn)確性。根據(jù)input對(duì)照分布，候選區(qū)域根據(jù)其平均input計(jì)數(shù)分為五層：第一層1（0~10）、第2層（10~20）、第3層（20~40）、第4層（40~80）和第5層（80~+∞). 以0.05值為標(biāo)稱(chēng)值顯著性水，所有方法在5個(gè)分層中的靈敏度和FDR如圖4所示。方法按各層的平均值排序，當(dāng)從較低分層轉(zhuǎn)到較高分層時(shí)，所有方法都提高了靈敏度（圖4A–C），低input區(qū)域通常容易受模擬噪聲影響。即使在第一層，DRME也具有較高靈敏度(?>0.75），且在所有區(qū)域中具有相對(duì)較好性能。隨著樣本量增加，DRME靈敏度仍在提高。隨著樣本量增加，所有方法都表現(xiàn)出增加和減少的可變靈敏度，且這種性能增益對(duì)于較低層非常顯著，表明大樣本量有助于更可靠預(yù)測(cè)，尤其是對(duì)于受高背景噪聲影響更大的區(qū)域。其中，exomePeak2從樣本量的增加中受益最大，從第七位上升到第四位。對(duì)于FDR的結(jié)果，更大樣本量不如靈敏度（圖4D–F）。TRESS和exomePeak2在所有分層和樣本量大小中顯示出較小且最一致的假發(fā)現(xiàn)率（FDR）。在小樣本量下（N=3），MeTDiff在較低input區(qū)域的FDR較差，而隨著樣本量增加，F(xiàn)DR得到很大的提升。exomePeak、FETHMM和DRME受較差FDR影響，即使在大樣本情況下也是如此（N=10）。
圖4：按平均input計(jì)數(shù)值分層分析靈敏度和FDR。靈敏度和FDR以BH調(diào)整后P值計(jì)算，以0.05為截止值確定顯著性。
A–C. 分層靈敏度，每組分別設(shè)置3個(gè)、5個(gè)和10個(gè)重復(fù)。
D–F. 分層FDR，每組分別設(shè)置3個(gè)、5個(gè)和10個(gè)重復(fù)。在10%DMR下進(jìn)行N=20模擬。
I型錯(cuò)誤和p值有效性（TYPE I ERROR AND VALIDITY OF P-VALUES）
為了研究8種方法的I型錯(cuò)誤和p值有效性，在null條件下進(jìn)行假設(shè)模擬，其中沒(méi)有(0%)候選區(qū)域存在差異甲基化。使用bh調(diào)整后的p值以0.05標(biāo)稱(chēng)值為顯著性水平獲得每種方法鑒定的DMR。在每組設(shè)置2、3、5、7和10個(gè)重復(fù)的情況下，計(jì)算經(jīng)驗(yàn)I型錯(cuò)誤率（表2）。在所有情況下，TRESS和FETHMM的I型錯(cuò)誤率都接近0.05，表明其I型錯(cuò)誤率接近標(biāo)稱(chēng)值。exomePeak2更為“保守”，小樣本下(N = 2、3)以低靈敏度獲得最佳FDR(圖2D、E)。DRME是最“自由”方法，與其高靈敏度和FDR相匹配（圖2D、E）。
P值有效性分析檢測(cè)了null條件下P值是否均勻分布在0和1之間，并在圖5中說(shuō)明了每組使用三個(gè)重復(fù)的結(jié)果。在QQ圖（Quantile-quantile plot）中，TRESS和exomePeak2產(chǎn)生的p值與預(yù)期值最為一致(圖5A，位于或接近對(duì)角線(xiàn)參考線(xiàn))。大多數(shù)方法生成自由p值(右下方區(qū)域)，而FETHMM在大多數(shù)區(qū)域過(guò)于保守(左上方區(qū)域)。由于小p值在DMR檢測(cè)中更具信息性，本研究還應(yīng)用-log10轉(zhuǎn)換，重點(diǎn)關(guān)注圖5B中小p值分布。TRESS、exomePeak2和RADAR表現(xiàn)最好，而其他方法提供的p值過(guò)小，表明I型錯(cuò)誤膨脹。且對(duì)樣本量不敏感（圖2E），其中TRESS、exomePeak2和RADAR產(chǎn)生的的FDR控制得最好，也最穩(wěn)定。
表2：在無(wú)真DMR的null假設(shè)下，八種方法檢測(cè)DMR的I型錯(cuò)誤（0.05標(biāo)稱(chēng)值顯著性水平計(jì)算，并在20次模擬中取平均值）
圖5：分析null條件下模擬觀(guān)察到的p值有效性。
QQ圖（Quantile–quantile plot）將p值分布與null下的期望分布U(0,1)進(jìn)行比較。
QQ圖進(jìn)行?log10轉(zhuǎn)換，重點(diǎn)關(guān)注小P值。在無(wú)DMR的null假設(shè)下進(jìn)行20次模擬。樣本量N=3 /組。
運(yùn)行時(shí)長(zhǎng)和內(nèi)存消耗（RUNTIME AND MEMORY CONSUMPTION）
BAM文件為默認(rèn)input評(píng)估每種方法的軟件運(yùn)行時(shí)間和計(jì)算內(nèi)存消耗?；谕还?jié)點(diǎn)、同一內(nèi)核和200 GB內(nèi)存的高性能計(jì)算（HPC），在不同樣本量下，五種方法的運(yùn)行時(shí)長(zhǎng)如圖6A所示。隨著樣本量增加，所有方法都顯示出更長(zhǎng)的運(yùn)行時(shí)間。與其他方法相比，TRESS和exomePeak2的運(yùn)行時(shí)長(zhǎng)都更短，且隨著樣本量增加更為明顯。exomePeak和MeTDiff在所有樣本量中具有相似的運(yùn)行時(shí)間。RADAR的運(yùn)行時(shí)間最慢。由于大多數(shù)方法都將BAM文件作為標(biāo)準(zhǔn)input，因此進(jìn)一步對(duì)計(jì)算內(nèi)存消耗進(jìn)行了基準(zhǔn)檢測(cè)（圖6B）。MeTDiff和exomePeak消耗內(nèi)存最少（分別為3.81 GB和4.62GB）。TRESS消耗的內(nèi)存略多于MeTDiff和exomePeak。exomePeak2利用了最多的內(nèi)存（170.28GB）。模擬在HPC中進(jìn)行，每個(gè)方法calling都有1個(gè)節(jié)點(diǎn)、40個(gè)內(nèi)核和200 GB可用內(nèi)存。
圖6：m6A DMR檢測(cè)方法的運(yùn)行時(shí)長(zhǎng)和內(nèi)存消耗比較。
五種不同方法在不同樣本量下的運(yùn)行時(shí)長(zhǎng)比較，以小時(shí)為單位。
計(jì)算五種不同方法的內(nèi)存消耗，單位為GB。
真實(shí)數(shù)據(jù)分析（REAL DATA ANALYSIS）
首先從一項(xiàng)研究METTL3-METTL14復(fù)合體介導(dǎo)哺乳動(dòng)物核RNA m6A甲基化的研究中獲得了真實(shí)數(shù)據(jù)集（GSE46705），將其標(biāo)記為“RD1”。在該研究中，人類(lèi)HeLa細(xì)胞系有四種樣品類(lèi)型：一種野生型（WT）樣品和三種處理過(guò)的樣品，這些處理對(duì)應(yīng)于復(fù)合體METTL3、METLL14和WTAP的敲除（KD）。每個(gè)樣品2個(gè)重復(fù)。將TRESS、exomePeak、exomePeak2、MeTDiff和RADAR方法應(yīng)用于該真實(shí)數(shù)據(jù)，以鑒定m6A差異甲基化。同時(shí)還采用了適用于分析MeRIP-seq數(shù)據(jù)的MACS3方法。MACS3已被先前的幾項(xiàng)研究表明其作為MeRIP-seq數(shù)據(jù)差異分析的有效工具的潛力。該研究只分析以BAM文件為input比較的方法，因此排除了QNB和DRME等以reads計(jì)數(shù)矩陣為input的方法。
原始FASTQ文件被比對(duì)到人類(lèi)參考基因組hg18，使用帶有默認(rèn)參數(shù)的STAR標(biāo)準(zhǔn)流程。比對(duì)后的BAM文件作為所有五種方法的input進(jìn)行比較，主要分析WT和METTL3樣品之間的差異甲基化。DMR calling在FDR<0.05的顯著性水平上進(jìn)行。在過(guò)濾掉短（寬度<150）和重疊區(qū)域后，TRESS、exomePeak、exomePeak2、MeTDiff和RADAR分別鑒定出1413、1397、5272、161和2924個(gè)DMR。exomePeak2鑒定出最多的DMR，而MeTDiff鑒定的DMR最少。
五種方法分析WT組與METTL3組真實(shí)數(shù)據(jù)的性能比較如圖7所示。使用ChIPseeker對(duì)DMR進(jìn)行注釋?zhuān)▓D7A）。結(jié)果顯示，除了RADAR以外的大多數(shù)方法都支持3'UTR的DMR。RADAR偏好基因下游外顯子區(qū)（即非第一外顯子）。所有方法的啟動(dòng)子和下游外顯子區(qū)均表現(xiàn)出相當(dāng)數(shù)量的組成基因組表征。圖7B顯示了五種方法的5個(gè)重疊區(qū)域。exomePeak2發(fā)現(xiàn)3348個(gè)特異性DMR，是所有DMR中最高的。兩種方法之間重疊區(qū)域的最高數(shù)量是由exomePeak和exomePeak2 calling的1038個(gè)重疊，而兩種方法間重疊區(qū)域的最少數(shù)量是由TRESS和MeTDiff calling的15個(gè)重疊。DMR的peaks寬分布（log scale）如圖7C所示。TRESS偏好150–400bp中長(zhǎng)區(qū)域，RADAR具有雙峰分布（bimodal distribution），覆蓋中長(zhǎng)和長(zhǎng)兩個(gè)區(qū)域。鑒定出1038個(gè)共有區(qū)域的FDR（圖7D）。與exomePeak2相比，exomePeak是一種更保守的方法。同時(shí)由TRESS、exomePeak、exomePeak2、MeTDiff和RADAR方法顯示了WT和METTL3樣品之間共有DMR的兩個(gè)示例（圖7E），這兩個(gè)區(qū)域覆蓋蛋白編碼基因TEX264（chr3）、PRICKLE4、TOMM6和USP49（chr6）。先前的研究表明，TEX264能夠激活信號(hào)受體活性，并參與蛋白-DNA共價(jià)交聯(lián)修復(fù)。USP46通過(guò)剪接體參與半胱氨酸型內(nèi)肽酶活性、組蛋白H2B保守的C-末端賴(lài)氨酸去泛素化和mRNA剪接。對(duì)exomePeak2和RADAR進(jìn)行GO(Gene Ontology)通路分析（圖7F），在exomePeak2的DMR中，前三個(gè)GO富集是“生長(zhǎng)因子受體和第二信使的信號(hào)通路疾?。―iseases of signal transduction by growth factor receptors and second messengers）”、“TP53轉(zhuǎn)錄調(diào)控（Transcriptional regulation by TP53）”和“I類(lèi)MHC介導(dǎo)的抗原處理和呈遞（Class I MHC mediated antigen processing & presentation）”。
同時(shí)在另外兩個(gè)真實(shí)數(shù)據(jù)集（GSE94613和GSE115105）中進(jìn)行檢測(cè)，并將它們標(biāo)記為“RD2”和“RD3”，其中， “RD2”包括12個(gè)METTL3敲低細(xì)胞系和對(duì)照的人類(lèi)樣本，“RD3”包括兩個(gè)Ythdf1敲低和對(duì)照的野生型小鼠骨髓來(lái)源樹(shù)突狀細(xì)胞（BMDC）。對(duì)這兩個(gè)數(shù)據(jù)集進(jìn)行相同的分析，根據(jù)DMR數(shù)量和3’UTR在三個(gè)真實(shí)數(shù)據(jù)集中的百分比對(duì)五種方法進(jìn)行排序（圖7G）。exomePeak2軟件calling了三個(gè)真實(shí)數(shù)據(jù)集中最多的DMR，其次是RADAR。在3’UTR方面，除了exomePeak2和MeTDiff之間的微小差異外，數(shù)據(jù)集之間再次觀(guān)察到一致結(jié)果（圖7H）。
圖7：真實(shí)數(shù)據(jù)的差異m6A甲基化方法。
條形圖顯示在已鑒定的DMR中各種基因組特征分布。TRESS、exomePeak、exomePeak2、MeTDiff、RADAR和MACS3采用相同的FDR 0.05截止值來(lái)calling顯著性。
維恩圖顯示通過(guò)五種方法鑒定的DMR重疊。
六種方法的峰寬分布密度圖（log scale）。
exomePeak和exomePeak2的1038共有區(qū)域的成對(duì)FDR值散點(diǎn)圖。
共有DMR的peaks差異分析可視化的兩個(gè)例子。所有差異peaks分析均在野生型（WT）組和METTL3組之間。
exomePeak2的DMR基因的GO分析，顯示DMR數(shù)最多。
三個(gè)真實(shí)數(shù)據(jù)集中DMR計(jì)數(shù)排序。
三個(gè)真實(shí)數(shù)據(jù)集中3’UTR百分比排序。
比較要點(diǎn)（Key Points）
新型表觀(guān)轉(zhuǎn)錄組測(cè)序技術(shù)能夠使用數(shù)據(jù)驅(qū)動(dòng)的方法評(píng)估RNA修飾。
差異表觀(guān)轉(zhuǎn)錄組分析需要對(duì)成對(duì)的input對(duì)照和IP樣本進(jìn)行適當(dāng)?shù)慕＃赃m應(yīng)技術(shù)和生物噪聲、peaks值檢測(cè)并解決小樣本量問(wèn)題。
TRESS和exomePeak2在基準(zhǔn)研究中表現(xiàn)出高TDR、低FDR和超高靈敏度。
檢測(cè)精確度可能會(huì)受低input表達(dá)影響，但受益于樣本量增加。
RADAR、TRESS和exomePeak2顯示了頂級(jí)嚴(yán)格的I型錯(cuò)誤控制和null下的有效p值分布。MeTDiff計(jì)算內(nèi)存消耗最少，TRESS運(yùn)行時(shí)間最快。
關(guān)于易基因RNA m6A甲基化測(cè)序(MeRIP-seq)技術(shù)
易基因MeRIP-seq技術(shù)利用m6A特異性抗體富集發(fā)生m6A修飾的RNA片段（包括mRNA、lncRNA等rRNA去除所有RNA），結(jié)合高通量測(cè)序，可以對(duì)RNA上的m6A修飾進(jìn)行定位與定量，總RNA起始量可降低至10μg，最低僅需1μg總RNA。廣泛應(yīng)用于組織發(fā)育、干細(xì)胞自我更新和分化、熱休克或DNA損傷應(yīng)答、癌癥發(fā)生與發(fā)展、藥物應(yīng)答等研究領(lǐng)域；可應(yīng)用于動(dòng)物、植物、細(xì)胞及組織的m6A檢測(cè)。
大樣本量m6A-QTL性狀關(guān)聯(lián)分析，傳統(tǒng)MeRIP單個(gè)樣品價(jià)格高，通常難以承擔(dān)。易基因開(kāi)發(fā)建立MeRIP-seq2技術(shù)，顯著提成IP平行性，實(shí)現(xiàn)不同樣本間相對(duì)定量，降低檢測(cè)成本。
易基因提供適用于不同科研需求的MeRIP技術(shù)：
m6A甲基化-常量mRNA 甲基化測(cè)序（MeRIP-seq）
m6A甲基化-常量mRNA +lncRNA甲基化測(cè)序（lnc-MeRIP-seq）
m6A甲基化-微量mRNA +lncRNA甲基化測(cè)序（Micro-lnc-MeRIP-seq）
高通量m6A甲基化-常量mRNA甲基化測(cè)序（MeRIP-seq2）
技術(shù)優(yōu)勢(shì)：
起始量低：樣本起始量可降低至10-20μg，最低僅需1μg總RNA；
轉(zhuǎn)錄組范圍內(nèi)：可以同時(shí)檢測(cè)mRNA和lncRNA；
樣本要求：可用于動(dòng)物、植物、細(xì)胞及組織的m6A檢測(cè)；
重復(fù)性高：IP富集重復(fù)性高，最大化降低抗體富集偏差；
應(yīng)用范圍廣：廣泛應(yīng)用于組織發(fā)育、干細(xì)胞自我更新和分化、熱休克或DNA損傷應(yīng)答、癌癥的發(fā)生與發(fā)展、藥物應(yīng)答等研究領(lǐng)域。
研究方向：
m6A甲基化目前主要運(yùn)用在分子機(jī)制的理論性研究
疾病發(fā)生發(fā)展：腫瘤、代謝疾?。ㄈ绶逝?糖尿病）、神經(jīng)和精神疾?。ㄈ绨柶澓ＤY/抑郁癥）、炎癥…
發(fā)育和分化：早期胚胎發(fā)育、個(gè)體/組織/器官生長(zhǎng)發(fā)育、干細(xì)胞分化與命運(yùn)決定、衰老
環(huán)境暴露與響應(yīng)：污染、抗逆、生活方式
關(guān)于m6A甲基化研究思路
（1）整體把握m6A甲基化圖譜特征：m6A peak數(shù)量變化、m6A修飾基因數(shù)量變化、單個(gè)基因m6A peak數(shù)量分析、m6A peak在基因元件上的分布、m6A peak的motif分析、m6A peak修飾基因的功能分析
（2）篩選具體差異m6A peak和基因：差異m6A peak鑒定、非時(shí)序數(shù)據(jù)的分析策略、時(shí)序數(shù)據(jù)的分析策略、差異m6A修飾基因的功能分析、差異m6A修飾基因的PPI分析、候選基因的m6A修飾可視化展示
（3）m6A甲基化組學(xué)&轉(zhuǎn)錄組學(xué)關(guān)聯(lián)分析：Meta genes整體關(guān)聯(lián)、DMG-DEG對(duì)應(yīng)關(guān)聯(lián)、m6A修飾目標(biāo)基因的篩選策略
（4）進(jìn)一步驗(yàn)證或后期試驗(yàn)
易基因科技提供全面的RNA甲基化研究整體解決方案，詳詢(xún)易基因0755-28317900。
參考文獻(xiàn)：
Duan D, Tang W, Wang R, Guo Z, Feng H. Evaluation of epitranscriptome-wide N6-methyladenosine differential analysis methods. Brief Bioinform. 2023 May 19;24(3) pii: 7111718.
