Bismark Bisulfite Mapper学习笔记（二）甲
(IV)Bismark甲基化提取
Bismark附带一个补充的bismark_methylation_extractor脚本，它对Bismark结果文件进行操作，并对每个C进行提取甲基化call。根据上下区域(CpG、CHG或CHH)，每个C的位置将被写入一个新的输出文件，因此甲基化的Cs将被标记为正read(+)，非甲基化的Cs将被标记为反read(-)。生成的文件可以导入到SeqMonk等基因组查看器中，然后就可以开始对甲基化数据的分析了。或者，甲基化提取的输出可以使用选项--bedGraph转换为一个bedGraph和coverage文件。这个步骤也可以使用独立脚本bismark2bedGraph从甲基化提取输出中完成。coverage文件也可以使用Import Data > Bismark (cov)直接导入到SeqMonk中。或者，bedGraph计数输出可用于生成全基因组胞嘧啶报告，报告基因组中每一个CpG(可选的每一个胞嘧啶)上的数字，不考虑是否被任何read覆盖。由于这种类型的报告对两条链上的胞嘧啶都具有信息性，因此输出的文件可能相当大。bedGraph到全基因组胞嘧啶报告转换也可以使用独立模块coverage2cytosine单独运行。
在Bismark版本0.6或更高版本中，bismark_methylation_extractor的默认输入格式是BAM/SAM。
bismark_methylation_extractor有几个可选项。如在甲基化call字符串里忽略第一个数字的位置；如删除一个限制性内切酶的位点(如果执行 RRBS非定向的BS-Seq库，可能需要删除每个read重组MspI位点，因为它会在第一个甲基化call里引入bias)。对于paired-end的reads，另一个有用的可选项是--no_overlap(默认情况下为on)，指定这个可选项将只对paired-end reads中间重叠部分的甲基化call提取一次(使用来自第一个read的calls，因为第一个read可能是错误率最低的)。
要获得可选项类型的完整列表:使用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
甲基化提取还可以自动检测比对模式，并自动设置上述选项。一个典型的代码包含可选的bedGraph输出的是这样的:
bismark_methylation_extractor --gzip --bedGraph --buffer_size 10G sample_bismark_bt2.bam
包括可选的全基因组胞嘧啶甲基化报告的典型代码是这样的:
bismark_methylation_extractor --gzip --bedGraph --buffer_size 10G --cytosine_report --genome_folder /path_to_genome_folder/ sample_bismark_bt2.bam
所以，对于这个例子，代码是：
$ 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
运行了大概40分钟（128G内存），生成以下文件：
-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链的互补链
OB – original bottom strand原始bottom链
CTOB – complementary to original bottom strand原始bottom链互补链
来自OT和CTOT的甲基化calls将提供原始top链上胞嘧啶甲基化位置的信息，来自OB和CTOB的甲基化calls将提供原始bottom链上胞嘧啶甲基化位置的信息。请注意，在Bismark比对步骤中指定--directional(默认模式)选项不会对CTOT或CTOB链生成任何报告，所以在这个例子中，我们只得到了OT和OB的甲基化报告信息。
由于胞嘧啶可以存在于三种不同的序列背景(CpG, CHG或CHH)中的任何一个中，bismark_methylation_extractor的默认输出将包含3个单独的输出文件，所以对于每一个OT/OB/CTOT/CTOB都会有各自的3个文件。在这个例子中，我们只有OT和OB，所以一共有6个文件生成。
如果对特定链的甲基化不感兴趣，那么所有可用的甲基化信息都可以合并到一个文件中(4条链的信息可以合并在一起)。使用--merge_non_CpG，将默认为三个输出文件(CpG-context、CHG-context和CHH-context)，或者两个输出文件(CpG-context和non-CpG-context)。(注意，对于non-CpG输出，会导致很大的文件大小)。
链特异性和背景相关的可选项都可以与--merge_non_CpG选项结合使用。
对于甲基化提取生成文件的解读：
（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
上面的内容里可以看到，每个文件里都是5列，这5列分别是：
- seq-ID （序列ID）
- methylation state （甲基化状态，“+”是甲基化的，“-”是非甲基化的）
- chromosome （染色体）
- start position (= end position) （起始/结束位置）
- 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，是一个基于0的基因组起始和基于1结束的坐标文件。它按染色体坐标排序，一共有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列分别是：
<染色体> <起始位置> <结束位置> <甲基化百分数>
由于甲基化百分数本身并不能提供在某一位置检测到的甲基化或非甲基化read的实际read覆盖率，bismark2bedGraph还会输出一个coverage文件(使用基于1的基因组坐标)，该文件具有两个附加列:
$ 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列的信息分别是：
<染色体> <起始位点> <结束位点> <甲基化百分数> <甲基化的count> <非甲基化的count>
（3）基因组范围内胞嘧啶甲基化报告
Bismark甲基化提取还可以选择性地输出全基因组胞嘧啶甲基化报告。coverage2cytosine模块也可以单独运行。该报告按染色体坐标进行排序，但也包含序列背景。在本例中输出的文件：
$ 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列的内容是：
<染色体> <位置> <链> <甲基化的count> <非甲基化的count> <C-context> <三核苷酸背景>
与bedGraph文件或coverage文件的主要区别在于，无论在实验中是否有任何read覆盖了它们，都将考虑top链和bottom链上的每个胞嘧啶。要使其工作，还必须使用--genome_folder <path>指定用于Bismark比对的基因组。对于bedGraph模式，默认情况下只考虑CpG背景中的胞嘧啶，但是可以通过使用--CX 扩展到任何序列背景中的胞嘧啶。但请注意，这可能意味着任何大型基因组的单个输出都超过11亿个胞嘧啶。
（4）M-bias plot
从Bismark v0.8.0开始，Bismark甲基化提取还可以生成甲基化bias图，显示read中每个可能位置的甲基化比例(Hansen et al., Genome Biology, 2012, 13:R83)。M-bias图的数据也写入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。同时还生成.png文件，该文件需要Perl模块GD::Graph(更具体地说，需要GD::Graph::lines和GD::Graph::colour两个模块)。如果系统上找不到GD::Graph，则只打印表。图中还包含每个位置的甲基化call的绝对数量(甲基化和非甲基化)。对于双端reads，将绘制两个单独的M-bias图。
Bismark自动生成的M-bias plot
(V) Bismark HTML 报告
脚本bismark2report使用Bismark比对报告，以及Bismark套件的其他报告(如去重、甲基化提取(splitting)或M-bias报告)来生成一个图形化的HTML报告页面。如果在同一文件夹中发现多个Bismark报告，那么将为每个报告生成一个单独的HTML报告，其中输出文件名来自Bismark比对报告文件。bismark2report尝试根据文件basename自动查找可选报告。
# --dir是指定输出文件夹
$ bismark2report --dir ./HTML_reports/
当然上面的代码还有其他可选的参数，比如：
--alignment_report FILE #你可以指定为比对报告生成HTML报告
--dedup_report FILE #指定为去重生成报告
--splitting_report FILE #为甲基化提取的splitting文件生成报告
--mbias_report FILE #为M-bias文件生成报告
--nucleotide_report #为核酸文件生成报告
上面的代码里，我没有指定为哪一个文件单独生成报告，所以bismark2report会为它找到的所有文件都生成一个报告，可以看一下bismark2report运行时弹出的信息：
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 总结报告
这个脚本使用几个样品(甚至数百个)的文件生成HTML报告。除非指定了某些BAM文件，否则bismark2summary首先识别文件夹中的Bismark BAM文件，然后根据输入文件basename自动检测Bismark比对、去重或甲基化提取(splitting)报告。
基本格式：
bismark2summary [options]
由于在这个例子中，只有一个样品，所以这一步就不练习了。具体的参数解释在：https://github.com/FelixKrueger/Bismark/tree/master/Docs
(VII) Bismark 核酸覆盖报告 (bam2nuc)
脚本bam2nuc读取BAM文件并计算read的单核苷酸和双核苷酸覆盖率(使用基因组序列而不是reads本身观察到的序列)，并将其与平均基因组序列组成进行比较。包含InDels的reads没有被考虑在内。含有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函数把这个表可视化：
总结
下面是一个对于不同甲基化测序文库类型的分析流程里要注意的步骤：
Bismark软件的学习先告一段落，具体的更多细节还请参照官网的参数介绍：这里。
