(IV)Bismark甲基化提取
Bismark附帶一個(gè)補(bǔ)充的bismark_methylation_extractor腳本，它對Bismark結(jié)果文件進(jìn)行操作，并對每個(gè)C進(jìn)行提取甲基化call蹦狂。根據(jù)上下區(qū)域(CpG允华、CHG或CHH)，每個(gè)C的位置將被寫入一個(gè)新的輸出文件袍榆，因此甲基化的Cs將被標(biāo)記為正read(+)廉羔，非甲基化的Cs將被標(biāo)記為反read(-)溉痢。生成的文件可以導(dǎo)入到SeqMonk等基因組查看器中，然后就可以開始對甲基化數(shù)據(jù)的分析了憋他∈适遥或者，甲基化提取的輸出可以使用選項(xiàng)--bedGraph轉(zhuǎn)換為一個(gè)bedGraph和coverage文件举瑰。這個(gè)步驟也可以使用獨(dú)立腳本bismark2bedGraph從甲基化提取輸出中完成。coverage文件也可以使用Import Data > Bismark (cov)直接導(dǎo)入到SeqMonk中蔬螟〈搜福或者，bedGraph計(jì)數(shù)輸出可用于生成全基因組胞嘧啶報(bào)告旧巾，報(bào)告基因組中每一個(gè)CpG(可選的每一個(gè)胞嘧啶)上的數(shù)字耸序，不考慮是否被任何read覆蓋。由于這種類型的報(bào)告對兩條鏈上的胞嘧啶都具有信息性鲁猩，因此輸出的文件可能相當(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庫男应，可能需要?jiǎng)h除每個(gè)read重組MspI位點(diǎn)，因?yàn)樗鼤诘谝粋€(gè)甲基化call里引入bias)娱仔。對于paired-end的reads沐飘，另一個(gè)有用的可選項(xiàng)是--no_overlap(默認(rèn)情況下為on)，指定這個(gè)可選項(xiàng)將只對paired-end reads中間重疊部分的甲基化call提取一次(使用來自第一個(gè)read的calls，因?yàn)榈谝粋€(gè)read可能是錯(cuò)誤率最低的)耐朴。
要獲得可選項(xiàng)類型的完整列表:使用bismark_methylation_extractor --help命令借卧。
基本格式：
bismark_methylation_extractor [options] <filenames>
對于單端測序的文件：
bismark_methylation_extractor --gzip sample_bismark_bt2.bam
對于雙端測序的文件：
bismark_methylation_extractor --gzip sample_bismark_bt2_pe.bam
甲基化提取還可以自動檢測比對模式，并自動設(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
所以铐刘，對于這個(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ǔ)鏈
來自O(shè)T和CTOT的甲基化calls將提供原始top鏈上胞嘧啶甲基化位置的信息，來自O(shè)B和CTOB的甲基化calls將提供原始bottom鏈上胞嘧啶甲基化位置的信息俯艰。請注意捡遍，在Bismark比對步驟中指定--directional(默認(rèn)模式)選項(xiàng)不會對CTOT或CTOB鏈生成任何報(bào)告，所以在這個(gè)例子中竹握，我們只得到了OT和OB的甲基化報(bào)告信息画株。
由于胞嘧啶可以存在于三種不同的序列背景(CpG, CHG或CHH)中的任何一個(gè)中，bismark_methylation_extractor的默認(rèn)輸出將包含3個(gè)單獨(dú)的輸出文件啦辐，所以對于每一個(gè)OT/OB/CTOT/CTOB都會有各自的3個(gè)文件谓传。在這個(gè)例子中，我們只有OT和OB芹关，所以一共有6個(gè)文件生成续挟。
如果對特定鏈的甲基化不感興趣，那么所有可用的甲基化信息都可以合并到一個(gè)文件中(4條鏈的信息可以合并在一起)侥衬。使用--merge_non_CpG诗祸，將默認(rèn)為三個(gè)輸出文件(CpG-context、CHG-context和CHH-context)轴总，或者兩個(gè)輸出文件(CpG-context和non-CpG-context)直颅。(注意，對于non-CpG輸出怀樟，會導(dǎo)致很大的文件大小)功偿。
鏈特異性和背景相關(guān)的可選項(xiàng)都可以與--merge_non_CpG選項(xiàng)結(jié)合使用。
對于甲基化提取生成文件的解讀：
（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ù)本身并不能提供在某一位置檢測到的甲基化或非甲基化read的實(shí)際read覆蓋率关拒，bismark2bedGraph還會輸出一個(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ū)別在于归露，無論在實(shí)驗(yàn)中是否有任何read覆蓋了它們洲脂，都將考慮top鏈和bottom鏈上的每個(gè)胞嘧啶。要使其工作剧包，還必須使用--genome_folder <path>指定用于Bismark比對的基因組恐锦。對于bedGraph模式，默認(rèn)情況下只考慮CpG背景中的胞嘧啶疆液，但是可以通過使用--CX 擴(kuò)展到任何序列背景中的胞嘧啶一铅。但請注意，這可能意味著任何大型基因組的單個(gè)輸出都超過11億個(gè)胞嘧啶堕油。
（4）M-bias plot
從Bismark v0.8.0開始潘飘，Bismark甲基化提取還可以生成甲基化bias圖，顯示read中每個(gè)可能位置的甲基化比例(Hansen et al., Genome Biology, 2012, 13:R83)掉缺。M-bias圖的數(shù)據(jù)也寫入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(更具體地說艰毒，需要GD::Graph::lines和GD::Graph::colour兩個(gè)模塊)。如果系統(tǒng)上找不到GD::Graph搜囱，則只打印表丑瞧。圖中還包含每個(gè)位置的甲基化call的絕對數(shù)量(甲基化和非甲基化)。對于雙端reads犬辰，將繪制兩個(gè)單獨(dú)的M-bias圖。
(V) Bismark HTML 報(bào)告
腳本bismark2report使用Bismark比對報(bào)告冰单，以及Bismark套件的其他報(bào)告(如去重幌缝、甲基化提取(splitting)或M-bias報(bào)告)來生成一個(gè)圖形化的HTML報(bào)告頁面。如果在同一文件夾中發(fā)現(xiàn)多個(gè)Bismark報(bào)告诫欠，那么將為每個(gè)報(bào)告生成一個(gè)單獨(dú)的HTML報(bào)告涵卵，其中輸出文件名來自Bismark比對報(bào)告文件。bismark2report嘗試根據(jù)文件basename自動查找可選報(bào)告荒叼。
# --dir是指定輸出文件夾
$ bismark2report --dir ./HTML_reports/
當(dāng)然上面的代碼還有其他可選的參數(shù)轿偎，比如：
--alignment_report FILE #你可以指定為比對報(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)告
上面的代碼里，我沒有指定為哪一個(gè)文件單獨(dú)生成報(bào)告被廓，所以bismark2report會為它找到的所有文件都生成一個(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文件用瀏覽器打開后是這樣的：
(VI)Bismark 總結(jié)報(bào)告
這個(gè)腳本使用幾個(gè)樣品(甚至數(shù)百個(gè))的文件生成HTML報(bào)告。除非指定了某些BAM文件，否則bismark2summary首先識別文件夾中的Bismark BAM文件昆婿，然后根據(jù)輸入文件basename自動檢測Bismark比對球碉、去重或甲基化提取(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沒有被考慮在內(nèi)。含有Ns的單核苷酸或雙核苷酸也會被忽略能庆。
bam2nuc處理Bismark單端和雙端文件(自動確定)施禾。BAM和CRAM文件都可以作為輸入，但是請注意相味，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è)對于不同甲基化測序文庫類型的分析流程里要注意的步驟：
Bismark軟件的學(xué)習(xí)先告一段落，具體的更多細(xì)節(jié)還請參照官網(wǎng)的參數(shù)介紹：這里丰涉。
