美文网首页BS-seq生信分析流程生信必备生物知识
DNA甲基化实战分析-----bismark 代码篇

DNA甲基化实战分析-----bismark 代码篇

作者: liu_ll | 来源:发表于2019-01-24 17:04 被阅读29次

      在上次的文章中简单的介绍了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

    相关文章

      网友评论

        本文标题:DNA甲基化实战分析-----bismark 代码篇

        本文链接:https://www.haomeiwen.com/subject/bjryjqtx.html