DNA甲基化测序数据处理(一)

作者: 面面的徐爷 | 来源:发表于2018-04-25 18:17 被阅读136次

    前言

    因为组里面出了一批甲基化测序数据,使用的技术为BS-seq,处理的时候顺带记录了学习过程,演示使用数据为官方提供的example.fastq。

    DNA甲基化及CpG岛定义了解

    DNA甲基化(英语:DNA methylation)为DNA化学修饰的一种形式,能在不改变DNA序列的前提下,改变遗传表现。为外遗传编码(epigenetic code)的一部分,是一种外遗传机制。DNA甲基化过程会使甲基添加到DNA分子上,例如在胞嘧啶环的5'碳上:这种5'方向的DNA甲基化方式可见于所有脊椎动物。
    人类细胞内,大约有1%的DNA碱基受到了甲基化。在成熟体细胞组织中,DNA甲基化一般发生于CpG双核苷酸(CpG dinucleotide)部位;而非CpG甲基化则于胚胎干细胞中较为常见[1] [2]。植物体内胞嘧啶的甲基化则可分为对称的CpG(或CpNpG),或是不对称的CpNpNp形式(C与G是碱基;p是磷酸根;N指的是任意的核苷酸)。
    特定胞嘧碇受甲基化的情形,可利用亚硫酸盐定序(bisulfite sequencing)方式测定。DNA甲基化可能使基因沉默化,进而使其失去功能。此外,也有一些生物体内不存在DNA甲基化作用。
    ——维基百科

    DNA甲基化作为基因组上的表观修饰(区别于组蛋白修饰),存在于各种生物中。

    简单来说,可以认为甲基化代表着基因的失活,去甲基化则标志基因的激活与表达

    虽然CpG序列出现的频率并不高,但是在某些基因区域内,CpG的密度很高,俗称CpG岛。这些CpG岛大多出现在基因的启动子区域(人类占到70%),长度达300-3000bp。目前的研究表明,大多数的管家基因都含有CpG岛,位于基因的5'端(其中的大多数CpG岛都是未甲基化的)。

    另外需要注意的是,目前的研究表明,肿瘤样本与正常样本的CpG岛甲基化差异大多不是发生CpG岛的内部而是位于CpG岛岸(CpG island shore)

    由于CpG位点的易甲基化导致胞嘧啶脱氨变成胸腺嘧啶,所以在漫长的进化过程中,CpG位点逐渐消失,但是又存在着对于基因表达的调控要求,所以CpG岛的出现也被理解为抵抗甲基化经常很,维持调控功能。

    DNA甲基化测序原理与方法

    此处略过,请自行了解(示例文件为WGBS单端测序文件)。

    其实是因为理解起来比较累,知道大致原理就可以了,以后再来填坑

    数据处理(使用Bismark软件处理)

    Bismark下载

    Bismark官网

    The less the people know about how sausages and our code are made, the better they sleep at night (untracable author)
    bismark_v0.19.0.tar.gz
    解压即用tar xzf bismark_v0.X.Y.tar.gz,常用的话将路径写入PATH:PATH=/PATH/TO/Bismark/:$PATH
    manual页面需要搭梯子才能看到,简书不支持附件,所以有需要的可以留言,我email pdf文件

    Dependencies

    需要用户已经装好bowtie1/bowtie2

    测试数据获取

    此处使用测试数据test.fastq
    (from SRR020138, Lister et al., 2009; trimmed to 50 bp; base call qualities are Sanger encoded Phred values (Phred33)).

    $ls -l
    total 2.1M
    -rw-r--r-- 1 sxj users 2.1M Apr 17  2018 test_data.fastq
    

    INDEX文件构建

    # under install folder
    ./bismark_genome_preparation --path_to_bowtie /PATH/to/bowtie2/ --verbose ~/PATH/to/GRCh38/
    

    比对

    # paired-end data
    ./bismark --genome ~/PATH/to/GRCh38/ -1 read1.fastq.gz -2 read2.fastq.gz -p 4 -o ./ 2>test.log
    
    # single-end data
    ./bismark --genome ~/PATH/to/GRCh38/ test_data.fastq -p 4 -o ./ 2>test.log
    

    甲基化位点提取

    ./bismark_methylation_extractor -s --gzip --bedGraph --buffer_size 10G --cytosine_report --comprehensive --genome_folder ~/PATH/to/GRCh38/ test_data_bismark_bt2.bam 2>extracor.log
    

    生成处理报告

    ./bismark2report
    

    所有结果文件

    ls -l
    -rw-r--r-- 1 sxj users  82K Apr 18 16:11 CHG_context_test_data_bismark_bt2.deduplicated.txt.gz
    -rw-r--r-- 1 sxj users 154K Apr 18 16:11 CHH_context_test_data_bismark_bt2.deduplicated.txt.gz
    -rw-r--r-- 1 sxj users 138K Apr 18 16:11 CpG_context_test_data_bismark_bt2.deduplicated.txt.gz
    -rw-r--r-- 1 sxj users 105K Apr 18 17:30 extracor.log
    -rw------- 1 sxj users  20K Apr 17 15:56 nohup.out
    -rw-r--r-- 1 sxj users 450K Apr 17 15:56 test_data_bismark_bt2.bam
    -rw-r--r-- 1 sxj users 449K Apr 18 10:34 test_data_bismark_bt2.deduplicated.bam
    -rw-r--r-- 1 sxj users  48K Apr 18 16:11 test_data_bismark_bt2.deduplicated.bedGraph
    -rw-r--r-- 1 sxj users  55K Apr 18 16:11 test_data_bismark_bt2.deduplicated.bismark.cov
    -rw-r--r-- 1 sxj users 217M Apr 18 17:30 test_data_bismark_bt2.deduplicated.CpG_report.txt.CpG_report.txt.gz
    -rw-r--r-- 1 sxj users 2.9K Apr 18 16:11 test_data_bismark_bt2.deduplicated.M-bias.txt
    -rw-r--r-- 1 sxj users  766 Apr 18 16:11 test_data_bismark_bt2.deduplicated_splitting_report.txt
    -rw-r--r-- 1 sxj users  260 Apr 18 10:34 test_data_bismark_bt2.deduplication_report.txt
    -rw-r--r-- 1 sxj users 347K Apr 19 23:18 test_data_bismark_bt2_SE_report.html
    -rw-r--r-- 1 sxj users 1.8K Apr 17 15:56 test_data_bismark_bt2_SE_report.txt
    -rw-r--r-- 1 sxj users 2.1M Apr 17 15:14 test_data.fastq
    -rw-r--r-- 1 sxj users 6.2K Apr 17 15:56 text.log
    

    结果解读

    --cytosine_report参数会根据当前目录下的信息文件生成一个HTML格式的报告文件,即test_data_bismark_bt2_SE_report.html文件,它包括了比对信息,甲基化信息,M-bias等,可以对数据有一个大概的认知(下图只展示了一部分):

    比对信息 M-bias

    同时因为使用了--comprehensive,所以结果合并正反链的数据后会输出CpG/CHG/CHH三种类型的甲基化文件,包含了胞嘧啶所有的组合形式,但实际上我们自然最关注的是CpG位点的甲基化。其中

    CpG_context_test_data_bismark_bt2.deduplicated.txt即CpG甲基化位点的文件。

    head CpG_context_test_data_bismark_bt2.deduplicated.txt
    
    Bismark methylation extractor version v0.19.0
    SRR020138.15024317_SALK_2029:7:100:1672:902_length=86   -   1   57333019    z
    SRR020138.15024319_SALK_2029:7:100:1672:31_length=86    +   2   10026473    Z
    SRR020138.15024331_SALK_2029:7:100:1673:1282_length=86  +   11  78025243    Z
    SRR020138.15024343_SALK_2029:7:100:1673:202_length=86   +   10  121617231   Z
    SRR020138.15024357_SALK_2029:7:100:1673:879_length=86   -   4   75173715    z
    SRR020138.15024361_SALK_2029:7:100:1673:235_length=86   -   2   130768889   z
    SRR020138.15024368_SALK_2029:7:100:1673:123_length=86   +   10  106402850   Z
    SRR020138.15024376_SALK_2029:7:100:1673:1908_length=86  -   12  124097382   z
    SRR020138.15024380_SALK_2029:7:100:1673:397_length=86   +   8   95038627    Z
    # 第一列为测序信息
    # 第二列为甲基化状态 + 代表甲基化 -代表未甲基化
    # 第三列代表chromosome
    # 第四列代表location
    # 第五列代表methylation call,简单来说大写的就是甲基化的(因为还有CHG,CHH的数据,分别对应x, X , h, H)
    

    test_data_bismark_bt2.deduplicated.bismark.cov文件则给了每个位点的甲基化比例,为下一步确定CpG岛提供了基础,其数据形式如下:

    $head test_data_bismark_bt2.deduplicated.bismark.cov
    1   975476  975476  100 1   0
    1   975488  975488  100 1   0
    1   975490  975490  100 1   0
    1   2224487 2224487 100 1   0
    1   2224489 2224489 100 1   0
    1   2224514 2224514 100 1   0
    1   2224520 2224520 100 1   0
    1   2313220 2313220 0   0   1
    1   9313902 9313902 100 1   0
    1   9313914 9313914 100 1   0
    
    # 第一列代表chromosome
    # 第二,三列代表location
    # 第四列代表甲基化百分比
    # 第五列代表甲基化数目
    # 第六列代表未甲基化数目
    

    test_data_bismark_bt2.deduplicated.CpG_report.txt.CpG_report.txt文件则是背景信息:

    $head test_data_bismark_bt2.deduplicated.CpG_report.txt.CpG_report.txt
    1   10469   +   0   0   CG  CGC
    1   10470   -   0   0   CG  CGA
    1   10471   +   0   0   CG  CGG
    1   10472   -   0   0   CG  CGC
    1   10484   +   0   0   CG  CGG
    1   10485   -   0   0   CG  CGG
    1   10489   +   0   0   CG  CGC
    1   10490   -   0   0   CG  CGG
    1   10493   +   0   0   CG  CGC
    1   10494   -   0   0   CG  CGG
    
    # 第一列为chromosome
    # 第二列为location
    # 第三列为strand
    # 第四,五列为甲基化和非甲基化的碱基数目
    # 第六列为CG
    # 第七列为背景(第一个C延伸两个碱基)
    

    其它参数

    bam2nuc
    bismark2summary
    coverage2cytosine
    NOMe_filtering
    filter_non_conversion
    # 有需要的可以自信探索 --help or manual
    

    结语

    此处根据测序数据得到了甲基化位点的信息,但是后续DML以及DMR的确定还需要R包的使用,以及后续的可视化还以探索以下包:


    挖坑待填

    Methy-Pipe: An Integrated Bioinformatics Pipeline for Whole Genome Bisulfite Sequencing Data Analysis
    2014年出的一个pipeline,分析作图一条龙,有兴趣的同学安排一下。

    相关文章

      网友评论

      • elaine0622:多谢分享,我最近也在做这个,有些地方不是特别明白,看到你的评论说“deduplication是要看是wgbs或者rrbs的”,这两种有什么区别吗?我的数据是rrbs
        面面的徐爷:@elaine0622 你好,你需要先从原理理解wgbs和rrbs的区别,然后再理解duplication的作用就没有什么问题了,相信你可以轻松解决。
      • 热衷组培的二货潜:Batmeth2,methylkit,bsseq,DSS,DMRcaller都还可以
        面面的徐爷:@二货潜 谢谢你的介绍,目前还在比较哪种方法最好,methykit DDS都在了解,下一篇文章会介绍关于这些包的学习。
      • 热衷组培的二货潜:deduplication那步忘了加吗?😂希望能多点可视化方面的,特别是那个基因整体上下游的分布等
        面面的徐爷:@二货潜 你好,deduplication是要看是wgbs或者rrbs的,所以这一步是可选择性的。可视化的一些包我还在学习,但是展开来讲就太大了……

      本文标题:DNA甲基化测序数据处理(一)

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