CN107622183B - 一种基于多重指标的胎儿染色体倍性检测分析方法 - Google Patents

Abstract

本发明公开的一种基于多重指标的胎儿染色体倍性检测分析方法，其包括如下步骤：(1)质量控制步骤；(2)胎儿DNA含量估计步骤；(3)计算Z、ZZ、BM统计指标步骤；(4)输出图表步骤；(5)结果评估步骤。本发明与现有技术相比，主要有以下几大优势：(1)质控更加严格。(2)检测指标多样化。(3)检测结果更全面。(4)数据结果可视化。

Description

一种基于多重指标的胎儿染色体倍性检测分析方法

技术领域

本发明涉及生物检测技术领域,特别涉及一种基于多重指标的胎儿染色体倍性检测分析方法。

背景技术

我国二孩政策于2016年1月正式实施后，据悉，2016年新生儿数量达到了1786万人。预计“十三五”期间年出生人口在1700万-2000万人之间。而我国是出生缺陷高发国家之一，出生缺陷总发生率约为5.6％，随着二孩政策放开，高龄高危孕产妇比例明显增高，无创产前基因检测(NIPT)市场有望迎来高速增长。

无创产前检测(NIPT)方法相比传统方法主要有以下几方面优势：(一)仅仅抽取母体血液进行检测，对母体和胎儿无伤害，大大降低了由于传统方法检测所可能造成的流产风险；(二)NIPT能在6周左右就检测出胎儿染色体是否正常从而更早地进行诊断及时干预，提高优生率。

目前市场上所有NIPT的应用中，应用最多最成熟的当属胎儿染色体非整倍性检测。对于胎儿染色体非整倍性检测大多采用Z值进行分析，但是由于生物学或者技术上的问题，现有的分析方法仍然存在一定比例的假阳性和假阴性问题。因此，仍需要对现有的检测分析方法进行改进，以提高检测的灵敏度，降低错误率。

NIPT胎儿染色体倍性分析目前存在以下几点问题：

(1)胎儿DNA含量：现有的分析方法中大多没有考虑胎儿DNA含量，有研究表明当胎儿DNA含量低于4％时候，将会大大增加假阴性结果。

(2)母源性干扰：由于母体的染色体拷贝数变异、嵌合、血癌等会掩盖或者模拟来自胎儿的真实染色体倍性变异，因此会一定程度上会干扰数据分析。

(3)检测指标：现有的分析方法中大多只通过一个Z值进行判断。这样会存在一定的假阳性和假阴性，且存在一定范围的灰区难以判定。

(4)报告结果:现有分析方法中，结果中只有一些简单的图表报告，还有很多数据信息没有有效的呈现。且现有的图表在可视化、形象化角度来说还可以进一步提升。

发明内容

本发明所要解决的技术问题在于针对上述现有技术所存在的不足而提供一种基于多重指标的胎儿染色体倍性检测分析方法。

本发明所要解决的技术问题可以通过以下技术方案来实现：

一种基于多重指标的胎儿染色体倍性检测分析方法，包括如下步骤：

(1)质量控制步骤；

(2)胎儿DNA含量估计步骤；

(3)计算Z、ZZ、BM统计指标步骤；

(4)输出图表步骤；

(5)结果评估步骤。

在本发明的一个优选实施例中，所述质量控制步骤是指将对原始测序数据质量进行评估，并降低由于测序因素造成的数据噪声，提高数据的有效性、分析的准确性，具体步骤包括：

(1.1)利用FastQC查看测序质量结果，查看序列测序得分值分布、GC分布和重复率；统计所有序列的总Total Reads，Total bases、Q20、Q30、GC含量、N字符的数量及相关比例；

(1.2)去除序列中含有的接头序列；

(1.3)对序列首尾的index或者测序质量比较差的reads进行trimming；

(1.4)利用滑动窗口，根据得分值对低质量的Reads进行过滤，如去除N或者低得分值的序列。

在本发明的一个优选实施例中，所述胎儿DNA含量估计步骤是指利用seqFF方法采用加权秩选择标准模型(Weighted Rank Selection Criterion,WRSC)和弹性网络模型(Elastic net,Enet)模型进行估算。

(2.1)通过弹性网络模型估算胎儿DNA含量，计算方法如公式(1)所示；

对n个常染色体的第ith个50kb基因组bin计算β，δ是弹性网络模型中的一个矫正参数。X是bin i的标准化计数；Ye是估算的胎儿DNA比例；

(2.2)通过加权秩选择标准模型估算胎儿DNA含量，计算方法如公式(2)所示；

Yw＝XCD+E (2)

Yw是m个样本和n个Y染色体bins组成的m×n矩阵，X是m个样本和p个常染色体bins组成的m×p矩阵，E为误差矩阵，C和D是p×r和r×n的满秩矩阵，其中r是未知的，需要通过训练数据集来获得；令A＝CD,则A为秩r的降秩矩阵，该矩阵的回归系数可由公式(3)计算：

其中Γ是对应协方差Y的n×n的权重矩阵，μ是一个调整参数，r(B)是矩阵B的秩；

(2.3)最终胎儿DNA(ff)的含量由公式(4)计算；

ff＝(Ye+Yw)/2 (4)。

在本发明的一个优选实施例中，所述计算Z、ZZ、iZ统计指标步骤：

(3.1)将基因组划分为等长的窗口，过滤掉包含N的窗口，并统计各窗口的GC含量；

(3.2)对比对结果进行过滤，获得无重复的比对结果，然后，根据基因组窗口信息，统计单个样本比对到每个窗口的reads数量及GC含量，并得到单个样本比对到每条染色体上的reads数量；

(3.3)根据基因组各个窗口的GC含量和reads数量，进行loess回归分析，校正每个窗口的reads数量，得到每条染色体校正后的reads数量；

Kmer覆盖度统计模块得到单个窗口的读段数，记做UR，对应窗口的GC含量记为GC_bin，使用LOESS回归分析(UR vs.GC_bin)得到的拟合值为UR_loess，对应染色体各窗口读段平均值为

则该窗口校正后的读段数UR_correct为kmer覆盖度：

每条染色体上覆盖度校正值为各窗口校正值之和：

j为染色体，b为j染色体窗口数量；然后计算每条染色体相对覆盖度，

(3.4)对同一个样品内每一条染色体的所有bin计算的Z值，用于检测每条染色内部的差异；通过公式(5)中计算的UR_loess计算各条染色体的所有bin的iZscore：

BMD(Bin Median):每条染色体的iZscore的中位数，该数值不受分布数据的极大或极小值影响，可以有效的代表染色体的iZscore值的分布；

BMA(Bin Mean):反映整条染色体的iZscore的平均值；

(3.5)得到各样本每条染色体的reads相对覆盖度，并根据各个样本的GC含量，对染色体覆盖度进行线性拟合，得到残差，用于z检验，计算z分值；

利用loess回归模型得到染色体相对覆盖度，可直接进行z检验，计算

同时还可以根据染色体相对覆盖度和GC含量进行线性拟合，得到斜率a和截距b，计算残差ε_i＝CR-(a×GC_i+b)，i为样本号，利用残差进行t检验得到z值，

(3.6)对同一样品内的不同染色体的Z值再次计算Z值，用于检测不同染色体间的差异；通过公式(6)得到所有染色体的拟合Z_score，进一步计算所有染色体的ZZ_score：

在本发明的一个优选实施例中，所述输出图表步骤：

(4.1)输出iZscore的分布的箱线图和折线图

(4.2)输出Zscore、染色体的相对覆盖度分布图

在本发明的一个优选实施例中，所述结果评估步骤是指:根据表(1)判断结果所属类型。

本发明与现有技术相比，主要有以下几大优势：

(1)质控更加严格：首先对测序结果进行严格的质控，其次对胎儿DNA含量进行评估，只有质控合格后才进行后续分析，以提高准确性。

(2)检测指标多样化：本方法结合Z、ZZ、BMD、BMA四个统计指标对结果进行判断，结果更加可靠准确。

(3)检测结果更全面：本方法不仅仅可以分析常规的13、18、21三条染色体的倍性检测，还可以对其它20条染色体的倍性进行检测。

(4)数据结果可视化:本软件分析结果中除了提供有效的数据表格文件外，还生成了形象直观的信息图，使得更多的数据信息得以展示，也使得结果一目了然。

附图说明

图1为本发明的主要分析流程概览示意图。

图2为本发明的iZscore的分布的箱线示意图。

图3为本发明的iZscore的分布的折线示意图。

图4为本发明的染色体的相对覆盖度分布示意图。

图5为本发明的Zscore分布示意图。

具体实施方式

为了实现本发明目的，本发明包括五大主要步骤:(1)质量控制步骤；(2)胎儿DNA含量估计步骤；(3)计算Z、ZZ、BM统计指标步骤；(4)输出图表步骤；(5)结果评估步骤。具体如图1所示。

(1)质量控制步骤

(1.1)利用FastQC查看测序质量结果，查看序列测序得分值分布、GC分布和重复率；统计所有序列的总包括总reads数目，总碱基数目、序列长度、Q20、Q30、GC含量、N数量、Q20％、Q30％、GC％含量、N％。

(1.2)利用AdapterRemoval去除序列5`和3`端含有的接头序列。

(1.3)本步骤利用程序Seqtk，对序列首尾的index或者测序质量比较差的序列进行trimming。

(1.4)利用滑动窗口，根据得分值对低质量的碱基序列进行过滤，如去除N或者低得分值的序列。

(2)胎儿DNA含量估计步骤

(2.1)通过弹性网络模型估算胎儿DNA含量，计算方法如公式(1)所示。

对n个常染色体的第ith个50kb基因组bin计算β，δ是弹性网络模型中的一个矫正参数。X是bin i的标准化计数。Ye是估算的胎儿DNA比例。

(2.2)通过加权秩选择标准模型估算胎儿DNA含量，计算方法如公式(2)所示。Yw是m个样本和n个Y染色体bins组成的m×n矩阵，X是m个样本和p个常染色体bins组成的m×p矩阵，E为误差矩阵，C和D是p×r和r×n的满秩矩阵，其中r是未知的，需要通过训练数据集来获得。令A＝CD,则A为秩r的降秩矩阵，该矩阵的回归系数可由公式(3)计算：

Yw＝XCD+E (2)

其中Γ是对应协方差Y的n×n的权重矩阵，μ是一个调整参数，r(B)是矩阵B的秩。

(2.3)最终胎儿DNA(ff)的含量由公式(4)计算。本实施案例中计算的ff值见表(2)。

ff＝(Ye+Yw)/2 (4)

(3)计算Z、ZZ、BM统计指标步骤

(3.1)将基因组划分为等长的窗口，过滤掉包含N的窗口，并统计各窗口的GC含量。

(3.2)对比对结果进行过滤，获得无重复的比对结果，然后，根据基因组窗口信息，统计单个样本比对到每个窗口的reads数量及GC含量，并得到单个样本比对到每条染色体上的reads数量。

(3.3)根据基因组各个窗口的GC含量和reads数量，进行loess回归分析，校正每个窗口的reads数量，得到每条染色体校正后的reads数量。

Kmer覆盖度统计模块得到单个窗口的读段数，记做UR，对应窗口的GC含量记为GC_bin，使用LOESS回归分析(UR vs.GC_bin)得到的拟合值为UR_loess，对应染色体各窗口读段平均值为

则该窗口校正后的读段数UR_correct为kmer覆盖度：

每条染色体上覆盖度校正值为各窗口校正值之和：

j为染色体，b为j染色体窗口数量。然后计算每条染色体相对覆盖度，

(3.4)对同一个样品内每一条染色体的所有bin计算的Z值，用于检测每条染色内部的差异。通过(3.3)中计算的UR_loess计算各条染色体的所有bin的iZscore，公式如下：

BMD(Bin Median):每条染色体的iZscore的中位数，该数值不受分布数据的极大或极小值影响，可以有效的代表染色体的iZscore值的分布。

BMA(Bin Mean):反映整条染色体的iZscore的平均值。

(3.5)得到各样本每条染色体的reads相对覆盖度，并根据各个样本的GC含量，对染色体覆盖度进行线性拟合，得到残差，用于z检验，计算z分值。

利用loess回归模型得到染色体相对覆盖度，可直接进行z检验，计算

同时还可以根据染色体相对覆盖度和GC含量进行线性拟合，得到斜率a和截距b，计算残差ε_i＝CR-(a×GC_i+b)，i为样本号，利用残差进行t检验得到z值，

(3.6)对同一样品内的不同染色体的Z值再次计算Z值，用于检测不同染色体间的差异。通过(3.4)得到所有染色体的拟合Z_score，进一步计算所有染色体的ZZ_score。

本实施案例中以上各统计指标结果示例见表1、表2。

表1

表2

(4)输出图表步骤；

(4.1)输出iZscore的分布的箱线图和折线图，结果示例见图2和图3。

(4.2)输出染色体的Zscore、相对覆盖度分布图，结果示例见图4和图5。

(5)结果评估步骤

在本实施例中，步骤(3)已得出相关的统计结果，见表3

表3

根据表3判断结果所属类型。如果Zscore在[-3,3]范围内，ZZscore在[-2.5,2.5],则可以初步判定该样本为正常样本，若BMD、BMA都小于0则更加能充分说明该样本为正常样本。若Zscore、ZZscore的值都大于3，则初步判定为Triomy，若BMA、BMD大于0则能进一步判断为Triomy。同理可判断Monosomy。其它条件则可判定为Undermined。本实施例中S1、S2、S3分别为13、18、21阳性。参见表2

表2

Samples	S1	S2	S3	S4	S5
						ff values(％)	12.08	10.82	11.83	11.31	12.13

Claims (2)

1.一种基于多重指标的胎儿染色体倍性检测分析方法，其特征在于，包括如下步骤：

(1)质量控制步骤；

(2)胎儿DNA含量估计步骤；

(3)计算Z、ZZ、iZ统计指标步骤；

(4)输出图表步骤；

(5)结果评估步骤；

所述质量控制步骤是指将对原始测序数据质量进行评估，并降低由于测序因素造成的数据噪声，提高数据的有效性、分析的准确性，具体步骤包括：

(1.1)利用FastQC查看测序质量结果，查看序列测序得分值分布、GC分布和重复率；统计所有序列的Total Reads、Total bases、Q20、Q30、GC含量、N字符的数量及相关比例；

(1.2)去除序列中含有的接头序列；

(1.3)对序列首尾的index或者测序质量比较差的reads进行trimming；

(1.4)利用滑动窗口，根据得分值对低质量的reads进行过滤，如去除N字符或者低得分值的序列；

所述胎儿DNA含量估计步骤是指利用seqFF方法采用加权秩选择标准模型和弹性网络模型进行估算；

(2.1)通过弹性网络模型估算胎儿DNA含量，计算方法如公式(1)所示；

对n个常染色体的第i个50kb基因组bin计算β，

δ是弹性网络模型中的一个矫正参数；

X是bin i的标准化计数；

Ye是估算的胎儿DNA比例；

(2.2)通过加权秩选择标准模型估算胎儿DNA含量，计算方法如公式(2)所示；

Yw＝XCD+E (2)

Yw是m个样本和n个Y染色体bins组成的m×n矩阵，

X是m个样本和p个常染色体bins组成的m×p矩阵，

E为误差矩阵，C和D是p×r和r×n的满秩矩阵，其中r是未知的，需要通过训练数据集来获得；令A＝CD,则A为秩r的降秩矩阵，该矩阵的回归系数可由公式(3)计算：

其中Γ是对应协方差Y的n×n的权重矩阵，μ是一个调整参数，r(B)是矩阵B的秩；

(2.3)最终胎儿DNA的含量(ff)由公式(4)计算；

ff＝(Ye+Yw)/2 (4)；

所述计算Z、ZZ、iZ统计指标步骤：

(3.1)将基因组划分为等长的窗口，过滤掉包含N的窗口，并统计各窗口的GC含量；

(3.2)对比对结果进行过滤，获得无重复的比对结果，然后，根据基因组窗口信息，统计单个样本比对到每个窗口的reads数量及GC含量，并得到单个样本比对到每条染色体上的reads数量；

(3.3)根据基因组各个窗口的GC含量和reads数量，进行loess回归分析，校正每个窗口的reads数量，得到每条染色体校正后的reads数量；

Kmer覆盖度统计模块得到单个窗口的读段数，记做UR，对应窗口的GC含量记为GC_bin，

使用LOESS回归分析得到的拟合值为UR_loess，对应染色体各窗口读段平均值为

则该窗口校正后的读段数UR_correct为kmer覆盖度：

每条染色体上覆盖度校正值为各窗口校正值之和：

j为染色体，b为j染色体窗口数量；

然后计算每条染色体相对覆盖度，

(3.4)对同一个样品内每一条染色体的所有bin计算的Z值，用于检测每条染色体内部的差异；

通过计算的UR_loess计算各条染色体的所有bin的iZscore：

BMD:每条染色体的iZscore的中位数，该数值不受分布数据的极大或极小值影响，可以有效的代表染色体的iZscore值的分布；

BMA:反映整条染色体的iZscore的平均值；

(3.5)得到各样本每条染色体的reads相对覆盖度，并根据各个样本的GC含量，对染色体覆盖度进行线性拟合，得到残差，用于z检验，计算z分值；

利用loess回归模型得到染色体相对覆盖度，可直接进行z检验，计算

同时还可以根据染色体相对覆盖度和GC含量进行线性拟合，得到斜率a和截距b，计算残差ε_i＝CR-(a×GC_i+b)，i为样本号，利用残差进行t检验得到z值，

(3.6)对同一样品内的不同染色体的z值再次计算z值，用于检测不同染色体间的差异；

通过所有染色体的拟合Z_score，进一步计算所有染色体的ZZ_score：

所述输出图表步骤包括：

(4.1)输出iZ_score的分布的箱线图和折线图；

(4.2)输出Z_score、染色体的相对覆盖度分布图。

2.如权利要求1所述的种基于多重指标的胎儿染色体倍性检测分析方法，其特征在于，所述结果评估步骤是指:根据步骤4输出的图表判断结果所属类型。

Images