CN108920904B

CN108920904B - 一种同源基因特异性甲基化时序数据的分析方法

Info

Publication number: CN108920904B
Application number: CN201810833596.1A
Authority: CN
Inventors: 吉冠玉; 吉红玉; 王君文; 胡琪; 高飞
Original assignee: Shenzhen E Gene Technology Co ltd
Current assignee: Shenzhen E Gene Technology Co ltd
Priority date: 2018-07-26
Filing date: 2018-07-26
Publication date: 2022-05-27
Anticipated expiration: 2038-07-26
Also published as: CN108920904A

Abstract

本发明公开了一种同源基因特异性甲基化时序数据的分析方法，属于基因组学和分子生物技术领域。该方法首先利用DNA甲基化测序数据序列特征和突变信息判定ASM区域，筛选出每个时期或者组别中特有的所述ASM区域，之后结合不同突变型对应的甲基化信息，利用某个时期的甲基化信息和其他不同时期样品进行比较，初步判定出原始时期特异性ASM区域，最后通过统计所有所述原始时期特异性ASM区域内的差异以及全部CpG位点数目，采用超几何检验，判定所述ASM区域是否为最终的时期特异性同源基因特异性甲基化区域。本发明进一步深化分析挖掘甲基化信息，可以检测多个不同阶段基于重亚硫酸盐处理的数据的ASM动态变化。

Description

一种同源基因特异性甲基化时序数据的分析方法

技术领域

本发明属于基因组学和分子生物技术领域，特别涉及一种同源基因特异性甲基化时序数据的分析方法。

背景技术

在DNA甲基化数据分析的过程中，特别是针对WGBS(全基因组甲基化)数据，很多人都觉得无从下手，常规的分析思路通常是找DMR(Differentially methylated regions，差异性甲基化区域)，并进行注释分析；但是此方案挖掘出来的信息非常有限，大量相关的DNA甲基化信息被遗漏。

ASM(Allele-Specific DNA Methylation，等位基因特异性甲基化)，除了我们通常关注的等位基因特异性甲基化与印记基因和女性X染色体的关系外，研究表明等位基因特异性甲基化影响等位基因间甲基化差异高达85％的完整CpG岛。等位基因间的甲基化差异与表型强烈相关。等位基因特异性甲基化能够导致甲基化的基因拷贝的等位基因抑制。同时，等位基因特异性甲基化可能影响了约10％的全部人类基因，并促进了等位基因特异性表达和单等位基因沉默。因此，等位基因特异性甲基化代表了遗传多态性如何能够导致表型变异的一条表观遗传途径。研究表明等位基因特异性甲基化是一种表观遗传漂变的结果，其方向由等位基因间的异常差异确定，同时获得一个等位基因高甲基化的倾向被遗传。

现有ASM确定的方法有两个类别，首先是通过针对BS数据进行SNP分型，选择突变频率接近50％的位置确定其不同突变型的DNA甲基化率差异，另外一种思路是通过DNA片段设计二项式分布模型，并结合贝叶斯方法判定ASM区间，要求区域内片段读长的甲基化分布表现为明显的两个类型的甲基化模式，并且两种模式的数量基本一致，如amrfinder软件，但是此种方法只针对单个样品；现阶段并没有适合判定多个阶段或者多个分组之间比较ASM特征以及ASM变化过程的方法。

发明内容

本发明的目的在于克服目前存在的技术缺失，提供一种同源基因特异性甲基化时序数据的分析方法。该方法在基于亚硫酸盐技术处理的DNA甲基化数据分析的基础上进一步深化分析挖掘甲基化信息，可以检测多个不同阶段基于重亚硫酸盐处理的数据的ASM动态变化，并可以用来检测对照组实验的ASM动态改变的数据。

本发明的目的通过下述技术方案实现：

一种同源基因特异性甲基化时序数据的分析方法，包括如下步骤：

首先利用DNA甲基化测序数据序列特征和突变信息判定ASM区域(同源基因特异性甲基化区域)，筛选出每个时期或者组别中特有的所述ASM区域，之后结合不同突变型对应的甲基化信息，利用某个时期的甲基化信息和其他不同时期样品进行比较，初步判定出原始时期特异性ASM区域(original-stage-specific-ASM)，最后通过统计所有所述原始时期特异性ASM区域内的差异以及全部CpG位点数目，采用超几何检验，判定所述ASM区域是否为最终的时期特异性同源基因特异性甲基化区域。

在其中的一个实施例中，所述同源基因特异性甲基化时序数据的分析方法，包括如下步骤：

(1)将重亚硫酸盐测序数据与参考基因组比对、排序，去除冗余测序数据；

(2)利用CpG位点的重亚硫酸盐转化C的个数信息得到CpG位点上覆盖到的C的个数R(m)以及经过重亚硫酸盐转化后的T的个数信息R(t)、以及甲基化率Methy_Rate＝R(m)/(R(m)+R(t))；

(3)利用甲基化数据得到每个样品的核苷酸突变信息，利用相应的突变信息将有关reads分为两类(A，B)，且每类的比例在0.3～0.7之间；分别统计区间内每条reads的对应位置的CpG位点的甲基化状态，定义读长中测序为C的位点为甲基化C(m)，测序为T的位点为非甲基化C(n)，并进一步统计测序读长的甲基化水平，methy_reads(i)＝C(m)/(C(m)+C(n))，i∈A U B；

同时判定所述两类reads中的甲基化率在CpG位点上的甲基化情况methy_rate(A_i)，methy_rate(B_i)；

利用t检验判定所述两类reads的甲基化率差异，若同时满足p(methy_reads(i))<0.05；p(methy_rate(i))<0.05，为保证统计有效性，该过程筛选CpG为点数>5；满足上述筛选条件则判定该区域为等位基因特异性甲基化区域(ASM)。

优选的，所述同源基因特异性甲基化时序数据的分析方法，还包括步骤(4)：利用步骤(3)中筛选出的所述方法在研究中多个时期的组别分别进行统计，筛选出多个时期中只有单独一个时期特有的ASM区域为待选时期特异性ASM区间U(ASM_i)，并筛选出所有位点的CpG位点的甲基化信息；进一步利用核苷酸突变信息分别统计各个时期的测序读长的甲基化信息reads_stage(I,A or B)，对应位点的甲基化信息Methy_rate_stage(I,A or B)，分别统计时期特异性的差异情况methy_reads(I,A or B)，对比reads_stage(I,A or B)；methy_rate(I,A or B)，对比Methy_rate_stage(I,A or B)；分别采用t检验p(reads_stage,A or B)<0.05；p(Methy_rate_stage,A or B)<0.05的结果确认疑似阶段特异性ASM区间U(ASM_i)；

优选的，所述同源基因特异性甲基化时序数据的分析方法，还包括步骤(5)：筛选步骤(4)中CpG甲基化信息，统计每个特异性ASM区间中的CpG位点的甲基化率在每个样品中的甲基化率，分别比较AB两种突变类型至少一种对应的reads的甲基化率有差异的位点个数k(ASM_i)，以及区间内的CpG个数m(CG_i)，总数为N，以及所有的疑似阶段特异性ASM区间U(ASM_i)中的总数n：

优选的，所述同源基因特异性甲基化时序数据的分析方法，还包括步骤(6)：采用超几何分布判别相应的区域内的富集情况，取p<0.05的ASM区域为显著时期特异性ASM区间：

优选的，步骤(3)中的CpG位点采用如下方法进行提取：首先合并有所样品的全部CpG位点，筛选出在每个样品中所述ASM区域内的CpG位点，对于同组中重复样品，筛选出重复样品中取交集，筛选出所有样品中都存在于ASM区域内的CpG位点。

本发明相对于现有技术具有如下的优点及效果：

本发明在基于亚硫酸盐技术处理的DNA甲基化数据分析的基础上进一步深化分析挖掘甲基化信息，可以检测多个不同阶段(>＝2)基于重亚硫酸盐处理的数据的ASM动态变化。

附图说明

图1是本申请的一个实施例的分析方法的流程图；

图2是实施例1的样品的前期的三个重复之间的交叠情况；

图3是实施例1的样品的中期的三个重复之间的交叠情况；

图4是实施例1的样品的后期的三个重复之间的交叠情况；

图5是实施例1的ASM对应基因区域在前中后三个时期的情况图；

图6是实施例1的合并后ASM区间长度与个数的分布情况图；

图7是实施例1中的ASM区域在基因上的位置分布情况图；

图8是实施例1中样品在三个时期特异性ASM的甲基化水平分布图。

具体实施方式

下面结合实施例对本发明作进一步详细的描述，但本发明的实施方式不限于此。

本申请的基于同源基因特异性甲基化分析时序甲基化数据的方法，包括如下步骤：

(1)将所检测的9例经过重亚硫酸盐测序数据与猪参考基因组比对、排序，去除冗余测序数据；此实验数据中9例猪全基因组甲基化数据共分为前中后三个时期，每个时期三个重复。

(2)利用CpG位点的重亚硫酸盐转化C的个数信息得到CpG位点上覆盖到的C的个数R(m)以及经过重亚硫酸盐转化后的T的个数信息R(t)；以及甲基化率Methy_Rate＝R(m)/(R(m)+R(t))；

利用t检验判定所述两类reads的甲基化率差异，若同时满足p(methy_reads(i))<0.05、p(methy_rate(i))<0.05，为保证统计有效性，该过程筛选CpG为点数>5；满足上述筛选条件则判定该区域为等位基因特异性甲基化区域(ASM)；

此处为保证实验数据的可靠性，结合本申请的研究目的为前中后三个时期特异性的ASM的变化，所以此处将三个时期的样品ASM结果取交集。对于样品间不同ASM边界取交集，最终得到每个时期的总的ASM的情况。

(4)利用步骤(3)中筛选出的所述方法在研究中多个时期的组别分别进行统计，筛选出多个时期中只有单独一个时期特有的ASM区域为待选时期特异性ASM区间U(ASM_i)，并筛选出所有位点的CpG位点的甲基化信息；进一步利用核苷酸突变信息分别统计各个时期的测序读长的甲基化信息reads_stage(I,A or B)，对应位点的甲基化信息Methy_rate_stage(I,A or B)，分别统计时期特异性的差异情况methy_reads(I,A or B)，对比reads_stage(I,A or B)；methy_rate(I,A or B)，对比Methy_rate_stage(I,A or B)；分别采用t检验p(reads_stage,A or B)<0.05；p(Methy_rate_stage,A or B)<0.05的结果确认疑似阶段特异性ASM区间U(ASM_i)；

(5)筛选步骤(4)中CpG位点的甲基化信息，统计每个特异性ASM区间中的CpG位点的甲基化率在每个样品中的甲基化率，分别比较AB两种突变类型至少一种对应的reads的甲基化率有差异的位点个数k(ASM_i)，以及区间内的CpG个数m(CG_i)，总数为N，以及所有的疑似阶段特异性ASM区间U(ASM_i)中的总数n：

(6)采用超几何分布判别相应的区域内的富集情况，取p<0.05的ASM区域为时期特异性ASM区间：

后续通过将上述ASM区间有重叠，或者间隔小于200bp区间合并，分别统计ASM的长度分布，并对区域注释到基因组，确定各基因元件分布的情况以及基因功能情况。其流程图如图1所示。

实施例1

(1)将所检测的9例经过重亚硫酸盐测序数据与猪参考基因组比对、排序，去除冗余测序数据；此实验数据中9例猪全基因组甲基化数据共分为前中后三个时期，每个时期三个重复；

此处为保证实验数据的可靠性，结合本申请的研究目的为前中后三个时期特异性的ASM的变化，所以此处将三个时期的样品ASM结果取交集；其交叠情况如图2-4所示；从图2-4可以看出，三个时期中的ASM区间在重复样品中大量重合，重合部分即为排除样品特异性时期相关ASM数据集。

(4)利用步骤(3)中筛选出的所述方法在研究中多个时期的组别分别进行统计，筛选出多个时期中只有单独一个时期特有的ASM区域为待选时期特异性ASM区间U(ASM_i)，并筛选出所有位点的CpG位点的甲基化信息；进一步利用核苷酸突变信息分别统计各个时期的测序读长的甲基化信息reads_stage(I,A or B)，对应位点的甲基化信息Methy_rate_stage(I,A or B)，分别统计时期特异性的差异情况methy_reads(I,A or B)，对比reads_stage(I,A or B)；methy_rate(I,A or B)，对比Methy_rate_stage(I,A or B)；分别采用t检验p(reads_stage,A or B)<0.05、p(Methy_rate_stage,A or B)<0.05的结果确认疑似阶段特异性ASM区间U(ASM_i)；

(5)筛选步骤(4)中CpG甲基化信息，统计每个特异性ASM区间中的CpG位点的甲基化率在每个样品中的甲基化率，分别比较A、B两种突变类型至少一种对应的reads的甲基化率有差异的位点个数k(ASM_i)，以及区间内的CpG个数m(CG_i)，总数为N，以及所有的疑似阶段特异性ASM区间U(ASM_i)中的总数n：

并对相应的区间做基因注释，最终取前中后三个时期独有的基因集用于后续研究，在此实施案例中，我们发现前期特异性ASM对应1539个特异性基因，中期对应1076个，后期对应基因3111个，结果见图5。

(7)对步骤(6)中的前期、中期、后期三个时期的相关基因进行GO和KEGG注释，确定所述相关基因中ASM对应的基因在整个发展过程中或者组间比较的过程中的功能层面的变化。

本实施例利用9例猪全基因组甲基化数据对方法进行测试，利用该方法判定的结果与实验预期的结果一致。

在整个研究过程中，ASM区域在基因组上大量的出现图5所示情形；结果也显示大量的基因区域有重合，如图5所示，但是在每个时期都有每个时期独自的ASM区域存在，说明在青春期发育过程中等位基因特异性甲基化的发生是一个动态的过程。

如图6所示，发生等位基因特异性甲基化的区域集中在3k以内。

针对这些各个时期特异性的ASM进行分析，我们从ASM区域在各个基因元件上的分布也可以看出，每个时期发生在各个基因元件上的区域分布情况也有明显的改变，结果见图7。

通过进一步考察这些区域的DNA甲基化，发现对应的甲基化率相较与其他时期都比较低，参见图8，这也从反面印证了我们检验区间的正确性。

针对印记相关的基因，我们发现印记相关的基因和对应的基因表达量并没有明显的规律，但是在前期和中期表现为表达下调的多，而后期趋势则相反，如图8所示。进一步，我们针对这些基因做功能注释，发现基因通路在前期和嗅觉传导(Olfactorytransduction)方面有密切的关系。在过程前期，生物过程层面和免疫应答等有密切的关系；在过程中期，细胞组分层面和细胞核有关系，而在后期，细胞组分层面和细胞外泌体、膜泡、细胞外细胞器以及细胞外区域有关系。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims

1.一种同源基因特异性甲基化时序数据的分析方法，其特征在于，包括如下步骤：

首先利用DNA甲基化测序数据序列特征和突变信息判定ASM区域，筛选出每个时期或者组别中特有的所述ASM区域，之后结合不同突变型对应的甲基化信息，利用某个时期的甲基化信息和其他不同时期样品进行比较，初步判定出原始时期特异性ASM区域，最后通过统计所有所述原始时期特异性ASM区域内的差异以及全部CpG位点数目，采用超几何检验，判定所述ASM区域是否为最终的时期特异性同源基因特异性甲基化区域。

2.根据权利要求1所述的同源基因特异性甲基化时序数据的分析方法，其特征在于，所述同源基因特异性甲基化时序数据的分析方法，包括如下步骤：

(3)利用甲基化数据得到每个样品的核苷酸突变信息，利用相应的突变信息将有关reads分为两类(A，B)，且每类的比例在0.3～0.7之间；分别统计区间内每条reads的对应位置的CpG位点的甲基化状态，定义读长中测序为C的位点为甲基化C(m)，测序为T的位点为非甲基化C(n)，并进一步统计测序读长的甲基化水平，methy_reads(i)＝C(m)/(C(m)+C(n))，i∈AUB；

利用t检验判定所述两类reads的甲基化率差异，若同时满足p(methy_reads(i))<0.05；p(methy_rate(i))<0.05，为保证统计有效性，该过程筛选CpG为点数>5；满足上述筛选条件则判定该区域为等位基因特异性甲基化区域。

3.根据权利要求2所述的同源基因特异性甲基化时序数据的分析方法，其特征在于，所述同源基因特异性甲基化时序数据的分析方法，还包括步骤(4)：利用步骤(3)中筛选出的所述方法在研究中多个时期的组别分别进行统计，筛选出多个时期中只有单独一个时期特有的ASM区域为待选时期特异性ASM区间U(ASM_i)，并筛选出所有位点的CpG位点的甲基化信息；进一步利用核苷酸突变信息分别统计各个时期的测序读长的甲基化信息reads_stage(I,A or B)，对应位点的甲基化信息Methy_rate_stage(I,A or B)，分别统计时期特异性的差异情况methy_reads(I,A or B)，对比reads_stage(I,A or B)；methy_rate(I,A orB)，对比Methy_rate_stage(I,A or B)；分别采用t检验p(reads_stage,A or B)<0.05；p(Methy_rate_stage,A or B)<0.05的结果确认疑似阶段特异性ASM区间U(ASM_i)。

4.根据权利要求3所述的同源基因特异性甲基化时序数据的分析方法，其特征在于，所述同源基因特异性甲基化时序数据的分析方法，还包括步骤(5)：筛选步骤(4)中CpG甲基化信息，统计每个特异性ASM区间中的CpG位点的甲基化率在每个样品中的甲基化率，分别比较AB两种突变类型至少一种对应的reads的甲基化率有差异的位点个数k(ASM_i)，以及区间内的CpG个数m(CG_i)，总数为N，以及所有的疑似阶段特异性ASM区间U(ASM_i)中的总数n：

5.根据权利要求4所述的同源基因特异性甲基化时序数据的分析方法，其特征在于，所述同源基因特异性甲基化时序数据的分析方法，还包括步骤(6)：采用超几何分布判别相应的区域内的富集情况，取p<0.05的ASM区域为显著时期特异性ASM区间：

6.根据权利要求2所述的同源基因特异性甲基化时序数据的分析方法，其特征在于，步骤(3)中的CpG位点采用如下方法进行提取：首先合并有所样品的全部CpG位点，筛选出在每个样品中所述ASM区域内的CpG位点，对于同组中重复样品，筛选出重复样品中取交集，筛选出所有样品中都存在于ASM区域内的CpG位点。