官方手册（2020 Sep更新）
https://github.com/FelixKrueger/Bismark/tree/master/Docs

（一）在开始前你需要了解的：

Bismark是用Perl编写的，并在命令行里运行。另外，你需要安装Bowtie 2 或者HISAT2 。 对于0.14或者更高的版本，Bismark可以使用并行运行，同时进行比对和甲基化提取步骤。你可以搜索--parallel/--multicore来看更多的细节。

首先你需要下载参考基因组，并把它放在一个文件夹里。参考基因组可以在 Ensembl或者NCBI 网站下载。对于下面讲到的例子，你需要下载Homo sapiens基因组。Bismark支持FastA格式的基因组序列，文件后缀可以是.fa或.fasta。可以输入单个fasta或者多个fasta文件。

教程的示例数据是test_dataset.fastq，下载地址here(it contains 10,000 reads in FastQ format, Phred33 qualities, 50 bp long reads, from a human directional BS-Seq library)。

（二）Bismark-主要信息

（1）什么是Bismark?

Bismark是一套工具，可高效率的分析亚硫酸氢盐序列(BS-Seq)数据。Bismark用参考基因组进行比对，同时进行胞嘧啶甲基化提取。Bismark是用Perl编写的，在命令行运行。用亚硫酸氢盐处理的reads用短Reads对准软件Bowtie 2或HISAT2进行比对。因此，需要在你的电脑上安装Bowtie 2(或HISAT2)。

（2）Bismark的安装

由于我使用的是服务器里的软件，所以直接调用就行了。如果想在自己的电脑上安装Bismark，可以在这个网站下载：here。
然后解压：

tar xzf bismark_v0.X.Y.tar.gz

（3）硬件要求

Bismark将参考基因组保存在内存中，除此之外，还运行了多达四个Bowtie 2的并行实例。内存的使用取决于参考基因组的大小。对于一个大的真核生物基因组(人类或小鼠)，我们的经验大约要用12GB的内存。因此，我们建议在具有5个CPU核和至少12 GB RAM的机器上运行Bismark。Bowtie 2的内存需求有点大(可能是为了允许gapped比对)。在使用Bowtie 2运行Bismark时，我们推荐至少有5个核和> 16GB RAM的系统。
比对速度很大程度上取决于所使用的read长度和比对参数。允许许多mismatches和使用较短的seed长度往往是相当慢的。

（4）这个软件支持什么样的BS-Seq文件？

Bismark支持在以下条件下亚硫酸氢盐处理的reads(全基因组shotgun BS-Seq (WGSBS)、reducated -representation BS-Seq (RRBS)或PBAT-Seq(后亚硫酸氢盐接头标签)进行比对:

1.序列格式可以是FastQ或FastA
2.单端或双端reads
3.输入文件可以被解压的或被gzip压缩的(以.gz结尾)
4.可变reads长度
5.定向或非定向的BS-Seq文库

下面是一个比对模式的完整列表：

注意：Bismark目前只支持Illumina平台的数据，尚不支持SOLiD平台的reads。

（5）Bismark是如何工作的?

序列reads首先由bisulfite进行转换(C->T)和反向(G->A转换正向链)版本，然后再与同样转换过的基因组进行比对(也包括C->T和G->A转换)。这样比对后就有4种结果，选取最佳比对的那一种后，然后将其与正常基因组序列进行比较，从而推断reads中所有胞嘧啶位置的甲基化状态。如果一个比对具有唯一的最佳比对得分(如as:i字段报告的那样)，则认为该read是惟一比对。如果一个read的比对产生了几个具有相同数量的 mismatches或相同比对分数(如AS:i字段)，那么这个read(或一个read-pair)将被完全丢弃。

光看文字可能比较难理解为什么会有4种结果，下面的图应该能比较直观的方便理解：

（6）Bismark比对和提取甲基化报告

完成后，Bismark生成一个包含以下信息的运行报告:

1.使用的比对参数总结
2.已分析的序列数
3.具有唯一最佳比对的序列数(mapping efficiency)

4.统计总结唯一最佳比对来自哪一条亚硫酸氢盐链
5.分析的胞嘧啶数目
6.甲基化和未甲基化胞嘧啶的数目

在CpG、CHG或CHH环境中胞嘧啶甲基化百分比(其中H可以是A、T或C)。该百分比根据公式分别计算:

% methylation (context) = 100 * methylated Cs (context) / (methylated Cs (context) + unmethylated Cs (context)).
需要强调的是，甲基化百分比值(背景)只是在映射步骤中直接执行的一个非常粗略的计算。经过后处理或过滤后的实际甲基化水平可能有所不同。

（三）运行Bismark

运行Bismark分为三个主要步骤:
首先，需要对感兴趣的基因组进行亚硫酸氢盐转换和并构建索引，以允许Bowtie进行比对。对于每个基因组，这一步只需要执行一次。注意，Bowtie 2和HISAT2需要不同的索引，因为它们的索引不兼容。

Bismark read比对步骤。只需指定一个要分析的文件，参考基因组和比对参数。Bismark将生成一个结合了比对/甲基化call的输出(默认是BAM格式)以及一个运行统计报告。

Bismark甲基化提取。这个步骤是可选的，它将从Bismark比对输出中提取甲基化信息。运行这个额外的步骤允许将甲基化信息分解到不同的背景中，获得特异链的甲基化信息并提供一些过滤选项。你还可以选择将甲基化信息排序到bedGraph/coverage文件中，或者进一步处理它们到全基因组胞嘧啶甲基化报告中。

下面几节将更详细地通过示例描述这些步骤。

(I) Bismark 基因组准备

这个脚本只需要运行一次，以准备感兴趣的亚硫酸氢盐比对基因组。你需要指定一个目录，其中包含你想要比对reads的基因组(请注意，bismark_genome_prepare脚本期望该文件夹中包含FastA文件，扩展名为.fa或. FastA，每个文件有单个或多个序列条目)。Bismark会在这个目录中创建两个单独的文件夹，一个用于C->T转换基因组，另一个用于G->A转换基因组。在创建了C->T和G->A版本的基因组后，它们将被使用bowtie2-build(或hisat2-build)并行构建索引。一旦C->T和G->A基因组索引都被创建，你就不需要再次使用该脚本。这一步根据基因组的大小，有可能需要几个小时来运行！

注：这一步我使用了服务器里的128G内存，运行大约2小时

请再次注意，Bowtie 2和HISAT2的索引是不兼容的!
这里我们使用的是bowtie2。
先看一下基本格式：

bismark_genome_preparation [options] <path_to_genome_folder>

代码：

$ bismark_genome_preparation --path_to_aligner /gpfs/share/apps/bowtie2/2.4.1/bin/ /gpfs/home/fangy04/Bismark/hg38_genome/WholeGenomeFasta

因为我是在服务器里运行的，所以我可以指定调用多少G的内存以及多少个cores（使用自己电脑运行的同学请忽略下面代码，直接跳到下面“(II) Bismark 比对步骤”，这里代码仅作为我自己的笔记使用）：

$ srun --mem=128G --cpus-per-task=1 --time=8:00:00 --pty bismark_genome_preparation --path_to_aligner /gpfs/share/apps/bowtie2/2.4.1/bin/ /gpfs/home/fangy04/Bismark/hg38_genome/WholeGenomeFasta

或：

$ sbatch genome_preparation.sh

#!/bin/bash
#SBATCH --job-name=bismark_genome_preparation           # Job name
#SBATCH --mail-type=END,FAIL                            # Mail events (NONE, BEGIN, END, FAIL, ALL)
#SBATCH --mail-user=yan.fang@***.org                    # Where to send mail
#SBATCH --ntasks=1                                      # Run on a single CPU
#SBATCH --mem=128gb                                     # Job memory request
#SBATCH --time=08:00:00                                 # Time limit hrs:min:sec
#SBATCH --output=/gpfs/home/fangy04/log_files/bismark_%j.log           # Standard output and error log
#SBATCH -p cpu_short

module load bismark/0.22.1
module load python/cpu/3.6.5
module load bowtie2/2.4.1
module load perl/5.28.0

bismark_genome_preparation --path_to_aligner /gpfs/share/apps/bowtie2/2.4.1/bin/ /gpfs/home/fangy04/Bismark/hg38_genome/WholeGenomeFasta

这一步运行完成后，看一下都有哪些输出文件：

$ tree #这是我存放原始基因组的文件夹
.
`-- WholeGenomeFasta
    |-- Bisulfite_Genome #这个文件夹是运行后生成的，可以看到里面有两个子文件夹
    |   |-- CT_conversion #CT转换的基因组文件夹
    |   |   |-- 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 #GA转换的基因组文件夹
    |       |-- 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
    |-- GenomeSize.xml
    |-- genome.dict
    |-- genome.fa
    `-- genome.fa.fai

4 directories, 18 files

(II) Bismark 比对步骤

这一步代表了亚硫酸氢盐比对和call甲基化的实际部分。Bismark要求用户只需指定两件事:

1.包含感兴趣的基因组的目录。这个文件夹必须包含未修改的基因组(如.fa或.fasta文件)，以及在Bismark基因组准备步骤中生成的两个亚硫酸氢盐基因组子目录。
2.要分析的序列文件(FastQ或FastA格式)。
所有其他参数都是可选的。

对于每个序列文件或每组配对的paired-end测序文件，Bismark生成一个比对和甲基化call输出文件，以及一个报告文件，详细说明比对和甲基化call统计信息，用于记录保存。

在运行Bismark之前，我们建议花些时间使用FastQC对原始序列文件进行质量控制。FastQC可能能够发现与你的BS-Seq文件相关的异常，如高碱基calling错误率或污染序列，如PCR引物或Illumina接头。许多误差来源会对比对效率和/或比对位置产生不利影响，因此也可能影响call甲基化和得出的结论。有必要的话，需要进行Trimming。但是这里由于是示例文件，这一步就不做了。在实际的数据处理中，这一步是要进行的。

如果没有指定其他选项，Bismark将使用一组默认值，其中一些是:

（1）使用Bowtie 2

1.使用Bowtie 2是默认模式
2.如果没有指定到Bowtie 2的路径，则假定bowtie2在PATH中
3.标准的比对使用multi-seed长度为20bp且0错配。可以使用选项-L和-N修改这些参数
4.标准比对使用默认的最小比对评分函数L,0，-0.2，即f(x) = 0 + -0.2 * x(其中x是读取长度)。对于75bp的read，这意味着在比对失效之前，read的最低比对分数可以是-15。这大致等于read中出现2次错配或1-2个indel(或两者的组合)。使用--score_min <func>设置严格度。

即使用户没有被要求指定额外的比对选项，但通常建议这样做(例如，当默认参数太严格时)。要查看选项的完整列表，请在命令行中键入bismark --help。

（2）定向的BS-Seq文库 (默认)

对于一个给定的locus，亚硫酸氢盐处理DNA和随后的PCR扩增可以产生四条(亚硫酸氢盐转换)链。根据使用的接头，BS-Seq库可以用两种不同的方式构建:

1.如果库是定向的，则只对原top链(OT)或原bottom链(OB)的(亚硫酸氢盐转换)版本进行测序。即使与OT (CTOT)互补链和OB (CTOB)互补链在BS-PCR步骤中产生，它们也不会被测序，因为它们在5 '端携带了错误的接头。默认情况下，Bismark只对OT和OB链执行2次read比对，因此忽略了来自互补链的比对，因为它们理论上不应该出现在BS-Seq库中。
2.或者，也可以构建BS-Seq文库，使在BS-PCR中产生的所有4条不同的链，以大致相同的可能性进入测序文库。在本例中，所有四条链(OT、CTOT、OB、CTOB)都可以产生有效的比对，这个库称为非定向的。指定--non_directional将指示Bismark使用所有四种比对的输出。

再总结一下:与原始top链的比对或与原始top链的互补链(OT和CTOT)都将产生top链上胞嘧啶的甲基化信息。与原始bottom链或与原始bottom链互补的链(OB和CTOB)比对，都会产生bottom链胞嘧啶的甲基化信息，也就是说，它们似乎会产生参考基因组top链G位置的甲基化信息。

有关如何提取四种不同比对链的甲基化信息的更多信息，请参阅下面关于Bismark methylation extractor的部分。

这一步的代码基本格式：

bismark [options] --genome <genome_folder> {-1 <mates1> -2 <mates2> | <singles>}

在我们这个例子里，示例数据是单端测序，所以代码是：

$ bismark --genome /gpfs/home/fangy04/Bismark/hg38_genome/WholeGenomeFasta sample.fastq.gz

（3）Bismark输出的文件是什么样的？

从0.6.x版本以后，输出的都是BAM/SAM格式。
运行后生成的文件：

（4）Bismark BAM/SAM输出文件

默认情况下，Bismark为所有比对模式生成SAM输出。请注意，报告的质量值是用Sanger格式(Phred 33)编码的。
每一列如下：

QNAME：序列ID
FLAG:这个FLAG试图考虑到亚硫酸氢盐读取的来源(这不同于普通的DNA比对FLAG!)
RNAME：染色体
POS：起始位置
MAPQ：为Bowtie 2和HISAT2计算
CIGAR
RNEXT
PNEXT
TLEN
SEQ
QUAL：Phred33
NM-tag
MD-tag：每个错配的碱基(14)XM-tag（甲基化calling字符）
XR-tag：read比对转换状态(16)XG-tag(基因组比对转换状态)

（5）数据可视化

要查看比对上的reads的位置，可以使用染色体、开始和结束位置将Bismark输出文件导入基因组查看器，比如SeqMonk(这对于识别基因组中显示大量比对上的Reads的区域非常有用)。比对的输出也可以用于后续步骤如去重(每个位置只允许1 条read，去除PCR artefacts)或过滤亚硫酸氢盐转换相关non-bisulfite错配的数量。

这里我就先不练习如何使用SeqMonk可视化了，如果后续需要我再深入学习

与亚硫酸氢盐转换相关的非亚硫酸氢盐错配是在基因组中有T，而在BS-read中有C的错配位置。这种错配可能是由于亚硫酸氢盐read比对的方式造成的。为了不引入甲基化reads的偏倚，包含这种错配的reads不会自动从比对输出中删除。值得注意的是，这些位置即使没有执行甲基化calls，含有亚硫酸氢盐转换相关的非亚硫酸氢盐错配可能导致错误的比对（如果使用特别宽松的比对参数）。

（6）Methylation call

看一下上一步输出的bam文件：

$ samtools view test_data_bismark_bt2.bam | head -n 5
SRR020138.15024317_SALK_2029:7:100:1672:902_length=86   16      chr1    57333005        42      50M     *       0   0         TTCTTTCCCATCCCATAAATCCTAAAAATAATAAAAAATCATCCCCAAAT       @@:AC@<=+@?+8)@BCCCA=6BCCCCCCCCCCCCCCCCACB=<88BCCA      NM:i:11       MD:Z:14G2G6G0G0G0G4G1G1G0G10G1   XM:Z:..............z..h......hhhh....h.h.hh..........h. XR:Z:CT XG:Z:GA
SRR020138.15024318_SALK_2029:7:100:1672:137_length=86   0       chr12   129289551       8       50M     *       0   0         AAAAAAAAAAAAAAGAAAAAAAAGAAAAAGAAAAGGAAAAGTAAAAAAAA       =@CAA=@B@CB=98%:AB?>@56/=3<=<)>B@:*=:=61%,<A@@1+12      NM:i:2        MD:Z:41C5G2      XM:Z:.........................................h........ XR:Z:CT XG:Z:CT
SRR020138.15024319_SALK_2029:7:100:1672:31_length=86    0       chr2    10026448        42      50M     *       0   0         ATTTTGTTATAGAGTGGGGTATTTTCGGGAAGAAGGAGGAGGAGTGTATT       BCCCCBCCCCA?:=ACCBCABCCCCCBCCA??5=9@4BB@;??B@BABBA      NM:i:8        MD:Z:1C1C5C9C2C0C22C1C1  XM:Z:.h.x.....x.........h..hh.Z....................h.x. XR:Z:CT XG:Z:CT
SRR020138.15024320_SALK_2029:7:100:1672:1164_length=86  16      chr5    28344365        8       50M     *       0   0         CACAAAATATCAACACCCCTAAACCCCACATTATTCAAAAATCAATTATA       @@@BBBA@A9=A@<?::2:<CB@?=:BBAC??CB@@BBBBC>:ACABCAB      NM:i:11       MD:Z:4G1G1G3G9G9G5G0G0G1T4G2     XM:Z:....x.h.h...x.........h.........h.....hhh......h.. XR:Z:CT XG:Z:GA
SRR020138.15024321_SALK_2029:7:100:1672:433_length=86   0       chr14   38241894        42      50M     *       0   0         TTTTGAGTAGAGAAGTTAGTATTTTAGGGAATTTTTGATTTTTTTAAGTT       BCCBB?B@@A>@-4BBB:7@BBBCBBC@@=A@BCACA;BCBBCBB@@@BB      NM:i:14       MD:Z:0C0C13C0C6C0C6C0C1C3C0C1C0C1C5      XM:Z:hh.............hx......hx......hh.x...hh.hh.h..... XR:Z:CT XG:Z:CT

甲基化call字符串对于BS-read中不涉及胞嘧啶的每个位置都用一个点"."代替，或包含以下三个不同胞嘧啶甲基化的字母(大写=甲基化，小写=未甲基化):

z - C in CpG context - unmethylated #C在CpG里，非甲基化
Z - C in CpG context - methylated #C在CpG里，甲基化
x - C in CHG context - unmethylated #C在CHG里，非甲基化
X - C in CHG context - methylated #C在CHG里，甲基化
h - C in CHH context - unmethylated #C在CHH里，非甲基化
H - C in CHH context - methylated #C在CHH里，甲基化
u - C in Unknown context (CN or CHN) - unmethylated #C在未知区域里，非甲基化
U - C in Unknown context (CN or CHN) - methylated #C在未知区域里，甲基化
. - not a C or irrelevant position #该位置不是胞嘧啶

所以根据这个信息我们就可以知道每一条read的甲基化信息了。

（7）Bowtie 2 或HISAT2模式的局部比对

(在上一个版本中提过:https://github.com/FelixKrueger/Bismark/releases/tag/0.22.0)

Wu等人发表的一篇文章进一步阐述了我们对单细胞的BS-seq(通常称为PBAT文库)可以生成嵌合reads（generate chimeric read pairs），该文章进一步详细描述了片段内嵌合体阻碍单细胞BS-seq文库的高效比对。在那篇文章里，作者描述了一个pipeline，首先使用paired-end比对，然后是单端比对步骤，使用局部比对来提高分子内嵌合体的mapping。为了允许对单细胞或PBAT库进行这种类型的改进，Bismark还允许使用局部对齐。

请注意，我们仍然不建议使用局部比对来大幅增加mapping效率(请参见here)，但是我们承认PBAT / scBSs-seq /scNMT-seq特殊应用局部比对可能的确有所不同。我们还没为局部比对设置更合适或更严格的默认值，也没有调查是否在甲基化提取步骤中需要--ignore标签。

(III) Bismark去重步骤

基本格式：

deduplicate_bismark --bam [options] <filenames>

deduplicate_bismark脚本是为了从Bismark输出中(包括单端和双端的SAM/BAM文件)去除基因组中同一位置的比对，这些重复由于PCR扩增过度而产生。排列在同一基因组位置但在不同链上的序列是单独计分的。

值得注意的是，RRBS不推荐去重、扩增子或其他靶标类型丰富的文库也不建议去重。

在默认模式下，与给定位置的第一次比对将被使用，而不管它的甲基化call。由于比对没有以任何方式进行排序，所以对于每个位置来说，这已经足够接近随机读取了。

对于单端比对，只有染色体、起始坐标和read链将用于去重。

对于双端比对，染色体、read链、第一条read的开始坐标以及第二条read的开始坐标将用于去重。该脚本期望Bismark输出为SAM/BAM格式。

请注意，对于paired-end的BAM文件，去重脚本要求Read 1和Read 2在连续的行中相互跟随!如果文件是按位置排序的，请确保根据read名称重新排序(例如使用samtools sort -n)。

代码：

$ deduplicate_bismark --bam test_data_bismark_bt2.bam

生成文件：

（1）使用UMI或者barcodes去重

在去重的时候，除了染色体，起始位置(和结束位置)和链的方向的选项--barcode，还要考虑潜在的barcodes或UMIs。barcode需要是read ID的最后一个元素，并且必须由一个:分隔，比如：

MISEQ:14:000000000-A55D0:1:1101:18024:2858_1:N:0:CTCCT

（2）对相同文库的多个文件进行去重

当使用选项--multiple时，所有指定的输入文件都被视为单个样本，并连接在一起进行去重。它对SAM文件使用Unix cat，对BAM文件使用samtools cat。
BAM文件的附加说明:尽管这可以在BAM或CRAM上工作，但所有输入文件必须彼此具有相同的格式。每个输入文件的序列字典必须是相同的。默认情况下，header取自要连接的第一个文件。

下一篇将进行甲基化信息提取，以及输出文件的解读。