Bismark软件使用入门
作者：马可菠萝 童蒙
编辑：angelica
一、甲基化文库建库原理
借助重亚硫酸盐（Bisulfite）转换的方法被认为是金标准，该化学试剂会将未甲基化的C进行转化，胞嘧啶（C）转换为尿嘧啶（U），后续的PCR扩增过程U会被胸腺嘧啶（T）替代，而，甲基化的C则不受影响。在了解如何使用Bismark之前，需要先了解DNA甲基化文库的结构。
通常，甲基化文库构建方法有两种：
-
一种是先片段化，然后末修，加上特殊处理的接头（甲基化的）后，进行CT转化，之后进行扩增，如下图：
该文库为链特异性文库，使用bismark默认参数比对即可。
-
另一种是先CT转换，然后扩增加接头，如下图。
该文库为非链特异性文库，需要使用PBAT文库参数。
通常而言，甲基化建库均为链特异性方式构建：
- 测序read1都是来自于原始的基因组，经历了C->T转化，来自于正链的叫OT，负链叫OB；
- 测序read2都是来自于原始基因组的互补链，经历了G->A转换，称之为CTOT或者CTOB。
二、bismark比对分析
bismark的分析分为三步：
基因组索引构建：需要将基因组进行相应的转换，模拟bisulfite处理后基因组，这样才能用于比对。在每个新物种比对之前，都需要进行该处理。bismark内部默认调用bowite软件进行比对。
序列比对：输入为下机序列、基因组和比对参数，bismark会输出比对结果和甲基化检测的结果，默认格式为BAM，同时还会有一个统计报告。
甲基化位点提取：这一步是可选的，会利用bismark的比对结果来获得甲基化信息。这一步会把甲基化分为不同的类别（CG，CHG和CHH），获得链特异性信息，并且提供过滤参数；也可以对甲基化信息进行调整，或者进行更多的深入分析。
01.基因组索引
使用bismark_genome_preparation来对基因组进行处理，输入是给定的目录，目录里面是.fa或者.fasta文件。
bismark会生成两个目录，一个是C->T转换的基因组，一个是G->A转换的基因组共两套基因组，之后可以用bowtie2-build来建立索引。
特别注意，索引路径中必需包含一个基因组文件。
运行命令示例如下：
bismark_genome_preparation
--path_to_bowtie /path/bowtie2/ 这里使用bowtie2
<path_to_genome_folder> 索引路径
--verbose 输出log信息
索引输出目录结果如下：
Bisulfite_Genome
├── CT_conversion
│ ├── BS_CT.1.bt2
│ ├── BS_CT.2.bt2
│ ├── BS_CT.3.bt2
│ ├── BS_CT.4.bt2
│ ├── BS_CT.rev.1.bt2
│ ├── BS_CT.rev.2.bt2
│ └── genome_mfa.CT_conversion.fa
└── GA_conversion
├── BS_GA.1.bt2
├── BS_GA.2.bt2
├── BS_GA.3.bt2
├── BS_GA.4.bt2
├── BS_GA.rev.1.bt2
├── BS_GA.rev.2.bt2
└── genome_mfa.GA_conversion.fa
02.比对
比对的过程需要提供两个信息
- 输入目录：已经建立好的索引目录，即基因组索引步骤中的索引目录
- 输入文件：测序reads，fastq格式
- 其余参数都是可选，自行调整使用
输出包括：
- 比对的结果，BAM文件
- 报告文件：包含各种统计信息和甲基化统计结果
2.2.1、比对原理
一般情况下，由于建库的时候会产生四种可能性的序列，原始的正链（OT），原始的互补链（OB），原始正链的互补链（CTOT），原始互补链的互补链（CTOB），因此bismark会进行进行四条路线的比对，如下图。
- 左边的C->T，是原始链的转化方式，然后跟两种参考基因组进行比对，得到OT和OB的结果；
- 右边的G->A，是扩增链的转化方式，然后跟两种参考基因组进行比对，得到CTOT和CTOB的结果；
上面示意图提到建库方式是有方向性的，因此bismark只需要进行两次比对即可，而不需要进行4次的比对，如果文库是非链特异性文库，需要加上--non_directional的参数，bismark会进行4次比对，选择最优的结果进行输出。
针对先转化后加接头的文库（PBAT），常见如单细胞甲基化文库，需要添加--pbat参数进行比对。
2.2.2、运行bismark
bismark
--bowtie2 使用bowtie2
-N 0 允许错配数目
--un 过滤多处匹配reads
-o 输出目录
<genome_folder> 索引路径
-1 read_1.fq.qz 测序read1
-2 read_2.fq.gz 测序read2
03.输 出
bam文件：可以利用SeqMonk进行可视化；也可以用于继续的去除dup等
比对输出路径示例如下：
# 样本：A
A_P_R1_bismark_bt2_pe.bam
A_P_R1_bismark_bt2_PE_report.txt
A_P_R2.fq.gz_unmapped_reads_1.fq.gz
A_P_R2.fq.gz_unmapped_reads_2.fq.gz
比对输出路径中包含一个比对后的文件（BAM格式）和一个比对报告文件以及未比对到参考基因组上的reads文件（--un参数）。
若测序数据为双端数据文件中会有PE或pe标签，若为单端数据则会有SE或se字符加以区别。
针对链特异性文库，若进行单端比对read1使用bismark单端参数即可，read2需要添加--pbat参数；而非链特异性文库--pbat参数的使用则相反。
通常，若某物种比对率较低，可先进行双端比对同时输出未必对的reads，随后使用未比对的read1和read2分别进行单端比对并合并输出的BAM文件，作为下游的输入。
比对报告文件示例如下：
Bismark report for: A_P_R1.fq.gz and A_P_R2.fq.gz
Bismark was run with Bowtie 2 against the bisulfite genome of /Index/ with the specified options: -q --score-min L,0,-0.2 --reorder --ignore-quals --no-mixed --no-discordant --dovetail --maxins 1000
Option '--directional' specified (default mode): alignments to complementary strands (CTOT, CTOB) were ignored (i.e. not performed)
Final Alignment report
======================
Sequence pairs analysed in total: 43930319
Number of paired-end alignments with a unique best hit: 26314484
Mapping efficiency: 60.0%
Sequence pairs with no alignments under any condition: 16147263
Sequence pairs did not map uniquely: 1468572
Sequence pairs which were discarded because genomic sequence could not be extracted: 2
Number of sequence pairs with unique best (first) alignment came from the bowtie output:
CT/GA/CT: 13199639 ((converted) top strand)
GA/CT/CT: 0 (complementary to (converted) top strand)
GA/CT/GA: 0 (complementary to (converted) bottom strand)
CT/GA/GA: 13114843 ((converted) bottom strand)
Number of alignments to (merely theoretical) complementary strands being rejected in total: 0
Final Cytosine Methylation Report
=================================
Total number of C's analysed: 1191141867
Total methylated C's in CpG context: 33527315
Total methylated C's in CHG context: 17696056
Total methylated C's in CHH context: 68376839
Total methylated C's in Unknown context: 97299
Total unmethylated C's in CpG context: 20813883
Total unmethylated C's in CHG context: 239796745
Total unmethylated C's in CHH context: 810931029
Total unmethylated C's in Unknown context: 301701
C methylated in CpG context: 61.7%
C methylated in CHG context: 6.9%
C methylated in CHH context: 7.8%
C methylated in unknown context (CN or CHN): 24.4%
Bismark completed in 0d 7h 30m 40s
报告展示了唯一比对上的read对数目等信息（Final Alignment report部分）以及有reads支持发生甲基化的C的数目统计（Final Cytosine Methylation Report部分），注意这里并非真正意义上的甲基化C，仍需要下游的分析。
04.去除duplicated序列
由于PCR的扩增，可能会存在比对到基因组相同位置的序列，重复测序的序列需要被去除。如果是RRBS文库，建议不要进行这一步骤。
对于单选测序序列，去dup的时候，会根据染色体、起始坐标和链的方向来确定是否是duplication；
对于双端测序序列，Read1的染色体、起始坐标，Read2的起始坐标会用于去除duplication。需要注意的是，去dup的时候需要注意read1和read2是相邻的，而不是按照坐标排序的。
去除dup的参数示例如下：
deduplicate_bismark
-p A_P_R1_bismark_bt2_pe.bam bismark比对输出的bam文件
--bam bam格式输出
--samtools_path <samtools_dir> SAMtools软件路径
--output_dir <out_dir> 输出路径
输出路径会生成一个去除dup后的bam文件和一个去除dup的统计文件。
三、提取甲基化位点
使用bismark_methylation_extractor来提取每个C位点的甲基化的情况，结果文件可以导入到SeqMonk中。
可以使用--bedGraph，来产生bedgraph和coverage的文件；--count可以产生count文件，输出每个CpG位点的深度文件。
bismark_methylation_extractor
--pair-end 指定双端数据
--comprehensive 输出CHG CHH CpG的甲基化信息
--no-overlap 去除reads重叠区域的bias
--bedGraph 输出bedGraph文件
--counts 每个C上甲基化reads和非甲基化reads的数目
--report 一个甲基化summay
--cytosine_report 输出全基因组所有CpG
--genome_folder <path_to_reference_genome>
input.bam 输入文件
-o >output_dir> 输出路径
bismark对不同状态的胞嘧啶进行了约定并记录到了BAM文件和中间文件中：
X 代表CHG中甲基化的C
x 代笔CHG中非甲基化的C
H 代表CHH中甲基化的C
h 代表CHH中非甲基化的C
Z 代表CpG中甲基化的C
z 代表CpG中非甲基化的C
U 代表其他情况的甲基化C(CN或者CHN)
u 代表其他情况的非甲基化C (CN或者CHN)
bismark_methylation_extractor的一个非常重要的输出文件为全基因组胞嘧啶报告文件（*CX_report.txt），该文件是一个非常核心的记录胞嘧啶深度信息的文件，可进行甲基化水平的计算、区域甲基化信号文件的转化、差异甲基化水平的计算等。
文件示例如下：
chr1 14716 + 1 6 CG CGT
chr1 14717 - 10 1 CG CGA
chr1 14741 + 2 1 CG CGG
chr1 14742 - 0 1 CG CGC
各列信息解释如下：
第1列 : 染色体编号
第2列 : 染色体起始位置
第3列 : 正负链信息
第4列 : 甲基化碱基数目
第5列 : 非甲基化碱基数目
第6列 : CG/CHG/CHH
第7列 : C背景序列
基于该文件可进行单个胞嘧啶位点的甲基化水平（C/C+T）计算，比如第一行：甲基化是水平为1/(1+6)；针对指定区域的甲基化水平计算则取区域内的胞嘧啶的总C/(总C+总T)即可。
四、bismark的优缺点
4.1、优点
（1）比对和mC Calling比较准确，受众广泛，使用率很高
（2）软件维护性好，更新快，bug修复及时
（3）软件成体系，集成了比对、去重、mC Calling分析
4.2、缺点
（1）比对和mC Calling速度慢，真的很慢，可以通过拆分输入数据为小份作为输出和增大计算资源解决。
（2）bismak相匹配的直接进行下游甲基化水平、图形可视化等分析的套件较少，需要自行计算或进行格式转化作为第三方软件的输入。
（3）说明文档不够系统和详细，在”外行“看来理解从建库到mC Calling有一定的挑战。
bismark的结果只是研究甲基化信息的第一步也是关键的一步。未来，期望bismark在速度上，和在对甲基化水平层面分析的支持上多一下改进。
