DNA甲基化與轉(zhuǎn)錄調(diào)控荐绝、基因組印記、干細胞分化避消、胚胎發(fā)育和炎癥等過程有關(guān)低滩。DNA甲基化異痴偌校可能揭示疾病狀態(tài)，包括癌癥和神經(jīng)系統(tǒng)疾病恕沫。因此监憎，人類基因組中5-甲基胞嘧啶（5mC）分布和位點是一個重要的研究方向。全基因組重亞硫酸鹽測序(Whole-genome bisulfite sequencing, WGBS)是一種用于分析DNA甲基化的高通量方法婶溯，本文綜述了WGBS的關(guān)鍵步驟鲸阔，總結(jié)了可用且最新的分析工具，比較了比對算法迄委，并分享了數(shù)據(jù)處理的示例代碼褐筛，為科研人員的DNA甲基化研究提供參考。
標題：Analysis and Performance Assessment of the Whole Genome BisulfiteSequencing Data Workflow: Currently Available Tools and a Practical Guide toAdvance DNA Methylation Studies
期刊：《small methods》
時間：2022年3月
影響因子：10.7
WGBS文庫制備
廣泛使用的亞硫酸鹽轉(zhuǎn)化文庫構(gòu)建方法：Accel-NGS Methyl-Seq（Accel）叙身、TruSeq DNA Methylation（TruSeq）和SPlinted ligation adapter tagging（SPLAT）渔扎。
WGBS分析流程
在處理原始測序數(shù)據(jù)時曲梗，生物信息學(xué)分析流程的可重復(fù)性和一致性至關(guān)重要赞警。由于WGBS數(shù)據(jù)通常很龐大，因此需要大量的計算資源虏两、內(nèi)存和存儲空間愧旦。這就要求分析流程不僅要穩(wěn)定，還要高效地使用內(nèi)存和時間定罢。
（1）質(zhì)量評估（Quality Assessment）
預(yù)比對質(zhì)量評估確保了原始數(shù)據(jù)的正確輸入并提高了可比對性与柑。多維度評估審查了低質(zhì)量reads、接頭序列蓄坏、污染序列和重復(fù)reads价捧。合成測序（增加測序的堿基數(shù)量）會降低質(zhì)量，尤其是Illumina平臺涡戳。在測序過程中结蟋，DNA簇中的單個分子可能無法成功合成。錯誤合成的分子聚類會導(dǎo)致motif calling錯誤渔彰，這些錯誤與片段長強正相關(guān)嵌屎。文庫構(gòu)建中長片段比例越高，可能會導(dǎo)致更多的calling錯誤恍涂、較低的Phred得分和更高的不匹配率宝惰，特別是對于成對末端對齊。建議保留質(zhì)量分數(shù)≥30的堿基再沧，這表示對motif calling準確性99.9%尼夺。DNA測序反應(yīng)讀長通常比DNA片段長，片段兩端的接頭序列可能被錯誤地測序炒瘸，這會引入構(gòu)成性的甲基化C淤堵，并在隨后的甲基化calling中引起偏差，因此建議提前檢查并修剪接頭顷扩。數(shù)據(jù)可能被污染拐邪，包括在文庫制備步驟中添加引物和載體序列，以及為校準測序反應(yīng)中的序列質(zhì)量而添加的Phi X噬菌體DNA隘截。過度表示的污染物通常呈現(xiàn)顯著高于核基因組的reads深度扎阶，導(dǎo)致在比對過程中失敗。當同一DNA片段被雙重計數(shù)時技俐，可能會發(fā)生重復(fù)reads（兩個具有相同基因位點的reads）乘陪。
長插入片段可能會從更高比例的高質(zhì)量reads、較少接頭污染和更高有效reads深度（增加基因組覆蓋效率）中受益雕擂。TruSeq文庫制備需要特別注意重復(fù)reads啡邑，其PCR擴增cycles數(shù)量大約是Accel的兩倍（10-12 cycles對比6-9 cycles）。
建議使用FastQC對原始FASTQ文件在修剪前后進行質(zhì)量評估井赌。QC結(jié)果包括一個.HTML格式的摘要和包含質(zhì)量圖表的.zip文件夾谤逼。檢測每個樣本的總reads數(shù)贵扰、每個堿基的序列質(zhì)量和接頭檢測圖，以確保良好的質(zhì)量比對和比對reads數(shù)量流部。在亞硫酸鹽轉(zhuǎn)化過程中觀察到GC比率和內(nèi)容分布的偏倚戚绕，修剪后的FastQC報告作為雙重保障，確保在修剪后（去除接頭）每個樣本的總reads數(shù)保持不變枝冀，并確保隨機引物引起的堿基序列變化舞丛。
（2）reads修剪（Trimming）
使用FastQC對原始FASTQ文件進行質(zhì)量評估后，可以利用Trim Galore軟件（具有簡化的命令行和參數(shù)果漾，能夠自動識別技術(shù)序列）和Trimmomatic軟件對低質(zhì)量reads球切、接頭序列、污染序列绒障、重復(fù)序列進行去除吨凑。這一步驟稱之為修剪(Trimming)。WGBS數(shù)據(jù)修剪主要有兩個方面：質(zhì)量修剪和接頭去除户辱。質(zhì)量修剪側(cè)重于原始reads端的質(zhì)量下降鸵钝，以及序列組成的頭部和尾部偏差。修剪方法包括剪除低于特定質(zhì)量閾值的區(qū)域（Phred默認分數(shù)為20庐镐，表示1/100的堿基可能是錯誤的）恩商，或者從reads開始和結(jié)束修剪自定數(shù)量堿基。對于非亞硫酸鹽測序的接頭去除焚鹊，依賴于“過度表示的序列”和“每條序列的GC含量”痕届。由于CG含量在WGBS數(shù)據(jù)中不適用，過度表示的序列作為接頭污染的指標末患，在Trim Galore中通常被檢測到研叫。修剪參數(shù)的其他考慮因素是文庫方向和輸入reads的成對/單端特性。Trim Galore對文庫的默認設(shè)置具有方向性璧针，在Trim Galore中應(yīng)指定成對末端的特性嚷炉，以保持成對reads，以便進一步比對探橱。兩個文件的reads數(shù)不匹配和reads名稱不一致會觸發(fā)警告申屹。
（3）比對（Alignment）
WGBS的原理是對未甲基化的胞嘧啶(Cs)通過亞硫酸鹽處理轉(zhuǎn)化為胸腺嘧啶(Ts)，同時保留甲基化的Cs隧膏。理想情況下哗讥，當reads序列比對到參考基因組時，可以識別未甲基化的Cs胞枕。然而亞硫酸鹽轉(zhuǎn)化帶來數(shù)據(jù)比對兩大計算挑戰(zhàn)：C-T比對不匹配和序列復(fù)雜性降低杆煞。C-T比對不匹配指的是在測序reads中，T可能與參考基因組中的C比對，反之則不然决乎。序列復(fù)雜性降低使得難以區(qū)分轉(zhuǎn)化后的Ts與系統(tǒng)錯誤队询，夸大了比對不準確。
為了解決這種變化帶來的reads與參考基因組不匹配的問題构诚，有兩種主要的比對策略：三字符策略和通配符策略蚌斩。三字符策略通過將參考基因組和序列reads中的所有Cs轉(zhuǎn)換為Ts，將基于四字符的基因組簡化為基于三字符的基因組范嘱。之后使用標準的比對工具處理reads序列送膳，如Bowtie1或Bowtie2、BWA mem和GEM3彤侍，主要采用Burrows-Wheeler變換(BWT)回溯算法肠缨。而通配符策略用基因組中的Cs轉(zhuǎn)換為Ys逆趋，這可以與reads序列中的Cs和Ts進行匹配盏阶。在選擇比對工具時垮抗，比對的準確性和計算時間是主要的考慮因素负甸。Bock(2012)建議，通配符比對策略實現(xiàn)了更高的基因組覆蓋率枣察，但增加了高甲基化水平評估偏差的可能魄眉，而三字符策略則相反砰盐。通配符比對軟件在BS轉(zhuǎn)化reads中保留Cs，并將序列復(fù)雜性提高到確保與參考基因組唯一比對水平坑律，而三字符比對軟件在BS轉(zhuǎn)化reads中去除Cs岩梳，降低序列復(fù)雜性，增加了模糊比對位點的機會晃择。但覆蓋差異和M偏差僅在基因組中的高相似區(qū)域中表現(xiàn)出來冀值，因此在如人類基因組的長序列reads比對不太相關(guān)。因此在比對軟件選擇中優(yōu)先考慮計算速度和內(nèi)存消耗宫屠。研究表明如Bismark和BWA-METH等三字符比對軟件在運行時間和峰值內(nèi)存使用方面優(yōu)于如BRAT_BW列疗、BSMAP和GSnap的通配符比對軟件。
在眾多比對軟件中浪蹂，BitMapperBS和FMtree相對更節(jié)省時間抵栈，但與Bismark、BWA-METH和gemBS比對軟件相比坤次，在1,000,000 bp成對末端模擬數(shù)據(jù)read上古劲，并沒有觀察到6-7倍的計算時間減少。對于處理大型哺乳動物測序項目的人來說缰猴，從大約24小時減少到7-8小時可能是有說服力的产艾，盡管在追求速度時可能會犧牲比對質(zhì)量。BitMapperBS可能無法保證在有多個不匹配的情況下獲得最高質(zhì)量比對，因此作者更為推薦Bismark胰舆、BWA-METH骚露、gemBS，這幾款軟件能節(jié)省約1/3的運行時間且能良好的平衡運行時間與比對質(zhì)量之間的關(guān)系缚窿。
用于亞硫酸鹽處理reads比對算法在運行時間和比對質(zhì)量上的差異會影響下游的甲基化calling棘幸。在Accel-NGS MethylSeq、SPLAT和TruSeq等WGBS框架以及TrueMethyl和EMSeq的氧化亞硫酸鹽測序中倦零，對廣泛應(yīng)用的比對軟件Bismark误续、BitMapperBS、BWA-METH和gemBS進行評估扫茅。通過使用Seqtk的相同數(shù)據(jù)檢測的五種文庫制備方法蹋嵌，從樣本數(shù)據(jù)中減去1,000,000 bp成對末端reads，比較其運行速度葫隙。比較結(jié)果顯示栽烂，BitMapperBS具有最高的對齊速度，平均每秒約550-650 reads（表1）恋脚。Bismark腺办、BWA-METH和gemBS顯示出相同的比對速度（約每秒200-300reads0；而Bismark最不穩(wěn)定糟描。
四個比對工具在甲基化calling后的比對質(zhì)量顯示，BWA-METH和gemBS有最高的唯一比對率和最少的未比對reads（圖3A）船响。在每個染色體的平均CpG覆蓋度≤×10時存在微小差異躬拢，在≥×20覆蓋度時，BWA-METH比其他方法有略好的覆蓋度（圖3B）见间。DNA片段兩端未甲基化的Cs數(shù)量對于從比對結(jié)果中獲得合格的DNA甲基化calling至關(guān)重要聊闯。使用M-bias圖對每個流程比對結(jié)果的影響規(guī)模非常相似，盡管在不同文庫之間存在較大差異（圖3C）缤剧。使用比對軟件在基因組上calling的CpGs一致性顯示馅袁，在每個注釋區(qū)域和轉(zhuǎn)錄起始位點(TSSs)周圍的平均甲基化水平上顯示出可比的基因組富集分布（圖3D-E）。與此同時荒辕，甲基化extraction水平不受比對軟件選擇的影響汗销，因為Bismark與其他三種比對軟件相比，甲基化calling相關(guān)性只略降低（圖4）抵窒。
（4）比對質(zhì)量評估（Alignment quality control）
比對后的QC在WGBS中至關(guān)重要，WGBS的內(nèi)在混雜變量會使甲基化估計偏向于過高或過低的估計李皇，主要通過M-bias圖來檢測削茁。在亞硫酸鹽處理過程中可能會發(fā)生部分（不完全）甲基化宙枷，此時可以觀察到C和T的peaks值，這通常會導(dǎo)致檢測過高茧跋。高于98.5%的比率可以確保沒有偏差慰丛。可以在所有背景中（CpG瘾杭、CHG和CHH）向甲基化樣本中添加帶有未甲基化C的Spike序列诅病，然后計數(shù)未甲基化C和T數(shù)量，并計算添加序列的轉(zhuǎn)化率粥烁。
同時贤笆，由于估計過低導(dǎo)致的偏差會捕獲到假陰性的甲基化位點。例如讨阻，通過酶介導(dǎo)的碎片化雙鏈DNA末端修復(fù)會在片段兩端引入未甲基化的C芥永，從而導(dǎo)致人為的甲基化水平低估。這在M-bias圖中反映為問題兩端的平均甲基化水平急劇下降钝吮，這應(yīng)在提取甲基化之前予以丟棄埋涧。亞硫酸鹽介導(dǎo)的降解是WGBS中偏差的主要來源，因為降解非隨機搀绣，發(fā)生在未甲基化的胞嘧啶上飞袋，這些胞嘧啶從文庫中舍棄。這導(dǎo)致許多隨后的序列偏差和整體甲基化的高估链患。
（5）甲基化信息提取（Methylation extraction）
在經(jīng)過BitMapperBS瓶您、BWA-METH麻捻、Bismark和GemBS等比對軟件進行比對后，推薦利用MethylDackel進行甲基化信息提取呀袱。例如用MethylDackel對BitMapperBS比對后的甲基化信息提取贸毕。甲基化評估通過比較測序reads和參考基因組進行，如果在參考基因組中某個位點顯示為C夜赵，在上述位點注意到C時明棍，就分配100%的甲基化，當指示為T時寇僧，則分配0%的甲基化摊腋。計算加權(quán)平均值，并在計算該位點的C和T數(shù)量后嘁傀，將其指定為最終的甲基化水平兴蒸。如10/10 Cs顯示完全胞嘧啶甲基化，6/10 Cs顯示部分甲基化（60%）细办，0/10
Cs代表未甲基化胞嘧啶橙凳。在提取之前，對兩個鏈中每個位點的平均甲基化水平進行M-bias分析，以識別提取reads時的基本技術(shù)偏差作為修剪參考岛啸。從理論上講钓觉，reads應(yīng)該是恒定的，但每對中的第一和第二reads通常在5'和3'端有偏倚坚踩。reads中的人為噪聲會在提取過程中引發(fā)錯誤的甲基化calling议谷。MethylDackel在頂部和底部線條上給出了修剪建議，這些建議可以作為后續(xù)提取的參數(shù)堕虹。MethylDackel通過比對得到的BAM文件生成bedGraph文件卧晓，記錄甲基化與未甲基化位點信息，這些可以用來進行數(shù)據(jù)過濾和進一步數(shù)據(jù)分析赴捞。
數(shù)據(jù)歸一化與統(tǒng)計分析
（1）CpG甲基化
文庫制備方法會顯著影響每個CpG位點的平均覆蓋度逼裆。與甲基化芯片數(shù)據(jù)不同，測序數(shù)據(jù)沒有標準化的歸一化方法赦政。但數(shù)據(jù)歸一化對下游差異甲基化檢測至關(guān)重要胜宇。降采樣(Downsampling)通過減少reads序列數(shù)量，使其與相似序列數(shù)據(jù)相匹配恢着，從而實現(xiàn)歸一化桐愉。比對產(chǎn)生的BAM文件和甲基化提取產(chǎn)生的bedGraph文件都可以降采樣。在比對階段降采樣可能在時間和內(nèi)存上要求很高掰派，而在提取階段降采樣則需要較少的時間和內(nèi)存从诲，同時保證了相似的甲基化calling數(shù)量、檢測到的CpG位點靡羡、reads數(shù)分布和平均覆蓋度的準確性系洛。因此，在進行進一步的數(shù)據(jù)比較之前略步，建議先對bedGraph文件進行降采樣描扯，特別是對于差異甲基化區(qū)域（DMRs）的檢測。
（2）DMR
DMR（差異甲基化區(qū)域）檢測是核心甲基化分析之一趟薄，涉及對多個樣本中的基因組區(qū)域進行分析绽诚。最常見的應(yīng)用是在癌癥和正常樣本之間尋找可能作為生物標志物或揭示疾病生物學(xué)的異常甲基化區(qū)域『技澹基于DMR統(tǒng)計分析方法因軟件而異恩够，以下是一些主要方法。
BSmooth：使用局部似然平滑方法（local likelihood smoothing approach）來鑒定樣本特異性甲基化信息中的DMR岔帽。應(yīng)用Welch's t-test（Student's t-test變體）比較多個樣本玫鸟。DMR是具有觀察到的P值高于預(yù)定義β值的CpG位點。然而預(yù)定義閾值可能導(dǎo)致II類錯誤（假陰性）犀勒，從而影響結(jié)果屎飘。
BiSeq：通過納入錯誤發(fā)現(xiàn)率和β二項分布模型來解決這一問題妥曲，充分考慮到生物學(xué)重復(fù)。然后通過triangular kernel模型調(diào)整分層過程引起的顯著變化钦购，計算目標區(qū)域的統(tǒng)計顯著性檐盟。P值被歸一化、轉(zhuǎn)換為z分數(shù)押桃、平均值進行比較葵萎。
MethylSig：類似于BiSeq，應(yīng)用β二項分布模型來考慮reads覆蓋度和生物學(xué)意義唱凯。
Metilene：結(jié)合了二項分割和多變量Kolmogorov-Smirnov擬合優(yōu)度檢驗(K-S 檢驗)羡忘。這種非參數(shù)方法使用逐步分割基因組區(qū)域，被穩(wěn)步地最小化到CpG數(shù)量少于預(yù)定義下限的區(qū)域磕昼，或在統(tǒng)計顯著性上沒有改善的區(qū)域卷雕。這種方法對累積樣本分布的差異性更敏感。
methylKit：對單樣本情況使用Fisher's精確檢驗票从，對多重復(fù)樣本使用基于logistical回歸的統(tǒng)計方法計算組間差異漫雕。
Defiant：是一個獨立的程序，使用加權(quán)Welch擴展鑒定DMR峰鄙。對于只有一個重復(fù)的兩個樣本浸间，使用Fisher's精確檢驗，對于有多個重復(fù)的樣本吟榴，使用Welch's t-test魁蒜，基于覆蓋度對無偏樣本方差進行加權(quán)。
Benjamini-Hochberg：應(yīng)用于調(diào)整DMR鑒定中多重t檢驗的P值煤墙。數(shù)據(jù)分布本質(zhì)上是二項的梅惯，因為大多數(shù)甲基化分布要么是完全甲基化的，要么是完全未甲基化的仿野，表明二項分布模型性能優(yōu)于其他模型。
由于reads覆蓋度變化和人口統(tǒng)計參數(shù)（如性別她君、年齡和種族）的共變異對DMR檢測有強相關(guān)脚作，因此數(shù)據(jù)的歸一化轉(zhuǎn)換和協(xié)變量調(diào)整至關(guān)重要。
案例研究
WGBS數(shù)據(jù)的計算分析具有挑戰(zhàn)性缔刹，包括分析 FASTQ 讀長球涛、甲基化估計、位點注釋校镐、DMR 檢測和可視化亿扁。以下為WGBS數(shù)據(jù)分析的綜合案例研究。
（1）分析工具
（2）數(shù)據(jù)分析
①Pre-alignment
使用以下代碼檢測 R1 （WGBS_R1.fastq） 和 R2 （WGBS_R2.fastq） 讀取的原始 FASTQ 質(zhì)量：
?結(jié)果包括一個.html格式的摘要和一個帶有質(zhì)量數(shù)字的.zip文件夾鸟廓。MultiQC 將具有多個 FastQC 結(jié)果的樣品合并到一份.html報告中从祝。
②reads修剪
③比對和甲基化calling
Bismark：
BitMapperBS:
DMR的注釋和分析
易小結(jié)：
DNA甲基化和其他表觀基因組分析的動態(tài)標記與不同人類疾病的診斷和預(yù)后相關(guān)聯(lián)襟己。對甲基化結(jié)果的深刻且低偏倚解釋是下游生物學(xué)機制研究的核心。本綜述討論了分析WGBS數(shù)據(jù)的計算方法牍陌，并介紹了使用現(xiàn)有工具從原始reads檢測DMR所需的基本QC分析步驟擎浴。此外，還提出了甲基化文庫制備和數(shù)據(jù)處理中固有的潛在混雜因素及其對策毒涧。希望潛在用戶能夠理解WGBS的基本概念贮预，從而加速人類疾病的發(fā)現(xiàn)。
參考文獻：
1契讲、Gong T,Borgard H, Zhang Z, Chen S, Gao Z, Deng Y. Analysis and Performance Assessmentof the Whole Genome Bisulfite Sequencing Data Workflow: Currently AvailableTools and a Practical Guide to Advance DNA Methylation Studies. Small Methods.2022 Mar;6(3):e2101251. doi: 10.1002/smtd.202101251. PubMed PMID: 35064762.
