1. 每種方法鑒定的排名靠前區(qū)域的真發(fā)現(xiàn)率(TDR)。TDR定義為調(diào)整后p值排名前靠前區(qū)域中真DMR占比。
2. DMR檢測(cè)方法的受試者工作特征(ROC)曲線。
3. 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）
接下來研究樣本量對(duì)DMR calling精確度的影響，因?yàn)闃颖玖客ǔＪ菍?shí)驗(yàn)設(shè)計(jì)中的主要參數(shù)。本研究兩組的模擬樣本量分別為2、3、5、7和10，每種條件下2、3、7、10個(gè)樣品的TDR分別如圖3A-D所示。幾乎所有方法在靠前排名 (如前100或前200)calling區(qū)域獲得高TDR(>0.8)，且在排名靠后時(shí)顯示出精確性下降。具體來說，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ū)域截止線標(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)確性高度依賴于表達(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)稱值顯著性水，所有方法在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è)模擬，其中沒有(0%)候選區(qū)域存在差異甲基化。使用bh調(diào)整后的p值以0.05標(biāo)稱值為顯著性水平獲得每種方法鑒定的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)稱值。exomePeak2更為“保守”，小樣本下(N = 2、3)以低靈敏度獲得最佳FDR(圖2D、E)。DRME是最“自由”方法，與其高靈敏度和FDR相匹配（圖2D、E）。
P值有效性分析檢測(cè)了null條件下P值是否均勻分布在0和1之間，并在圖5中說明了每組使用三個(gè)重復(fù)的結(jié)果。在QQ圖（Quantile-quantile plot）中，TRESS和exomePeak2產(chǎn)生的p值與預(yù)期值最為一致(圖5A，位于或接近對(duì)角線參考線)。大多數(shù)方法生成自由p值(右下方區(qū)域)，而FETHMM在大多數(shù)區(qū)域過于保守(左上方區(qū)域)。由于小p值在DMR檢測(cè)中更具信息性，本研究還應(yīng)用-log10轉(zhuǎn)換，重點(diǎn)關(guān)注圖5B中小p值分布。TRESS、exomePeak2和RADAR表現(xiàn)最好，而其他方法提供的p值過小，表明I型錯(cuò)誤膨脹。且對(duì)樣本量不敏感（圖2E），其中TRESS、exomePeak2和RADAR產(chǎn)生的的FDR控制得最好，也最穩(wěn)定。
表2：在無真DMR的null假設(shè)下，八種方法檢測(cè)DMR的I型錯(cuò)誤（0.05標(biāo)稱值顯著性水平計(jì)算，并在20次模擬中取平均值）
 


圖5：分析null條件下模擬觀察到的p值有效性。

1. QQ圖（Quantile–quantile plot）將p值分布與null下的期望分布U(0,1)進(jìn)行比較。
2. QQ圖進(jìn)行−log10轉(zhuǎn)換，重點(diǎn)關(guān)注小P值。在無DMR的null假設(shè)下進(jìn)行20次模擬。樣本量N=3 /組。

 
運(yùn)行時(shí)長和內(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í)長如圖6A所示。隨著樣本量增加，所有方法都顯示出更長的運(yùn)行時(shí)間。與其他方法相比，TRESS和exomePeak2的運(yùn)行時(shí)長都更短，且隨著樣本量增加更為明顯。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í)長和內(nèi)存消耗比較。

1. 五種不同方法在不同樣本量下的運(yùn)行時(shí)長比較，以小時(shí)為單位。
2. 計(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”。在該研究中，人類HeLa細(xì)胞系有四種樣品類型：一種野生型（WT）樣品和三種處理過的樣品，這些處理對(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ì)到人類參考基因組hg18，使用帶有默認(rèn)參數(shù)的STAR標(biāo)準(zhǔn)流程。比對(duì)后的BAM文件作為所有五種方法的input進(jìn)行比較，主要分析WT和METTL3樣品之間的差異甲基化。DMR calling在FDR<0.05的顯著性水平上進(jìn)行。在過濾掉短（寬度<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)行注釋（圖7A）。結(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中長區(qū)域，RADAR具有雙峰分布（bimodal distribution），覆蓋中長和長兩個(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通過剪接體參與半胱氨酸型內(nèi)肽酶活性、組蛋白H2B保守的C-末端賴氨酸去泛素化和mRNA剪接。對(duì)exomePeak2和RADAR進(jìn)行GO(Gene Ontology)通路分析（圖7F），在exomePeak2的DMR中，前三個(gè)GO富集是“生長因子受體和第二信使的信號(hào)通路疾�。―iseases of signal transduction by growth factor receptors and second messengers）”、“TP53轉(zhuǎn)錄調(diào)控（Transcriptional regulation by TP53）”和“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ì)照的人類樣本，“RD3”包括兩個(gè)Ythdf1敲低和對(duì)照的野生型小鼠骨髓來源樹突狀細(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ù)集之間再次觀察到一致結(jié)果（圖7H）。

圖7：真實(shí)數(shù)據(jù)的差異m6A甲基化方法。

1. 條形圖顯示在已鑒定的DMR中各種基因組特征分布。TRESS、exomePeak、exomePeak2、MeTDiff、RADAR和MACS3采用相同的FDR 0.05截止值來calling顯著性。
2. 維恩圖顯示通過五種方法鑒定的DMR重疊。
3. 六種方法的峰寬分布密度圖（log scale）。
4. exomePeak和exomePeak2的1038共有區(qū)域的成對(duì)FDR值散點(diǎn)圖。
5. 共有DMR的peaks差異分析可視化的兩個(gè)例子。所有差異peaks分析均在野生型（WT）組和METTL3組之間。
6. exomePeak2的DMR基因的GO分析，顯示DMR數(shù)最多。
7. 三個(gè)真實(shí)數(shù)據(jù)集中DMR計(jì)數(shù)排序。
8. 三個(gè)真實(shí)數(shù)據(jù)集中3’UTR百分比排序。

• 新型表觀轉(zhuǎn)錄組測(cè)序技術(shù)能夠使用數(shù)據(jù)驅(qū)動(dòng)的方法評(píng)估RNA修飾。
• 差異表觀轉(zhuǎn)錄組分析需要對(duì)成對(duì)的input對(duì)照和IP樣本進(jìn)行適當(dāng)?shù)慕�模，以適應(yīng)技術(shù)和生物噪聲、peaks值檢測(cè)并解決小樣本量問題。
• 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í)間最快。

 
參考文獻(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.