分析DNA甲基化的手段有很多,除了甲基化芯片外,还有WGBS和RRBS等实验与高通量测序相结合的手段,不管是哪种策略,都需要对DNA进行亚硫酸氢盐处理。

在​​WGBS​​​和​​RRBS​​​的分析过程中,最重要就是通过比对识别甲基化位点。原始序列由于经过了亚硫酸氢盐处理,在比对时,为了正确区分甲基化位点和SNP位点,就需要特殊的比对软件。​​bismark​​ 就是甲基化测序常用的比对软件之一。

​bismark​​ 用于将亚硫酸氢盐处理的reads与参考基因组进行比对,并识别甲基化位点。该软件具有以下特点:

  1. 一个步骤同时完成比对和methylcation calling;
  2. 支持单端和双端数据
  3. 支持gapped 和 ungapped 比对
  4. seed length 和 错配个数都是允许调整的
  5. 可以识别CpG, CHG, CHH 等甲基化位点

软件官网

​https://www.bioinformatics.babraham.ac.uk/projects/bismark/​

软件安装

wget https://github.com/FelixKrueger/Bismark/archive/0.19.0.tar.gz
tar xzvf 0.19.0.tar.gz
cd Bismark-0.19.0/
export PATH=$PWD:$PATH

​bismark​​首先对基因组进行两种转换

C->T

将原始序列中的C全部变成T

> chr1
ACAGTACTGGCGGATTATAGGGAAACACCCGGAGCATATGCTGTTTGGTC
> CT
ATAGTATTGGTGGATTATAGGGAAATATTTGGAGTATATGTTGTTTGGTT

G->A

将原始序列中的G全部变成A

> chr1
ACAGTACTGGCGGATTATAGGGAAACACCCGGAGCATATGCTGTTTGGTC
> GA ACAATACTAACAAATTATAAAAAAACACCCAAAACATATACTATTTAATC

具体的转换过程如下

my $sequence = uc$_;
my $CT_sequence = $sequence;
my $CT_transliterations_performed = ($CT_sequence =~ tr/C/T/);
my $GA_sequence = $sequence;
my $GA_transliterations_performed = ($GA_sequence =~ tr/G/A/);

然后调用bowtie后者bowite2 对转换后的序列建立索引,bowtie1和bowtie2最大的区别就是支不支持gapped  align, 只有bowite2支持,也就是说bowtei2比对时支持插入缺失,bowite1不支持。对于目前主流的测序平台产生的数据,一般选择bowtie2就可以了。

对基因组建立索引

将基因组fasta 文件放在1个目录下 ,后缀为​​.fa​​​ 或者 ​​.fasta​​, 然后运行以下命令

bismark_genome_preparation hg19/hg19.fasta

运行成功后,会在fasta 所在目录生成如下的目录结构

├── 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

可以看到,分别对CT和GA转换的基因组建立了bowtie2的索引。