文摘

DNA序列数据现在被用来研究人类的祖先历史。这种合并的现有方法推理用递推公式来计算概率的数据。这些方法在实际应用中非常有用,但计算复杂。这里我们首先调查的渐近行为这样的推理;结果表明,广泛的,估计合并的时间将一致的一个有限的限制。然后我们研究一个相对简单的计算方法来分析和说明如何使用它。

1。介绍

在过去的几十年里,已经取得了相当大的进展在群体遗传学领域。的一个主要目标是推断聚结时间下的人口研究,也就是说,推断出从他们最近的共同祖先(MRCA)及其分布基于观测数据。

在遗传学、合并的理论是一个回顾性的群体遗传学的痕迹中所有基因样本人口单一祖先副本由所有成员共享的人口。人口的合并的时间的时候他们最近的共同祖先。基因之间的继承关系,通常表示为一个基因家谱,类似于系统发育树。联合分析的目的是来推断的合并时间的样本 个人独立抽样从人口规模 根据他们的观察,DNA序列的多样性。与参数推理为独立同分布(iid)数据,为渐近极限可以方便地用于描述估计当数据规模很大,各种现有的研究表明,估计MRCA单位 代,还不清楚是否会随着数据样本量的增加集中没有绑定。相比之下,估计的突变速率相同的设置,估计是一致的和渐近正态的1),虽然速度慢得多 相比,的速度 i.i.d.数据。同时,不同于通常的参数,MRCA变化 的序列。这提示我们调查联合估计的渐近行为。我们想知道是否这样的估计将渐近一致和如果它在何种意义上。调节隔离总数的网站,我们发现这种估计收敛或不后的一些非负有限的限制意味着,根据行为突变的数量在所有的分支的树由观测数据。这个问题的同时,分析这种类型的数据通常是计算广泛而复杂;我们研究这个问题的相对简单的模拟方法。我们首先研究这种方法的渐近行为3,然后描述和说明我们的方法这个问题的部分4

在聚结推理,线粒体DNA (mtDNA)数据扮演重要的角色。线粒体是为数不多的现有基因在细胞核外,对哺乳动物来说,这仅仅是母体遗传来的。人类mtDNA是一个双链分子序列长度约16500个碱基对。在细胞外核,众所周知,mtDNA的突变率是大约10倍的核基因、线粒体的一个部分,其控制地区,高突变率甚至一个秩序。简单的继承模式和高可变性使mtDNA人类进化历史研究的一个重要来源。DNA链上的每个网站的四个基地、C、G和t .分子进化,基因突变发生在基地的形式替换。之间的变化嘌呤(G)或嘧啶(C、T)过渡;嘌呤和嘧啶之间是颠换。前类型替换比后者更为常见。

我们专注于线粒体的控制区域数据在格里菲斯和Tavare [2),这是部分的数据在病房等。3]。他们从一段控制区域,有352个碱基对(网站),其中159是嘧啶嘌呤网站和193网站。这个数据包含63序列采样从北美印第安部落,Nuu-Chah-Nulth,从温哥华岛。消除序列与多个突变后在一些单一的网站,这样的假设一个突变每个站点,其余数据有55个序列,与14个不同的序列数据(称为血统)。网站,并不是所有的观察序列有相同的基础是一个隔离的网站。整个序列很长,但是只有隔离网站信息的分析;其他网站都忽略了。提到数据18隔离网站和提出了表1,每个家族的频率(或多重性)。

合并是一个随机样本的系谱树模型 从大量的DNA序列。样本大小的树的一个例子 图中给出了1

这个话题的更详细的评论,看到哈德逊(4]和唐纳利Tavare [5]。

在聚结推理有以下之一。

基本假设。人口规模 很大,不变对许多代过去,是已知的,或从其他来源可以估计;人口是一个随机样本的数据;的数量在每一代生育之前Wright-Fisher模型(自常数大小,人口死亡的人数也遵循类似的模式);在任何核苷酸突变(替换)网站只会发生一次的祖先,是不可逆转的;突变发生在不同的时间间隔是独立的;突变发生的时间点是毒药和速度分布 定义,独立在家谱树的每个分支, 是已知的,或从其他方法可以估计或来源。

聚结时间的推断 样本的人口规模 有两个步骤。第一步是建模的分布 没有任何数据,predata分布;然后在第二步中,更新predata分布,利用观测数据,postdata分布,在此基础上进行正式的推理。predata分布是由金曼开创6,7];他表明,在单位的时间 代, 在哪里 是独立的等待时间。 是时候从 样本的共同祖先 共同的祖先。可以找到一个快速参考Tavare [8]。在这里 分布是指数 , 。的 年代可以以图形方式表示作为一个联合树如图1;然后 树的高度。树的长度定义为 然后(金曼) 年的时间单位转换的关系 ,在那里 每一代的平均年,通常是20 - 25。在这里我们看到,作为一个初步分析没有观察到的数据,合并的时间随机样本的大小 从人口规模 大概是 一代又一代,只要 是比较大的。因此,从分组人口样本的合并的时间大致相同的人口(只要样本容量比较大的)。这种现象进一步研究了沃特森(9显示的) 在哪里 祖先的数量,在哪里 代以前,人口的规模 从数据样本的大小 是画的。这里的样品必须是随机从人口;否则可能不可靠的结果。例如,样本的大小 是从一个族群的尺寸吗 从人口规模 ;然后由(3),合并的predata估计时间 这个示例是约 代,但也差不多 代以来样本也从整个人口。抽样方案的悖论的出现。如果样本来自的族群大小 ,一个只能使用 随着时间尺度,不是 ,因为样本来自分组人口预计将有较小的遗传变异,而不是整个人口。

突变,常见的假设是,突变发生的时间遵循一个毒药与恒定速率过程 的,因此,在任何分支长度 从树上,突变的数量分支有毒药分布的意思 ,独立于其他分支上的突变。前面提到的时间尺度,通常 ,在那里 发生突变的概率是每序列生成。DNA序列, 是序列长度(基地)乘以突变率/网站/生成和通常可以从其他来源。合并的时间以来与温和的大小大约是样本 代, 大约可以解释为累积(MRCA自从)突变速率(每序列突变的数量)。同时,由于人口规模 , 也可以理解为整个人口每一代的突变率。

因此,鉴于突变率 和树长度 的突变 在一个样本 个人从给定的人口是阿宝的毒药分布 (10] 注意,这并不取决于概率 ,但在 , , 。为什么 ?。取 ;然后 ,这是预期数量的累积突变 代前的一个序列。所以 是一个毒药分布参数的合理选择。但是如果我们模型数量 累积突变的序列 代之前,适度 ,我们应该使用

在聚结推理的关键是评估的postdata分布 ,这是比它更涉及predata分布,它在很大程度上依赖于突变分布数据。例如,如果突变发生在早期的家谱树,然后估计 将会更大。尽管假设突变最多只能发生一次在每个站点和变异是不可逆转的,突变在观测数据的总数是隔离站点的数量。但突变如何分发的谱系树的分支是未知的。这种分布的推理是至关重要的 ,这取决于所使用的数据信息和实际方法。这是我们的重点。表示由 观测数据,估计合并的时间 的样本是由postdata分布的均值 作为 可以视为一个贝叶斯推理过程,与predata postdata分布对应于先验和后验分布在贝叶斯框架。但与普通贝叶斯设置,这里的参数 随样本容量 建模,数据不能i.i.d.这个参数。这是原因的推断 不能任意精确,在某种意义上,不能任意postdata分布的方差小,随着样本量的增加没有绑定。一般postdata分布也不是在封闭的形式,由抽样评估方法。Tavare et al。10)派生postdata分布仅基于隔离网站在样本的数量。这个方法很方便使用,但是不使用DNA序列的结构信息。众所周知的方法在格里菲斯和Tavare [2),以后GT,是基于完整的数据信息由一组根树。这种方法的一个基本工具合并推理使用完整的数据信息,但计算复杂。

评估postdata联合分布,GT使用概率的递归公式,推导出在Ethier和格里菲思(11]。方法不容易完全理解和正确使用许多遗传学家。也这些概率计算限制性;的postdata分布 由马尔可夫链蒙特卡罗采样和计算是很复杂的。

在这里,我们研究一个相对简单的近似方法充分利用数据信息;在这种方法中,而不是计算树概率在GT,我们只是将post数据树概率的均匀 根树和使用模拟方法计算联合分布;因此,绕过复杂的评价树的概率,很容易理解和更简单的计算。

根树分析中扮演着重要的角色,这不是唯一确定的数据。相当于一个拔起树的数据,相当于一套拔起树。每个根树有一个0 - 1的矩阵表示一些计算方便,但并不是任何价值0 - 1矩阵对应于一个根树。在下面,我们提供更多细节关于他们和他们的关系。

根树。根树的一个分支系统,支行,等等。每个分行或支行的尖端代表一个已知的血统。观察到的样本被表示为点突变的分行,支行等等在指定位置。每个家族的观察到的多样性表现为叶子在每个分行或支行的,等等。

根树的表示是独一无二的相对位置的分支,支行,等等。根树有几个级别的随机性。如果我们只知道样本的大小 ,然后根树共有 叶子;除此之外,这棵树的形状,如何分割,如何分配树叶,有多少突变,突变的分布都是随机的。如果数据和突变的数量,然后树可以只需要几形状。不同于GT和其他相关文献,我们将观察到的血统频率(多样性)的叶子在相应的分支机构,支行等等的树。

不同于接合的树,有一个完整的时间排序的分割点树枝,扎根树只有局部时间序的分裂和变异。我们只知道分歧分支(es)发生前支行的,但是不知道分裂的排序不同的分支。我们知道突变(s)在树枝上发生过那些在其支行(es),但是不知道订单的同一分支,同样支行(es),或在不同的支行。对于一个给定的序列数据,它可能对应多个不同的树。观测数据的表1隔离网站的所有列,没有颠换。假设突变最多只能发生一次在每个站点和变异是不可逆的,在每个隔离的网站,一个且只有一个基本类型的突变;另一种类型是祖传的。如果我们知道突变状态在每个隔离的网站,据说突变状态标签,我们可以使用一个0 - 1的矩阵 表示观测数据, ,如果血统的基类型 在网站 突变体, 否则。这样的0 - 1矩阵表示数据方便分析。很容易看到,每一根树唯一确定一个0 - 1的矩阵 ,但一个任意的0 - 1矩阵价值可能不对应于一个根树。它必须满足一些条件对应一个根深蒂固的树。有丰富的方法和算法如何判断一个给定值0 - 1矩阵是一个有效的根树,如果是如何构建的树(例如,12- - - - - -16])。我们发现的方法出现在一系列的文章和规定为引理1 Gusfiled [16很容易使用。给出一个有效的价值0 - 1矩阵 (代表树)意味着它满足条件,可以唯一地对应,画一根树。这里,独特性意味着家谱关系,包括哪些血统在同一分行或支行,等等和哪些突变网站部分分行或支行等等,决心,但树的特殊形状,比如一些树枝放在左边或者右边,分支机构的角度,他们的长度,等等,都是无关紧要的。因此,根树之间有一个一对一的对应关系,一个有效的价值0 - 1矩阵。考虑到观测数据,网站的突变状态通常是未知的。数据与 隔离网站,有 不同方法标签突变雕像,但大多数标记矩阵不符合根树的表示;众所周知,只有 不同的树,因此 不同的标签(矩阵)对应数据,还有现有算法构造的树和它们相应的矩阵(例如,16,17])。然而,我们发现该方法在GT是方便。用这种方法,一分之一需要构造一根树的数据或其有效价值0 - 1矩阵。例如,从最开始共享突变标记,在每一列(网站)的数据,标签不太常见的基类型的突变(其他祖先)。很容易检查其有效性使用引理1上面提到的条件。构造的树对应于这个矩阵,将它转换成一个拔起树GT;即吸收无突变支行到分行(es),然后伸直分行,支行,等等。拔起树是唯一确定的 根树。

然后,在此基础上拔起树,一个人可以获得所有的其他的树木在格里菲斯和Tavare [18];,或者把树根附近点的每个顶点顶点伸出,然后安排分行,支行,等等成所需的形状;如果有不止一个两个相邻顶点之间的变异,把树根两个相邻的点突变,或者等所有成对的突变,和形状上面的树。这样我们得到的所有根树拔起树。事实上,给出任何根树,所有其他的 根树上面可以以同样的方式构造的,不用拔起树。一旦扎根树构造,相应的矩阵表示。

3所示。MRCA估计的渐近行为

参数推理与独立同分布数据和样本大小 ,众所周知,估计是渐近一致和渐近正态的 。但MRCA推理,数据 不是独立同分布,现有研究表明,分布的估计MRCA吗 不能集中注意力,即使 。估算的突变率与相同的数据,估计是发现一致性和渐近正态的 (1]。这激励我们调查的渐近性态 作为常用的点估计量的联合。我们想知道这是否估计也有类似的突变速率估计量的渐近行为。我们发现这样的估计几乎肯定是不一致的。描述结果,我们认为数据集在三个不同的常用形式。第一种形式的数据,我们考虑的是合并树如图1。这种类型的数据通常是不实际的,至于最真实的数据我们没有信息构建这样的树。但作为起点会为我们提供一些指导的结果。有 树中的节点(分割点)编号2 他们的时间顺序。召回的定义 th合并的时间 。之间的 th和 节点有确切 段,表示它们 从左到右,每个人都有长度 。假设突变的数量 上段 是已知的。让 , 相对应的突变分布 。这种类型的数据是完全为代表 。当我们没有 , 不是唯一确定。但考虑到每一根树 , 和位置信息的突变,突变向量 可以由一个随机的方式(详细节4)对应 。表示 是我们之前在根树 的,没有额外的知识我们平等对待每个根树可能从观测数据。在这里我们 从概率的有不同的意义 的年代,GT(后者不总结,但获得的概率拔起树从观察到的数据)。我们有(附录)

常用数据表的形式1,相当于 根树;在这里

最后一个方法是估计 只有通过突变的数量 在根的树木,不使用信息。

我们有以下结果(证明在附录)。

命题1。(我)一个 因此,上述估计会发散几乎可以肯定的是,如果 看作是随机的。
(2) 将收敛与否取决于上述系列。
(3) 和上面的估计量的渐近行为取决于上述系列。

备注2。上面的结果告诉我们 不能由一个渐近的确定性特征量,即使对于大数据的大小。估计量是由突变的数量在最初几个合并的时代。因此,唯一可行的方法来推断出合并的时间是通过数值方法,随着postdata联合分布甚至没有关闭形式渐近。相比之下,predata的意思 收敛,但不准确的估计量的聚结时间下的人口研究。

4所示。该方法

方法是构造变异向量 和计算概率的数据直接从家谱树扎根 的,假设有 序列中的隔离网站数据,这正是历史上发生突变的总数 个体取样,然后 不同的树 兼容的数据。的每一根树是一个固定的家谱结构,以多样性为叶子,但是树片段中突变的数量是随机的,突变的总数 。树的结构由分行,支行在每个分支,sub-subbranches、等等,和树叶。这些都是固定的特征根树。考虑到数据,根树的显示 突变分布沿血统,但是没有时间尺度的树,(5)不能用于计算变异概率。每一根树告诉我们的部分点突变。例如,在根树,我们知道突变在网站4,6,14发生在分裂之前的血统 , , , ,从而发生突变在网站1日之前5和10。但是我们不知道哪一个4,6,14日发生。我们知道突变1发生前10,但是我们不知道1和5的顺序,等等。如果我们有完整的数据 对应所有的树木, 的年代,我们可以计算 在命题1(二)。但是 s不是直接可用;然而, 可以很容易地模拟的指数分布之前,每一根树 有一个初始突变分布在其分支段。表示由 th段(顺序是任意的,例如,我们可以从上层标签他们到低,从左到右位置),并让 是突变的数量(很多都是零;我们可以专注于非零变异的片段)。表示 。鉴于 , 可以从采样 后者(详细)。让 期望对 。上述激励我们估计 通过

上述预期是不容易直接计算,因为我们不知道的联合分布 。相反,我们使用仿真方法。为此,我们样品 独立并生成 (见下文)对应 为每一个 然后近似 作为

现在我们考虑生成 。后 分配的分支 ,我们只需要考虑每一部分 非零号 的突变。每一个的长度 的总和吗 博览会和符号。为了简单起见,假设 ;然后考虑到 突变 和使用(5),很容易发现突变的数量 在每一个 时间间隔 , 遵循了多项分布 在哪里 。在所有的非零 在相应的时间间隔,分配 求和是为所有在哪里 的下降

具体来说,仿真方法如下。为 ,请执行以下操作步骤。(我)样本 从联合分布(1);也就是说, 是独立的, 。或者说,样本 并设置 (2)对于每一个固定 ,分配 合并的事件 根据每一根树序列 。有关详细信息,请参阅下面的说明。(3)分配 突变在相应的部分根据(13)。然后得到 的年代(14)。

在所有的 迭代,计算(12),直到收敛,可以评估相对误差,例如。

说明:分配 合并的事件 基于根树的序列。我们使用落后的方法;也就是说,第一次分配 ,然后 , ,最后 。例如,考虑的树。有 序列,与频率 的血统 。注意序列(叶)在每个血统(分支)只有聚结在每个分支(如果有多个分支的叶子),和分支与单个叶只在MRCA聚结 。我们首先决定 去哪个部门或对单一的分支。因为它是最新的合并的时候,它只能去一对叶子在一些分支与多个叶子。因为分支 只有1片;它是排除在这一步。剩下的树枝 都有多个叶子共有54。我们分配 其中的一个分支与重量成正比的叶子,也就是说,与概率 。假设 被分配给 ;我们需要决定去哪个支行。我们有三个候选人支行 与数量的叶子 。我们随机分配 他们与权重 。假设它是分配给分公司 ;然后 将一对在这个分支,这是无关紧要的。但两人将被视为一个叶分配 ”年代,所以这一步后,我们重新分配的叶子 为18。

现在我们分配 。上面的程序是一样的;唯一的区别是现在 有18个叶子。候选人仍分支 与重量 。假设 也分配给 的分支 ;然后 现在有17个叶子。

现在我们分配 给候选人 与重量 。假设 被分配给 ;我们需要决定哪些4支店去。 只有1叶和排除在外。所以我们分配支行 与重量 。假设它进入 ,因为它只有一条叶子, 这个任务后合并为一片树叶。

继续这样,直到 分配。那么这个根树的分支长度的 分配给他们。在这一步中,每段的长度 总结了一些吗 以来的年代。 被从每个 ,我们可以分配的每一个 的年代(13),然后得到 通过接下来的公式。然后计算(12)。

假设人口规模 是常数可以放松一样在GS和Tavare et al。10]。

附录

证明的命题

(我)召回的 是独立的, 是独立的, 是独立的 , , , ,所以 。观察 右边上面的密度 分布,一个正常化常数。它的意思是 和方差 。自 ,我们有 在哪里 , 的是我。维与 , , , 。请注意, 收敛序列有 还请注意, ,也就是说, 为每个固定权重序列 。自 固定 通过i.i.d.随机变量加权和的结果(见[19),审查的结果),一个必要条件 收敛(主导者——) 对所有固定 ,或者没有任期将占主导地位 。因为这个条件不满意,是占主导地位,我们的前几项 部分的证明(2)和(3)是类似的(我)和省略的部分。

利益冲突

作者宣称没有利益冲突有关的出版。