全基因組DNA甲基化測(cè)序數(shù)據(jù)工作流程分析和性能評(píng)估之分析軟件比較

DNA甲基化與轉(zhuǎn)錄調(diào)控、基因組印記、干細(xì)胞分化、胚胎發(fā)育和炎癥等過(guò)程有關(guān)。DNA甲基化異常可能揭示疾病狀態(tài)，包括癌癥和神經(jīng)系統(tǒng)疾病。因此，人類(lèi)基因組中5-甲基胞嘧啶（5mC）分布和位點(diǎn)是一個(gè)重要的研究方向。全基因組重亞硫酸鹽測(cè)序(Whole-genome bisulfite sequencing, WGBS)是一種用于分析DNA甲基化的高通量方法，本文綜述了WGBS的關(guān)鍵步驟，總結(jié)了可用且最新的分析工具，比較了比對(duì)算法，并分享了數(shù)據(jù)處理的示例代碼，為科研人員的DNA甲基化研究提供參考。
 
標(biāo)題：Analysis and Performance Assessment of the Whole Genome Bisulfite Sequencing Data Workflow: Currently Available Tools and a Practical Guide to Advance DNA Methylation Studies
期刊：《small methods》
時(shí)間：2022年3月
影響因子：10.7
 
WGBS文庫(kù)制備
廣泛使用的亞硫酸鹽轉(zhuǎn)化文庫(kù)構(gòu)建方法：Accel-NGS Methyl-Seq（Accel）、TruSeq DNA Methylation（TruSeq）和SPlinted ligation adapter tagging（SPLAT）。
   
WGBS分析流程
在處理原始測(cè)序數(shù)據(jù)時(shí)，生物信息學(xué)分析流程的可重復(fù)性和一致性至關(guān)重要。由于WGBS數(shù)據(jù)通常很龐大，因此需要大量的計(jì)算資源、內(nèi)存和存儲(chǔ)空間。這就要求分析流程不僅要穩(wěn)定，還要高效地使用內(nèi)存和時(shí)間。
   
四個(gè)比對(duì)工具在甲基化calling后的比對(duì)質(zhì)量顯示，BWA-METH和gemBS有最高的唯一比對(duì)率和最少的未比對(duì)reads（圖3A）。在每個(gè)染色體的平均CpG覆蓋度≤×10時(shí)存在微小差異，在≥×20覆蓋度時(shí)，BWA-METH比其他方法有略好的覆蓋度（圖3B）。DNA片段兩端未甲基化的Cs數(shù)量對(duì)于從比對(duì)結(jié)果中獲得合格的DNA甲基化calling至關(guān)重要。使用M-bias圖對(duì)每個(gè)流程比對(duì)結(jié)果的影響規(guī)模非常相似，盡管在不同文庫(kù)之間存在較大差異（圖3C）。使用比對(duì)軟件在基因組上calling的CpGs一致性顯示，在每個(gè)注釋區(qū)域和轉(zhuǎn)錄起始位點(diǎn)(TSSs)周?chē)�的平均甲基化水平上顯示出可比的基因組富集分布（圖3D-E）。與此同時(shí)，甲基化extraction水平不受比對(duì)軟件選擇的影響，因?yàn)锽ismark與其他三種比對(duì)軟件相比，甲基化calling相關(guān)性只略降低（圖4）。
   
（4）比對(duì)質(zhì)量評(píng)估（Alignment quality control）
比對(duì)后的QC在WGBS中至關(guān)重要，WGBS的內(nèi)在混雜變量會(huì)使甲基化估計(jì)偏向于過(guò)高或過(guò)低的估計(jì)，主要通過(guò)M-bias圖來(lái)檢測(cè)。在亞硫酸鹽處理過(guò)程中可能會(huì)發(fā)生部分（不完全）甲基化，此時(shí)可以觀察到C和T的peaks值，這通常會(huì)導(dǎo)致檢測(cè)過(guò)高。高于98.5%的比率可以確保沒(méi)有偏差�？梢栽谒�有背景中（CpG、CHG和CHH）向甲基化樣本中添加帶有未甲基化C的Spike序列，然后計(jì)數(shù)未甲基化C和T數(shù)量，并計(jì)算添加序列的轉(zhuǎn)化率。

同時(shí)，由于估計(jì)過(guò)低導(dǎo)致的偏差會(huì)捕獲到假陰性的甲基化位點(diǎn)。例如，通過(guò)酶介導(dǎo)的碎片化雙鏈DNA末端修復(fù)會(huì)在片段兩端引入未甲基化的C，從而導(dǎo)致人為的甲基化水平低估。這在M-bias圖中反映為問(wèn)題兩端的平均甲基化水平急劇下降，這應(yīng)在提取甲基化之前予以丟棄。亞硫酸鹽介導(dǎo)的降解是WGBS中偏差的主要來(lái)源，因?yàn)榻到夥请S機(jī)，發(fā)生在未甲基化的胞嘧啶上，這些胞嘧啶從文庫(kù)中舍棄。這導(dǎo)致許多隨后的序列偏差和整體甲基化的高估。
 
（5）甲基化信息提�。∕ethylation extraction）
在經(jīng)過(guò)BitMapperBS、BWA-METH、Bismark和GemBS等比對(duì)軟件進(jìn)行比對(duì)后，推薦利用MethylDackel進(jìn)行甲基化信息提取。例如用MethylDackel對(duì)BitMapperBS比對(duì)后的甲基化信息提取。甲基化評(píng)估通過(guò)比較測(cè)序reads和參考基因組進(jìn)行，如果在參考基因組中某個(gè)位點(diǎn)顯示為C，在上述位點(diǎn)注意到C時(shí)，就分配100%的甲基化，當(dāng)指示為T(mén)時(shí)，則分配0%的甲基化。計(jì)算加權(quán)平均值，并在計(jì)算該位點(diǎn)的C和T數(shù)量后，將其指定為最終的甲基化水平。如10/10 Cs顯示完全胞嘧啶甲基化，6/10 Cs顯示部分甲基化（60%），0/10 Cs代表未甲基化胞嘧啶。在提取之前，對(duì)兩個(gè)鏈中每個(gè)位點(diǎn)的平均甲基化水平進(jìn)行M-bias分析，以識(shí)別提取reads時(shí)的基本技術(shù)偏差作為修剪參考。從理論上講，reads應(yīng)該是恒定的，但每對(duì)中的第一和第二reads通常在5'和3'端有偏倚。reads中的人為噪聲會(huì)在提取過(guò)程中引發(fā)錯(cuò)誤的甲基化calling。MethylDackel在頂部和底部線條上給出了修剪建議，這些建議可以作為后續(xù)提取的參數(shù)。MethylDackel通過(guò)比對(duì)得到的BAM文件生成bedGraph文件，記錄甲基化與未甲基化位點(diǎn)信息，這些可以用來(lái)進(jìn)行數(shù)據(jù)過(guò)濾和進(jìn)一步數(shù)據(jù)分析。
 
數(shù)據(jù)歸一化與統(tǒng)計(jì)分析
（1）CpG甲基化
文庫(kù)制備方法會(huì)顯著影響每個(gè)CpG位點(diǎn)的平均覆蓋度。與甲基化芯片數(shù)據(jù)不同，測(cè)序數(shù)據(jù)沒(méi)有標(biāo)準(zhǔn)化的歸一化方法。但數(shù)據(jù)歸一化對(duì)下游差異甲基化檢測(cè)至關(guān)重要。降采樣(Downsampling)通過(guò)減少reads序列數(shù)量，使其與相似序列數(shù)據(jù)相匹配，從而實(shí)現(xiàn)歸一化。比對(duì)產(chǎn)生的BAM文件和甲基化提取產(chǎn)生的bedGraph文件都可以降采樣。在比對(duì)階段降采樣可能在時(shí)間和內(nèi)存上要求很高，而在提取階段降采樣則需要較少的時(shí)間和內(nèi)存，同時(shí)保證了相似的甲基化calling數(shù)量、檢測(cè)到的CpG位點(diǎn)、reads數(shù)分布和平均覆蓋度的準(zhǔn)確性。因此，在進(jìn)行進(jìn)一步的數(shù)據(jù)比較之前，建議先對(duì)bedGraph文件進(jìn)行降采樣，特別是對(duì)于差異甲基化區(qū)域（DMRs）的檢測(cè)。
 
（2）DMR
DMR（差異甲基化區(qū)域）檢測(cè)是核心甲基化分析之一，涉及對(duì)多個(gè)樣本中的基因組區(qū)域進(jìn)行分析。最常見(jiàn)的應(yīng)用是在癌癥和正常樣本之間尋找可能作為生物標(biāo)志物或揭示疾病生物學(xué)的異常甲基化區(qū)域�；�于DMR統(tǒng)計(jì)分析方法因軟件而異，以下是一些主要方法。

BSmooth：使用局部似然平滑方法（local likelihood smoothing approach）來(lái)鑒定樣本特異性甲基化信息中的DMR。應(yīng)用Welch's t-test（Student's t-test變體）比較多個(gè)樣本。DMR是具有觀察到的P值高于預(yù)定義β值的CpG位點(diǎn)。然而預(yù)定義閾值可能導(dǎo)致II類(lèi)錯(cuò)誤（假陰性），從而影響結(jié)果。

BiSeq：通過(guò)納入錯(cuò)誤發(fā)現(xiàn)率和β二項(xiàng)分布模型來(lái)解決這一問(wèn)題，充分考慮到生物學(xué)重復(fù)。然后通過(guò)triangular kernel模型調(diào)整分層過(guò)程引起的顯著變化，計(jì)算目標(biāo)區(qū)域的統(tǒng)計(jì)顯著性。P值被歸一化、轉(zhuǎn)換為z分?jǐn)?shù)、平均值進(jìn)行比較。

MethylSig：類(lèi)似于BiSeq，應(yīng)用β二項(xiàng)分布模型來(lái)考慮reads覆蓋度和生物學(xué)意義。

Metilene：結(jié)合了二項(xiàng)分割和多變量Kolmogorov-Smirnov擬合優(yōu)度檢驗(yàn)(K-S 檢驗(yàn))。這種非參數(shù)方法使用逐步分割基因組區(qū)域，被穩(wěn)步地最小化到CpG數(shù)量少于預(yù)定義下限的區(qū)域，或在統(tǒng)計(jì)顯著性上沒(méi)有改善的區(qū)域。這種方法對(duì)累積樣本分布的差異性更敏感。

methylKit：對(duì)單樣本情況使用Fisher's精確檢驗(yàn)，對(duì)多重復(fù)樣本使用基于logistical回歸的統(tǒng)計(jì)方法計(jì)算組間差異。

Defiant：是一個(gè)獨(dú)立的程序，使用加權(quán)Welch擴(kuò)展鑒定DMR。對(duì)于只有一個(gè)重復(fù)的兩個(gè)樣本，使用Fisher's精確檢驗(yàn)，對(duì)于有多個(gè)重復(fù)的樣本，使用Welch's t-test，基于覆蓋度對(duì)無(wú)偏樣本方差進(jìn)行加權(quán)。

Benjamini-Hochberg：應(yīng)用于調(diào)整DMR鑒定中多重t檢驗(yàn)的P值。數(shù)據(jù)分布本質(zhì)上是二項(xiàng)的，因?yàn)榇蠖鄶?shù)甲基化分布要么是完全甲基化的，要么是完全未甲基化的，表明二項(xiàng)分布模型性能優(yōu)于其他模型。

由于reads覆蓋度變化和人口統(tǒng)計(jì)參數(shù)（如性別、年齡和種族）的共變異對(duì)DMR檢測(cè)有強(qiáng)相關(guān)，因此數(shù)據(jù)的歸一化轉(zhuǎn)換和協(xié)變量調(diào)整至關(guān)重要。
 
案例研究
WGBS 數(shù)據(jù)的計(jì)算分析具有挑戰(zhàn)性，包括分析 FASTQ 讀長(zhǎng)、甲基化估計(jì)、位點(diǎn)注釋、DMR 檢測(cè)和可視化。以下為WGBS數(shù)據(jù)分析的綜合案例研究。

（1）分析工具
（2）數(shù)據(jù)分析
① Pre-alignment
   
使用以下代碼檢測(cè) R1 （WGBS_R1.fastq） 和 R2 （WGBS_R2.fastq） 讀取的原始 FASTQ 質(zhì)量： 結(jié)果包括一個(gè).html格式的摘要和一個(gè)帶有質(zhì)量數(shù)字的.zip文件夾。MultiQC 將具有多個(gè) FastQC 結(jié)果的樣品合并到一份.html報(bào)告中。  
② reads修剪
③ 比對(duì)和甲基化calling
Bismark：
BitMapperBS:
DMR的注釋和分析
易小結(jié)：
DNA甲基化和其他表觀基因組分析的動(dòng)態(tài)標(biāo)記與不同人類(lèi)疾病的診斷和預(yù)后相關(guān)聯(lián)。對(duì)甲基化結(jié)果的深刻且低偏倚解釋是下游生物學(xué)機(jī)制研究的核心。本綜述討論了分析WGBS數(shù)據(jù)的計(jì)算方法，并介紹了使用現(xiàn)有工具從原始reads檢測(cè)DMR所需的基本QC分析步驟。此外，還提出了甲基化文庫(kù)制備和數(shù)據(jù)處理中固有的潛在混雜因素及其對(duì)策。希望潛在用戶能夠理解WGBS的基本概念，從而加速人類(lèi)疾病的發(fā)現(xiàn)。

參考文獻(xiàn)：
1、Gong T, Borgard H, Zhang Z, Chen S, Gao Z, Deng Y. Analysis and Performance Assessment of the Whole Genome Bisulfite Sequencing Data Workflow: Currently Available Tools and a Practical Guide to Advance DNA Methylation Studies. Small Methods. 2022 Mar;6(3):e2101251. doi: 10.1002/smtd.202101251. PubMed PMID: 35064762.