DNA甲基化实战分析-----bismark 代码篇
在上次的文章中简单的介绍了bismark的基本知识DNA甲基化分析----------甲基化比对软件专题（Bismark），bismark的功能非常强大，分析甲基化的数据非常好用。今天我们来看看如何用bismark来分析甲基化的数据。
----------------------------------------------------分割线----------------------------------------------------
1：我们都知道bismark的首先是要解决mapping的问题，那么mapping肯定需要设计到序列比对的问题，那么设计到了文件比较基因序列比对的时候，第一步必要做的就是：建立 index
Q:问题来了，我们都知道bismark 可以调用2中比对的模式，一种是bowtie一种是bowtie2，那么我们用哪个呢？
A:用bowtie2比较适用，因为bowtie2可以支持插入缺失，bowite1不支持。对于目前主流的测序平台产生的数据，一般选择bowtie2。
bismark_genome_preparation --bowtie2 /data/genomes/homo_sapiens/GRCh37/
------bismark_genome_preparation 是建立index的命令
------bowtie2 调用bowtie2来建立index
------/data/genomes/homo_sapiens/hg19/ 参考基因的路径，该路径包含了ref的文件
(ps: 我的是已经安装好了bowtie2,bowtie,bismark等，并且加入了bashrc
给出官方文档的教程代码：
bismark_genome_preparation --path_to_bowtie /usr/local/bowtie/ --verbose
/data/genomes/homo_sapiens/GRCh37/ )
--verbose 输出log信息
建立的索引文件信息
2： 建立完索引信息了，接下来开始比对。
在这里需要特别的注意一下！！！！如果是利用bowtie2建立的index，那么在比对的时候也是需要利用bowtie2来进行比对的，不能用bowtie！
cd /liull/test/bs-seq-test/4_TZctDNAme_in_action/3_bismark
bismark_path= /liull/test/bs-seq-test/4_TZctDNAme_in_action/3_bismark
mkdir ./reference
mkdir ./lamdaDNA
bismark --bowtie2 -N 0 -L 20 --quiet --un --ambiguous --sam \
-o ${bismark_path}/reference \
/liull/Database/hg19/ \
-1 /liull/test/bs-seq-test/4_TZctDNAme_in_action/2_QC/after_QC/TZCTDNA40_H7HLCCCXY_L2_1_trimmed.fq.gz \
-2 /liull/test/bs-seq-test/4_TZctDNAme_in_action/2_QC/after_QC/TZCTDNA40_H7HLCCCXY_L2_2_trimmed.fq.gz
ps:
1：这里分别将基因组和参考基因组和lambdaDNA进行了比较，因为在建库的时候，为了得到BS转化率掺入了一段lambdaDNA序列，这段序列是完全已知的，所以通过lambdaDNA的转化率可以得到这次BS的转化率。
2：--bowtie2 利用bowtie2进行的操作
3：-N 0 在比对的时候允许0个错配信息。
4：-L 20 是在比对的时候建立的seed的大小是20
5： --quiet 是说不输出在比对的时候的比对流程信息
6：--ambiguous 是说如果有一个序列匹配到了多个地方把这个序列记录下来，它会保存在 _ambiguous_reads_1.txt/_ambiguous_reads_2.txt.fq
7：-un 是说这些多处匹配的reads信息不会写在输出的fq中
8：--sam 指输出的格式为sam
3：接下来对比对后的sam文件进行去重复，用到的是deduplicate_bismark
deduplicate_bismark \
-p \
/liull/test/bs-seq-test/4_TZctDNAme_in_action/3_bismark/Me-1/reference/Me-1_1_head100000.fastq_bismark_bt2_pe.sam
deduplicate_bismark \
-p \
/liull/test/bs-seq-test/4_TZctDNAme_in_action/3_bismark/Me-1/lamdaDNA/Me-1_1_head100000.fastq_bismark_bt2_pe.sam
-p：对 paired-end 的 bismark 输出文件 (default format: SAM) 去重复
4：去掉了测序带来的重复片段，接下来我们看一下如何去解读甲基化信息，这里用到的是bismark_methylation_extractor
bismark_methylation_extractor \
-p \
--comprehensive \
--no_overlap \
--bedGraph \
--counts \
--buffer_size 20G \
--report \
--cytosine_report \
--genome_folder ${ref_PATH}/hg.fa \
${total_PATH}/3_bismark/Me-1/reference/Me-1_1_head100000.fastq_bismark_bt2_pe.deduplicated.sam \
-o ${total_PATH}/3_bismark/Me-1/reference
bismark_methylation_extractor \
-p \
--comprehensive \
--no_overlap \
--bedGraph \
--counts \
--buffer_size 20G \
--report \
--cytosine_report \
--genome_folder ${lamdaDNA_path}/lambda_dna.fa \
${total_PATH}/3_bismark/Me-1/lamdaDNA/Me-1_1_head100000.fastq_bismark_bt2_pe.deduplicated.sam \
-o ${total_PATH}/3_bismark/Me-1/lamdaDNA
ps:
-p : 输入的文件是pair-end
--comprehensive 将可能的甲基化信息都输出，包括CHG,CHH,CpG
--no_overlap：对于双端读取，read_1和read_2理论上是可能重叠的。这个选项避免了重复的甲基化计算了2次。虽然这消除了对序列片段中心更多甲基化的计算偏差。
--bedGraph:输出bedGraph文件
--cytosine_report：指报道全基因组所有的CpG。
--counts:指在bedGraph中有每个C上甲基化reads和非甲基化reads的数目(在--cytosine_report指定的情况下。)
--buffer_size 指定的内存去计算甲基化信息
--report ：会有一个甲基化的summary
--genome_folder:后跟着参考基因组的位置
-o:输出的文件路径
5：解读甲基化信息
cd /share_bio/unisvx1/sunyl_group/liull/test/bs-seq-test/4_TZctDNAme_in_action/3_bismark/Me-1/reference
bismark2report
cd /share_bio/unisvx1/sunyl_group/liull/test/bs-seq-test/4_TZctDNAme_in_action/3_bismark/Me-1/lamdaDNA
bismark2report
ps：直接利用bismark2report这个命令可以看到报告
6：unix_sort
对生成的sam文件进行排序，然后方便methylkit进行可视化
#####6: unix_sort
grep -v '^[[:space:]]*@' /share_bio/unisvx1/sunyl_group/liull/test/bs-seq-test/4_TZctDNAme_in_action/3_bismark/Me-1/reference/Me-1_1_head100000.fastq_bismark_bt2_pe.deduplicated.sam | \
sort -k3,3 -k4,4n | \
perl -F"\t" -lane 'print if $F[2]=~ /^(?:\d|X|Y)\d*$/' >/share_bio/unisvx1/sunyl_group/liull/test/bs-seq-test/4_TZctDNAme_in_action/3_bismark/Me-1/reference/Me-1_1_head100000.fastq_bismark_bt2_pe.deduplicated.sorted.chr1-22-XY.sam
-------------------------------------------------分割线--------------------------------------------------------
后续的可视化留到下篇的笔记，这篇主要是利用bismark来分析甲基化的数据。走一遍分析的流程，后续的统计分析没有介绍，有兴趣的朋友可以关注methykit等R包。这里先做一个小预告。
(PS:大神@二货潜 之前写过一篇踩坑的小文章，有兴趣的可以看看：bismark_methylation_extractor提取甲基化信号的问题)
Reference:
1:http://www.bioinformatics.babraham.ac.uk/projects/bismark/
2:bismark官方文档
3:http://www.bioinformatics.babraham.ac.uk/training/Methylation_Course/BS-Seq%20theory%20and%20QC%20lecture.pptx
