查看原文
其他

Bismark软件使用入门

生信阿拉丁 生信阿拉丁 2022-05-16

Bismark软件使用入门



01


甲基化文库建库原理

借助重亚硫酸盐(Bisulfite)转换的方法被认为是金标准,该化学试剂会将未甲基化的C进行转化,胞嘧啶(C)转换为尿嘧啶(U),后续的PCR扩增过程U会被胸腺嘧啶(T)替代,而,甲基化的C则不受影响。在了解如何使用Bismark之前,需要先了解DNA甲基化文库的结构。

通常,甲基化文库构建方法有两种:
  • 一种是先片段化,然后末修,加上特殊处理的接头(甲基化的)后,进行CT转化,之后进行扩增,如下图:

该文库为链特异性文库,使用bismark默认参数比对即可。

  • 另一种是先CT转换,然后扩增加接头,如下图。


该文库为非链特异性文库,需要使用PBAT文库参数。

通常而言,甲基化建库均为链特异性方式构建:

  • 测序read1都是来自于原始的基因组,经历了C->T转化,来自于正链的叫OT,负链叫OB;

  • 测序read2都是来自于原始基因组的互补链,经历了G->A转换,称之为CTOT或者CTOB。


02


bismark比对分析

bismark的分析分为三步:

1

基因组索引构建:需要将基因组进行相应的转换,模拟bisulfite处理后基因组,这样才能用于比对。在每个新物种比对之前,都需要进行该处理。bismark内部默认调用bowite软件进行比对。


2

序列比对:输入为下机序列、基因组和比对参数,bismark会输出比对结果和甲基化检测的结果,默认格式为BAM,同时还会有一个统计报告。


3

甲基化位点提取:这一步是可选的,会利用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的统计文件。


03


提取甲基化位点


使用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)即可。


04


bismark的优缺点


4.1、优点

(1)比对和mC Calling比较准确,受众广泛,使用率很高

(2)软件维护性好,更新快,bug修复及时

(3)软件成体系,集成了比对、去重、mC Calling分析

4.2、缺点

(1)比对和mC Calling速度慢,真的很慢,可以通过拆分输入数据为小份作为输出和增大计算资源解决。

(2)bismak相匹配的直接进行下游甲基化水平、图形可视化等分析的套件较少,需要自行计算或进行格式转化作为第三方软件的输入。

(3)说明文档不够系统和详细,在”外行“看来理解从建库到mC Calling有一定的挑战。

bismark的结果只是研究甲基化信息的第一步也是关键的一步。未来,期望bismark在速度上,和在对甲基化水平层面分析的支持上多一下改进。


05


参考资料

1.https://github.com/FelixKrueger/Bismark

2.https://github.com/FelixKrueger/Bismark/tree/master/Docs

3.https://www.epigenesys.eu/images/stories/protocols/pdf/20120720103700_p57.pdf


作者:马可菠萝 童蒙

编辑:angelica


往期精彩回顾

转座子分类软件deepTE简介

单细胞ATAC高级分析

Motif可视化——从PFM矩阵到sequence logo

pycharm使用入门

CopyKat——基于高通量单细胞测序方法鉴定肿瘤细胞拷贝数变异和亚克隆结构

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存