全基因組DNA甲基化測序數(shù)據(jù)工作流程分析和性能評估之分析軟件比較
DNA甲基化與轉(zhuǎn)錄調(diào)控、基因組印記、干細胞分化、胚胎發(fā)育和炎癥等過程有關。DNA甲基化異�？赡芙沂炯膊顟B(tài)，包括癌癥和神經(jīng)系統(tǒng)疾病。因此，人類基因組中5-甲基胞嘧啶（5mC）分布和位點是一個重要的研究方向。全基因組重亞硫酸鹽測序(Whole-genome bisulfite sequencing, WGBS)是一種用于分析DNA甲基化的高通量方法，本文綜述了WGBS的關鍵步驟，總結(jié)了可用且最新的分析工具，比較了比對算法，并分享了數(shù)據(jù)處理的示例代碼，為科研人員的DNA甲基化研究提供參考。
標題：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年3月
影響因子：10.7
WGBS文庫制備
廣泛使用的亞硫酸鹽轉(zhuǎn)化文庫構建方法：Accel-NGS Methyl-Seq（Accel）、TruSeq DNA Methylation（TruSeq）和SPlinted ligation adapter tagging（SPLAT）。
WGBS分析流程
在處理原始測序數(shù)據(jù)時，生物信息學分析流程的可重復性和一致性至關重要。由于WGBS數(shù)據(jù)通常很龐大，因此需要大量的計算資源、內(nèi)存和存儲空間。這就要求分析流程不僅要穩(wěn)定，還要高效地使用內(nèi)存和時間。
表1：WGBS比對軟件在速度上的比較。
四個比對工具在甲基化calling后的比對質(zhì)量顯示，BWA-METH和gemBS有最高的唯一比對率和最少的未比對reads（圖3A）。在每個染色體的平均CpG覆蓋度≤×10時存在微小差異，在≥×20覆蓋度時，BWA-METH比其他方法有略好的覆蓋度（圖3B）。DNA片段兩端未甲基化的Cs數(shù)量對于從比對結(jié)果中獲得合格的DNA甲基化calling至關重要。使用M-bias圖對每個流程比對結(jié)果的影響規(guī)模非常相似，盡管在不同文庫之間存在較大差異（圖3C）。使用比對軟件在基因組上calling的CpGs一致性顯示，在每個注釋區(qū)域和轉(zhuǎn)錄起始位點(TSSs)周圍的平均甲基化水平上顯示出可比的基因組富集分布（圖3D-E）。與此同時，甲基化extraction水平不受比對軟件選擇的影響，因為Bismark與其他三種比對軟件相比，甲基化calling相關性只略降低（圖4）。
（4）比對質(zhì)量評估（Alignment quality control）
比對后的QC在WGBS中至關重要，WGBS的內(nèi)在混雜變量會使甲基化估計偏向于過高或過低的估計，主要通過M-bias圖來檢測。在亞硫酸鹽處理過程中可能會發(fā)生部分（不完全）甲基化，此時可以觀察到C和T的peaks值，這通常會導致檢測過高。高于98.5%的比率可以確保沒有偏差。可以在所有背景中（CpG、CHG和CHH）向甲基化樣本中添加帶有未甲基化C的Spike序列，然后計數(shù)未甲基化C和T數(shù)量，并計算添加序列的轉(zhuǎn)化率。
同時，由于估計過低導致的偏差會捕獲到假陰性的甲基化位點。例如，通過酶介導的碎片化雙鏈DNA末端修復會在片段兩端引入未甲基化的C，從而導致人為的甲基化水平低估。這在M-bias圖中反映為問題兩端的平均甲基化水平急劇下降，這應在提取甲基化之前予以丟棄。亞硫酸鹽介導的降解是WGBS中偏差的主要來源，因為降解非隨機，發(fā)生在未甲基化的胞嘧啶上，這些胞嘧啶從文庫中舍棄。這導致許多隨后的序列偏差和整體甲基化的高估。
（5）甲基化信息提�。∕ethylation extraction）
在經(jīng)過BitMapperBS、BWA-METH、Bismark和GemBS等比對軟件進行比對后，推薦利用MethylDackel進行甲基化信息提取。例如用MethylDackel對BitMapperBS比對后的甲基化信息提取。甲基化評估通過比較測序reads和參考基因組進行，如果在參考基因組中某個位點顯示為C，在上述位點注意到C時，就分配100%的甲基化，當指示為T時，則分配0%的甲基化。計算加權平均值，并在計算該位點的C和T數(shù)量后，將其指定為最終的甲基化水平。如10/10 Cs顯示完全胞嘧啶甲基化，6/10 Cs顯示部分甲基化（60%），0/10 Cs代表未甲基化胞嘧啶。在提取之前，對兩個鏈中每個位點的平均甲基化水平進行M-bias分析，以識別提取reads時的基本技術偏差作為修剪參考。從理論上講，reads應該是恒定的，但每對中的第一和第二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ù)歸一化對下游差異甲基化檢測至關重要。降采樣(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ū)域進行分析。最常見的應用是在癌癥和正常樣本之間尋找可能作為生物標志物或揭示疾病生物學的異常甲基化區(qū)域�；贒MR統(tǒng)計分析方法因軟件而異，以下是一些主要方法。
BSmooth：使用局部似然平滑方法（local likelihood smoothing approach）來鑒定樣本特異性甲基化信息中的DMR。應用Welch's t-test（Student's t-test變體）比較多個樣本。DMR是具有觀察到的P值高于預定義β值的CpG位點。然而預定義閾值可能導致II類錯誤（假陰性），從而影響結(jié)果。
BiSeq：通過納入錯誤發(fā)現(xiàn)率和β二項分布模型來解決這一問題，充分考慮到生物學重復。然后通過triangular kernel模型調(diào)整分層過程引起的顯著變化，計算目標區(qū)域的統(tǒng)計顯著性。P值被歸一化、轉(zhuǎn)換為z分數(shù)、平均值進行比較。
MethylSig：類似于BiSeq，應用β二項分布模型來考慮reads覆蓋度和生物學意義。
Metilene：結(jié)合了二項分割和多變量Kolmogorov-Smirnov擬合優(yōu)度檢驗(K-S 檢驗)。這種非參數(shù)方法使用逐步分割基因組區(qū)域，被穩(wěn)步地最小化到CpG數(shù)量少于預定義下限的區(qū)域，或在統(tǒng)計顯著性上沒有改善的區(qū)域。這種方法對累積樣本分布的差異性更敏感。
methylKit：對單樣本情況使用Fisher's精確檢驗，對多重復樣本使用基于logistical回歸的統(tǒng)計方法計算組間差異。
Defiant：是一個獨立的程序，使用加權Welch擴展鑒定DMR。對于只有一個重復的兩個樣本，使用Fisher's精確檢驗，對于有多個重復的樣本，使用Welch's t-test，基于覆蓋度對無偏樣本方差進行加權。
Benjamini-Hochberg：應用于調(diào)整DMR鑒定中多重t檢驗的P值。數(shù)據(jù)分布本質(zhì)上是二項的，因為大多數(shù)甲基化分布要么是完全甲基化的，要么是完全未甲基化的，表明二項分布模型性能優(yōu)于其他模型。
由于reads覆蓋度變化和人口統(tǒng)計參數(shù)（如性別、年齡和種族）的共變異對DMR檢測有強相關，因此數(shù)據(jù)的歸一化轉(zhuǎn)換和協(xié)變量調(diào)整至關重要。
案例研究
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)標記與不同人類疾病的診斷和預后相關聯(lián)。對甲基化結(jié)果的深刻且低偏倚解釋是下游生物學機制研究的核心。本綜述討論了分析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 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.
標題：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年3月
影響因子：10.7
WGBS文庫制備
廣泛使用的亞硫酸鹽轉(zhuǎn)化文庫構建方法：Accel-NGS Methyl-Seq（Accel）、TruSeq DNA Methylation（TruSeq）和SPlinted ligation adapter tagging（SPLAT）。
圖1：WGBS的亞硫酸鹽轉(zhuǎn)后文庫制備方法。
WGBS分析流程
在處理原始測序數(shù)據(jù)時，生物信息學分析流程的可重復性和一致性至關重要。由于WGBS數(shù)據(jù)通常很龐大，因此需要大量的計算資源、內(nèi)存和存儲空間。這就要求分析流程不僅要穩(wěn)定，還要高效地使用內(nèi)存和時間。
圖2：WGBS數(shù)據(jù)處理和分析流程概述和原則。
a）質(zhì)量評估：預比對質(zhì)控(QC)以進行測序質(zhì)量評估、接頭檢測和序列偏倚。然后對原始數(shù)據(jù)reads進行質(zhì)量過濾和接頭修剪。
b）數(shù)據(jù)比對：處理后的數(shù)據(jù)使用專門設計的軟件與參考基因組進行比對，生成標準的SAM/BAM格式結(jié)果。比對軟件會考慮到胞嘧啶和胸腺嘧啶的不匹配分布變化，并進行相應的調(diào)整。
c）比對后的質(zhì)量評估：可以通過一些比對軟件或使用外部軟件如Samtools或Qualimap來獲得比對后的質(zhì)控結(jié)果。甲基化偏差(M-bias)圖表示reads中每個位點的平均甲基化水平，應該是恒定的。如果不是，則需要在甲基化calling中修剪偏倚reads。
d）甲基化calling：甲基化calling過程會在參考基因組中迭代，為胞嘧啶背景(CpG、CHG和CHH)生成原始甲基化要求。在某些情況下，比對軟件的選擇可能會影響甲基化分析結(jié)果。最終的甲基化calling報告會列出每個位點的甲基化百分比和覆蓋率。
e）甲基化calling結(jié)果的注釋和分析：甲基化calling結(jié)果可以進行注釋或進行差異甲基化區(qū)域(DMR)分析，并可以通過多種方法（如BedGraph）進行可視化。處理大量數(shù)據(jù)，如人類基因組，可能需要降維以提高速度并避免內(nèi)存“耗盡”。
（1）質(zhì)量評估（Quality Assessment）
預比對質(zhì)量評估確保了原始數(shù)據(jù)的正確輸入并提高了可比對性。多維度評估審查了低質(zhì)量reads、接頭序列、污染序列和重復reads。合成測序（增加測序的堿基數(shù)量）會降低質(zhì)量，尤其是Illumina平臺。在測序過程中，DNA簇中的單個分子可能無法成功合成。錯誤合成的分子聚類會導致motif calling錯誤，這些錯誤與片段長強正相關。文庫構建中長片段比例越高，可能會導致更多的calling錯誤、較低的Phred得分和更高的不匹配率，特別是對于成對末端對齊。建議保留質(zhì)量分數(shù)≥30的堿基，這表示對motif calling準確性99.9%。DNA測序反應讀長通常比DNA片段長，片段兩端的接頭序列可能被錯誤地測序，這會引入構成性的甲基化C，并在隨后的甲基化calling中引起偏差，因此建議提前檢查并修剪接頭。數(shù)據(jù)可能被污染，包括在文庫制備步驟中添加引物和載體序列，以及為校準測序反應中的序列質(zhì)量而添加的Phi X噬菌體DNA。過度表示的污染物通常呈現(xiàn)顯著高于核基因組的reads深度，導致在比對過程中失敗。當同一DNA片段被雙重計數(shù)時，可能會發(fā)生重復reads（兩個具有相同基因位點的reads）。
長插入片段可能會從更高比例的高質(zhì)量reads、較少接頭污染和更高有效reads深度（增加基因組覆蓋效率）中受益。TruSeq文庫制備需要特別注意重復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ù)，能夠自動識別技術序列）和Trimmomatic軟件對低質(zhì)量reads、接頭序列、污染序列、重復序列進行去除。這一步驟稱之為修剪(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對文庫的默認設置具有方向性，在Trim Galore中應指定成對末端的特性，以保持成對reads，以便進一步比對。兩個文件的reads數(shù)不匹配和reads名稱不一致會觸發(fā)警告。
（3）比對（Alignment）
WGBS的原理是對未甲基化的胞嘧啶(Cs)通過亞硫酸鹽處理轉(zhuǎn)化為胸腺嘧啶(Ts)，同時保留甲基化的Cs。理想情況下，當reads序列比對到參考基因組時，可以識別未甲基化的Cs。然而亞硫酸鹽轉(zhuǎn)化帶來數(shù)據(jù)比對兩大計算挑戰(zhàn)：C-T比對不匹配和序列復雜性降低。C-T比對不匹配指的是在測序reads中，T可能與參考基因組中的C比對，反之則不然。序列復雜性降低使得難以區(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，并將序列復雜性提高到確保與參考基因組唯一比對水平，而三字符比對軟件在BS轉(zhuǎn)化reads中去除Cs，降低序列復雜性，增加了模糊比對位點的機會。但覆蓋差異和M偏差僅在基因組中的高相似區(qū)域中表現(xiàn)出來，因此在如人類基因組的長序列reads比對不太相關。因此在比對軟件選擇中優(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ì)量之間的關系。
用于亞硫酸鹽處理reads比對算法在運行時間和比對質(zhì)量上的差異會影響下游的甲基化calling。在Accel-NGS MethylSeq、SPLAT和TruSeq等WGBS框架以及TrueMethyl和EMSeq的氧化亞硫酸鹽測序中，對廣泛應用的比對軟件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)定。
b）數(shù)據(jù)比對：處理后的數(shù)據(jù)使用專門設計的軟件與參考基因組進行比對，生成標準的SAM/BAM格式結(jié)果。比對軟件會考慮到胞嘧啶和胸腺嘧啶的不匹配分布變化，并進行相應的調(diào)整。
c）比對后的質(zhì)量評估：可以通過一些比對軟件或使用外部軟件如Samtools或Qualimap來獲得比對后的質(zhì)控結(jié)果。甲基化偏差(M-bias)圖表示reads中每個位點的平均甲基化水平，應該是恒定的。如果不是，則需要在甲基化calling中修剪偏倚reads。
d）甲基化calling：甲基化calling過程會在參考基因組中迭代，為胞嘧啶背景(CpG、CHG和CHH)生成原始甲基化要求。在某些情況下，比對軟件的選擇可能會影響甲基化分析結(jié)果。最終的甲基化calling報告會列出每個位點的甲基化百分比和覆蓋率。
e）甲基化calling結(jié)果的注釋和分析：甲基化calling結(jié)果可以進行注釋或進行差異甲基化區(qū)域(DMR)分析，并可以通過多種方法（如BedGraph）進行可視化。處理大量數(shù)據(jù)，如人類基因組，可能需要降維以提高速度并避免內(nèi)存“耗盡”。
（1）質(zhì)量評估（Quality Assessment）
預比對質(zhì)量評估確保了原始數(shù)據(jù)的正確輸入并提高了可比對性。多維度評估審查了低質(zhì)量reads、接頭序列、污染序列和重復reads。合成測序（增加測序的堿基數(shù)量）會降低質(zhì)量，尤其是Illumina平臺。在測序過程中，DNA簇中的單個分子可能無法成功合成。錯誤合成的分子聚類會導致motif calling錯誤，這些錯誤與片段長強正相關。文庫構建中長片段比例越高，可能會導致更多的calling錯誤、較低的Phred得分和更高的不匹配率，特別是對于成對末端對齊。建議保留質(zhì)量分數(shù)≥30的堿基，這表示對motif calling準確性99.9%。DNA測序反應讀長通常比DNA片段長，片段兩端的接頭序列可能被錯誤地測序，這會引入構成性的甲基化C，并在隨后的甲基化calling中引起偏差，因此建議提前檢查并修剪接頭。數(shù)據(jù)可能被污染，包括在文庫制備步驟中添加引物和載體序列，以及為校準測序反應中的序列質(zhì)量而添加的Phi X噬菌體DNA。過度表示的污染物通常呈現(xiàn)顯著高于核基因組的reads深度，導致在比對過程中失敗。當同一DNA片段被雙重計數(shù)時，可能會發(fā)生重復reads（兩個具有相同基因位點的reads）。
長插入片段可能會從更高比例的高質(zhì)量reads、較少接頭污染和更高有效reads深度（增加基因組覆蓋效率）中受益。TruSeq文庫制備需要特別注意重復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ù)，能夠自動識別技術序列）和Trimmomatic軟件對低質(zhì)量reads、接頭序列、污染序列、重復序列進行去除。這一步驟稱之為修剪(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對文庫的默認設置具有方向性，在Trim Galore中應指定成對末端的特性，以保持成對reads，以便進一步比對。兩個文件的reads數(shù)不匹配和reads名稱不一致會觸發(fā)警告。
（3）比對（Alignment）
WGBS的原理是對未甲基化的胞嘧啶(Cs)通過亞硫酸鹽處理轉(zhuǎn)化為胸腺嘧啶(Ts)，同時保留甲基化的Cs。理想情況下，當reads序列比對到參考基因組時，可以識別未甲基化的Cs。然而亞硫酸鹽轉(zhuǎn)化帶來數(shù)據(jù)比對兩大計算挑戰(zhàn)：C-T比對不匹配和序列復雜性降低。C-T比對不匹配指的是在測序reads中，T可能與參考基因組中的C比對，反之則不然。序列復雜性降低使得難以區(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，并將序列復雜性提高到確保與參考基因組唯一比對水平，而三字符比對軟件在BS轉(zhuǎn)化reads中去除Cs，降低序列復雜性，增加了模糊比對位點的機會。但覆蓋差異和M偏差僅在基因組中的高相似區(qū)域中表現(xiàn)出來，因此在如人類基因組的長序列reads比對不太相關。因此在比對軟件選擇中優(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ì)量之間的關系。
用于亞硫酸鹽處理reads比對算法在運行時間和比對質(zhì)量上的差異會影響下游的甲基化calling。在Accel-NGS MethylSeq、SPLAT和TruSeq等WGBS框架以及TrueMethyl和EMSeq的氧化亞硫酸鹽測序中，對廣泛應用的比對軟件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)定。
表1：WGBS比對軟件在速度上的比較。
四個比對工具在甲基化calling后的比對質(zhì)量顯示，BWA-METH和gemBS有最高的唯一比對率和最少的未比對reads（圖3A）。在每個染色體的平均CpG覆蓋度≤×10時存在微小差異，在≥×20覆蓋度時，BWA-METH比其他方法有略好的覆蓋度（圖3B）。DNA片段兩端未甲基化的Cs數(shù)量對于從比對結(jié)果中獲得合格的DNA甲基化calling至關重要。使用M-bias圖對每個流程比對結(jié)果的影響規(guī)模非常相似，盡管在不同文庫之間存在較大差異（圖3C）。使用比對軟件在基因組上calling的CpGs一致性顯示，在每個注釋區(qū)域和轉(zhuǎn)錄起始位點(TSSs)周圍的平均甲基化水平上顯示出可比的基因組富集分布（圖3D-E）。與此同時，甲基化extraction水平不受比對軟件選擇的影響，因為Bismark與其他三種比對軟件相比，甲基化calling相關性只略降低（圖4）。
圖3：四種比對算法的比對質(zhì)量比較。
圖4：所有 CpG 位點的全基因組甲基化水平的相關性矩陣
（4）比對質(zhì)量評估（Alignment quality control）
比對后的QC在WGBS中至關重要，WGBS的內(nèi)在混雜變量會使甲基化估計偏向于過高或過低的估計，主要通過M-bias圖來檢測。在亞硫酸鹽處理過程中可能會發(fā)生部分（不完全）甲基化，此時可以觀察到C和T的peaks值，這通常會導致檢測過高。高于98.5%的比率可以確保沒有偏差。可以在所有背景中（CpG、CHG和CHH）向甲基化樣本中添加帶有未甲基化C的Spike序列，然后計數(shù)未甲基化C和T數(shù)量，并計算添加序列的轉(zhuǎn)化率。
同時，由于估計過低導致的偏差會捕獲到假陰性的甲基化位點。例如，通過酶介導的碎片化雙鏈DNA末端修復會在片段兩端引入未甲基化的C，從而導致人為的甲基化水平低估。這在M-bias圖中反映為問題兩端的平均甲基化水平急劇下降，這應在提取甲基化之前予以丟棄。亞硫酸鹽介導的降解是WGBS中偏差的主要來源，因為降解非隨機，發(fā)生在未甲基化的胞嘧啶上，這些胞嘧啶從文庫中舍棄。這導致許多隨后的序列偏差和整體甲基化的高估。
（5）甲基化信息提�。∕ethylation extraction）
在經(jīng)過BitMapperBS、BWA-METH、Bismark和GemBS等比對軟件進行比對后，推薦利用MethylDackel進行甲基化信息提取。例如用MethylDackel對BitMapperBS比對后的甲基化信息提取。甲基化評估通過比較測序reads和參考基因組進行，如果在參考基因組中某個位點顯示為C，在上述位點注意到C時，就分配100%的甲基化，當指示為T時，則分配0%的甲基化。計算加權平均值，并在計算該位點的C和T數(shù)量后，將其指定為最終的甲基化水平。如10/10 Cs顯示完全胞嘧啶甲基化，6/10 Cs顯示部分甲基化（60%），0/10 Cs代表未甲基化胞嘧啶。在提取之前，對兩個鏈中每個位點的平均甲基化水平進行M-bias分析，以識別提取reads時的基本技術偏差作為修剪參考。從理論上講，reads應該是恒定的，但每對中的第一和第二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ù)歸一化對下游差異甲基化檢測至關重要。降采樣(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ū)域進行分析。最常見的應用是在癌癥和正常樣本之間尋找可能作為生物標志物或揭示疾病生物學的異常甲基化區(qū)域�；贒MR統(tǒng)計分析方法因軟件而異，以下是一些主要方法。
BSmooth：使用局部似然平滑方法（local likelihood smoothing approach）來鑒定樣本特異性甲基化信息中的DMR。應用Welch's t-test（Student's t-test變體）比較多個樣本。DMR是具有觀察到的P值高于預定義β值的CpG位點。然而預定義閾值可能導致II類錯誤（假陰性），從而影響結(jié)果。
BiSeq：通過納入錯誤發(fā)現(xiàn)率和β二項分布模型來解決這一問題，充分考慮到生物學重復。然后通過triangular kernel模型調(diào)整分層過程引起的顯著變化，計算目標區(qū)域的統(tǒng)計顯著性。P值被歸一化、轉(zhuǎn)換為z分數(shù)、平均值進行比較。
MethylSig：類似于BiSeq，應用β二項分布模型來考慮reads覆蓋度和生物學意義。
Metilene：結(jié)合了二項分割和多變量Kolmogorov-Smirnov擬合優(yōu)度檢驗(K-S 檢驗)。這種非參數(shù)方法使用逐步分割基因組區(qū)域，被穩(wěn)步地最小化到CpG數(shù)量少于預定義下限的區(qū)域，或在統(tǒng)計顯著性上沒有改善的區(qū)域。這種方法對累積樣本分布的差異性更敏感。
methylKit：對單樣本情況使用Fisher's精確檢驗，對多重復樣本使用基于logistical回歸的統(tǒng)計方法計算組間差異。
Defiant：是一個獨立的程序，使用加權Welch擴展鑒定DMR。對于只有一個重復的兩個樣本，使用Fisher's精確檢驗，對于有多個重復的樣本，使用Welch's t-test，基于覆蓋度對無偏樣本方差進行加權。
Benjamini-Hochberg：應用于調(diào)整DMR鑒定中多重t檢驗的P值。數(shù)據(jù)分布本質(zhì)上是二項的，因為大多數(shù)甲基化分布要么是完全甲基化的，要么是完全未甲基化的，表明二項分布模型性能優(yōu)于其他模型。
由于reads覆蓋度變化和人口統(tǒng)計參數(shù)（如性別、年齡和種族）的共變異對DMR檢測有強相關，因此數(shù)據(jù)的歸一化轉(zhuǎn)換和協(xié)變量調(diào)整至關重要。
案例研究
WGBS 數(shù)據(jù)的計算分析具有挑戰(zhàn)性，包括分析 FASTQ 讀長、甲基化估計、位點注釋、DMR 檢測和可視化。以下為WGBS數(shù)據(jù)分析的綜合案例研究。
（1）分析工具
（2）數(shù)據(jù)分析
① Pre-alignment
圖5：FASTQ 文件元素
使用以下代碼檢測 R1 （WGBS_R1.fastq） 和 R2 （WGBS_R2.fastq） 讀取的原始 FASTQ 質(zhì)量：
圖6：FastQC的質(zhì)控報告
② reads修剪
③ 比對和甲基化calling
Bismark：
BitMapperBS:
DMR的注釋和分析
易小結(jié)：
DNA甲基化和其他表觀基因組分析的動態(tài)標記與不同人類疾病的診斷和預后相關聯(lián)。對甲基化結(jié)果的深刻且低偏倚解釋是下游生物學機制研究的核心。本綜述討論了分析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 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.
標簽：
DNA甲基化
Copyright(C) 1998-2024 生物器材網(wǎng) 電話：021-64166852;13621656896 E-mail：info@bio-equip.com