导读DNA甲基化测序数据处理（二）：差异分析答衔接上一篇数据比对后的结果，使用R包DSS进行处理。我们先来复习一下上一节课得到的数据结果：*.bismark.cov.gz 文件这里我们使用的R包为DSS，...

今天运困体育就给我们广大朋友来聊聊山西甲基化分析,希望能帮助到您找到想要的答案。

DNA甲基化测序数据处理（二）：差异分析

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会报错，详情解决方案可见

以下为简单解决方案

甲基化套路

答和甲基化有关的。

可以先了解下甲基化：

450k甲基化基础

450K甲基化芯片数据处理传送门

450k甲基化芯片常用工具包：ChAMP和minfi等。

甲基化的一些预备知识

甲基化程度的量化

DMP（或DML，差异甲基化位点）与 DMR（差异甲基化区域）的关系。如何定义DMR？

一般来说，DMR是通过统计bump来计算出来的，可以参考： ChAMP 分析甲基化芯片数据-差异分析下篇

一般来说，我们还会关注两个方面的信息：DMR与CpG岛的关系，DMR与基因的关系。

DMR与CpG岛的关系：图片来自 ShengXinRen

关于DMR或DMP与基因的关系（笔者特别关注甲基化位点的功能注释），简要总结如下。

一般而言，启动子区域的甲基化程度影响基因的转录（但也有报道说第一外显子等位置的甲基化也与基因的转录相关）。如何描述一个基因的转录相关的甲基化程度呢？

有个论坛上是这么说的（ ShengXinRen ）：

也有人总结如下：

以及这种说法（ ）：

另外，补充一个知识（ “启动子预测”技能 ）：

感觉暂时并没有统一的标准。

可以自己尝试各种界定标准：

1、 TSS上游1500bp、2000bp、5000bp内的甲基化位点的平均值；

2、 TSS上游及下游1500bp、2000bp、5000bp内的甲基化位点的平均值；

3、 TSS上游1500bp、2000bp、5000bp内或5-UTR或第一外显子的甲基化位点平均值。

考虑5000是因为CpG岛的加上两边的Shore一般可达到6kb左右。注意到，5-UTR是第一外显子的一部分。有时候甚至还可以加上是否为CpG岛（或Shore）这个限定。

下图来自： 彻底搞清楚promoter, exon, intron, and UTR

3.4分，简单的肠癌甲基化分析。主要涉及差异分析、关联分析、功能注释。

解读： 如何从甲基化入手，轻松整篇预后标志物的文章

1、 数据质控 ：共485,577个基因座的DNA甲基化数据，在预处理数据和质量控制后，保留了467,971个探针。

2、 差异筛选 ：minfi包筛选DMR（差异化甲基化区域），这一步类似于RNA-Seq的筛选差异基因。结果：最终得到675个差异甲基化区域，其中654个上调。

3、 注释和功能

3-1 DMR的注释 ：这些DMR区域与基因的关系是什么呢？我们利用这些差异甲基化区域的位置与基因的各个元件位置的关系，观察这些差异甲基化区域主要分布在基因的哪些位置上。结果：上调的甲基化区域大多数位于基因的第一外显子，5'UTR，TSS200，TSS150和基因体中，而只有少数UMR位于基因间和3'UTR中区域，同样的下调的甲基化区域也有相同的现象。

3-2 DMR与CpG岛的关系 ：差异甲基化区域与CpG岛的关系如图，从中可以看出上调的差异甲基化区域主要聚集在CpG岛区域，而下调的差异甲基化区域主要聚集在低CpG岛密度区域。

总结：

在本研究中，在大量COAD样品中进行了DNA甲基化谱的综合分析，以研究COAD中存在的改变的DNA甲基化模式。COAD样品和邻近组织样品之间的DNA甲基化谱的比较揭示了COAD样品中异常的DNA甲基化变化，并导致675个DMR的鉴定，包括654个高甲基化和21个低甲基化DMR。这些结果与先前的研究结果一致，即DNA高甲基化是结直肠癌的常见特征。

此外，这些DMR可用于有效区分COAD样品和相邻组织样品，这表明DMR可能在COAD的形成中具有致病作用。基因组分析显示，DMR主要位于启动子区域（包括第1 外显子，5'UTR和TSS）和体区，这与之前在其他类型癌症中的观察结果一致。在基因间和 3'UTR 区域中仅发现了一小部分DMR。此外，大多数高甲基化DMR位于CpG岛中，而大多数低甲基化DMR不位于CpG岛或注释基因中。

DNA甲基化数据分析(二)

答DSS是一个R包，对基于计数的测序数据进行差异分析。 它可以检测来自RNA-seq的差异表达基因（DEG），以及来自亚硫酸氢盐测序（BS-seq）的差异甲基化位点或区域（DML / DMR）。

DSS (Dispersion Shrinkage for Sequencing data)，为基于高通量测序数据的差异分析而设计的Bioconductor包。主要应用于BS-seq（亚硫酸氢盐测序）中计算不同组别间差异甲基化位点（DML）和差异甲基化区域（DMR）即Call DML or DMR。

使用R包DSS多种方式检验差异甲基化信号区域。

这里我使用的R版本为3.6.0,安装“DSS”包

1.准备输入文件

DSS包要求输入数据格式：每一行代表一个CpG位点，格式如下：

第一列为染色体

第二列为位置

第三列为总reads数

第四列为甲基化的reads数

我们看一下上一步我们得到的数据test_R1_bismark_bt2_pe.deduplicated.bismark.cov.gz

在file.cov.gz中

第一列数据为染色体

第二、三列数据为甲基化位置

第四列为甲基化百分比

第五列为甲基化数目

第六列为未甲基化数目

在这里,需要对我们的数据利用R，linux或者python整理成DSS包所需输入文件格式。这里我使用python整理的输入文件格式。

2.对不同组别之间DNA甲基化进行差异分析

这里我们使用DSS包自带的数据

2.1 加载DSS包，并读取甲基化数据

2.2构建BSobj对象

可以看出我们的数据包括4个样本，34739个甲基化位点。

2.3利用DMLtest函数call DML

For whole-genome BS-seq data, perform DML test with smoothing

smoothing: 用于指示是否在估计平均甲基化水平时应用平滑的标志,这里我理解的smoothing有点类似滑动窗口的意思。我们也可以选择smoothing=False

2.4利用callDML函数可以找出差异甲基化位点

delta:定义DML的阈值。在DML检验程序中，在每个CpG位点进行两组均值相等的假设检验。这里如果指定了delta，函数将计算均值差大于delta的后验概率，然后基于此调用DML。

p.threshold: 当未指定delta时，这是定义DML的p值阈值，例如p值小于该阈值的位点将被视为DML。当指定delta时，后验概率大于1-p阈值的CpG位点被视为DML。

2.5利用callDMR函数可以找出差异甲基化区域

当不同组别间CpG位点区域具有显著的统计学差异时这段差异区域被定义为DMRs。DMR也是基于DML被检测出来的。

如果想把甲基化位点和甲基化区域信息输出，可通过write.csv()输出。

2.6 #使用showOneDMR函数可视化DMR

给定一个DMR和一个BSseq对象，此函数将生成一个多面板图，每个图对应一个示例文件，以可视化甲基化水平。每个CpG上都有一个条，灰色线条表示总覆盖率，蓝色线条表示甲基化水平。

showOneDMR唯一缺陷只能展示一个区域甲基化水平，有点丑，感兴趣的同学可以自己利用脚本画一下全基因组甲基化水平。

参考：

1.

methylkit差异甲基化分析

答典型的文件应该是这样的

一般呈现两边高中间低的特点。要么高甲基化，要么低甲基化。

如果数据受到了PCR扩增的影响，那么右边还会出现一个峰。

得到在所有样本中都有cover的甲基化位点or甲基化区域（这里是甲基化区域）。

得到相关系数图

We can also do a PCA analysis on our samples. The following function will plot a scree plot for importance of components.

We can also plot PC1 and PC2 axis and a scatter plot of our samples on those axis which will reveal how they cluster.

The calculateDiffMeth() function is the main function to calculate differential methylation. Depending on the sample size per each set it will either use Fisher’s exact or logistic regression to calculate P-values. P-values will be adjusted to Q-values using SLIM method (Wang, Tuominen, and Tsai 2011). If you have replicates, the function will automatically use logistic regression.

calculateDiffMeth()函数是计算差异甲基化的主要函数。根据每个组的样本多少，它将使用Fisher精确检验或逻辑回归检验来计算p值。使用SLIM方法将p值调整为q值(Wang, Tuominen, and Tsai 2011)。如果有重复样本，函数将自动使用逻辑回归检验。

这个原理我也不是很懂。

可视化每条染色体的低甲基化/高甲基化碱基/区域的分布：

读取CpG岛（CpGi）和CpG海岸注释时报错了，不过好像没有什么关系文档里也这么用。

a GenomicRangesList contatining one GRanges object for flanks and one for GRanges object for the main feature.

NOTE: This can not return a GRangesList at the moment because flanking regions do not have to have the same column name as the feature. GRangesList elements should resemble each other in the column content. We can not satisfy that criteria for the flanks

读入甲基化调用文件之后会返回一个methylRawList对象。

但是当样本很多而且样本的规格很大的时候，可以返回不一样的对象

methylRawListDB, methylRawDB, methylBaseDB and methylDiffDB统称为methylDB objects。

Sometimes, when dealing with multiple samples and increased sample sizes coming from genome wide bisulfite sequencing experiments, the memory of your computer might not be sufficient enough.

Therefore methylKit offers a new group of classes, that are basically pendants to the original methylKit classes with one important difference: The methylation information, which normally is internally stored as data.frame, is stored in an external bgzipped file and is indexed by tabix (Li 2011), to enable fast retrieval of records or regions. This group contains methylRawListDB, methylRawDB, methylBaseDB and methylDiffDB, let us call them methylDB objects.

We can now create a methylRawListDB object, which stores the same content as myobj from above. But the single methylRaw objects retrieve their data from the tabix-file linked under dbpath.

今天的内容先分享到这里了，读完本文《西藏甲基化分析--山西甲基化分析》之后，是否是您想找的答案呢？想要了解更多，敬请关注www.zuqiumeng.cn,您的关注是给小编最大的鼓励。