测序数据的处理方法和处理装置的制造方法

测序数据的处理方法和处理装置的制造方法
【专利摘要】本发明提供了一种测序数据的处理方法和处理装置。该处理方法包括：通过高通量测序获取来源于母体外周血样本的核苷酸序列信息；将参考基因组划分成多个特异性区域，每个特异性区域内NRSc值相等；将来源于母体外周血样本的所有染色体的核苷酸序列信息分配至参考基因组的多个特异性区域，统计样本在每个特异性区域内的NRSs值；利用GC含量修正样本每个特异性区域内的NRSs值，记为NRSs'值；基于NRSs'值分别统计目标染色体和对照染色体上所有特异性区域的NRSs'值的均值，记为第一均值和第二均值；将第一均值与第二均值进行差异性检验，根据差异检验结果确定染色体是否存在非整倍性。该处理方法提高对测序数据处理的准确性。
【专利说明】
测序数据的处理方法和处理装置
技术领域
[0001] 本发明涉及测序数据处理领域，具体而言，涉及一种测序数据的处理方法和处理 装置。
【背景技术】
[0002] 染色体异常可能是数目上的或结构上的。数量异常，包括三体性(多一个染色体）、 单体性(丢一个染色体)和多倍性(整个多一套染色体）。结构异常包括由染色体断裂等引起 的结构性重排，如易位、翻转、缺失和插入。
[0003] 染色体数量异常，如非整倍性和多倍性，与多种疾病包括出生缺陷及癌症有关。我 国每年新生儿近2000万，其中约4 %~6 %存在出生缺陷，其中胎儿染色体异常是临床最常 见的出生缺陷类型之一，据统计约160例新生儿中就有1例是染色体异常患者。染色体三体 综合征是染色体疾病中发病率最高的一类，当细胞内某染色体的数目不是正常的两条而是 三条，也即染色体总数目为47条时就会导致三体综合征。最常见的染色体三体综合征有:21 三体综合征（T21 )、18三体综合征（T18)和13三体综合征（T13)。为减少出生缺陷婴儿的比 例，对染色体非整倍性的快速、准确的检测是必要的。
[0004] 超声波扫描技术或生物化学标记物筛查的无创方法，已用于进行染色体异常的风 险判断，但该方法准确率较低，仅为60-80%，且受孕龄等生理因素的影响。而常规的产前诊 断方法则需通过侵入性方法如羊膜穿刺或绒毛膜绒毛取样，因此存在流产风险，且检测周 期较长。1997年，在母体血浆中发现了循环的无细胞胎儿DNA(Lancet.l997 Aug 16;350 (9076):485-7·Presence of fetal DNA in maternal plasma and serum.Lo YMl, Corbetta N，Chamberlain PF，Rai V,Sargent IL,Redman CW,Wainscoat 几.）。1999年，发 现怀有21号染色体三体胎儿的妇女血浆中循环胎儿DNA的浓度明显高于怀有整倍体胎儿妇 女血浆中循环胎儿DNA的浓度(Lo，Y.M.D.et al.，Clin Chem 45:1747-1751(1999) ;Zhong， X.Y.et al^Prenat Diagn 20:795-798(2000))。上述发现为无创产前诊断提供了新的可 能性。在此基础上，产前无创领域取得了诸多进展，如通过甲基化敏感酶富集胎儿DNA以降 低母体背景干扰(PCT/US2004/033175 2004.10.08);通过PCR比对特异基因片段的Ct值以 筛查21号三体(CN200610003103.9，2006.02.10);通过基于RNA-SNP的等位基因扩增检测推 断胎儿染色体非整倍性（0呢00680007354.2,2006.03.17)。然而对胎儿0财的富集耗时耗 力，且扩增技术要求序列的特异性或基因的杂合性，使其难以成为一种通用的技术。
[0005] 2008年，Rossa W.K.Chiu等人提出测序手段可获得外周血中核酸分子的大量信息 (Rossa W.K.Chiu,et al.Noninvasive prenatal diagnosis of fetal chromosomal aneuploidy by massively parallel genomic sequencing of DNA in maternalplasma·PNAS，2008，105:20458-20463)，并发现，在临床意义染色体上有异常的样 本中，其临床意义染色体有异常的核酸分子的量与背景染色体的核酸分子的量的比值参 数，与由正常样本所建得的一个或多个正常对照值的参数存在差异。由此，基于高通量测序 的方法可用来检测染色体异常，并消除了对特异序列扩增的依赖。但现有的基因组测序检 测方法需要将待测样本与多个样本或标准正常样本进行比对分析，耗时长，对样本需求量 大（如，申请号为CN200880108377.1的中国专利申请），且对各批次样本实验条件的一致性 有严格要求，制约了其便捷化与高通量的应用。
[0006] 因此，仍需要对现有的处理测序数据的方法进行改进，以提高数据处理的准确性。

【发明内容】

[0007] 本发明的主要目的在于提供一种测序数据的处理方法和处理装置，以提高对测序 数据处理的准确性。
[0008] 为了实现上述目的，根据本发明的一个方面，提供了一种测序数据的处理方法，该 处理方法包括:通过高通量测序获取来源于母体外周血样本的所有染色体的核苷酸序列信 息;将参考基因组划分成多个特异性区域，每个特异性区域内非重复序列的数目NRSc相等； 将来源于母体外周血样本的所有染色体的核苷酸序列信息分配至参考基因组的多个特异 性区域，统计样本在每个特异性区域内的NRSs值;利用GC含量修正样本在每个特异性区域 内的NRSs值，记为NRSs '值;基于NRSs '值，分别统计目标染色体和对照染色体上所有特异性 区域的NRSs'值的均值，分别对应记为第一均值和第二均值;将第一均值与第二均值进行差 异性检验，根据差异检验结果确定染色体是否存在非整倍性。
[0009] 进一步地，利用GC含量修正样本在每个特异性区域内的NRS值的步骤包括:利用修 正公式NRSs ' =NRSs Xa修正样本在每个特异性区域内的NRSs值，其中，a=NRSs/NRSs"， _为所有特异性区域NRSs值的中位数值，NRSs"为利用样本的每个特异性区域的GC含量 与NRSs值进行多项式样条拟合后获得的拟合值。
[0010]进一步地，在利用样本的每个特异性区域的GC含量与NRSs值进行多项式样条拟合 之前，处理方法还包括从样本的所有特异性区域中去除NRSs值异常的特异性区域的步骤， 优选采用线性拟合或者局部多项式回归拟合的方法去除NRSs值异常的特异性区域。
[0011] 进一步地，NRSc值为10000~50000中的任意整数。
[0012] 进一步地，目标染色体选自以下任意一条或几条的组合：13号染色体、18号染色 体、21号染色体、X染色体以及Y染色体;对照染色体选自以下任意一条或几条的组合:1号染 色体、2号染色体、3号染色体、4号染色体、5号染色体、6号染色体、7号染色体、8号染色体、9 号染色体、10号染色体、11号染色体以及12号染色体;优选地，对照染色体选自以下任意一 条或几条的组合：1号染色体、2号染色体、3号染色体、6号染色体、7号染色体、11号染色体、 12号染色体以及16号染色体。
[0013] 为了实现上述目的，根据本发明的另一个方面，提供了一种测序数据的处理装置， 该处理装置包括:测序模块，用于通过高通量测序获取来源于母体外周血样本的所有染色 体的核苷酸序列信息;特异性区域划分模块，用于根据NRSc值相等的原则将参考基因组划 分成多个特异性区域;分配模块，用于根据与参考基因组进行序列比对的原则，将来源于母 体外周血样本的所有染色体的核苷酸序列信息分配至参考基因组的多个特异性区域;第一 统计模块，用于统计样本在每个特异性区域内的NRSs值;修正模块，用于利用GC含量修正样 本在每个特异性区域内的NRSs值，记为NRSs '值;第二统计模块，用于基于NRSs '值，分别统 计目标染色体和对照染色体上所有特异性区域的NRSs '值的均值，记为第一均值和第二均 值;检验模块，用于将第一均值与第二均值进行差异性检验;确定模块，用于根据差异检验 结果确定染色体是否存在非整倍性。
[0014] 进一步地，修正模块包括:第一计算单元，用于计算所有特异性区域NRSs值的中位 数值拟合单元，用于利用样本的每个特异性区域的GC含量与NRSs值进行多项式样 条拟合，得到拟合曲线;获取单元，用于根据拟合曲线获取每个特异性区域的拟合值NRSs"； 第二计算单元，用于根据 (r ?^/NRSs"公式计算修正系数α;修正单元，用于根据修正公式 NRSs ' = NRSs Xa修正样本在每个特异性区域内的NRSs值。
[0015] 进一步地，拟合单元在执行利用样本的每个特异性区域的GC含量与NRSs值进行多 项式样条拟合，得到拟合曲线的步骤之前，拟合单元还包括过滤子单元，过滤子单元用于执 行从样本的所有特异性区域中去除NRSs值异常的特异性区域的步骤，优选过滤子单元为线 性拟合子单元或者局部多项式回归拟合子单元。
[0016] 进一步地，NRSc值为10000~50000中的任意整数。
[0017] 进一步地，目标染色体选自以下任意一条或几条的组合：13号染色体、18号染色 体、21号染色体、X染色体以及Y染色体;对照染色体选自以下任意一条或几条的组合:1号染 色体、2号染色体、3号染色体、4号染色体、5号染色体、6号染色体、7号染色体、8号染色体、9 号染色体、10号染色体、11号染色体以及12号染色体;优选地，对照染色体选自以下任意一 条或几条的组合：1号染色体、2号染色体、3号染色体、6号染色体、7号染色体、11号染色体、 12号染色体以及16号染色体。
[0018] 应用本发明的技术方案，通过以测序数据为基础，通过以相等条数的非重复序列 为原则来划分特异性区域，避免了各个特异性区域内非重复性序列数目不均一而导致的数 据波动，进而优化了染色体间核酸数据参数的相关性，利用与生物样本中临床相关染色体 的参数与其他非临床相关染色体区的参数比对，从而确定待测样本中染色体非整倍性是否 存在。该方法实现了单样本检测，可以不需要标准的正常样本，消除了对实验条件的依赖 性，加快了分析速度，是种简易、快速、准确的检测手段，其常染色体检测的准确率在99%以 上，假阳性率小于1%。
【附图说明】
[0019] 构成本申请的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示 意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
[0020] 图1示出了根据本发明的一种优选的实施例1中的SOOl样品（阴性样本)测序序列 中非重复序列在基因组上各特异性区域内的分布示意图；以及
[0021] 图2示出了图1中SOOl样品测序序列中过滤异常值后的非重复序列在基因组上各 特异性区域内的分布示意图；
[0022] 图3示出了图2中SOOl样品测序序列中过滤异常值后的非重复序列在基因组上各 特异性区域内的样条曲线拟合图；
[0023]图4a和图4b分别示出了实施例1中SOOl样品的各条常染色体的修正前和修正后的 特异性区域内的非重复序列的数目；其中，图4a显示修正前，图4b显示修正后；
[0024]图5a和图5b分别示出了另一种优选的实施例中S002样品的各条常染色体的修正 前和修正后的特异性区域内的非重复序列的数目；其中，图5a显示修正前，图5b显示修正 后；
[0025]图6a和图6b分别示出了又一优选的实施例中S007样品的各条常染色体的修正前 和修正后的特异性区域内的非重复序列的数目；其中，图6a显示修正前，图6b显示修正后； [0026]图7a和图7b分别示出了再一种优选的实施例中S006样品的各条常染色体的修正 前和修正后的特异性区域内的非重复序列的数目；其中，图7a显示修正前，图7b显示修正 后；
[0027]图8a、图8b和图8〇分别示出了本申请的实施例2中384例在线数据样本中第13号染 色体、第18号染色体和第21号染色体的Z值分布图，其中，图8a显示第13号染色体，图8b显示 第18号染色体，图8c显示第13号染色体。
【具体实施方式】
[0028]需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相 互组合。下面将参考附图并结合实施例来详细说明本发明。
[0029]术语解释：
[0030] 测序数据:指待测样本经过高通量测序获得的核苷酸序列信息。
[0031] kmer:将序列以逐个移动碱基的方式进行连续切割，得到的序列长度为k的核苷酸 序列，比如以下这条序列:ATCGTTGCTTAATGACGTCAGTCGAAT，如果是13-mer分析的话，k-mer 为ATCGTTGCTTAAT、TCGTTGCTTAATG、CGTTGCTTAATGA、GTTGCTTAATGAC ……。
[0032] 非重复性序列（non-repeated sequence，简称NRS);通过将待测样本测序获得的 序列与正常人类基因组进行比对，获得的全基因组水平上的唯一的kmer即为非重复性序 列。本申请中，按照等条数的非重复序列来划分特异性区域时，划分的条数是按照参考基因 组序列来进行划分的，因而，划分得到的每个特异性区域中非重复序列的条数记为NRSc，而 待测样本的测序序列实际在上述各特异性区域内的非重复序列的条数记为NRSs。
[0033]特异性区域（specified region，简称SR)，根据本发明中所描述的特异性区域的 划分方法所得到的基因组各条染色体上的特定区域。
[0034]染色体:既可以指整个染色体，也可以指染色体的一部分。处理一条染色体片段的 数学推导和处理所有染色体片段的数学推导是一致的，本领域的技术人员知悉相应的改变 方法。对照染色体是健康个体中的染色体或推定正常的染色体，包括统计学推定正常，这里 的染色体是单个染色体或染色体组(大于等于2条染色体，或者说是非13，18,214，￥的染色 体或其任意组合）。
[0035] "非整倍性"和"多倍性"是细胞具有的染色体数与通常的单倍体数η或双倍体数2n 不同的情况。非整倍体细胞可为具有三倍体的细胞，即具有一个染色体的三个拷贝数的细 胞;或为单倍体，即具有一个染色体的单拷贝的细胞。染色体非整倍性，改变了对应染色体 的表达量，可以通过新一代测序平台(NGS)结合生物信息学分析方法，根据测序比对结果统 计各条染色体的表达量可以判定待测样本是否存在该染色体的非整倍性变异。
[0036] 样本为细胞、组织或体液，可选自：母体全血(外周血）、血浆、血清、尿液、唾液、生 殖道冲洗液;胎儿细胞或胎儿细胞残留物、胚胎植入前的活检材料;羊水、绒膜绒毛样品等。 样品可来自任何动物，优选哺乳动物，更优选人。
[0037] 对DNA测序文库进行测序可以是双端短序列测序、单端长序列测序或单端短序列 测序。其中双端短序列是指紧接着5'端链接引物的小于50bp的序列和紧接着3'端链接引物 的小于50bp的序列。优选地，双端短序列是指紧接着5'端链接引物的不大于36bp的序列和 紧接着3 '端链接引物的不大于36bp的序列。
[0038]单端短序列是指紧接着5'端链接引物的小于50bp的序列或紧接着3'端链接引物 的小于50bp的序列。优选地，单端短序列是指紧接着5'端链接引物的不大于36bp的序列或 紧接着3'端链接引物的不大于36bp的序列。单端长序列是指紧接着5'端链接引物的大于 99bp的序列或紧接着3'端链接引物的大于99bp的序列。双端测序是指分别测试位于序列两 端的序列。单端测序是指对位于序列一端的序列进彳T测序。
[0039] 由于现有的染色体非整倍性的检测方法在准确性和便利性方面仍存在欠缺，为了 改善这一状况，在本申请一种典型的实施方式中，提供了一种测序数据的处理方法，该处理 方法包括:通过高通量测序获取来源于母体外周血样本的所有染色体的核苷酸序列信息； 将参考基因组划分成多个特异性区域，每个特异性区域内非重复序列的数目（记为NRSc)相 等;将来源于母体外周血样本的所有染色体的核苷酸序列信息分配至参考基因组的多个特 异性区域，统计样本在每个特异性区域内的NRSs值;利用GC含量修正样本在每个特异性区 域内的NRSs值，记为NRSs '值;基于NRSs '值，分别统计目标染色体和对照染色体上所有特异 性区域的NRSs'值的均值，分别对应记为第一均值和第二均值;将第一均值与第二均值进行 差异性检验，根据差异检验结果确定染色体是否存在非整倍性。
[0040] 本申请的上述处理方法，通过以测序数据为基础，通过以相等条数的非重复序列 为原则来划分特异性区域，避免了各个特异性区域内非重复性序列数目不均一而导致的数 据波动，进而优化了染色体间核酸数据参数的相关性，利用与生物样本中临床相关染色体 的参数与其他非临床相关染色体区的参数比对，从而确定待测样本中染色体非整倍性是否 存在。该方法实现了单样本检测，可以不需要标准的正常样本，消除了对实验条件的依赖 性，加快了分析速度，是种简易、快速、准确的检测手段，其常染色体检测的准确率在99%以 上，假阳性率小于1%。
[0041 ]具体地，上述差异性检验的方法可以为现有的各种差异性检验，比如，Z检验（Z-test)，u检验或t检验等。本申请优选Z检验。
[0042]上述处理方法中，利用GC含量修正样本在每个特异性区域内的NRS值的步骤可以 采用现有的GC校正方法也能提高检测的准确性。为了使检测准确性更高，在本申请一种优 选的实施例中，上述修正方法包括:利用修正公式NRSs ' =NRSs X α修正样本在每个特异性 区域内的NRSs值，其中，护.^gi/NRSs",?β为所有特异性区域NRSs值的中位数值，NRSs" 为利用样本的每个特异性区域的GC含量与NRSs值进行多项式样条拟合后获得的拟合值。修 正后的NRSs '更服从正态分布，因而使得后续差异性检验结果更准确。
[0043]拟合是根据已知离散的点(GC含量与NRSs值为X、Y轴的坐标）出32，'"，作}，通过 调整拟合函数中若干待定系数f(h，X2，···，λη)，使得该函数与已知点集的差别(最小二乘意 义）最小。已知点（X i，Y i ) ; X I〈 X 2〈…〈X η，i e Z是一系列观测值，符合某关系式 η
[(K 构建拟合函数身，使得:Yi=y(Xi)最小。如果拟合函数是非线 /=1 '1 性函数，则称之为非线性拟合，也叫作样条拟合。相应的，如果拟合函数为多项式，则可以称 之为多项式样条拟合。本发明优选多项式样条拟合，样条曲线是光滑的三次曲线。
[0044]三次样条曲线给定η个数据点，共有n-1个区间，每个区间的方程是:fi = ai+bi(x- 1〇+(^(11)2+(1办1)3，需确定4(11-1)个未知系数，通过连续性、节点处一阶导数相等、二 阶导数相等，可得到4n-6个方程，再人为添加2个边界条件。通过R软件系统的函数 smooth.spline完成样条拟合（http://www .stat.wisc.edu ~xie/smooth_spl ine tutorial.html)。
[0045]在利用样本的每个特异性区域的GC含量与NRSs值进行多项式样条拟合之前，上述 处理方法还包括从样本的所有特异性区域中去除NRSs值异常的特异性区域的步骤，可以采 用GC线性拟合的方法或者通过人工筛选的方式去除异常值，比如删除GC值为0、非重复序列 数目为〇或者非重复序列数目显著过多的窗口。在本申请中，优选采用局部多项式回归拟 合的方法去除NRSs值异常的特异性区域，该方法有利于排出部分非特异性区域因染色体结 构特异性而导致的内部非重复序列的数目过高或过低的异常特异性区域。此外，也可以采 用线性拟合拟合方法。拟合方法为统计学或生物信息学领域常用的去除异常值的方法，具 体方法在此不再赘述。
[0046]上述处理方法中，划分特异性区域时是按照NRSc值相等的原则来进行划分的，具 体NRSc值可以根据待测样本的基因组大小、序列复杂度等方式来进行确定。优选NRSc值为 10000~50000中的任意整数。
[0047]上述处理方法中，目标染色体和对照染色体可以根据待测样本的组织、细胞来源 不同、或者物种的不同或者实际检测需求的不同进行合理选择。当待测样本为人类时，优选 目标染色体选自以下任意一条或几条的组合:13号染色体、18号染色体、21号染色体、X染色 体以及Y染色体;对照染色体选自以下任意一条或几条的组合：1号染色体、2号染色体、3号 染色体、4号染色体、5号染色体、6号染色体、7号染色体、8号染色体、9号染色体、10号染色 体、11号染色体以及12号染色体;更优选地，对照染色体选自以下任意一条或几条的组合：1 号染色体、2号染色体、3号染色体、6号染色体、7号染色体、11号染色体、12号染色体以及16 号染色体。
[0048]在本申请另一种典型的实施方式中，还提供了一种测序数据的处理装置，该处理 装置包括:测序模块，用于通过高通量测序获取来源于母体外周血样本的所有染色体的核 苷酸序列信息;特异性区域划分模块，用于根据NRSc值相等的原则将参考基因组划分成多 个特异性区域;分配模块，用于根据与参考基因组进行序列比对的原则，将来源于母体外周 血样本的所有染色体的核苷酸序列信息分配至参考基因组的多个特异性区域;第一统计模 块，用于统计样本在每个特异性区域内的NRSs值;修正模块，用于利用GC含量修正样本在每 个特异性区域内的NRSs值，记为NRSs '值;第二统计模块，用于基于NRSs '值，分别统计目标 染色体和对照染色体上所有特异性区域的NRSs'值的均值，记为第一均值和第二均值;检验 模块，用于将第一均值与第二均值进行差异性检验;确定模块，用于根据差异检验结果确定 染色体是否存在非整倍性。
[0049]上述检测装置通过以测序模块获取的测序数据为基础，运用改进的特异性区域划 分模块以相等条数的非重复序列为原则来划分特异性区域，优化了染色体间核酸数据参数 的相关性，然后通过依次执行分配模块、第一统计模块、修正模块、第二统计模块、检验模 块，利用与生物样本中临床相关染色体的参数与其他非临床相关染色体区的参数比对，确 定模块最终通过检验模块的差异性检验结果确定待测样本中染色体非整倍性是否存在。该 装置实现了单样本测序数据的检测，且不需要标准的正常样本，消除了对实验条件的依赖 性，使分析速度加快、总体改进了对染色体异常的评估。是一种简易、快速、准确的染色体非 整倍性检测装置，其常染色体检测的准确率在99%以上，假阳性率小于1%。
[0050]具体地，上述检验模块可以为现有的各种差异性检验模块，比如，Z检验(Z- test) 模块，u检验模块或者t检验模块等。本申请优选Z检验模块。
[0051 ]上述修正模块采用现有的GC修正模块也能提高检测的准确性。在本申请一种优选 的实施例中，上述修正模块包括:第一计算单元，用于计算所有特异性区域NRSs值的中位数 值拟合单元，用于利用样本的每个特异性区域的GC含量与NRSs值进行多项式样条 拟合，得到拟合曲线;获取单元，用于根据拟合曲线获取每个特异性区域的拟合值NRSs"；第 二计算单元，用于根据α-?^/NRSs"公式计算修正系数α ;修正单元，用于根据修正公式 NRSs ' = NRSs Xa修正样本在每个特异性区域内的NRSs值。
[0052]上述优选的实施例中，通过利用多项式样条拟合的拟合单元具有拟合准确度高的 优点，以便于更准确地获取拟合值，相应地，通过第二计算单元计算得到的修正系数也更准 确，进而能够通过修正单元更准确地获得待测样本在每个特异性区域内的NRSs值，即得到 准确度更高的NRSs '值。
[0053]上述处理装置中，拟合单元在执行利用样本的每个特异性区域的GC含量与NRSs值 进行多项式样条拟合，得到拟合曲线的步骤之前，拟合单元还包括过滤子单元，过滤子单元 用于执行从样本的所有特异性区域中去除NRSs值异常的特异性区域的步骤，这样能够进一 步提高拟合单元在进行多项式样条拟合过程中拟合准确度。优选过滤子单元采用常规的线 性拟合子单元或者局部多项式回归拟合子单元进行异常值过滤即可。
[0054] 优选地，上述处理装置中NRSc值为10000~50000中的任意整数。
[0055]上述处理装置中，目标染色体和对照染色体可以根据待测样本的组织、细胞来源 不同、者物种的不同或者实际检测需求的不同进行合理选择。当待测样本为人类时，优选目 标染色体选自以下任意一条或几条的组合：13号染色体、18号染色体、21号染色体、X染色体 以及Y染色体;对照染色体选自以下任意一条或几条的组合：1号染色体、2号染色体、3号染 色体、4号染色体、5号染色体、6号染色体、7号染色体、8号染色体、9号染色体、10号染色体、 11号染色体以及12号染色体;更优选地，对照染色体选自以下任意一条或几条的组合：1号 染色体、2号染色体、3号染色体、6号染色体、7号染色体、11号染色体、12号染色体以及16号 染色体。
[0056]本申请的上述方法及其装置可与其它已知方法、装置或组合物联用，优选能改善 染色体异常检测技术的方法、装置或组合物。例如，母体生化指标的数学模型分析。
[0057]本申请所提供的上述方法，其具有高通量、低成本、简便、准确度与灵敏度高的优 势。现有的方法需要将待测样本与多个样本或标准正常样本进行比对分析，耗时长，且对样 本需求量大。本申请实现了单样本检测，可以不依赖于标准的正常样本，避免了对实验条件 的依赖性，加快了分析速度且提高了检测准确率。
[0058]本申请提供的上述方案是将DNA测序手段与生物信息分析的方法相结合，通过Z值 检验等差异性检验方法来判断染色体是否存在异常。如Z值在4.5之外，可确定为存在染色 体非整倍性。染色体异常优选是21号染色体三体，13号染色体三体，18号染色体，X染色体和 Y染色体的异常。
[0059]本申请方法尤其适用于检测染色体数量异常，优选染色体非整倍性数量异常，更 优选常染色体非整倍性异常。
[0060] 下面将结合具体的实施例来进一步说明本申请的有益效果。
[0061] 实施例1对待测样本测序数据的处理方法
[0062](一)对待测样本孕妇外周血中游离DNA片段进行高通量测序 [0063] (1)采集孕妇全血，经过预处理得到血浆；
[0064] 在获准同意通知书后，通过静脉穿刺从怀孕22周妇女（即后续表2中的样品S001) 取血采血量5-10ml，加入到乙二胺四乙酸(EDTA)管中，血液样品经高速离心后得到除去了 血细胞的血浆样品，每个样品血浆量约为700u 1。
[0065] (2)提取血浆 DNA;
[0066] 使用Magen公司生产的DNA抽提试剂盒HiPure Circulating DNA Kits来提取血衆 中的DNA(产品号为D3180-02)。
[0067] (3)将从血浆抽提得到的DNA制备成可供新一代高通量测序平台测序的文库
[0068] 血浆DNA使用T4 DNA聚合酶、T4 PNK和Klenow酶进行末端修复和加A处理，使用T4 DNA连接酶和测序接头进行加接头处理。最后使用加有标签的文库引物进行PCR，并使用磁 珠进行纯化筛选，最终得到上机用的测序文库。
[0069] (4)对制备好的文库进行DNA测序
[0070] 测序文库在Illumina的cBot仪器上扩增，DNA单端测序文库制成DNA簇，得到海量 测序读长为36bp的序列。
[0071](二)确定血浆中的DNA片段的序列信息
[0072] 1.对正常人类参考基因组进行特异性区域划分与统计
[0073] (1)筛选非重复性序列
[0074]将人类参考基因组（hgl9 GRCh37http : //www · ncbi · nlm · nih · gov/pro jects/ genome/assembly/grc/)，切分成长度为35bp，偏移量为Ibp的海量kmer集合;从中筛选得到 全基因组上唯一的km er，即非重复性序列，并记录对应的位置坐标信息。
[0075] (2)特异性区域划分
[0076]从1号染色体的第一条非重复性序列开始记录起始位置，直至累积至20000条时记 录其终止位置，将此定义为1号染色体上的第一个特异性区域，每个特异性区域之间不存在 重叠。
[0077]针对1号染色体直至Y染色体均重复上边的处理步骤，获得所有染色体特异性区域 的位置信息和GC含量(对正常人类参考基因组进行特异性区域划分只需进行一次，后续每 个待测样品按照参考基因组划分的特异性区域进行处理即可）。
[0078] (3)特异性区域统计
[0079]统计各条染色体上的特异性区域数量及区域内所有非重复性序列的GC含量分布 情况。
[0080] 2.样品DNA序列比对
[0081 ] 通过生物信息序列比对软件BWA(Burrows_Wheeler Aligner)，将测序所得的DNA 序列与正常人类参考基因组(hgl9，GRCh37)进行不容错比对(完全匹配，不允许有碱基错 配），确定所有测序DNA序列在基因组上的详细位置信息，包括染色体来源、染色体上的坐标 以及在基因组特异性区域分布情况等(表2中SOOl样品测序序列中非重复序列在基因组上 各特异性区域内的分布情况见图I)。
[0082](三)确定待测染色体的表达量 [0083] 1、过滤异常值
[0084]将待测样本的基因组特异性区域的GC含量以及该区域内非重复性序列的数目 (NRSs)通过loess函数进行局部多项式回归拟合(线性拟合亦可），将NRSs数在拟合值正负3 倍标准差之外(P〈〇.005)的定义为异常值，将异常值过滤之后的分布如图2所示。
[0085] 2、加权修正
[0086] 按照GC含量将待测样本的基因组的所有特异性区域进行分类后，进行样条曲线拟 合得到每个GC含量对应的NRSs的拟合值，记为NRSs"，其对应的分布情况如图3所示。
[0087] 其中，具体拟合步骤为：以NRSs的中位数值为基线，将NRS的拟合值NRSs"与 基线值相比得到修正系数α，计算公式如下；
[0088] ?= NRSs/NRSs" (1)
[0089] NRSs，=NRSsXa (2)
[0090] 上述公式是针对待测样本基因组上的每个特异性区域进行计算的，其中，_指 的是基因组上所有特异性区域上NRS数的中位数值，NRSs"为拟合值，NRSs '是修正后的非重 复序列数。
[0091] 从以下图4a与图4b、图5a与图5b、图6a与图6b以及图7a与图7b的修正前和修正后 可以看出，未经修正的数据波动比较大，直接进行染色体之间的差异性比较容易导致假阴 性或假阳性的检测结果。而修正后各染色体的特异性区域内的非重复序列数目分布情况趋 于平稳，数据差异性更显著，更容易判断异常值，表明本申请的方法能消除GC结构差异，并 避免GC偏好性问题。可用于染色体非整倍性异常的检测，减少假阴性结果的出现，如下图7a 和图7b的chr21对应的NRS数同其它常染色体比较明显高出，对应的检测结果是该样品21号 染色体存在非整倍性异常的风险极高。
[0092](四）Z值检验判断染色体表达量是否存在显著性差异
[0093]以NRSs经GC修正后的NRSs'，将目标检测染色体(chr21、chrl8、chrl3、XSYW^；f 有特异性区域的NRSs '的均值，与对照染色体组合(chrl、chr2……chrl2)的所有NRSs '的均 值进行差异性比较，得到检测值Z(Z-Score)，根据Z值判断目前目标染色体是否存在非整倍 性变异。当 Z-score 2 4.5或Z-score < -4.5，即检测结果为三体变异高风险，或单体变异高 风险；当-4.5<Z-sc〇re<4.5，即检测结果为非整倍性变异低风险。
[0094] 或通过管家基因的分布情况，筛选出对照染色体组合，包括chrl、chr2、chr3、 chr6、chr7、chr11、chr12、chr16 〇
[0095]实施例2有效性评价
[0096](一)利用在线数据样本进行评价
[0097] 实施例1中所示的处理方法中的步骤可以模块或单元的形式通过计算装置来实 现。为了评价实施例1的方法的有效性，以下用能够执行上述步骤的模块或单元形成的处理 装置来进行测试。该处理装置包括：
[0098] 测序模块，用于通过高通量测序获取来源于母体外周血样本的所有染色体的核苷 酸序列信息；
[0099] 可选地，上述模块包括 Il Iumina的 cBot 仪器、I Ilumina 的Genome Analyzer、 HiSeq2000/2500、Hiseq3000/4000、NextseqCN500等配套型号测序仪或者 Life Technologies公司的SOLiD等配套测序仪中执行测序功能的模块。
[0100]特异性区域划分模块，调用特异性区域划分模块程序，根据NRSc值相等的原则将 参考基因组划分成多个特异性区域;可以按照10000~50000(优选20000)中任一整数条的 非重复序列为单位进行划分，以克服按长度如20Kb或50Kb划分的特异性区域内存在非重复 序列数目差异较大而数据均一性差的缺陷。分配模块，运行分配模块，将测序模块输出的结 果与参考基因组进行序列比对，将来源于母体外周血样本的所有染色体的核苷酸序列信息 分配至特异性区域划分模块产生的特异性区域内；
[0101] 可选地，能够执行序列比对原则的模块如BWA模块、BOWTIE模块或者NOVOALIGr^t 块用来进行待测样本测序数据的分配，
[0102]第一统计模块，用于统计样本在每个特异性区域内的NRSs值;可选地的统计模块 有SAMT00LS模块；
[0103]修正模块，用于利用GC含量修正样本在每个特异性区域内的NRSs值，记为NRSs' 值；
[0104] 优选地，修正模块包括:第一计算单元，用于计算所有特异性区域NRSs值的中位数 值拟合单元，用于利用样本的每个特异性区域的GC含量与NRSs值进行多项式样条 拟合，得到拟合曲线;获取单元，用于根据拟合曲线获取每个特异性区域的拟合值NRSs"；第 二计算单元，用于根据 a=i^/NRSS"公式计算修正系数α ;修正单元，用于根据修正公式 NRSs ' = NRSs Xa修正样本在每个特异性区域内的NRSs值。
[0105] 第二统计模块，用于基于NRSs'值，分别统计目标染色体和对照染色体上所有特异 性区域的NRSs'值的均值，记为第一均值和第二均值；
[0106] 检验模块，用于将第一均值与第二均值进行差异性检验;可选地，采用Z检验模块 来进行差异性分析；
[0107] 确定模块，用于根据差异检验结果确定染色体是否存在非整倍性；
[0108] 优选地，当目标染色体为常染色体，而-4.5 < Z值< 4.5时，用于确定目标染色体不 存在非整倍性，否则，确定存在非整倍性。
[0109] 以源自不同实验室、不同NGS平台的数据（从NCBI的SRA数据库http:// WWW. ncbi . nlm. nih. gov/sra/中下载的其他机构上传的无创产前基因检测项目临床研究孕 妇外周血的高通量测序数据，其中包含384例样品数据)为样本来进一步说明本申请处理装 置的有效性及通用性。
[0110]其中，针对该384例样品中的第21号、18号和13号染色体的检测结果如下表1所示： [0111]表1.384例NCBI在线数据阳性样本检出结果。
LU114J 附：上表l中，··chr"表不染β体；··gc"表不GC?、重；··ZV"表不ZValue，Z值 ；"TEST" 表示通过该方法得到的染色体非整倍性异常检测结果。
[0115] 从上述表1中及图8a、图8b和图8c可知，检出1例T13阳性样品SRR358477,其余样品 的13号染色体的Z值均稳定分布在（-4.5，4.5)区间内；检出5例T18阳性样品SRR357943、 51^357972、51^358089、51^358257、51^358325，其余样品18号染色体的2值均稳定分布在(-4.5,4.5)区间内；检出7例了21阳性样品51^357843、51^358020、51^358126、51^358144、 31^358322、51^358352、51^358353，其余样品21号染色体的2值均稳定分布在(-4.5,4.5)区 间内。
[0116] (二)血液样品进行评价
[0117] 以68例样品（由卫生部临检中心和北京人民医院提供）的检测结果为例来评价上 述处理方法，评价结果见表2,表2中仅显示前30个样本的结果，并用核型结果进行验证。
[0118] 表2.血液样品检测结果
[0120] 注：NRSs -(Η3,18:，21)代表的是目标检测染色体未经修正前，所有特异性区域中 NRS的平均值。HRSs'chri(：i=13，18,21)代表的是经过特异性区域GC含量加权修正后，该染色体 上所有特异性区域NRS的平均值。Ζν^( i = 13，18，21)代表的是该染色体通过与对照染色 体进行显著性差异分析得到的Z值。TEST代表的是通过该方法得到的染色体非整倍性异常 检测结果，N(Negative)表示检测结果为阴性，未检出明显异常。T13/T18/T21表示检测结果 显示目标检测染色体存在非整倍性异常。核型(Karyotype)代表的是临床核型分析结果，即 金标准结果(46，XN代表的是核型正常样品的染色体数目以及性染色体情况，47，XN，+21代 表的是该样品核型分析显示有47条染色体，比正常核型多出一条21号染色体，即唐氏综合 征)。
[0121] 表2数据表明，根据显著性差异检验结果：S0002和S0013两个样品ZVchri3均大于等 于4.5，判定13号染色体非整倍性异常存在高风险;S0007和S0012两个样品ZV ehrl8均大于等 于4.5，判定18号染色体非整倍性异常存在高风险，S0003、S0006和SOOl 1这三个样品ZVchr2I 都大于等于4.5,判定21号染色体非整倍性异常存在高风险。针对21号、18号和13号染色体， 本申请的检测结果均与染色体核型分析结果一致，本方法检测结果判定为低风险的样品， 即ZV值在-4.5到4.5之间的样品，其核型分析结果也均为正常。表明本方法在用于染色体非 整倍性异常的检测时检测准确性较高。
[0122] 实施例3稳定性与数据量研究
[0123] (一)样品稳定性
[0124] 利用上述方法，对8002、8006、8007、8008这四个样品（对应的核型结果分别为113 阳性、T21阳性、T18阳性和正常)分别重复测8次，统计染色体相对表达量(记为CR)及Z检验 值(记为ZV)研究数据波动情况，以评价该检测方法的稳定性，评价结果见表3。
[0125]
[0126] 表3 · s002、s006、s007、s008重复性检测数据总表

L〇129」上表3中，Mean代表平均值，SD代表标准差，CV代表离散系数。从表3可知，该4个样 品重复检测8次对应的CR值的CV(离散值）均小于O . 01，且ZV的波动（SD值)也均在± 1.1之 内，数据波动较小，表明本方法的稳定性较好。
[0130] (二)数据量研究
[0131] 在测序数据量在0.25M(raw reads)至15M范围，研究基因组特异性区域上NRS数目 的波动情况。针对s002、s006、s007、s008这四个样品的测序数据(对应的核型结果分别为 T13阳性、T21阳性、T18阳性和正常），随机截取2M至15M的数据量，以此进行基因组比对及 统计ZV与CV(该样品在基因组所有特异性区域内非重复序列数的离散系数值。统计结果见 表4。
[0132] 表4.不同测序数据量对应的CV值(离散系数)和ZV值(Z值）

L〇135」从上表4?知，本万法迠亘厂谮数据量的染色体检测，X其在数据量为IM及IM以上 时，数据的稳定性与Z检验的结果都较好。
[0136] 从以上的描述中，可以看出，本发明上述的实施例实现了如下技术效果:通过以测 序数据为基础，通过以相等条数的非重复序列为原则来划分特异性区域，避免了各个特异 性区域内非重复性序列数目不均一而导致的数据波动，进而优化了染色体间核酸数据参数 的相关性，利用与生物样本中临床相关染色体的参数与其他非临床相关染色体区的参数比 对，从而确定待测样本中染色体非整倍性是否存在。该方法实现了单样本检测，且不需要标 准的正常样本，消除了对实验条件的依赖性，使分析速度加快、总体改进了对染色体异常的 评估。提供了一种简易、快速、准确的染色体非整倍性检测手段，其常染色体检测的准确率 在99%以上，假阳性率小于1%。本方法相对多样本方法降低了假阴性检出率;相对现有的 单样本方法，对测序数据量的要求更小。
[0137] 显然，本领域的技术人员应该明白，上述的本申请的一些模块、元件或一些步骤可 以用通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算装 置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将 它们存储在存储装置中由计算装置来执行，或者将它们分别制作成各个集成电路模块，或 者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本申请不限制于任 何特定的硬件和软件结合。
[0138]以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技 术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修 改、等同替换、改进等，均应包含在本发明的保护范围之内。
【主权项】
1. 一种测序数据的处理方法，其特征在于，所述处理方法包括： 通过高通量测序获取来源于母体外周血样本的所有染色体的核巧酸序列信息； 将参考基因组划分成多个特异性区域，每个所述特异性区域内非重复序列的数目NRSc 值相等； 将所述来源于母体外周血样本的所有染色体的核巧酸序列信息分配至所述参考基因 组的多个所述特异性区域，统计所述样本在每个所述特异性区域内的NRSs值； 利用GC含量修正所述样本在每个所述特异性区域内的NRSs值，记为NRSs '值； 基于所述NRSs '值，分别统计目标染色体和对照染色体上所有特异性区域的NRSs '值的 均值，分别对应记为第一均值和第二均值； 将所述第一均值与第二均值进行差异性检验，根据差异检验结果确定染色体是否存在 非整倍性。2. 根据权利要求1所述的处理方法，其特征在于，利用GC含量修正所述样本在每个所述 特异性区域内的NRS值的步骤包括： 利用修正公式NRSs'=NRSsXa修正所述样本在每个所述特异性区域内的NRSs值，其 中，所述α= NRSs/ NRSs", NRSs为所有特异性区域NRSs值的中位数值，NRSs"为利用所述样 本的每个所述特异性区域的GC含量与NRSs值进行多项式样条拟合后获得的拟合值。3. 根据权利要求2的处理方法，其特征在于，在利用所述样本的每个所述特异性区域的 GC含量与NRSs值进行多项式样条拟合之前，所述处理方法还包括从所述样本的所有特异性 区域中去除NRSs值异常的特异性区域的步骤，优选采用线性拟合或者局部多项式回归拟合 的方法去除NRSs值异常的特异性区域。4. 根据权利要求1的处理方法，其特征在于，所述NRSc值为10000~50000中的任意整 数。5. 根据权利要求1至4中任一项所述的处理方法，其特征在于， 所述目标染色体选自W下任意一条或几条的组合：13号染色体、18号染色体、21号染色 体、X染色体W及Y染色体； 所述对照染色体选自W下任意一条或几条的组合：1号染色体、2号染色体、3号染色体、 4号染色体、5号染色体、6号染色体、7号染色体、8号染色体、9号染色体、10号染色体、11号染 色体W及12号染色体； 优选地，所述对照染色体选自W下任意一条或几条的组合：1号染色体、2号染色体、3号 染色体、6号染色体、7号染色体、11号染色体、12号染色体W及16号染色体。6. -种测序数据的处理装置，其特征在于，所述处理装置包括： 测序模块，用于通过高通量测序获取来源于母体外周血样本的所有染色体的核巧酸序 列信息； 特异性区域划分模块，用于根据NRSc值相等的原则将参考基因组划分成多个特异性区 域； 分配模块，用于根据与所述参考基因组进行序列比对的原则，将所述来源于母体外周 血样本的所有染色体的核巧酸序列信息分配至所述参考基因组的多个所述特异性区域； 第一统计模块，用于统计所述样本在每个所述特异性区域内的NRSs值； 修正模块，用于利用GC含量修正所述样本在每个所述特异性区域内的NRSs值，记为 NRSs'值； 第二统计模块，用于基于所述NRSs '值，分别统计目标染色体和对照染色体上所有特异 性区域的NRSs'值的均值，分别对应记为第一均值和第二均值； 检验模块，用于将所述第一均值与第二均值进行差异性检验； 确定模块，用于根据差异检验结果确定染色体是否存在非整倍性。7. 根据权利要求6所述的处理装置，其特征在于，所述修正模块包括： 第一计算单元，用于计算所有特异性区域NRSs值的中位数值而?京； 拟合单元，用于利用所述样本的每个所述特异性区域的GC含量与NRSs值进行多项式样 条拟合，得到拟合曲线； 获取单元，用于根据所述拟合曲线获取每个所述特异性区域的拟合值NRSs"; 第二计算单元，用于根据α= NRSs/ NRSs'，公式计算修正系数α; 修正单元，用于根据修正公式NRSs ' =NRSs X α修正所述样本在每个所述特异性区域内 的NRSs值。8. 根据权利要求7的处理装置，其特征在于，所述拟合单元在执行利用所述样本的每个 所述特异性区域的GC含量与NRSs值进行多项式样条拟合，得到拟合曲线的步骤之前，还包 括过滤子单元，所述过滤子单元用于执行从所述样本的所有特异性区域中去除NRSs值异常 的特异性区域的步骤，优选所述过滤子单元为线性拟合子单元或者局部多项式回归拟合子 单元。9. 根据权利要求6的处理装置，其特征在于，所述NRSc值为10000~50000中的任意整 数。10. 根据权利要求6至9中任一项所述的处理装置，其特征在于， 所述目标染色体选自W下任意一条或几条的组合：13号染色体、18号染色体、21号染色 体、X染色体W及Y染色体； 所述对照染色体选自W下任意一条或几条的组合：1号染色体、2号染色体、3号染色体、 4号染色体、5号染色体、6号染色体、7号染色体、8号染色体、9号染色体、10号染色体、11号染 色体W及12号染色体； 优选地，所述对照染色体选自W下任意一条或几条的组合：1号染色体、2号染色体、3号 染色体、6号染色体、7号染色体、11号染色体、12号染色体W及16号染色体。
【文档编号】G06F19/22GK105844116SQ201610161767
【公开日】2016年8月10日
【申请日】2016年3月18日
【发明人】张必良, 曹亮
【申请人】广州市锐博生物科技有限公司