DNA甲基化测序数据处理（二）：差异分析
衔接上一篇数据比对后的结果，使用R包DSS进行处理。 我们先来复习一下上一节课得到的数据结果： *.bismark.cov.gz 文件 这里我们使用的R包为DSS，使用 Bioconductor 进行安装。 这个包可以对甲基化数据做两件事： DSS包对差异甲基化的检测基于β负二项分布的严格沃尔德检验。 输入数据的格式如下: DML是甲基化差异位点，DMR为甲基化差异区域。 使用DSS包自带的数据进行演示 1. 载入两个包DSS和bsseq（需要该包构建obj对象） 2. 利用DMLtest函数call DML，分为一下几个步骤： 根据甲基化水平进行loci的差异分析 3. 根据dmlTest来call DML 当然，用户也可以指定差异的阈值，只有差异大于阈值的才会被call出来。 4. 根据dmlTest Call DMR 我们可以发现，smoothing前后得到的结果差异还是很大，所以针对不同的实验类型我们需要注意是否使用smoothing。 同理，也可以使用的delta参数以及调整p.threshold得到合适的结果。 5. 可视化 DSS包提供了一个不是很美观的可视化函数，用户其实可以使用coverage结果在R里面作图。 分析到此就告一段落了，随后就是自行对差异甲基化区域的注释以及可视化文献作图。 后续在处理数据的时候，发现有一个情况，多核跑DMLtest会报错，详情解决方案可见 https://github.com/haowulab/DSS/issues/13 以下为简单解决方案
甲基化不同测序环境cg、chg、chh是什么意思？
甲基化C碱基在基因组上的分布包含三种形式（CG, CHG和CHH，其中H代表A 或T或C碱基）。就这个意思在植物中的甲基化分三种类型：CG，CHG，CHH，其中C是加上甲基的胞嘧啶，H是任何除了G之外的核苷酸残基。这几种甲基化的不同在于它们建立，维护以及遗传的原理；而且它们在基因组调节过程中的作用不同，这样就会对机体的特性—或者称为表型产生不同的影响。以前在几种不同植物中通过对全基因组范围的甲基化特性分析确认了许多可变的DNA甲基化区域。限制性内切酶 它能识别特定的核苷酸序列 并在特定的切点上切割DNA分子。一种限制性内切酶只能识别一种特定的核苷酸序列。
甲基化的甲基化检测方法
DNA甲基化是最早发现的基因表观修饰方式之一，真核生物中的甲基化仅发生于胞嘧啶，即在DNA甲基化转移酶（DNMTs）的作用下使CpG二核苷酸5’-端的胞嘧啶转变为5’-甲基胞嘧啶。DNA甲基化通常抑制基因表达，去甲基化则诱导了基因的重新活化和表达。这种DNA修饰方式在不改变基因序列前提下实现对基因表达的调控。脊椎动物DNA的甲基化状态与生长发育调控密切相关，比如在肿瘤发生时，抑癌基因CpG岛以外的CpG序列非甲基化程度增加，CpG岛中的CpG则呈高度甲基化状态，导致抑癌基因表达的下降。1、甲基化特异性的PCR（Methylation-specific PCR，MSP）用亚硫酸氢盐处理基因组DNA，所有未发生甲基化的胞嘧啶被转化为尿嘧啶，而甲基化的胞嘧啶不变；随后设计针对甲基化和非甲基化序列的引物进行PCR。通过电泳检测MSP扩增产物，如果用针对处理后甲基化DNA链的引物能得到扩增片段，则说明该位点存在甲基化；反之，说明被检测的位点不存在甲基化。2、亚硫酸氢盐测序法（Bisulfite sequencing PCR，BSP）用亚硫酸氢盐处理基因组DNA，则未发生甲基化的胞嘧啶被转化为尿嘧啶，而甲基化的胞嘧啶不变。随后设计BSP引物进行PCR，在扩增过程中尿嘧啶全部转化为胸腺嘧啶，最后对PCR产物进行测序就可以判断CpG位点是否发生甲基化称为BSP-直接测序方法。将PCR产物克隆至载体后进行测序，可以提高测序成功率，这种方法称为BSP-克隆测序法。3、高分辨率熔解曲线法（High Resolution Melting，HRM）在非CpG岛位置设计一对针对亚硫酸氢盐修饰后的DNA双链的引物，这对引物中间的片段包含感兴趣的CpG岛。若这些CpG岛发生了甲基化，用亚硫酸氢盐处理后，未甲基化的胞嘧啶经PCR扩增后转变成胸腺嘧啶，而甲基化的胞嘧啶不变，样品中的GC含量发生改变，从而导致熔解温度的变化（图1）。其中，样品要求：细胞（≥106 个）、组织（≥300mg）、血液（≥1ml）、血清（≥1.5ml）等样品材料，基因组DNA（体积≥20μl，浓度≥50 ng/μl）。
甲基化的甲基化检测技术
曾经是2009年最值得关注的技术。遗传上除了ATGC这四种碱基，人们对第五种碱基-甲基化的胞嘧啶的兴趣也日益增加。甲基化修饰的存在对DNA转录的调控起了重要作用，异常的甲基化往往是许多疾病的起因。既然如此重要，系统绘制甲基化组的方法自然也多了起来。甲基化组（methylome）这个词还不太常见，它指的是全基因组范围内的甲基胞嘧啶。美国Salk生物研究院的Joseph Ecker及其同事刚刚通过高通量测序的方法，展现了一张人胚胎干细胞中所有甲基胞嘧啶的完整图谱。这是第一张单碱基分辨率的哺乳动物甲基化图谱，并伴随着mRNA和小RNA以及一些组蛋白修饰的比较分析。他们发现胚胎干细胞中近四分之一的甲基化是在非CG背景下，暗示胚胎干细胞采用不同的甲基化机制来影响基因调控。随着ES细胞的诱导分化，非CG甲基化消失，而在iPSC中还原。这张图谱为未来人类疾病和发育中的表观遗传学修饰研究打下了基础。美国Whitehead研究院的Meissner等也曾绘制了类似的图谱。他们利用高通量的亚硫酸氢盐测序和单分子测序，产生了覆盖大部分CpG岛的DNA甲基化图谱。他们发现，DNA甲基化模式与组蛋白甲基化模式的相关性比基因组序列更好；CpG的甲基化是动态的表观遗传标志，在细胞分化期间经历大量的改变。他们还发现，胚胎干细胞和原代细胞中与发育调控相关的“弱”CpG岛在体外增殖期间经历了异常的超甲基化，让人想起某些原发肿瘤。另外，两个独立的研究小组，分别为哈佛大学的George Church等，以及加州大学的Kun Zhang连同弗吉尼亚联邦大学的Yuan Gao等，也将传统的甲基化工具如DNA的重亚硫酸盐转化与目标基因组捕获技术和高通量测序相结合，定量测定人基因组中的甲基化。尽管这些甲基化图谱的绘制方法略有不同，但他们都采用了亚硫酸氢盐转化，将未甲基化的胞嘧啶转化成尿嘧啶，并在随后的扩增步骤中转化成胸腺嘧啶。虽然很有效，但这种方法需要一些手工操作来确保完全的转化，并通过计算分析来绘制图谱。在2012年2月份，英国Oxford Nanopore Technologies的科学家提出，单分子纳米孔测序可替代这种劳动密集型的技术，测序仪能直接分辨出未修饰的胞嘧啶和甲基化胞嘧啶。当核酸外切酶消化单链DNA后，单个碱基落入孔中，它们瞬间与环式糊精相互作用，并阻碍了穿过孔中的电流。每个碱基ATGC以及甲基胞嘧啶都有自己特有的电流振幅，因此很容易转化成DNA序列。这样，纳米孔技术就能直接读出这第五种碱基。由于纳米孔测序仅限于短的寡核苷酸，因此全基因组测序还有一段很长的路要走。此外，还有一些技术障碍需要克服，比如确保碱基以正确的次序进入纳米孔，然后从另一侧离开。不过，一旦成功，第五种碱基的直接测序将会产生重大影响。
DNA甲基化测序数据处理（一）：数据比对
因为组里面出了一批甲基化测序数据，使用的技术为BS-seq，处理的时候顺带记录了学习过程，演示使用数据为官方提供的example.fastq。 DNA甲基化作为基因组上的表观修饰（区别于组蛋白修饰），存在于各种生物中。 虽然CpG序列出现的频率并不高，但是在某些基因区域内，CpG的密度很高，俗称CpG岛。这些CpG岛大多出现在基因的启动子区域（人类占到70%），长度达300-3000bp。目前的研究表明，大多数的管家基因都含有CpG岛，位于基因的5'端（其中的大多数CpG岛都是未甲基化的）。 另外需要注意的是，目前的研究表明， 肿瘤样本 与正常样本的CpG岛甲基化差异大多不是发生CpG岛的内部而是位于 CpG岛岸（CpG island shore） 。 由于CpG位点的易甲基化导致胞嘧啶脱氨变成胸腺嘧啶，所以在漫长的进化过程中，CpG位点逐渐消失，但是又存在着对于基因表达的调控要求，所以CpG岛的出现也被理解为抵抗甲基化经常很，维持调控功能。 此处略过，请自行了解(示例文件为WGBS单端测序文件)。 Bismark官网 需要用户已经装好bowtie1/bowtie2 此处使用测试数据 test.fastq （from SRR020138, Lister et al., 2009; trimmed to 50 bp; base call qualities are Sanger encoded Phred values (Phred33)). --cytosine_report 参数会根据当前目录下的信息文件生成一个HTML格式的报告文件，即 test_data_bismark_bt2_SE_report.html 文件，它包括了比对信息，甲基化信息，M-bias等，可以对数据有一个大概的认知（下图只展示了一部分）： 同时因为使用了 --comprehensive ，所以结果合并正反链的数据后会输出CpG/CHG/CHH三种类型的甲基化文件，包含了胞嘧啶所有的组合形式，但实际上我们自然最关注的是CpG位点的甲基化。其中 CpG_context_test_data_bismark_bt2.deduplicated.txt 即CpG甲基化位点的文件。 test_data_bismark_bt2.deduplicated.bismark.cov 文件则给了每个位点的甲基化比例，为下一步确定CpG岛提供了基础，其数据形式如下： test_data_bismark_bt2.deduplicated.CpG_report.txt.CpG_report.txt 文件则是背景信息： 此处根据测序数据得到了甲基化位点的信息，但是后续DML以及DMR的确定还需要R包的使用，以及后续的可视化还以探索以下包：
