CN111724860B - 一种基于测序数据识别染色质开放区域的方法及装置 - Google Patents

Abstract

本申请公开了一种基于测序数据识别染色质开放区域的方法及装置。本申请方法包括，测序数据处理步骤，根据cfDNA全基因组测序数据在参考基因组的位置，计算窗口化保护评分和测序深度；初步筛选步骤，定位WPS的波峰和波谷，将相邻波峰的间距大于阈值的区域作为候选染色质开放区域；波峰参数筛选步骤，根据候选染色质开放区域两侧各3‑9个WPS的波峰面积、波峰夹角，波峰间距、波峰高度和波峰宽度进一步筛选；染色质开放区域识别步骤，通过cfDNA平均标准化测序深度验证，准确识别染色质开放区域。本申请方法，无需借助参考图谱，不依赖任何试剂或试验，仅通过数据分析即可简单、有效的获得染色质开放区域，提高了检测质量和效率。

Description

一种基于测序数据识别染色质开放区域的方法及装置

技术领域

本申请涉及核酸检测技术领域，特别是涉及一种基于测序数据识别染色质开放区域的方法及装置。

背景技术

核小体作为真核生物染色质的基本单位，是一种典型的DNA与组蛋白结合的生物大分子，是表观遗传学的重要内容。在不同细胞状态下，核小体定位是动态变化的，其位置的改变会影响转录因子与DNA的结合，从而调控基因特异性表达。

当DNA和蛋白结合形成核小体，DNA就会受到蛋白的保护而不被酶降解，在高通量测序中该区域测序深度高；相反，当DNA复制转录时，需要将转录调控区的紧密结构打开，从而允许一些调控因子，如转录因子、其他调控因子等与DNA相结合。这部分打开的染色质，被称为开放染色质，即染色质开放区域。染色质开放区域裸露的DNA会被酶降解，因此测序深度低。

开放染色质的研究方法主要包括：ATAC-Seq，以及传统的DNase-Seq和FAIRE-Seq等。其中，DNase-Seq主要是使用DNaseI内切酶识别开放染色质区域；FAIRE-Seq是先进行超声裂解，然后用酚-氯仿富集。ATAC-Seq是用Tn5转座酶，随后对酶捕获到的DNA序列进行富集和扩增；该方法所需细胞量少，实验简单，可以在全基因组范围内检测染色质的开放状态，目前已经成为研究染色质开放性的首选技术方法。

以上几种方法都是从实验手段获得开放染色质区域，对实验操作人员的技术要求高，建库的准备工作较为繁琐；并且，试剂比较昂贵、对酶的依赖性强，存在酶切割偏好性等问题。

因此，亟需研发一种更简单有效的染色质开放区域检测方法，以减小对试剂和技术的依赖和要求。

发明内容

本申请的目的是提供一种新的基于测序数据识别染色质开放区域的方法、装置和存储介质。

为了实现上述目的，本申请采用了以下技术方案：

本申请的第一方面公开了一种基于测序数据识别染色质开放区域的方法，其中，测序数据为cfDNA全基因组测序数据，该方法包括以下步骤：

测序数据处理步骤，包括根据cfDNA全基因组测序数据在参考基因组上的位置，结合区间划分和滑动窗口，计算窗口化保护评分(windowed protection score，缩写WPS)和测序深度；

初步筛选步骤，包括定位窗口化保护评分的波峰和波谷，将相邻波峰的间距大于阈值的区域作为候选染色质开放区域，其中，阈值大于或等于核小体大小；

波峰参数筛选步骤，包括对每个候选染色质开放区域两侧各3-9个窗口化保护评分的波峰面积、波峰夹角、波峰间距、波峰高度和波峰宽度，总计五个波峰参数进行计算；同时，对持家基因转录起始区域的窗口化保护评分的相同的五个波峰参数进行计算，并将其作为波峰参数阈值；筛选候选染色质开放区域的波峰参数在波峰参数阈值范围内的区域，用于后续判断；

染色质开放区域识别步骤，包括对波峰参数筛选步骤获得的候选染色质开放区域，使用该区域的cfDNA平均标准化测序深度进行验证，判断各候选染色质开放区域是否出现cfDNA覆盖度降低的现象，如果出现cfDNA覆盖度降低的现象，则将其识别为染色质开放区域。

其中，波峰参数筛选步骤，选取候选染色质开放区域两侧各3-9个窗口化保护评分，主要原因是3-9个WPS波峰占据的测序范围能够获得可靠的结果。进一步优选的，对每个候选染色质开放区域两侧各5-7个WPS的波峰面积、波峰夹角、波峰间距、波峰高度和波峰宽度进行计算；因为，选取3或4个WPS波峰的参考范围略小，在一些试验中存在一定的误差，而选取8个或以上的WPS波峰，区域范围过大，核小体排列出现部分波峰不整齐的情况，会影响识别效果，并且计算复杂度也有所提升。因此，综合考虑可靠性和计算效率优选5-7个WPS波峰。本申请的一种实现方式中，更优选的，取5个WPS，其波峰占据约2000bp左右的测序范围，已经能够较好地获得可靠结果。

需要说明的是，本申请的方法，利用核小体缺失时，WPS波峰间距会呈现一定的特征，基于此对可能的染色质开放区域进行筛选，即初步筛选步骤；然后，采用持家基因转录起始区域WPS的波峰面积、波峰夹角、波峰间距、波峰高度和波峰宽度，作为波峰参数阈值，进一步根据染色质开放区域的波峰参数特征筛选符合要求的区域，即波峰参数筛选步骤；最后，利用cfDNA平均标准化测序深度进行验证，将符合条件的区域识别为染色质开放区域，即染色质开放区域识别步骤。本申请的方法，只需要对cfDNA全基因组测序数据进行分析，通过染色质开放区域cfDNA覆盖深度低和WPS波峰消失两个特性查找染色质开放区域，无需借助任何参考图谱，也无需依赖任何试验技术或酶；并且，结合波形的参数特征，包括波峰的面积、夹角、波峰间距、波峰高度、波峰宽度多个指标准确识别WPS波峰消失的区间，降低假阳性和假阴性染色质开放区域。

还需要说明的是，现有利用WPS波形进行核小体峰识别的方法，仅仅适用于查找两组数据间的核小体峰差异，无法实现单一样本的染色质开放区域的识别；本申请创造性的利用相邻波峰间距、波峰参数和cfDNA平均标准化测序深度，三个参数特征，逐步对单一样本的染色质开放区域进行分析，提高了染色质开放区域的检测效率和质量。

优选的，染色质开放区域识别步骤中，判断各候选染色质开放区域是否出现cfDNA覆盖度降低的现象，具体包括，计算全基因组随机等长区域内cfDNA的测序深度，当cfDNA平均标准化测序深度小于全基因组随机等长区域内cfDNA的测序深度的n倍时，则将判断该区域出现了cfDNA覆盖度降低的现象，其中n＜0.6。

其中，n＜0.6是本申请通过对数据采样、实验分析之后得到的值。通过对数据进行采样分析本申请得到的一个比较合适的范围为0.15<n<0.6，一般染色质开放区域的标准化测序深度均在这个范围之内。n<0.6的范围内都能较好地识别出可能的染色质开放区域，特别是，当n的取值贴近0.35左右的区间，查准率最高。当n的取值向上增加，查全率增加，即更多的染色质开放区域被检测出来，但是差准率下降，即假阳性的占比增高。因此，为了降低识别错误的概率，可以调低n的取值，通过调整n的取值在查全率和查准率之间进行平衡。本申请的一种实现方式中，优选的n＜0.5，可以在查全率和查准率两者之间有较好的平衡效果。

需要说明的是，如本申请背景技术中提到的，染色质开放区域裸露的DNA会被酶降解，因此导致测序深度低；本申请正是利用这个特征，最终通过判断各候选染色质开放区域是否出现cfDNA覆盖度降低的现象，来确认该区域是否为染色质开放区域。可以理解，将cfDNA平均标准化测序深度与全基因组随机等长区域内cfDNA的测序深度进行比较，只是本申请的一种实现方式采用的判断是否出现cfDNA覆盖度降低的方法，不排除还可以采用其它方式进行覆盖度判断，在此不作具体限定。

优选的，测序深度标准化的方法，具体包括，将参考基因组按m kb大小划分为不同的区间，计算m kb区间cfDNA的覆盖度，并进行标准化，具体公式如式一所示，

式一

其中，x^*为标准化后的值，x为初始值，x_min为区间最小值，x_max为区间最大值。

优选的，10kb≤m≤50kb。

优选的，测序数据处理步骤中，计算窗口化保护评分和测序深度具体包括，将参考基因组按m kb大小划分为不同的区间，每个区间再以nbp大小为滑动窗口，计算每个窗口的cfDNA读长覆盖度和窗口化保护评分。

优选的，10kb≤m≤50kb。

优选的，120bp≤n≤150bp。

需要说明的是，本申请通过划分区间，可以使用多线程技术加快计算速度，还可以使用分布式计算技术更快加速WPS的计算；至于具体的每个区间的大小，可以根据基因组的大小和实际的处理能力或速度要求而定，一般10kb≤m≤50kb即可。

优选的，测序数据处理步骤还包括使用Savitzky-Golay滤波器(缩写SG滤波器)消除窗口化保护评分的波形噪声，Savitzky-Golay滤波器的输入数据为测序数据处理步骤获得的一维的窗口化保护评分波形数组，数组的索引为参考基因组位置，数组存储的内容为当前基因组位置计算得到的窗口化保护评分数值；Savitzky-Golay滤波器通过kbp大小的滑动窗口，通过滑动窗口内的所有点得到局部多项式回归方程，最后计算窗口中点的预测值，作为每个碱基位点的窗口化保护评分的校正值。

优选的，25bp≤k≤50bp。

其中，局部多项式回归方程一般是通过这k个点以最小二乘拟合为优化的目标函数，拟合出一个多项式函数，最后求出窗口中心位置的函数值作为该中心位置的预测值。例如，对应k＝31，中心位置的坐标为60，需要选取31个点，这31个点的坐标区间为[45，75]，通过预设的多项式次数(选取为2)以及这31个点的值，以最小二乘拟合为优化的目标函数，计算得出一个二阶多项式函数，通过这个获得的函数将窗口中心坐标60代入进去，就可以计算出该位置的预测值，或称预估值、拟合值。具体到本申请，SG滤波器为基于最小二乘的局部多项式回归算法，局部指的是该算法采用的是一个个小窗口，即窗口大小为k，范围为25bp≤k≤50bp，多项式指的是本申请选择的拟合该区域的函数为多项式函数，本申请预先给定该多项式的次数，例如取值为2，基于最小二乘指的是我们使用最小二乘法为优化目标来获得拟合效果比较好的函数；回归指的是本申请通过拟合得到的函数获得所需要预测的中心点的值；对于每一个中心点，SG滤波器算法都会取一个大小为k的窗口，计算得到多项式函数，得到预测值，因此需要滑动中心点的位置，即滑动窗口。总的来说，k是本申请选取的SG滤波器的窗口大小，取值范围为25bp≤k≤50bp；而SG滤波器的输入数据的大小即本申请的参考基因组m，m取值范围为10kb≤m≤50kb，需要的滤波的范围也就是这个范围。通过滑动窗口对改该区间的每一个点进行滤波操作。另外，SG滤波器还有一个参数为多项式的次数，本申请选择为2，用二次多项式来拟合函数。

需要说明的是，由于WPS波形存在一定的由测序带来的噪声干扰，并且其波形数据的局部相似性较高，因此需要使用能很好地保留WPS波峰形状和宽度信息的滤波方法消除WPS波形噪音；本申请的一种实现方式中优选采用Savitzky-Golay滤波器，不排除还可以采用其它的类似功能的滤波器。至于滑动窗口的大小和每次滑动窗口的取点数可以根据实际情况而定。

优选的，初步筛选步骤中，定位窗口化保护评分的波峰和波谷，具体包括，使用三个大小一致的小窗口，分别记为窗口1、窗口2、窗口3，三个窗口依次间隔5bp，滑动三个窗口，并计算三个窗口的窗口化保护评分数值之和，分别记为sum_win1、sum_win2、sum_win3；若sum_win1<sum_win2，且sum_win3<sum_win2，则为波峰；若sum_win1>sum_win2，且sum_win3>sum_win2，则为波谷。

其中，三个窗口间隔5bp，是因为波峰位置的跨度一般在15bp左右，经过试验分析认为，选择在5bp左右比较合适。

优选的，小窗口的大小为4-7bp。

其中，小窗口的大小优选为4-7bp，当窗口小于4bp(不包含4bp)时，因为窗口较小，容易受到测序误差的影响，波峰的定位结果产生抖动，将一些局部最大值或最小值识别为波峰或者波谷；而当窗口大于7bp(不包含7bp)时，因为窗口较大，可能直接将宽度较小的波峰或者波谷跨越过去，造成波峰定位的不准确。

需要说明的是，以上方法只是本申请的一种实现方式中具体判断波峰波谷的一种方法，不排除还可以采用其它判断方法，在此不作具体限定。其中，小窗口的大小可以根据具体项目分析需求而定。

优选的，本申请的方法中，在初步定位了波峰和波谷后，还包括对波峰或波谷的可靠性进行判断，具体包括计算波峰波谷的波峰高度、半波形宽度、相邻波峰宽度、相对波谷宽度、相邻波峰相对高度，总计五个波峰数值；同时，对持家基因转录起始区域的窗口化保护评分的相同的五个波峰数值进行计算，并将其作为波峰数值阈值；波峰波谷的五个波峰数值在所述波峰数值阈值范围内，则认为该位置为可靠的波峰或者波谷。最终去掉可靠性差的波峰波谷，获得全基因组范围的波峰波谷集合，用于后续分析。

优选的，初步筛选步骤中，阈值为250bp，即将相邻波峰的间距大于250bp的区域作为候选染色质开放区域。

需要说明的是，初步筛选步骤主要是根据核小体缺失导致的相邻波峰间距特征进行筛选；一般来说，核小体的大小为167bp左右，因此，WPS波峰间距应该大于200bp；为了是筛选更加准确，本申请将相邻波峰间距大于250bp的区域判定为可能的染色质开放区域。

优选的，波峰参数筛选步骤中，计算波峰面积的公式如式二所示，计算波峰夹角的公式如式三所示；

式二

式三

式二中，x为波谷在基因组上的坐标；i为对应的第几个波谷的编号，范围为1-9；y为该碱基位点对应的窗口化保护评分数值，d_x为x的微分；式三中，A为波峰位点，B和C分别为波峰位点A两侧相应的波谷位点。

本申请的第二方面公开了一种基于测序数据识别染色质开放区域的装置，该装置用于基于cfDNA全基因组测序数据进行染色质开放区域识别，该装置包括测序数据处理模块、初步筛选模块、波峰参数筛选模块和染色质开放区域识别模块；

测序数据处理模块，包括用于根据cfDNA全基因组测序数据在参考基因组上的位置，结合区间划分和滑动窗口，计算窗口化保护评分和测序深度；

初步筛选模块，包括用于定位窗口化保护评分的波峰和波谷，将相邻波峰的间距大于阈值的区域作为候选染色质开放区域，阈值大于或等于核小体大小；

波峰参数筛选模块，包括用于对每个候选染色质开放区域两侧各3-9个窗口化保护评分的波峰面积、波峰夹角、波峰间距、波峰高度和波峰宽度，总计五个波峰参数进行计算；同时，对持家基因转录起始区域的窗口化保护评分的相同的五个波峰参数进行计算，并将其作为波峰参数阈值；筛选候选染色质开放区域的波峰参数在波峰参数阈值范围内的区域，用于后续判断；

染色质开放区域识别模块，包括用于对波峰参数筛选模块获得的候选染色质开放区域，使用该区域的cfDNA平均标准化测序深度进行验证，判断各候选染色质开放区域是否出现cfDNA覆盖度降低的现象，如果出现cfDNA覆盖度降低的现象，则将其识别为染色质开放区域。

需要说明的是，本申请的装置，实际上就是通过各模块实现本申请的方法的各个步骤，因此，本申请装置中各模块的具体实现方式或参数条件可以参考本申请的方法。例如，染色质开放区域识别模块中，具体如何判断各候选染色质开放区域是否出现cfDNA覆盖度降低的现象，如何计算全基因组随机等长区域内cfDNA的测序深度；测序数据处理模块中，具体如何计算窗口化保护评分和测序深度，使用Savitzky-Golay滤波器消除WPS的波形噪声；初步筛选模块中，具体如何定位窗口化保护评分的波峰和波谷，如何对波峰或波谷的可靠性进行判断，初步筛选模块的阈值；以及波峰参数筛选模块中的波峰面积和波峰夹角的计算公式等；都可以参考本申请的基于测序数据识别染色质开放区域的方法。

本申请的第三方面公开了一种基于测序数据识别染色质开放区域的装置，该装置包括存储器和处理器；存储器，用于存储程序；处理器，用于通过执行存储器存储的程序以实现本申请的基于测序数据识别染色质开放区域的方法。

本申请的第四方面公开了一种计算机可读存储介质，其包括程序，该程序能够被处理器执行以实现本申请的基于测序数据识别染色质开放区域的方法。

由于采用以上技术方案，本申请的有益效果在于：

本申请基于测序数据识别染色质开放区域的方法，仅仅通过对cfDNA全基因组测序数据进行分析，通过染色质开放区域cfDNA覆盖深度低和WPS波峰消失两个特性查找染色质开放区域，无需借助任何参考图谱，也不依赖于任何试剂、酶或试验，能够简单、有效的从单一样本的cfDNA全基因组测序数据中检测获得染色质开放区域，且降低了假阳性和假阴性，提高了全基因组范围搜索染色质开放区域的质量和效率。

附图说明

图1是本申请实施例中基于cfDNA测序数据识别染色质开放区域的方法的流程框图；

图2是本申请实施例中基于cfDNA测序数据识别染色质开放区域的装置的结构框图；

图3是本申请实施例中Savitzky-Golay滤波器效果图和滤波器算法原理分析图，其中，(a)图为滤波之前的WPS波形，(b)图为滤波后的WPS波形，(c)图为Savitzky-Golay滤波器原理分析图；

图4是本申请实施例中滑动窗口方法寻找WPS波峰波谷的示意图；

图5是本申请实施例中WPS波形示意图，包括波峰、波谷、波峰高度、半波形宽度、相邻波峰宽度、相对波谷宽度、相邻波峰相对高度等七个参数的示意图；

图6是本申请实施例中候选染色质开放区域两侧各5个WPS的波峰面积、波峰夹角、波峰间距、波峰高度、波峰宽度的示意图；

图7是本申请实施例中检测出的染色质开放区间中心和持家基因的转录起始点的距离的分布服从均值为0的类正态分布；

图8是本申请实施例中检测出的染色质开放区域与不同组织特异性基因TSS区域的重叠百分比。

具体实施方式

下面通过具体实施方式结合附图对本申请作进一步详细说明。在以下的实施方式中，很多细节描述是为了使得本申请能被更好的理解。然而，本领域技术人员可以毫不费力的认识到，其中部分特征在不同情况下是可以省略的，或者可以由其他元件、材料、方法所替代。在某些情况下，本申请相关的一些操作并没有在说明书中显示或者描述，是为了避免本申请的核心部分被过多的描述所淹没，而对于本领域技术人员而言，详细描述这些相关操作并不是必要的，他们根据说明书中的描述以及本领域的一般技术知识即可完整了解相关操作。

现有的染色质开放区域研究都需要采用特殊的酶或试剂，通过试验技术，进行检测和判断。利用WPS波形进行核小体峰识别，主要是用于对比分析两组数据，用于查找两组数据间的核小体峰差异。因此，目前尚未有直接根据测序数据特征分析识别染色质开放区域的方案。

本申请研究分析认为，cfDNA主要是人体各组织细胞凋亡后产生的，在凋亡过程中，不受核小体保护的DNA会被酶分解；因此，原理上可以通过cfDNA测序数据的测序深度和WPS波峰进行染色质开放区域识别。具体的，在染色质开放区域的DNA无核小体保护，片段末端定位随机分布，WPS波峰不明显，在染色质开放区域两侧有良好定位的核小体；因此，染色质开放区域两侧片段末端定位具有规律性，WPS波峰形态良好，排列整齐有规律。基于以上研究和认识，本申请创造性的提出可以通过染色质开放区域cfDNA覆盖深度低和WPS波峰消失两个特性，直接从单一样本的测序数据中查找染色质开放区域；使得染色质开放区域的识别无需依赖任何试剂或试验技术，只需要对cfDNA全基因组测序数据进行相应的分析即可，更加简单有效。具体的，本申请基于测序数据识别染色质开放区域的方法，如图1所示，包括测序数据处理步骤11、初步筛选步骤12、波峰参数筛选步骤13和染色质开放区域识别步骤14。

其中，测序数据处理步骤11，包括根据cfDNA全基因组测序数据在参考基因组上的位置，结合区间划分和滑动窗口，计算WPS和测序深度。

本申请的一种实现方式中，对cfDNA全基因组测序数据利用比对软件比对到参考基因组，将比对到参考基因组上唯一位置的序列，进行后续生物信息分析。WPS(windowedprotection score)是一个通过滑动窗口计算而来的窗口保护分数，即通过cfDNA片段末端的分布来推导人基因组中核小体的分布位置。具体而言，预期cfDNA片段末端应聚集在核小体边界附近，基本不分布在核小体所在位置。基因组上每个窗口内每个碱基位点的WPS指标的计算方法是通过计算跨越nbp大小窗口的长度在120-180bp范围内的片段数量减去有一端落入窗口的长度在120-180bp范围内片段数量获得WPS的值。WPS值与强烈定位的阵列内核小体的位置关联。高WPS值代表该区域DNA被核小体保护，低WPS值代表DNA未被保护。

初步筛选步骤12，包括定位窗口化保护评分的波峰和波谷，将相邻波峰的间距大于阈值的区域作为候选染色质开放区域，该阈值大于或等于核小体大小。

本申请的一种实现方式中，定位WPS波峰和波谷的算法分为两个主要步骤，第一步是通过滑动窗口方法获取可能的波峰波谷位置，第二步是计算波峰高度、半波形宽度、相邻波峰宽度、相对波谷宽度、相邻波峰相对高度这五个参数，参数数值在阈值范围，即根据持家基因转录起始区域的WPS波形计算五个参数的阈值范围，则认为该位置为可靠的波峰或者波谷。其中，通过滑动窗口方法获取可能的波峰波谷位置，是使用三个大小一致的小窗口来分析实现的。初步筛选步骤主要是根据当缺失一个核小体时，WPS波峰间距的阈值进行筛选。核小体大小为167bp，当缺失一个核小体时，WPS波峰间距应该大于200bp；因此，本申请的一种实现方式中，将相邻波峰间距大于250bp时，将该区域判断为可能是一个染色质开放区域。

波峰参数筛选步骤13，包括对每个候选染色质开放区域两侧各3-9个窗口化保护评分的波峰面积、波峰夹角、波峰间距、波峰高度和波峰宽度，总计五个波峰参数进行计算；同时，对持家基因转录起始区域的窗口化保护评分的相同的五个波峰参数进行计算，并将其作为波峰参数阈值；筛选候选染色质开放区域的波峰参数在波峰参数阈值范围内的区域，用于后续判断。

本申请的一种实现方式中，不仅计算了每个候选染色质开放区域两侧各3-9个WPS的波峰面积、波峰夹角、波峰间距、波峰高度和波峰宽度；而且，计算了该区域两侧各5个WPS的波峰面积、波峰夹角、波峰间距、波峰高度、波峰宽度等数值的均值和方差，其中，平均值的计算公式如式四所示，方差的计算公式如式五所示：

式四

式五

式四和式五中，

为平均值，n＝10，S_d ²为方差。

根据持家基因转录起始区域的WPS波形计算波峰的面积、夹角、波峰间距、波峰高度、波峰宽度等的数值作为参考阈值，计算指标数值满足参考阈值范围时，则认为该区域周边WPS波峰规则，中间缺失至少一个核小体，为候选的开放染色质区域。

染色质开放区域识别步骤14，包括对波峰参数筛选步骤获得的候选染色质开放区域，使用该区域的cfDNA平均标准化测序深度进行验证，判断各候选染色质开放区域是否出现cfDNA覆盖度降低的现象，如果出现cfDNA覆盖度降低的现象，则将其识别为染色质开放区域。

本申请的一种实现方式中，通过将参考基因组按m kb大小(10kb≤m≤50kb)划分不同的区间，计算m kb区间cfDNA的覆盖度，并进行标准化。当标准化后的cfDNA覆盖度小于全基因组随机等长区域内cfDNA的测序深度的n(n<0.5)倍时，则认为该区域确实出现了cfDNA覆盖度降低的现象，符合开放染色质区域的筛选条件。

基于测序数据识别染色质开放区域的方法：

(1)传统地从实验手段获得开放染色质区域，对实验操作人员的技术要求高，建库的准备工作较为繁琐，试剂比较昂贵、对酶的依赖性强，存在酶切割偏好性等问题，本申请方法无需依赖特殊实验操作；

(2)本申请仅依靠cfDNA全基因组测序数据，通过染色质开放区域cfDNA覆盖深度低和WPS波峰消失两个特性查找染色质开放区域，无需借助任何参考图谱；而现有的利用WPS波形进行核小体峰识别的方法仅适用于查找两组数据间的核小体峰差异，无法实现单一样本的染色质开放区域的识别；

(3)本申请同时结合波形的参数特征，包括波峰的面积、夹角、波峰间距、波峰高度、波峰宽度多个指标准确识别WPS波峰消失的区间，降低假阳性和假阴性染色质开放区域；

(4)本申请方法通过划分区间，可以使用多线程技术加快计算速度，还可以使用分布式计算技术更快加速WPS的计算，提高全基因组范围搜索染色质开放区域的效率。

本领域技术人员可以理解，上述方法的全部或部分功能可以通过硬件的方式实现，也可以通过计算机程序的方式实现。当上述方法中全部或部分功能通过计算机程序的方式实现时，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘、光盘、硬盘等，通过计算机执行该程序以实现上述功能。例如，将程序存储在设备的存储器中，当通过处理器执行存储器中程序，即可实现上述全部或部分功能。另外，当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序也可以存储在服务器、另一计算机、磁盘、光盘、闪存盘或移动硬盘等存储介质中，通过下载或复制保存到本地设备的存储器中，或对本地设备的系统进行版本更新，当通过处理器执行存储器中的程序时，即可实现上述实施方式中全部或部分功能。

因此，基于本申请的方法，本申请提出了一种基于测序数据识别染色质开放区域的的装置，如图2所示，包括测序数据处理模块21、初步筛选模块22、波峰参数筛选模块23和染色质开放区域识别模块24。

其中，测序数据处理模块21，包括用于根据cfDNA全基因组测序数据在参考基因组上的位置，结合区间划分和滑动窗口，计算窗口化保护评分和测序深度；初步筛选模块22，包括用于定位窗口化保护评分的波峰和波谷，将相邻波峰的间距大于阈值的区域作为候选染色质开放区域，阈值大于或等于核小体大小；波峰参数筛选模块23，包括用于对每个候选染色质开放区域两侧各3-9个窗口化保护评分的波峰面积、波峰夹角、波峰间距、波峰高度和波峰宽度，总计五个波峰参数进行计算；同时，对持家基因转录起始区域的窗口化保护评分的相同的五个波峰参数进行计算，并将其作为波峰参数阈值；筛选候选染色质开放区域的波峰参数在波峰参数阈值范围内的区域，用于后续判断；染色质开放区域识别模块24，包括用于对波峰参数筛选模块23获得的候选染色质开放区域，使用该区域的cfDNA平均标准化测序深度进行验证，判断各候选染色质开放区域是否出现cfDNA覆盖度降低的现象，如果出现cfDNA覆盖度降低的现象，则将其识别为染色质开放区域。

本申请的装置，利用各模块相互协调作用，能够实现基于测序数据的染色质开放区域识别，特别是通过本申请装置的各模块能够实现本申请基于测序数据识别染色质开放区域的方法中的相应的各个步骤，使染色质开放区域识别能够实现自动化处理。

本申请的另一实现方式中还提供了一种基于测序数据识别染色质开放区域的装置，该装置包括存储器和处理器；存储器，包括用于存储程序；处理器，包括用于通过执行存储器存储的程序以实现以下方法：测序数据处理步骤，包括根据cfDNA全基因组测序数据在参考基因组上的位置，结合区间划分和滑动窗口，计算WPS和测序深度；初步筛选步骤，包括定位窗口化保护评分的波峰和波谷，将相邻波峰的间距大于阈值的区域作为候选染色质开放区域，其中，阈值大于或等于核小体大小；波峰参数筛选步骤，包括对每个候选染色质开放区域两侧各3-9个窗口化保护评分的波峰面积、波峰夹角、波峰间距、波峰高度和波峰宽度，总计五个波峰参数进行计算；同时，对持家基因转录起始区域的窗口化保护评分的相同的五个波峰参数进行计算，并将其作为波峰参数阈值；筛选候选染色质开放区域的波峰参数在波峰参数阈值范围内的区域，用于后续判断；染色质开放区域识别步骤，包括对波峰参数筛选步骤获得的候选染色质开放区域，使用该区域的cfDNA平均标准化测序深度进行验证，判断各候选染色质开放区域是否出现cfDNA覆盖度降低的现象，如果出现cfDNA覆盖度降低的现象，则将其识别为染色质开放区域。其中，各步骤的具体实现方式可以参考本申请的基于测序数据识别染色质开放区域的方法。

本申请另一种实现方式中还提供一种计算机可读存储介质，该存储介质中包括程序，该程序能够被处理器执行以实现如下方法：测序数据处理步骤，包括根据cfDNA全基因组测序数据在参考基因组上的位置，结合区间划分和滑动窗口，计算WPS和测序深度；初步筛选步骤，包括定位窗口化保护评分的波峰和波谷，将相邻波峰的间距大于阈值的区域作为候选染色质开放区域，其中，阈值大于或等于核小体大小；波峰参数筛选步骤，包括对每个候选染色质开放区域两侧各3-9个窗口化保护评分的波峰面积、波峰夹角、波峰间距、波峰高度和波峰宽度，总计五个波峰参数进行计算；同时，对持家基因转录起始区域的窗口化保护评分的相同的五个波峰参数进行计算，并将其作为波峰参数阈值；筛选候选染色质开放区域的波峰参数在波峰参数阈值范围内的区域，用于后续判断；染色质开放区域识别步骤，包括对波峰参数筛选步骤获得的候选染色质开放区域，使用该区域的cfDNA平均标准化测序深度进行验证，判断各候选染色质开放区域是否出现cfDNA覆盖度降低的现象，如果出现cfDNA覆盖度降低的现象，则将其识别为染色质开放区域。同样的，各步骤的具体实现方式可以参考本申请的基于测序数据识别染色质开放区域的方法。

下面通过具体实施例和附图对本申请作进一步详细说明。以下实施例仅对本申请进行进一步说明，不应理解为对本申请的限制。

实施例1

本例对100个健康人cfDNA全基因组测序数据进行合并，共180G数据，利用BWA比对软件比对到人参考基因组GRCH37版本，将比对到参考基因组上唯一位置的序列，进行后续生物信息分析，用于染色质开放区域识别，具体如下：

测序数据处理步骤，包括根据cfDNA全基因组测序数据在参考基因组上的位置，结合区间划分和滑动窗口，计算WPS和测序深度。

cfDNA在转录起始位点、转录因子结合位点、其它调控因子结合部位等染色质开放区域的分布与核小体分布区存在明显的差异性片段特征信号。染色质开放区域核小体消失，其周边区域有排列整齐的核小体，其中开放区域cfDNA的测序深度降低，深度波峰消失或波峰相对高度降低；有核小体保护的区域，可以被测序的cfDNA片段较多，cfDNA测序深度较高，波峰信号强烈，深度波峰排列规则。

本例计算WPS、测序深度指标步骤为：将参考基因组按20kb划分不同的区间，计算每20kb区间cfDNA的覆盖度，并进行标准化，具体公式如式一所示，

式一

其中，x^*为标准化后的值，x为初始值，x_min为区间最小值，x_max为区间最大值。

每个20kb区间再以120bp为滑动窗口，计算每个窗口的WPS分值。

WPS(windowedprotection score)是一个通过滑动窗口计算而来的窗口保护分数。基因组上每个窗口内每个碱基位点的WPS指标的计算方法是通过计算跨越120bp窗口的长度在120-180bp范围内的片段数量减去有一端落入窗口的长度在120-180bp范围内片段数量获得WPS的值。WPS值与强烈定位的阵列内核小体的位置关联。高WPS值代表该区域DNA被核小体保护，低WPS值代表DNA未被保护。

进一步地，由于WPS波形存在一定的由测序带来的噪声干扰，其波形数据的局部相似性较高，因此需要使用能很好地保留WPS波峰形状和宽度信息的滤波方法消除WPS波形噪音。本例使用Savitzky-Golay滤波器(缩写SG滤波器)消除WPS波形噪声。SG滤波器原理具体如图3的(c)图所示，SG滤波器的输入数据为一维的WPS波形数组，数组的索引为参考基因组位置，数组存储的内容为当前基因组位置计算得到的WPS值。SG滤波器通过35bp滑动窗口截取局部数据，以最小二乘拟合的残差为优化目标计算35个点的局部多项式回归方程，最后计算窗口中点的预测值，作为每个碱基位点的WPS校正值。图3中的(a)图为滤波之前的WPS波形，(b)图为SG滤波之后的结果，黑色散点为原始WPS数据。图3的结果显示，SG滤波器成功的消除了波形噪声。

初步筛选步骤，包括定位窗口化保护评分的波峰和波谷，将相邻波峰的间距大于阈值的区域作为候选染色质开放区域，阈值大于或等于核小体大小。

为了找到染色质开放区域，首先需要定位WPS波形的波峰和波谷。定位波峰波谷的算法分为两个主要步骤，第一步是通过滑动窗口方法获取可能的波峰波谷位置，第二步是计算波峰高度(height)、半波形宽度(width)、相邻波峰宽度(Peak Distance)、相对波谷宽度(Troughwidth)、相邻波峰相对高度(Relative peak height)这五个参数，参数数值在阈值范围，根据持家基因转录起始区域的WPS波形计算五个参数的阈值范围，则认为该位置为可靠的波峰或者波谷。算法首先采用滑动窗口方法寻找波峰波谷，使用三个5bp的小窗口，如图4所示，三个窗口依次间隔5bp，滑动这三个窗口，并计算着三个窗口的WPS值之和，分别记为sum_win1、sum_win2、sum_win3。若sum_win1<sum_win2，且sum_win3<sum_win2，则为波峰；若sum_win1>sum_win2，且sum_win3>sum_win2，则为波谷。如图5所示，计算波峰高度、半波形宽度、相邻波峰宽度、相对波谷宽度、相邻波峰相对高度五个参数。根据1000个持家基因转录起始区域的WPS波形计算五个参数的阈值范围如下。窗口w大小为5bp，波峰高度height需大于0.58，半波形宽度width范围为25-100bp，相邻波峰宽度Peakwidth和相对波谷宽度Troughwidth的范围均为50-220bp，相邻波峰相对高度Relative peak height需小于0.2。如果波形参数满足阈值条件，则认为该位置为可靠的波峰或波谷。继续滑动窗口计算下一个可能的波峰或者波谷。最后获得全基因组范围的波峰波谷集合。

本例根据核小体大小为167bp，当缺失一个核小体时，WPS波峰间距应该大于200bp的特征，基于计算所得的波峰位置，当相邻波峰间距大于250bp左右时，认为该区域可能是一个染色质开放区域。

波峰参数筛选步骤，包括对每个候选染色质开放区域两侧各5个窗口化保护评分的波峰面积、波峰夹角、波峰间距、波峰高度和波峰宽度，总计五个波峰参数进行计算；同时，对持家基因转录起始区域的窗口化保护评分的相同的五个波峰参数进行计算，并将其作为波峰参数阈值；筛选候选染色质开放区域的波峰参数在所述波峰参数阈值范围内的区域，用于后续判断。

本例计算候选染色质开放区域两侧各5个WPS的波峰面积、波峰夹角、波峰间距、波峰高度、波峰宽度等的均值和方差。首先提取区域两侧各5个波峰对象，根据每一个波峰对象的波峰位置、以及其两侧的波谷位置，如图6所示计算该区域两侧各5个WPS的波峰面积、波峰夹角、波峰间距、波峰高度、波峰宽度。其中，计算波峰面积的公式如式二所示，计算波峰夹角的公式如式三所示；

式二

式三

式二中，x为波谷在基因组上的坐标；i为对应的第几个波谷的编号，范围为1-9；y为该碱基位点对应的窗口化保护评分数值，d_x为x的微分；式三中，A为波峰位点，B和C分别为相应波峰位点两侧的波谷位点。

计算该区域两侧各5个WPS的波峰面积、波峰夹角、波峰间距、波峰高度、波峰宽度等数值的均值和方差，其中，平均值的计算公式如式四所示，方差的计算公式如式五所示：

式四

式五

式四和式五中，

为平均值，n＝10，S_d ²为方差。

本例根据1000个持家基因转录起始区域的WPS波形计算阈值范围，当波峰面积的方差小于300、波峰宽度的方差小于450、波峰夹角的方差小于200、波峰间距的夹角小于350以及波峰高度的夹角小于1200这些条件均满足时，认为该区域周边WPS波峰规则，中间缺失至少一个核小体，符合对于开放染色质区域的定义。

染色质开放区域识别步骤，包括对波峰参数筛选步骤获得的候选染色质开放区域，使用该区域的cfDNA平均标准化测序深度进行验证，判断各候选染色质开放区域是否出现cfDNA覆盖度降低的现象，如果出现cfDNA覆盖度降低的现象，则将其识别为染色质开放区域。

本例判断各候选染色质开放区域是否出现cfDNA覆盖度降低的现象，具体包括，计算全基因组随机等长区域内cfDNA的测序深度，当所述cfDNA平均标准化测序深度小于全基因组随机等长区域内cfDNA的测序深度的n倍时，则将判断该区域出现了cfDNA覆盖度降低的现象，其中n＜0.5。

本例再使用20kb区间范围的cfDNA平均标准化测序深度进行验证，标准化的公式如式一所示，当覆盖度小于全基因组随机等长区域内cfDNA的测序深度的n(n<0.5)倍时，则认为该区域确实出现了cfDNA覆盖度降低的现象，符合开放染色质区域的筛选条件，即识别为染色质开放区域。

本例最终从100个健康人的总计180G的cfDNA全基因组测序数据中识别得到57730个染色质开放区域，与预期相符。

实施例2

本例从Encode数据库下载人体造血细胞的DNase-seq数据，这些数据是目前判断染色质开放区域的金标准。本例采用实施例1相同的染色质开放区域识别方法对下载的人体造血细胞的DNase-seq数据进行染色质开放区域识别，以此评价实施例1的染色质开放区域识别方法的准确性和可靠性。DNase-seq数据用于表示染色体是否处于开放状态，处于开放状态的区域被称作DNase酶高敏位点(DHSs)。以该数据为基础，在所有细胞系中根据DHS信号强度，挑选出代表性DHSs。

造血细胞的DNase-seq数据测序结果给出了203360个的500bp-1500bp的大区间和约227205个150bp的小区间，这些区间为可能的染色质开放区域，由于参考集合较大，本例根据DHS信号强度进行筛选，当认为DNase-seq测序数据的峰值高度-log₁₀(Pvalue)≥800时对应的染色质开放区域比较可信。

本例采用实施例1的方法检测出的染色质开放区域有47.3％和实验结果重复。说明实施例1的染色质开放区域识别方法适用于正常人血浆cfDNA的60×以上的测序数据，能够识别50％的造血谱系的染色质开放区域，在测序深度提高的情况能识别更多的染色质开放区域。

实施例3

基因的转录起始位点周边区域一般都存在启动子、增强子以及转录因子结合位点等与表达相关的DNA区域，当基因处于活跃表达状态时，会出现染色质开放区域。

持家基因是在所有细胞中都会稳定表达的基因，因此本例首先分析持家基因的转录起始点(TSS)周边区域和实验结果的重复程度。

本例从持家基因数据库http://www.tau.ac.il/～elieis/HKG/下载3444个持家基因的位置信息，对比实施例1检测出的57730个染色质开放区域集合，结果显示有57％的持家基因TSS区域和实施例1检测出的染色质开放区域有交集，并且实施例1检测出的染色质开放区间中心和持家基因的TSS的距离的分布服从均值为0的类正态分布，如图7所示。本例分析了全体基因的TSS区域和我们的结果的重合程度，结果显示，实施例1检测出32.7％的人体基因的TSS区域含有染色质开放区域，可能处于基因表达状态。

实施例4

本例从数据库http://bioinf.xmu.edu.cn/PaGenBase/advance.jsp下载组织特异性基因TSS区域坐标信息，对比分析了实施例1方法搜索的染色质开放区域和组织特异性基因TSS区域的重叠百分比。

结果如图8所示，造血系统的组织特异性基因的TSS重叠率明显高于其他组织，这与正常人群的cfDNA主要来自造血系统的认识是一致的。因此，可以解释实施例1检测到的染色质开放区域具有一定的组织特异性，并且具有作为组织追踪方法的能力。

以上内容是结合具体的实施方式对本申请所作的进一步详细说明，不能认定本申请的具体实施只局限于这些说明。对于本申请所属技术领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干简单推演或替换。

Claims (11)

1.一种基于测序数据识别染色质开放区域的方法，其特征在于：所述测序数据为cfDNA全基因组测序数据，所述方法包括以下步骤，

测序数据处理步骤，包括根据所述cfDNA全基因组测序数据在参考基因组上的位置，结合区间划分和滑动窗口，计算窗口化保护评分和测序深度；

初步筛选步骤，包括定位窗口化保护评分的波峰和波谷，将相邻波峰的间距大于阈值的区域作为候选染色质开放区域，所述阈值大于或等于核小体大小；定位窗口化保护评分的波峰和波谷，具体包括，使用三个大小一致的小窗口，分别记为窗口1、窗口2、窗口3，三个窗口依次间隔5bp，滑动三个窗口，并计算三个窗口的窗口化保护评分数值之和，分别记为sum_win1、sum_win2、sum_win3；若sum_win1<sum_win2，且sum_win3<sum_win2，则为波峰；若sum_win1>sum_win2，且sum_win3>sum_win2，则为波谷，所述小窗口的大小为4-7bp；

波峰参数筛选步骤，包括对每个候选染色质开放区域两侧各3-9个窗口化保护评分的波峰面积、波峰夹角、波峰间距、波峰高度和波峰宽度，总计五个波峰参数进行计算；同时，对持家基因转录起始区域的窗口化保护评分的相同的五个波峰参数进行计算，并将其作为波峰参数阈值；筛选候选染色质开放区域的波峰参数在所述波峰参数阈值范围内的区域，用于后续分析；

染色质开放区域识别步骤，包括对所述波峰参数筛选步骤获得的候选染色质开放区域，使用该区域的cfDNA平均标准化测序深度进行验证，判断各候选染色质开放区域是否出现cfDNA覆盖度降低的现象，如果出现cfDNA覆盖度降低的现象，则将其识别为染色质开放区域。

2.根据权利要求1所述的方法，其特征在于：所述染色质开放区域识别步骤中，判断各候选染色质开放区域是否出现cfDNA覆盖度降低的现象，具体包括，计算全基因组随机等长区域内cfDNA的测序深度，当所述cfDNA平均标准化测序深度小于全基因组随机等长区域内cfDNA的测序深度的n倍时，则将判断该区域出现了cfDNA覆盖度降低的现象，其中n＜0.6。

3.根据权利要求2所述的方法，其特征在于：测序深度标准化的方法，具体包括，将参考基因组按m kb大小划分为不同的区间，计算m kb区间cfDNA的覆盖度，并进行标准化，具体公式如式一所示，

式一

其中，x^*为标准化后的值，x为初始值，x_min为区间最小值，x_max为区间最大值；10kb≤m≤50kb。

4.根据权利要求1所述的方法，其特征在于：所述测序数据处理步骤中，计算窗口化保护评分和测序深度具体包括，将参考基因组按mkb大小划分为不同的区间，每个区间再以nbp大小为滑动窗口，计算每个窗口的cfDNA读长覆盖度和窗口化保护评分；10kb≤m≤50kb，120bp≤n≤150bp。

5.根据权利要求1所述的方法，其特征在于：所述测序数据处理步骤，还包括使用Savitzky-Golay滤波器消除窗口化保护评分的波形噪声，Savitzky-Golay滤波器的输入数据为所述测序数据处理步骤获得的一维的窗口化保护评分波形数组，数组的索引为参考基因组位置，数组存储的内容为当前基因组位置计算得到的窗口化保护评分数值；Savitzky-Golay滤波器通过kbp大小的滑动窗口，通过滑动窗口内的所有点得到局部多项式回归方程，最后计算窗口中点的预测值，作为每个碱基位点的窗口化保护评分的校正值；25bp≤k≤50bp。

6.根据权利要求1所述的方法，其特征在于：还包括对波峰或波谷的可靠性进行判断，具体包括计算波峰波谷的波峰高度、半波形宽度、相邻波峰宽度、相对波谷宽度、相邻波峰相对高度，总计五个波峰数值；同时，对持家基因转录起始区域的窗口化保护评分的相同的五个波峰数值进行计算，并将其作为波峰数值阈值；波峰波谷的五个波峰数值在所述波峰数值阈值范围内，则认为该位置为可靠的波峰或者波谷。

7.根据权利要求1所述的方法，其特征在于：所述初步筛选步骤中，阈值为250bp，即将相邻波峰的间距大于250bp的区域作为候选染色质开放区域。

8.根据权利要求1-7任一项所述的方法，其特征在于：所述波峰参数筛选步骤中，计算波峰面积的公式如式二所示，计算波峰夹角的公式如式三所示；

式二 波峰面积

式三 波峰夹角

式二中，x为波谷在基因组上的坐标；i为对应的第几个波谷的编号，范围为1-9；y为波峰范围的碱基位点对应的窗口化保护评分数值，d_x为x的微分；式三中，A为波峰位点，B和C分别为相应波峰位点两侧的波谷位点。

9.一种基于测序数据识别染色质开放区域的装置，其特征在于：所述装置用于基于cfDNA全基因组测序数据进行染色质开放区域识别，所述装置包括测序数据处理模块、初步筛选模块、波峰参数筛选模块和染色质开放区域识别模块；

所述测序数据处理模块，包括用于根据所述cfDNA全基因组测序数据在参考基因组上的位置，结合区间划分和滑动窗口，计算窗口化保护评分和测序深度；

所述初步筛选模块，包括用于定位窗口化保护评分的波峰和波谷，将相邻波峰的间距大于阈值的区域作为候选染色质开放区域，所述阈值大于或等于核小体大小；定位窗口化保护评分的波峰和波谷，具体包括，使用三个大小一致的小窗口，分别记为窗口1、窗口2、窗口3，三个窗口依次间隔5bp，滑动三个窗口，并计算三个窗口的窗口化保护评分数值之和，分别记为sum_win1、sum_win2、sum_win3；若sum_win1<sum_win2，且sum_win3<sum_win2，则为波峰；若sum_win1>sum_win2，且sum_win3>sum_win2，则为波谷，所述小窗口的大小为4-7bp；

所述波峰参数筛选模块，包括用于对每个候选染色质开放区域两侧各5个窗口化保护评分的波峰面积、波峰夹角、波峰间距、波峰高度和波峰宽度，总计五个波峰参数进行计算；同时，对持家基因转录起始区域的窗口化保护评分的相同的五个波峰参数进行计算，并将其作为波峰参数阈值；筛选候选染色质开放区域的波峰参数在所述波峰参数阈值范围内的区域，用于后续判断；

所述染色质开放区域识别模块，包括用于对所述波峰参数筛选模块获得的候选染色质开放区域，使用该区域的cfDNA平均标准化测序深度进行验证，判断各候选染色质开放区域是否出现cfDNA覆盖度降低的现象，如果出现cfDNA覆盖度降低的现象，则将其识别为染色质开放区域。

10.一种基于测序数据识别染色质开放区域的装置，其特征在于：所述装置包括存储器和处理器；

所述存储器，包括用于存储程序；

所述处理器，包括用于通过执行所述存储器存储的程序以实现权利要求1-8任一项所述的方法。

11.一种计算机可读存储介质，其特征在于：所述存储介质中包括程序，所述程序能够被处理器执行以实现权利要求1-8任一项所述的方法。

Images