(IV)Bismark甲基化提取
Bismark附帶一個(gè)補(bǔ)充的bismark_methylation_extractor腳本，它對(duì)Bismark結(jié)果文件進(jìn)行操作，并對(duì)每個(gè)C進(jìn)行提取甲基化call。根據(jù)上下區(qū)域(CpG、CHG或CHH)，每個(gè)C的位置將被寫(xiě)入一個(gè)新的輸出文件，因此甲基化的Cs將被標(biāo)記為正read(+)，非甲基化的Cs將被標(biāo)記為反read(-)。生成的文件可以導(dǎo)入到SeqMonk等基因組查看器中，然后就可以開(kāi)始對(duì)甲基化數(shù)據(jù)的分析了?；蛘?，甲基化提取的輸出可以使用選項(xiàng)--bedGraph轉(zhuǎn)換為一個(gè)bedGraph和coverage文件。這個(gè)步驟也可以使用獨(dú)立腳本bismark2bedGraph從甲基化提取輸出中完成。coverage文件也可以使用Import Data > Bismark (cov)直接導(dǎo)入到SeqMonk中?；蛘撸琤edGraph計(jì)數(shù)輸出可用于生成全基因組胞嘧啶報(bào)告，報(bào)告基因組中每一個(gè)CpG(可選的每一個(gè)胞嘧啶)上的數(shù)字，不考慮是否被任何read覆蓋。由于這種類型的報(bào)告對(duì)兩條鏈上的胞嘧啶都具有信息性，因此輸出的文件可能相當(dāng)大。bedGraph到全基因組胞嘧啶報(bào)告轉(zhuǎn)換也可以使用獨(dú)立模塊coverage2cytosine單獨(dú)運(yùn)行。
在Bismark版本0.6或更高版本中，bismark_methylation_extractor的默認(rèn)輸入格式是BAM/SAM。
bismark_methylation_extractor有幾個(gè)可選項(xiàng)。如在甲基化call字符串里忽略第一個(gè)數(shù)字的位置；如刪除一個(gè)限制性內(nèi)切酶的位點(diǎn)(如果執(zhí)行 RRBS非定向的BS-Seq庫(kù)，可能需要?jiǎng)h除每個(gè)read重組MspI位點(diǎn)，因?yàn)樗鼤?huì)在第一個(gè)甲基化call里引入bias)。對(duì)于paired-end的reads，另一個(gè)有用的可選項(xiàng)是--no_overlap(默認(rèn)情況下為on)，指定這個(gè)可選項(xiàng)將只對(duì)paired-end reads中間重疊部分的甲基化call提取一次(使用來(lái)自第一個(gè)read的calls，因?yàn)榈谝粋€(gè)read可能是錯(cuò)誤率最低的)。
要獲得可選項(xiàng)類型的完整列表:使用bismark_methylation_extractor --help命令。
基本格式：
bismark_methylation_extractor [options] <filenames>
對(duì)于單端測(cè)序的文件：
bismark_methylation_extractor --gzip sample_bismark_bt2.bam
對(duì)于雙端測(cè)序的文件：
bismark_methylation_extractor --gzip sample_bismark_bt2_pe.bam
甲基化提取還可以自動(dòng)檢測(cè)比對(duì)模式，并自動(dòng)設(shè)置上述選項(xiàng)。一個(gè)典型的代碼包含可選的bedGraph輸出的是這樣的:
bismark_methylation_extractor --gzip --bedGraph --buffer_size 10G sample_bismark_bt2.bam
包括可選的全基因組胞嘧啶甲基化報(bào)告的典型代碼是這樣的:
bismark_methylation_extractor --gzip --bedGraph --buffer_size 10G --cytosine_report --genome_folder /path_to_genome_folder/ sample_bismark_bt2.bam
所以，對(duì)于這個(gè)例子，代碼是：
$ bismark_methylation_extractor --gzip --bedGraph --buffer_size 10G --cytosine_report --genome_folder /gpfs/home/fangy04/Bismark/hg38_genome/WholeGenomeFasta test_data_bismark_bt2.deduplicated.bam
運(yùn)行了大概40分鐘（128G內(nèi)存），生成以下文件：
-rw------- 1 fangy04 fangy04 80735 Dec 4 23:17 CHH_OT_test_data_bismark_bt2.deduplicated.txt.gz
-rw------- 1 fangy04 fangy04 13255 Dec 4 23:17 CpG_OT_test_data_bismark_bt2.deduplicated.txt.gz
-rw------- 1 fangy04 fangy04 44091 Dec 4 23:17 CHG_OT_test_data_bismark_bt2.deduplicated.txt.gz
-rw------- 1 fangy04 fangy04 80815 Dec 4 23:17 CHH_OB_test_data_bismark_bt2.deduplicated.txt.gz
-rw------- 1 fangy04 fangy04 13464 Dec 4 23:17 CpG_OB_test_data_bismark_bt2.deduplicated.txt.gz
-rw------- 1 fangy04 fangy04 42658 Dec 4 23:17 CHG_OB_test_data_bismark_bt2.deduplicated.txt.gz
-rw------- 1 fangy04 fangy04 780 Dec 4 23:17 test_data_bismark_bt2.deduplicated_splitting_report.txt
-rw------- 1 fangy04 fangy04 2945 Dec 4 23:17 test_data_bismark_bt2.deduplicated.M-bias.txt
-rw------- 1 fangy04 fangy04 12606 Dec 4 23:17 test_data_bismark_bt2.deduplicated.M-bias_R1.png
-rw------- 1 fangy04 fangy04 15750 Dec 4 23:17 test_data_bismark_bt2.deduplicated.bedGraph.gz
-rw------- 1 fangy04 fangy04 14627 Dec 4 23:17 test_data_bismark_bt2.deduplicated.bismark.cov.gz
-rw------- 1 fangy04 fangy04 211870669 Dec 4 23:59 test_data_bismark_bt2.deduplicated.CpG_report.txt.gz
上面生成的文件中，文件名里OB,OT的意思是：
OT – original top strand原始top鏈
CTOT – complementary to original top strand原始top鏈的互補(bǔ)鏈
OB – original bottom strand原始bottom鏈
CTOB – complementary to original bottom strand原始bottom鏈互補(bǔ)鏈
來(lái)自O(shè)T和CTOT的甲基化calls將提供原始top鏈上胞嘧啶甲基化位置的信息，來(lái)自O(shè)B和CTOB的甲基化calls將提供原始bottom鏈上胞嘧啶甲基化位置的信息。請(qǐng)注意，在Bismark比對(duì)步驟中指定--directional(默認(rèn)模式)選項(xiàng)不會(huì)對(duì)CTOT或CTOB鏈生成任何報(bào)告，所以在這個(gè)例子中，我們只得到了OT和OB的甲基化報(bào)告信息。
由于胞嘧啶可以存在于三種不同的序列背景(CpG, CHG或CHH)中的任何一個(gè)中，bismark_methylation_extractor的默認(rèn)輸出將包含3個(gè)單獨(dú)的輸出文件，所以對(duì)于每一個(gè)OT/OB/CTOT/CTOB都會(huì)有各自的3個(gè)文件。在這個(gè)例子中，我們只有OT和OB，所以一共有6個(gè)文件生成。
如果對(duì)特定鏈的甲基化不感興趣，那么所有可用的甲基化信息都可以合并到一個(gè)文件中(4條鏈的信息可以合并在一起)。使用--merge_non_CpG，將默認(rèn)為三個(gè)輸出文件(CpG-context、CHG-context和CHH-context)，或者兩個(gè)輸出文件(CpG-context和non-CpG-context)。(注意，對(duì)于non-CpG輸出，會(huì)導(dǎo)致很大的文件大小)。
鏈特異性和背景相關(guān)的可選項(xiàng)都可以與--merge_non_CpG選項(xiàng)結(jié)合使用。
對(duì)于甲基化提取生成文件的解讀：
（1）胞嘧啶在CpG/CHG/CHH背景的文件
用zcat查看壓縮的txt文件：
#top鏈在CpG背景下的甲基化信息
$ zcat CpG_OT_test_data_bismark_bt2.deduplicated.txt.gz | head -5
Bismark methylation extractor version v0.22.1
SRR020138.15024319_SALK_2029:7:100:1672:31_length=86 + chr2 10026473 Z
SRR020138.15024331_SALK_2029:7:100:1673:1282_length=86 + chr11 78025243 Z
SRR020138.15024338_SALK_2029:7:100:1673:1625_length=86 + chr3 197545315 Z
SRR020138.15024338_SALK_2029:7:100:1673:1625_length=86 + chr3 197545329 Z
#Bottom鏈在CHG背景下的甲基化信息
$ zcat CHG_OB_test_data_bismark_bt2.deduplicated.txt.gz | head -5
Bismark methylation extractor version v0.22.1
SRR020138.15024320_SALK_2029:7:100:1672:1164_length=86 - chr5 28344377 x
SRR020138.15024320_SALK_2029:7:100:1672:1164_length=86 - chr5 28344369 x
SRR020138.15024326_SALK_2029:7:100:1672:1418_length=86 - chr5 126882694 x
SRR020138.15024326_SALK_2029:7:100:1672:1418_length=86 - chr5 126882662 x
#top鏈在CHH背景下的甲基化信息
$ zcat CHH_OT_test_data_bismark_bt2.deduplicated.txt.gz | head -5
Bismark methylation extractor version v0.22.1
SRR020138.15024318_SALK_2029:7:100:1672:137_length=86 - chr12 129289592 h
SRR020138.15024319_SALK_2029:7:100:1672:31_length=86 - chr2 10026449 h
SRR020138.15024319_SALK_2029:7:100:1672:31_length=86 - chr2 10026467 h
SRR020138.15024319_SALK_2029:7:100:1672:31_length=86 - chr2 10026470 h
上面的內(nèi)容里可以看到，每個(gè)文件里都是5列，這5列分別是：
- seq-ID （序列ID）
- methylation state （甲基化狀態(tài)，“+”是甲基化的，“-”是非甲基化的）
- chromosome （染色體）
- start position (= end position) （起始/結(jié)束位置）
- methylation call （甲基化call）
z - C in CpG context - unmethylated
Z - C in CpG context - methylated
x - C in CHG context - unmethylated
X - C in CHG context - methylated
h - C in CHH context - unmethylated
H - C in CHH context - methylated
u - C in Unknown context (CN or CHN) - unmethylated
U - C in Unknown context (CN or CHN) - methylated
（2）bedGraph輸出文件
Bismark甲基化提取的選擇性輸出文件bedGraph，是一個(gè)基于0的基因組起始和基于1結(jié)束的坐標(biāo)文件。它按染色體坐標(biāo)排序，一共有4列：
$ zcat test_data_bismark_bt2.deduplicated.bedGraph.gz | head -5
track type=bedGraph
chr2 35131 35132 0
chr2 2744216 2744217 0
chr2 3292734 3292735 100
chr2 5367628 5367629 0
4列分別是：
<染色體> <起始位置> <結(jié)束位置> <甲基化百分?jǐn)?shù)>
由于甲基化百分?jǐn)?shù)本身并不能提供在某一位置檢測(cè)到的甲基化或非甲基化read的實(shí)際read覆蓋率，bismark2bedGraph還會(huì)輸出一個(gè)coverage文件(使用基于1的基因組坐標(biāo))，該文件具有兩個(gè)附加列:
$ zcat test_data_bismark_bt2.deduplicated.bismark.cov.gz | head -5
chr2 35132 35132 0 0 1
chr2 2744217 2744217 0 0 1
chr2 3292735 3292735 100 1 0
chr2 5367629 5367629 0 0 1
chr2 6471076 6471076 0 0 1
6列的信息分別是：
<染色體> <起始位點(diǎn)> <結(jié)束位點(diǎn)> <甲基化百分?jǐn)?shù)> <甲基化的count> <非甲基化的count>
（3）基因組范圍內(nèi)胞嘧啶甲基化報(bào)告
Bismark甲基化提取還可以選擇性地輸出全基因組胞嘧啶甲基化報(bào)告。coverage2cytosine模塊也可以單獨(dú)運(yùn)行。該報(bào)告按染色體坐標(biāo)進(jìn)行排序，但也包含序列背景。在本例中輸出的文件：
$ zcat test_data_bismark_bt2.deduplicated.CpG_report.txt.gz | head -5
chr2 10001 + 0 0 CG CGT
chr2 10002 - 0 0 CG CGN
chr2 10215 + 0 0 CG CGA
chr2 10216 - 0 0 CG CGG
chr2 10270 + 0 0 CG CGT
7列的內(nèi)容是：
<染色體> <位置> <鏈> <甲基化的count> <非甲基化的count> <C-context> <三核苷酸背景>
與bedGraph文件或coverage文件的主要區(qū)別在于，無(wú)論在實(shí)驗(yàn)中是否有任何read覆蓋了它們，都將考慮top鏈和bottom鏈上的每個(gè)胞嘧啶。要使其工作，還必須使用--genome_folder <path>指定用于Bismark比對(duì)的基因組。對(duì)于bedGraph模式，默認(rèn)情況下只考慮CpG背景中的胞嘧啶，但是可以通過(guò)使用--CX 擴(kuò)展到任何序列背景中的胞嘧啶。但請(qǐng)注意，這可能意味著任何大型基因組的單個(gè)輸出都超過(guò)11億個(gè)胞嘧啶。
（4）M-bias plot
從Bismark v0.8.0開(kāi)始，Bismark甲基化提取還可以生成甲基化bias圖，顯示read中每個(gè)可能位置的甲基化比例(Hansen et al., Genome Biology, 2012, 13:R83)。M-bias圖的數(shù)據(jù)也寫(xiě)入coverage文本文件:
$ head -8 test_data_bismark_bt2.deduplicated.M-bias.txt
CpG context
===========
position count methylated count unmethylated % methylation coverage
1 45 13 77.59 58
2 31 9 77.50 40
3 26 10 72.22 36
4 22 16 57.89 38
5 35 11 76.09 46
上面5列分別代表：
<read位置> <甲基化的count> <非甲基化的count> <甲基化%> <總coverage>
你可以用其他方法生成漂亮的圖表，例如使用R或Excel。同時(shí)還生成.png文件，該文件需要Perl模塊GD::Graph(更具體地說(shuō)，需要GD::Graph::lines和GD::Graph::colour兩個(gè)模塊)。如果系統(tǒng)上找不到GD::Graph，則只打印表。圖中還包含每個(gè)位置的甲基化call的絕對(duì)數(shù)量(甲基化和非甲基化)。對(duì)于雙端reads，將繪制兩個(gè)單獨(dú)的M-bias圖。
(V) Bismark HTML 報(bào)告
腳本bismark2report使用Bismark比對(duì)報(bào)告，以及Bismark套件的其他報(bào)告(如去重、甲基化提取(splitting)或M-bias報(bào)告)來(lái)生成一個(gè)圖形化的HTML報(bào)告頁(yè)面。如果在同一文件夾中發(fā)現(xiàn)多個(gè)Bismark報(bào)告，那么將為每個(gè)報(bào)告生成一個(gè)單獨(dú)的HTML報(bào)告，其中輸出文件名來(lái)自Bismark比對(duì)報(bào)告文件。bismark2report嘗試根據(jù)文件basename自動(dòng)查找可選報(bào)告。
# --dir是指定輸出文件夾
$ bismark2report --dir ./HTML_reports/
當(dāng)然上面的代碼還有其他可選的參數(shù)，比如：
--alignment_report FILE #你可以指定為比對(duì)報(bào)告生成HTML報(bào)告
--dedup_report FILE #指定為去重生成報(bào)告
--splitting_report FILE #為甲基化提取的splitting文件生成報(bào)告
--mbias_report FILE #為M-bias文件生成報(bào)告
--nucleotide_report #為核酸文件生成報(bào)告
上面的代碼里，我沒(méi)有指定為哪一個(gè)文件單獨(dú)生成報(bào)告，所以bismark2report會(huì)為它找到的所有文件都生成一個(gè)報(bào)告，可以看一下bismark2report運(yùn)行時(shí)彈出的信息：
Writing Bismark HTML report to >> ./HTML_reports/test_data_bismark_bt2_SE_report.html <<
==============================================================================================================
Using the following alignment report: > test_data_bismark_bt2_SE_report.txt <
Processing alignment report test_data_bismark_bt2_SE_report.txt ...
Complete
Using the following deduplication report: > test_data_bismark_bt2.deduplication_report.txt <
Processing deduplication report test_data_bismark_bt2.deduplication_report.txt ...
Complete
Using the following splitting report: > test_data_bismark_bt2.deduplicated_splitting_report.txt <
Processing splitting report test_data_bismark_bt2.deduplicated_splitting_report.txt ...
Complete
Using the following M-bias report: > test_data_bismark_bt2.deduplicated.M-bias.txt <
Processing M-bias report test_data_bismark_bt2.deduplicated.M-bias.txt ...
Complete
No nucleotide coverage report file specified, skipping this step
==============================================================================================================
生成的html文件用瀏覽器打開(kāi)后是這樣的：
(VI)Bismark 總結(jié)報(bào)告
這個(gè)腳本使用幾個(gè)樣品(甚至數(shù)百個(gè))的文件生成HTML報(bào)告。除非指定了某些BAM文件，否則bismark2summary首先識(shí)別文件夾中的Bismark BAM文件，然后根據(jù)輸入文件basename自動(dòng)檢測(cè)Bismark比對(duì)、去重或甲基化提取(splitting)報(bào)告。
基本格式：
bismark2summary [options]
由于在這個(gè)例子中，只有一個(gè)樣品，所以這一步就不練習(xí)了。具體的參數(shù)解釋在：https://github.com/FelixKrueger/Bismark/tree/master/Docs
(VII) Bismark 核酸覆蓋報(bào)告 (bam2nuc)
腳本bam2nuc讀取BAM文件并計(jì)算read的單核苷酸和雙核苷酸覆蓋率(使用基因組序列而不是reads本身觀察到的序列)，并將其與平均基因組序列組成進(jìn)行比較。包含InDels的reads沒(méi)有被考慮在內(nèi)。含有Ns的單核苷酸或雙核苷酸也會(huì)被忽略。
bam2nuc處理Bismark單端和雙端文件(自動(dòng)確定)。BAM和CRAM文件都可以作為輸入，但是請(qǐng)注意，CRAM文件需要Samtools 1.2或更高版本。
基本格式：
bam2nuc [options] --genome_folder <path> [input.(bam|cram)]
代碼：
# --dir輸出文件夾
# --genome_folder 基因組文件夾
$ bam2nuc --dir ./Nucleotide_Coverage_report --genome_folder /gpfs/home/fangy04/Bismark/hg38_genome/WholeGenomeFasta test_data_bismark_bt2.deduplicated.bam
生成文件：
$ cat test_data_bismark_bt2.deduplicated.nucleotide_stats.txt
(di-)nucleotide count sample percent sample count genomic percent genomic coverage
A 80512 32.84 866420001 29.52 0.000
C 40381 16.47 598683433 20.40 0.000
G 56325 22.98 600854940 20.47 0.000
T 67932 27.71 868918077 29.61 0.000
AA 28791 11.98 286827281 9.77 0.000
AC 10702 4.45 147939993 5.04 0.000
AG 20855 8.68 205417428 7.00 0.000
AT 18616 7.75 226234956 7.71 0.000
CA 15612 6.50 212726856 7.25 0.000
CC 8419 3.50 151371046 5.16 0.000
CG 2020 0.84 29303965 1.00 0.000
CT 13529 5.63 205281198 6.99 0.000
GA 18284 7.61 175501875 5.98 0.000
GC 9234 3.84 124675968 4.25 0.000
GG 14579 6.07 152424560 5.19 0.000
GT 13024 5.42 148252316 5.05 0.000
TA 15598 6.49 191363599 6.52 0.000
TC 10970 4.57 174696228 5.95 0.000
TG 18727 7.79 213708670 7.28 0.000
TT 21287 8.86 289149328 9.85 0.000
我們可以使用前面提到的bismark2report函數(shù)把這個(gè)表可視化：
總結(jié)
下面是一個(gè)對(duì)于不同甲基化測(cè)序文庫(kù)類型的分析流程里要注意的步驟：
Bismark軟件的學(xué)習(xí)先告一段落，具體的更多細(xì)節(jié)還請(qǐng)參照官網(wǎng)的參數(shù)介紹：這里。
