一种无创产前检测染色体异常的装置
阅读说明：本技术 一种无创产前检测染色体异常的装置 (Noninvasive prenatal chromosome abnormality detection device ) 是由 张静波 曲丽 王伟伟 徐冰 伍启熹 王建伟 刘倩 唐宇 于 2020-12-10 设计创作，主要内容包括：本发明涉及生物信息学技术领域,具体公开了一种无创产前检测染色体异常的装置。所述装置包括：检测模块、数据质控模块、数据预处理模块、计算模块和判断模块；计算模块包括：染色体杂合比计算单元、染色体自适应校正标准差计算单元和染色体Z值计算单元；染色体自适应校正标准差计算单元,根据参考集样本的测序深度与染色体杂合比的标准差的关系进行逻辑回归得到拟合函数,基于所述拟合函数代入待测孕妇样本的测序深度,得到待测样本染色体自适应校正标准差。本装置可以提高染色体异常真阳样本的检出率。(The invention relates to the technical field of bioinformatics, and particularly discloses a noninvasive prenatal detection device for chromosome abnormality. The device comprises: the device comprises a detection module, a data quality control module, a data preprocessing module, a calculation module and a judgment module; the calculation module comprises: a chromosome heterozygosis ratio calculation unit, a chromosome adaptive correction standard deviation calculation unit and a chromosome Z value calculation unit; and the chromosome adaptive correction standard deviation calculating unit is used for performing logistic regression according to the relation between the sequencing depth of the reference set sample and the standard deviation of the chromosome hybridization ratio to obtain a fitting function, and substituting the fitting function into the sequencing depth of the pregnant woman sample to be detected to obtain the chromosome adaptive correction standard deviation of the sample to be detected. The device can improve the detection rate of the chromosome abnormality true positive sample.)
技术领域
本发明涉及生物信息学技术领域，具体地说，涉及一种无创产前检测染色体异常的装置。
背景技术
产前筛查是产科护理的重要组成部分。现行的产前检查方法包括羊膜穿刺术和绒毛膜绒毛取样，这些有创性的胎儿取样都有一定的流产风险。无创产前筛查(NIPT)可以解决有创产前检查的流产风险问题。通常，无创产前筛查利用孕妇外周血中胎儿游离DNA(cfDNA)检测胎儿是否染色体异常。生物信息学算法和工具的快速发展为胎儿染色体异常的检测开辟了新的方向。
基于高通量测序技术方法检测染色异常时，当一个待测孕妇样品经过提取、建库、测序和数据预处理后，需要用生物信息方法对结果数据进行解读。通常，传统NIPT方法使用Z值方法来检测样本胎儿是否染色体异常，即通过比较待检测孕妇样品中校正后的染色体读段数目与参考集样本中同一染色体的读段数目的均值和标准差计算Z值，使用得到的Z值判断待测孕妇样品的胎儿是否有染色体异常。然而，真实临床实践中常常因为参考集测序深度与待测样本测序深度的差异导致假阳性与假阴性的发生。具体地说，由于构建参考集的测序深度与待检测孕妇样品的测序深度存在差异，而构建参考集测序深度的不同会导致染色体读段数标准差变化，从而影响Z值。传统的NIPT方法在计算Z值时没有考虑到参考集的测序深度不同而带来的染色体读段数标准差的差异，而对不同测序深度的待测孕妇样品使用相同的标准差，会导致假阳性与假阴性结果，有很明显的缺点。
如果训练集用测序深度较高的样本训练，而待测样本测序深度较低时，会导致所使用的染色体读段数的标准差偏小，Z值变大，出现假阳(FP)；如果训练集用测序深度较低的样本训练，而待测样本测序深度较高时，会导致所使用的染色体读段数的标准差偏大，Z值变小，出现假阴(FN)。
为解决该问题，现有技术中存在针对不同测序深度的待测样本训练多个相应测序深度的参考集数据，即每次针对特定的测序样本重新训练针对该样本的参考集的染色体杂合比的标准差参数，但该方案存在一些缺陷，主要包含以下几个方面：(1)对于深度很高的样本，如需要训练出相应深度样本的参考集数据，将浪费很多测序成本；(2)每次在检测待测样本的染色体Z值时，如果该样本对应的测序深度没有相应的参考集被训练过，就需要重新构建一个新的参考集，将会花费大量的时间，并且构建出的参考集也只用于一个测试样本的标注差，因此效率很低。
因此，需要提供一种新的无创产前检测染色体异常的装置以解决现有技术的问题。
发明内容
针对上述技术问题，本发明的目的是提供一种可以提高染色体异常真阳样本(TP)的检出率，降低阳性样本的假阴率(FP)的无创产前检测染色体异常的装置。
为了实现本发明的发明目的，本发明的技术方案如下：
一种无创产前检测染色体异常的装置，其特征在于，所述装置包括：检测模块、数据质控模块、数据预处理模块、计算模块和判断模块；
所述计算模块包括：染色体杂合比计算单元、染色体自适应校正标准差计算单元和染色体Z值计算单元；
所述染色体自适应校正标准差计算单元，根据参考集样本的测序深度与染色体杂合比的标准差的关系进行逻辑回归得到拟合函数，基于所述拟合函数代入待测孕妇样本的测序深度，得到待测样本染色体自适应校正标准差，所述拟合函数为R函数：fit＝lm(as.vector(y)～x+I(x^2))，predict(y(test)[[1]]data.frame(x(test)＝num))，其中，fit代表本发明中最终采纳的拟合函数，lm代表R语言中的拟合线性模型，x代表所有参考集样本的测序深度，y代表所有参考集样本的染色体杂合比标准差，x(test)代表待测孕妇样本的测序深度，y(test)代表根据所述拟合函数计算出来的待测孕妇样本的染色体杂合比的标准差。其中R语言中的fit函数用于将参考集样本的总读段数和样本的染色体杂合比的标准差的关系拟合成二次函数((y(fit))～A0+A1x(fit)+A2(x(fit)^2)，其中y(fit)为使用线性回归拟合出来的函数因变量，x(fit)为拟合出来的函数的自变量)并得到该拟合函数的多项式系数(A0，A1，A2)，该拟合函数用于后续的测试样本代入该函数后能够得到自适应读段数标注差；R语言中的predict函数为将上述得到的拟合函数的x值(x(fit))代入测试样本的测序深度(即读段数x(test))得到拟合函数的因变量，即测试样本的染色体杂合比的标准差y(test)。
本发明是一种基于高通量测序数据检测染色体异常的装置，通过全基因组测序技术从DNA中获取文库产物，然后通过高通量测序方法得到其序列，通过预处理和特定的自适应校正确定样本染色体的异常状态。
该装置可根据不同待测样本自身的测序深度信息，利用本发明的线性拟合函数，获得检测结果。本发明只需利用一个参考集样本训练出拟合函数，对于不同的深度的待测样本，无需重新训练参考集样本，只需代入拟合函数即可，在保证染色体杂合比的标准差参数准确的同时节省了大量时间，另外，相比于复杂的函数(例如，多项式抽取理论的随机模型)，本发明仅应用数学模型中简单易懂的线性拟合函数，既降低了过拟合风险，又提高了检测速度和准确性，避免在训练参考集样本环节花费大量时间。
本装置不应用于疾病的诊断和治疗。
本发明中，所述数据预处理模块：将通过质控的待测孕妇样本基因组划分成20kb～100kb大小的滑窗，两个相邻窗口之间有10kb～50kb的重叠序列，计算每个滑窗中的读段数RC；
对待测孕妇样本基因组所有滑窗中的读段数RC进行平均数标准化校正；旨在消除各个样本的总读段数不同而引起的每个滑窗中的读段数的波动。
对待测孕妇样本基因组所有滑窗中通过标准化校正后的读段值RNor使用平滑线条法进行GC校正，得到GC校正后的读段数RGc；旨在消除由于基因组上的GC高低所带来的的读段数的偏差。
对待测孕妇样本基因组所有滑窗中通过GC校正后的读段数RGc进行基线校正。旨在消除批次效应带来的读段数偏差。
优选，将通过质控的待测孕妇样本基因组划分成100kb大小的滑窗，两个相邻窗口之间有50kb的重叠序列。
窗宽或重叠序列过小时会导致检测时间过长，窗宽或重叠序列过大时会导致算法检出效果不稳定，50k的重叠序列既有利于增强算法检出效果的稳定性又可以节省算法运行的时间成本。
本发明中，标准化校正后的每个滑窗中的读段数RNor＝RC/(total/BinT)，其中，待测孕妇样本基因组读段总数为total，滑窗总数为T。
本发明中，GC校正通过R命令实现：RGc＝smooth.spline(RNor[ord])。
本发明中，染色体杂合比计算单元：根据一条染色体中所有滑窗中经过数据预处理模块进行预处理后的读段数值计算均值，记为该条染色体杂合比hh；
染色体Z值计算单元：根据待测样本染色体自适应校正标准差，计算每条染色体Z值：
其中，Zscore:待测孕妇样本某条染色体Z值；
Sample hh：待测孕妇样本的某条染色体杂合比；
Reference hh：参考集样品的某条染色体杂合比；
Adaptive sd：基于逻辑回归得到的某条染色体的自适应校正标准差。
本发明中，检测模块：用于对孕妇外周血游离DNA进行高通量测序，获得待测孕妇样本基因组。
本发明中，数据质控模块：用于将测序得到的待测孕妇样本基因组与人类基因组hg19进行比对，大于2M读段数的文件通过质控。
本发明中，判断模块：根据Z值判断染色体是否异常，当Z值大于3时，判定该条染色体为异常样本，当Z值小于3时，判定该条染色体为正常样本。
本发明的有益效果至少在于：
1、本发明基于逻辑回归的方法，根据待检测孕妇样品的实际测序深度自适应学习校正后染色体读段标准差参数，解决了由于参考集样本和待测孕妇样品的测序深度不同导致的染色体读段数标准差的差异问题。且所用函数简单，降低了过拟合风险，提高了检测速度和准确性；
2、本发明提高了无创产前检测中染色体异常真阳样本的检出率，减少了假阴样本的数量。
附图说明
图1为本发明装置进行检测的流程示意图。
图2、图3、图4为22条染色体在不同测序深度下的22条染色体读段数标准差变化图。其中，每幅图代表一条染色体，各图中横坐标为参考集样本的读段数，纵坐标为该条染色体在此读段数下的标准差。
图5为本发明实施例3中PN00G18DE1061测试时，每个窗中的读段数，其中，上图为进行GC校正前每个窗中的读段数，下图为进行GC校正后每个窗中的读段数。
具体实施方式
下面将结合实施例对本发明的优选实施方式进行详细说明。需要理解的是以下实施例的给出仅是为了起到说明的目的，并不是用于对本发明的范围进行限制。本领域的技术人员在不背离本发明的宗旨和精神的情况下，可以对本发明进行各种修改和替换。
下述实施例中所使用的实验方法如无特殊说明，均为常规方法。下述实施例中所用的材料、试剂等，如无特殊说明，均可从商业途径得到。
实施例1拟合函数建立
本实施例对22条染色体在不同测序深度下的22条染色体读段数标准差进行计算，其变化图参见图2、图3、图4。
其中，每幅图代表一条染色体，各图中横坐标为参考集样本的读段数，纵坐标为该条染色体在此读段数下的标准差。从图2-图4中可以看到，不同测序深度下的染色体读段数的标准差变化明显，随着读段数增加(也就是测序深度的增加)，标准差会变小，这说明对于一个样本来说，应该使用与其测序深度相匹配的标准差参数，而不应对于所有样本都使用同一个标准差。
每一幅图中，根据所有参考集样本的读段数(训练集样本中的22条染色体的读段数)和标准差的对应关系(图中的多个点)，本发明线性拟合出读段数和标准差的线性拟合函数(图中的曲线)，即根据参考集样本得到本发明的拟合函数曲线。fit＝lm(as.vector(y)～x+I(x^2))，predict(y(test)[[1]]data.frame(x(test)＝num))，其中，fit代表本发明中最终采纳的拟合函数，lm代表R语言中的拟合线性模型，x代表所有参考集样本的测序深度，y代表所有参考集样本的染色体杂合比标准差，x(test)代表待测孕妇样本的测序深度，y(test)代表根据所述拟合函数计算出来的待测孕妇样本的染色体杂合比的标准差。
本发明根据这个曲线函数来预测待测样本应该使用的标准差。同时可以从图2-图4中看出，本发明的拟合函数很好地拟合了相关的点，该函数是合理的。
具体过程如下：
(1)如图2-图4所示的样本中，本实施例共采用了69,667个NIPT模拟健康样本(无染色体异常或大于10M微缺失重复变异，孕周大于10周)作为参考集样本以用来拟合出图2-图4中所示的线性拟合函数。所用的测序平台为MGI2000平台，通过测序平台测序可以得到这些样本的fq文件，读段长度为50bp，测序深度为0.1～1X，这些样本的读段数在7M～76M之间，对于上述所有的69,667个样本中的每一个样本都进行(2)～(6)步骤中的预处理流程；
(2)对fq文件进行质控筛选：去除PCR重复片段，去除含有连续N碱基(或者连续5个核苷酸的Phred score平均评分小于20)的低质量读段，保留去除上述读段后读段数大于2M的样本；
(3)将每条染色体划分成100kb大小的窗，窗重叠区为50kb，计算读段在每个窗中的读段数，即读段的起始碱基位置均落在滑窗的起始碱基位置内；
(4)对窗中的读段数进行标准化校正，即每个滑窗中的读段数RC除以该样本所有窗中读段数的均值，标准化校正后的每个滑窗中的读段数RNor＝RC/(total/BinT)，其中，待测孕妇样本基因组读段总数为total，滑窗总数为T；
(5)统计每个窗中读段的GC含量，对上述标准化校正后的读段数根据GC含量使用R语言中的平滑样条法进行GC校正；
(6)计算每条染色体上所有窗的经过GC校正后的读段数的均值，记为该条染色体杂合比hh，以此计算22条染色体的杂合比，记为Reference hh；本实施例中针对13、18、21号染色体计算而来的杂合比均值分别为0.0000034，0.0000052，0.0000015，此均值作为参考集样本的杂合比；进一步计算所有样本(69,667个NIPT模拟健康样本)针对13、18、21号染色体的杂合比标准差为0.06009、0.09525、0.012189，此标准差为传统检测方法中所用的SD值，即依据传统检测方法，13号染色体的标准差固定为0.06009、18号染色体的标准差固定为0.09525、21号染色体的标准差固定为0.012189。
(7)根据69,667个参考集样本中的每个样本所有窗中的总读段数，将上述的69,667个样本分为138个区间，间距为0.5M，分别为(7M,7.5M)，(7.5M,8M)…(75.5M,76M)，例如，若某个样本的总读段数为7.25M，则该样本属于区间(7M，7.5M)统计138个区间中每个区间内所有样本，并计算在该区间内的样本的每条染色体杂合比的标准差，这样对于每条染色体可以得到138个标准差。然后以区间的左边界作为横坐标，此区间对应的标准差为纵坐标，组成138个pair，即(7，SD1)，(7.5，SD2)…(75.5，SD138)，即138个点；
(8)对于22条染色体中的每一条染色体来说，根据步骤(7)中的138个点(以划分的读段数区间的左边界为X坐标，以该区间内所有样本的染色体杂合比的标准差为纵坐标)，采用R包中二次函数的线性回归对这些点进行拟合。然后，利用最小二乘方法构造此二次函数多项式系数，记为A0，A1，A2，即可画出图2-图4中各条染色体的线性拟合曲线，所拟合出的二次函数记为(y(fit))～A0+A1x(fit)+A2(x(fit)^2)，其中y(fit)为使用线性回归拟合出来的函数因变量，x(fit)为拟合出来的函数的自变量，A0为拟合函数的多项式系数中的常数项，A1为拟合函数的多项式系数中的一次项系数，A2为拟合函数的多项式系数中的二次项系数；
到此为止，参考集样本的训练过程完成。
本实施例中针对13、18、21号染色体拟合构造出的二次函数多项式系数见表1。
表1
A0
A1
A2
13号染色体
0.0498219876
-0.0016572671
0.0001552737
18号染色体
0.0614557842
-0.0021351633
0.0002116426
21号染色体
0.0164491534
-0.0028875793
0.0002817179
实施例2染色体异常判断
本实施例提供一种采用本发明的无创产前检测染色体异常的装置检测染色体异常的方法。检测流程示意图见图1。具体如下：
(1)以检测模块对待测孕妇样本外周血游离DNA数据进行测序，所用的测序平台为MGI2000平台，通过测序平台测序可以得到fq文件(包含测序序列及序列的测序质量)，读段长度为50bp；
(2)以数据质控模块对fq文件进行质控筛选：去除PCR重复片段，去除含有连续N碱基(或者连续5个核苷酸的Phred score平均评分小于20)的低质量读段，若去除上述片段剩余读段数小于2M，则重新建库测序；
(3)以数据预处理模块将每条染色体划分成100kb大小的窗，窗重叠区为50kb，即(1,100k)，(50k，150k)…,计算读段在每个窗中的读段数，即读段的起始碱基位置均落在滑窗的起始碱基位置内，如果一个读段的起始碱基位置均在窗的起始位置坐标内，则该窗的读段数加1，对于每一个读段重复这个操作，直至所有读段数归属哪个窗都完成；
(4)以数据预处理模块对窗中的读段数进行标准化校正，即每个滑窗中的读段数RC除以该样本所有窗中读段数的均值，标准化校正后的每个滑窗中的读段数RNor＝RC/(total/BinT)，其中，待测孕妇样本基因组读段总数为total，滑窗总数为T；
(5)以数据预处理模块统计每个窗中读段的GC含量，对上述标准化校正后的读段数根据GC含量使用R语言中的平滑样条法进行GC校正；GC校正通过R命令实现：RGc＝smooth.spline(RNor[ord])。
(6)以数据预处理模块对GC校正后的读段数使用加权线性回归函数进行基线校正，加权线性回归函数的权重为该批次的所有样本在各个窗中的读段数标准差的倒数；
(7)以计算模块中的染色体杂合比计算单元计算每条染色体上所有窗的经过预处理后的读段数的均值，记为该条染色体杂合比hh，以此计算22条染色体的杂合比；
(8)以计算模块中的染色体自适应校正标准差计算单元根据实施例1中得到的22条拟合函数及该样本中所有窗的总读段数x(text)，计算出待测样本染色体自适应校正标准差，Adaptivesd＝A0+A1x(text)+A2(x(text)^2))，其中x(text)为待测样本的总读段数(拟合出来的函数的自变量x(fit))，A0，A1和A2为实施例1中拟合出的多项式系数，Adaptivesd为自适应标准差参数(使用线性回归拟合出来的函数因变量y(fit))；
(9)以计算模块中的染色体Z值计算单元计算每条染色体Z值：
其中，Zscore:待测孕妇样本某条染色体Z值；
Sample hh：待测孕妇样本的某条染色体杂合比；
Reference hh：实施例1中参考集样品的某条染色体杂合比；
Adaptive sd：基于逻辑回归得到的某条染色体的自适应校正标准差。
(10)以判断模块判断染色体是否异常：对于每一条染色体，判断步骤(9)中计算出的Z值是否大于3，如果Z值大于3，则判定为染色体异常，若Z值小于3，则判定为染色体正常。
实施例3
本实施例利用实施例2的方法和传统方法对135例染色体异常的真阳样本和3例阴性样本进行了验证。该138例样本的孕周均大于10周，抽取其外周血数据进行验证。
传统方法与本发明实施例2的方法检测流程只在最后一步Z值计算的过程不同，其他处理流程包括测序、计算各个窗宽中的读段数、GC校正、基线校正等都相同(如实施例2所述流程)。
传统方法：
本发明方法：
从上述的传统方法和本发明方法的Z值计算的方法可以看到，传统方法的分母使用固定标准差(本实施例中13号染色体的固定标准差为0.06009、18号染色体的固定标准差为0.09525、21号染色体的固定标准差为0.012189，此标准差来源于实施例1中步骤(6))，即无论测试样本的深度如何，都采用同一个SD；而本发明方法所使用的分母使用自适应SD，即采用实施例1中的参考集样本所训练出来的线性拟合函数根据待测样本的实际测序深度计算出相应的自适应的标准差参数，进而再计算Z值。
样品来源：在无创产前检测中138例孕妇样本外周血数据。
以编号PN00G18DE1061的样本在检测21号染色体是否异常的检测为例具体说明检测过程如下：
(1)对编号为PN00G18DE1061的待测孕妇样本外周血游离DNA数据进行测序，通过测序平台MGI2000测序可以得到fq文件(包含测序序列及序列的测序质量)，读段长度为50bp，样本的总读段数为6,438,936；
(2)对上述生成的fq文件进行质控：去除PCR重复片段，去除含有连续N碱基(或者连续5个核苷酸的Phred score平均评分小于20)的低质量读段，剩余保留下来总读段数为5,129,088，即约为5.1M，读段数大于2M，则无需重新测序，直接进行下面步骤的分析。
(3)将21号染色体划分成100kb大小的窗，窗重叠区为50kb，21号染色体共被划分出963个窗，计算样本在每个窗中的读段数，即读段的起始碱基位置均落在滑窗的起始碱基位置内。以这种方式，可以得到963个窗的读段数，例如该待测样本在窗(46650000，46750000)中的读段数为936，其他窗的读段数不一一列举。
(4)对上述窗中的读段数进行标准化校正，即每个滑窗中的读段数除以该样本所有窗中读段数的均值，例如，该样本的读段数均值为892，对于窗(46650000，46750000)的读段数标准化校正后为936/892＝1.04；
(5)统计每个窗中读段的GC含量，对上述标准化校正后的读段数根据GC含量使用R语言中的平滑样条法进行GC校正，对于该测试样本，每个窗GC校正后的读段数如图5下图所示；每个窗GC校正前的读段数如图5上图所示。
(6)对上述步骤(5)中GC校正后的每个窗中的读段数使用加权线性回归函数进行基线校正，使用的权重为该批次所有样本的窗中读段数标准差倒数；
(7)计算21号染色体上所有窗(共963个窗)中经过上述预处理后读段数的均值，为0.030172，即21号染色体的杂合比为0.030172；
(8)根据实施例1，得到的21号染色体的拟合函数为：(y(fit))～0.0164491534-0.0028875793x(fit)+0.0002817179(x(fit)^2)，将该函数的x(fit)代入该样本的总读段数5.1，计算出待测样本染色体自适应校正标准差为0.00905；
(9)根据实施例1得到的21号染色体的参考集样本的21号染色体的杂合比为0.0000015，根据步骤(7)中得到的21号染色体的杂合比为0.030172，根据步骤(8)中得到的该样本的21号染色体的自适应标准差为0.00905计算21号染色体Z值：Zscore＝(Samplehh-Reference hh)/(Adaptive sd)＝(0.030172-0.0000015)/0.00905＝3.333；
(10)根据步骤(9)中计算出的Z值为3.333，判定编号为PN00G18DE1061的样本的21号染色体异常。
以本发明上述方法获得所有样本的杂合比、自适应标准差、Z值等信息，从而实现染色体异常的判断。具体如下编号样本的杂合比、自适应标准差、Z值见表2。
表2
样本编号
杂合比
自适应标准差
Z值
PN00G18EE6259
0.1803556(Chr13)
0.04550987(Chr13)
3.963(Chr13)
NE00G19EB0338
0.5252057(Chr13)
0.05958767(Chr13)
8.814(Chr13)
PN00G18PA0136
0.2037483(Chr18)
0.05633075(Chr18)
3.617(Chr18)
PN00G18DD9338
0.2300391(Chr18)
0.05694036(Chr18)
4.04(Chr18)
PN00G16AD0047
0.1102264(Chr21)
0.00941864(Chr21)
11.703(Chr21)
PN00G16AC0579
0.2422204(Chr21)
0.00948322(Chr21)
25.542(Chr21)
全部样本使用传统检测方法(固定标准差FSD)和本发明方法(自适应标准差ASD)得到的13号染色体、18号染色体，21号染色体的Z值情况和具体判断结果信息见表3。
表3
根据表3的数据，两种方法样品验证统计结果见表4。
表4
假阴样本数
假阳样本数
传统方法
8
3
本发明方法
5
0
从表3、表4中的结果可知，3例在传统检测方法中被判定为假阴性样本，使用本发明检测判定为真阳性样本。另外3例在传统方法中被判定为假阳性的样本，使用本发明检测被判定为染色体异常真阴样本。综上所述，本发明提出的基于逻辑回归自适应学习参数的模型可以提高染色体异常真阳样本(TP)的检出率，降低阳性样本的假阴率(FP)样本数。
虽然，上文中已经用一般性说明及具体实施方案对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。