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