CN109979535B - 一种胚胎植入前遗传学筛查装置 - Google Patents

Abstract

本发明涉及胚胎植入前遗传学筛查的装置。该装置包含数据获取模块、序列比对模块、标准化模块、分类模块、GC校正模块、背景库校正模块、染色体异常分析模块以及结果报告模块。

Description

一种胚胎植入前遗传学筛查装置

技术领域

本发明涉及早期产前筛查领域，尤其涉及一种用于胚胎植入前遗传学筛查的装置及方法。

背景技术

胚胎植入前遗传学筛查(Preimplantation Genetic Screening，PGS)是指胚胎植入着床之前，对早期胚胎进行染色体数目(核型)和结构异常(CNV)的检测，通过一次性检测胚胎23对染色体的结构和数目，分析胚胎是否有遗传物质异常的一种早期产前筛查方法。通过胚胎植入前遗传学筛查，可以挑选正常的胚胎植入子宫，以期获得正常的妊娠，提高患者的临床妊娠率，降低多胎妊娠。

PGS遗传学筛查是针对胚胎所有染色体的筛查，可以查看染色体的对数是否有缺失、染色体的形态结构是否正常等。染色体在细胞分裂之前才形成，因此PGS会在受精卵形成胚胎(培育第3天)，或形成囊胎(培育第5天)后进行筛查。染色体有问题的胚胎很难自然发育到成熟，一般常见的情况是会在第5、6个月停育流产。即便胚胎能够存活到自然生产，未来生育出来的婴儿也极有可能发生健康问题。比如，智力低下、头小、眼距宽、耳位低、短颈、鼻塌而短、外生殖器发育不良、腭裂、肌张低下或亢进、颠痫、肛门闭锁、发育迟缓、眼裂小、持续性新生儿黄疸及明显的青斑、眼睑下垂、心脏畸形、肾脏畸形、虹膜或视网膜缺损等等。

使用基因芯片来进行试管婴儿PGS检测是目前比较常见的方法，市面上检测PGS的基因芯片产品如Illumina的24Sure Arrays,Perkin Elmer的KaryoLiteBoBs^TMPGS检测芯片和Oxford Gene Technology的CytoSure^TMArrays。随着二代测序的普及，实验机构与生物公司也在尝试使用低通量的NGS方法来进行PGS的检测。

基于低通量NGS的核型检测方法中，GC是影响数据波动的主要来源。市面上有类似产品采用了LOWESS回归结合PCA的方法。该方法的优点在于：既打消了GC的影响又能用背景库抵消其余无法解释的数据波动，其缺点在于：LOWESS非参数回归对于数据波动的抗干扰能力较弱，同时也无法适应高度偏好的数据。

基于NGS平台的CNV(Copy number variations，拷贝数变异)检测，可以在保证检测性能的前提下一次性给出多个基因的CNV检测结果，同时对于低纯度肿瘤样本有更好的检测效果。常规的NGS平台CNV检测使用全基因组低深度测序数据，该种技术成本较高，对于CNV断点的确定也较为模糊。随着NGS技术的不断进步，基于目标区域捕获的高深度测序技术在临床检测的应用场景下逐渐表现出了的优势。但相比于全基因组测序数据，此种技术所产生的数据具有一定的捕获偏好性，数据波动较大。目前尚无公认、成熟且稳定的用于此平台的CNV检测方法。

CNV检测方法中，断点的确定直接影响到检测结果的判定。CBS是目前检测CNV断点的主流算法之一。其优点是：分段结果稳定，对低浓度的断点检测灵敏度较高。但也存在如下缺点：对波动较大的数据分段过于琐碎；无法对CNV进行定性判断；参数适应性差，需针对数据调试参数阈值。

发明内容

就检测标准而言，使用基于二代测序数据来分析核型和CNV都是基于深度的分析。就数据的表现形式而言，拷贝数变异和CNV表现在分析结果上都是数量比的差异。因此，有可能在基于NGS平台的胚胎植入前遗传学筛查中，同时进行染色体数目(核型)和结构异常(CNV)的检测。

本发明的在于提供一种对于数据波动的抗干扰能力强、且能够适应高度偏好的数据的同时检测染色体数目和结构异常的装置及方法。

即，本发明包括：

1.一种用于胚胎植入前遗传学筛查的装置，其包含下述模块:

数据获取模块，用于获取对测试样本进行全基因组低深度测序而得到目的测序数据；

序列比对模块，其与所述数据获取模块连接，用于将所述数据获取模块获取的测序数据与参考基因组序列进行比对，其中，将参考基因组(例如HG19)划分为多个窗口，统计每个所述窗口内部的比对到的测序数据，所述比对到的测序数据包括样本的碱基序列(reads)的数量和样本的唯一匹配碱基序列(unique reads)的数量；

标准化模块，其与所述序列比对模块相连接，用于以背景库为参照，对每个所述窗口内部的样本的唯一匹配碱基序列做标准化处理，得到标准化窗口；

分类模块，其与标准化模块相连接，用于根据阈值将所述标准化窗口划分为第一类窗口和第二类窗口，所述第一类窗口的计算值在阈值以上，所述第二类窗口的计算值小于阈值，

其中，所述计算值为每个所述标准化窗口内部的参考基因组的唯一匹配碱基序列的比对率(即，所述计算值为每个所述标准化窗口内部全参考基因组唯一的碱基序列的比例；唯一匹配的参考基因组碱基序列长度与唯一匹配的样本碱基序列的长度相同)；

GC校正模块，其与所述分类模块相连接，包括第一类窗口校正子模块、第二类窗口校正子模块以及GC校正结果合并模块，

所述第一类窗口校正子模块与所述分类模块相连接，用于将所述第一类窗口内部的样本的唯一匹配碱基序列分别按照全染色体和分染色体进行GC校正，得到第一类窗口GC校正结果，

所述第二类窗口校正子模块与所述第一类窗口校正子模块相连接，用于采用所述第一类窗口GC校正结果对所述第二类窗口内部的样本的唯一匹配碱基序列进行模拟，得到第二类窗口GC校正结果，

所述GC校正结果合并模块用于合并所述第一类窗口GC校正结果和所述第二类窗口GC校正结果，得到样本GC校正结果，

其中，所述GC校正采用广义线性模型(GLM模型,General Liner Model)；

背景库校正模块，其与所述GC校正模块相连接，用于消除系统固有的数据波动，获得背景库矫正结果；

染色体异常分析模块，其与所述背景库校正模块和所述GC校正模块相连接，用于分析染色体的异常；

以及

结果报告模块，其与所述染色体异常分析模块相连接，用于报告染色体数目和结构的变异情况。

2.根据项1所述的装置，其中，所述染色体异常分析模块包括染色体数目分析子模块和染色体结构变异分析子模块，

所述染色体数目分析子模块，其与所述背景库校正模块相连接，用于根据所述样本系统校正结果确定染色体数目变异的情况，

所述染色体结构变异分析子模块，其与所述GC校正模块相连接，用于根据所述样本GC校正结果确定存在染色体结构变异的区间。

3.根据项1或2所述的装置，所述GLM模型采用下述式(1)所示公式，

D＝beta0+beta1*GC

D＝log(UR-min(UR)+1)……(1)

ΔD＝D-hat_D

式(1)中，UR表示每个所述第一类窗口中样本的唯一匹配碱基序列的数量，函数min表示取对象的最小值，beta0和beta1为回归系数，GC表示碱基序列中G碱基C碱基的含量，ΔD表示回归方程的残差，hat_D表示D的回归估计值，

4.根据项3所述的装置，其中，GC校正模块的第一类窗口GC校正结果为ΔD，ΔD即D_chri，通过式(2)获得，

ΔD＝{ΔD_chr1，ΔD_chr2…ΔD_chr22，ΔD_chrX，ΔD_chrY}……(2)。

式(2)表示ΔD是一个集合，集合内部的元素是每条染色体的校正结果值。

5.根据项1～4中任一项所述的装置，所述背景库校正模块中采用的所述背景库为阴性背景库，通过建立多元线性回归模型进行背景库校正，获得校正残差ΔY，

所述多元线性回归模型如式(3)所示：

Y＝f(X₁,X₂,X₃,X₄,...,X_n)……(3)

这里，Y是待处理的样本，X₁,X₂,...，X_n是同流程处理得到的阴性样本，残差ΔY如式(4)所示，

ΔY＝{ΔY_{chr1,...,chr22}，ΔY_chrX，ΔY_chrY}……(4)

式(4)中，ΔY是一个集合，集合内部的元素是每条染色体的校正结果值。

其中，将全染色体窗口分为{chr 1，chr 2,...,chr 22},{chrX},{chrY}三部分，分别使用上述多元线性回归模型获得其残差，

ΔY_{chr1,...,chr22}＝ΔY_{chr1,...,chr22}+median(D _chr1_,...,chr22)

ΔY_chrX＝ΔY_chrX+median(D _chrX)

ΔY_chrY＝ΔY_chrY+median(D _chrY)

将分染色体的ΔY的中位数作为估计的该染色体的拷贝数，但是，如果判断不存在Y染色体，则令Y染色体的ΔY＝ΔY-median(ΔY))

ΔY＝ΔY_{chr1,...,chr22}+ΔY_chrX+ΔY_chrY(意为合并1:22号，X，Y染色体四组分类计算的结果值)。

6.根据项2～5中任一项所述的装置，所述染色体结构变异分析子模块包括分段元件、偏离度计算元件、结构变异确认元件；

所述分段元件采用小波分段法，其与所述背景库校正模块相连接，用于分染色体对所述ΔY进行小波变换，得到Seg_i值，下标i取值为所属的染色体号；其中，小波分段法是一种数学算法，能够将随机数分类，用于区分是否发生了拷贝数变异；

所述偏离程度计算元件与所述分段元件相连接，用于基于所述Seg_i值计算每条染色体各窗口的偏离程度(Deviation)；这里，偏离程度可以理解为CNV的信号强弱，越大越可能是个dup，越小越可能是个del，趋于0则认为是正常区间；

所述结构变异确定元件与所述偏离程度计算元件相连接，用于基于所述每条染色体各窗口的偏离程度确定染色体结构变异的区间；这里，偏离度越大越可能是一段CNV。

7.根据项6所述的装置，其中，所述分段元件通过下述公式(5)进行分段：

Seg_i＝haarSeg(ΔY_chri)……(5)

该公式(5)中，haarSeg表示基于R程序语言haarSeg函数的运算结果，ΔY_chri表示打消GC偏好性、消除系统固有波动之后的计算结果，此计算结果用于检测染色体结构变异(如，CNV)。

8.根据项6所述的装置，其中，所述偏离程度计算元件通过下述公式(9)计算中心偏离程度，用于表示CNV信号的强弱，

9.根据项6所述的装置，其中，染色体变异确定元件确定的染色体变异的区间满足下述约束条件：

这里，segCut是一个常数，例如可以为0.3，中心偏离程度Deviation大于segCut认为是CNV，

A是一个常数，例如可为20，minWin＝A表示该数据信号波动的最少窗口数为A个，

B是一个常数，例如可为10，minRun＝B表示不含缺口的最少窗口数为B个，

C是一个常数，例如可为3，maxGap＝C表示相邻两连续窗口最大缺口数为C个。

10.根据项1～9中任一项所述的装置，其中，所述标准化模块通过下述公式(7)将所述唯一匹配测序片段进行标准化:

该公式(10)中，depth为固定的参数，表示预期达到的标准深度，例如可以为200，UR表示测序数据文件中每个窗口中唯一匹配测序片段的计数。式(7)完成了UR的标准化处理，并把处理结果赋值于原始变量UR。

11.根据项1～10中任一项所述的装置，其中，所述比对率mapability可以通过下述公式(9)计算：

M_k(x)＝1/F_k(x)

公式(9)中，F_k(x)表示某长度为k的DNA片段在全基因组中出现的频率。

因此mappability介于0,1之间，越接近于1则表示窗口内唯一程度越高。

具体的，mapability可以通过以下公式(8)计算：

mapa blity＝UR/(maprate+1)……(8)

该公式(8)中，maprate为测序数据中获得的比对率，UR表示测序数据文件中每个窗口中唯一匹配测序片段的计数。

发明效果

根据本发明，能够提供一种同时检测染色体数目和结构异常的装置及方法，该装置及方法具有对于数据波动的抗干扰能力强、且能够适应高度偏好的数据的优点，能够平稳地处理具有较高偏好度的单细胞扩增数据，提高染色体数目和结构异常的检测准确率。

附图说明

图1为显示本发明的用于胚胎植入前遗传学筛查的装置的一例的示意图。

图2为显示实施例步骤2的拆分结果的一例的示意图。图中以Project开头的都是批次数据，分别代表一单分析任务。

图3为显示实施例步骤3获得以窗口为单位的序列数据的一例的示意图。图中方框内标识的就是划分窗口后的文件。

图4为实施例1步骤4获得的核型结果图。

图5为实施例1步骤4获得的CNV结果图。

图6为实施例2对样本2的核型分析结果图。

图7为实施例2对样本2的CNV分析结果图。

发明的具体实施方式

本说明书中提及的科技术语具有与本领域技术人员通常理解的含义相同的含义，如有冲突以本说明书中的定义为准。

实施例1

1、下机文件是以bcl为拓展名的illumina hiseq4000测序仪所获得的测序数据文件；

2、bcl测序文件经过拆分后获得以文库号为单位的批次数据，拆分结果的一例如图2所示。一个批次内部一般会有多个文库，每个文库中包含了这个带检测样本的测序信息。

3、将拆分数据接入预处理模块，划分窗口，获得以窗口为单位的序列数据，其一例如图3所示。

4、读取上述划分窗口的文档的数据，文档数据结构如下。将窗口分类成为第一类窗口以及第二类窗口，用tag表示，是一个逻辑值，由0，1组成。

打消GC偏好性，获得模型残差值一项，用residD表示，第二类窗口的residD值暂时空缺，有第一类窗口的residD值进行模拟。用第一类窗口的residD模拟值作为第二类窗口的residD值，此时residD中的每一个元素都对应一个实数。

用背景库中的阴性样本对待检测样本进行校正，获得结果，记为result，用于计算核型与CNV。result的值即为该窗口染色体拷贝数估计值，某条染色体的拷贝数即为该染色体所有窗口的result值的中位数。核型结果图如图4所示。

CNV计算：

用haarseg函数对result进行分段，获得分段的结果segRes。比较Segmented值与阈值的关系，若高于或低于阈值，则认为可能是CNV。

CNV的Segmented即为偏离程度。

对上述获得的CNV进行筛选，筛选条件包含CNV长度，窗口数等参数。CNV结果如图5所示。常规CNV筛选方法通常将阈值基线设定在通常将阈值基线分别设定单倍体基线和三倍体基线的位置(本实施例中2.0附近)，本实施例的阈值基线与样本检测结果相适应。图5在横坐标约1400位置有一处离群窗口，在常规CNV筛选方法会判定为CNV。通过本实例的染色体结构变异分析子模块，由于其阈值自适应的特性，用于判断CNV的阈值基线适应分段原件根据局部片段的特异性调节阈值的范围，能够有效降低由固定阈值基线引起的CNV的假阳性。

最后，将所获的CNV用OMIM数据库进行注释，获得有意义的，并且长度大于某一标准的CNV，总和到结果报告中。

实施例2

采用与实施例1相同的方法对样本2进行检测。核型和CNV的检测结果分别如图6、图7所示。样本数据集中在二倍体基准线，约2.15附近。若按照常规CNV筛选方法，通常将阈值基线分别设定单倍体基准线和三倍体基准线的位置(本实施例中在2.23和2.05附近)。本实施例的阈值基线与样本检测结果相适应。图7在横坐标约0到150位置有一处离群窗口，在常规CNV筛选方法和本实施例中均判定为CNV；在横坐标约950到970位置有一处离群窗口，在常规CNV筛选方法会判定为CNV，在本实施例中判定为正常波动区域。

工业实用性

根据本发明，提供了一种对于数据波动的抗干扰能力强、且能够适应高度偏好的数据的同时检测染色体数目和结构异常的装置及方法。

Claims (10)

1.一种用于胚胎植入前遗传学筛查的装置，其包含下述模块:

数据获取模块，用于获取对测试样本进行全基因组低深度测序而得到目的测序数据；

序列比对模块，其与所述数据获取模块连接，用于将所述数据获取模块获取的测序数据与参考基因组序列进行比对，其中，将参考基因组划分为多个窗口，统计每个所述窗口内部的比对到的测序数据，所述比对到的测序数据包括样本的碱基序列reads的数量和样本的唯一匹配碱基序列unique reads的数量；

标准化模块，其与所述序列比对模块相连接，用于以背景库为参照，对每个所述窗口内部的样本的唯一匹配碱基序列做标准化处理，得到标准化窗口；

分类模块，其与标准化模块相连接，用于根据阈值将所述标准化窗口划分为第一类窗口和第二类窗口，所述第一类窗口的计算值在阈值以上，所述第二类窗口的计算值小于阈值，其中，所述计算值为每个所述标准化窗口内部的参考基因组的唯一匹配碱基序列的比对率；

GC校正模块，其与所述分类模块相连接，包括第一类窗口校正子模块、第二类窗口校正子模块以及GC校正结果合并模块，

所述第一类窗口校正子模块与所述分类模块相连接，用于将所述第一类窗口内部的样本的唯一匹配碱基序列分别按照全染色体和分染色体进行GC校正，得到第一类窗口GC校正结果，

所述第二类窗口校正子模块与所述第一类窗口校正子模块相连接，用于采用所述第一类窗口GC校正结果对所述第二类窗口内部的样本的唯一匹配碱基序列进行模拟，得到第二类窗口GC校正结果，

所述GC校正结果合并模块用于合并所述第一类窗口GC校正结果和所述第二类窗口GC校正结果，得到样本GC校正结果，

其中，所述GC校正采用广义线性模型GLM,General Liner Model；

背景库校正模块，其与所述GC校正模块相连接，用于消除系统固有的数据波动，获得背景库矫正结果；

染色体异常分析模块，其与所述背景库校正模块和所述GC校正模块相连接，用于分析染色体的异常；

以及

结果报告模块，其与所述染色体异常分析模块相连接，用于报告染色体数目和结构的变异情况。

2.根据权利要求1所述的装置，其中，所述染色体异常分析模块包括染色体数目分析子模块和染色体结构变异分析子模块，其中，

所述染色体数目分析子模块，其与所述背景库校正模块相连接，用于根据所述样本GC校正结果确定染色体数目变异的情况，

所述染色体结构变异分析子模块，其与所述GC校正模块相连接，用于根据所述样本GC校正结果确定存在染色体结构变异的区间。

3.根据权利要求1所述的装置，所述广义线性模型GLM采用下述式(1)所示公式，

式(1)中，UR表示测序数据文件中每个窗口中唯一匹配测序片段的计数，函数min表示取对象的最小值，beta0和beta1为回归系数，GC表示碱基序列中G碱基C碱基的含量，ΔD表示回归方程的残差，hat_D表示D的回归估计值。

4.根据权利要求3所述的装置，其中，所述ΔD为ΔD_chri，通过式(2)获得，

ΔD＝{ΔD_chr1，ΔD_chr2…ΔD_chr22，ΔD_chrX，ΔD_chrY}……(2)

式(2)中，ΔD是一个集合，集合内部的元素是每条染色体的校正结果值。

5.根据权利要求1所述的装置，所述背景库校正模块中采用的所述背景库为阴性背景库，通过建立多元线性回归模型进行背景库校正，获得校正残差ΔY，所述多元线性回归模型如式(3)所示：

Y＝f(X₁,X₂,X₃,X₄,…,X_n)……(3)

这里，Y是待处理的样本，X₁,X₂,X₃,X₄,…,X_n是同流程处理得到的阴性样本，残差ΔY如式(4)所示，

ΔY＝{ΔY_{chr1,...,chr22}，ΔY_chrX，ΔY_chrY}……(4)

式(4)中，ΔY是一个集合，集合内部的元素是每条染色体的校正结果值。

6.根据权利要求5所述的装置，所述染色体结构变异分析子模块包括分段元件、偏离度计算元件、结构变异确定元件；

所述分段元件采用小波分段法，其与所述背景库校正模块相连接，用于分染色体对所述ΔY进行小波变换，得到Seg_i值，下标i取值为所属的染色体号；

所述偏离程度Deviation计算元件与所述分段元件相连接，用于基于所述Seg_i值计算每条染色体各窗口的偏离程度Deviation；

所述结构变异确定元件与所述偏离程度Deviation计算元件相连接，用于基于所述每条染色体各窗口的偏离程度Deviation确定染色体结构变异的区间。

7.根据权利要求6所述的装置，其中，所述分段元件通过下述公式(5)进行分段：

Seg_i＝haarSeg(ΔY)……(5)

该公式(5)中，haarSeg表示基于R程序语言haarSeg函数的运算结果。

8.根据权利要求6所述的装置，其中，所述结构变异确定元件确定的染色体变异的区间满足下述约束条件：

这里，segCut是一个常数，中心偏离程度Deviation大于segCut认为是CNV，

A是一个常数，minWin＝A表示数据信号波动的最少窗口数为A个，

B是一个常数，minRun＝B表示不含缺口的最少窗口数为B个，

C是一个常数，maxGap＝C表示相邻两连续窗口最大缺口数为C个。

9.根据权利要求1所述的装置，其中，所述标准化模块通过下述公式(7)将所述唯一匹配测序片段进行标准化:

该公式(7)中，depth为固定的参数，表示预期达到的标准深度，UR表示测序数据文件中每个窗口中唯一匹配测序片段的计数。

10.根据权利要求1所述的装置，其中，所述比对率mapability通过下述公式(8)计算：

mapability＝UR/(maprate+1)……(8)

该公式(8)中，maprate为测序数据中获得的比对率，UR表示测序数据文件中每个窗口中唯一匹配测序片段的计数。

Images