treatment参数指定样本的分组，0代表control组，1代表treatment组bam⽂件
2. 合并所有样本的数据
将所有样本的甲基化情况合并，得到所有样本的甲基化表达谱，⽤法如下
meth=unite(myobj, destrand=FALSE)
meth 中的内容如下，其实就是之前的methylation calling⽂件的合并
在合并的过程中，默认情况下，只有所有的样本都包含该位点时，才会保留，本质就是取的所有样本的交集，如果你想要取并集，可以修改up参数的值，该参数的值代表每组中⾄少有多少个样本覆盖该位点时才保留，如果设置为1，就是取并集。
meth.min=unite(myobj,up=1L)
3. 执⾏差异分析
通过calculateDiffMeth函数来执⾏差异甲基化分析，⽤法如下
myDiff=calculateDiffMeth(meth)
根据甲基化C是变多了还是变少了，可以将差异甲基化的结果分为两⼤类：
1. hypermethylated
2. hypomethylated
hypermethylated表⽰相⽐control组，treatment组中的甲基化C更多；hypomethylated则相反，表⽰treatment组中的甲基化C⽐control组中少。
采⽤getMethylDiff函数提取差异分析的结果，⽤法如下
difference函数表明差异的阈值，只有差异⼤于该阈值时，才会保留，起始就是meth.diff的值，注意是绝对值⼤于difference的值。
除了difference阈值之外，还有qvalue阈值，⼩于该阈值的结果保留。在methylKit中，校正p值采⽤的是SILM算法，和我们常规的BH算法不同。type参数定义差异的类型，如果你只关注hypermethylated或者hypomethylated，可以设置type 参数的值，单独筛选。
diff函数
在methylKit中，它的差异分析总是针对合并后的甲基化表达谱，如果你的甲基化表达谱每⼀⾏是⼀个
甲基化位点，那么差异分析的结果就是差异甲基化位点；如果你的表达谱每⼀⾏是⼀个甲基化区域，那么差异分析的结果就是差异甲基化区域。上⾯的例⼦都是针对差异甲基化位点的，下⾯看下差异甲基化区域的分析。
⾸先遇到的问题就是甲基化区域如何界定，在methylKit中，按照滑动窗⼝的⽅式定义甲基化区域，默认窗⼝⼤⼩为10000 bp ，步长为10000bp，通过tileMethylCounts函数实现。
完整的差异甲基化区域分析的代码如下：

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。