本申请依照35U.S.C§119(e)的规定,要求下述美国临时专利申请的权利:系列号No.60/739882,申请日2005年11月26日;系列号No.60/742305,申请日2005年12月6日;系列号No.60/754396,申请日2005年12月29日;系列号No.60/774976,申请日2006年2月21日;系列号No.60/789506,申请日2006年4月4日;系列号No.60/817741,申请日2006年6月30日;系列号No.11/496982,申请日2006年7月31日;以及系列号No.60/846610,申请日2006年12月22日;上述公开的文献的全部内容在此引入作为参考。
具体实施方式
对该体系的概念性描述
本发明公开的体系的一个目标在于能够提供高准确度的遗传数据,用来进行遗传诊断。当来自于个体的遗传数据中含有大量的干扰信息或者错误信息时,本发明公开的体系能够利用相关个体的遗传数据间的相似性,以及二级遗传数据中包含的信息,对目标染色体组中的干扰信息进行清除。这通过下述方式来完成:确定在配偶子的形成过程中哪些染色体片段参加了反应,在减数分裂过程中在哪里发生了交联,进而能够认定哪些次级染色体组片段是与目标染色体组的片段几乎等同的。在某些情况中,这种方法可以被用来进行干扰性碱基对测量值的清除,但其也可以被用来对未被测量的个体碱基对或者DNA的全部区域的同一性进行推断。除此之外,可以对每个重新构建的分型(call)的置信度进行计算。首先使用一种高度简化的解释,进行非真实存在的设想,以阐述本发明所涉及的概念。随后再公开可以应用于当今技术中的详细的统计学方法。
本发明体系的另外一个目的是检测异常染色体的数量,染色体的片段,以及染色体的来源。在遗传样本中,异倍体、不平衡易位、单亲二体性、或者其他高水平的染色体异常、在各个基因座上存在的遗传物质的量,这些都可以用来测定所述样本的染色体状态。实现这一方法有许多种途径,本发明描述了其中的几种途径。在某些途径中,存在于样本中的遗传物质的剂量足够用来直接检测异倍体。在其他的途径中,可以使用清除遗传物质的方法,以增强对染色体不平衡的检测的有效性。可以对每个染色体的分型(call)的置信度进行计算。
本发明体系的另外一个目的是提供一种从遗传数据中有效、高效的分离出最简单、并且可理解的统计学模型的方式,该方式需要探测一系列宽范围的条目,用来模仿与遗传数据相关的变量所产生的效果。更具体的,当前可以利用的基于遗传数据的表型模仿方法或者表型易感性模仿方法中的绝大部分或者全部都具有下述缺点:(i)这类方法不使用凸集合最佳技术,因此不能保证在给定的一组训练数据中找到模型参数的局部最小解;(ii)这类方法不使用能够减化模型的复杂性的技术,因此,当结果的数量小于自变量的数量时,这类方法不能建立出具有优秀的代表性的模型;(iii)在逻辑回归当中,如果不对正常分布的数据进行简化的假设,这类方法就不能够从数据中提取出最简单的、可理解的标准(rules);(iv)这类方法不能平衡某些先验信息(a-priori information),因而不能对表型或者表型易感性做出最可能的预测,其中所述的先验信息是关于基因-基因间的相互关系、基因-表型间的相互关系、以及基因-疾病间的相互关系的信息;(v)这类方法只提供一种模型,因此不能够将不同的模型与训练数据进行交叉验证(cross-validating),以提供一种选择最可能的数据的常规途径。在基于对大量与遗传信息以及表型信息相关的数据组的分析而进行的结果预测中,这些缺点是致命的。总的来说,当前可利用的方法不能够有效的告知个体:当给定基因型时,特定的表型特征出现的可能性,或者当给定父母的基因型特征时,在其后代中特定的表型特征出现的可能性。
值得注意的是,下面给出的某些说明包括了本文作者先前公开发表的文章。这些文章被作为背景信息,以帮助理解和进一步说明本发明公开的内容。
可以把基因型-表型预测模型分为三个种类:i)已知遗传缺陷或者等位基因以100%的可能性引发疾病表型;ii)遗传缺陷以及等位基因增加了疾病表型的可能性,其中预测值的数量足够小,因而能够使用列联表对表型的可能性进行模仿;以及iii)利用多维线性回归模型或者多维非线性回归模型,使用遗传标记物的复杂组合对表型进行预测。在具有当前已知序列的359个基因(参见表1第2行)以及人类孟德尔遗传在线数据库(OMIM)中的疾病表型当中,绝大部分属于种类(i);并且剩余的部分主要属于种类(ii)。然而,随着时间的发展,期望获得多种属于种类(iii)的基因型-表型模型,其中需要对多个等位基因或者变异间的相互反应进行模仿,用以评价特定表型的可能性。例如,种类(iii)在现金被用在预测人类免疫缺陷病毒对抗逆转录治疗产生的应答性的情况中,其中所述的预测以人类免疫缺陷病毒的遗传数据为基础。
对于种类(i)来说,通常基于专家规则直接预测表型的发生。在一方面,描述了一类统计学技术,所述的统计学技术可以被用来针对种类(ii)进行准确的表型预测。在另一方面,描述了一类统计学技术,所述的统计学技术可以被用来针对种类(iii)进行准确的表型预测。在另外一个方面,描述了一类方法,通过所述的方法能够针对特定的表型、特定的一组整合数据、以及特定的个体数据选择出最佳的模型。
本发明公开的方法的某些实施方式利用列联表针对种类(ii)进行准确的预测。这类技术对涉及基因-基因间的相互关系以及基因-疾病间的相互关系的先验信息进行平衡,用以加强对表型的预测或者对表型易感性的预测。这类技术能够对先前研究得出的数据进行平衡(leverage),其中所述的先前的研究并没有对所有相关的自变量进行取样。这类技术并没有因为先前的结果缺少数据而将其丢弃,相反,它们对来自于HapMap(人类基因组单体型图)计划的数据进行平衡(leverage),并且在其他方面对这些先前的研究结果加以利用,其中所述的先前的研究结果仅仅对相关自变量的子集进行了测定。通过这种方式,可以根据所有的整合数据对预测模型进行训练,而不是简单的对来自于宿主的数据进行整合,其中针对该宿主的所有相关自变量均进行了测定。
本发明中描述的某些方法使用凸集合最佳技术来生成稀疏模型,所述的稀疏模型能够被用来对种类(iii)进行准确的预测。基因型-表型模型问题经常被完善,或者产生不适定(ill-posed),这是由于潜在的预测值——基因,蛋白,变异以及它们间的相互反应——的数量大于测定结果的数量。通过运用一个类似于Occam’s Razor(奥卡姆剃刀)的原理:当可以使用许多种可能的理论来解释一种现象时,最简单的理论最有可能是正确的,上述数据组仍然可以被用来训练具有准确的归纳性的稀疏参数模型。在上文中讨论的关于在种类(iii)中构建基因型-表型模型的方面具体描述了这一哲学体系。这里描述的应用于遗传数据的技术包括针对无法检测出的基因型-表型数据组或者不适定(ill-conditioned)的基因型-表型数据组建立稀疏参数模型。在对稀疏参数组进行选择时,运用类似于Occam’s Razor(奥卡姆剃刀)的原理,并且因此能够建立准确的模型,即使当潜在的预测值的数量大于测定结果的数量时。除此之外,这里描述的用于在种类(iii)中构建基因型-表型模型的技术的某些实施方式使用了凸集合最佳技术,对于一组给定的训练数据而言,该技术能够确保找到其模型参数的全局最小解。
对于一组给定的整合数据以及一组针对于个体的可利用数据而言,不能明确得知哪种预测途径最适合用来为该个体做出最佳的表型预测。除了对能够进行准确的表型预测的一组方法进行了描述之外,本发明公开的实施方式还公开了一种体系,对于一个给定的表型预测、一组给定的整合数据、以及针对于被预测个体的一组可利用的数据而言,该体系能够对多种方法进行检测,并且选择出最佳的方法。这里公开的方法以及体系利用多种模型以及多个调整参数,在一组给定的数据中检查不同的自变量与依变量的所有组合,随后,对自变量与依变量的组合进行选择,并且使用试验数据对那些调整参数进行测定,从而选择出能够达到最佳的模仿准确度的调整参数。在针对种类(i)的情况中,可以使用(draw)专家规则;在具有较少的自变量的其他情况中,例如在针对种类(ii)的情况中,列联表能够提供最佳的表型预测;在其他的情况例如在针对种类(iii)的情况中,可以利用线性回归技术或者非线性回归技术提供最佳的预测方法。值得说明的是,在阅读过本发明之后,本领域技术人员能够明确,用于选择针对个体进行预测的最佳模型的方法也可以被用来在许多模型技术中进行选择,其中所述的模型技术是超过本发明公开的技术范围的。
本发明所述技术的某些实施方式在几个方面被进行了验证。首先,在预测患有阿尔茨海默氏症的可能性方面对其进行了验证,该方法使用列联表以及从许多临床研究中整合得出的不完整的数据组进行预测,其中所述的临床研究以遗传标记物为基础,着眼于对阿尔茨海默氏症(早老性痴呆症)进行预测。接着,在模仿I型人类免疫缺陷病毒(HIV-1)的药物应答性的方面对本发明所述的体系进行了验证,该方法利用了回归分析以及病毒染色体中的遗传标记物的知识。最后,在对由它莫西芬以及伊立替康引起的副作用的预测方面对该体系进行验证,其中它莫西芬以及伊立替康是分别用在针对各种乳腺癌病例以及各种结肠癌病例的治疗当中的,该方法利用回归分析以及个体中的遗传标记物的不完整数据和相关癌症的实验室信息和临床信息的不完整数据。
由于基因定型测试的成本逐渐降低,能够从遗传数据中可靠的预测病毒药物应答、癌症药物应答、以及其他表型应答或者表型结果的统计学模型成为重要的工具,用于对适当的反应过程进行选择,这些反应过程可以是疾病的治疗、生活方式的决策或者生活习惯的决策、或者其他反应。这里描述的最佳技术能够被应用于许多基因型-表型模型问题当中,达到促进临床决策的目的。
本发明体系的技术性描述
清除数据:一个简单化的实施例
附图1描述了在减数分裂过程中发生的重组反应的过程,用于在亲本中生成配偶子。橙色部分(或者灰色部分)表示来自于个体的母亲的染色体101。白色部分表示来自于个体的父亲的染色体102。在这一被称作双线期的时间间隔内,在减数分裂的细胞分裂前期I中,可以看见由四个染色单体103组成的一个四分体。在被称为重组节的点104上,发生了同系对中的非姊妹染色单体间的交换。为了进行描述,本实施例将着眼于单个的染色体,以及三种单核苷酸多态性(SNPs),假定它们能够表示三种基因的等位基因特征。为了进行这个论题,假定可以分别测定母亲的染色体以及父亲的染色体的单核苷酸多态性。这一概念可以应用于多种单核苷酸多态性、具有多种单核苷酸多态性特征的许多等位基因、许多染色体中,也可以应用于当前的基因定型技术,在所述的基因定型技术中,在基因定型之前不能各自分离出母亲的染色体以及父亲的染色体。
需要注意能够发生感兴趣的单核苷酸多态性间的可能的交换的基因座。可以根据对应的单核苷酸多态性(SNP1,SNP2,SNP3),将三个母体基因的等位基因组描述为(am1,am2,am3)。将三个父亲基因的等位基因组描述为(ap1,ap2,ap3)。考虑到在附图1中形成的重组节,并且假定每一对重组染色单体仅发生一次重组。在这个过程中形成的配偶子组将具有下述等位基因:(am1,am2,ap3),(am1,ap2,ap3),(ap1,am2,ap3),(ap1,ap2,am3)。在不发生染色单体的交换的情况中,所述的配偶子将具有下述等位基因:(am1,am2,am3),(ap1,ap2,ap3)。当在相关区域内具有两个交换点时,所述的配偶子将具有下述等位基因:(am1,ap2,am3),(ap1,am2,ap3)。这八种等位基因的不同组合将被作为针对上述特定亲本所假设的等位基因组而被引证。
从胚胎DNA中测定的等位基因将是干扰性的。为了进行讨论,从胚胎DNA中分离出一个单个的染色体,并且假定该染色体来自于附图1中所描述的发生了减数分裂的亲本。可以使用指示变量的向量(vector)来描述对于该染色体上的等位基因的测定值:A=[A1 A2 A3]T,其中,当胚胎染色体中被测定的等位基因是am1时,A1=1,当胚胎染色体中被测定的等位基因是ap1时,A1=-1,并且当胚胎染色体中被测定的等位基因既不是am1也不是ap1时,A1=0。根据对假定亲本的等位基因组的假设,能够生成八个向量,分别对应于上述描述的所有可能的配偶子。对于上述描述的等位基因来说,这些向量可以是a1=[1 1 1]T,a2=[11 -1]T,a3=[1 -1 1]T,a4=[1 -1 -1]T,a5=[-1 1 1]T,a6=[-1 1 -1]T,a7=[-1 -1 1]T,a8=[-1 -1 -1]T。在这个高度简化的体系应用中,可以通过在假设组以及被测定的向量之间进行一个简单的相关性分析,来确定晶胚具有的可能的等位基因:
i*=argmaxiATai,i=1......8 (1)
一旦确定了i
*,则假定的
就被认定为是胚胎DNA中最可能的等位基因组。之后,使用不同的假设、假定上述胚胎染色体来自于母亲或者来自于父亲,将上述过程重复进行两次。能够产生最大的
相关性的假设就被认定为是正确的假设。每次利用一组假定的等位基因组,基于分别对母亲DNA或者父亲DNA的测量进行假定。值得注意的是,在本发明公开的方法的一个典型的实施方式中,在那些与特定的疾病表型相关的重要的单核苷酸多态性之中进行大量的单核苷酸多态性测定——其中所述的与特定的疾病表型相关的重要的单核苷酸多态性也被称为表型相关性单核苷酸多态性或者PSNP。通过从NCBI(美国国立生物技术信息中心)的单核苷酸多态性数据库(dbSNP)中选择出那些在个体之间具有本质区别的对照性单核苷酸多态性(RefSNP),可以通过推理(a-priori)从表型相关性单核苷酸多态性中选择出非表型相关性单核苷酸多态性(NSNP)(例如,通过建立具体的基因定型阵列)。或者,在表型相关性单核苷酸多态性中,可以针对一对具体的父母选择出非表型相关性单核苷酸多态性,这是因为在不同的父母中,这种单核苷酸多态性是不同的。在表型相关性单核苷酸多态性中的其他单核苷酸多态性可以用于进行具有高水平置信度的测定,测定在表型相关性单核苷酸多态性中是否发生了交换。很重要的一点是,应当注意,当使用这种表示方法表示不同的“等位基因”时,仅仅是出于方便的目的;其单核苷酸多态性不与编码蛋白的基因发生关联。
当前技术中的体系
在另外一个更加复杂的实施方式中,使用给定的一种特定的测量方式来计算一组等位基因的后验概率,该测量方式考虑了特定交换发生的可能性。除此之外,进行了scenario typical微阵列以及其他的基因定型技术,其中,通过针对一对染色体进行单核苷酸多态性的测量,而不是针对单个染色体进行测量。可以使用表示几对单核苷酸多态性测量值的随机变量(e1,i,e2,i),(p1,i,p2,i)以及(m1,i,m2,i),分别对胚胎、父亲染色体以及母亲染色体的基因座I处的基因型测量结果进行表征。如果所有的测量都是成对进行的,那么就无法确定交换反应发生在母亲染色体中还是父亲染色体中,因为这一原因,需要对该方法进行下述改良:除了对受精的晶胚、父亲二倍体组织以及母亲二倍体组织进行基因定型之外,还需要对来自于父母各方的一个单倍体细胞进行基因定型,其中所述的单倍体细胞是精子细胞以及卵细胞。使用p1,i、i=1......N表示精子细胞中测定的等位基因,使用p2,i表示从父亲的二倍体组织中测定的互补等位基因。类似的,使用m1,i表示从卵细胞中测定的等位基因,使用m2,i表示从母亲的二倍体细胞中测定的互补等位基因。上述测定方法不能够提供亲本染色体在生成标准的精子细胞以及卵细胞时发生亲本染色体发生交换的基因座。然而,可以假定卵细胞或者精子细胞上的N等位基因的序列是由亲本染色体通过发生少量的交换或者不发生交换而产生的。这一信息足可以应用本发明公开的运算法则。某些的误差概率与父亲的单核苷酸多态性分型以及母亲的单核苷酸多态性分型是存在联系的。对这种误差概率的估算需要依据对于(p1,i,p2,i)和(m1,i,m2,i)的测定以及所使用的技术的信号干扰比(信噪比)而发生变化。尽管可以针对每个基因座唯一的计算出这些误差概率且不影响本发明公开的方法,但可以通过下述方式将这一代数学问题在此进行简化:假定父亲的单核苷酸多态性分型的准确概率以及母亲的单核苷酸多态性分型的准确概率分别是常数pp以及pm。
假定在胚胎DNA上进行了M测量法。除此之外,将符号进行了少许的修改,此时的A表示一个组而不是一个向量:A表示一个对于来自于每个亲本的等位基因的组合(或者组)的特定的假设。来自于双亲的等位基因的所有可能的组合的组A用SA来表示。其目标在于利用M测量法确定具有最高后验概率的等位基因的组合(或者假设)A∈SA:
利用条件概率的规则,P(A|M)=P(M|A)P(A)/P(M)。由于P(M)在所有不同的A组中是共用的,因而最优的研究结果可以被重新修改为:
现在考虑对P(M/A)的计算。从单个的基因座i开始,并且假设位于晶胚上的该位点来自于双亲的单核苷酸多态性pt,1,i以及mt,1,i,其中下脚标t用来表示这些双亲单核苷酸多态性的真实值,该真实值与测量得到的p1,i以及m1,i相反,所述的p1,i以及m1,i可以是准确的或者是不准确的。胚胎的单核苷酸多态性的真实值用(et,1,i,et,2,i)来表示。如果假设A是真实值,那么(et,1,i,et,2,i)=(pt,1,i,mt,1,i)或者(mt,1,i,pt,1,i)。由于无法区分哪个测量值(e1,i,e2,i)来自于哪个亲本,因而两种顺序都需要考虑,因而假设的组A=[(pt,1,i,mt,1,i),(mt,1,i,pt,1,i)]。特定的测量值M取决于亲本的单核苷酸多态性的真实值或者其潜在的状态,即(pt,1,i,pt,2,i)以及(mt,1,i,mt,2,i)。由于具有四种单核苷酸多态性,pt,1,i,pt,2,I,mt,1,i,mt,2,i,并且上述每一个值都能够假定四个核苷碱基对A,C,T,G的值,因此存在44或者256个可能的状态。针对一种状态s1来描述上述运算法则,在该种状态中,假定pt,1,i≠pt,2,i≠mt,1,i≠mt,2,i。从这种解释中,可以明确怎样将这一方法应用到所有的256种可能的状态sk,k=1......256中。假定对胚胎的单核苷酸多态性(e1,i,e2,i)进行了M法测定,并且获得e1,i=p1,i,e2,i=m1,i的结果。对于给定的假设A以及状态s1而言,通过下述等式计算这种测量值的先验概率:
P(e1,i=p1,i,e2,i=m1,i|A,s1)=
P(et,1,i=pt,1,i,et,2,i=mt,1,i|A,s1)P(e1,i=p1,i|et,1,i=pt,1,i)P(e2,i=m1,i|et,2,i=mt,1,i)
+P(et,1,i=mt,1,i,et,2,i=pt,1,i|A,s1)P(e1,i=p1,i|et,1,i=mt,1,i,pt,1,i≠mt,1,i)P(e2,i=m1,i|et,2,i=pt,2,i,pt,2,i≠mt,1,i)
(4)
因为假定的A=[(pt,1,i,mt,1,i),(mt,1,i,pt,1,i)]使得胚胎单核苷酸多态性的两种顺序的几率相同,因而将第一个条目(term)以及第二个条目(term)间的第一个等式认为是:P(e1,i=p1,i,e2,i=m1,i|A,s1)=P(e1,i=m1,i,e2,i=p1,i|A,s1)=0.5。现在考虑第一个条目的第二个等式,P(e1,i=p1,i|et,1,i=pt,1,i),测定e1,i=p1,i的概率给出了下述假设:胚胎的单核苷酸多态性et,1,i实际上是来自于父亲的单核苷酸多态性pt,1,i的。准确测定父亲的单核苷酸多态性、母亲的单核苷酸多态性、以及胚胎的单核苷酸多态性的概率是pp、pm以及pe。给定假定(et,1,i=pt,1,i),对(e1,i=p1,i)的测量需要在准确测量胚胎的单核苷酸多态性以及父亲的单核苷酸多态性的前提下进行,或者在上述两者都不是被准确测量的,并且是由于具有相同的核苷(A,C,T或者G)才导致它们产生不准确测量性的前提下进行。因此,P(e1,i=p1,i|et,1,i=pt,1,i)=pepp+(1-pe)(1-pp)/3,其中为了简单化而假设对四种核苷的不正确分型的概率是相同的——这一运算法则可以被容易的进行改进,从而适应一个特定核苷(A,C,T,G)分型的不同的概率,在另外一个特定核苷上进行测定。相同的方法可以被应用到第一个条目的第三个等式中,从而得到P(e2,i=m1,i|et,2,i=mt,1,i)=pepm+(1-pe)(1-pm)/3。现在考虑第二个条目的第二个等式,P(e1,i=p1,i|et,1,i=mt,1, i,mt,1,i≠pt,1,i)需要e1,i或者p1,i是不准确的测量值,或者需要两者都是不准确的测量值,这样一来测量的值就是相等的:P(e1,i=p1,i|et,1,i=mt,1,i,mt,1,i≠pt,1,i)=pe(1-pp)/3+(1-pe)pp/3+(1-pe)(1-pp)2/9。相同的论据可以被应用到第二个条目的最后一个等式中去,得到P(e2,i=m1,i|et,2,i=pt,2,i,mt,1,i≠pt,2,i)=pe(1-pm)/3+(1-pe)pm/3+(1-pe)(1-pm)2/9。现在,将这些所有的条目进行组合,并且进行下述假设——仅仅是为了简化上述代数方法——假定pe=pp=pm=p,因而可以计算:
尽管计算方法不同,但是这里描述的计算采用的相似的概念性方法可以被用于全部的256种可能的状态sk,k=1......256中。计算P(e1,i=p1,i,e2,i=m1,i|A,s1),得出全部256种状态下的si,将每个si相加求和,得到P(e1,i=p1,i,e2,i=m1,i|A)。换句话说,即为:
为了计算每种状态si的概率P(si),必须将构成一种状态的所有的单独的等位基因看做是分别的事件,这是因为这些等位基因分布在不同的染色体中,换句话说,即为:P(si)=P(pt,1,i,pt,2,i,mt,1,i,mt,2,i)=P(pt,1,i)P(pt,2,i)P(mt,1,i)P(mt,2,i).可以将Bayesian(贝叶斯定理)技术用来评估个体测量值的概率分布。可以将对位于母亲染色体基因座i上或者父亲染色体基因座i上的等位基因的每次测量看成是一次抛硬币试验,来测定每个等位基因是特定值(A,C,T或者G)的概率。在成人组织样本中进行这类测量,并且可以将其看成是完全可信的,即使针对每种单核苷酸多态性对等位基因对进行了测量,并且无法确定哪个等位基因来自于哪个染色体。使wp,1,i=P(pt,1,i),根据父亲染色体上的单核苷酸多态性i的概率为值pt,1,i。在下面的解释中,用w来替代wp,1,i。将父亲染色体上的单核苷酸多态性i的测量值表示为收集数据D(collecting data D)。可以为w建立一个概率分布p(w),并且根据贝叶斯定理p(w|D)=p(w)p(D|w)/p(D)在测量了数据之后对其进行更新。假定观察到了单核苷酸多态性i的等位基因n,并且这一特定的等位基因依照w出现了h次——换句话说,heads被观察到h次。可以通过二项分布来描述这种观察的概率
在收集数据之前,假定具有一个先验分布p(w),该值均匀分布在0至1之间。通过利用贝叶斯定理,能够直接表明p(w|D)的结果分布是一种下述形式的β分布:
其中c是正常化常数(normalizing constant)。然而,许多时候使用贝叶斯定理以及新的测量值对p(w|D)进行更新,它将继续具有如上所述的β分布。每次收集到一个新的测量值时,p(w)的估计值将获得更新。值得注意的是,对于不同的种类以及不同的性别都有一个不同的p(w)函数,这是因为在特定单核苷酸多态性中的不同等位基因的概率取决于这些种族以及性别的分组,其中所述的种类以及性别的分组与Hapmap(人类基因组单体型图)计划中的分组规则是相同的。为了计算P(si),将每个染色体之上的每个等位基因与估算的概率分布相联系起来,即pp,1,i(Wp,1,i),pp,2,i(Wp,2,i),pm,1,i(Wm,1,i)以及pm,2,i(Wm,2,i)。随后,可以针对每个个体分布计算出P(si)的最大后验概率估计。例如,使用wp,1,i *表示最大化的pp,1,i(wp,1, i)。可以根据下述等式计算出P(si)的最大后验概率估算值:
P(si)MAP=wp,1,i *wp,2,i *wm,1,i *wm,2,i *(9)
由于针对每个w都具有概率分布,因而也可以计算任意特定置信度水平上的P(si)值的保守估计值,这种计算通过对概率分布进行整合,而不是简单的利用最大后验概率估计值。例如,可以通过上述方法在某种置信度水平上对p(M|A)进行保守性的估计。无论使用保守性估计值还是使用最大后验概率估计值,都对P(si)值进行了不断的精确,用于计算p(M|A)。在随后的步骤中,可以将参考的假定状态排除,以简化该符号,假定使用状态s1来对详细的计算过程进行解释。需要记住的是,这些计算实际上应当针对256种状态中的每一种状态来进行,随后将每种状态的概率进行求和。
现在,将计算p(M|A)的方法扩展到多重单核苷酸多态性基因座中,假定M表示的是晶胚上的N对单核苷酸多态性的测量值的组合,M=[M1,......,MN]。同样假定A表示的是具有上述单核苷酸多态性的亲本染色体的每种单核苷酸多态性的假设值的组合,A=[A1,......,AN]。用SA’表示其他可能的假设值的组合,所述的假设值的组合不同于A或者在A’组合中。可以通过下式对P(M|A)以及P(M|A’)进行计算:
考虑对P(A)的计算。实质上,这是根据形成晶胚的配偶子形成过程中特定的交换反应发生的可能性来计算的。特定的等位基因组的概率取决于两个因素,即胚胎染色体来自于母亲或者父亲的概率,以及特定的交换组合的概率。对于一组不存在异倍体的正常的胚胎染色体组来说,胚胎染色体来自于母亲或者父亲的先验概率是大约50%,因而在所有的A中都是同样的。现在,考虑一组特定的重组节的概率。相关的重组位点R的数量取决于测定的单核苷酸多态性的数量:R=N-1。由于构成感兴趣的表型相关性单核苷酸多态性(PSNP)周围的非表型相关性单核苷酸多态性(NSNP)的DNA片段相对较短,交换干扰使得在一个区域内的相同染色体上发生两次交换是高度不可能的。由于计算效率的原因,假定在每个区域的每个相关染色体中仅发生一次交换反应,且该交换反应可以发生在R可能性位点上。本领域技术人员能够显而易见的是,可以将这一方法扩展到包括在一个给定区域内发生多次交换反应的可能性。
将每个区域内在单核苷酸多态性之间发生的交换反应的概率用Pr来表示,r=1......N-1。在第一种顺序中,在区域r中的两个单核苷酸多态性之间形成重组节的概率与这些单核苷酸多态性之间的遗传距离(以cMorgans为测量单位)是成比例的。然而,近期的很多研究都能够对两个单核苷酸多态性基因座之间的重组概率进行精确的模仿(modeling)。对精子研究以及遗传变异类型的观察表明,重组率以千碱基的规模进行广泛的变化,并且这类重组发生在重组热点区,且引起连锁不平衡(linkage disequilibrium),从而显示出似块状结构(block-like structure)。可以通过UCSC(加州大学圣克鲁斯分校)基因组注释数据库公开获得NCBI(美国国立生物技术信息中心)关于人类基因组重组率的数据。
可以单独使用各种数据组,也可以将其组合使用。最常用的两种数据组来自于Hapmap(人类基因组单体型图)计划以及来自于Perlegen Human Haplotype(Perlegen人类单体型图)计划。后者具有较高的密度,而前者具有较高的质量。参见附图2,为染色体1的位置1038423至位置4467775间的区域重组率,依据Hapmap Phase I(第一代人类基因组单体型图)的数据,释放(release)16a。这些速率是利用可逆转跳跃的(reversible-jump)的Markov Chain Monte Carlo(马尔可夫链蒙特卡尔理论)(MCMC)方法进行估算的,所述的MCMC方法可以在LDHat软件包中获得。状态空间(state-space)模型被认为是分段不变的重组率图表的分布。马尔可夫链除了研究每个片段201的速率之外,还研究速率改变位点的数量以及位置的分布。这些结果可以被用来进行对Pr的估算,通过使用单核苷酸多态性间的每个恒定片段的长度对重组率次数进行整合。在核苷202上发生的累计的重组率在附图2中用红色部分来表示。
用C表示一组指示变量Cr,这样一来,如果在区域r中发生了一次交换反应,则Cr=1,否则,Cr为零。如果没有发生交换反应,C0=1,否则,C0为零。由于假设在非表型相关性单核苷酸多态性区域内仅发生一次交换反应,因而在C组中仅有一个元素是非零的。因此,由C组表示的交换的概率被认为是:
在关于单核苷酸多态性1......N的A假设中,具有4个相关性的可能的交换。即,发生在下述中的可能的交换:i)形成晶胚的父亲染色体(由指示变量组Cpe表示),ii)形成排序的精子的父亲染色体(组Cps),iii)形成晶胚的母亲染色体(组Cme)以及iv)形成排序的卵细胞的母亲染色体(组Cee)。其他两个另外的假设是v)是否第一个父亲的胚胎单核苷酸多态性来自于pt,1, i或者pt,2,i,以及vi)是否第一个母亲的胚胎单核苷酸多态性来自于mt,1,i或者mt,2,i。由于发现了在单核苷酸多态性间的交换发生的概率在不同的种族以及不同的性别中是不相同的,用Pp,r来表示父亲染色体的不同的交换概率,用Pm,r来表示母亲染色体的不同的交换概率。因此,包含了Cpe、Cps、Cme、Cee的特定假设A的概率被表示为:
(12)
现在结合用于测定P(A)以及P(M/A)的方程式,上述方程式3中计算A*所必需的元素都已经被定义。因此,从具有高度错误倾向的胚胎单核苷酸多态性的测量当中确定交换反应繁盛的位置已经成为可能,并且因此能够以高置信度对胚胎的测量值进行清除。还可以在最佳的假设A*中确定置信程度。为了确定置信程度,必须找到几率比(odds ratio)P(A*|M)/P(A*’|M)。用于下述计算的工具已经全部在上文中进行了描述:
随后,用
来表示A
*的置信度。这一算式可以表明特定的假设A
*的置信度,但不能表明特定的单核苷酸多态性测定值的置信度。为了计算胚胎表型相关性单核苷酸多态性n测定值的置信度,必须建立所有的假设A组,所述的A组不改变这一单核苷酸多态性值。这个组被表示为
其对应于能够导致晶胚上的表型相关性单核苷酸多态性n具有与假设A
*预测的相同的值的所有假设。类似的,建立一个组
其对应于能够导致晶胚上的表型相关性单核苷酸多态性n具有与假设A
*预测的不同的值的所有假设。现在,可以计算正确分型的单核苷酸多态性的概率与不正确分型的单核苷酸多态性的概率相比而产生的几率比:
根据几率比
而计算的胚胎单核苷酸多态性n的特定分型的置信度可以进行如下计算:
值得注意的是,这项技术同样可以被用来检测例如单亲二体性(UPD)这样的缺陷,其中两个相同的染色体来自于同一个亲本,而不存在来自于另外一个亲本的染色体。在试图推导亲本染色体中的交换反应时,不存在能够以高置信度对数值进行充分解释的假设,如果允许进行的可以选择的假设包括对单亲二体性的可能性的假设,则更有可能进行上述推导。
对不确定的重组率以及单核苷酸多态性测定可靠性的限制
这里公开的方法取决于:对特定单核苷酸多态性间的重组概率进行的假定;对胚胎、精子、卵细胞、父亲以及母亲染色体上每种单核苷酸多态性的正确测量的概率进行的假定;以及对不同种族中的某些等位基因的可能性进行的假定。这些假定中的每一个进行考虑:重组的机理并没有被充分理解和模仿,并且交换概率会根据个体的基因型而发生改变。而且,测定重组率的技术表现出了本质上的可变性。例如,用来执行可逆转跳跃的马尔可夫链蒙特卡尔理论(MCMC)方法的LDHat软件包进行了一组假设,并且需要一组用户输入重组的机理以及特征。这些假设能够影响对单核苷酸多态性间的重组率的预测,这也正如各种研究得到的不同的结果所表明的。
可以预期的是,除了上述列出的假设之外,对于重组率的假设将对方程式15具有最大的影响。上述描述的计算应当是以对单核苷酸多态性间的交换概率的最佳估算Pr为基础的。此后,可以对Pr使用保守性估算,利用例如具有95%的置信度界限的重组率的值,朝着减小置信测量值P(正确分型的单核苷酸多态性n)的方向。所述的95%置信度界限可以从置信度数据中得出,其中所述的置信度数据来自于针对重组率的各种研究中,并且该值可以通过查看使用不同的方法从不同的分组中得到的公开数据之间的不一致水平来进行确证。
类似的,上述的95%置信度界限可以用来估算每个正确分型的单核苷酸多态性的概率:pp,pm,pe。这些数值可以根据包括在基因定型检测输出文件中的实际测量阵列强度来计算,并且结合测量技术可靠性的经验数据。值得注意的是,对于非表型相关性单核苷酸多态性而言,那些没有确立好的参数pp,pm,pe可以被忽略。例如,由于二倍体亲本数据是经过可靠测量的,因而可以忽略亲本单倍体细胞的非表型性单核苷酸多态性的测量值以及晶胚非表型性单核苷酸多态性的测量值中那些与亲本二倍体组织上的相关单核苷酸多态性中的任何等位基因对应不上的测量值。
最后,对不同种族的组中某些等位基因的可能性的假设进行考虑,这引起了对P(si)的计算。这类假设同样不会给本发明公开的方法带来大的影响,因为对亲本二倍体数据的测量是可靠的,即,对亲本样本的状态si进行的直接测量通常会获得具有高置信度的数值。虽然如此,可以利用如等式8所描述的针对每个w的概率分布来计算每种状态的概率的置信度界限P(si)。如上所述,可以计算出每个P(si)的95%置信度界限,朝着能够减小置信度测量值P(准确分型的单核苷酸多态性n)的保守性方向。
对P(准确分型的单核苷酸多态性n)的确定能够获知下述信息:了解在每个表型相关性单核苷酸多态性周围具有多少需要进行测量的非表型相关性单核苷酸多态性,以达到期望的置信度水平。
值得注意的是,有许多不同的途径来实现本发明公开的方法的观点,即双亲DNA的测量方法的组合,对一个或者一个以上晶胚的DNA的测量方法,以及对于减数分裂过程的先验知识,其目的在于得到针对胚胎单核苷酸多态性的较好的估算。本领域技术人员应该明了,当不同子集的先验知识是已知的或者是未知的,或者是已知的内容具有较高程度或者较低程度的确定性时,可以怎样应用类似的方法进行计算。例如,可以使用多个晶胚的测量值来增加测量的确定性,通过这种方法,可以对特定的晶胚的单核苷酸多态性进行分型,或者提供来自于双亲的丢失的数据。值得注意的是,通过这种测量技术,不需要对感兴趣的表型相关性单核苷酸多态性进行测量。即使这种测量体系没有对非表型相关性单核苷酸多态性进行测定,也仍然能够通过本发明公开的方法以高置信度的条件对其进行重新构建。
还需要注意的是,一旦确定了减数分裂过程中发生交换反应的位点,并且目标基因组的区域被确定在亲本DNA的相关区域内,那么不但能够推断出感兴趣的个体单核苷酸多态性的同一性,而且可以推断出由于等位基因脱失或者其他测量中发生的错误而导致的目标基因组中的DNA全区域的缺失。同样可以测量出双亲DNA中的插入以及缺失,以及使用本发明公开的方法来推断出它们存在于目标DNA中。
可以利用许多技术来改善如上所述的本发明公开的运算法则的计算复杂性。例如,可以仅仅选择或者主要选择那些在母亲以及父亲之间存在差别的非表型性单核苷酸多态性。另外一种方法是仅仅使用在空间上与表型相关性单核苷酸多态性相近的非表型性单核苷酸多态性,使得在感兴趣的非表型性单核苷酸多态性以及表型相关性单核苷酸多态性之间发生交换的机会最小化。也可以利用在染色体之上存在的非表型性单核苷酸多态性,这样一来能够使多种表型性单核苷酸多态性发生最大限度的覆盖。另外的方法是最初仅使用少量的非表型性单核苷酸多态性对交换反应发生的位置进行大致的确定,并且以有限程度的确定性进行确定。随后,可以利用其他的非表型性单核苷酸多态性对交换模型进行精制,并且增加对表型相关性单核苷酸多态性的准确分型的概率。用于大致确定规模的交换组合的数量被表示为NC,其中N是单核苷酸多态性的数量,C是交换反应的最大数量。因此,当C=4时,针对每个表型相关性单核苷酸多态性容纳大致N=100值,这在Pentium-IV(奔腾IV)处理器中仍然属于计算机可处理的。使用上面描述的方法以及其他用于加强计算效率的方法,可以容易的容纳N>100,C>4的情况。下面将描述这类方法中的一种。
值得注意的是,还有许多其他的方法能够对表型相关性单核苷酸多态性进行分型,并且生成对正确分型的表型相关性单核苷酸多态性的概率进行的估算,以特定的胚胎数据、双亲数据组,以及使用的运算法则为依据,不需要改变根本的观点。这一概率可以用来进行个体的决策,以及用来在体外受精或者早期非入侵性产前遗传诊断(NIPGD)情况中提供可靠的服务。
遗传数据清除运算法则中的递归解法
这里描述了本发明的另外一种实施方式,该实施方式涉及一种线性测量的运算法则。假定运算量是有限的,那么在使用本发明公开的方法时,运算长度可能是一个重要的因素。当进行计算时,任何的运算法都必须计算某些值,其中,计算的数量需要指数性的上涨,因而单核苷酸多态性的数量将变得难以处理。从时间的立场出发,随着单核苷酸多态性的数量线性增长的许多运算解法通常都是优选的,这是由于单核苷酸多态性的数量是巨大的。下面描述了这种方式。
针对所有可能的假设进行考虑的一种简单的方式必须将运算时间作为单核苷酸多态性数量的指数函数进行处理。如前所述,假设测定的数据是一组针对于晶胚、父亲染色体以及母亲染色体上的k单核苷酸多态性的测量值,即,M={M1,......Mk},其中Mi=(e1i,e2i,p1i,p2i,m1i,m2i)。如前所述,设想的空间是SH={H1,......,Hq}={所有的假设组},其中每个假设具有Hj=(Hj 1,......Hj k)的形式,其中Hj 1是对于snip(剪切)i的“小型”假设,具有Hj 1=(pi *,mi *)的形式,其中pi *∈{p1i,p2i}并且mi *∈{m1i,m2i}。存在4种不同的“小型(mini)”假设Hj 1,具体是:
Hj 11:(e1i,e2i)={(p1i,m1i)或者(m1i,p1i)}
Hj 12:(e1i,e2i)={(p1i,m2i)或者(m2i,p1i)}
Hj 13:(e1i,e2i)={(p2i,m1i)或者(m1i,p2i)}
Hj 14:(e1i,e2i)={(p2i,m2i)或者(m2i,p2i)}
其目的在于选择最具可能性的假设H*,如:
H*=argmaxH∈SHP(H|M)=argmaxH∈SHF(M,H),其中函数F(M,H)=P(H|M)
在空间SH中具有4k种不同的假设。用尽一切办法探测整个空间SH以试图找到最佳的假设,必须使用的运算法则应当是k的指数级O(exp(k)),其中k是涉及的单核苷酸多态性的数量。对于一个较大的k来说,即使k>5,这种运算将是非常缓慢并且不实用的。因此,更实用的办法是采取递归解法,在该解法中,将对大小k的问题处理为在常量时间内对大小(k-1)的问题的函数。这里表示的解决方法是按照k的线性级O(k)。
单核苷酸多态性数量的线性递归解法
从F(M,H)=P(H|M)=P(M|H)*P(H)/P(M)开始计算。随后,argmaxHF(M,H)=argmaxHP(M|H)*P(H),其目标在于在线性时间内解决P(M|H)*P(H)。假设M(s,k)=在单核苷酸多态性s至单核苷酸多态性k之间的测量值,H(s,k)=在单核苷酸多态性s至单核苷酸多态性k之间进行的假设,并且将表示方法进行简化为M(k,k)=Mk,H(k,k)=Hk=针对单核苷酸多态性k(k个单核苷酸多态性)的测量值以及假设值。如下所述:
同样的:
其中
并且PC(Hi-1,Hi)=在Hi-1,Hi之间发生交换的概率最后,计算k个单核苷酸多态性:
F(M,H)=P(M|H)*P(H)=P(M(1,k)|H(1,k))*P(H(1,k))
=P(M(1,k-1)|H(1,k-1))*P(H(1,k-1))*P(Mk|Hk)*PF(Hk-1,Hk)
因而化简为F(M,H)=F(M(1,k),H(1,k))=F(M(1,k-1),H(1,k-1))*P(Mk|Hk)*PF(Hk-1,Hk)
即,我们可以将对于k个单核苷酸多态性中的F值的计算简化为对于k-1个单核苷酸多态性的F值的计算。
在H=(H1,......,Hk)中,对于单核苷酸多态性k(k个单核苷酸多态性)的假设而言:
其中
将其总结为:
其中可以发现G是递归的:当n=2,......k时,
并且G(M(1,1),H1)=0.25*P(M1|H1)
其运算法则如下:
当n=1时,生成4个假设的H1i,用i=1,......,4来计算G(M1|H1i)。
当n=2时,生成4个假设的H2i,在常量时间内用i=1,......,4来计算G(M(1,2)|H2i)
使用下述公式:
当n=k时,生成4个假设的Hki,用i=1,...,4来计算G(M(1,k)|Hki),i,利用
在任何时候,都只存在4种假设需要记忆存储以及恒定数量的运行。因此,上述运算法则与指数级相反,是关于k的线性的,其中k是单核苷酸多态性的数量。
以线性时间解答P(M)
不是必须对P(M)进行解答来获得最佳的假设,因为对于所有的H来说它是恒定不变的。但是为了获得真实的、有意义的条件概率P(H|M)=P(M|H)*P(H)/P(M)的数值,则该结果必须得自P(M)。如上所述,我们可以这样来写:
其中
我们可以通过递归来解答W(M,H):
因此为了简便,将大小k的问题简化为大小(k-1)的问题,通过
并且
如前所述,当n=2......k时,递归的生成W(2),...,W(K)=W(M(1,k-1)|Hk),直到最后,可能能够得自
对于每个水平来说,仅存在四种不同的假设Hk,因此,上述运算法则也是关于单核苷酸多态性的数量k的线性运算。
线性时间的个体单核苷酸多态性置信度
一旦我们计算出最佳的假设H*=(H1 *,......,Hk *),我们就能够在最后的答案中得到对于每个单核苷酸多态性的置信度,即当i=1,......,k时的P(Hi *|M)。如前所述,P(Hi *|M)=P(M|Hi *)P(Hi *)/P(M)=W(Hi *,M)/P(M),其中P(M)是已知的。
即,当针对i-1个单核苷酸多态性、i个单核苷酸多态性、以及i+1个单核苷酸多态性至k个单核苷酸多态性的假设的假设值H被分解。如前所述:
并且
因此
P(H(1,k)=1/4*T(H(1,k))=1/4*T(H(1,i-1))*PF(Hi-1,Hi *)*PF(Hi-1,Hi *)*T(H(i+1,k))
其中
由此,可以表明:
再一次,将大小k的情形简化为两组较小的大小,尽管这样一来比未简化之前稍微复杂了一点。可以对每一组进行如下的计算
因此,当n=1,......,k时,m=k,......,1时,上述运算法则将针对4种不同的Hn、Hm中的每一种来计算W(M(1,n),Hn),W(M(m,k),Hm),随后在需要的时候将它们进行组合,从而计算W(M(1,k),Hi *),其中i=1,......,k。运算的数量仍然是线性的k。
在可以获得较小的数据组或者不同的数据组的情况中,本发明公开的方法在胚胎数据中的应用
在本发明体系的一种实施方式中,仅仅必须利用来自于双亲中的一方(假定是来自于母亲)的二倍体数据,可以使用或者不使用来自于双亲中的一方或者双方的单倍体数据,并且在上述情况中,这一数据具有已知的高确定性或者低确定性。例如,通常期望考虑到卵细胞捐献的疲劳特性(grueling nature),将存在不容易获得母亲单倍体数据的情形。在阅读过本说明书之后,本领域技术人员应当明了,对于有限的数据组来说,计算特定的单核苷酸多态性的可能性的统计学方法将被进行怎样的改进。
另外一种可以选择的方法利用来自于关系更加远的亲属的数据来弥补双亲中的一方或者双方的丢失了的二倍体数据或者单倍体数据。例如,由于已知一组个体染色体来自于他的或者她的双亲中的每一位,因此可以利用来自于外祖父外祖母的二倍体数据对丢失的母亲二倍体数据或者错误测量的母亲二倍体数据进行部分的重建。
注意这一方法的递归性质:知晓了双亲的单细胞单倍体数据的自然存在的干扰测量值,以及来自于适当的祖父祖母的二倍体数据和/或单倍体数据,就可以利用本发明公开的方法对双亲的单倍体数据进行清除,反之将提供更加准确的晶胚的基因定型。如何对这一方法进行改进以使其应用于这类情况中,对于本领域技术人员来说是显而易见的。
优选使用更多的信息而不是更少的信息,因为这样能够在给定单核苷酸多态性的情况下增加正确分型的几率,并且能够增加这些分型的置信度。但随之而来的是增加了反应体系的复杂性,因为使用到了其他的技术以及其他的数据来源。额外信息的来源有许多,利用这些信息来增加数据的可利用的技术也有许多。例如,具有许多以信息学为基础的方法,这类方法能够利用从Hapmap(人类基因组单体型图)数据中、或者从遗传数据的其他知识库中发现的相关性。除此之外,还存在许多生物学方法,它们能够对遗传数据进行直接测量,不然则需要使用计算机数据分析方法对数据进行再创造。例如,通过从二倍体细胞中分离出个体染色体的方法将原本不可利用的单倍体数据变为可测量的数据,上述方法是利用流式细胞术来分离荧光标记的染色体。或者,可以利用细胞融合术生成单等位基因杂交细胞,来影响从二倍体到单倍体的转化。
应用本发明公开的方法来选择哪个晶胚可能被植入
在一种实施方式中,本发明所述的体系可以被用来确定晶胚植入到母体内并发育成婴儿的可能性。从某个方面来说,晶胚植入的可能性是由晶胚的单核苷酸多态性、和/或其与母亲的单核苷酸多态性间的相互关系来决定的,本发明公开的方法在帮助进行晶胚的选择方面是很重要的,该方法以清除单核苷酸多态性数据为基础,对哪个晶胚能够被成功的植入进行了可靠的预测。为了对这种可能性进行最佳的预测,必须考虑到确定的晶胚基因型,或者结合考虑晶胚中的基因表达水平、母体中的基因表达水平、和/或确定的母体基因型。
除此之外,本领域已知异倍体晶胚的植入可能性较小,导致成功受孕的可能性较小,导致发育成健康儿童的可能性较小。因此,筛选异倍体是选择最有可能获得成功结果的晶胚的一个重要的方面。关于这一方法的更加详细的描述将在下面给出。
双亲单倍体数据的推导
在本方法的一种实施方式中,可能需要推导双亲的单体型,给出双亲二倍体数据的详细信息。可以通过多种方式完成这一过程。在最简单的情况中,可以通过对直系亲属(母亲,父亲,儿子或者女儿)对单个单倍体细胞进行分子检测,从而推断出单体型。这种情况对于本领域技术人员来说是很容易处理的事情,可以通过从分子检测测定出的二倍体基因型中减去已知的单体型,从而推导出其姊妹单体型。例如,如果一个特定的基因座是杂合的,一个未知的亲本单体型就是与已知的亲本单体型相反的等位基因。
在另外一种情况中,可以通过对个体的亲本单倍体细胞的分子生物学单体定型(单体型分析)来获知亲本的干扰性单倍体数据,其中所述的亲本单倍体细胞可以是例如精子细胞,或者是通过个体染色体来获知,所述的个体染色体可以通过各种方法包括磁性珠以及流式细胞术而被分离。在这种情况中,可以使用如上所述的相同步骤,不同之处在于确定的单体型与测定的单体型都是干扰性的。
还存在一类方法,能通过二倍体数据直接推导出单倍体数据组,该方法借助统计学方法并利用普通大众的已知的单体域(haplotype block)(例如在公开的人类基因组单体型图计划中建立的那些单体型)进行推导。单体域本质上是一系列相关的等位基因,这些等位基因重复出现在各种不同的人口当中。由于这些单体域是古老的并且是常见的,因此它们可以被用来从二倍体基因型中预测单体型。随后,在本发明公开的方法中将双亲推断的单体域作为输入数据进行输入,用来从晶胚中清除干扰数据。能够完成这项任务的公共可利用的运算法则包括不完美种系发生方法,基于共轭先验(conjugate priors)以及来自于人类遗传学的先验知识的贝叶斯方法。这些运算法则中的某些利用隐马尔可夫模型。一项研究利用public trio以及不相关的个体数据来证明这类运算法则在1MB的序列内的误差率低至0.05%。然而,正如所预料到的,当单体域稀少时,针对个体的准确率就较低。在一项估算中,计算方法在多达5.1%的基因座上定相(phase)失败,并具有20%的较小的等位基因频率。
在本发明的一种实施方式中,利用在体外受精周期中的不同晶胚中分离出的多个卵裂球中得到的遗传数据对双亲的单体域进行推断,具有更高的可靠性。
利用高生产量基因定型以及中等生产量基因定型筛选异倍体的技术
在本发明所述体系的一种实施方式中,测量到的遗传数据可以被用来检测个体中是否存在异倍体和/或镶嵌性。这里公开的几种方法利用中等生产量的基因定型或者高生产量的基因定型来检测染色体数量或者DNA片段的拷贝数,这种检测是在组织样本中的扩增的或者未扩增的DNA中进行的。其目的在于对利用不同的定量基因定型平台和/或定性基因定型平台检测某些种类的异倍体以及检测镶嵌性水平中表现出的可靠性进行评价,其中所述的定量基因定型平台和/或定性基因定型平台是例如ABI Taqman,MIPS(分子倒置探针),或者来自于Illumina、Agilent以及Affymetrix的微阵列。在许多这类情况当中,在与探针杂交之前,遗传物质在基因定型阵列中通过聚合酶链式反应进行扩增,用以检测特定等位基因的存在。这类检测在基因定性中的使用方法将在本发明的其他内容中进行描述。
下面描述筛选异常数量的DNA片段的几种方法,其中所述的异常是由缺失、异倍体和/或镶嵌性所引发的。这类方法被分成下列几组:(i)不进行等位基因分型的定量技术;(ii)平衡(leverage)等位基因分型的定性技术;(iii)平衡(leverage)等位基因分型的定量技术;(iv)利用遗传数据在每个基因座上的扩增的概率分布函数的技术。所有的技术都涉及到对给定的染色体的给定片段上的多个基因座进行的测量,从而确定上述给定的片段在目标个体的染色体组中出现的次数。除此之外,这类方法包括建立由一个或者一个以上假设组成的组,所述的假设是针对给定片段的出现次数而进行的;测量给定片段上的多个基因座上的遗传数据的数量;确定每个假设的相对概率,给出目标个体的遗传数据的测量值;并且使用与每种假设相关的相对概率来确定给定片段出现的次数。并且,这类方法都包括了建立一种组合的测量值M,M是一个关于多个基因座上的遗传数据的数量的测量值的计算函数(computed function)。在所有这类方法中,根据测量值M确定阈值,用来对每种假设Hi进行选择,并且估算需要进行测量的基因座的数量,以获得特定水平的针对每种假设的错误检验。
给定测量值M的每种假设的概率是P(H
i|M)=P(M|H
i)*P(H
i)/P(M)。由于P(M)与H
i是无关的,我们通过P(M|H
i)*P(H
i)来确定每种假设给出M的概率。在接下来的步骤中,为了简化分析并且对不同的技术进行比较,我们假定对于所有的{H
i}来说P(H
i)都是相同的,因而我们能够仅通过考虑P(M|H
i)来计算所有P(H
i|M)的相对概率。因此,我们对于阈值的确定以及需要测量的基因座的数量的确定都是基于选择错误的假设的特定概率而进行的,前提是在假设对于所有的{H
i}来说P(H
i)都是相同的情况下。本领域技术人员在阅读过本说明书之后应该知晓,怎样对这一方法进行改进,以适应在{H
i}组中不同的假设会产生不同的P(H
i)这一事实。在某些实施方式中,将阈值设定为使其选择假设
在这一假设中,对于所有的i值而言,都能得到最大的P(H
i|M)。然而,不必需将阈值设定为获得最大化的P(H
i|M),只是为了达到在{H
i}组中的不同假设间的误差检验概率的特定比率。
非常重要的一点是,要注意这里使用的检测异倍体的技术同样可以很好的用来检测单亲二体性、不平衡易位、以及进行染色体的性别鉴定(男性或者女性;XY或者XX)。上述所有观点关注的都是给定样品的特性(identity)以及给定样品中存在的染色体(或者染色体片段)的数量,因此也都是本发明描述的方法所运用的观点。怎样将本发明所述的任意一种方法进行扩展以用来对这些异常中的任意一种进行检测对于本领域技术人员来说将是显而易见的。
匹配滤过(matched filtering)的概念
这里所应用的方法与数字信号的最优化检测中应用的方法是类似的。可以利用施瓦茨不等式(Schwartz inequality)来表示该方法,其中,当存在正常分布的干扰时,使信号干扰比(信噪比)(SNR)最大化的最佳方式是建立一个理想化的匹配信号或者匹配滤过,所述的匹配信号或者匹配滤过与每种可能的无干扰的信号相对应,并且将这种匹配信号与接收到的干扰信号进行相互关联。这种方法需要已知可能的信号组,以及干扰的统计学分布——平均偏差以及标准偏差(SD)。这里描述的是在样本中检测染色体、或者DNA片段存在与否的常规方法。该方法在对全染色体进行检查的过程与对插入或者缺失的染色体片段进行检查的过程是没有区别的。它们都被叫做DNA片段。在阅读过说明书之后,应当明了这类技术应当被进行怎样的扩展,以用来进行多个种类的异倍体的检测以及性别确定,或者在晶胚、胎儿或者出生婴儿的染色体中检测插入以及缺失。可以将这一方法应用于宽范围的定量基因定型平台以及定性基因定型平台,包括Taqman,定量聚合酶链式反应(qPCR),Illumina阵列,Affymetrix阵列,Agilent阵列,分子倒置探针(MIPS)试剂盒等等。
一般性问题的描述
假定在发生了两种等位基因变异x以及y的单核苷酸多态性处放置探针。在每个基因座i中,i=1......N,相应于上述两个等位基因中的遗传物质的量来收集数据。在Taqman检测中,这类测量值可以是例如周期时间,Ct,在其中,每个等位基因特异性染料的水平都超过了(cross)阈值。本领域技术人员能够明确,这种方法可以被进行怎样的扩展,以用来进行每个基因座的遗传物质的剂量的不同测量,或者进行在一个基因座上的每个等位基因相对应的遗传物质的剂量的不同测量。对遗传物质剂量的定量测量可能是非线性的,在这种情况中,因为目标片段的存在而导致的特定基因座上的测量值的变化将取决于该样本中具有的上述基因座的其他拷贝的数量,其中所述的上述基因座的其他拷贝是来自于其他DNA片段的。在某些情况中,所使用的技术需要进行线性测量,这样一来,因为目标片段的存在而导致的特定基因座上的测量值的变化将不取决于该样本中具有的上述基因座的其他拷贝的数量,其中所述的上述基因座的其他拷贝是来自于其他DNA片段的。下面将描述一种方法,该方法可以将来自于Taqman或者定量聚合酶链式反应检测法的测量值线性化,但是用于使非线性测量值线性化的其他技术也有很多,它们可以被应用到不同的检测当中。
用数据dx=[dx1......dxN]来表示在基因座1......N上的等位基因x中的遗传物质的剂量的测量值。对于等位基因y,也使用类似的表示方法,dy=[dy1......dyN]。假定每个片段j上都具有等位基因aj[aj1......ajN],其中每个元素aji是x或者是y。用dx=sx+vx来描述等位基因x中的遗传物质的剂量的测量数值,其中sx是信号而vx是干扰。信号sx=[fx(a11,......,aj1)......fx(aJN,......,aJN)],其中fx是来自于等位基因组的测量值的图谱(mapping),而J是DNA片段拷贝的数量。干扰矢量(disturbance vector)vx是由错误测量而引起的,在非线性测量的情况中,是由除目标DNA片段之外的其他遗传物质的存在而引起的。假定错误测量是呈正态分布(normally distributed)的,并且大于由非线性而导致的干扰(参见线性测量部分),因而vxi≈nxi,其中nxi存在变异σxi 2,且矢量nx是呈正态分布的,~N(0,R),R=E(nxnx T)。现在,假定将某些滤波h应用到这一数据中,从而得到测量值mx=hTdx=hTsx+hTvx。为了使信号干扰比(信噪比)(hTsx/hTnx)最大化,利用匹配滤过h=μR-1sx来表示h,其中μ是度量常数。可以将用于等位基因x的方法重复的应用于等位基因y中。
方法1a:当每个基因座的平均偏差以及标准偏差已知时,利用定量技术对异倍体或者性别进行测量,其中并不进行等位基因的分型
在这个环节中,假定关于某个基因座上的遗传物质的剂量的数值与等位基因值无关(例如,使用定量聚合酶链式反应),或者上述数值仅仅针对的是在人群中具有100%外显率的等位基因,或者将每个基因座上的多个等位基因的数值进行组合(参见线性测量部分)用以测量在该基因座上的遗传物质的剂量。因此,在这一环节中,可以研究数值dx而忽略dy。同样假定存在两种假设:h0表示具有两个DNA片段的拷贝(这些拷贝通常不是相同的),h1表示仅仅具有一个拷贝。对于每种假设来说,可以分别将数值表示为dxi(h0)=sxi(h0)+nxi以及dxi(h1)=sxi(h1)+nxi,其中sxi(h0)表示当存在两个DNA片段时在基因座i上的遗传物质的预期预测值(预期信号),而sxi(h1)表示当存在一个DNA片段时在基因座i上的遗传物质的预期预测值。通过将假设h0的预期信号进行差分(differencing out),来构建每个基因座的测量值:mxi=dxi-sxi(h0)。如果h1是真实的,那么该测量的预期值是E(mxi)=sxi(h1)-sxi(h0)。利用如上所述的匹配滤过概念,将h设定为h=(1/N)R-1(sxi(h1)-sxi(h0))。其测量值被表示为m=hTdx=(1/N)∑i=1......N((sxi(h1)-sxi(h0))/σxi 2)mxi。
如果h1是真实的,预期值E(m|h1)=m1=(1/N)∑i=1......N((sxi(h1)-sxi(h0))2/σxi 2,并且m的标准偏差是σm|h1 2=(1/N2)∑i=1......N((sxi(h1)-sxi(h0))2/σxi 4)σxi 2=(1/N2)∑i=1...... N((sxi(h1)-sxi(h0))2/σxi 2。
如果h0是真实的,m的预期值是E(m|h0)=m0=0,并且m的标准偏差仍然是σm|h0 2=(1/N2)∑i=1......N((sxi(h1)-sxi(h0))2/σxi 2。
附图3描述的是确定假阴性检测的概率以及假阳性检测的概率的方法。假定将阈值t设定为m1和m0的中间值,这样能够得到相等的假阴性概率以及假阳性概率(这不需要进行如下所述的情形)。通过下述比例(m1-t)/σm|h1=(m1-m0)/(2σm|h1)能够确定假阴性的概率。可以利用“5-Sigma”统计学,这样一来假阴性的概率即为1-normcdf(5,0,1)=2.87e-7。在此种情况下,其结果为(m1-m0)/(2σm|h1)>5或者5或10sqrt(平方根)((1/N2)∑i=1......N((sxi(h1)-sxi(h0))2/σxi 2<(1/N)∑i=1......N((sxi(h1)-sxi(h0))2/σxi 2或者sqrt(平方根)(∑i=1......N(sxi(h1)-sxi(h0))2/σxi 2)>10。为了计算N的大小,可以通过整合数据计算出平均的信号干扰比(平均信噪比):MSNR=(1/N)∑i=1......N((sxi(h1)-sxi(h0))2/σxi 2。随后,可以从上述不等式中得出N:sqrt(N).sqrt(MSNR)>10或者N>100/MSNR。
可以在Taqman检测法进行数据测量中应用到这一方法,所述的Taqman检测是来自于Applied Biosystems(应用生物体系)的,利用X染色体上的48个单核苷酸多态性。对每个基因座的测量值是时间,Ct,它利用对应于这个基因座的孔中释放的模板(die)来超过阈值。样本0由每个孔中的大约0.3纳克(50个细胞)的总体DNA(total DNA)组成,其中所述的总体DNA来自于混合的女性提供者,这些宿主具有两个X染色体;样本1由每个孔中的大约0.3纳克的DNA组成,其中所述的DNA来自于混合的男性提供者,这些宿主具有一个X染色体。附图4以及附图5表示出了样本1以及样本0的测量值的柱状图。这些样本的分布被表示为m0=29.97;SD0=1.32,m1=31.44,SD1=1.592。由于这个数据来自于混合的男性样本以及混合的女性样本,因而由此得出的SD取决于在混合样本中的每个单核苷酸多态性中的不同等位基因的频率。除此之外,得到的某些SD将取决于不同的检测在每个单核苷酸多态性中的不同的效率,以及每个孔中吸入的染料的不同的剂量。附图6提供了针对男性样本以及女性样本的每个基因座上的测量值的差别的柱状图。男性样本与女性样本间的平均差为1.47,这种差值的SD是0.99。由于这个SD值仍将受到混合的男性样本以及混合的女性样本中的不同的等位基因频率的支配,因而它将不再受到在每个基因座上进行的每种检测的不同效率的影响。由于其目的在于利用大致相似的SD来对两个测量值中的每一种进行区分,因而可以将SD值调整至适合所有基因座上的每个测量值,为0.99/sqrt(2)=0.70。针对每个基因座进行两次试验,目的在于为在该基因座上进行的检测估算σxi,从而能够应用匹配滤过。将σxi的下限设置为0.2,是为了避免因为仅经过两次试验而计算出的σxi所导致的统计学异常。只有那些在两个等位基因中、在两次试验的运行中、在男性样本以及女性样本中都没有发生等位基因脱失的基因座(总计37个)才能够被用于上述方案以及计算中。将上述描述的方法应用于该数据中,从而得出MSNR=2.26,因此,N=2252/2.26^2=17个基因座。
方法1b:当每个基因座的平均偏差以及标准偏差未知时或者是相同时,利用定量技术对异倍体或者性别进行测量,其中并不进行等位基因的分型
当每个基因座的特点没有被很好的了解时,可以做出一个最简单的假设:在每个基因座上进行的所有的检测都将表现出类似的结果,即所有的基因座i的E(mxi)以及σxi都是常数,这样一来,可以仅仅考虑E(mx)以及σx。在这种情况下,匹配滤过的方法m=hTdx就变为计算dx分布的平均值。这种方法将被用来进行平均值的比较,并且该方法可以被用来估算利用真实数据的不同种类的检测方法所需要的基因座的数量。
如上所述,当在样本中存在两个染色体(假设h0)或者存在一个染色体(假设h1)时,需要对其类型进行考虑。对于h0来说,其分布为N(μ0,σ0 2),对于h1来说,其分布为N(μ1,σ1 2)。分别使用N0样本以及N1样本对每种分布进行测量,其中利用到了测量的样本的平均值以及SD值:m1,m0,s1,以及s0。可以使用呈正态分布的随机变量M0,M1对上述平均值进行模仿(modeled),M0~N(μ0,σ0 2/N0)并且M1~N(μ1,σ1 2/N1)。假定N1以及N0足够大(>30),因此可以假定M1~N(m1,s1 2/N1)并且M0~N(m0,s0 2/N0)。为了检验这些分布是否是不相同的,可以对平均值间的区别进行检验,其中d=m1-m0。随机变量D的变异是σd 2=σ1 2/N1+σ0 2/N0,其可以近似于σd 2=s1 2/N1+s0 2/N0。给定h0时,E(d)=0;给定h1时,E(d)=μ1-μ0。随后将对进行h1和h0分型的不同技术进行讨论。
通过不同次Taqman检测测定到的数值可以被用来进行校准工作,其中所述的检测利用了在X染色体上的48个单核苷酸多态性。样本1由每个孔中的大约0.3纳克的DNA组成,其中所述的DNA来自于混合的男性提供者,这些提供者具有一个X染色体;样本0由每个孔中的大约0.3纳克的DNA组成,其中所述的DNA来自于混合的女性提供者,这些提供者具有两个X染色体。N1=42并且N0=45。附图7以及附图8表示出了样本1以及样本0的柱状图。这些样本的分布被表示为m1=32.259,s1=1.460,σm1=s1/sqrt(N1)=0.225;m0=30.75,s0=1.202,σm0=s0/sqrt(N0)=0.179。这些样本的d=1.509并且σd=0.2879。
由于这个数据来自于混合的男性样本以及混合的女性样本,因而其大部分标准偏差取决于在混合样本中的每个单核苷酸多态性中的不同等位基因的频率。在多次试验中的同一时间中,对一个单核苷酸多态性的Ct(周期时间)的变化进行考虑,从而估算出SD。这一数值在附图9中有所表示。这个柱状图是关于0点对称的,因为每个单核苷酸多态性的Ct(周期时间)都是在两次试验中进行测量的,并且从其中减去了每个单核苷酸多态性的周期时间的平均值。在经过了两次试验的混合的男性样本的20个单核苷酸多态性中,平均的标准偏差是s=0.597,这一SD值将被保守的应用在男性样本和女性样本中,因为女性样本的SD值将小于男性样本的SD值。除此之外,值得注意的是,仅使用了一种染料进行测量,因为上述混合的样本被假定为对于所有的单核苷酸多态性而言是杂合的。如果使用两种染料,则需要将一个基因座上的每个等位基因的测量值进行组合,这是一个较为复杂的过程(参见线性测量部分)。将两种染料得到的测量值进行组合将使信号振幅加倍,并且使干扰振幅增大至大致sqrt(2),导致信号干扰比(信噪比)(SNR)加大大致sqrt(2)或者3dB。
假定没有镶嵌性以及没有参考样本的情况下的检测
假定已经通过许多试验很清楚的已知了m0,并且每次试验仅仅针对一种样本而计算出m1,并且将其与m0进行比较。N1是检测的次数,并且假定每次检测针对的是一个不同的单核苷酸多态性基因座。将阈值t设置成m0与m1的中间值,使得假阳性的可能性与假阴性的数量相等,并且如果一种样本超过了阈值,其将被标记为异常。假定s1=s2=s3=0.597,并且利用5-sigma方法(五西格玛方法),因而假阴性或者假阳性的概率为1-normcdf(5,0,1)=2.87e-7。其目标为5s1/sqrt(N1)<(m1-m0)/2,因此N1=100s1 2/(m1-m0)2=16。现在,还可以使用一种方法,在该方法中,允许假阳性的概率高于假阴性的概率,这是一种有害的情形。如果测量到阳性结果,那么试验将重新进行。因此,可以这样认为,假阴性的概率应当等于假阳性的概率的平方。对附图3进行考虑,让t=阈值,并且假定Sigma_0=Sigma_1=s。因此,(1-normcdf((t-m0)/s,0,1))2=1-normcdf((m1-t)/s,0,1)。解决了这一问题,则可以表示为t=m0+0.32(m1-m0)。因此,其结果为5s/sqrt(N1)<(m1-m0)-0.32(m1-m0)=(m1-m0)/1.47,因此N1=(52)(1.472)s2/(m1-m0)2=9。
在具有镶嵌性但没有参考样本的情况下进行的检测
假定的情况与前述的相同,不同之处在于其目的在于以97.7%的概率对镶嵌性进行检测(即,二西格玛方法)。该方法要优于羊水诊断的标准方法,在所述的标准方法中,需要提取大约20个细胞并且对这些细胞进行拍照。如果假定在这20个细胞中有一个细胞是异倍体,并且以100%的可靠度对其进行了检测,利用标准方法检测到的的具有至少一组异倍体的概率是1-0.9520=64%。如果上述细胞中有0.05%的细胞是异倍体(将其称之为样本3),那么m3=0.95m0+0.05m1并且var(m3)=(0.95s0 2+0.05s1 2)/N1。因此,std(m3)2<(m3-m0)/2=>sqrt(0.95s0 2+0.05s1 2)/sqrt(N1)<0.05(m1-m2)/4=>N1=16(0.95s2 2+0.05s1 2)/(0.052(m1-m2)2)=1001。值得注意的是,使用一西格玛统计学方法的结果,该结果仍然比常规方法达到更好的效果(即,以84.1%的可靠度进行检测),该结果可以使用相类似的方法进行表示为N1=250。
不具有镶嵌性并且使用参考样本的情况下进行的检测
尽管这一方法不是必须的,假定每次试验针对两个样本,从而将得到的m1与真实样本m2进行比较。假定N=N1=N0。计算d=m1-m0并且,假定σ1=σ0,设置阈值t=t=(m0+m1)/2,这样一来,使得假阳性的概率与假阴性的概率相等。为了得到2.87e-7的假阴性概率,必须具有下式表示的情形:(m1-m2)/2>5sqrt(s1 2/N+s2 2/N)=>N=100(s1 2+s2 2)/(m1-m2)2=32。
在具有镶嵌性并且使用参考样本的情况下进行的检测
如上所述,假定假阴性的概率是2.3%(即,二西格玛方法)。如果这些细胞中有0.05%的细胞是异倍体(将其称之为样本3),那么m3=0.95m0+0.05m1并且var(m3)=(0.95s0 2+0.05s1 2)/N1。d=m3-m2并且σd 2=(1.95s0 2+0.05s1 2)/N。得到的结果必然是std(m3)2<(m0-m2)/2=>sqrt(1.95s2 2+0.05s1 2)/sqrt(N)<0.05(m1-m2)/4=>N=16(1.95s2 2+0.05s1 2)/(0.052(m1-m2)2)=2002。再次使用一西格玛方法,其结果可以使用相类似的方式进行表示为N=500。
在一种情形中,其目的仅在于以64%的概率检测5%的镶嵌性,其中上述概率也表示了现有技术目前的状态,对该种情形进行考虑。随后,假阴性的概率应当是36%。换句话说,必须计算出x,1-normcdf(x,0,1)=36%。因此,对于二西格玛方法来说,N=4(0.36^2)(1.95s2 2+0.05s1 2)/(0.052(m1-m2)2)=65,对于一西格玛方法来说,N=33。值得注意的是,这种方法将产生一个非常高水平的假阳性,这种情况需要被进行处理(addressed),因为这种水平的假阳性在当前而言不是一种可行的选择。
同样需要注意的是,如果将N限定为384(即,每个染色体使用一个具有384个孔的Taqman培养板),并且其目的在于以97.72%的概率检测镶嵌性,那么,利用一西格玛方法对8.1%的镶嵌性进行检测将成为可能。为了以84.1%的概率检测镶嵌性(或者以15.9%的假阴性率),那么利用一西格玛方法对5.8%的镶嵌性进行检测将成为可能。为了以97.72%的置信度对19%的镶嵌性进行检测,可能需要大致70个基因座。因此,应当在每个培养板中筛选5个染色体。
在表2中提供了关于这些不同类型中的每一个的概况。在该表中还包括了由定量聚合酶链式反应以及SYBR检测生成的结果。使用如上描述的方法,并且进行简单化的假设,假定定量聚合酶链式反应在每个基因座上的检测是相同的。附图10以及附图11表示出了如上所述的样本1以及样本0的柱状图。N0=N1=47。这些样本的测量值的分布被表示为m1=27.65,s1=1.40,σm1=s1/sqrt(N1)=0.204;m0=26.64;s0=1.146,σm0=s0/sqrt(N0)=0.167。对于这些样本而言,d=1.01并且σd=0.2636。附图12表示的是男性样本以及女性样本的每个基因座上的周期时间(Ct)间的差别,全部这些基因座的差别的标准偏差是0.75。在男性样本以及女性样本中,每个基因座的每个测量值的SD值近似为0.75/sqrt(2)=0.53。
方法2:利用等位基因分型的定量技术
在这种情形中,无需假设所述的检测是定量的。相反的,假定等位基因分型是定性的,并且从该检测中无法获得有意义的定量数据。这一方法适合用于进行等位基因分型的任何检测中。附图13描述了在减数分裂过程中不同的单倍体配偶子是怎样形成的,并且可以用来描述与这种情形相关的不同种类的异倍体。最佳的运算法则取决于被检测的异倍体的类型。
考虑这样一种情形:其中异倍体的产生是由第三个片段所引起的,所述的第三个片段不含有其他两个片段中的任何一个或者两个的拷贝的部分。从附图13中可以看出,例如当p1以及p4,或者p2以及p3,都出现在婴儿细胞以及来自于另外一方亲本的一个片段中时,这种情况将会出现。这一情况非常普遍,阐明了形成异倍体的机制。一种方法是从假设h0开始,假设在细胞中存在两个片段,并且这些片段都是什么。为了进行阐述,假定h0针对的是附图13中的p3以及m4。在一个优选的实施方式中,该假设的想法来自于本文其他内容中描述的运算法则。假设h1是这样的情形:具有一个另外的片段,该片段不包括其他片段的拷贝部分。例如当p2或者m1也存在时,这种情形将会出现。可以识别出p3以及m4中的所有的纯合的基因座。可以通过寻找基因座上的杂合的基因型分型来对异倍体进行检测,其中所述的杂合的基因型分型被期望是纯合的。
假定每个基因座具有两个可能的等位基因,x以及y。将等位基因x以及等位基因y的概率分别概括为px以及py,并且px+py=1。如果h1是真实的,那么对于每个基因座i来说,p3以及m4是纯合的,那么非纯合性分型的概率是py或者px,这取决于分别在x以及y中的基因座是否是纯合的。需要注意:基于对双亲数据的认识,即p1,p2,p4以及m1,m2,m3,可以进一步的确定在每个基因座上具有非纯合型等位基因x或者y的概率。通过上述方法,利用相同数量的单核苷酸多态性可以对每种假设进行更加可靠的测量,但需要利用复杂的表示符号,因而这一扩展方法还没有被明确的研究过。本领域技术人员应该明确,怎样利用这一信息来增加假设的可靠性。
用pd来表示等位基因脱失的概率。在假设h0中,将基因座i上出现杂合基因型的概率表示为p0i,在假设h1中,将基因座i上出现杂合基因型的概率表示为p1i。
给定h0:p0i=0
给定h1:p1i=px(1-pd)or p1i=py(1-pd),这取决于在x以及y中的基因座是否是纯合的。
建立一个测量值m=1/Nh ∑i=1...Nh Ii,其中Ii是指示变量,并且当具有杂合分型时,Ii等于1,否则,Ii等于0。Nh是纯合基因座的数量。假定px=py并且对于所有的基因座来说,p0i,p1i都表示相同的两个值p0i以及p1i。在假设h0中,E(m)=p0=0并且σ2 m|h0=p0(1-p0)/Nh。在假设h1中,E(m)=p1 andσ2 m|h1=p1(1-p1)/Nh。利用五西格玛统计学方法,并且使假阳性的概率与假阴性的概率相等,则可以表示为(p1-p0)/2>5σm|h1,因此Nh=100(p0(1-p0)+p1(1-p1))/(p1-p0)2。使用二西格玛置信度替代五西格玛置信度,则可以表示为Nh=4.22(p0(1-p0)+p1(1-p1))/(p1-p0)2。
必须对足够的基因座N进行取样,并且这些基因座必须具有足够的可以利用的纯合基因座Nh-avail,这样一来可以达到至少97.7%的置信度(二西格玛方法)。定义Nh-avail=∑i=1...N Ji,其中Ji是指示变量,并且当所述的基因座是纯合型时,Ji等于1,否则,Ji等于0。所述的基因座为纯合型的概率是px 2+py 2。因此,E(Nh- avail)=N(px 2+py 2)并且σNh-avail 2=N(px 2+py 2)(1-px 2-py 2)。为了保证N是足够大的并且具有97.7%的置信度,则必须使E(Nh-avail)-2σN h-avail=Nh,其中Nh是从前述内容中得到的。
例如,如果假定pd=0.3,px=py=0.5,可以得出五西格玛置信度下的Nh=186以及N=391。相类似的,在假阴性以及假阳性中,可以在二西格玛置信度即97.7%的置信度下表示出Nh=30并且N=68。
值得注意的是,可以使用相类似的方法检测片段的缺失,其中将h0假设为存在两个已知的染色体片段,将h1假设为上述染色体片段之一发生了丢失。例如,可以寻找到所有那些本应是杂合型但实际为纯合型的基因座,在前述的等位基因脱失的影响下起作用。
同样需要注意的是,即使该检测是定性的,等位基因脱失率也可以被用来提供一种类型的定量测量,所述的定量测量针对的是存在的DNA片段的数量。
方法3:利用参考序列中已知的等位基因,进行定量的等位基因测量
在这里,假定一组常规的片段组中的等位基因或者一组预期的片段组中的等位基因是已知的。为了检查三个染色体,第一步需要对数据进行清除,假定其中的两个染色体。在本发明的一个优选的实施方式中,利用本文其他内容中描述的方法完成第一步骤数据的清除。然后,从测量数据中减去与预期的两个片段相关的信号。之后,可以在剩余的信号中寻找另外一个片段。使用匹配滤过的方法,并且用来表征另外一个片段的信号是基于确定每个片段以及这些片段的互补染色体都存在的情况下产生的。例如,针对附图13进行考虑,如果PS的结果表明了片段p2以及m1的存在,那么可以使用这里描述的技术检查另外一个染色体上的p2,p3,m1以及m4的存在。如果存在另外一个片段,则可以确保多于50%的等位基因与这些测试信号中的至少一种是共同的。值得注意的是,在这里没有进行详细描述的另外一种方法利用的是本文其他内容中描述的运算法则来进行数据的清除,假定具有异常数量的染色体,即1个、3个、4个以及5个染色体,随后再应用这里讨论的方法。在阅读过本说明书之后,这一方法的详细内容对于本领域技术人员来说应当是显而易见的。
在假设h0中,具有两个染色体,其等位基因向量是a1,a2。在假设h1中,还具有第三个染色体,其等位基因向量是a3。利用本发明中描述的方法或者其他技术对遗传数据进行清除,可以对h0中预期的两个片段的等位基因进行确定:a1=[a11...a1N]并且a2=[a21...a2N],其中每个元素aji是x或者是y。在假设h0中产生的预期信号为:s0x=[f0x(a11,a21)...fx0(a1N,a2N)],s0y=[fy(a11,a21)...fy(a1N,a2N)],其中fx,fy描述的是从等位基因组至每个等位基因测量值的图谱(mapping)。在假设h0中,可以将数据表示为dxi=s0xi+nxi,nxi~N(0,σxi 2);dyi=s0yi+nyi,nyi~N(0,σyi 2)。对数据和参考信号进行微分(differencing),从而生成一个测量值:mxi=dxi-sxi;myi=dyi-syi。完整的测量向量是m=[mx T my T]T。
现在,生成目标片段的信号——所述的片段的存在是假设性的,并且将在剩余物中对其进行寻找——根据该片段的假定的等位基因:a3=[a31...a3N]。将剩余物的信号表示为:sr=[srx T sry T]T,其中srx=[frx(a31)...frx(a3N)],并且sry=[fry(a31)...fry(a3N)],当a3i=x时,frx(a3i)=δxi,否则,该值为零,当a3i=y时,fry(a3i)=δyi,否则,该值为零。这项分析假定所使用的测量值已经被线性化(参见下面的部分),因而在基因座i处的等位基因x的一个拷贝的出现将生成数据δxi+nxi,在基因座i处的等位基因x的κx个拷贝的出现将生成数据κxδxi+nxi。然而,值得注意的是,这种假设对于本发明描述的常规方法而言不是必需的。在假设h1中,如果等位基因a3i=x,那么mxi=δxi+nxi,myi=nyi,如果a3i=y,那么mxi=nxi,myi=δyi+nyi。因此,可以得到匹配滤过器h=(1/N)R-1sr,其中R=diag([σx1 2...σxN 2σy1 2...σyN 2])。测量值m=hTd。
h0:m=(1/N)∑i=1..N srxinxi/σxi 2+sryinyi/σyi 2
h1:m=(1/N)∑i=1..N srxi(δxi+nxi)/σxi 2+sryi(δyi+nyi)/σyi 2
为了估算所需的单核苷酸多态性的数量,进行下述简单化的假设:对于所有的等位基因以及所有的基因座所作的检测都具有相类似的性质,即δxi=δyi=δ并且σxi=σyi=σfor i=1...N。之后,可以得出如下的平均偏差以及标准偏差:
h0:E(m)=m0=0;σm|h0 2=(1/N2σ4)(N/2)(σ2δ2+σ2δ2)=δ2/(Nσ2)
h1:E(m)=m1=(1/N)(N/2σ2)(δ2+δ2)=δ2/σ2;σm|h1 2=(1/N2σ4)(N)(σ2δ2)=δ2/(Nσ2)
现在计算这个检测方法的信号干扰比(信噪比)(SNR),并且将h1与h0的情况进行比较。其信号是m1-m0=δ2/σ2,并且这一测量的干扰变化是σm|h0 2+σm|h1 2=2δ2/(Nσ2)。因此,这项检测的信号干扰比(信噪比)(SNR)是(δ4/σ4)/(2δ2/(Nσ2))=Nδ2/(2σ2)。
将这一信号干扰比(信噪比)与下述情形进行比较,在所述的情形中,不进行基于等位基因分型的匹配滤过方法,仅仅简单的将每个基因座上的遗传信息进行相加。假定h=(1/N)1,其中1是N的向量,并且进行如上所述的简单化的假设,即δxi=δyi=δ并且σxi=σyi=σfor i=1...N。对于这种情形,如果m=hTd,则可以直接表示为:
h0:E(m)=m0=0;σm|h0 2=Nσ2/N2+Nσ2/N2=2σ2/N
h1:E(m)=m1=(1/N)(Nδ/2+Nδ/2)=δ;σm|h1 2=(1/N2)(Nσ2+Nσ2)=2σ2/N
因此,这项检测的信号干扰比(信噪比)是Nδ2/(4σ2)。换句话说,通过使用仅仅将预期的片段a3的等位基因的测量值进行求和的匹配滤过方法,使得所需要的单核苷酸多态性的数量减少了两个因素(factor)。它忽略了由匹配滤过的使用而导致的单核苷酸多态性的增加(gain),从而解决了在每个基因座上进行的检测的效率不同的问题。
值得注意的是,如果我们不能够准确的定义参考信号sxi以及syi,那么得到的测量信号mxi以及myi中的噪声或者干扰的SD值将会增加。当δ<<σ时,这个现象是无关紧要的,否则,它将导致假性检测概率的增加。因此,这一技术非常适合应用于检验下述假设:存在三个片段,并且其中的两个片段被假定为互为对方的准确拷贝。在这种情况中,可以使用数据清除技术从而获得可靠的sxi以及syi,其中所述的数据清除技术是以其他内容中描述的定性等位基因分型为基础的。在一种实施方式中,将方法3与方法2联合使用,其中使用了定性基因定型并且,除了等位基因脱失的定量测量之外,不能够检测到某个片段的第二个准确拷贝的存在。
现在,我们描述另外一种利用等位基因分型的定量技术。该方法包括对一个给定的等位基因的四个位点(registers)中的每一个所产生的信号的相对剂量进行比较。可以设想,在理想化的情况中,包括了一个单个的、正常的细胞,在该细胞中发生了同质扩增,(或者扩增的相对剂量是规格化的)可以发生下述四种可能的状况:(i)如果是一个杂合的等位基因,那么其四个位点处的相对强度大概为1:1:0:0,并且信号的绝对强度将与一个碱基对相符;(ii)如果是纯合的等位基因,那么其相对强度大概为1:0:0:0,并且信号的绝对强度将与两个碱基对相符;(iii)如果等位基因脱失发生在等位基因中的一个之上,那么该等位基因的相对强度大概为1:0:0:0,并且信号的相对强度将与一个碱基对相符;以及(iv)如果等位基因的脱失同时发生在两个等位基因中,那么该等位基因的相对强度大概为0:0:0:0,并且信号的绝度强度将与不含碱基对的情况相符。
然而,如果是异倍体,将观察到不同的情况。例如,如果是三倍体,且不产生等位基因的脱失(ADO),则将会发生下述三种情形中的一种:(i)如果是三重杂合的等位基因,那么其四个位点处的相对强度大概为1:1:1:0,并且信号的绝对强度将与一个碱基对相符;(ii)如果其中有两个等位基因是纯合的,那么其相对强度大概为2:1:0:0,并且信号的绝对强度将分别与两个碱基对以及一个碱基对相符;(iii)如果等位基因都是纯合的,那么该等位基因的相对强度大概为1:0:0:0,并且信号的相对强度将与三个碱基对相符。如果等位基因的脱失发生在等位基因所在的细胞具有三倍体的情况下,将会观察到预期为一个正常细胞的情形中的一种。如果是单倍体(单体性),那么其四个位点处的相对强度大概为1:0:0:0,并且信号的绝对强度将与一个碱基对相符。这种情况与在等位基因之一上发生了等位基因脱失(ADO)的正常细胞的情况相符,然而,在正常细胞的情况中,这种情形只能在少部分的等位基因中观察到。如果是其中存在两个相同染色体的单亲二体性,其四个位点处的相对强度大概是1:0:0:0,并且信号的绝对强度将与两个碱基对相符。如果是其中存在来自于一个亲本的两个不同染色体的单亲二体性,该方法将表明这一细胞是正常的,尽管利用本专利中描述的其他方法对数据进行的进一步的分析能够揭示这种单亲二体性。
在所有这些情况中,无论在正常细胞、具有异倍体的细胞、或者具有单亲二体性的细胞中,来自于一种单核苷酸多态性的数据都不足以用来确定细胞所呈现的状态。然而,如果对上述假设中的每一种的概率进行计算,并且将位于一个给定染色体中的足够数量的单核苷酸多态性的概率进行结合,一种假设将占主导地位,因而将有可能在高置信度条件下确定染色体的存在状态。
线性定量测量的方法
可以利用很多方法在一个特定基因座上对遗传物质的剂量进行线性化的测量,因而可以容易的对来自于不同等位基因的数据进行求和或者微分。我们首先讨论一种遗传方法,随后讨论一种针对特定类型的检测而设计的方法。
假定数值dxi指的是对基因座i上的等位基因x的遗传物质的剂量的非线性测量值。使用N个测量值建立一组训练数据,其中对于每个测量值而言,根据数值dxi可以估计或者得知遗传物质的剂量是βxi。对训练数据组βx,i=1...N进行选择,使其包含了在操作中可能遇到的所有不同剂量的遗传数据。使用标准的回归技术对函数进行训练,所述的函数描绘了由非线性测量值dxi至线性测量的预期值E(βxi)的情况。例如,可以使用线性回归来训练次序(order)P的多项式函数,即E(βxi)=[1 dxi dxi 2...dxi P]c,其中c是系数向量,c=[c0 c1 ...cP]T。为了对这个线性函数进行训练,建立一个N次测量的遗传物质剂量的向量βx=[βx1...βxN]T并且建立一个测量数值自乘(raise to powers)0......P的矩阵:D=[[1 dx1 dx1 2...dx1 P]T[1 dx2 dx2 2...dx2 P]T...[1 dxN dxN 2...dxN P]T]T。之后,使用最小平方适配法(least squares fit)对系数进行计算,c=(DTD)-1DTβx。
除了依靠例如适配多项式之类的常规函数,我们还可以建立特殊函数用以描述一个特定的检测,其优于常规函数。例如,我们对Taqman检测或者定量聚合酶链式反应检测进行研究。等位基因x以及某些基因座i的剂量作为时间的函数,当其达到超过某些阈值的点时,可以被描述为一个带有偏离补偿(bias offset)的指数曲线:gxi(t)=αxi+βxiexp(γxit),其中αxi是偏离补偿,γx i是指数生长速率,并且βxi依照的是遗传物质的剂量。为了根据βxi计算测量值,通过该曲线的渐进极限gxi(-∞)来计算参数αxi,随后,可以通过选取该曲线的对数来得出βxi以及γxi,从而获得log(gxi(t)-αxi)=log(βxi)+γxit,并且进行标准的线性回归。一旦我们得到αxi以及γxi的值,便可以利用另外一种方法通过时间tx来计算βxi,其中所述的时间tx指的是超过阈值gx的时间处,得出βxi=(gx-αxi)exp(-γxitx)。这将得到特定等位基因遗传数据的真实剂量的干扰测量值。
无论使用何种技术,我们都可以模拟线性测量结果βxi=κxδxi+nxi,其中κx是等位基因x的拷贝数,δxi是等位基因x以及基因座i的常数,并且nxi~N(0,σx 2),其中σx2是可以凭借经验来测量的。
方法4:对每个基因座上的遗传数据的扩增使用概率分布函数
一个特定的单核苷酸多态性的遗传物质的量将取决于该细胞中存在这种单核苷酸多态性的初始片段的数量。然而,由于扩增过程以及杂交过程的随意性,特定单核苷酸多态性的遗传物质的量将不与片段的初始数量成直接的比例关系。用qs,A,qs,G,qs, T,qs,C来表示构成等位基因的四个核酸(A,C,T,G)中的每一个的特定单核苷酸多态性的遗传物质的扩增数量。需要注意的是,这些数量可以恰好为零,这取决于扩增所使用的技术。同样需要注意的是,这些数量通常是通过来自于特定杂交探针的信号强度而测量得到的。这种对于强度的测量可以用来替代对于数量的测量,或者可以通过标准的技术在不改变本发明的性质的情况下被转化为数量的估计值。用qs来表示从特定的单核苷酸多态性的所有等位基因中得到的全部遗传物质的数量的总和:qs=qs,A+qs,G+qs,T+qs,C。用N来表示在一个细胞中含有上述单核苷酸多态性的片段的数量。N通常为2,但也可以是0,1,或者3,或者更大。对于这里讨论的任何高生产力或者中度生产力的基因定型方法而言,都可以使用qs=(A+Aθ,s)N+θs来表示得到的遗传物质的数量,其中A表示全部的扩增,其可以是凭借先验知识估算得到的,也可以是凭借经验简单测量得到的,Aθ,s是针对单核苷酸多态性而言估算值A中的错误,并且θs是导入到单核苷酸多态性的扩增、杂交以及其他反应过程中的附加的干扰。干扰项Aθ,s以及θs通常是足够大的,以至于qs将不是对于N的可靠测量值。然而,通过对染色体上的多个单核苷酸多态性进行测量,能够减轻这些干扰项带来的影响。用S来表示在一个特定染色体上测量的单核苷酸多态性的数量,所述的特定染色体是例如染色体21。可能得到针对一个特定染色体的全部单核苷酸多态性中的遗传物质的平均数量,如下:
假定A
θ,s以及θ
s是正态分布的随机变量,其平均值为0并且具有变量
以及
那么可以建立模型
其中
是正态分布的随机变量,其平均值为0并且具有变量
因此,如果对该染色体上的足够数量的单核苷酸多态性进行了测量,即
那么则可以准确的估算出N=q/A。
在另外一种实施方式中,假定所述的扩增是依照下述模型而发生的:其中从一个单核苷酸多态性中发出的信号的水平是s=a+α,其中(a+α)具有的分布与附图14中左图的分布相类似。在0点处的δ函数模拟了大约为30%的等位基因脱失率,其平均值为a,并且如果不发生等位基因的脱失,则该扩增在0至a0之间是均匀分布的。根据这一分布的平均值,可以计算出a0,a0=2.86a。现在,通过附图14中的右图来模拟α的概率密度函数。用sc来表示从基因座c中产生的信号;用n来表示片段的数量;用αi来表示根据附图14分布的随机变量,其构成来自基因座i的信号;并且用σ来表示全部{αi}的标准偏差。sc=anc+∑i=1..nc αi;平均值(sc)=anc;std(sc)=sqrt(nc)σ。如果根据附图14的右图中的分布来计算σ,则可以得出σ=0.907a2。我们可以从n=sc/(ac)中得出片段的数量,并且对于“五西格玛统计学方法”而言,我们需要std(n)<0.1so std(sc)/(ac)=0.1=>0.95a.sqrt(nc)/(ac)=0.1so c=0.952n/0.12=181。
另外一个模型对分型的置信度进行评估,并且估算必须对多少个基因座或者单核苷酸多态性进行测量才能确保得到一个给定程度的置信度,该模型将随机变量作为扩增的乘数因子(multiplier)进行合并,而不是作为额外的干扰(噪声)来源,即s=a(1+α)。取对数,log(s)=log(a)+log(1+α)。现在,建立一个新的随机变量γ=log(1+α),并且可以将这个变量假设为是呈正态分布的~N(0,σ)。在这个模型中,扩增反应可以从很小的范围到很大的范围,取决于σ,但绝对不是负数。因此α=eγ-1;并且sc=∑i=1... cna(1+αi)。对于表示符号而言,平均值(sc)与预期值E(sc)是可以相互交换的
E(sc)=acn+aE(∑i=1...cnαi)=acn+aE(∑i=1..cnαi)=acn(1+E(α))
为了得出E(α),必须得到可能的α的概率密度函数(pdf),因为α是γ的函数,是一个已知的高斯(Gaussian)概率密度函数pα(α)=pγ(γ)(dγ/dα)。因此:
并且
并且:
当σ=1时,其具有如附图15所示的形式。现在,可以通过整合该概率密度函数从而得到E(α), 这可以根据多种不同的σ数字来完成。通过这种方式得到了作为σ的函数的E(sc)或者平均值(sc)。现在,同样能够利用这个概率密度函数来计算var(sc):
var(sc)=E(sc-E(sc))2=E(∑i=1..cna(1+αi)-acn-aE(∑i=1...cnαi))2
=E(∑i=1..cnaαi-aE(∑i=1...cnαi))2
=a2E(∑i=1..cnαi-cnE(α))2
=a2E((∑i=1..cnαi)2-2cnE(α)(∑i=1..cnαi)+c2n2E(α)2)
=a2E(cnα2+cn(cn-1)αiαj-2cnE(α)(∑i=1..cnαi)+c2n2E(α)2)
=a2c2n2(E(α2)+(cn-1)E(αiαj)-2cnE(α)2+cnE(α)2)
=a2c2n2(E(α2)+(cn-1)E(αiαj)-cnE(α)2)
同样可以利用pα(α)通过多种不同的σ数字来解决,从而得到作为σ的函数的var(sc)。之后,我们可以从具有已知数量c的基因座以及已知数量n的片段的样本中选取一系列的测量值,并且从这一数据中得到std(sc)/E(sc)。这将使我们可以计算出σ的值。为了对n进行估算,E(sc)=nac(1+E(α)),因此 可以测量出,从而
当对足够多数量的独立随机变量进行求和,其分布接近于高斯形态,其中所述的独立随机变量的平均值为0,那么因此可以把sc(以及)看做是正态分布的,并且如前所述,我们可以使用五西格玛统计学方法:
其目的在于得到2normcdf(5,0,1)=2.7e-7的错误概率。由此,可以计算出基因座的数量c。
性别鉴定
在本发明所述体系的一种实施方式中,可以使用遗传数据来确定目标个体的性别。利用本发明公开的方法来确定来自于双亲的哪些染色体中的哪些片段构成了目标个体的遗传物质,之后,可以通过查找哪些性染色体是来自于父亲,从而确定目标个体的性别:X表示女性,而Y表示男性。本领域技术人员应当知晓怎样利用这一方法来确定目标个体的性别。
对假设的验证
在本发明所述体系的某些实施方式中,存在一个缺点:为了以最可能的置信度对正确的遗传状态进行预测,需要针对每个可能的状态进行假设。然而,由于遗传状态的可能性数量是异常庞大的,而计算时间是有限的,则不可能对每项假设进行检验。在这些情况中,一个另外的方法是利用假设验证的概念。这一概念包括:如果某些假设、或者某些类假设是真实的,对某些数值、某些数值组的限值进行估算,并且对期望在测量数据中出现的性质或者类型进行估算。之后,检验这些测量数值是否落入那些预期的限值范围之内,和/或对某些预期的性质或者类型进行检验,如果没有出现预期的情况,那么该运算法则能够表征那些需要进一步研究的测量值。
例如,在目标DNA中,如果一条染色体臂的一端被折断(broken off),在该种情形中,最可能的假设可以被计算为“正常”(相反的,可以是例如“异倍体”)。这是因为与遗传物质的真实状态即染色体的一端被折断相符的特定假设没有被进行检验,因为这种情况的可能性是非常低的。如果使用了假设验证的概念,那么这种运算法则将注意到相当多的数值将在预期的测量值的限值之外,其中这些数值与等位基因位于染色体上被折断的部分相符合。随机将产生表征(flag),标志着对这一情形进行进一步的调查,增加了揭开遗传物质的真实状态的可能性。
本领域技术人员应当知晓,怎样对本发明公开的方法进行修改,以使其包括上述验证技术。值得注意的是,使用本发明公开的方法被认为非常难以检测的一种异常情况是平衡易位。
含有污染DNA的方法的应用
在本发明所述体系的一种实施方式中,也可以使用本发明公开的方法对来自于目标DNA的遗传数据进行清除,其中所述的遗传数据已经明确的或者可能的被异源DNA所污染。上面概括的假设验证的观点可以被用来识别落入预期限值之外的遗传样本;在样本受到污染的情况中,期望这一验证能够提供一个表征(flag),并且能够识别出被污染了的样本。
由于可以从亲本遗传数据中得知目标DNA的大片段,并且假若污染的程度相当低并且测量了足够的单核苷酸多态性,则由于异源遗传物质而引起的欺骗性的数据将可以被识别出来。本发明公开的方法还能够进行目标基因组的重新构建,虽然其具有较低的置信度水平。假如污染的水平足够低,被计算为最具可能性的假设仍然期望与目标DNA样本中的遗传物质的真实状态相一致。
本领域技术人员应当知晓,怎样对这些方法进行优化,使其能够达到清除遗传数据的目的,其中所述的遗传数据被异源DNA所产生的欺骗性信号所污染。
付诸实施的减少实施例(Example of Reduction to Practice)
在本发明所述的体系的一种实施方式中,可以利用一组运算法则来完成如上所述的方法,这组运算法则将计算一列相关的单核苷酸多态性中的每个单核苷酸多态性间的最具可能的同一性,以及每个单核苷酸多态性分型的置信度水平。这里描述的是完成本专利中所述方法的一种可能的方式。附图16以及附图17真实的表示出了本发明公开的方法的执行过程的分解图(breakdown),输入的请求信息以及输出的格式。
附图16着重表示了输入数据(1601)及其格式和请求信息,以及输出数据(1605)及其格式。在该运算法则中输入的内容由测量的数据(1602)以及在数据库中存储的已有数据(1603)组成,其中所述的测量数据包括由使用者输入的数据,并且其中所述的已有数据因而被新收集的数据进行了更新。所述的测量数据(MD,1602)由遗传数据组成,所述的遗传数据来自于对晶胚中的目标单核苷酸多态性的测量、以及对父亲等位基因和母亲等位基因的测量,也包括上述每一个等位基因的测量的准确性或者置信度。所述的已有数据(1603)由种群频率数值(population frequency data)(FD),测量偏离数值(BD)以及交叉数值(CD)组成。
种群频率数值(FD)包含每个可获得的单核苷酸多态性的等位基因频率(每个A,C,T,G的值)。这些数据可以是先前已知的,或者是测量得出的,并且可以利用本文其他内容中描述的新的收集数据对其进行更新。
测量偏离数值(BD)捕获的是测量过程朝向某些数值的偏离。例如,假定等位基因的真实值是X=A,并且正确测量的概率是px,则测量值x的分布为:
|
A |
C |
T |
G |
概率 |
px |
PC |
PT |
PG |
不发生偏离的概率 |
px |
(1-px)/3 |
(1-px)/3 |
(1-px)/3 |
其中pX+pC+pT+pG=1。如果测量过程没有发生朝向任意数值的偏离,那么pC=pT=pG=(1-pX)/3。这一信息可以从对测量过程及其相关手段的机理的经验知识以及理论知识中得以认识。
交叉数值(CD)由遗传距离数据库以及剪切对(pairs of snips)之间的交叉概率组成,该数据是从HAPMAP(人类基因组单体型图)数据中收集得到的。
测量数值(MD)、种群频率数值(FD)、测量偏移数值(BD)、交叉数值(CD)共同组成了用于本发明公开方法(称为“双亲支持”,1604)的运算法则需要输入的必要内容。之后,该运算法则(1604)对输入数据进行处理,从而生成输出数据(1605),所述的输出数据描述了目标遗传数据给定的测量值中最有可能是“真实”的数值,以及依照双亲等位基因而判断的每个单核苷酸多态性的最有可能的来源。
附图17着重表示了运算法则(称为“亲本支持”)本身的结构,以及该运算法则是怎样运用这些输入数据中的每一个的。后台工作:为了找到最有可能的假设,必须计算出P(H|M)1707,即给出测量值的假设的概率,对于所有可能的假设H而言。
如前所述:
为了得出P(H|M)(1707),首先需要得出P(M|H)(1707),以及P(H)(1708),上述数值针对的是所有的假设H。之后就可以通过上述列出的等式计算出P(M)。假设P(H)的概率1708取决于其中假定了多少交叉反应,以及这些交叉反应(CD,1704)中的每一个发生的可能性,其中的交叉反应在前述内容中有所解释。
可以使用下述等式对P(M|H)进行计算:
P(t)1706表示的是父亲等位基因以及母亲等位基因的特定值t出现的频率,其来自于种群频率数据(FD,1703)。P(M|H&t)1705表示的是准确测量晶胚、父亲、以及母亲的等位基因值的概率,其中假定一个特定的“真实”值t。测量数值、使用者输入的准确性(MD,1701)、以及测量偏离数据库(BD,1702)是计算P(M|H&t)1705所必需的输入内容。
下面将立即给出对于上述方法的更加详细的说明。从单核苷酸多态性R={r1,...,rk}(k个单核苷酸多态性组成的一组)、以及相应测量的父亲与晶胚的同一性M=(e1,e2,p1,p2,m1,m2)开始,对于k个单核苷酸多态性而言,使用id’s s1,...,sk对每个单核苷酸多态性进行识别,其中:
e1=(e11,e12,...,e1k)表示的是针对所有的单核苷酸多态性,在晶胚的染色体之一(它们不是全部必须来自于相同的亲本染色体的)中进行的测量
e2=(e21,e22,...,e2k)表示的是在晶胚的另外一个染色体中进行的测量
p1=(p11,p12,...,p1k)表示的是在父亲的第一染色体中进行的测量(全部来自于相同的染色体)
p2=(p21,p22,...,p2k)表示的是在父亲的第二染色体中进行的测量(全部来自于相同的染色体)
m1=(m11,m12,...,m1k)表示的是在母亲的第一染色体中进行的测量(全部来自于相同的染色体)
m2=(m21,m22,...,m2k)表示的是在母亲的第二染色体中进行的测量(全部来自于相同的染色体)
还可以表示为M={M1,...,Mk},其中Mi=(e1i,e2i,p1i,p2i)
本方法的目的在于确定“真实的”晶胚数值T=(E1,E2),即,给定测量值M的最有可能的情形,其中:
E1=(E11,E12,...,E1k)表示的是在晶胚的第一染色体中进行的测量,其中所述的染色体与父亲的染色体相符合,E1i∈{p1i,p2i}
E2=(E21,E22,...,E2k)表示的是在晶胚的第二染色体中进行的测量,其中所述的染色体与母亲的染色体相符合,E2i∈{m1i,m2i}
还可以表示为T={T1,...,Tk},其中Ti=(E1i,E2i)。
有效的,将双亲的染色体值(p1,p2,m1,m2)作为支持,用于进行对(e1,e2)的测量值进行检查、验证以及修正,因此,该运算法则得名为“双亲支持运算法则”。
为了达到这一目的,建立关于晶胚值的来源的所有可能的假设,并且从其中选出最具可能性的一个,计算测量值M。设想的空间是SH={H1,......,Hq}={所有的假设组},其中每个假设具有Hj=(Hj 1,......Hj k)的形式,其中Hj 1是对于单核苷酸多态性i的“小型”假设,具有Hj 1=(pi *,mi *)的形式,其中pi *∈{p1 i,p2i}并且mi *∈{m1j,m2i}。存在4种不同的“小型(mini)”假设Hj 1,具体是:
Hj 11:(e1i,e2i)={(p1i,m1i)或者(m1i,p1i)}
Hj 12:(e1i,e2i)={(p1i,m2i)或者(m2i,p1i)}
Hj 13:(e1i,e2i)={(p2i,m1i)或者(m1i,p2i)}
Hj 14:(e1i,e2i)={(p2i,m2i)或者(m2i,p2i)}
从理论上说,SH可以具有q=4k个不同的成员以供选择,尽管在后来的步骤中这一空间将被父亲染色体以及母亲染色体的交叉反应的最大数量进行限制。
选择最具可能性的假设H*,如:
从中选择出的最具可能性的假设H*=argmaxH∈SHP(H|M)
对于特定的H而言:
因此,从每种假设中可以得出:
(1)P(M/H)是测量值M给出特定的假设H的概率
(2)P(H)是特定的假设H的概率
(3)P(M)是测量值M的概率
在计算过所有的H的P(H|M)之后,选择出具有最高概率的假设。
计算P(M|H)
由于对于每个单核苷酸多态性的测量都是独立的,M=(M1,...Mk)并且对于所有的k个单核苷酸多态性的特定假设H=(H1,...Hk),那么:
P(M|H)=P(M1|H1)*......*P(Mk|Hk)
对于特定的单核苷酸多态性r来说,计算P(Mr|Hr)。因为Ω={A,C,T,G}X{A,C,T,G}X={A,C,T,G}X{A,C,T,G},通过贝叶斯公式计算出的“真实的”双亲值(P1r,P2r,M1r,M2r)的所有可能的值的空间为:
计算P(Mr/Hr&(P1r,P2r,M1r,M2r)=t)
Mr=(e1r,e2r,p1r,p2r,m1e,m2r)表示的是在这个单核苷酸多态性中得到的测量值。
T=(E1r,E2r,P1r,P2r,M1r,M2r)表示的是假定的“真实”值,通过假设,t=(P1r,P2r,M1r,M2r)并且(E1r,E2r)是由T而得来的。(E1r是P1r和P2r中的一种,E2r是M1r和M2r中的一种)
P(Mr=(e1r,e2r,p1r,p2r,m1r,m2r)/T=(E1r,E2r,Pr1,P2r,M1r,M2r))=
P(e1r/E1r)*P(e2r/E2r)*P(p1r/Pr1)*P(p2r/P2r)*P(m1r/M1r)*P(m2r/M2r)
得出:
peri=P(准确测量的晶胚的值i,在单核苷酸多态性r中)
ppri=P(准确测量的父亲的值i,在单核苷酸多态性r中)
pmri=P(准确测量的母亲的值i,在单核苷酸多态性r中)
如果不存在测量偏离,则其中p(e1r,E1r,r)=1/3,否则,该值可以通过试验数据进行确定,例如通过来自于Hapmap(人类基因组单体型图)计划中的数据进行确定。
计算P((P1r,P2r,M1r,M2r)=t)
让t=(t1,t2,t3,t4):
P((P1r,P2r,M1r,M2r)=(t1,t2,t3,t4))=P(P1r=t1)*P(P2r=t2)*P(M1r=t3)*P(M2r=t4)
假定具有n个样本(P1,P2,M1,M2),假定所有的父亲的值以及母亲的值都是独立的,并且对于{A,C,T,G}中的ti而言,t=(t1,t2,t3,t4)
当t1=A时,为了得到一个特定的p1A=P(P1=t1),假定当缺少任意数据时,这一概率可以是0至1间的任意值,因此它被指定为一个值U(0,1)。由于获得了数据,该值被新值进行了更新,并且这一参数的分布成为一种β分布。假定在P1的n种观察结果中,有h个值的P1=A,并且w=(事件P1=A)并且D=(给定的数据)。在前面的部分描述了p(w|数据)呈一种β分布,其中α=h+1,β=n-h+1(参见等式(8))。预期值以及变量X~B(α,β)的分布如下:
因此,下述参数的后验平均值为p1rA=P(P1r=A|Data)=(h+1)/(n+2)
类似的,p1rB=(#(p1r=B)+1)/(n+2),...m2rG=(#(m2r=G)+1)/(n+2),等等。因此,能够计算出所有的值p1rA,...,m2rG,并且:
计算P(H)
假设H的概率H=(H1,...,Hk)取决于染色体交叉的量,其中所述的概率H具有Hi=(pi *,mi *)。例如,
当P(交叉)=0时,则P(H)=1/4,并且如果p*在{(p11,p21,...ps1),(p12,p22,...,ps2)},m*在{(m11,m21,...,ms1),(m12,m22,...,ms2)}中,则H=(p*,m*),否则,H等于零。
当P(交叉)>0时,将每个单核苷酸多态性之间的交叉的概率结合进来是很重要的。
假设H是由对父亲染色体以及母亲染色体中的每个单核苷酸多态性进行的假设而组成的pi *∈{p1i,p2i}并且mi *∈{m1i,m2i},即H=(Hp,Hm),其中Hp=(p1 *,...pk *),并且Hm=(m1 *,...mk *),它们都是各自独立的。
P(H)=P(Hp)*P(Hm)。假定单核苷酸多态性被增加的位置所排序(ordered),
其中PCi=P(交叉(ri-1,ri)),即,在单核苷酸多态性ri-1,ri之间的某处发生交叉反应的概率,并且,如果pi *,pi-1 *全部来自于p1或者p2,则Ii=1,否则,Ii=0。
计算P(交叉(a,b))
给定单核苷酸多态性a,b,在碱基位置la,lb处(给定的碱基),交叉反应的概率大约为:
P(la,lb)=0.5(1-exp(-2G(la,lb)))
其中G(la,lb)等于位置la,lb之间的遗传距离,以Morgans为单位。不存在关于G的精确的闭合形状函数,但是其可以被松散的估算为G(la,lb)=|la-lb|*1e-8。通过利用HapMap(人类基因组单体型图)数据库中的碱基位置si以及距离G(si,si+1),可以使用更好的近似值,其中i的范围包括了所有的位置。特别的, 因此,可以将其用在交叉概率中。
计算P(M)
一旦已知了P(M|H),则可以为SH中所有不同的H计算出P(H),
计算出具有最大概率的假设的一种更加有利的方法
给定计算时间的限制,由于单核苷酸多态性的数量增加,上述方法的复杂性成指数缩放(exponential scaling),在某些情况中,必须利用更加有利的方法来确定具有最大概率的假设,并且因此进行相关单核苷酸多态性的分型。下面描述了一种完成上述方法的更加快速的方式:
从前述内容可知,P(H|M)=P(M|H)*P(H)/P(M),argmaxHP(H|M)=argmaxH并且P(M|H)*P(H)=argmaxHF(M,H),其目的在于计算H,取F(M,N)的最大值。
假定M(s,k)等于在剪切(snips)s至剪切k间的测量值,H(s,k)等于对剪切s至剪切k间的假设,简略表示形式为M(k,k)=Mk,H(k,k)=Hk=针对剪切k的测量值以及假设。如前所示:
并且同样的
其中
并且PC(Hi-1,Hi)=在Hi-1,Hi之间发生交叉反应的概率因而,最后,对于剪切n而言:
F(M,H)=P(M|H)*P(H)=P(M(1,n)|H(1,n))*P(H(1,n))
=P(M(1,n-1)|H(1,n-1))*P(H(1,n-1))*P(Mn|Hn)*PF(Hn-1,Hn)
从而:F(M,H)=F(M(1,n),H(1,n))=F(M(1,n-1),H(1,n-1))*P(Mn|Hn)*PF(Hn-1,Hn)
因此,可以将对n个剪切(snips)的计算简化为对n-1个剪切的计算。
对于关于n个剪切的假设H=(H1,...Hn)而言:
其中
概括得出:
其中,可以递归的计算出G;当i=2,......,n时,
并且G(M(1,1),H1)=0.25*P(M1|H1)
通过进行下述运算符法则的运算,可以找出最佳的假设:
步骤1:当I=1时,为H1建立四种假设,为这四种假设建立G(M1|H1),并且记录G1,G2,G3,G4
步骤2:当I=2时,为H2建立四种假设,使用上面的规则建立G(M(1,2)|H2):
并且记录这四个新的Gn。
用I=k来重复步骤2,其中k
i=k
i-1+1直至k=n:生成四种关于H
k的假设,得出G(M
(1,k)|H
k),
并且记录这四个G
n。
由于在任何时间都只能记录四种假设,操作的数量是恒定的,因而该运算法则是线性的。
计算P(M):P(H|M)=P(M|H)*P(H)/P(M)=F(M,H)/P(M))
如上所述:
其中
可以使用递归式来解答W(M,H):
因此:
并且
该运算法则与前述描述的情形相类似,其中i=2∶n,并且在每个步骤中,都会建立一组新的W(i),直至最后的步骤产生了最优化的W。
从d1,d2,h,pd1,pd2,ph中计算p1,p2,pp1,pp2
为了进行解释,这一部分将着重涉及父亲的二倍体数据以及单倍体数据,很重要的一点是,需要注意,相同的运算法则也可以应用于母亲。使用:
○d1,d2——二倍体测量值中的等位基因分型
○h——单倍体测量值中的等位基因分型
○pd1,pd2——在每个二倍体测量值中的正确等位基因分型的概率
○ph——在单倍体测量值中的正确等位基因分型的概率
这些数据应当被测绘(map)至这里公开的运算法则的下述输入参数中:
○p1——单倍体细胞以及二倍体细胞之一所对应的等位基因
○p2——余下的二倍体细胞所对应的等位基因
○pp1,pp2——正确等位基因分型的概率
由于h对应于d1,那么为了计算p1的值,必须要用到h以及d1。之后,p2将自动对应于d2。类似的,如果h对应于d2,那么为了计算p1的值,必须要用到h以及d2。之后,p2将自动对应于d1。
在这里使用术语“对应于”,是因为它既可以具有“相同”的含义,也可以具有“以更高的概率来自于”的含义,其含义取决于不同的测量结果以及种群概率。
上述运算法则的目的在于,计算出隐藏在原始测量值h,d1,d2,ph,pd1,pd2以及种群频率背后的“真实”等位基因值的概率。
基础的运算法则步骤如下:
(i)根据h,d1,d2,ph,pd1,pd2值以及种群频率数据,确定h对应于d1还是对应于d2。
(ii)将等位基因分型指定为p1以及p2;基于步骤(1),而计算出概率pp1以及pp2。
将h指定为d1或者d2
建立两种假设:
H1:h对应于d1(h来自于d1)
H2:h对应于d2(h来自于d2)
其任务在于计算出这两种假设得到测量值M的概率:
P(H1/M(h,d1,d2,ph,pd1,pd2))并且P(H2/M(h,d1,d2,ph,pd1,pd2))
(为了对内容进行简化,此后将使用P(H1/M)以及P(H2/M)来表示)
为了计算出这些概率,应用贝叶斯定理:
其中,P(M)=P(M/H1)*P(H1)+P(M/H2)*P(H2)。由于假设H1与假设H2具有相等的可能性,P(H1)=P(H2)=0.5,因此:
并且
为了计算P(M/H1)以及P(M/H2),必须考虑到二倍体结果d1以及d2的所有可能值的组,Ω={AA,AC,...,GG},即,A、C、T、G的任意组合,即所谓的潜在状态。当针对潜在状态进行假设时(即,在值d1以及d2中,伴有基于假设H1或者H2的假定值h),可以分别生成关于h,d1以及d2的“真实值”H,D1以及D2的所有可能的组合(状态S={s1,s2,...,s16}),如下表所示:
由于“真实值”H,D1以及D2是未知的,并且仅仅已知原始测量结果h,d1,d2,ph,pd1,pd2,则必须通过下述方式来计算整个组的Ω的P(M/H1)以及P(M/H2):
如果,为了进行计算,假定d1=d2,并且pd1以及pd2都是自变量,则可以表示为:
计算上述等式中最后一个求和项中的前三项:P(M(x)/X),其中x在{h,d1,d2}中。
为了计算正确的等位基因分型概率(命中“真实的等位基因值”),需要基于结果的测量值x来算出等位基因的真实值X。如果测量值x与真实值X相等,则该概率为px(正确测量的概率)。如果x与X是不相等的,则该概率为(1-px)/3。例如,在X=C,并且测量值x=A的情况下,计算“真实值”C的概率。得到A的概率是px。得到C、T或者G的概率是(1-px)。因而,命中C的概率是(1-px)/3,因为可以假定C、T以及G具有相同的可能性。
如果将指示变量Ix也包括在该计算当中,其中,当x=X时,Ix=1,并且当x≠X时,Ix=0,得到的概率如下:
P(M(x)/X)=I{x=X}*px+(1-I{x=X})*(1/3)*(1-px),x在{h,d1,d2}中选取
现在,考虑在P(M|H1)中的最后两项。P(D1)以及P(D2)表示的是等位基因A、C、T以及G的种群频率,所述的种群频率可以从现有技术中获知。
针对特定的状态2,计算上述表达式,给出特定的测量值M(h=A,d1=G,d2=C):
P(M(h)|H)*P(M(d1)|D1)*P(M(d2)|D2)*P(D1)*P(D2)=
=P(M(h)=A|H=A)*P(M(d1)=G|D1=A)*P(M(d2)=C|D2=C)*P(D1=A)*P(D2=C)=
=ph*((1-pd1)/3)*pd2*f(D1=A)*f(D2=C)
类似的,针对剩余的15种状态及其总和Ω组,计算(1)给出的特定测量值(在这种情况中,M(h=A,d1=G,d2=C))。
现在,已经计算出P(M/H1)以及P(M/H2)。最后,按照如下描述的方法计算P(H1/M)以及P(H1/M):
指定(Assigning)等位基因分型以及相应的概率
现在建立下述四种不同的假设:
Hp2A:p2的“真实值”是A
Hp2C:p2的“真实值”是C
Hp2T:p2的“真实值”是T
Hp2G:p2的“真实值”是G
并且计算P(Hp2A/M),P(Hp2C/M),P(Hp2T/M),P(Hp2G/M)。其中最大的值决定了特定的等位基因分型以及相应的概率。
由于p2的来源是未知的(其来自于d1的概率是P(H2/M),其来自于d2的概率是P(H1/M)),因而必须考虑到p2等位基因来源于d1或者来源于d2的两种情况。针对假设HA,应用贝叶斯定理,得出:
P(Hp2A|M)=P(Hp2A|M,H1)*P(H1|M)+P(Hp2A|M,H2)*P(H2|M)
在步骤1中已经计算出了P(H1/M)以及P(H1/M)。通过贝叶斯定理:
由于在假设H1中意味着p2来自于d2:
如前所述:P(H1,M/Hp2A)=P(M(d2)/D2=A)=I{d2=D2}*pd2+(1-I{d2 =D2})*(1/3)*(1-pd2)。
P(Hp2A)=P(D2=A)=fd2(A),其中fd2(A)是从种群频率数据中得到的。
P(H1,M)=P(H1,M/Hp2A)*P(Hp2A)+P(H1,M/Hp2C)*P(Hp2C)+P(H1,M/Hp2T)*P(Hp2T)+P(H1,M/Hp2G)*P(Hp2G)。
类似的,计算P(Hp2A&H2/M)。
P(Hp2A/M)=P(Hp2A&H1/M)+P(Hp2A&H2/M),因此,可以计算出p2等于A的概率。针对C、T、以及G重复进行上述计算。得到的最大值将是p2等位基因分型的答案,及其相应的概率。
指定(Assigning)等位基因分型为p1(相应于单倍体细胞以及二倍体细胞之一的等位基因)
如前所述,我们建立四种不同的假设:
Hp1A:p1的“真实值”是A
Hp1C:p1的“真实值”是C
Hp1T:p1的“真实值”是T
Hp1G:p1的“真实值”是G
并且计算P(Hp1A/M),P(Hp1C/M),P(Hp1T/M),P(Hp1G/M)
下面对假设Hp1A进行详细的细节描述。在“真实情况”的情况下,仅当单倍体细胞以及相应的二倍体细胞等于A时,p1将等于A。因此,为了计算p1以及pp1,必须考虑到单倍体细胞与相应的二倍体细胞相同的情况。因此,在假设Hp1A中,p1的“真实值”是A并且得到HhdA:单倍体细胞以及相应的二倍体细胞的“真实值”是A。
由于h的来源是未知的(其来自于d1的概率是P(H1/M),其来自于d2的概率是P(H2/M)),因而必须考虑到h等位基因来源于d1或者来源于d2的两种情况,并且在确定p1时应用这两种情况。这意味着,利用贝叶斯定理:
P(HhdA|M)=P(HhdA|M,H1)*P(H1|M)+P(HhdA|M,H2)*P(H2|M)
如前所述,可以从先前的计算中得知P(H1/M)以及P(H2/M)。
P(H1,M/HhdA)=P(M(h)/H=A)*P(M(d1)/D1=A)=
=[I{h=H}*ph+(1-I{h=H})*(1/3)*(1-ph)]*[I{d1=D1}*pd1+(1-I{d1=D1})*(1/3)*(1-pd1)],
因为在假设H1中意味着p1来自于d1。P(HhdA)=P(h=A)*P(D1=A)=fh(A)*fd1(A),其中fh(A)以及fd2(A)是从种群频率数值中得到的。P(H1,M)=
P(H1,M/HhdA)*P(HhdA)+P(H1,M/HhdC)*P(HhdC)+P(H1,M/HhdT)*P(HhdT)+P(H1,M/HhdG)*P(HhdG)
类似的,计算出P(HhdA&H2/M)。
P(HhdA/M)=P(HhdA&H1/M)+P(HhdA&H2/M),并且现在,我们计算出了p1等于A的概率。针对C、T、以及G重复进行上述计算。得到的最大值将是p1等位基因分型的答案,及其相应的概率。
输入实施例
下面给出了两个输入实施例。第一个实施例是一组具有低倾向性共分离的单核苷酸多态性,其中,单核苷酸多态性分散在染色体中,其输入数据如表3中所示。第二个实施例是一组具有高倾向性共分离的单核苷酸多态性,其中,单核苷酸多态性聚集在染色体中,其输入数据如表4中所示。这两组数据都包括个体的测量的单核苷酸多态性数据,个体双亲的单核苷酸多态性数据,以及相应的置信度数值。值得注意的是,这一数据是从真实的人体中测量得到的真实的数据。每一横行表示的是一个特定单核苷酸多态性位点上的测量值。每一列中包含的数据的含义在列头处有所指示。在列头中使用的缩略语的主要意思如下:
○家族_id=每个人具有的独特的id(由于抄写原因而被包括在内)
○单核苷酸多态性_id=单核苷酸多态性的识别号码
○e1,e2=晶胚的单核苷酸多态性值
○p1,p2=父亲的单核苷酸多态性值
○m1,m2=母亲的单核苷酸多态性值
○pe1,pe2=对e1,e2的测量准确性
○pp1,pp2=对p1,p2的测量准确性
○pm1,pm2=对m1,m2的测量准确性
输出实施例
两个输出实施例的数据在表5以及表6中有所表示,并且这些输出数据分别与表3以及表4中列出的输入数据是相互对应的。两个表格中都给出了个体测量的单核苷酸多态性数据,个体双亲的单核苷酸多态性数据,个体单核苷酸多态性数据中最具可能性的真实值,以及相应的置信度。每一横行表示的是一个特定单核苷酸多态性位点上的测量值。每一列中包含的数据的含义在列头处有所指示。在列头中使用的缩略语的主要意思如下:
○单核苷酸多态性_id=单核苷酸多态性的识别号码
○真实_值(true_value)=设想的e1,e2的核苷值
○真实_假设(true_hyp)=对e1,e2的来源的假设
○ee=针对e1,e2的测定的单核苷酸多态性核苷值
○pp=针对p1,p2的测定的单核苷酸多态性核苷值
○mm=针对m1,m2的测定的单核苷酸多态性核苷值
○假设可能性(HypProb)=最终假设的概率。对于该项输出数据来说,仅存在一个数值,但由于excel列表的结构,导致该值重复出现在所有的行中。
值得注意的是,可以人工执行这项运算法则,或者通过计算机执行。表3以及表4表示的是由计算机完成该方法的版本中使用的输入数据的实施例。表5表示的是表3中所示的输入数据的相应输出数据。表6表示的是表4中所示的输入数据的相应输出数据。
模拟运算法则
下面描述的是第二个模拟方法,其作用在于确保本发明体系的完整性,并且对上述运算法则在各种不同的情况中的实际效率进行评价。为了达到上述目的,需要运行1000次全体系的模拟。这包括随机建立双亲遗传数据,利用计算机数据分析方法(in silico)仿效减数分裂从而生成胚胎数据,模拟对胚胎数据的不完全测量,随后运行本发明公开的方法对模拟的胚胎测量数据进行清除,之后,将“清除”过的数据与“真实”数据进行比较。下面的内容中对该模拟方法进行了更加详细的解释,并且在附图18中描绘出了该方法中所有事件的流程图。对同一理论的两种不同的执行方法进行了检验。将在下面的内容中给出更加全面的解释。
针对DH以及PS的模拟运算法则及其结果
对于这两个运算法则而言,其初始的输入变量为:
(i)需要进行检测的单核苷酸多态性列表
(ii)母亲染色体以及父亲染色体的种群频率
(iii)单倍体测量中的正确等位基因分型的概率(ph,pe)以及无序的二倍体测量中的正确等位基因分型的概率(pd)。
应当根据从来自于关于相关单核苷酸多态性的经验数据(种群频率)、以及来自于计量仪表的显示(ph,pd,pe)中得到的结果来确定上述数值。针对几个类型进行这种模拟,这些类型例如是最具可能性的(告知的)、普通的(未被告知的)、以及非常不可能的(极限情况)。
一旦确定了上述的统计学参数,对于所有的模拟来说,给出特定单核苷酸多态性的交叉概率都是相同的,并且将在给出剪切位点数据库(SNIPLOC_NAME_MAT)以及遗传距离数据库(HAPLOC_NAME_MAT)之前得出。
[交叉概率,剪切]=
GetCrossProb(snips,SNIPLOC_NAME_MAT,parameters,HAPLOC_NAME_MAT)
初步模拟循环
初步模拟循环是用来证明将被用于整个模拟过程中的遗传数据是真实的。将步骤1至5重复进行10000次。值得注意的是,这项模拟可以应用于双亲中的任何一方或者双方;其步骤是相同的。在这种情况下,为了进行阐述,将使本次模拟应用于父亲的情况中,并且,在附图18中所示的相关内容也将包括附图18的圆括号内所表示的相应的母亲的值输入入口。
步骤1:建立初始的父亲二倍体细胞(P1,P2)
[P1,P2]=建立初始的染色体(剪切;父亲染色体的种群频率);1801(1802)
建立初始的父亲细胞取决于父亲细胞中的每个单核苷酸多态性的种群频率。
步骤2:针对DHAlgo建立单倍体数据以及无序的二倍体数据
模拟父亲染色体的交叉1803,从而得到两组交叉的染色体:P1C1,P2C1以及P1C2,P2C2;1804(1805)。在交叉反应1806发生之后,从父亲的等位基因中选取一个(从第一组中)作为P1情况中的单倍体等位基因HP1807(1808)(因为选择哪个等位基因都没有区别),并且在二倍体等位基因中打乱顺序,从而得到(D1P,D2P)1807(1808)。
HP=选取的(P1C1,P2C1);
[D1P,D2P]=杂乱的(P1,P2)。
步骤3:向初始数据组中导入误差,从而模拟测量值
基于给定的正确测量的概率(ph-单倍体测量值,pd-二倍体测量值),向这些测量值中导入误差,从而得到模拟的双亲测量数据1811(1812)。
hp=导入的误差(HP,ph);
d1p=导入的误差(D1P,pd);
d2p=导入的误差(D2P,pd)。
步骤4:应用DH运算法则来计算(p1,p2)(pp1,pp2)
DH运算法则从单倍体细胞中选取等位基因,并且从二倍体细胞中选取无序的等位基因,并还原出引起上述情形的最具可能性的有序的二倍体等位基因。DH运算法则试图重建(P1,P2),并且还原对于父亲细胞的估算误差(pp1,pp2)。为了进行比较,同样进行经验性的运算法则,该运算法则进行简单的等位基因配对。这样做的目的在于比较本发明公开的运算法则究竟有多么优越,这种比较针对的是简单的经验性运算法则。
[p1,p2,pp1,pp2]=DHAlgo(hp,d1p,d2p,ph,pd,snips,popfreqlistPP,′DH′);[p1s,p2s,pp1s,pp2s]=DHAlgo(hp,d1p,d2p,ph,pd,snips,popfreqlistPP,′ST′);
步骤5:收集操作所需的统计学参数
比较(P1,P2)从而得出(p1,p2)。
[P1cmp(:,i),P2cmp(:,i),P1prob(:,i),P2prob(:,i),P1mn(i),P2mn(i)]=DHSimValidate(P1,P2,p1,p2,pp1,pp2);
需要注意的是:(P1Si,P2Si,P1Pi,P2Pi,P1Ai,P2Ai)=(I{P1=p1},I{P2=p2},pp1,pp2,p1acc,p2acc),其中I{P1=p1}表示的是二元指示阵列(binary indicator array),用来评估DH运算法则针对所有的单核苷酸多态性的运算准确性。类似的,对于I{P2=p2}来说,pp1,pp2是从该运算法则中得到正确的等位基因分型的概率,并且p1acc=mean(I{P1=p1}),即,对于p1而言,此次操作的平均准确率与p2acc(对于p2的操作的准确率)相类似。
初步的模拟结果
进行一万次的模拟,来评估该项运算法则的准确性。DH准确性P1=mean(P1Ai),DH准确性P2=mean(P2Ai),这表明了DH运算法则对于P1,P2的总体的准确性。对于个体的单核苷酸多态性而言,对于每个单核苷酸多态性的准确性:单核苷酸多态性准确性P1=mean(P1Si)应当与正确测量该单核苷酸多态性的估算概率的平均值相适合,即与单核苷酸多态性概率P1(SNPProb.P1)=mean(P1Si)相适应,即,如果该运算法则的运算是正确的,则单核苷酸多态性准确性P1(SNPAcc.P1)应当与单核苷酸多态性概率P1紧密相关。上述两者之间的关系是通过它们的相关性反映出来的。
上述10000个模拟循环是在不同的设置类型下运行的:
(1)利用更加真实的已有的基因型数据而得到的潜在的种群频率,以及均一的种群频率,其中在每个单核苷酸多态性中,A,C,T,G具有相同的概率。
(2)对于单倍体的测量(ph)以及无序的二倍体的测量(pd)的几种测量准确性的组合。给出各种假设;假设这些测量都是非常准确的(0.95,0.95),不太准确(0.75,0.75)以及不准确或者随机(0.25,0.25),以及不均衡的组合(0.9,0.5),(0.5,0.9)。最接近真实的情况应当是准确性大约为0.6至0.8的情况。
(3)模拟是针对所有的情况同时使用DH运算法则以及简单匹配的ST运算法则,这是为了评价本发明公开的运算法则的表现。
上述所有操作的结果在表7当中有所概况。
在这些模拟中,本发明公开的运算法则的表现要优于现有的经验性的运算法则,特别是在不均一的种群频率的真实情况中、以及在正确测量的不均衡概率或者降低的概率的真实情况中。同样被证明的是,在这些情况中,我们对于个体单核苷酸多态性的运算法则的准确性的估算是很好的,因为在平均比为1的情况下,对于正确的等位基因分型的估算的准确率与模拟的平均准确率的相关性是大约99%。
在最真实的情况中,当数据种群频率以及(ph,pd)=(0.6,0.8)时,在执行方法1中,能够准确的重新获得(P1,P2)的单核苷酸多态性的平均百分比是(0.852,0.816),在执行方法2中,能够准确的重新获得(P1,P2)的单核苷酸多态性的平均百分比是(0.601,0.673)。
需要注意的是,在表7以及表8中,以“数据”开头的行中使用的种群频率来自于经验结论,而以“均一”开头的行中使用的种群频率是假定的均一的种群。
很重要的一点是,需要注意,在表7以及表8中的准确性被定义为具有正确的单核苷酸多态性分型以及被正确的识别了染色体的来源的单核苷酸多态性的平均百分比。同样重要的是,需要注意,这些模拟反应反映的是该运算法则的两种可能的执行方法。还可能存在有其他执行该运算法则并可能取得更好效果的方法。这一模拟仅仅意在证明该方法是可以付诸实施的。
完全模拟循环
将步骤1-8重复进行10000次。这项模拟是为了检验利用完全公开的方法使用从相关个体中测量得到的遗传数据对目标个体的测量得到的遗传数据进行的清除,在这个案例中,所述的相关个体是双亲。
步骤1:建立初始的双亲二倍体细胞(P1,P2),(M1,M2)
[P1,P2]=建立初始的染色体(剪切,父亲染色体的种群频率);(1801)
[M1,M2]=建立初始的染色体(剪切,母亲染色体的种群频率);(1802)
建立初始的双亲细胞取决于父亲细胞以及母亲细胞中的每个单核苷酸多态性的种群频率。
步骤2:交叉的双亲细胞(P1C,P2C),(M1C,M2C)(1803)
建立两组交叉的父亲细胞:首先使用DH运算法则得到(P1C1,P2C1),接下来使用PS运算法则来得到(P1C2,P2C2)(1804)
建立两组交叉的母亲细胞:首先使用DH运算法则得到(M1C1,M2C1),接下来使用PS运算法则来得到(M1C2,M2C2)(1805)
[P1C1,P2C1]=Cross(P1,P2,snips,fullprob);
[P1C2,P2C2]=Cross(P1,P2,snips,fullprob);
[M1C1,M2C1]=Cross(M1,M2,snips,fullprob);
[M1C2,M2C2]=Cross(M1,M2,snips,fullprob)
步骤3:针对DH运算法则而建立单倍体细胞以及无序的二倍体细胞(1806)
为单倍体细胞HP而选取一组父亲细胞(1804,第一组),并且在二倍体细胞中打乱顺序从而得到(D1P,D2P)(1807)。为母亲细胞进行同样的操作(1805,第一组),从而得到MH(D1M,D2M)(1808)。
HP=选取(P1C1,P2C1);
HM=选取(M1C1,M2C1);
[D1P,D2P]=杂乱的(P1,P2);
[D1M,D2M]=杂乱的(M1,M2);
步骤4:进行二倍体晶胚分型(1809)
为晶胚细胞选取父亲细胞之一(1804,第二组)以及母亲细胞之一(1805,第二组)。为了进行测量而打乱顺序。
E1=选取(P1C2,P2C2);
E2=选取(M1C2,M2C2);
[E1J,E2J]=杂乱的(E1,E2);(1810)
步骤5:向测量值中导入误差(1811,1812,1813)
为了给出测量值误差(ph-单倍体细胞,pd-无序的二倍体细胞,pe-晶胚细胞),向测量值中导入误差。
hp=导入的误差(HP,ph);(1811)
d1p=导入的误差(D1P,pd);(1811)
d2p=导入的误差(D2P,pd);(1811)
hm=导入的误差(HM,ph);(1812)
d1m=导入的误差(D1M,pd);(1812)
d2m=导入的误差(D2M,pd);(1812)
e1=导入的误差(E1J,pe1);(1813)
e2=导入的误差(E2J,pe2);(1813)
步骤6:运用DH运算法则计算出(p1,p2),(m1,m2),(pp1,pp2),(pm1,pm2)
DH运算法则选取一个单倍体细胞以及一个无序的二倍体细胞,并且还原出引起上述情形的最具可能性的有序的二倍体等位基因。DH运算法则试图针对父亲染色体重建(P1C1,P2C1),针对母亲染色体重建(M1C1,M2C1),并且还原对于父亲细胞的估算误差(pp1,pp2)以及对于母亲细胞的误差估算(pm1,pm2)。
[p1,p2,pp1,pp2]=DH运算法则(hp,d1p,d2p,snips,popfreqlistPP);(1814)
[m1,m2,pm1,pm2]=DH运算法则(hm,d1m,d2m,snips,popfreqlistMM);(1815)
步骤7:运用PS运算法则计算出(DE1,DE2)(1816)
PS运算法则选取重建的双亲细胞(p1,p2,m1,m2)以及无序测量的晶胚细胞(e1,e2),并还原出最具可能性的有序的真实的晶胚细胞(DE1,DE2)。PS运算法则试图重建(E1,E2)。
[DE1,DE2,alldata]=PS运算法则(snips,e1,e2,p1,p2,m1,m2,pe,pp1,pp2,pm1,pm2,parameters,crossprob,popfreqlistPP,popreqlistMM);
步骤8:从这次模拟操作中收集想要的统计学信息
为本次操作建立统计学信息:
simdata=Sim Validate(alldata,DE1,DE2,P1,P2,M1,M2,E1,E2,p1,p2,m1,m2,e1,e2,pe,pe,pp1,pp2,pm1,pm2);
模拟结果
进行一万次的模拟,对于该运算法则的准确性的最终评估PS准确性E1(PSAccuracy.E1)=mean(E1Ai),PS准确性E2(PSAccuracy.E2)=mean(E2Ai),这表明,通过E1、E2能够计算出PS运算法则的总体的准确性。对于一个个体的单核苷酸多态性而言,每个单核苷酸多态性的准确性:单核苷酸多态性准确性E1(SNPAcc.E1)=mean(E1Si)应当与正确测量该单核苷酸多态性的估算概率的平均值相适合,即与单核苷酸多态性概率E1(SNPProb.E1)=mean(E2Pi)相适应,即,如果该运算法则的运算是正确的,则单核苷酸多态性准确性E1(SNPAcc.E1)应当与单核苷酸多态性概率E1相关。上述两者之间的关系是通过它们的相关性反映出来的。
上述10000个模拟循环是在不同的设置类型下运行的:
(1)利用更加真实的已有的基因型数据而得到的潜在的种群频率,以及均一的种群频率,其中在每个单核苷酸多态性中,A,C,T,G具有相同的概率。
(2)对于单倍体的测量(ph)、无序的二倍体的测量(pd)、以及晶胚的测量(pe)的几种测量准确性的组合。模拟下述各种不同的准确性;非常准确(0.95,0.95,0.95),不太准确(0.75,0.75,0.75)以及不准确或者随机(0.25,0.25、0.26),以及不均衡的组合(0.9,0.5,0.5),(0.5,0.9,0.9)。最接近真实的情况应当是准确性大约为(0.6,0.8,0.8)的情况。
(3)我们进行的模拟是针对所有的情况同时使用PS运算法则以及简单匹配的STPS运算法则,这是为了评价本发明公开的运算法则的表现。
上述所有操作的结果在表8当中有所概况。
在这些模拟中,本发明公开的运算法则的表现要优于现有的经验性的运算法则,特别是在不均一的种群频率的真实情况中、以及在正确测量的不均衡概率或者降低的概率的真实情况中。同样被证明的是,在这些情况中,我们对于个体单核苷酸多态性的运算法则的准确性的估算是很好的,因为在平均比为1的情况下,对于正确的等位基因分型的估算的准确率与模拟的平均准确率的相关性是大约99%。
在最真实的情况中,当数据种群频率以及(ph,pd,pe)=(0.6,0.8,0.8)时,在执行方法1中,能够准确的重新获得(E1,E2)的单核苷酸多态性的平均百分比是(0.777,0.788),在执行方法2中,能够准确的重新获得(E1,E2)的单核苷酸多态性的平均百分比是(0.835,0.828)。如前所述,用来表示运算法则的平均准确性的数字不仅仅涉及正确的单核苷酸多态性分型,还涉及对上述单核苷酸多态性的双亲来源的准确识别。为了有效,与简单的接纳测量数据的运算法则相比,该运算法则必须还原出更加好的结果。可以惊讶的发现,在某些情况中,上述运算法则的准确性低于列出的测量的准确性。很重要的一点是,需要记住,为了达到这种模拟的目的,仅仅当一个单核苷酸多态性具有了正确的分型,并且其亲本以及染色体来源也被正确的识别之后,才能认为该单核苷酸多态性分型是准确的。偶然的获得这种正确性的几率被认为低于测量的准确性。
获得产前遗传数据以及胚胎遗传数据所必需的实验室技术
有许多可以利用的技术,能够用来进行细胞和DNA片段的分离,用以进行基因定型。可以将这里公开的体系以及方法应用于这些技术中的任何一项中,特别是那些包括从母亲血液中分离胎儿细胞或者DNA片段的技术,或者是那些在体外受精过程中从晶胚中分离胚泡的技术。这里公开的体系以及方法同样可以被应用于遗传数据的计算机数据分析方法中,即,不直接从遗传数据中进行测量。
在本发明体系的一种实施方式中,可以通过下面描述的方法获得这类数据。
细胞的分离
成人的二倍体细胞可以从大块儿组织或者血液样本中获得。成人的二倍体单个细胞可以通过流式细胞术(FACS)、或者荧光激活细胞分拣法从全血样本中获得。成人单倍体单个精子细胞同样可以通过流式细胞术从精子样本中分离出来。成人单倍体单个卵细胞可以在体外受精过程中在卵细胞收集的步骤中分离出来。
可以按照体外受精临床操作中普遍的技术完成目标单个胚泡从人类晶胚中的分离。可以使用单克隆抗体、或者其他技术例如流式细胞术或者密度梯度离心来完成目标胎儿细胞从母亲血液中的分离。
在本申请中,DNA的提取也可能需要使用到非标准的方法。对不同的DNA的提取方法进行比较的文献记载,在某些情况中,发现新的方案例如使用额外添加的N-十二烷基肌氨钠将更加有效并且产生最少的假阳性。
染色体组DNA的扩增
染色体组的扩增可以通过多种方法来完成,这些方法包括:连接介导的聚合酶链式反应(LM-PCR),简并寡核苷酸引物聚合酶链式反应(DOP-PCR),以及多重置换扩增(MDA)。在上述三种方法中,简并寡核苷酸引物聚合酶链式反应可靠的通过少量的DNA生成大量的DNA,其中也包括了单拷贝染色体;这个方法也许是最适合用于对双亲二倍体数据进行基因定型的方法,因为在这类基因定型的过程中数据的忠实度是至关重要的。多重置换扩增是最快速的方法,在几个小时之内生成了几百倍的DNA扩增;这个方法也许最适合用于对胚胎细胞进行基因定型的方法,因为在这类基因定型中,时间是实质因素。
背景扩增对于上述这些方法中的每一个来说都是个问题,因为每一个方法都可能扩增污染的DNA。非常微小量的污染都能够不可逆的破坏检测并且给出错误数据。因此,使用干净的实验室环境是至关重要的,因为在该环境中,扩增前以及扩增后的工作流程是完全被物理分离的。干净的、无污染的工作流程是工业化分子生物学中的一种新的途径,仅需要对细节进行仔细的关注。
基因定型检测以及杂交
对扩增的DNA的基因定型可以通过许多方法来完成,这些方法包括分子倒置探针(MIPs),例如Affymetrix’s Genflex标记阵列,微阵列,例如Affymetrix’s 500K阵列或者Illumina磁珠阵列,或者单核苷酸多态性基因定型检测,例如AppliedBioscience’s Taqman检测。Affymetrix 500K阵列、分子倒置探针/GenFlex,Taqman以及Illumina检测都需要微克数量级的DNA,因而利用任何的工作流程对单个细胞进行基因定型都需要某种类型的扩增。这些技术中的每一项都需要在以下方面进行各种不同的权衡:成本方面,数据质量方面,数据的数量质量比方面,用户自定义(customizability)方面,完成检测的时间方面以及可测量出的单核苷酸多态性的数量方面,以及其他方面。500K以及Illumina阵列的一个优点在于它们可以在大量的单核苷酸多态性的基础上收集数据,大约250000个单核苷酸多态性,它们与分子倒置探针相反,后者能够对大致(on the order of)10000个单核苷酸多态性进行检测,而Taqman检测能够检测的数量更少。分子倒置探针、Taqman以及Illumina检测优于500K阵列的一个优点在于它们固有的用户自定义化,允许使用者对单核苷酸多态性进行选择,而500K阵列则不允许这种用户自定义化。
在体外受精过程中的胚胎植入前的诊断情况中,固有的时间限制是很重要的;在这种情况下,牺牲数据的质量来换取时间可能是有利的。尽管还具有其他的明显的优势,但标准的分子倒置探针检测方案是一个相对而言时间集约型(time-intensive)的操作过程,通常需要2.5天到3天的时间来完成。在分子倒置探针方法中,探针从目标DNA处的退火以及扩增后的杂交是格外的时间集约,对于这些时间的任何背离都将导致数据质量的下降。探针从DNA样本中退火过夜(12-16个小时)。扩增后的杂交从阵列中退火过夜(12-16个小时)。在退火以及扩增的前后发生的其他许多步骤使得该方案的总体标准时间期限达到2.5天。对于分子倒置探针检测的速度的优化能够潜在的将这一过程的时间减少至少于36个小时。500K阵列以及Illumina检测具有一个更加快速的轮回(turnaround):在标准的操作方案中,需要大约1.5天至2天时间来生成高度可靠的数据。这些方法都是可以被优化的,据估计,对于500K阵列和/或Illumina检测而言,其进行基因定型检测的一个轮回的时间可以被减少至少于24个小时。更加快速的是Taqman检测,其可以在三个小时内完成。对于上述所有的方法而言,检测时间的减少将导致数据质量的降低,然而,这正是本发明公开的内容所要从事的。某些可利用的技术虽然快速,但不是特别的高产率,因此,对于此处的高度相似的产前遗传诊断而言是不可行的。
自然的,在时间是至关重要的情况中,例如在体外受精过程中对胚泡进行基因定型的情况中,快速的检测相对于较慢的检测而言具有明显的优势,而当不存在这种时间压力时,例如在体外受精之前对双亲DNA的基因定型已经开始时,在选择合适的方法时,其他国素将占主导地位。例如,在一项技术以及另外一项技术间的另外一种权衡针对的是价格与数据质量的比。下述做法是有意义的:当测量数据的质量更为重要时,使用能够获得高质量的测量数据的相对昂贵的技术,而当数据的忠诚度不是至关重要的时候,使用具有相对低质量的测量数据的相对便宜的技术。任何能够得到足够快速的高生产力的基因定型技术都可以用来与这种方法一同进行遗传物质的基因定型。
本发明所述方法的上下文实施例
这里描述了本发明公开的方法应用于体外受精实验室操作的一个实施例,在该操作中,能够在体外受精过程的时间限制内对所有可存活的晶胚进行完全的基因定型。在一次体外受精的实验室操作中,从卵细胞的受精到晶胚的植入所需的一个轮回的时间在三天之内。这意味着相关的实验室工作、数据的清除、以及表型的预测必须在这个时间内完成。这一体系的示意图在附图19中有所表示,并且在下面进行描述。这一体系可以由下述步骤组成:从体外受精对象(母亲)1902以及体外受精对象(父亲)1903中提取双亲遗传样本1901,利用基因定型体系在体外受精实验室操作下1904对其进行分析。该体系可以包括从母亲1902处收集的多个卵细胞,并且使用从父亲1903处获得的精子对上述卵细胞进行受精,从而得到多个受精的胚胎1905。该体系可以包括一个实验室技师,该技师为每个晶胚提取胚泡,对每个胚泡进行DNA的扩增,并且使用高生产力的基因定型体系1906对其进行分析。该体系可以包括将来自于双亲的遗传数据以及来自于胚泡的通传数据传送至安全的数据处理体系1907中,其中所述的安全的数据处理体系验证并且清除胚胎的遗传数据。该体系可以包括使用表型运算法则1909对清除过的胚胎数据1908进行操作,从而为每个晶胚预测表型易感性。该体系可以包括将这些预测值与其相关的置信度水平一同传送至内科医师1910处,其中所述的内科医师帮助体外受精对象1902以及1903来选择以供植入母体1901内的晶胚。
涉及到遗传数据的清除的杂项
很重要的一点是,应当注意,这里描述的方法涉及的是遗传数据的清除,并且由于所有的存活生物都含有遗传数据,因而这些方法可以相同的应用于从双亲中遗传染色体的任何人类、动物、或者植物。这些动物以及植物的列表可以包括,但不局限于:大猩猩,黑猩猩,倭黑猩猩,猫,狗,熊猫,马,牛,绵羊,山羊,猪,印度豹,老虎,狮子,鲑鱼,鲨鱼,鲸,骆驼,野牛,海牛,麋鹿,旗鱼,海豚,犰狳,黄蜂,蟑螂,蠕虫,秃鹫,鹰,麻雀,蝴蝶,美洲杉,玉米,小麦,大米,矮牵牛花,兰花苕子,向日葵,豚草,橡树,栗树,以及头虱。
对遗传数据的测量不是一个完美的过程,特别是当遗传物质的样本很小时。这些测量值通常包含错误的测量值、不清楚的测量值、伪造的测量值、以及丢失的测量值。这里公开的方法的目的在于检测并且校正所有这些错误或者其中的某些错误。使用这个方法能够提高置信度,对遗传数据进行很大程度上的了解。例如,使用当前的技术,从单个细胞中扩增得到的DNA中不清楚的测量得到的遗传数据中可能含有20%至50%的未测量的区域,或者等位基因的脱失。在某些情况中,这些遗传数据中可能含有1%至99%的未测量的区域,或者等位基因的脱失。除此之外,给出单核苷酸多态性测量值的置信度也受到错误的影响。
当不清楚的数据具有大约50%的等位基因脱失率时,在这种情况下,期望的是在应用了这里公开的方法之后,在至少90%的情况中清楚的数据将具有正确的等位基因分型,在理想状态下,这一数值将升至99%甚至更高。当不清楚的数据具有大约80%的等位基因脱失率时,在这种情况下,期望的是在应用了这里公开的方法之后,在至少95%的情况中清楚的数据将具有正确的等位基因分型,在理想状态下,这一数值将升至99.9%甚至更高。当不清楚的数据具有大约90%的等位基因脱失率时,在这种情况下,期望的是在应用了这里公开的方法之后,在至少99%的情况中清楚的数据将具有正确的等位基因分型,在理想状态下,这一数值将升至99.99%甚至更高。当一个特定的单核苷酸测量值的置信度接近90%时,在这类情况中,所期望的是,清楚的数据具有的等位基因分型的置信度超过95%,在理想的状态下,该置信度超过99%甚至更高。当一个特定的单核苷酸测量值的置信度接近99%时,在这类情况中,所期望的是,清楚的数据具有的等位基因分型的置信度超过99.9%,在理想的状态下,该置信度超过99.99%甚至更高。
同样重要的是,需要注意,通过测量来自于一个卵裂球的扩增DNA而得到的胚胎遗传数据可以被用于很多目的中。例如,可以将其用来检测异倍体、单亲二体性、对个体进行性别鉴定、以及进行各种表型的预测。当前,在体外受精的实验室操作中,由于所使用的技术的关系,通常产生下述情形:一个卵裂球仅能提供检测一种疾病所需要的足够的遗传物质,所述的疾病例如是异倍体,或者是具体的单性生殖疾病。由于这里公开的方法具有一个常见的第一步骤:对来自于卵裂球的一大组单核苷酸多态性进行检测,不考虑需要进行的预测的类型,则内科医师或者家长不需要被迫选择需要筛选的有限数量的疾病。相反,允许对医学知识状态所已知的许多基因和/或表型的筛选进行选择。使用本发明公开的方法,在对卵裂球进行基因定型之前识别特定的筛选条件的仅有的优点在于:如果确定了某些表型相关性单核苷酸多态性(PSNPs)是特别相关的,那么可以选择一组更加适合的非表型相关性单核苷酸多态性(NSNPs),这类非表型相关性单核苷酸多态性更可能与目标表型相关性单核苷酸多态性发生共分离,因此,增加了目标等位基因分型的置信度。值得注意的是,即使当单核苷酸多态性没有提早的被个人化(personalized),在这种情况中,其置信度也被预期为足够适合用于这里描述的各种目的当中。
表型预测以及临床预测
可以利用许多模型通过基因型信息以及临床信息来预测表型数据。不同的模型适合于不同的情况中,以可利用的数据的量以及类型为依据。为了选择表型预测的最适合方法,通常最好使用一组检验数据对多种方法进行检验,并且决定哪种方法提供了最佳的预测准确性,这种预测准确性是基于与检验数据的测量结果相比较而得到的结论。这里描述的某些实施方式中包括了一组方法,当将这组方法组合使用并且基于检验数据的表现而进行选择,这组方法将提供获得准确的表型预测值的高度可能性。首先,描述了使用列联表在种类(scenario)(ii)中进行基因型-表型模拟的技术。其次,描述了使用由凸集合最佳技术建立的回归模型在种类(iii)中进行基因型-表型模拟的技术。之后,描述了用来选择最佳模型的技术,其中所述的最佳模型给定了被预测的特定表型,特定的患者数据,以及用于训练并且检验模型的特定的数据组。
如今的数据:基于列联表的模拟表型结果
如果存在已知的遗传缺陷以及能够增加疾病表型的概率的等位基因,并且用以进行预测的对象(predictors)的数量足够的少,在这种情况下,可以通过列联表来模拟表型的概率。如果仅存在一个相关的遗传等位基因,可以使用A+/A-来描述一个特定等位基因的存在/缺失,并且使用D+/D-来描述一个疾病表型的存在/缺失。包含(f1,N1,f2,N2)的列联表是:
|
D+ |
D- |
# |
G+ |
f1 |
1-f1 |
N1 |
G- |
f2 |
1-f2 |
N2 |
其中f1以及f2表示测量频率或者不同结果的概率,并且这些条目的总数是N=N1+N2。从这个表中可以看出,在具有自变量(IV)G+或者G-的两种情况中,存在疾病状态D+的概率的让步比(odds ratio)可以被表示为OR=f1(1-f2)/f2(1-f1),其具有95%的置信区间(confidence interval):OR1±1.96/S,其中S表示标准偏差。例如,利用一项对10000个个体进行的乳腺癌研究,其中M+表示BRCA1或者BRCA2等位基因的存在:
|
D+ |
D- |
# |
M+ |
.563 |
.437 |
1720 |
M- |
.468 |
.532 |
8280 |
这一数据得出了一个让步比,OR=1.463,其置信区间为[1.31;1.62],它可以被用来预测具有给定突变的乳腺癌发生的增加的概率。值得注意的是,大于2*2的列联表可以被用来供给更多的自变量或者结果变量。例如,在乳腺癌的情形中,列联M+以及M-可以被下述四种列联所替代:BRCA1以及BRCA2,BRCA1以及非BRCA2,非BRCA1以及BRCA2,以及最后的非BRCA1以及非BRCA2。具有本领域知识的人应当很清楚怎样确定多于2*2的列联表的置信区间。这项技术将被用在当没有足够的自变量以及足够的数据来建立具有低标准偏差的模型的情况中,通过对不同组中的患者进行计算来进行该项技术,其中所述的不同的组被自变量的不同列联所定义。这种方法避免了设计数学模型的困难,其中所述的数学模型将不同的自变量与需要被模拟的结果相关起来,在构建回归模型的过程中需要使用到这种数学模型。
值得注意的是,来自于特定的单核苷酸多态性的遗传数据也可以被投射(project)至其他的自变量空间中,特别是例如在HapMap(人类基因组单体型图)计划被识别的不同类型的单核苷酸多态性。HapMap(人类基因组单体型图)计划将个体集合在库(bin)中,并且使用一个特定类型的单核苷酸多态性对每个库进行定义。例如,一个库(B1)具有的单核苷酸多态性类型中含有BRCA1以及BRCA2,另外一个库(B2)具有的单核苷酸多态性类型中含有BRCA1以及非BRCA2,第三个库(B3)具有的单核苷酸多态性的类型关系到了突变的所有其他组合。与建立一个表示这些单核苷酸多态性的全部不同组合的列联表相比,更胜一筹的做法是建立一个表示列联B1,B2以及B3的列联表。
此外还需要注意的是,如HapMap(人类基因组单体型图)计划中所述,某些单核苷酸多态性同时出现的趋向可以被用来建立模型,所述的模型利用多种单核苷酸多态性作为预测因素,即使那时所得的数据是由不同组的患者数据所组成的,其中每个组仅具有一个被测量了的单核苷酸多态性。这一问题通常出现在从公众可以获得的研究论文中建立模型的情况中,其中所述的研究论文例如是那些从OMIM(人类孟德尔遗传在线)中获得的论文,其中每篇论文中含有的数据都是同一类的,仅对一个相关的单核苷酸多态性进行了测量,尽管有多个单核苷酸多态性对于表型来说都是具有预言性的。为了对这一方面的内容进行阐述,其中这些内容能够有效的用于通过当今可获得的数据建立预测模型的方法中,利用阿尔茨海默氏症作为具体的参考,对于该病症的预测模型可以依据下述自变量而建立:阿尔茨海默氏症的家族历史,性别,种族,年龄,以及三个基因的不同等位基因,所述的三个基因是APOE(血浆载脂蛋白E),NOS3(一氧化氮合成酶-3),以及ACE(血管紧张素转换酶)。在该疾病的情况中,对除阿尔茨海默氏症之外还可以存在于其他疾病中的普遍深入的问题进行了讨论:尽管在确定一种特定表型的倾向的过程中涉及到许多的基因,但绝大多数的历史研究仅仅对一个特定基因的等位基因进行了取样。对于阿尔茨海默氏症来说,几乎所有被研究的同类(cohorts)都仅仅进行了一个基因的取样,所述的基因是APOE(血浆载脂蛋白E),NOS3(一氧化氮合成酶-3),或者ACE(血管紧张素转换酶)。然而,输入多个遗传等位基因来进行模型的建立是很重要的,即使大部分可利用的数据是来自于仅仅调查了一个基因的研究中的。这一问题可以通过以下方面进行解决:对一种简单化的情形进行考虑,在该情形中存在两种表型状态,以及表示两个相关基因的两个自变量,每个自变量都具有两种状态。给定一个用以描述疾病表型的随机变量D∈[D+,D-],以及用以描述上述基因的两个随机变量A∈[A+,A-]以及B∈[B+,B],其目的在于找到最可能的估计值P(D/A,B)。可以通过应用贝叶斯定理来实现上述目的,利用P(D/A,B)=P(A,B/D)P(D)/P(A,B)。P(D)以及P(A,B)可以从公共数据中获得。具体的,P(D)指的是该疾病在人群中的总体普及率,这一数值可以从公众可获得的统计学信息中得到。除此之外,P(A,B)指的是基因A与基因B在一个个体中同时发生的特定状态的普及率,这一数值可以从公开的数据库中得到,例如从HapMap(人类基因组单体型图)计划中得到,该计划在不同种族的多个个体中对各种不同的单核苷酸多态性进行了测量。需要注意的是,在一个优选的实施方式中,所有这些概率都是针对特定的种族群体以及特定的性别来计算的,其中存在概率偏差,其结果要优于针对整个人类群体进行计算的结果。一旦确定了这些概率,挑战来自于对P(A,B/D)的准确估计,因为大部分的同类数据(cohort data)提供的是P(A/D)以及P(B/D)的估计值。相关的信息可以在各种公开的数据库中找到,例如HapMap(人类基因组单体型图)计划,这些信息是关于不同的遗传等位基因间的统计学关联的信息,即,关于P(A/B)的信息。然而,仅仅给出P(A/B)、P(A/D)、P(B/D)对于得到P(A,B/D)来说不起任何作用,因为存在不受拘束的自由度(unconstrained degree of freedom)。然而,如果已知关于P(A,B/D)的某些信息,并且这些信息是从一类对基因A以及基因B同时进行了取样的同类中得到的,即使仅仅是一个单一的列联例如(A-,B-),那么,与P(A/D)、P(B/D)、P(A/B)有关的信息的价值将被进行平衡,从而促进对P(A,B/D)的估计。将使用列联表对这一概念进行描述。
针对下面的两个列联表进行考虑,这两个列联表表示的是受遗传状态A+以及A-所影响的结果概率D+以及D-。这项研究针对的是状态A。使用f来表示A的测量频率,并且使用p来表示我们试图进行评估的实际概率。
A |
D+ |
D- |
A+ |
f1 |
f2 |
A- |
f3 |
f4 |
A |
D+ |
D- |
A+ |
p1 |
p2 |
A- |
p3 |
p4 |
其中f3=1-f1,f4=1-f2并且p3=1-p1,p4=1-p2。用K1来表示在状态A的情况中条目(subjects)的数量,即,产生结果D+的条目的数量。用K2来表示状态A的对照组中条目的数量,即,产生结果D-的条目的数量。
相类似的,针对下面的两个列联表进行考虑,这两个列联表表示的是受遗传状态B+以及B-所影响的结果概率D+以及D-。这项研究针对的是状态B。使用f来表示B的测量频率,并且使用p来表示我们试图进行评估的实际概率。
B |
D+ |
D- |
B+ |
f5 |
F6 |
B- |
f7 |
F8 |
其中f7=1-f5,f8=1-f6并且p7=1-p5,p8=1-p6。使用K3来表示在状态B的情况中的数量,并且使用K4来表示状态B的对照组中的数量。上述的列联表表示的是分别测量遗传状态A以及遗传状态B的试验。然而,理想情况中选取的列联表应当包括不同状态的A以及B的组合。下面描述的列联表是一种假设的研究,针对的是AB组合的状态,其中使用f来表示测量频率,并且使用p来表示实际概率。
AB |
D+ |
D- |
A+B+ |
f9 |
f10 |
A+B- |
f11 |
f12 |
A-B+ |
f13 |
f14 |
A-B- |
f15 |
f16 |
AB |
D+ |
D- |
A+B+ |
p9 |
p10 |
A+B- |
p11 |
p12 |
A-B+ |
p13 |
p14 |
A-B- |
p15 |
p16 |
其中f15=1-f9-f11-f13,f16=1-f10-f12-f14并且p15=1-p9-p11-p13,p16=1-p10-p12-p14。使用K5来表示在状态AB的情况下的数量,并且使用K6来表示状态AB的对照组中的数量。
为了使用符号进行表示,记作K7=K9=K5并且K8=K10=K6。因此,实际上各组的规模为:
# |
D+ |
D- |
A |
K1 |
K2 |
B |
K3 |
K4 |
AB |
K5 |
K6 |
可以使用统计学的基本规则来加强假设的列联表AB中的细胞间的相关性。在这个实施例中,对于符合状态D+的细胞而言,下述相关性可以被加强:
P(A+B-/D+)=P(A+/D+)-P(A+B+/D+)
P(A-B+/D+)=P(B+/D+)-P(A+B+/D+)
P(A-B-/D+)=1-P(A+/D+)-P(B+/D+)+P(A+B+/D+)
并且类似的,对于符合状态D-的细胞而言:
P(A+B-/D-)=P(A+/D-)-P(A+B+/D-)
P(A-B+/D-)=P(B+/D-)-P(A+B+/D-)
P(A-B-/D-)=1-P(A+/D-)-P(B+/D-)+P(A+B+/D-)
使用上述列联表中的符号,不考虑多余存在的相关性,这些相关性可以被转化为:
p11=p1-p9
p13=p5-p9
p12=p2-p10
p14=p6-p10
或者等价变换为
p1=p9+p11
p2=p10+p12
p5=p9+p13
p6=p10+p14
为了对所有这些相关性进行概括,下面给出了p1,...,p16与p9,...,p16的所有相关性。为了得出这些值之间的相关性,在横行中的概率是在纵列中具有数值1的概率的总和,例如,第一横行给出p1=p9+p11。
|
p9 |
p10 |
p11 |
p12 |
p13 |
p14 |
p15 |
P16 |
p1 |
1 |
|
1 |
|
|
|
|
|
p2 |
|
1 |
|
1 |
|
|
|
|
p3 |
|
|
|
|
1 |
|
1 |
|
p4 |
|
|
|
|
|
1 |
|
1 |
p5 |
1 |
|
|
|
1 |
|
|
|
p6 |
|
1 |
|
|
|
1 |
|
|
p7 |
|
|
1 |
|
|
|
1 |
|
p8 |
|
|
|
1 |
|
|
|
1 |
p9 |
1 |
|
|
|
|
|
|
|
p10 |
|
1 |
|
|
|
|
|
|
p11 |
|
|
1 |
|
|
|
|
|
p12 |
|
|
|
1 |
|
|
|
|
p13 |
|
|
|
|
1 |
|
|
|
通过频率与概率的相关性,可以建立等式fi=pi+ni,其中n=9......16,其中ni是干扰项,表示的是基于fi的出现频率的概率pi的不完美测量。将这一方法应用于上面描述的相关性中,并且假定已经对列联表AB中的所有细胞进行了测量(这样做的目的仅仅在于描述,将在以下的内容中进行讨论),这10个观测结果可以被表示为:
这些测量等式可以通过矩阵符号(matrix notation)来呈现为:
F=XP+N
其中F=[F1,...,F16]T,P=[p9,...,p16]T并且N=[n9,...,n16]T,其中X是前述表格中所表示的矩阵。这个矩阵方程式可以用来计算8个未知的系数p9...p16。在这个特定的情况中,我们计算了所有的参数p9...p16。如果我们不能获得组合的A、B基因的所有测量值,那么我们至少需要知道一个D+的测量值以及一个D-的测量值。给定了上述的相关性,之后,我们可以填写表格的其余部分。换句话说,为了能够填写针对假设的研究AB的列联表,期望具有至少一种样本,对该样本中存在的A以及B的特定状态的某些条目同时进行测量,其中这些条目同时具有D+以及D-的结果。这使得我们能够得到矩阵X的满秩(full rank),表示的是进行的测量,因此可以得到p9...p16的值,并且可以将这些值填入AB的列联表中。如果还存在更多的研究数据,可以在矩阵X的底部添加更多的横行,这些横行的结构与上面给出的横行的结构相类似。
为了完成精确的回归,需要进行一个权重回归,这一回归对于每一个由组间样本的大小来决定的观测值(observation)fi进行权重分析,因而,具有更多观测值的研究以及细胞将获得更多的权重。对于测量等式fi=pi+ni而言,ni并不是总具有相同的变量,并且上述回归不是同方差性的(homoscedastic)。具体的,fi=1/Ki *二项式(pi,Ki)~N(pi,pi(1-pi)/Ki),其中二项式(pi,Ki)表示的是一个二项分布,其中每项检验的结果概率是pi,并且进行了Ki项检验。这个二项分布可以由N(pi,pi(1-pi)/Ki)来近似计算,这是一个正态分布,其具有平均值pi,以及变量pi(1-pi)/Ki。因此,干扰(噪声)可以被模拟为一个正态变量ni~N(0,pi(1-pi)/Ki),其具有理论变量Vi=pi*(1-pi)/Ki。这一变量可以通过样本频率vi=fi*(1-fi)/Ki来近似计算。
进行一种权重回归,对每项观测值i进行权重分析,所述的观测值与变量vi成反比。可以将干扰矩阵N的分布描述为~N(0,V),其中V是一个具有对角元素[v9,...,v16]的矩阵,并且其他所有的元素都是0。该矩阵可以被表示为V=diag([v9,...,v16])。相类似的,使W=diag([1/v9,...,1/v16])。现在,可以使用权重回归对P进行计算:
P=(X’WX)-1X’WY
可以直接看出变量P可以是
Var(P)=(X’WX)-1
其可以被用来得出P值计算的置信度。
为了进行概括,我们使用了来自于个体基因的数据(A:f1,...,f4,B:f5,...,f8),以及来自于A和B的组合的数据(AB:f9,...,f16),从而有助于估算A和B组合的概率(p9,...,p16)以及它们的变量(v9,...,v16)。最后,在我们的研究中,我们最常处理的是对数让步比而不是概率,因此,我们需要将这些概率转化为LORs(对数让步比)。一般的,如下给出了事件H的概率以及变量。
|
D+ |
D- |
H+ |
p1 |
p2 |
H- |
1-p1 |
1-p2 |
V |
v1 |
v2 |
对数让步比的公式为LOR=[log(p1)-log(1-p1)]-[log(p2)-log(1-p2)],变量的公式为(通过δ方法)V=[(p1)-1+(1-p1)-1]2*V(p1)+[(p2)-1+(1-p2)-1]2*V(p2)。下面的表格所列出的概率与对数让步比以及变量相符合,针对的是A,B组合的情况
|
D+ |
D- |
LOR |
Var |
A+B+ |
P9 |
p10 |
lor1 |
V1=[1/p9+1/(1-p9)]2v9+[1/p10+1/(1-p10)]2v10 |
A+B- |
P11 |
p12 |
lor2 |
V2=[1/p11+1/(1-p1)]2v1+[1/p12+1/(1-p12)]2v12 |
A-B+ |
P13 |
p14 |
lor3 |
V3=[1/p13+1/(1-p13)]2v3+[1/p14+1/(1-p14)]2v14 |
A-B- |
P15 |
p16 |
lor4 |
V4=[1/p15+1/(1-p15)]2v15+[1/p16+1/(1-p16)]2v16 |
该表提供了对数让步比的估算值以及相应的变量的估算值。
作为对该方法的说明,该技术可以被用来获得改进的P(A,B/D)的估算值,其中D表示的是存在阿尔茨海默氏症的状态,而A以及B分别表示的是APOE(血浆载脂蛋白E)基因以及ACE(血管紧张素转换酶)基因的两种不同的状态。表9表示的是三次不同的研究:由Alvarez于1999年进行的研究,其中仅仅对基因A进行了取样;由Labert于1998年进行的研究,其中仅仅对基因B进行了取样;以及由Farrer于2005年进行的研究,其中对基因A以及基因B都进行了取样。从这些研究中获得了两组结果,这两组结果被表示在表10当中。第一组(参见表10,第2列,第3列,第4列以及第5列)对所有的同类(cohorts)进行了分析,并且利用这里公开的方法对P(A,B/D)进行了改进的估算,给出P(A/D)以及P(B/D)。第二组(参见表10,第6列,第7列,第8列以及第9列)仅仅利用了由Farrer(于2005年)研究的现代同类而生成的那些结果,对P(A,B/D)进行了估算,在该方法中两种基因都被进行了取样。在前一种情况中,预测值的置信界限明显的减小。值得注意的是,可以使用来自于公共来源的描述P(A/B)的数据对这些预测值进行进一步的改进——可以将这些测量值加入到如上所述的X矩阵中。同样需要注意的是,这里描述的技术也可以被用来对单独的A概率以及B概率的估算值进行改进,例如对P(A+/D+),P(A+/D-),P(B+/D+)以及P(B-/D-)进行改进,所述的改进需要利用如上所述的相关性,例如p1=p5+p7。
需要注意的是,虽然仅仅使用了两个变量A以及B对该方法进行了说明,但应当注明的是上述列联表可以包括许多不同的自变量,例如那些在阿尔茨海默氏症的预测方法中提到的自变量:阿尔茨海默氏症的家族历史,性别,种族,年龄,以及三个基因的不同等位基因,所述的三个基因是APOE(血浆载脂蛋白E),NOS3(一氧化氮合成酶-3),以及ACE(血管紧张素转换酶)。可以对例如年龄之类的连续变量进行分类处理,将其归为不同数值的库(bins),以使其适应列联表信息。在一个优选的实施方式中,使用最大数量的自变量来模拟一种结果的概率,上述概率的标准偏差通常低于某些具体的阈值。换句话说,给定一个特定患者可利用的自变量,可能生成最具体的列联,而保留与该列联足够相关的训练数据,能使相关概率的估算值变得有意义。
值得注意的是,在阅读过本发明公开的内容之后,本领域技术人员应当明了,可以怎样利用类似的技术来提高多变量线性回归模型、多变量非线性回归模型、以及逻辑回归模型的准确性,其中所述的技术利用到了关于疾病-基因相关性的数据,基因-基因相关性的数据,和/或种群中的基因频率数据。进一步的,在阅读过本发明公开的内容之后,本领域技术人员应当明了,可以怎样利用类似的技术来提高多变量线性回归模型、多变量非线性回归模型、以及逻辑回归模型的准确性,其中所述的技术利用到了关于疾病-基因相关性的数据,基因-基因相关性的数据,和/或种群中的基因频率数据,并且该技术是通过对结果数据进行平衡来达到训练模型的目的,在所述的模型中,并不是对所有的与该模型相关的自变量中的所述结果数据进行测量。进一步的,在阅读过本发明公开的内容之后,本领域技术人员应当明了,可以怎样利用类似的技术来提高列联表模型的准确性,其中所述的技术利用到了关于疾病-基因相关性的数据,基因-基因相关性的数据,和/或种群中的基因频率数据,并且所述的列联表模型是通过其他技术来构建的,例如通过本领域熟知的期望值最大化(EM)运算法则来构建的。这些技术对于平衡来自于HapMap(人类基因组单体型图)计划的数据以及被包含在公共数据库中的其他数据而言是特别相关的,其中所述的公共数据库是例如美国国立生物技术信息中心(NCBI),人类孟德尔遗传在线(OMIM)以及单核苷酸多态性(dbSNP)数据库。
同样需要注意的是,在本专利中,当我们谈到来自于个体或者宿主的数据时,同时也假定了该数据可以指代该宿主中感染的任何病原体,或者该宿主患有的任何癌症。所述的个体数据或者宿主数据也可能指的是来自于人类晶胚的数据,来自于人类卵裂球的数据,来自于人类胎儿的数据,来自于某些其他细胞或者细胞组的数据,或者是来自于任意种类的动物或者植物的数据。
未来的数据:利用回归模型模拟多因子表型
随着越来越多的数据累积不断的将基因型与多因子表型联系起来,则如上所述的类型(scenario)(iii)将成为主要的类型,即需要对遗传标记的复杂组合进行考虑,从而准确的预测表型,并且将要调用多维线性回归模型或者多维非线性回归模型。通常,在针对这一类型进行模型训练时,可能的预测因素(predictors)的数量将会大于测量结果的数量。这里描述的体系以及方法的例子包括一项新的技术,该项技术可以为未确定的基因型-表型数据组或者不适定的(ill-conditioned)基因型-表型数据组建立稀疏参数模型。通过着重模拟人类免疫缺陷病毒(HIV)/获得性免疫缺陷综合症(AIDS)对抗逆转录治疗(ART)产生的应答来对该项技术进行说明,其中具有很多可以用来进行比较的模拟操作,并且可以利用包含了许多可能的遗传预测因素的数据。当利用真实的实验室测量值通过交叉验证进行检验时,这些模型预测的药物应答表型比先前在文献中讨论的模型以及这里描述的其他规范技术的结果更加准确。
下面对两种回归技术进行描述以及说明,针对预测病毒表型对抗逆转录治疗的应答这一方面来进行描述以及说明,其中所述的预测是根据遗传序列数据来进行的。两种技术都使用了凸集合最佳技术对一组稀疏的模型参数进行了连续子集选取。第一种技术利用最小的绝对缩减和变量选择算子(Least Absolute Shrinkage and Selection Operator)(LASSO),该方法应用l1常态损失函数(norm loss function)来建立一个稀疏的线性模型;第二种技术利用支持向量机(SVM)与径向基核函数(radial basis kernel functions),该方法应用ε-不敏感损失函数来建立一个稀疏的非线性模型。上述技术被用来预测人类免疫缺陷病毒-1病毒对十种逆转录酶抑制剂(RTIs)药物以及七种蛋白酶抑制剂(PIs)药物的应答。所述的遗传数据来自于人类免疫缺陷病毒为逆转录酶以及蛋白酶编码的序列。使这类方法具有这种特性的主要特征在于损失函数趋向于生成简单模型,在该简单模型中许多参数为零,而价值函数的凸性(convexity)确保了可以找到模型参数,该参数能够使对于一组特定的训练数据的价值函数全局最小化。
最小的绝对缩减和变量选择算子(LASSO)以及L1选择函数
当预测因素(predictors)M的数量超过了训练样本N的数量时,模拟的问题就变得冗余(overcomplete),或者不适定(ill-posed),因为任何任意的N个预测因素子集足可以在训练数据中建立一个具有零误差的线性模型,只要X矩阵中相关的纵列是线性独立的。因此,我们不愿意信任一个由线性回归方法还原出的N-预测因素模型。然而,设想一下,一个具有明显少于N个变量的模型将具有低的训练误差。模型越稀疏,由偶然的人为因素而引起的低的训练错误的可能性越小;因此预测因素有原因的与依变量相关的可能性越大。在针对逆转录酶抑制剂治疗数据的情形中,这成为稀疏解在冗余问题中的重要性的基础。相类似的论点可以被应用于不适定(ill-conditioned)的问题当中,在针对蛋白酶抑制剂数据的情形中,用矩阵X
TX中用一个大的条件数(condition number)来表征该不适定问题。在这种情况中,估算的参数
对模型误差以及测量干扰具有高度的易感性,并且因此而不可能被准确的概括出。冗余问题以及不适定问题是遗传数据中典型的问题,在上述情形中,可能的预测因素的数量——基因、蛋白、或者在我们所列举的情况中为突变位点——要大于测量结果的数量。
针对上述情况的一种规范的做法是子集的选取。例如,通过进行逐步选取,在每个步骤中将一个预测因素添加到模型中,这种添加是以获得最高的F检验结果为基础的,所述的F检验能够从统计学上表明该变量与预测误差的相关性的显著性水平。在每个变量被添加之后,可以对其余的所有变量进行检查,以确保这些变量中没有任何一个变量与该模型的预测误差的相关性低于统计学显著性阈值。这项技术已经被成功的运用于对药物应答预测的问题中。然而,由于选取过程的离散特性,在数值上发生的小的变化可以相当大的改变预测因素的选取组(chosen set)。一个变量的存在或者不存在可能影响到其与另外一个变量的相关性的统计学显著性,并且影响到是否将该变量导入到模型中或者离开该模型。这将影响归纳概括的准确性,特别是对于不适定的问题而言。
另外一种方法是依靠缩减函数来约束估算参数
的值。一种规范的缩减函数是参数的平方和,将这种缩减函数应用于岭回归,依照下述计算式得出该参数:
(17)
其中λ是一个调整参数(tuning parameter),通常由交叉验证来确定其值。这种方法是非稀疏的,并且不用将参数设置为零。这种方法趋向于破坏归纳概括(generalization)的准确性,并且使所得的解很难被翻译。
可以通过最小的绝对缩减和变量选择算子(LASSO)技术来克服上述问题。与子集的选取相反,最小的绝对缩减和变量选择算子(LASSO)不必进行离散的预测因素的接受或者拒绝;其允许我们通过一个连续子集优化的过程一同进行选取,一同选取的变量组便是最有效的预测因素。该技术使用l1常态损失函数:
(18)
其中λ通常是由交叉验证来设定的。最小的绝对缩减和变量选择算子(LASSO)将趋向于把许多参数设置为零。附图20表明了对最小的绝对缩减和变量选择算子(LASSO)的这一特性的洞察,称之为选择性。假定利用下述训练数据X=[1 0;01]T,y=[2 1]T来建立一个仅基于两个变异的模型,并且x轴以及y轴分别表示两个参数b1以及b2。对l1缩减函数以及l2缩减函数的作用进行比较,在两种函数中,得出一个非常相等的符合训练数据的值,如此一来||y-Xb||2=2。大圆(2001),小圆(2002),以及正方形(2003)分别表示价值函数||y-Xb||2,l2常态||b||2,以及l1常态|b1|+|b2|的层级曲线(level curve)。在两个圆的相交处(2004)得出岭回归(l2)的值;在正方形与大圆的相交处(2005)得出最小的绝对缩减和变量选择算子(LASSO)(l1)的值。由于l1常态的层级曲线的“尖端(pointiness)”得出的值位于轴b1之上,因而是稀疏的。将这一论点扩展到更多维的空间中,能够解释最小的绝对缩减和变量选择算子(LASSO)生成稀疏解的趋向性,并且暗示了得到的结果的可测量程度优于现有技术文献中报道的可测量程度的原因。
l1常态可以被认为是最具选择性的缩减函数,尽管其保留有凸集合。凸集合性确保能够针对一个给定的数据组得出一个全局解(global solution)。一项被称为最小角回归(Least Angle Regression)的近来新出现的高效的运算法则能够确保在M步骤中(M steps)集中至最小的绝对缩减和变量选择算子(LASSO)的全局解。
需要注意的是,在阅读过本发明公开的内容之后,本领域技术人员应当清楚,怎样将l1常态也应用于逻辑回归的情况中,从而模拟分类变量的每个状态的概率。在逻辑回归中,可以依照一组测量值的后验概率的相反值来建立一个凸集合价值函数。所述的后验概率指的是当假定模型估算的是每种结果的可能性时,观察到的训练数据的概率。通过向凸集合价值函数中加入l1常态,可以是得到的凸集合价值函数最小化,从而得到一个稀疏参数模型,用以模拟特定结果的概率。当测量结果的数量小于预测因素的数量时,使用逻辑函数的l1常态将是特别适合的。
支持向量机以及L1-常态
可以对支持向量机(SVM)进行设定,使其很好的模拟药物应答以及其他表型,特别是在所述的模型中的自变量之间存在复杂的相互作用的情况中。支持向量机的训练法则隐含的使用了l1常态选择函数。支持向量机是一类学习算法(learning algorithms),能够完成对实值函数(real-valued function)的近似计算,并且能够对样本数据进行准确的归纳概括,即使所估计的问题是阿达玛意义(Hadamard sense)上的不适定。支持向量机进行准确的归纳概括的能力通常会被支持向量机模型以及训练法则中的两个选择性的特征所影响。第一个特征是对价值函数的选择,或者是对在训练过程中被最小化了的函数的选择。第二个特征是对支持向量机的核函数(kernels)的选择,或者是对那些能够使支持向量机映射出复杂的非线性函数的函数的选择,包括自变量间的相互作用,其中所述的映射(map)使用的是一组相对小的线性回归参数。这些特征将在以下的内容中进行讨论。
利用线性函数近似值: 对宿主i的表型yi进行模拟。首先,通过使一个价值函数最小化来估算b,所述的价值函数由一个针对参数的l2缩减函数以及“ε-不敏感损失”函数组成,该函数不会降低误差至低于某些ε>0。可以对支持向量回归进行下述最优化:
(19)
受制于下述约束:
yi-bTxi≤ε+ξi +,i=1...N
(20)
bTxi-yi≤ε+ξi -,i=1...N
(21)
ξi +≥0,ξi -≥0,i=1...N
(22)
上述价值函数的第二项对模拟误差的绝对值进行最小化,使其超越了“不敏感”阈值ε。参数C能够对误差与缩减之间的相对重要性的权重进行衡量。为了满足Kuhn-Tucker约束,可以使用寻找拉格朗日(Lagrangian)鞍点(saddle-point)的标准方法来解答这种约束最佳化。所述的用于调节上述价值以及约束的拉格朗日函数为:
(23)
最小化针对的是参数向量b,ξ-,ξ+,而最大化针对的是拉格朗日乘子向量α-,α+,λ-,λ+.。需要注意的是,依照Kuhn-Tucker约束,拉格朗日乘子被期望是正的。因此,可以根据下式来得出最佳的参数组合:
(24)
其中
αi +,αi -,λi +,λi -≥0,i=1...N
(25)
由于最小化/最大化的顺序可以互相改变,首先对变量b,ξi +,ξi -进行最小化,通过将L的偏微商(partial derivatives)中的那些变量设置为0来进行上述最小化。从得到的等式中,可以发现权向量可以被表示为
(26)
同样,从得到的等式中,消去拉格朗日函数中的变量,通过对得到的二次形式(quadratic form)进行最大化,可以得出下列系数αi +,αi -,i=1...N:
其中
0≤αi +≤C,i=1...N
(29)
0≤αi -≤C,i=1...N
(30)
这样能够计算出向量b,并且将支持向量机模型充分定义为ε-不敏感损失函数。从等式(11)中可以看出,该模型可以被定义为
(31)
其中βi=αi +-αi -。得到的模型将趋向于是稀疏的,并且在{βi,i=1...M}中的许多参数将为零。那些与非零值βi相符的向量xi是已知的该模型的支持向量。支持向量的数量取决于可调整参数C的值、训练数据、以及该模型的适合度。在下面描述的内容中,体现了怎样利用核函数对该模型进行扩展,使其适应复杂的非线性函数。接着,将表示出ε-不敏感损失函数与l1常态缩减函数是相关的,它们在实质上是相同的,即依靠l1常态对一个稀疏参数组进行了一同选取。
为了模拟一个复杂函数,所述的复杂函数的变量之间存在可能的耦合,使用一个核函数来取代等式(17)中简单的内积(inner product),其中所述的核函数能够计算向量之间的更加复杂的相互作用。插入核函数之后,等式(17)得到的函数近似值具有下述形式:
(32)
其中将K(x,xo)定义为等于1。为了计算这些参数,使用与如上所述的最优化方法完全相同的方法,并且使用K(x,xi)来取代所有的xTxi项。如前所述,根据βi=αi +-αi -计算参数组,这种计算是基于发现了其最大化的受制于如上所述的相同的约束的论点。
对于如上所述的支持向量机结果,选取径向基核函数。
现在,为了说明l1常态的隐含使用:针对上述内容进行考虑,以取代尝试对等式(17)进行最优化,从下述最优化入手:
(34)
其中所述的l1缩减已经被明确的用来约束β的值以及上述数据的拟合误差,不同于使用训练数据的离散样本进行定义,其通过被模拟的假设函数的域来定义:βi=αi +-αi -;αi +,αi -≥0,αi +αi -≥0,i=1...N.。之后,可以通过下式重新进行最优化:
其受制于下述约束
αi +,αi -≥0
(36)
αi +αi -=0
(37)
尽管这个解具有不同的约束,虽然如此,但其将于ε-不敏感损失函数的解相符,当支持向量方法中的值C被选取的足够大并且约束0≤αi +,αi -≤C可以被容易的转变为约束(21)以及(22)并且基底函数之一是常数时,正如我们所说的情形中的等式(17)所示。在这种情况中,不需要支持向量方法中使用到的额外的约束 需要注意的是,约束(25)已经暗含在等式(15)中,因为约束(8)以及约束(9)是不能够同时起效的,因此,拉格朗日乘子之一αi +或者αi -应当是松弛的(slack),或者是零。
在这类情形中,我们能够发现,当暗含使用了l1缩减函数方法时,ε-不敏感损失函数可以达到稀疏函数的近似值。
多因子表型预测的实施例:模拟人类免疫缺陷病毒-1(HIV-1)的药物应答
用于对补救性抗逆转录病毒药物(ART)的表型结果进行预测的当前的方法没有表现出良好的预测能力,这很大程度上是因为缺乏统计学显著的结果数据,以及药物治疗方案间的许多不同的改变以及遗传变异。这一领域急切的需要对多种异源数据组进行整合,并且对药物应答预测进行加强。
这里公开的模型使用到了来自于斯坦福人类免疫缺陷病毒数据库逆转录酶以及蛋白酶药物抗性数据库(Stanford HIVdb RT and Protease Drug Resistance Database)中的数据,用来进行训练以及检验。这一数据由6644个人类免疫缺陷病毒-1病毒的体外表型检测结果所组成,在这些被检验病毒中,已经对逆转录酶或者蛋白酶编码的片段进行了测序。对十种逆转录酶抑制剂(RTI)以及七种蛋白酶抑制剂(PI)进行了检验。其中所述的逆转录酶抑制剂包括拉米夫定(3TC),阿巴卡韦(ABC),齐拉夫定(AZT),司他夫定(D4T),扎西他滨(DDC),去羟基苷(DDI),德拉维拉丁(DLV),依法韦恩茨(EFV),耐维拉平(NVP)以及特洛福韦(TDF)。所述的蛋白酶抑制剂包括ampranavir(APV),阿扎纳瓦(ATV),耐非纳瓦(NFV),瑞托纳瓦(RTV),赛科纳瓦(SQV),洛匹纳瓦(LPV)以及英迪纳瓦(IDV)。
对于每种药物而言,其数据被构建成下述成对的形式(xi,yi),i=1...N,其中N表示构成训练数据的样本的数量,yi是测量得到的药物的抗性倍数(fold resisitance)(或者表型),而xi是加上了一个常数项的变异向量,xi=[1xi1,xi2...xiM]T,其中M表示在相关的酶上存在的可能的变异的数量。假定当在ith样本中存在mth变异时,则元素xim=1,否则,将设定xim=0。同时使用一个密码子位点以及一个取代的氨基酸来定义每个变异。未对氨基酸序列产生影响的变异忽略不计。需要注意的是,对于每种药物而言,仅仅当多于1%的样本中存在变异时,其才能被包括在用于模型的一组可能的预测因素中,因为与抗性相关的变异不可能发生的如此不频繁。测量值yi表示的是与野生型相比,变异病毒对药物的抗性倍数。具体的,yi是变异病毒的IC50(药物将病毒复制抑制50%所需的浓度)的比例的对数,,与野生型病毒的IC50做比较。其目的在于为每种药物建立一个模型,该模型能够通过xi准确的预测出yi。为了对数据进行分批式最优化,利用M+1矩阵将自变量堆叠(stack)为N,X=[x1,x2...xN]T,并且将所有的观察值堆叠为一个向量y=[y1,y2...yN]T。
使用交叉验证方法对每种运算法则的表现进行测量。对于每种药物来说,可以通过该模型预测的表型应答与实际测量的检验数据的体外表型应答的比值来计算一阶相关系数(first-order correlation coefficient)R。
(38)
其中向量
表示的是对表型y的预测值,y表示的是向量y中的所有元素的平均值,而
表示的是所有元素的向量。对于每种药物以及每个方法而言,这些数据以9∶1的比例分别被随机的进行划分,分别用来进行训练以及检验。在一个实施例中,进行了十种不同的划分,目的是为了在不存在训练数据与检验数据的任何交迭的情况下计算向量
以及R。之后,可以将这个完整的过程重复进行十次,从而得出十个不同的R值。取上述十个不同R值的平均数,从而生成R的报告值(R reported)。计算每个模型的R值标准偏差,对多于十次的不同试验进行测量,以确保这些模型是在一种统计学显著的水平下进行比较的。
表11表示的是针对蛋白酶抑制剂药物的上述提及的模型的结果;表12表示的是针对十个逆转录酶抑制剂药物的模型的结果。这些结果是按照相关系数R来陈列的,取的是十次划分的训练数据以及检验数据的平均值。该表还列出了通过取样方差(sample variance)而计算出的R值平均值的标准偏差估算值。在最后一行列出了每种药物可以利用的样本的数量。为了提高平均性能,在上述方法中检测了:i)RR-岭回归,ii)DT-决策树;iii)NN-神经网络,iv)PCA-首要组分分析,v)SS-逐步筛选,vi)SVM_L-带有线性核函数的支持向量机,vii)LASSO-最小的绝对缩减和变量选择算子,以及viii)SVM-带有径向基核函数的支持向量机。在表11以及表12的最后一列中列出的信息在附图21中有所描述。附图21中的圆形表示的是针对每种蛋白酶抑制剂的十次不同试验得到的相关系数R的平均值,其中使用的是七种不同的蛋白酶抑制剂的平均值。附图21中的菱形表示的是针对每个逆转录酶抑制剂的十次不同试验得到的相关系数R的平均值,其中使用的是十种不同的逆转录酶抑制剂的平均值。该附图中还表示出了标准偏差的误差图。
当模拟技术中涉及调整参数时,对这些调整参数进行调整,使该技术取得最佳表现,所述的最佳表现是通过交叉验证来测量的,其中对参数的调整使用的是格点搜索法(grid search approach)。在所有的情况中,格点量子化都足够精细,所述格点中的最佳表现参数在实践中不能与给定数据的最佳参数进行区别,因为在预测中存在的差别是因为格点量子化在实验性底噪(noise floor)之下。
尽管这些数据具有明显的趋向性,但应当注意的是,由于样本的数量不同,各种药物间的潜在的遗传预测因素的相互作用以及数据的其他特性也不同,每种运算法则得出的R值因为药物的不同而不同。可以通过对表11中不同的药物纵列(第3列至第9列)以及表12中的不同的药物纵列(第3列至第12列)进行研究来发现这种变化。
在上述所有的方法中,支持向量机(SVM)的表现最佳,其些许的胜过最小的绝对缩减和变量选择算子(LASSO)(对于逆转录酶抑制剂而言,P<0.001;对于蛋白酶抑制剂而言,P=0.18)。经过ε-不敏感损失函数训练的支持向量机的表现明显的优于先前报道的基于支持向量机的方法。使用非线性核函数的支持向量机的表现优于使用线性核函数的支持向量机SVM_L,其中,使用非线性核函数的支持向量机也经过了ε-不敏感损失函数的训练(对于逆转录酶抑制剂而言,P=0.003;对于蛋白酶抑制剂而言,P<0.001)。支持向量机的表现明显的优于其他的非线性技术,其中所述的非线性技术使用了神经网络并且没有建立凸集合价值函数(对于逆转录酶抑制剂以及蛋白酶抑制剂而言,P<0.001)。最小的绝对缩减和变量选择算子技术使用凸集合价值函数以及连续的子集选取对线性回归模型进行训练,其表现明显的优于逐步筛选(SS)技术(对于逆转录酶抑制剂以及蛋白酶抑制剂而言,P<0.001)。最佳的五种方法,即逐步筛选(SS),首要组分分析(PCA),带有线性核函数的支持向量机(SVM_L),最小的绝对缩减和变量选择算子(LASSO),带有径向基核函数的支持向量机(SVM_R),都趋向于生成稀疏模型,或者趋向于生成具有有限数量的非零参数的模型。
为了对被选取为预测因素的变异子集进行说明,这里公开的某些实施方式着眼于具有次佳表现的模型,即最小的绝对缩减和变量选择算子(LASSO),其生成的线性回归模型不同于支持向量机,并不试图仿效在预测因素之间产生的非线性耦合或者逻辑耦合。因此,该模型能够直接的表明选取了多少个预测因素。表13表示的是由最小的绝对缩减和变量选择算子选取的变异的数量,这些被选取的变异与变异的数量(表13,第3行),样本的总数(表13,第2行)一同被用来作为对于每种蛋白酶抑制剂药物的预测因素(表13,第4行),对每个模型进行训练。在一个同样的表格中表示出了针对逆转录酶抑制剂的情况(表14,与表13具有相同项的行)。
这些被选取的变异同样能够加强对药物抗性的引起原因的理解。附图22、附图23以及附图24表示的是由最小的绝对缩减和变量选择算子选取的参数值,分别用于预测对蛋白酶抑制剂的应答,对核苷类逆转录酶抑制剂(NRTIs)的应答,以及对非核苷类逆转录酶抑制剂(NNRTIs)的应答。在这些附图中,每一行代表一种药物;每一列代表一个变异。对于蛋白酶抑制剂药物来说,相关的变异发生在蛋白酶上,而对于核苷类逆转录酶抑制剂药物以及非核苷类逆转录酶抑制剂药物来说,相关的变异发生在逆转录酶上。每个方块中的阴影表示的是该种药物中与该变异相关的参数的值。如右侧的彩条中所表示的(分别为2201,2301以及2401),阴影颜色较深的预测因素与增加的抗性相关;阴影颜色较浅的参数与增加的易感性相关。将这些变异从左至右进行排列,是为了减少相关参数的平均值的量级(magnitude)。所述的相关参数取的是所有行的平均数,或者是属于一类的所有药物的平均数。图中表示出了与四十个最大参数量级相关的突变。值得注意的是,对于一个具体的变异来说,或者对于一个具体的列来说,所有行中的参数值相差很大,或者属于同一类的不同药物中的参数值相差很大。
在岭回归(RR)、决策树(DT)、神经网络(NN)、以及逐步筛选(SS)这些运算法则中,不针对模型中的所有遗传变异进行训练,而是仅仅针对在某些位点处发生的变异的子集对模型进行训练,其中所述的位点处发生的变异被认为是影响抗性,这是由美国卫生和人类服务部(Department of Health and Human Services)(DHHS)确定的。已经发现,减少自变量的数量能够提高这些运算法则的表现。在使用线性核函数的支持向量机运算法则中,仅仅使用美国卫生和人类服务部(DHHS)确定的变异子集,能够针对逆转录酶抑制剂而取得最佳的表现,而针对所有的变异对模型进行训练,能够针对蛋白酶抑制剂而取得最佳的表现。对于其他所有的运算法则而言,通过针对所有的变异对模型进行训练,可以取得最佳的总体表现。
附图22、附图23以及附图24中表示的是由最小的绝对缩减和变量选择算子选取的作为预测因素的变异组,这些变异组通常不与由美国卫生和人类服务部(DHHS)确定的影响抗性的基因座相关,这些变异组是:对于蛋白酶抑制剂而言-19P,91S,67F,4S,37C,11I,14Z;对于核苷类逆转录酶抑制剂而言-68G,203D,245T,208Y,218E,208H,35I,11K,40F,281K;对于非核苷类逆转录酶抑制剂而言-139R,317A,35M,102R,241L,322T,379G,292I,294T,211T,142V。值得注意的是,在某些情况中,例如在最小的绝对缩减和变量选择算子以及支持向量机中,与仅针对模型中那些由美国卫生和人类服务部(DHHS)确定的影响抗性的基因座的情形相比(R=81.72,标准偏差=0.18),对于某类特定的药物例如洛匹纳瓦(LPV)来说,针对模型中所有的变异进行训练(R=86.78,标准偏差=0.17)将显著的提高其表现(P<0.001)。这说明除了那些由美国卫生和人类服务部(DHHS)确定的变异之外,其他的变异也可能在药物的抗性中起到作用。
这里已经证明,为了对表型预测模型进行训练以使其能够进行准确的归纳概括,使用凸集合最佳技术能够完成对稀疏参数组的连续子集选取。最小的绝对缩减和变量选择算子利用l1常态缩减函数来生成一个线性回归参数的稀疏组。带有径向基核函数并且使用ε-不敏感损失函数进行训练的支持向量机能够生成稀疏的非线性模型。这些技术的出众表现可以解释为它们的价值函数的凸集合性,以及它们趋向于生成稀疏模型的性质。凸集合性确保了当存在许多可能的预测因素时,我们可以针对一组特定的训练数据组而找到全局最优参数。稀疏模型趋向于进行很好的归纳概括,特别是在欠定数据或者不适定数据为主要的遗传数据的情况中时。可以将l1常态认为是最具选择性的凸集合函数。使用一个选择性缩减函数对一个稀疏参数组进行选取时,其运用的法则与奥卡姆剃刀法则(Occam’s Razor)相类似:当存在许多可能的理论来解释观察到的数据时,最简单的理论最可能是正确的。同时使用了l2缩减函数以及ε-不敏感损失函数的支持向量机能够趋向于产生与l1常态缩减函数应用于与支持向量相关联的参数中所起到的明确的作用而相类似的作用。
当自变量的数量很大,并且其数据是欠定的或者是不适定的时候,使用l1缩减函数的技术通常能够准确的进行归纳概括。因此,向模型中加入自变量的非线性组合或者逻辑组合成为可能,并且期望这些作为优秀的预测因素的组合能够在训练中被选取。使用了非线性核函数的支持向量机能够模拟自变量之间的相互作用,其中所述的非线性核函数是例如径向基函数,其表现明显的优于线性核函数。因此,在不改变这里公开的基本观点的情况下,可以通过向模型中加入自变量的逻辑组合来提高最小的绝对缩减和变量选择算子的表现。逻辑项可以来自于那些由决策树生成的逻辑项、来自于由专家法则所描述的逻辑相互作用、来自于逻辑回归技术、或者甚至是来自于逻辑项的一组随机排列。最小的绝对缩减和变量选择算子的一个优点在于其得到的模型很容易被翻译(interpret),因为其中的参数直接与自变量相结合,或者与包括自变量的表达式相结合,而不是与支持向量相结合。最小的绝对缩减和变量选择算子对于模型中大量的自变量的强固性(robustness)是由于l1常态的选择特性以及该方法的凸集合性的原因。
还存在其他的技术,在这些技术中使用缩减函数比使用l1常态更具选择性。例如,对数缩减回归使用到的缩减函数是来自于编码理论的,其对存在于模型参数组中的信息的量进行测量。这项技术使用对数函数作为缩减函数来替代l1常态,因而是非凸集合性的。当提供一种理论上有兴趣的方法来寻找一个稀疏参数组时,其补偿参数(penalty function)的非凸集合性表明了:与最小的绝对缩减和变量选择算子相比,该方法对于相应的回归的计算机计算仍然是较难处理的,并且对于一个给定的数据而言,大量的预测因素组可能仅仅产生一个局部最小解,而不是一个全局最小解。
这里描述的技术可以被应用于很大范围的表型预测问题当中,用于建立线性回归模型以及非线性回归模型。当潜在的遗传预测因素的数量大于测量结果的数量时,这类技术是格外适用的。
通过将遗传自变量映射(Mapping)至不同的空间中的方法对回归模型进行简化
如上所述,值得注意的是,对存在遗传标记的复杂组合的情况进行考虑,为了对分析进行简化,可能将单核苷酸多态性变量投射(project)至另一个变量空间中。这个变量空间可以代表已知的变异类型,例如在HapMap(人类基因组单体型图)计划中描述的群(clusters)或者库(bins)。换句话说,除了使用向量xi表示如上所述的特定的单核苷酸多态性变异,还可以表示出该个体是否落入特定的人类基因组单体型图群或者库中。例如,使用如上所述的表示符号,假设具有一个向量xi=[xi1,xi2...xiB]T,其中B表示相关的人类基因组单体型图库的数量。我们可以设定元素xib=1,当上述个体的单核苷酸多态性类型落入bth库中时,否则,将其设定为0。或者,如果个体的单核苷酸多态性与特定的库之间的交迭是不完善的(incomplete),并且不希望简单的将所述的个体归结为“其他”种类,那么我们可以设定每个xib,使其等于该单核苷酸多态性类型与该库b间的交迭的片段。在不改变这里公开的观点的情况下,许多其他的技术都可能被用来制定上述回归问题。
通过交叉验证对模型进行选取,用于进行结果的预测
在先于本发明之前的文献中,描述了各种表型预测技术,包括专家法则,列联表,线性回归以及非线性回归。现在,描述一种从一组模型技术中进行选择的常规方法,这种选择针对的是一个具体的宿主,该方法基于训练数据的使用来进行选择,目的在于选取能够最佳的模拟特定的无条件结果或者非无条件结果的模型技术。附图25提供了一个对于该体系的说明性的流程图。附图25中描述的过程是选择最佳模型的一个常规方法,所述的最佳模型能够给出特定的患者可以利用的数据,能够对表型进行模拟,并且给出一组检验数据以及训练数据,并且这一过程不依赖于具体的模型技术。在一个优选的实施方式中,可以使用的模型技术组包括专家法则,列联表,使用最小的绝对缩减和变量选择算子或者简单的最小平方进行训练的线性回归模型,以及使用支持向量机的非线性回归模型,其中所述的使用最小的绝对缩减和变量选择算子或者简单的最小平方进行训练的线性回归模型中不存在不适定的数据。
上述过程通过对需要被模拟的特定的宿主以及特定的依变量(DV)进行选取,从而开始2501步骤,或者——如果这是一个无条件变量——可以对该变量的概率进行模拟。随后,这一体系确定2502中与该宿主的记录相关的自变量(IV)组,所述的自变量组可能与模拟所述的依变量的结果是相关的。该体系中的人类使用者也可以对自变量的子集进行选取,选取那些使用者认为可能与模型相关的自变量子集。之后,该体系对2503a进行检查,查看模型是否已经进行了训练,并且是否已经为需要被模拟的给定的自变量的组合以及给定的依变量进行了选取。如果已经进行了训练以及选取,并且用于训练的数据以及用于检验已经建立的模型的数据不是过时的,该体系将直接使用该模型生成一个预测2519。否则,该体系将从数据库中提取其他所有的记录,其中所述的记录具有特定的感兴趣的依变量,并且其具有或者不具有与特定的感兴趣的宿主相同的自变量组。通过这种做法,该体系确定了2503b数据是否可以用于训练模型以及检验模型。如果答案是否定的,该体系对2515进行检查,查看是否具有任何可以利用的专家法则,用来基于该宿主可以利用的自变量的子集对结果进行预测。如果不存在可以利用的专家法则,那么该体系从2504处退出,并且表明无法做出有效的预测。如果存在一种或者一种以上可以利用的专家法则,那么该体系将选取2505一个专家法则子集,所选取的专家法则子集最相称于特定宿主的数据。在一个优选的实施方式中,对应用于宿主的专家法则的选取将基于该专家法则的估计值的置信度水平来进行。如果没有这样的置信度估算可以利用,可以按照它们的特异性水平对这些专家法则进行分类,即通过计算这些专家法则在预测过程中使用到的该目标宿主可以利用的自变量的数量进行分类。随后,被选取出的专家法则的子集可以被用来进行预测2506。
如果在2503b处确定了该数据是可以利用的,该体系将检查2516,确认在检验数据以及训练数据中是否存在任何任何数据的丢失。换句话说,对于所有那些包含相关的依变量的记录,该体系都将对其进行检查,查看所有这些记录是否确实具有与目标患者可利用的相同的自变量组,其中这些自变量组可能是模型中潜在的预测因素。通常,这一答案将会是“否”,因为不同的患者将利用不同的信息。如果存在丢失的数据,该体系将进行一个步骤,用来找到应当被用来对该宿主进行最可能的预测的自变量的子集。这一步骤是耗费时间的,因为它包括了多轮的模型训练以及交叉验证。因此,在这个过程中,第一步骤就是根据可以利用的计算时间来减少2507自变量组,使其减少到一个可以操作的数量。在一个优选的实施方式中,根据自变量中含有的针对某些同样具有可利用的依变量的部分宿主的数据来进行自变量组的减少。可以使用本领域已知的其他技术来进一步的对自变量组进行减少,例如使用逐步筛选技术,该技术假定一个简单的线性回归模型,并且依据自变量与模型误差的相关性程度来选取自变量。之后,该体系进入一种循环,在该循环中,对剩余的自变量的每种组合进行了检查。在一个优选的实施方式中,对每个自变量以及依变量的下述状态进行了考虑:在模型中可以包括或者不包括每个自变量,并且对于所有的宿主来说,如果一个自变量或者一个依变量的数值数据是正的,该数据可以或者不必进行取对数的预处理。对于每种特定的包含/不包含与自变量和依变量的预处理的组合,应用一组模拟技术2510。
大多数的模拟技术都具有某种调整参数,调整参数可以基于格点搜索法通过交叉验证使用检验数据进行最优化或者调整。例如,对于如上所述的最小的绝对缩减和变量选择算子而言,需要为可变参数λ而探测许多值。对于每个λ值而言,可以对回归参数进行训练,并且模型的预测值将与检验数据的测量值进行比较。类似的,对于如上所述的支持向量机方法而言,需要使用格点搜索法进行最优化的调整参数包括C参数、ε参数以及描述核函数的特征的其他参数。对于基于列联表的技术而言,其中的调整参数可以与最高的标准偏差相一致,其中所述的最高标准偏差可以是一个列联表模型中认可的,尽管如上所述,对于给定的宿主而言,尽可能的将所述的列联建立的更加具体。
可以利用许多不同的计量学方法,用来将模型的预测值与检验数据进行比较,其目的在于对调整参数进行最优化,并且在模型中进行选取。在一个优选的实施方式中,使用到了误差的标准偏差。在另外一个实施方式中,可以使用预测结果与测量结果间的相关系数R。在逻辑回归或者列联表的方法中,还可以使用最大后验概率,即给定的检验数据组假定每个检验结果的可能性的模型预测值的概率。不论使用了怎样的计量学方法,对调整参数的值进行选取,使其对计量学值进行最优化,例如如果使用预测误差的标准偏差作为检验计量值时,则需要对预测误差的标准偏差进行最小化。由于模型的训练以及交叉验证是一个缓慢的过程,在2510过程中用来定义需要进行检查的不同的调整参数的格点分布的很粗糙,根据可以利用的时间的多少,因而仅仅能够获得最佳模型以及最佳调整参数的大体想法。
一旦使用这种方式2511对不同的自变量/依变量的所有组合进行了检查,该体系将对自变量/依变量的组合、模型以及调整参数进行选取,选取那些能够生成检验计量学最佳值的组合、模型以及调整参数。值得注意的是,如果其中不存在丢失数据,那么该体系将省略检查所有的自变量/依变量的组合的步骤。相反的,该体系将检查不同的模型技术以及调整参数2508,并且将对模拟方法以及调整参数的设置进行选取,选取那些能够使检验计量值最大化的模型以及设置。之后,该体系使用更加细微的空间格点对最佳回归模型进行精确的调整,并且为每组调整参数值确定其与检验数据的相关性。被选取的调整参数组应当生成检验计量学的最佳值。之后,该体系确定2518所述的检验计量值例如预测误差的标准偏差是否是低的或者是否低于一个选取的阈值,这样一来,得到的预测值被认为是有效的。例如,在一个实施方式中,对于一个预测值而言,如果其相关系数R>0.5,则该预测值被认为是有效的。如果作为结果的检验计量值不满足阈值,那么将无法进行预测2517。如果该检验计量值满足了必需的阈值,则可以得到一个表型预测值,以及用于该预测的自变量的组合,以及该模型与所述的检验数据之间形成的相关系数。
在癌症群组中使用丢失数据通过交叉验证进行模型选取方法的说明
为了对这个方面进行证明,我们着重利用与结肠癌相关的遗传数据以及表型数据,这些数据可以从遗传药理学和基因组药理学数据库(PharmGKB)中得到,所述的遗传药理学和基因组药理学数据库是健康基因组药理学研究网络国立研究中心(the National Institute of Health’s Pharmacogenomic Research Network)的一个组成部分,这些数据的任务是探索个体遗传变异对不同的药物应答具有怎样的影响。对于这个数据集而言,其主要的挑战是丢失信息。在理想状态下,我们愿意使用如上所述的回归技术,从特定患者可以获得的所有自变量中为模型自动选取一个自变量子集。然而,这限制了可以从其他患者处获得的、用于对模型进行训练以及检验的数据的数量。因此,对于含有足够少的自变量的数据集而言,在自变量的所有可能的子集中进行寻找成为可能。如上所述,对于每个子集而言,我们可以从该组患者中提取已经被测量了的所需的结果,因而可以获得相关的自变量组。如上所述,还可以寻找对包括在内的自变量进行预处理的可能的方式的空间,例如对呈正值的自变量取对数。在每种包括在内的自变量的组合以及自变量的预处理技术中,使用检验数据通过交叉验证对模型进行训练以及检验。选取那个与检验数据发生了最佳的交叉验证的模型。一旦针对给定设置的自变量组建立了一个模型,该模型即可以被应用于新的患者数据中,只要所述的新的患者数据服从于同样的自变量设置即可,不需要进行详尽的模型研究。
这项技术被用来预测结肠直肠肿瘤药物伊立替康的临床副作用。在接受伊立替康治疗的癌症患者的体内通常会观察到严重的毒性。其中包括了描述伊立替康的药物动力学与副作用间的关系的数据,以及编码伊立替康代谢酶以及推定相关的转运蛋白的基因的等位基因变异。对患者进行基因定型分析,对基因中的变异进行基因定型,其中所述的基因编码MDR1 P-糖蛋白(ABCB1),多种药物抗性相关蛋白MRP-1(ABCC1)以及MRP-2(ABCC2),乳腺癌抗性蛋白(ABCG2),细胞色素P450同功酶(CYP3A4,CYP3A5),羧酸酯酶(CES1,CES2),尿苷二磷酸葡萄糖醛酸转移酶(UGT1A1,UGT1A9),以及肝脏转录因子TCF1。该项研究中与遗传序列数据相关的表型数据在表15中有所描述。
附图26描述了对于使用伊立替康的结肠癌治疗的预测结果的模型,给定了可以利用的遗传药理学和基因组药理学数据库(PharmGKB)数据,该数据是通过基因组药理学翻译引擎翻译提供的。在附图26中,该模型表示出了相关的遗传基因座(2601),在这种情况中,指示剂利用0-24小时内的CPT-11浓度曲线下区域(AUC)的对数(2602)以及0-24小时内的SN-38浓度曲线下区域(AUC)的对数来预测从第12天至第14天的嗜中性粒细胞绝对计数的最低点(Nadir of Absolute Neutrophil Count)的对数(2604)。使用检验数据对上述模型进行交叉验证,得到R=64%的相关系数(2605)。表示出了模型预测值的经验标准偏差(2606),该标准偏差依靠用来对模型进行训练的结果的柱状图(2607)形成叠加。这些统计学信息可以被用来进行一个通用的治疗决策,例如完全放弃伊立替康治疗或者施用第二种药物,例如粒细胞结肠刺激因子,用来防止较低的嗜中性粒细胞绝对计数(ANC)以及由此引起的感染。
加强的诊断报告
在疾病治疗的情况中,生成的表型数据多为临床医师所使用,临床医师利用这些数据帮助选择治疗方案。在一个方面,可以将表型的预测值融入或者安排入临床医师或者患者的报告中。在另外一个方面,可以将这里公开的体系以及方法作为一个更大的体系的一部分(参见附图27)来使用,在该体系中,诊断实验室2703对来自于实验室检验的数据2701以及医学报告的数据2702进行验证,并且将其送至一个数据中心2704,在该数据中心,这些数据被整合为一个标准的本体论(ontology),使用本发明公开的方法对其进行分析,可以生成一个加强了的诊断报告2705,并且将该报告送至内科医师2706处。
在一种可能的情况中,生成的报告涉及的是使用伊立替康的结肠癌患者的临床结果的预测。该报告可能对例如治疗的禁忌症候、剂量进度表、副作用曲线之类的观点纳入了考虑范围。这类副作用的范例包括骨髓抑制以及晚期发作的腹泻(late-onset diarrhea),上述两种症状是伊立替康治疗中常见的两种剂量限制性副作用,需要进行紧急的医疗救护。除此之外,严重的嗜中性白血球减少症以及严重的腹泻分别影响到28%以及31%的患者。某些尿苷二磷酸葡萄糖醛酸转移酶(UGT1A1)等位基因,肝功能检测,吉尔伯特氏综合症(Gilbert’s Syndrome)的既往病史,以及对于能够诱发细胞色素P450的患者的药物治疗的识别,例如抗痉挛药物以及某些抗呕吐药物,都可以作为指示剂,用来确保对伊立替康剂量的调整。
附图28是一个针对使用伊立替康的结肠直肠肿瘤治疗的加强了的报告的实体模型(mock-up),其中使用到了表型的预测。在治疗之前,该报告中对患者的肿瘤阶段、既往病史、当前的药物治疗、以及尿苷二磷酸葡萄糖醛酸转移酶(UGT1A1)的基因型进行了考虑,借此推荐药物剂量。在首次药物剂量后的大致一个数据中,该报告包括了对于大致两周时间内期望的患者的嗜中性粒细胞绝对计数的最低点进行的预测,这种预测是基于尿苷二磷酸葡萄糖醛酸转移酶(UGT1A1)基因中的变异以及从患者血液中测量到的代谢物(例如,SN-38,CPT-11)而进行的。根据这种预测,医生可以决定是否给患者施用结肠刺激因子药物,或者改变伊立替康的剂量。所述患者的血球计数以及腹泻程度也同样被进行了监控。该模型还提供了用于推荐药物剂量的数据来源以及推荐理由。
多方面的组合
如前所述,在给定了本发明公开内容的益处的条件下,其他的方面、特征以及实施方式也可以用来完成这里公开的一种或者一种以上的方法以及体系。下面是一些实施例的简短列表,这些实施例描述的是本发明公开的各种不同的方面通过各种方式进行组合所产生的情形。很重要的一点是,需要注意,这个列表并不意在进行全面的概括,本发明的某些方面、特征以及实施方式的许多其他的组合也是可能的。
在一个实施例中,以一种方式利用各种不同的基因型测量技术,使得每项技术的值都得到最优化。例如,当信号较低时,实验室可以利用比较昂贵但是能够得到高质量的数据的技术,例如AppliedBioscience’s Taqman检测,来测量目标DNA,并且使用比较便宜但是需要大量的遗传物质才能获得高质量的数据的技术,例如Affymetrix’s 500K基因芯片,或者分子倒置探针,来测量双亲DNA。
在另外一个实施例中,一对夫妇接受了体外受精治疗,从妻子处获得了卵细胞,并且使用从丈夫处获得的精子对该卵细胞进行受精,生成了八个存活的晶胚。从每个晶胚中获得一个胚泡,并且使用Taqman基因定型检测对每个胚泡中的遗传数据进行测量。同时,从双亲中分别提取组织,并且使用分子倒置探针测量该组织中的二倍体数据。同样使用分子倒置探针对来自于丈夫的精子之一的单倍体数据进行测量,以及对来自于妻子的卵细胞之一的单倍体数据进行测量。使用双亲的遗传数据对八个胚泡中的单核苷酸多态性数据进行清除。随后,上述被清除过的遗传数据可以被用来进行预测,所述的预测是关于晶胚的可能的表型。选取具有最具希望的曲线的两个晶胚,并且将其植入母亲的子宫当中。
在另外一个实施例中,一名怀孕的妇女其丈夫具有家族性白痴病(Tay-Sachs disease)的家族病史,该妇女想知道她怀有的胎儿是否具有遗传易感性,但该妇女不希望接受羊水诊断,因为羊水诊断会带来严重的流产的危险。该妇女进行了抽血,从她的血液中分离出了胎儿的DNA,并且使用分子倒置探针对胎儿的DNA进行了分析。该妇女及其丈夫先前已经对它们的全基因组数据进行了分析,该分析可以通过计算机数据分析方法(in silico)来获得。医生可以利用双亲基因组的计算机数据分析方法知识以及这里公开的方法来进行胎儿DNA数据的清除,并且检查在胎儿的基因组中是否存在导致家族性白痴病的关键基因。
在另外一个实施例中,一名44岁的怀孕妇女担心她怀有的胎儿可能患有唐氏综合症(Downs Syndrome)。该妇女抗拒使用入侵性技术来进行产前诊断会产生流产的个人病史,因而她选择对其血液进行分析。健康医疗从业者也可以从母亲的血液样本中找到胎儿的细胞,并且使用这里公开的方法,结合这位妇女自己的遗传数据知识,对异倍体进行诊断。
在另外一个实施例中,一对夫妇接受了体外受精治疗;从妻子处获得了卵细胞,并且使用从丈夫处获得的精子对该卵细胞进行受精,生成了九个存活的晶胚。从每个晶胚中获得一个胚泡,并且使用Illumina磁珠检测对每个胚泡中的遗传数据进行测量。同时,从双亲中分别提取组织,并且使用分子倒置探针测量该组织中的二倍体数据。使用同样的方法对来自于丈夫的精子中的单倍体数据进行测量。由于不能从妻子处获取额外的卵细胞,因而从妻子自己的父亲以及母亲体内提取了大块二倍体组织样本,并且从其父亲体内提取了一个精子样本。使用分子倒置探针对上述所有的样本进行分析,并且使用这里公开的方法对妻子的染色体组进行遗传分析。随后,将分析得到的数据与丈夫的二倍体数据以及单倍体数据一同使用,对每个胚泡中的遗传数据进行高度准确的分析。根据表型的预测值,这对夫妇选择了三个晶胚以供植入体内。
在另外一个实施例中,一位赛马饲养者希望增加它的冠军赛马的后代马驹也能成为冠军的可能性。他安排了一匹选定的母马,使该马通过体外受精而得到受孕,并且利用种马以及母马的遗传数据对从存活的晶胚中测量到的遗传数据进行清除。清除过的胚胎遗传数据能够帮助饲养者找到相关的基因型-表型相关性,并且可以对植入的晶胚进行选择,选出最有可能生成期望的种马的晶胚。
在另外一个实施例中,一名怀孕的妇女想知道她怀有的胎儿是否具有患有任何严重疾病的倾向性。因为该妇女的丈夫已经去世,因此使用来自于其丈夫的兄弟以及其丈夫的父亲的单倍体数据以及二倍体数据来帮助对胎儿的遗传数据进行清除,其中所述的遗传数据是从胎儿血液取样的过程中收集到的胎儿细胞中测量到的。健康医疗从业者利用清除过的胎儿遗传数据提供出了一系列胎儿可能表现出来的表型,以及每个预测值的置信度。
在另外一个实施例中,一个羊水诊断实验室有时候需要与污染了的胎儿遗传数据进行斗争,其中所述的污染是由于较差的实验室技术而导致的。可以利用本发明公开的方法,使用母亲以及父亲的遗传数据对污染了的胎儿的遗传数据进行清除。我们可以设想下述情形:如果已知本发明公开的方法能够用来补偿增加了的DNA的污染率的时候,实验室就可以通过减少无菌操作来降低成本。
在另外一个实施例中,一名四十岁的妇女通过体外受精方式而获得受孕。她希望通过对晶胚进行筛选,选取那个/那些具有最小的患有遗传疾病的可能性、并且具有最大的被植入并且足月生产(carry to term)的可能性的晶胚。该名妇女所接受治疗的体外受精诊所从每个存活的晶胚中提取了一个胚泡,并且利用标准技术对其中的DNA进行了扩增,并且对主要的单核苷酸多态性进行了测定。随后,技术人员使用这里公开的方法对染色体不平衡进行了筛选,也找到并且清除了晶胚的遗传数据,从而对每个晶胚的表型素因(phenotypic predispositions)进行了预测。
在另外一个实施例中,一名怀孕的妇女进行了羊水诊断,并且使用这里公开的方法以及该妇女的血液样本中的胎儿细胞的遗传物质对异倍体以及其他染色体异常进行筛选。
在一个实施例中,一个非线性模型使用带有径向基核函数的支持向量机以及一个常态损失函数利用来自于一个人类成人的遗传数据以及表型数据预测患有早期发作的阿尔茨海默氏症(早老性痴呆症)的可能性,并且为可能的生活方式的改变以及锻炼方案的制定提供建议,其中所述的可能的生活方式的改变以及锻炼方案的制定能够延缓上述疾病的发作。
在另外一个实施例中,一个使用最小的绝对缩减和变量选择算子技术的线性模型利用一个患有肺癌的成年妇女的遗传数据以及表型数据,以及该肿瘤的遗传数据为该妇女的内科医师建立了一份报告,该报告预测了哪种药物对于延缓上述疾病的恶化是最有效的。
在另外一个实施例中,使用由遗传数据、表型数据以及临床数据组成的集合数据对许多模型进行了检验,所述的这些数据均来自于克罗恩氏病(Crohn’s disease)患者,随后,找出最准确的非线性回归模型,利用来自于一个成年男子的表型数据以及临床数据建立一份报告,该报告给出了某些营养补充方案的建议,这些方案具有缓解该患者的克罗恩氏病的症状的可能性。
在另外一个实施例中,一个利用列联表以及遗传信息的模型被用来对一个婴儿具有的可能的表型进行预测,其中所述的列联表是通过从Hapmap(人类基因组单体型图)计划中获得的数据而建立的,其中所述的遗传信息是从一个晶胚的胚泡中收集得到的,其中所述的婴儿是由上述的晶胚被植入之后发育而成的。
在另外一个实施例中,线性回归模型利用一个新生儿所感染的人类免疫缺陷病毒(HIV)菌株中的遗传信息为该新生儿的内科医师建立了一份报告,该报告建议了施用哪类抗逆转录病毒药物能够给这个新生儿带来生长至成人期的最大机会。
在另外一个实施例中,公布了一项新的研究,该研究揭示了心肌梗塞在中年妇女中的流行性与某些遗传标记以及表型标记之间的某些相关性。之后,这项研究促使使用一个非线性回归模型对中年妇女的数据的集合数据以及来自于被确诊的个体的遗传数据以及表型数据进行重新检查,其中所述的被确诊的个体的数据是该体系中已知的,之后,该模型将鉴别出那些最具有患心肌梗塞危险性的妇女,并且为其建立报告,该报告将被送至这些妇女各自的内科医师处,通知医师所存在的预测的危险。
在另外一个实施例中,使用来自于患有结肠癌的患者的集合数据对各种模型进行检验,包括曾经尝试过的各种药物干预。使用能够提供出最佳预测值的模型来鉴别最有可能从试验性的新药物中获益的患者,拥有该种新药物的公司将利用这些结果来帮助它们进行临床试验。
定义
SNP(单核苷酸多态性):趋向于表现个体与个体间的差异的一个染色体上的一个具体的基因座。
单核苷酸多态性分型(To call a SNP):对一个特定的碱基对的同一性进行询问(interrogate),考虑直接证据以及间接证据。
等位基因分型:单核苷酸多态性分型。
遗传数据的清除:利用相关个体的遗传数据以及这里公开的方法,对不完美的遗传数据进行选取,并且修正其中的某些误差或者全部的误差。
不完美的遗传数据:具有下述情形中的任何之一的遗传数据:等位基因脱失,不清楚的碱基对测量值,不正确的碱基对测量值,伪造信号,或者丢失的测量值。
置信度:分型的单核苷酸多态性、等位基因、或者等位基因组正确的表示出个体的真实遗传状态的统计学可能性。
多基因的(Multigenic):由多个基因或者等位基因所影响的。干扰遗传数据/噪声遗传数据:不完全的遗传数据,也叫做不完全的遗传数据。
未清除的遗传数据:在不使用任何方法对存在于原始遗传数据中的干扰(噪声)的存在进行修正的情况下测量的遗传数据;也叫做天然遗传数据(crude genetic data)。
直系亲属:母亲,父亲,儿子,或者女儿。
染色体区域:一个染色体片段,或者是一个全染色体。
双亲支持(Parental Support):有些时候本发明公开的清除遗传数据的方法所使用的名称。
染色体片段(Section of a chromosome):一个染色体的一部分(section),其大小可以从一个碱基对到整个染色体。
表格
表1
表2
表3
表4
表5
表6
表7
|
PS运算法则1 |
PS运算法则2 |
种群频 ph pd pe率 |
P1准确 P2准确性 性 |
P1准确 P2准确性 性 |
数据 0.95 0.95 0.95数据 0.75 0.75 0.75数据 0.25 0.25 0.25数据 0.5 0.9 0.9数据 0.9 0.5 0.5数据 0.6 0.8 0.8 |
0.834 0.8150.797 0.7690.711 0.6820.849 0.8380.792 0.8090.777 0.788 |
0.928 0.9310.819 0.8190.703 0.6870.866 0.8640.756 0.7520.835 0.828 |
一致性 0.95 0.95 0.95一致性 0.75 0.75 0.75一致性 0.25 0.25 0.25一致性 0.5 0.9 0.9一致性 0.9 0.5 0.5一致性 0.6 0.8 0.8 |
0.673 0.6310.549 0.4970.239 0.2490.601 0.6110.459 0.3910.544 0.511 |
0.898 0.9010.635 0.650.252 0.250.814 0.8180.449 0.4680.672 0.679 |
表8
Farrer(2005) |
|
|
|
D+ |
D- |
A+B+ |
56/151 |
25/206 |
A+B- |
26/151 |
69/206 |
A-B+ |
43/151 |
20/206 |
A-B- |
26/151 |
20/206 |
Alvarez(1999) |
D+ |
D- |
A+ |
161/350 |
176/400 |
A- |
189/350 |
224/400 |
Labert(1998) |
D+ |
D- |
B+ |
334/573 |
104/509 |
B- |
239/573 |
405/509 |
表9
表10
表11
表12
药物 |
APV |
阿扎纳瓦 |
耐非纳瓦 |
瑞托纳瓦 |
赛科纳瓦 |
洛匹纳瓦 |
英迪纳瓦 |
#样本 |
577 |
142 |
617 |
510 |
598 |
253 |
579 |
#变异 |
320 |
320 |
320 |
320 |
320 |
320 |
320 |
#预测因素 |
120 |
50 |
90 |
120 |
120 |
100 |
110 |
表13
药物 |
拉米夫定 |
阿巴卡韦 |
齐多夫定 |
司他夫定 |
扎西他滨 |
去羟基苷 |
德拉维拉丁 |
EFV |
耐维拉平 |
特洛福韦 |
#样本 |
366 |
364 |
363 |
366 |
319 |
364 |
395 |
384 |
401 |
96 |
#变异 |
791 |
791 |
791 |
791 |
791 |
791 |
791 |
791 |
791 |
791 |
#预测因素 |
50 |
70 |
100 |
110 |
80 |
100 |
170 |
120 |
120 |
60 |
表14
表15