CN104053789A

CN104053789A - 确定胎儿基因组中预定区域碱基信息的方法、***和计算机可读介质

Info

Publication number: CN104053789A
Application number: CN201280067404.1A
Authority: CN
Inventors: 陈盛培; 葛会娟; 李旭超; 易赏; 汪建; 王俊; 杨焕明; 张秀清
Original assignee: BGI Shenzhen Co Ltd
Current assignee: Shenzhen Huada Medical Laboratory
Priority date: 2012-05-14
Filing date: 2012-05-14
Publication date: 2014-09-17
Anticipated expiration: 2032-05-14
Also published as: JP6045686B2; PL2851431T3; KR20140146193A; RU2014150655A; AU2012380221A1; ES2656023T3; EP2851431B1; AU2012380221B2; KR101770460B1; RU2597981C2; US20150094210A1; WO2013170429A1; EP2851431A1; CN104053789B; HK1196401A1; US20180320235A1; JP2015525062A; EP2851431A4; SG11201407515RA

Abstract

提供了确定胎儿基因组中预定区域碱基信息的方法、***和计算机可读介质。其中，确定胎儿基因组中预定区域碱基信息的方法，包括下列步骤：针对胎儿基因组DNA样本，构建测序文库；对测序文库进行测序，以便获得胎儿的测序结果，该胎儿的测序结果由多个测序数据构成；基于胎儿的测序结果，结合胎儿遗传相关个体的遗传信息，根据隐马尔可夫模型，确定预定区域的碱基信息。

Description

确定胎儿基因组中预定区域碱基信息的方法、 ***和计算机可读介质优先权信息

无技术领域

本发明涉及确定胎儿基因组中预定区域碱基信息的方法、 ***和计算机可读介质。背景技术

遗传性疾病是由于遗传物质发生改变而造成的疾病，具有先天性、家族性、终身性和遗传性的特点。遗传性疾病可分为 3个大类：单基因遗传病、多基因遗传病及染色体异常。其中单基因病多由于单个致病基因的显性或隐性遗传所致基因功能异常；而多基因遗传病则是由多个基因变化影响所致的疾病，会在一定程度上受到外界环境因素的影响；染色体异常包括数目异常和结构异常，最为多见的是由于第 21号染色体三体所致的唐氏综合症，患儿表现为先天愚型和肢体形状异常等其他先天性特征。由于目前对遗传性疾病尚无有效的治疗方式，只能针对性地进行支持治疗或者药物緩解，费用昂贵，给社会和家庭带来沉重经济和精神负担。因此，在孩子出生前就对孩子的患病状态进行检测， ^故好预防工作，以达到优生优育的目的，是十分必要的。

然而，目前的相关检测手段仍有待改进。发明内容

本发明旨在至少解决现有技术中存在的技术问题之一。

在本发明的一个方面，本发明提出了一种确定胎儿基因组中预定区域碱基信息的方法。根据本发明的实施例，该方法包括下列步骤：针对胎儿基因组 DNA样本，构建测序文库；对所述测序文库进行测序，以便获得胎儿的测序结果，所述胎儿的测序结果由多个测序数据构成；基于所述胎儿的测序结果，结合胎儿遗传相关个体的遗传信息，根据隐马尔可夫模型，确定所述预定区域的碱基信息。子代的基因组形成过程，相当于亲代基因组的一次随机重组（即连锁互换单倍体型重组，以及配子的随机组合）。对于孕期血浆，假若我们将胎儿的单倍型（父母单倍型的重组型）作为隐含状态（ hidden states ), 可以将血浆的测序数据当做观察序列（ observations ),借助先验数据推算出状态转移概率（ transition probabilities ), 观察序列概率分布（observation symbol probabilities ) 和初始状态概率分布（initial state distribution ), 我们则可以通过诸如惠特比算法（ Viterbi algorithm )根据隐马尔可夫模型推断出最可能的胎儿单倍型组合，从而获得更多胎儿的信息。因而，根据本发明的实施例，借助隐马尔可夫模型，例如可以通过利用惠特比算法（ Viterbi algorithm ), 参考胎儿遗传相关个体的遗传信息，可以确定胎儿基因组中特定区域的核酸序列，由此，可以有效地对胎儿基因组的遗传信息进行产前检测。在本发明的又一方面，本发明提出了一种用于确定胎儿基因组中预定区域碱基信息的 ***。才艮据本发明的实施例，该***包括：文库构建装置，所述文库构建装置适于针对胎儿基因组 DNA样本，构建测序文库；测序装置，所述测序装置与所述文库构建装置相连，并且适于对所述测序文库进行测序，以便以便获得胎儿的测序结果，所述胎儿的测序结果由多个测序数据构成；分析装置，基于所述胎儿的测序结果，结合胎儿遗传相关个体的遗传信息，根据隐马尔可夫模型，确定所述预定区域的碱基信息。利用该***，能够有效地实施前面所述的确定胎儿基因组中预定区域碱基信息的方法，可以借助隐马尔可夫模型，例如可以通过利用惠特比算法（Viterbi algorithm ), 参考胎儿遗传相关个体的遗传信息，可以确定胎儿基因组中特定区域的核酸序列，由此，可以有效地对胎儿基因组的遗传信息进行产前检测，从而可以有效地对胎儿基因组的遗传信息进行产前确定。

在本发明的另一方面，本发明还提出了一种计算机可读介质。根据本发明的实施例，该计算机可读介质上存储有指令，所述指令适于被处理器执行以便基于胎儿的测序结果，结合胎儿遗传相关个体的遗传信息，根据隐马尔可夫模型，确定预定区域的碱基信息。利用本发明的计算机可读介质，能够有效地被处理器执行其存储的指令，以便借助隐马尔可夫模型，例如可以通过利用惠特比算法（Viterbi algorithm ), 基于胎儿的测序结果，参考胎儿遗传相关个体的遗传信息，可以确定胎儿基因组中特定区域的核酸序列，由此，可以有效地对胎儿基因组的遗传信息进行产前检测。

本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。附图说明

本发明的上述和 /或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图 1为根据本发明一个实施例的利用隐马尔可夫模型进行分析的流程示意图；以及图 2为根据本发明的一个实施例的用于确定胎儿基因组中预定区域核酸序列的***的结构示意图。发明详细描述

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。

需要说明的是，术语 "第一" 、 "第二" 仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有 "第一"、 "第二" 的特征可以明示或者隐含地包括一个或者更多个该特征。进一步地，在本发明的描述中，除非另有说明， "多个" 的含义是两个或两个以上。

确定胎儿基因组中预定区域信息的方法在本发明的第一方面，本发明提出了一种确定胎儿基因组中预定区域碱基信息的方法。根据本发明的实施例，该方法包括下列步骤：

首先，针对胎儿基因组 DNA样本，构建测序文库。根据本发明的实施例，胎儿基因组 DNA样本的来源并不受特别限制。根据本发明的一些实施例，可以釆用任何含有胎儿核酸的孕妇样本。例如，根据本发明的实施例，可以釆用的孕妇样本为孕妇乳汁、尿液和外周血。其中，优选孕妇外周血。釆用孕妇外周血作为胎儿基因组 DNA样本的来源，可以有效地实现无创取样方式获得胎儿基因组 DNA, 从而可以在不影响胎儿正常发育的前提下，对胎儿的基因组进行有效监测。关于针对核酸样本，构建测序文库的方法和流程，本领域技术人员可以根据不同的测序技术进行适当选择，关于流程的细节，可以参见测序仪器的厂商例如 Illumina公司所提供的规程，例如参见 Illumina公司 Multiplexing Sample Preparation Guide ( Part#l 005361; Feb 2010 )或 Paired-End SamplePrep Guide ( Part#l 005063； Feb 2010 ), 通过参照将其并入本文。根据本发明的实施例，从生物样本提取核酸样本的方法和设备，也不受特别限制，可以釆用商品化的核酸提取试剂盒进行。

在构建测序文库后，将测序文库应用于测序仪器，对测序文库进行测序，并获得相应的测序结果，该测序结果是由多个测序数据构成的。根据本发明的实施例，可以用于进行测序的方法和设备并不受特别限制，包括但不限于双脱氧链终止法；优选高通量的测序方法，由此，能够利用这些测序装置的高通量、深度测序的特点，进一步提高测序效率。从而，能够提高后续对测序数据进行分析，尤其是统计检验分析时的精确性和准确度。所述高通量的测序方法包括但不限于第二代测序技术或者是单分子测序技术。所述第二代测序平台 ( Metzker ML. Sequencing technologies-the next generation. Nat Rev Genet. 2010 Jan;ll(l):31-46 ) 包括但不限于 Illumina-Solexa ( GA^TM,HiSeq2000™等）、 ABI-Solid 和 Roche-454 (焦碑酸测序）测序平台；单分子测序平台（技术）包括但不限于 Helicos公司的真实单分子测序技术 ( True Single Molecule DNA sequencing ), Pacific Biosciences公司单分子实时测序（ single molecule real-time (SMRT™) ), 以及 Oxford Nanopore Technologies 公司的纳米孔测序技术等（Rusk, Nicole (2009-04-01). Cheap Third-Generation Sequencing. Nature Methods 6 (4): 244-245 )„ 随着测序技术的不断进化，本领域技术人员能够理解的是还可以釆用其他的测序方法和装置进行全基因组测序。根据本发明的具体示例，可以利用选自 Illumina-Solexa、 ABI-SOLiD、 Roche-454和单分子测序装置的至少一种对所述全基因组测序文库进行测序。

任选地，在得到测序结果之后，可以将所述测序结果与参照序列进行比对，以便确定与所述预定区域对应的测序数据。在本文中所使用的术语 "预定区域" 应作广义理解，是指任何包含可能发生预定事件位点的核酸分子的区域。对于 SNP分析而言，可以是指包含 SNP位点的区域。对于分析染色体非整倍性，则预定区域指的是所要分析的染色体的全长或者部分，即选择所有来自该染色体的测序数据。从测序结果中选择来自相应区域的测序数据的方法可以不受特别限制。根据本发明的实施例，可以通过将所得到的所有测序数据与已知的核酸参照序列进行比对，从而得到来自于预定区域的测序数据。另外，根据本发明的实施例，预定区域也可以是基因组上不连续的多个分散点。根据本发明的实施例，可以使用的参照序列的类型并不受特别限制，可以为任何含有感兴趣区域的已知序列。根据本发明的实施例，可以釆用已知的人类参考基因组作为参照序列。例如，根据本发明的实施例，釆用的人类参考基因组为 NCBI 36.3 , HG18。另外，根据本发明的实施例，进行比对的方法并不受特别限制。根据本发明的具体实施例，可以釆用 SOAP进行比对。

接下来，基于与预定区域对应的测序数据，确定预定区域中的部分核酸序列；以及基于所确定的预定区域中的部分核酸序列，按照惠特比算法，确定预定区域的其他核酸序列，以便获得预定区域的核酸序列。根据本发明的实施例，可以通过基于所述胎儿的测序结果，结合胎儿遗传相关个体的遗传信息，根据隐马尔可夫模型，确定所述预定区域的碱基信息。根据本发明的实施例，可以利用惠特比算法（Viterbi algorithm ), 借助隐马尔可夫模型，确定胎儿基因组中特定区域的碱基信息。由此，可以有效地对胎儿基因组的遗传信息进行产前检测。

下面参考图 1 , 对利用惠特比算法借助隐马尔可夫模型进行分析的原理进行详细描述：在本文中所使用的术语 "胎儿遗传相关个体" 指的是在遗传意义上，与胎儿之间具有亲缘关系的个体，例如根据本发明的实施例，可以釆用的 "胎儿遗传相关个体" 为胎儿的亲代例如父母。由此，子代的基因组形成过程，相当于亲代基因组的一次随机重组（即连锁互换单倍体型重组，以及配子的随机组合）。对于孕期血浆，假若我们将胎儿的单倍型（父母单倍型的重组型）作为隐含状态（hidden states ), 可以将血浆的测序数据当做观察序列 ( observations ), 借助先验数据推算出状态转移概率（transition probabilities )、观察序列概率分布 ( observation symbol probabilities )和初始状态概率分布 ( initial state distribution ), 我们则可以通过惠特比算法（ Viterbi algorithm )推断出最可能的胎儿单倍型组合，从而获得更多胎儿的信息。

详细分析步骤如下：

记号：

I.需要检测的位点数为 N。

II. 父母的单倍型分别记为 FH = [fh。,fh 和 MH =

其巾，

mK = {m_u,...,m ,...,m^ } , fh_k = { ".., ""， ^} , VK ≡{A,C,G,T} ^ k e {0, l] , i = 2, 3, ...,N。

III. 将未知的胎儿单倍型记为 H = { ¾， ¾ } , 特别地，。和/ ¾分别遗传自母亲和父亲。

= ^，…，"^,，…，） ' Κ = } , 其中， χ {0，1}， y_{ e{0,l}₅ 下标 ,.和 χ.组成的序列对， q_t = {χ }^组成了我们需要解码的隐藏状态，而所有可能出现的隐藏状态组成了集合 β。

IV. 测序数据记为， S= ..， }, 其中， = .,_A, 代表此位点的测序信息，包含了 ACGT四种碱基的数量。

V. 平均胎儿浓度和平均测序错误率分别记为 f和 ί？。第一步，构建初始状态概率分布向量，以及单倍体重组转移矩阵：

I. 初始状态概率分布记为； τ={ ( j eQ )₀ 根据本发明的实施例，在没有参考数据的情况下，可以设^ ？！^^ 全，即每种隐藏状态在第一个位点出现的可能性相等。

II. 根据本发明的实施例，记单倍体重组概率为 p_r = re/N , 其中 re代表人类配子基因组重组平均次数，为先验数据在 25到 30之间。

III. 根据本发明的实施例，记单倍体重组转移矩阵记为 A = ( j,k eQ ), 其中为隐藏

胎儿单倍型 h₀ = {m _i,...,m,. ...,m^}, = (/ ，…，^,，…，^^ }的下标和 y,组成的序列对， = {χ _;νί}组成了我们需要解码的隐藏状态。举个例子， = 0代表， "母源性染色体上, 对应基因座上等位基因型为 w¾。 "。

第二步，构建观察序列概率矩阵：

根据本发明的实施例，记观察序列概率矩阵为 B = {b_j (_Si)} ( i = 1,2,3,...,N , j eQ ), 其中 .(5_;)代表 "在位点 , 考虑母亲单倍型和胎儿单倍型（状态 j={x_h yi} )时，观测 il 这种测序信其中^ ^代表 "在位点考虑母亲单倍型和胎儿单倍型（状态 j={Xi, yi} ) 时，该碱基出现的可能

其中指示函

此步是进行 HMM参数，每个位点的观察序列概率分布 .(5_;)计算，即计算每个位点上不同胎儿单倍型（隐藏状态）下，血浆出现当前测序数据（观察序列）的可能性。

第三步，构建局部概率矩阵，和逆向指针 (下面以一维局部概率矩阵构建为例；)：

定义局部概率 · (q_t ) = (max^. (q · α_{¾ ι¾} ) · _b^ (_5; )

定义逆向指针 Ψ； ( ) = arg max ^(<?；) ·α_{¾ ι¾} 这里所使用的术语局部概率和逆向指针 Ψ_;( )都是沿用 Viterbi 算法的经典定义。关于该参数的定义的详细描述，可以参见 Lawrence R. Rabiner, PROCEEDINGS OF THE IEEE, Vol.77, No.2, 1989年 2月，通过参照将其全文并入本文。

第四步，确定最终状态，并回溯最优路径：

确定最终状态， =01¾11^ ^( ) 按照逆向指针回溯最优路径，即最可能胎儿单基因型 = _;( ) = l,2,3,...,N- 1)。第五步，输出结果。

由此，能够有效地对胎儿基因组的序列进行分析。相比其他已有的产前检测技术方法，本方法有以下技术优势，主要体现在准确性和可获得的遗传信息量上：

1) 根据本发明实施例，检测的位点不仅针对父源性位点，对母源性位点，即母亲杂合位点，同样可以很好地检测出来胎儿是否遗传了母亲的致病位点，检测准确率可高达 95% 以上，且可以检测多种变异类型，扩大了疾病检测的范围。

2) 根据本发明实施例，不但可以通过一次测序获得多个位点、多种疾病的信息，对于一些在孕妇血浆中覆盖程度较低，单纯提高测序深度依然无法准确判定的基因序列，可以通过本方法推断得到，结果准确可靠。

3) 根据本发明实施例，可进行遗传疾病作图，对于一些连锁相关疾病，可通过其他位点的信息直接推断出来，一次可获得的信息量大，对临床检测更加具有指导意义。

另外，根据本发明实施例，本发明的确定胎儿基因组中预定区域碱基信息的方法，不仅限于 SNP或者 STR等某一种遗传多态性位点，对所有的遗传多态性位点均可适用，且可以多种位点同时使用，以便互相验证。除了可进行产前无创检测胎儿基因组信息，达到疾病检测的目的，还可以进行无创产前亲子鉴定，在孩子出生前判定孩子父亲身份，为一些涉及抚养责任和义务、财产纠纷、性侵案等协助侦破。用于确定胎儿基因组中预定区域信 , 的***

在本发明的又一方面，本发明提出了一种用于确定胎儿基因组中预定区域核酸序列的 ***。根据本发明的实施例，参考图 2, 该*** 1000可以包括：文库构建装置 100、测序装置 200以及分析装置 400。

根据本发明的实施例，文库构建装置 100适于针对胎儿基因组 DNA样本，构建测序文库。才艮据本发明的实施例，测序装置 200与文库构建装置 100相连，并且适于对所构建的测序文库进行测序，以便获得测序结果，所得到的测序结果由多个测序数据构成。根据本发明的实施例，还可以进一步包括 DNA样本分离装置，该 DNA样本分离装置适于从孕妇外周血中提取胎儿基因组 DNA样本。由此，该***可以适用于进行无创产前检测。

根据本发明的实施例，任选地，还可以包括比对装置 300。根据本发明的实施例，比对装置 300与测序装置 200相连，并且适于将所得到的测序结果与参照序列进行比对，以便确定与预定区域对应的测序数据。根据本发明的实施例，可以用于进行测序的方法和设备并不受特别限制，包括但不限于双脱氧链终止法；优选高通量的测序方法，由此，能够利用这些测序装置的高通量、深度测序的特点，进一步提高测序效率。从而，提高后续对测序数据进行分析，尤其是统计检验分析时的精确性和准确度。所述高通量的测序方法包括但不限于第二代测序技术或者是单分子测序技术。所述第二代测序平台（Metzker ML. Sequencing technologies-the next generation. Nat Rev Genet. 2010 Jan;ll(l):31-46 ) 包括但不限于 Illumina-Solexa ( GA™,HiSeq2000™等）、 ABI-Solid和 Roche-454 (焦碑酸测序）测序平台；单分子测序平台（技术）包括但不限于 Helicos公司的真实单分子测序技术（ True Single Molecule DNA sequencing ) , Pacific Biosciences 公司单分子实时测序 ( single molecule real-time (SMRT™) ), 以及 Oxford Nanopore Technologies公司的纳米孔测序技术等（Rusk, Nicole (2009-04-01). Cheap Third-Generation Sequencing. Nature Methods 6 (4): 244-245 )。随着测序技术的不断进化，本领域技术人员能够理解的是还可以釆用其他的测序方法和装置进行全基因组测序。根据本发明的具体示例，可以利用选自 Illumina-Solexa、 ABI-SOLiD、 Roche-454和单分子测序装置的至少一种对所述全基因组测序文库进行测序。根据本发明的实施例，可以使用的参照序列的类型并不受特别限制，可以为任何含有感兴趣区域的已知序列。根据本发明的实施例，可以釆用已知的人类参考基因组作为参照序列。例如，根据本发明的实施例，釆用的人类参考基因组为 NCBI 36.3 , HG18。另外，根据本发明的实施例，进行比对的方法并不受特别限制。根据本发明的具体实施例，可以釆用 SOAP进行比对。

根据本发明的实施例，分析装置 400适于基于所述胎儿的测序结果，结合胎儿遗传相关个体的遗传信息，根据隐马尔可夫模型，确定所述预定区域的碱基信息。

根据本发明的实施例，惠特比算法釆用 0.25作为初始状态概率分布，釆用 re/N作为重组概率，其中 re=25~30, 优选 25 , N为所述预定区域的长度，

(i - p>- ) ' Pr ^x> = ¾ -; ≠ i-_; or x,≠ x__; = ¾_；

作为重组转移矩阵， _ΡΓ=Γε/Ν。

根据本发明的实施例，将所述测序结果与参照序列进行比对，以便确定与所述预定区域对应的测序数据进一步包括按照下列公式确定概率最高的碱基： base = ∑ ^ (1 - Δ (base, m_k ) + ^ ε - Α (base, m ) + 1 ^ · Δ (base, /_¾ )

Ae{0,lj l — e 二

关于数据分析部分，前面已经进行了详细描述，也当然地适用于确定胎儿基因组中预定区域核酸序列的***。不再赘述。

由此，利用该***，能够有效地实施前面所述的确定胎儿基因组中预定区域核酸序列的方法，可以通过例如惠特比算法（Viterbi algorithm ), 借助隐马尔可夫模型，确定胎儿基因组中特定区域的碱基信息，由此，可以有效地对胎儿基因组的遗传信息进行产前检测。

此外，根据本发明的实施例，预定区域为已知存在遗传多态性的位点，而遗传多态性为选自单核苷酸多态性和 STR的至少一种。

在本文中所述的术语 "相连" 应作广义理解，既可以是直接相连，也可以是间接相连，只要能够实现上述功能上的衔接即可。

需要说明的是，本领域技术人员能够理解，在前面所描述的确定胎儿基因组中预定区域核酸序列的方法的特征和优点也适合于确定胎儿基因组中预定区域核酸序列的***，为描述方便，不再详述。计算机可读介质

在本发明的又一方面，本发明提出了一种计算机可读介质。根据本发明的实施例，计算机可读介质上存储有指令，所述指令适于被处理器执行以便基于胎儿的测序结果，结合胎儿遗传相关个体的遗传信息，根据隐马尔可夫模型，确定所述预定区域的碱基信息。由此，利用该计算机可读介质，能够有效地实施前面所述的方法，从而可以通过例如惠特比算法（Viterbi algorithm ), 借助隐马尔可夫模型，确定胎儿基因组中特定区域的碱基信息，由此，可以有效地对胎儿基因组的遗传信息进行产前检测。

根据本发明的实施例，指令适于按照惠特比算法，根据隐马尔可夫模型，确定所述预定区域的碱基信息。根据本发明的实施例，在所述惠特比算法中，釆用 0.25作为初始状态概率分布，釆用 re/N作为重组概率，其中 re=25~30, 优选 25, N为所述预定区域的长度，釆用 t i ~ Pr ) ^Χϊ = ^X;-i -) ~ ϊ-!

ajk ■Pr( = (1 ^ ,., , ≠j i or x,≠x_i_ _ 作为重组转移矩阵，其中， =re/N

根据本发明的实施例，所述指令将所述测序结果与参照序列进行比对以便确定与所述预定区域对应的测序数据进一步包括按照下列公式确定概率最高的碱基：

P_{i se} = ∑ (1 - Δ (base, m_k ) + ^ ε - Α (base, m ) + 1 ^ · Δ (base, /_¾ ) ，关于数据分析部分，前面已经进行了详细描述，也当然地适用于确定胎儿基因组中预定区域核酸序列的***。不再赘述。

就本说明书而言， "计算机可读介质" 可以是任何可以包含、存储、通信、传播或传输程序以供指令执行***、装置或设备或结合这些指令执行***、装置或设备而使用的装置。计算机可读介质的更具体的示例（非穷尽性列表）包括以下：具有一个或多个布线的电连接部（电子装置），便携式计算机盘盒（磁装置），随机存取存储器（RAM ), 只读存储器 ( ROM ), 可擦除可编辑只读存储器（EPROM 或闪速存储器），光纤装置，以及便携式光盘只读存储器（CDROM )。另外，计算机可读介质甚至可以是可在其上打印所述程序的纸或其他合适的介质，因为例如可以通过对纸或其他介质进行光学扫描，接着进行编辑、解译或必要时以其他合适方式进行处理来以电子方式获得所述程序，然后将其存储在计算机存储器中。

应当理解，本发明的各部分可以用硬件、软件、固件或它们的组合来实现。在上述实施方式中，多个步骤或方法可以用存储在存储器中且由合适的指令执行***执行的软件或固件来实现。例如，如果用硬件来实现，和在另一实施方式中一样，可用本领域公知的下列技术中的任一项或他们的组合来实现：具有用于对数据信号实现逻辑功能的逻辑门电路的离散逻辑电路，具有合适的组合逻辑门电路的专用集成电路，可编程门阵列（PGA ), 现场可编程门阵列（FPGA )等。

本技术领域的普通技术人员可以理解实现上述实施例方法携带的全部或部分步骤是可以通过程序来指令相关的硬件完成，所述的程序可以存储于一种计算机可读存储介质中，该程序在执行时，包括方法实施例的步骤之一或其组合。

此外，在本发明各个实施例中的各功能单元可以集成在一个处理模块中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个模块中。上述集成的模块既可以釆用硬件的形式实现，也可以釆用软件功能模块的形式实现。所述集成的模块如果以软件功能模块的形式实现并作为独立的产品销售或使用时，也可以存储在一个计算机可读取存储介质中。

下面将结合实施例对本发明的方案进行解释。本领域技术人员将会理解，下面的实施例仅用于说明本发明，而不应视为限定本发明的范围。实施例中未注明具体技术或条件的，按照本领域内的文献所描述的技术或条件（例如参考 J.萨姆布鲁克等著，黄培堂等译的《分子克隆实验指南》，第三版，科学出版社）或者按照产品说明书进行。所用试剂或仪器未注明生产厂商者，均为可以通过市购获得的常规产品，例如可以釆购自 Illumina公司。

一般方法

本发明实施例的主要步骤包括：

1) 无创釆取含有胎儿遗传物质的孕妇样品，提取其中含有的遗传物质。

2) 胎儿家庭成员如父母和外祖父母等基因组 DNA提取和纯化。

3) 各遗传物质根据不同测序平台测序要求进行文库构建。

4) 测序获得的数据进行过滤，过滤条件根据质量值、接头污染等来设定。

5) 获得的高质量序列根据需要进行组装处理，组装结果与人类基因组参考序列进行比对。

获得唯一比对的序列，带入模型进行分析。分析模型：

记号：

I.需要检测的位点数为 N。

II. 父母的单倍型分别记为 FH = 和 MH = {m , mh_x \ mh_k = {m_lk,...,m_lk,...,m_Nk}, fh_k = {f_lk,...,f_ik,...,f_Nk} ^ ^_ik,mh_ik e{A,C,G,T} ^ ke{0,l], i = 2,3,...,N。

III. 将未知的胎儿单倍型记为 H = {¾， ¾}, 特别地，。和/ ¾分别遗传自母亲和父亲。

K = {m _l,...,m_Xi,...,m_N^} , = {/₁ ，…， ·, ···，/ 其中， x_; e{0，l}， y_{ e{0,l}₅ 下标 ,.和组成的序列对， q_t = {χ }^组成了我们需要解码的隐藏状态，而所有可能出现的隐藏状态组成了集合 β。

IV. 测序数据记为， S={s₁,...,s_i,...,s_N} 其中，； ^ = .,_Α， C^G^.G}, 代表此位点的测序信息，包含了 ACGT四种碱基的数量。

V. 平均胎儿浓度和平均测序错误率分别记为 f和 ί？。

第一步，构建初始状态概率分布向量，以及单倍体重组转移矩阵：

I. 初始状态概率分布记为； τ={ ( j eQ ), 根据本发明的实施例，在没有参考数据的情况下，可以设^ ？！^^ 全，即每种隐藏状态在第一个位点出现的可能性相等。

III. 根据本发明的实施例，记单倍体重组转移矩阵记为 A = ( j,keQ ), 其中为隐藏

第二步，构建观察序列概率矩阵:

根据本发明的实施例，记观察序列概率矩阵为 Β = [b_j (s,.)} ( i = 1,2,3, N , j eQ ), 其中 .(5_;)代表 "在位点 " 考虑母亲单倍型和胎儿单倍型（状态 j) 时，观测到这种测序信息的可能性"，即

(ⁿi + ⁿi,c + ⁿi,c + ^ηί,τ ) ! ίρ γ (p (p γ 其中代表 '在位点 , 考虑母亲单倍型和胎儿单倍型（状态 ) 时，该碱基出现的可能性"，即

^丄 ( 1 _ Δ [base, m_k) +— ε- (base, m_x^ +— s- {base, f_y> )

Ae{0,l}

其中指示函数 l-e x=y

第三步，构建局部概率矩阵，和逆向指针 (下面以一维局部概率矩阵构建为例；):

定义局部概率 ) = (max^. (¾) · α_{¾ ι¾} ) · (_Si ) 定义逆向指针 Ψ； ( ) = arg max S^q^-a 第四步，确定最终状态状态，并回溯最优路径:

确定最终状态， q_N*= argmax^ (¾ 根据逆向指针回溯最优路径，即最可能胎儿单基因型第五步，输出结果。实施例 1

样品收集及处理：

( 1 )所取样品包括一个家庭中父亲和母亲孕期的外周血，胎儿出生后取脐带血，以 EDTA抗凝管收集， Oragene®DNA唾液釆集 /DNA纯化试剂盒 OG-250釆集祖父母和外祖父母唾液。

(2)祖父母和外祖父母唾液 DNA提取后用 Infinium® HD Human610-Quad BeadChip 基因芯片进行分型。

(3)取母亲孕期外周血， 1600g, 4°C离心 10分钟，将血细胞和血浆分开，血浆再以 16000g, 4°C离心 10分钟，进一步去除残留的白细胞。孕妇血浆用 TIANamp Micro DNA Kit

(TIANGEN)提取 DNA, 得到母亲和胎儿基因组 DNA混合物，并且从全血分离血浆后剩余的白细胞中提取母亲基因组 DNA。将所得到的血浆 DNA根据 Illumia®公司 HiSeq2000^TM 测序仪的上机要求进行建库，构建好的文库经 Agilent®Bioanalyzer 2100检测片段分布范围符合要求，再经过 Q-PCR方法对两个文库进行定量，合格后 Illumina® HiSeq2000™测序仪测序，测序循环数为 PElOlindex (即双向 lOlbp index测序），其中仪器的参数设置及操作方法;^按照 Illumina®操作手册 (可由 http：〃 www.illumina.com/support/documentation.ilmn 获取）。

( 4 ) 父亲外周血、母亲外周血白细胞和胎儿脐带血则直接用 TIANamp Micro DNA Kit ( TIANGEN )提取试剂盒提取基因组 DNA。

除血浆 DNA样品外，将所获得的所有 DNA样品，需用 Covaris™打断仪打断至 500bp 大小的片段。将获得的 DNA片段以及血浆 DNA样品根据 Illumia®公司 HiSeq2000^TM测序仪的上机要求进行建库，具体步骤如下：

末端修复反应体系：

10 x T4 多核苷酸激酶緩冲液 10 μΐ

dNTPs(lOmM) 4 μ1

T4 DNA聚合酶 5 μΐ

Klenow 片段 1 μΐ

T4 多核苷酸激酶 5 μΐ

DNA片段 30μ1

ddH₂0补齐至 100 μΐ

20 °C反应 30分钟后，使用 PCR Purification Kit(QIAGEN)回收末端修复产物。将所得到的产物最后溶于 34μ1的 ΕΒ緩冲液中。

末端添加碱基 Α反应体系：

10 X Klenow緩冲液 5μ1

dATP(lmM) ΙΟμΙ

Klenow (3'-5' exo— ) 3μ1

DNA 32μ1

37 °C温育 30分钟后，经 MinElute® PCR Purification Kit(QIAGEN)纯化并溶于 12μ1的 ΕΒ 中。

接头连接反应体系：

2x Rapid DNA连接緩冲液 25μ1

PEI Adapter oligomix(20uM) ΙΟμΙ

Τ4 DNA连接酶 5μ1

添加碱基 Α的产物 ΙΟμΙ

20°C反应 15分钟后，使用 PCR Purification Kit(QIAGEN)回收连接产物。将所得到的产物最后溶于 32μ1的 ΕΒ緩冲液中。

PCR反应体系：

接头连接反应产物 10 μΐ Phusion DNA Polymerase Mix 25 μΐ

PCR 引物（lO pmol/μΙ) 1 μΐ

Index N(10 pmol/μΐ) 1 μΐ

超纯水 13 μΐ

反应程序如下：

98 °C 30 s

98 °C 10 s

65 °C 30 s [ . 10个循环

72 °C 30 s J

72 °C 5 min

4°C Hold

使用 PCR Purification Kit(QIAGEN)回收 PCR产物。样品最后溶于 50μ1的 ΕΒ緩冲液中。构建好的文库经 Agilent®Bioanalyzer 2100检测片段分布范围符合要求，再经过 Q-PCR 方法对两个文库进行定量，合格后，用 Illumina® HiSeq2000™测序仪测序，测序循环数为 PElOlindex (即双向 lOlbp index测序；)，其中仪器的参数设置及操作方法都按照 Illumina® 操作手册 (可由 http：〃 www.illumina.com/support/documentation.ilmn获取 )。

( 5 ) 父母基因组测序分型：

a. 使用 SOAP2将测序数据比对到人类参考基因组（版本为 NCBI 36.3 , HG18 )。

b. 使用 SOAPsnp (南方汉族 ( CHS ) 家系数据使用的是千人计划数据）对数据进行一致序歹¹ J ( consensus sequence, CNS )构建。

c. 提取出标记位点的基因型。

( 6 ) 父母单倍体型推断：

a.构建含祖辈与父母基因型的群体基因型矩阵，即提取父母、祖辈和南方汉族家系在标 i己位点的基因型。

b. 使用 BEAGLE对父母的单倍型进行推断。

( 7 )胎儿单倍体型推断：

a. 用 SOAP2将血浆测序数据比对到人类参考基因组（版本为 NCBI 36.3 , HG18 )。 b. 构建初始状态概率向量，以及单倍体重组转移矩阵。

构建初始状态概率向量：釆取无参考数据模式，即各个初始状态概率相等，均为 0.25。单倍体重组转移矩阵：保守地，我们取 re = 25 (其余按一般方法所述）。 c 统计每个位点的测序信息，并构建观察序列概率矩阵（其余按一般方法所述）。 d. 构建局部概率矩阵，和逆向指针（其余按一般方法所述）。

e. 确定最终状态状态，并回溯最优路径。

f. 输出。

根据胎儿出生后的脐带血基因分型结果，我们的分类准确性统计如下：母亲

纯合 ^：合

位点致准靈准确举位点数准确:数准. i 率位点数 m 准确率

199,552 199,552 63^968 96,57% 265J90 263,520 99.15% 常染色体杂合 65,409 64,735 39,944 95.45% 107,258 104,679 97,60%

合计 264,961 264,287 99.75% 1OS₍0S7 103₍912 96.14% 373舞 9870%

X染色体 4,881 4,881 100.00% 1₍71S 1478 86.03% 6,599 6,359 96.36% 工业实用性

本发明的确定胎儿基因组中预定区域碱基信息的方法、用于确定胎儿基因组中预定区域碱基信息的***以及计算机可读介质，能够有效地应用于对胎儿基因组中预定区域的核酸序列进行分析。尽管本发明的具体实施方式已经得到详细的描述，本领域技术人员将会理解。根据已经公开的所有教导，可以对那些细节进行各种修改和替换，这些改变均在本发明的保护范 o <

围之内。本发明的全部范围由所附权利要求及其任何等同物给出。

在本说明书的描述中，参考术语 " 个实施例" 些实施例"、 "示意性实施例"、 "示例"、 "具体示例"、或 "一些示例" 等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

Claims

权利要求书

1、一种确定胎儿基因组中预定区域碱基信息的方法，其特征在于，包括下列步骤：针对胎儿基因组 DNA样本，构建测序文库；

对所述测序文库进行测序，以便获得胎儿的测序结果，所述胎儿的测序结果由多个测序数据构成；以及

基于所述胎儿的测序结果，结合胎儿遗传相关个体的遗传信息，根据隐马尔可夫模型，确定所述预定区域的碱基信息。

2、根据权利要求 1所述的方法，其特征在于，所述胎儿基因组 DNA样本是从孕妇外周血中提取的。

3、根据权利要求 1所述的方法，其特征在于，所述测序是利用选自 Illumina-Solexa、 ABI-Solid、 Roche-454和单分子测序装置的至少一种对所述测序文库进行的。

4、根据权利要求 1所述的方法，其特征在于，进一步包括将所述胎儿的测序结果与参照序列进行比对，以便确定来自于所述预定区域的测序结果。

5、根据权利要求 4所述的方法，其特征在于，所述参照序列为人类参考基因组。

6、根据权利要求 1所述的方法，其特征在于，所述胎儿遗传相关个体是所述胎儿的父母。

7、根据权利要求 1所述的方法，其特征在于，按照惠特比算法，根据隐马尔可夫模型，确定所述预定区域的碱基信息。

8、根据权利要求 7所述的方法，其特征在于，在所述惠特比算法中，釆用 0.25作为初始状态概率分布，釆用 re/N作为重组概率，其中 re=25~30 , 优选 25 , N为所述预定区域的长度，

釆 J_k γ, =

作为重组转移矩阵，其中， /^=re/N。

9、根据权利要求 4所述的方法，其特征在于，将所述胎儿的测序结果与参照序列进行比对，以便确定来自于所述预定区域的测序结果进一步包括按照下列公式确定概率最高的碱基：

P_{i se} = ∑ (1 - Δ (base, m_k ) + ^ ε - Α (base, m ) + 1 ^ · Δ (base, /_¾ ) 其巾， e

3 X≠ γ

10、根据权利要求 1 所述的方法，其特征在于，进一步包括：所述预定区域为已知存在遗传多态性的位点。

11、根据权利要求 10所述的方法，其特征在于，所述遗传多态性为选自单核苷酸多态性和 STR的至少一种。

12、一种用于确定胎儿基因组中预定区域碱基信息的***，其特征在于，包括：文库构建装置，所述文库构建装置适于针对胎儿基因组 DNA样本，构建测序文库；测序装置，所述测序装置与所述文库构建装置相连，并且适于对所述测序文库进行测序，以便获得胎儿的测序结果，所述胎儿的测序结果由多个测序数据构成；

分析装置，基于所述胎儿的测序结果，结合胎儿遗传相关个体的遗传信息，根据隐马尔可夫模型，确定所述预定区域的碱基信息。

13、根据权利要求 12所述的***，其特征在于，进一步包括 DNA样本分离装置，所述 DNA样本分离装置适于从孕妇外周血中提取胎儿基因组 DNA样本。

14、根据权利要求 12所述的***，其特征在于所述测序装置为选自 Illumina-Solexa、 ABI-Solid、 Roche-454和单分子测序装置的至少一种。

15、根据权利要求 12所述的***，其特征在于，进一步包括比对装置，所述比对装置与所述测序装置相连，用于将所述胎儿的测序结果与参照序列进行比对，以便确定来自于所述预定区域的测序结果。

16、根据权利要求 12所述的***，其特征在于，所述分析装置适于按照惠特比算法，根据隐马尔可夫模型，确定所述预定区域的碱基信息。

17、根据权利要求 16所述的***，其特征在于，所述惠特比算法釆用 0.25作为初始状态概率分布，釆用 re/N作为重组概率，其中 re=25~30, 优选 25, N为所述预定区域的长度, 釆用 x; 一 ί-ι 一

作为重组转移矩阵，其中， p_r=re/N。

18、根据权利要求 15所述的***，其特征在于，将所述胎儿的测序结果与参照序列进行比对，以便确定来自于所述预定区域的测序结果进一步包括按照下列公式确定概率最高的碱基： P_{i se} = ∑ - (1 - Δ (base, m_k ) + -s - A (base, m ) + - ^ · Δ (base, f_y> )

Ae{0,lj 2 2 2

其巾，

19、一种计算机可读介质，其特征在于，所述计算机可读介质上存储有指令，所述指令适于被处理器执行以便基于胎儿的测序结果，结合胎儿遗传相关个体的遗传信息，根据隐马尔可夫模型，确定预定区域的碱基信息。

20、根据权利要求 19所述的计算机可读介质，其特征在于，所述指令适于按照惠特比算法，根据隐马尔可夫模型，确定所述预定区域的碱基信息。

21、根据权利要求 20所述的计算机可读介质，其特征在于，在所述惠特比算法中，釆用 0.25作为初始状态概率分布，釆用 re/N作为重组概率，其中 re=25~30 , 优选 25 , N为所述预定区域的长度，

釆用

¾ = ^Pr(¾ = = γ,__χ

作为重组转移矩阵，其中， /^=re/N。

22、根据权利要求 19所述的计算机可读介质，其特征在于，所述指令将所述胎儿的测序结果与参照序列进行比对，以便确定来自于所述预定区域的测序结果。

23、根据权利要求 22所述的计算机可读介质，其特征在于，所述指令将所述胎儿的测序结果与参照序列进行比对，以便确定来自于所述预定区域的测序结果进一步包括按照下列公式确定概率最高的碱基：

P_{i se} = ∑ (1 - Δ (base, m_k ) + ^ ε - Α (base, m ) + 1 ^ · Δ (base, /_¾ ) 其巾，

₍, .. I l - e ,¾二 v

A(j，_ )二

24、根据权利要求 19所述的计算机可读介质，其特征在于，进一步包括：所述预定区域为已知存在遗传多态性的位点。 25、根据权利要求 24所述的计算机可读介质，其特征在于' 所述遗传多态性为选自单核苷酸多态性和 STR的至少一种。