bismark对参考基因组建立索引
分析DNA甲基化的手段有很多，除了甲基化芯片外，还有WGBS和RRBS等实验与高通量测序相结合的手段，不管是哪种策略，都需要对DNA进行亚硫酸氢盐处理。
在WGBS和RRBS的分析过程中，最重要就是通过比对识别甲基化位点。原始序列由于经过了亚硫酸氢盐处理，在比对时，为了正确区分甲基化位点和SNP位点，就需要特殊的比对软件。bismark 就是甲基化测序常用的比对软件之一。
bismark 用于将亚硫酸氢盐处理的reads与参考基因组进行比对，并识别甲基化位点。该软件具有以下特点：
- 一个步骤同时完成比对和methylcation calling;
- 支持单端和双端数据
- 支持gapped 和 ungapped 比对
- seed length 和 错配个数都是允许调整的
- 可以识别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的索引。