第六章 溯祖

模拟 DNA 序列样本

为了准确推断影响分子变异的进化力量，我们必须能够准确模拟 DNA 序列。理想情况下，我们希望能够在各种模型、不同的人口历史、选择形式和强度以及任何其他感兴趣的参数下模拟种群。此类模拟的结果应该是 DNA 序列样本，相当于我们从自然种群中收集的样本。这个过程重复了数千次，然后可以用来找到最适合我们观察到的数据的模型，或者作为零分布，我们可以用它来测试特定的假设。

我们如何模拟这样的种群？最直接的方法是，从大量二倍体个体（相当于我们想要模拟的种群规模）开始，这些个体根据我们选择的种群模型进化。我们将对适当长度的序列（可能等于我们希望与模拟进行比较的 DNA 序列的长度）应用突变和重组，同时进行交配和任何可能的自然选择。然后，该系统必须运行多代，至少直到种群达到平衡，最后需要从种群中抽取少量染色体以创建一个模拟数据集。为了生成数千个独立模拟的样本，这个过程必须运行数千次。这个过程通常被称为正向模拟，因为它从最初相同的种群开始向前进行。

从上面的描述中可以明显看出，这种方法非常低效：我们只想知道样本的属性，但我们正在跟踪整个种群。为了解决这种低效率问题，几位数学遗传学家独立提出了后来被称为合并过程的方法（Kingman 1982a、1982b、1982c；Hudson 1983a；Tajima 1983）。合并过程的核心是，我们可以根据生成谱系间关系集（系谱）的简单规则来生成 DNA 序列样本。这些规则使我们能够生成具有几乎任何人口历史的样本，而无需模拟整个种群。由于该过程从现在开始并向后运行，因此被称为向后模拟。如今，有许多软件包可用于执行合并模拟（表 6.1）。

由于合并过程仅考虑样本而不是整个种群，因此它可以成为生成大量模拟数据集的一种非常有效的方法，但它也存在自身的局限性，无论是生物学还是计算上的。从生物学上讲，合并是模型的模型：也就是说，它是 Wright-Fisher 或 Moran 种群（尽管非常好）的近似值，而 Wright-Fisher 或 Moran 种群本身是自然种群的近似值。合并的一个关键假设是样本大小 n 比有效种群大小 Ne 小得多。违反这一假设可能会导致对种群过程的误导性推断（例如，Wakeley 和 Takahashi，2003)。虽然可以修改标准合并器以更精确地匹配 Wright-Fisher 模型并放宽 n << Ne 假设，但这种修改会导致算法速度降低 (Fu 2006)。从生物学角度来看，更重要的是，合并器仅从没有自然选择的种群中生成少量随机染色体样本。它不是对整个种群的模拟，因此无法提供除样本属性之外的任何结果。尽管模拟合并器的方法 w尽管选择方法一直在改进（例如，Spencer 和 Coop 2004；Teshima 和 Innan 2009；Ewing 和 Hermisson 2010；Kern 和 Schrider 2016），但可以建模的选择形式仍然有限。

从计算角度来看，当计算机处理器速度和内存资源有限时，合并器提供的优势更为重要，而这已不再是一个问题。技术发展使得合并器的效率变得不那么重要，除非在使用贝叶斯采样方法的应用中（参见第 9 章）。然而，最重要的是，当模拟基因组的非常大的区域或具有非常高重组率的区域时，合并器会变得非常低效，以至于无法使用它来建模当前正在收集的数据集。解决这个问题的一个办法是发明近似合并方法，该方法基于在沿着 DNA 片段移动时生成相关谱系，而不是为整个区域生成一个大谱系（Wiuf 和 Hein 1999）；因此，它们是横向模拟。这些方法更正式地称为顺序马尔可夫合并 (SMC) 模型（McVean 和 Cardin 2005；Marjoram 和 Wall 2006），并且有多个程序可以执行它们（表 6.1）。SMC 模型已经变得非常有用，但它们与自然种群的距离甚至更大：它们是模型的模型的模型！

在未来，正向模拟可能会成为此类方法中最广泛使用的方法。越来越多的快速灵活的正向模拟器可以在许多不同的自然选择形式和许多不同的人口历史下对非常大的基因组区域进行建模（表 6.1）。计算技术的进步使这些方法更加可行，而 DNA 测序技术的进步可能使它们成为必要。

虽然到目前为止，我只介绍了合并算法作为生成样本的算法，但这种方法有许多重要用途（请参阅 Hein、Schierup 和 Wiuf 2005 和 Wakeley 2009 以获得良好的概述）。合并算法为我们提供了一种概率建模谱系的方法，因此为分子群体遗传学中的各种结果提供了预期。这些预期中的许多都用于本书的后面章节。合并算法还为我们提供了一个从谱系角度思考的框架，因此可以让我们深入了解我们采样的 DNA 序列之间的潜在关系。这个框架对于直观地了解种群遗传过程非常有帮助，因此在本章的其余部分，我将讨论中性联合谱系的基础知识，并解释联合谱系的一些最重要的方面。在后面的章节中，我们将研究自然选择和非平衡人口历史对某个基因座谱系的影响。

模拟溯祖

溯祖谱系

来自一个种群的任何同源染色体样本都以某种方式相关：也就是说，它们在最近的过去都有共同的祖先。一些序列可能比其他序列更紧密相关，但随着时间的推移，它们最终都拥有共同的祖先。如果我们将上述向前时间模拟想象为 Wright-Fisher 种群模型，那么在最近一代中采样的任何染色体集合都将在某个基因座处具有一组特定的关系，我们称之为谱系（图 6.1；有时这些被称为基因谱系，即使该基因座不是基因）。 回溯到过去，两个在特定世代中拥有共同祖先的序列被称为融合成一个谱系； 向前看，这个融合事件代表一个染色体产生两个子染色体。 目前，我们只考虑单个种群中没有重组和选择的单倍体个体，Hudson（1990）称之为“典型的普通单倍体物种”。个体中的 N 个二倍体可以通过随机连接 2N 对单倍体形成。最终，所有谱系合并为一个序列 — 即我们样本最近的共同祖先的染色体。

合并是一种生成与我们的样本染色体相关的谱系的方法。通过反复运行向前时间的 Wright-Fisher 模拟形成的谱系可以用一组非常简单的规则来近似，这些规则将概率放在合并事件之间的时间长度上，然后可以将其转换为一组由谱系表示的关系。我们可以在生成的谱系树上随机放置突变，并在底部读取相当于 DNA 序列的内容。这两个步骤（生成树、添加突变）实际上就是为我们提供一组近似于从 Wright-Fisher 群体中采样的序列所需的全部步骤。接下来我将更详细地描述这些步骤。生成合并谱系的一种直接方法基于一个简单的近似值。如果我们用连续近似值替换 Wright-Fisher 模型中使用的离散代，那么距离下一个合并事件的时间就是从指数分布中抽取的随机变量。该指数分布仅由可能合并的剩余谱系数量参数化，这意味着当谱系数量较多时（例如，朝向树的尖端），合并事件之间的时间较短，而当仅剩下最后两个谱系时，合并事件之间的时间最长。使用指数分布，我们现在可以写下五个简单的步骤来为样本生成谱系

1. 从 i = n 条染色体开始。

2. 从参数为 x = i(i − 1)/2 的指数分布中选择下一次合并的时间。

3. 随机选择两条染色体

4. 合并所选的两个谱系

并设置 i → i − 1。

5. 如果 i > 1，则转到步骤 2；如果不是，则停止。

例如，图 6.2 显示了按照上述步骤生成的一个可能的谱系，其中 n = 5。从 i = n = 5 开始，该过程选择两个谱系进行合并，此时 i = 4，发生另一个合并事件，i = 3，依此类推。每次运行此过程时，都会创建一个略有不同的谱系，如图 6.3 所示。要生成具有非平衡种群历史的树，我们需要调整合并事件之间的时间，但这些转换相对简单（参见第 9 章），并且生成的树可以以相同的方式解释。这棵树中的时间以 N 代为单位（单倍体模型）、2N 代为单位（二倍体模型），有时甚至以 4N 代为单位（为方便起见）。

图 6.3 大小为 n = 5 的合并谱系的高度和拓扑变化。每棵树都是从完全相同的种群中独立生成的谱系。模拟是在 discoal 中进行的（Kern 和 Schrider 2016）。

图 6.2 样本大小为 n = 5 的合并谱系示例。合并事件之间的时间 Ti 相当于谱系中有 i 个谱系的时间量。

合并谱系的构建完全独立于突变的存在与否。无论是否存在区分它们的突变，染色体之间都存在一组关系。这意味着我们可以在第二步将突变应用于我们的谱系，而这不依赖于我们的谱系构建方式的细节，尽管突变的数量和位置可能取决于树。这种独立性还使我们能够在任何突变模型（无限位点、无限等位基因等）下使用任何类型的突变（SNP、indel、转座因子等）。

有两种方法可用于将突变置于谱系中。第一种方法称为固定 S 法 fixed-S method（Hudson 1993），通过在树上放置预定数量S的突变来进行。此方法可用于直接比较具有一定数量分离位点的样本和具有相同数量变体的模拟谱系。当研究除多态性水平以外的某些变异方面（例如等位基因频谱）时，这些类型的比较最为常见。由于分离位点的数量是预先指定的，我们只需将 S 突变随机放置在树的分支上，概率与分支长度成比例。

当最初未指定 S 时，使用第二种方法：在这种情况下，我们需要使用种群突变参数θ（二倍体中常染色体基因座 =4Nμ）生成要放置在树上的突变数。使用 t 表示谱系中任何特定分支上的时间，每个分支的突变数呈泊松分布，平均值为 t * θ/2。在这里，我们将突变的生成与它们在树上的位置结合起来：使用这种方法，较长的分支将再次具有更多的突变，因为更多的突变将在较长的时间段内出现。为了计算简单，我们还可以对树中的所有时间 Ttotal 求和，首先从泊松分布中绘制突变ution，其平均值为 Ttotal* θ/2 ，然后像在固定 S 方法中所做的那样，将它们全部随机“扔”到树上。

图 6.4 显示了具有四条染色体和六个突变的联合谱系。如前所述，这些突变可以代表任何突变模型下的任何事件；然而，为了简化讨论，让我们将它们视为无限位点模型下产生的单核苷酸多态性。为了将这些突变变成类似于 DNA 序列的东西，我们从四条染色体中每条染色体的长度 L = 6 的序列开始，最初将每个位置设置为 0（祖先状态）。然后，我们将每个突变沿树向下传播，这样，如果所有后代谱系继承了等位基因的突变版本，它们都会收到 1。图 6.4 中每个样本下方的字符串显示了结果序列，位置对应于树中从左到右、从上到下编号的突变的字符串。请记住，我们仍在处理非重组位点，请注意，比对中的突变顺序没有任何暗示，尽管在本章后面我们将考虑重组的影响，此时顺序很重要。可以看出，发生在谱系外部分支上的突变会导致单例突变，而靠近根部的突变通常会在样本中更多的个体中发生。然后可以像任何 DNA 或蛋白质序列一样处理此过程产生的“序列”，每个 0/1 位点代表祖先/衍生的等位基因。

图 6.4 将突变置于合并谱系中。六个突变被随机放置在左侧的合并谱系中。从由序列 [0, 0, 0, 0, 0, 0] 表示的祖先开始，衍生的等位基因由从 0 → 1 的变化表示。底部的序列代表结果单倍型，也显示在右侧的比对中。比对中突变的顺序没有意义。

谱系和系统发育

虽然合并谱系在许多方面类似于典型的系统发育树，但这两种类型的树之间存在非常重要的区别。从一个种群中采集的样本彼此之间通过一些未知的谱系相关。事实上，描述样本间关系的树通常是未知的，几乎不可知：即使没有突变区分我们的样本，仍然有一个谱系树。我们绝不会试图从等位基因的身份推断出树，就像在系统发育学中所做的那样，或者在为物种内的线粒体 DNA 或 Y 染色体单倍型构建“基因树”时所做的那样。此外，正如下面将更详细讨论的那样，在重组的情况下，基因座中的每个核苷酸位置实际上可能存在不同的谱系，甚至排除了推断单个树的可能性。相反，上面提出的合并算法会生成这些染色体之间的一组可能的关系，这些关系可能与我们的样本完全匹配，也可能不完全匹配。合并的明确目的之一是生成许多不同的谱系，有效地将树视为必须平均的干扰参数。染色体间的关系集本身就是数据集中变异的原因，这一观点突出了合并过程中以及任何分子数据集中的两个随机性来源：由关系集（即谱系）采样引入的方差和由具有任何给定关系集的染色体组采样突变引入的方差。回顾我们用于进行合并模拟的两个步骤（生成树、添加突变），这些变异来源显而易见，每个步骤都具有重要的随机性元素。通常，这两个随机性来源分别称为进化方差和抽样方差。使用显式假设检验时，了解我们必须考虑哪些变异来源以及哪些可以忽略非常重要（参见第 7 章和第 8 章）。

理解合并

重要的谱系数量

到目前为止，讨论还没有触及合并的数学基础，也没有触及合并谱系的不同方面究竟如何影响分子多样性的测量。虽然对合并的深入理论理解超出了本书的范围，但下面我将讨论四种数学谱系描述，这些描述至少应该能阐明合并和种群遗传过程之间的联系。

合并事件的概率 在任何给定的一代中，显然要求每条染色体在前一代中都有一个祖先（图 6.1）。如果没有自然选择，当前代的每个染色体都可以随机选择一个祖先；或者，我们可以说，上一代的每个染色体都有相同的概率成为当前代中给定染色体的父代。在 Wright-Fisher 群体中，大小恒定为 2N，这意味着有 2N 个可能的祖先可供选择，每个祖先被选中的概率为 1/2N。回想一下，当两个个体选择同一个祖先时，就会发生合并事件，对于任何特定的一对染色体，其概率为 1/2N：也就是说，以一对染色体中的一个选择祖先序列为条件，第二个染色体有 1/2N 的概率选择相同的序列。如果我们有一个大小为 n 的样本，那么就有 n(n − 1)/2 对可能的序列可以选择上一代中的同一个祖先。一般而言，在过去的任何一代中，对于 i 个谱系，有 i(i − 1)/2 个可能的配对可以以 1/2N 的概率选择相同的祖先。因此，在一代中发生任何合并事件的概率（即从 i 到 i − 1 个谱系的概率）是这两个项的乘积：

此概率假设在任何一代中最多有两个染色体可以选择相同的祖先，这一假设与 n << N 的要求密切相关。我们将在下一节中看到单个合并事件的概率如何导致进一步的洞察，但现在请注意，公式 6.1 意味着在给定的一代中，样本谱系的数量越多，发生合并事件的可能性就越大。

合并事件之间的时间 i 个谱系之间发生合并事件的概率可以直接告诉我们 i 个谱系的预期持续时间。任何特定代中，在剩余 i 个谱系的情况下，均不会发生合并事件，概率为 1 – P(i → i − 1)。如果没有发生合并事件，则发生合并事件的概率为 P(i → i − 1)，前一代中不发生合并事件的概率为 1 − P(i → i − 1)，依此类推。因此，合并发生的预期平均等待时间（即从 i 到 i − 1 个谱系所需的时间）——此处称为 Ti（图 6.2）——由指数分布近似，其平均值是公式 6.1 的倒数：正如谱系数量越多，发生合并事件的概率越高一样，公式 6.2 表明合并事件之间的时间越短。这意味着，例如，从 10 个谱系到 9 个谱系的时间将比从 5 个谱系到 4 个谱系的时间要小得多（例如，图 6.3）。平均而言，合并事件之间的时间在树的尖端最短，而在等待最后两个谱系合并为整个样本的最近共同祖先 (MRCA) 时最长。

合并谱系的高度我们可能对样本中每个谱系合并为单个谱系（MRCA）之前的预期时间感兴趣。因为我们知道从

i → i − 1、i − 1 → i − 2、i − 2 → i − 3、…、2 → 1 的预期等待时间为 Ti

通过将这些时间相加来确定平均树高：

对于两个染色体的样本，公式 6.3 表示树高（到最近共同祖先的时间）平均为 2N 代。随着样本量变大，树的高度接近 4N。

树的总长度树的总长度是树中所有分支的总和，因此该值将决定在样本中观察到多少突变。我们可以得出预期的总树长 Ttotal

注意公式 6.2 给出了预期存在

i 个谱系的时间量，因此在此期间所有谱系之间的总分支长度为 i * E(Ti

)（见图 6.2）。将存在

i 个谱系、i – 1 个谱系、i – 2 个谱系等的总时间量相加，可得出：

对于 n = 2，预期的总树长为 4N（是 TMRCA 值的两倍，因为只有两个分支），并且随着样本量变大，树会变大。但请注意，样本中每增加一条新染色体，树的长度就会越来越短：n = 5 时，预期长度为 4N * 2.08，n = 10 时为 4N * 2.83，n = 20 时仅为 4N * 3.55。树的长度（和高度）大部分仅由两条染色体贡献，而增加更多个体可能只会增加更多在树的尖端发生的融合。

融合和多态性测量

我们现在可以将上述谱系测量与从染色体样本中收集的实际数量联系起来。我们还可以非正式地证明，用于测量样本多样性的几个统计数据（第 3 章）实际上是 Wright-Fisher 平衡种群下种群突变参数 q 的估计量。回想一下，我们早期r 将统计量 p 定义为任意两个序列之间核苷酸差异的平均数。从公式 6.3 可以看出，任意两个随机采样序列的 MRCA 平均时间为 2N 代。假设突变以每代 m 的概率发生，并且种群以有效大小 (Ne) 来衡量，则预计两个序列之间会有 2 * 2Ne

m 个差异。这

表明在我们的中性平衡假设下，p 有望成为 q 的估计量。

虽然统计量 p 基于序列之间成对差异的平均数量，但我们根据分离位点的数量 S 定义统计量 qW

，因此：

其中 a 等于：

公式 6.6 中的术语应该很熟悉 - 它用于定义合并树的总长度 Ttotal

（公式 6.4）。很明显，收集更多序列将使我们能够找到更多分离位点，但 qW

tic 对我们每个额外样本的惩罚越来越少。如上所示，每个额外的序列也会越来越少地增加总长度合并树，这当然会导致添加的新多态性更少。我们可以通过重新排列公式 6.5 来明确地看到这种模式，得到 S = qW

a。再次假设 (1) 突变以每代 m 的概率发生，并且在整个树中均匀随机地发生，以及 (2) 无限位点模型，则样本中分离位点的预期数量为：

( ) *4 1 4 *

ES N i N a a

使用公式 6.4 和 6.6 中给出的关系。这有助于我们看到在我们的中性平衡假设下，预计是 q 的估计量。合并谱系的结构为我们提供了一种图形表示等位基因频率和等位基因频谱的方法。我们之前定义了 i 条染色体中存在衍生等位基因、n-i 条染色体中存在祖先等位基因的分离位点数，即 Si。从树的结构中我们可以看出，对于 i 个序列上存在的任何特定衍生等位基因，突变必定发生在具有恰好 i 个后代的分支上（图 6.4）。例如，在样本中出现三次的等位基因必定出现在具有恰好三个后代的分支上。

这里我们假设一个无限位点模型，因为在有限位点模型下，同一位点可能在后代谱系中恢复到祖先状态，或者在树的完全不同的分支上突变为相同状态。此类事件将使单个突变与其频率之间的对应关系复杂化。

单例突变（仅存在于单个染色体上的突变）发生在直接导致当今样本的分支上。这些有时被称为外部分支，与所有导致两个或更多后代的分支相对，这些分支被称为内部分支（Fu 和 Li 1993a）。并非所有等位基因频率都可能出现在每个家谱中，因为每棵树中不一定都有 2、3、... n - 1 个后代的内部分支（总是有 n 个分支有一个后代）。家谱的确切拓扑结构将决定哪些等位基因频率是可能的，从而决定可能的等位基因频谱。但是，对所有家谱进行平均，我们可以通过以下公式预测任何频率（即等位基因频谱）的预期分离位点数：其中 i 从 1 变为 n - 1（Fu 1995；Griffiths 和 Tavaré 1998）。这意味着

有 q 个单重态位点、q/2 个双重态位点、q/3 个三重态位点等等，并且

单重态位点将是最常见的位点类型。公式 6.8 还说明了为什么统计量 qe

；第 3 章）是 Wright-Fisher 平衡种群中 q 的估计量，并且一般来说，iSi

是

1 和 n – 1 之间所有 i 的 q 估计量。作为 q 的估计量，对于任何样本，我们可以将每个频率的中性等位基因的预期数量写为：

其中 S 表示样本中分离位点的总数，a

定义如公式 6.6 所示。图 6.5A 显示了 S = 28 和 n = 11 的基因座上预期的中性等位基因频谱的示例。

以上所有讨论都涉及我们所谓的展开频谱，我们能够为其分配祖先和派生等位基因（第 3 章）。对于折叠等位基因频谱，我们不知道 i 条染色体上存在的等位基因是由于在具有 i 个后代的分支上发生的新突变还是在具有 n - i 个后代的分支上发生突变后残留的祖先等位基因。无论如何，我们

也可以写下频率为 1

到 n/2 的分离位点的预期数量（对于奇数大小的样本向下舍入）：

我们可以再次用 qW

代替 q 的估计量，得到：

图 6.5B 显示了 S = 28 和 n = 11 的基因座的预期折叠等位基因频率谱，与图 6.5A 相同，但具有h 没有分配衍生和祖先等位基因。

确定偏差对等位基因频谱的影响

等位基因频谱的形状是分子群体遗传学中的一个重要指标。在没有选择影响的情况下，了解其在平衡群体中的预期形状意味着我们能够检测到

合并的

1 2 3 4 5 6 7 8 9 10

具有衍生等位基因的个体数 (i)

具有次要等位基因的个体数 (i)

图 6.5 根据样本中的频率预期分离位点数。给定一个 S = 28 和 n = 11 的基因座，频率 i [=E(Si

)] 处的多态性预期数量显示为 (A) 展开谱和 (B) 折叠谱。展开谱显示衍生等位基因的预期频率，而折叠谱显示次要等位基因的预期频率。与此预期的偏差，无论是由于非平衡人口历史（参见第 9 章）还是连锁选择（参见第 8 章）。然而，如第 2 章所述，如果我们对先前从较小样本中确定的大样本中的一组变异进行基因分型，我们会在等位基因频率谱中引入偏差。初始发现样本规模小的一个重要结果是，我们更有可能发现中频变异，而不太可能发现任何特定的低频变异（即次要等位基因频率较低的变异）。这种趋势意味着从预先确定的多态性推断出的等位基因频率谱与理论预期相比可能存在很大偏差（图 6.6）。正如关于预先确定的多态性和连锁不平衡（第 4 章）和 FST（第 5 章）的测量所讨论的那样，各个分离位点的样本频率没有任何不准确之处。这些估计值预计是正确的频率，并且可能

图 6.6 确定偏差对等位基因频谱的影响。给定一个大小为 k 的发现样本，显示了第二个独立样本大小为 n = 41 的预期折叠等位基因频率谱。无偏样本代表没有先前发现过程时的预期频谱。

（基于 Marth 等人，2004 年。）

1 3 5 7 9 1113

具有

次要等位基因（i）的个体数

频率 i 处的 SNP 的预期比例

甚至在发现样本和第二个（通常更大）样本之间相同（如果它们来自同一人群）。但是，如果没有预先确定，基因分型样本中每个频率的多态性预期比例将与预期值不匹配。事实上，即使发现样本大于基因分型样本，也会产生轻微的偏差，主要是因为第二个样本独有的低频等位基因将不会出现。如果确定方案已知，我们可以对其进行校正以恢复底层（无偏）等位基因频谱（例如，Nielsen 和 Signorovitch 2003；Nielsen、Hubisz 和 Clark 2004；Clark 等人 2005）。但是，没有完美的校正，在解释这些更新的光谱时必须谨慎（Clark 等人 2005）。

扩展合并

具有重组的合并

上面的讨论集中在足够大（或突变率足够低）的基因座的属性上，以满足无限位点模型的假设，但足够小以至于不会发生基因座内重组。这些假设对于大多数数据集来说显然是不现实的，所以我们现在必须考虑重组对合并谱系的影响。如前所述，序列中的每个核苷酸都有自己的谱系，描述样本中染色体之间的关系。如果一个核苷酸和它紧邻的核苷酸之间没有重组，那么它们都具有相同的谱系（图 6.7A）。如果有自由重组——如不同染色体上的两个核苷酸之间——那么这两个位点将有完全独立的谱系。更有趣的情况发生在位点部分链接时：当存在有限的重组时，合并重组重组位置 1 位置 240 位置 241 位置 327 位置 328 祖先重组图图 6.7 重组对合并谱系的影响。（A）假设基因座上的谱系。从位置 1 到 240 没有重组事件，因此只有一个谱系。位置 240 和 241 之间发生重组事件，从而产生一棵新树。同样，位置 241 和 327 之间没有重组。位置 327 和 328 之间发生第二次重组事件，并创建一个向右移动的新树。（B）祖先重组图总结了三个独特的谱系s 表示在 (A) 中的位置。相对于位置 1 的树，重组引入的新合并事件用虚线表示。它们之间的重组事件数。就合并谱系而言，两个位点之间的每个重组事件都会导致树中的重新排列，从而导致合并事件的位置在重组事件之前发生变化（图 6.7A）。请注意，这两棵树仍然保留了上述单个树的所有属性 - 重组不会改变树的平均拓扑结构（Kaplan 和 Hudson 1985）。由单个重组事件产生的两个相邻谱系将不再相同，但将高度相似；许多属性，例如它们的高度和总长度，将高度相关（Griffiths 1981；Hudson 1983b；Kaplan 和 Hudson 1985）。这意味着许多变异度量（例如 p 和 qW）也将在关联的基因座之间相关联。越来越多的重组事件将导致谱系发生越来越多的变化，直到最终两个被许多重组事件分隔开的谱系看起来完全不同，实际上完全独立。

考虑具有有限重组事件数 R 的单个基因座，我们现在可以将其想象为一组 R + 1 个非重组片段，每个片段都有自己的谱系（如图 6.7A 所示）。我们将每个片段的单个谱系称为该基因座的边缘谱系。给定所有边缘谱系，我们可以构建一个祖先重组图 (ARG)，它表示所有单独非重组片段的组合历史（图 6.7B）。合并事件在 ARG 中再次表示为分叉的谱系，而重组事件导致谱系结合在一起（向前看）。尽管这种对 ARG 的描述使它们似乎仅限于一组边缘树的简单图形表示，但这些图在群体遗传学中有许多重要应用（例如，Song 和 Hein 2005）。对 ARG 的讨论还应该清楚地表明，对于具有重组的常染色体基因座，没有单一的树来描述样本之间的关系。相反，有一系列树，每个片段在其历史上没有重组事件。这意味着不可能推断出重组区域的单一、分叉的基因树，尽管使用 Y 染色体或 mtDNA 基因座仍然可以进行这种推断。正如上一节中介绍的联合谱系与核苷酸变异摘要之间的关系一样，考虑重组对谱系的影响可以更深入地理解总结样本中重组事件历史的统计数据，以及与核苷酸多态性的多个测量值相关的方差。在第 4 章中，我们考虑了配子连锁不平衡的两个测量值，D′ 和 r2。回想一下，“完整” LD 意味着 D′ = 1，“完美” LD 意味着 r2 = 1，联合谱系让我们直观地以图形方式了解这两个值的含义。因为两个多态性只要它们之间没有（可检测的）重组就可以有 D′ = 1，所以在联合上下文中，这个值只是意味着两个分离位点共享相同的谱系。突变不一定发生在谱系的同一分支上，因此样本中可能存在三种单倍型。相反，只有当两个多态性完全相关时，r2 = 1，因此样本中只有两种单倍型。这意味着位点必须具有相同的拓扑结构，并且突变必须发生在同一分支上或根节点的任一侧。我们还在第 4 章中看到了群体重组参数 r = 4Ne 的几个估计量。正如群体突变参数 q 与样本中的多态性位点数量成正比一样，r 与样本中的重组事件数量 R 成正比。如第 4 章所述，估计样本中重组事件的真实数量可能很困难，因为大多数此类事件是无法检测到的。然而，

合并谱系的结构仍然允许我们看出 r 和 R 之间的关系。使用 c 表示每代每个核苷酸的重组率，我们可以再次想象重组事件在空间和时间上均匀随机地发生——也就是说，重组可以在基因座的任何地方以相等的概率发生，我们可以像处理突变一样将重组事件“抛”到我们的合并谱系上（尽管重组改变了谱系的基本形状）。给定公式 6.4 中树的总长度的期望，因此，我们样本中重组事件的预期数量为（Hudson 和 Kaplan 1985）：

（）*4 1 4 *

ER c N i Nc a a

这与公式 6.7 中给出的分离位点数量的预期完全类似，其中 q 替换为 r。但是，公式 6.7 和 6.12 中给出的关系有三个重要区别。首先，正如多次提到的，大多数重组事件是不可检测的，因此使用可检测事件的数量（例如）不会得到 r 的良好估计量。因为大多数突变事件在无限位点模型下都是可检测到的——并且可以针对有限位点场景进行调整（第 3 章）——S 是 q 的良好估计量，至少在 Wright-Fisher 模型的均衡假设下是如此。其次，与单独模拟突变的情况不同，我们现在必须明确选择序列中发生单个重组事件和突变的位置。重组事件的顺序是生成某个基因位点边缘谱系的关键。第三，我们将重组事件放到哪个确切的树上可能并不明确，因为将重组事件放到树上必然会改变拓扑结构。考虑这种情况的一种方法是考虑研究区域左侧第一个核苷酸处的固定谱系。然后，我们可以想象样本中的每个重组事件都被放到这个固定的树上，新的谱系会在这个位点的右侧，也就是我们感兴趣的区域内生成。这种方法不是标准后向合并算法模拟重组事件的方式（参见 Hudson 1983b、1990），但它可能有助于更直观地理解该过程。

最后，回想一下，随着重组率的上升，种群突变参数估计量的方差（例如 p 和 qW）会下降（公式 3.4 和 3.7）。合并算法提供了这种关联的谱系视角。在没有重组的基因座内，有一个单一的谱系描述染色体之间的关系。这个谱系的高度和总长度本身是从方差较大的分布中提取的，而我们只有一棵树这一事实意味着实际上从这些分布中只提取了一个。在发生重组的位点中，更高的重组率意味着有越来越多的边缘谱系，而这些个体谱系彼此之间的独立性越来越强。实际上，这意味着我们从间接决定多态性水平（树高和树长）的分布中抽取了大量的独立样本；因此方差大大降低。在大量位点之间发生自由重组的极端情况下，进化方差基本上没有影响，因为我们对大量独立树进行了平均。唯一的方差是通过在树上抽样突变而引入的，因此对于 qW，它相当于我们选择突变数量的泊松过程的方差（对于泊松分布，平均值等于方差）。这些结果的要点是，准确估计重组是产生有关感兴趣统计数据的准确置信区间的关键。

合并基因组和参考基因组

一个物种的参考基因组通常来自对单个个体的测序，无论是近交还是杂交（第 2 章）。因为这个个体没有什么特别之处——它不代表该物种的祖先状态，也不是“共识”个体——参考基因组将包含祖先和衍生状态，以独立分离多态性，就像从自然界采集的任何染色体一样。在短读测序的背景下，参考基因组被用作所有读数映射的主干，因此在种群样本中，它往往发挥着比它应得的更重要的作用。例如，当对替代等位基因的置信度不高时，从短读数据集调用多态性的保守方法有时会默认使用参考碱基。这可能对估计衍生等位基因频谱产生重要影响（请注意，因此默认为 N 更好）。因此，合并谱系的一个有趣（且有用）的用途是问：样本中分离位点的比例是多少，参考基因组中会有衍生等位基因？如果样本大小为 n = 2，答案应该是显而易见的：预计任何两条染色体之间平均有 q 个差异，突变在祖先的每个染色体的谱系中均等发生。因此，两个分支中任一个（任意指定一个作为参考）的衍生等位基因比例将是1/2。

对于大于 2 的样本，需要稍微复杂一些的论证。

实际上有两种方法可以回答这个问题，每种方法都只是对核心问题的重新陈述。为了找到任意选择的染色体上具有衍生等位基因的分离位点的比例，我们可以找到 (1) 样本中衍生等位基因的预期平均频率，

或 (2) 单个染色体相对于整个样本的历史中发生的突变的预期比例。

样本中衍生等位基因的平均频率给出了任何染色体上衍生等位基因的预期比例，平均所有分离位点。为了看到这一点，想象一下从自然界中随机选择一条染色体。对于任何给定的具有衍生等位基因频率 p 的多态性，选择包含衍生等位基因的染色体的概率也将是 p，因此以频率 p = 0.8 采样衍生等位基因的概率将是 0.8，以频率 p = 0.01 采样等位基因的概率将是 0.01，依此类推。因此，对许多分离位点进行平均，采样衍生等位基因的概率将只是平均衍生等位基因频率；该度量将相当于采样衍生等位基因的位点比例。从公式 6.9 中，我们得到了在 i 条染色体上发现的具有衍生等位基因的分离位点的预期数量 Si。现在我们想知道——在所有分离位点中——大小为 n 的样本中每个位点包含衍生等位基因的平均染色体数量，我们将其表示为 m。设 q(m) 为样本中 m 的分布（即衍生等位基因频谱），Griffiths 和 Tavaré (1998) 表明，在平衡状态下，Wright-Fisher 群体中 m 的预期值为：∑ = = − = −Em i i a a n n，使得样本中衍生等位基因的平均频率为 m/n。对于 n = 2 的情况，公式 6.13 再次得出 m = 1，平均衍生等位基因频率为 1/2。对于大小为 n = 11 的样本（如图 6.5 所示），m = 3.59，平均衍生等位基因频率为 0.326。

为了从谱系角度直观地了解这一结果，请考虑另一种表述问题的方式：在整个样本中，单个染色体历史中发生的突变的预期比例是多少？因为我们以看到分离位点为条件，所以我们只需要将染色体的历史追溯到样本的 MRCA；根据定义，在此点之前发生的所有突变都是固定的。

认识到任何单个染色体上的衍生突变都必须在 MRCA 之后发生，我们需要知道追溯到这个祖先的时间。

幸运的是，公式 6.3 给出了从谱系的任何尖端追溯到 MRCA 的预期时间。为了得到只发生在这个分支上的所有突变的比例，我们将公式 6.3 除以公式 6.4，即树中的总时间（图 6.8）。如果突变再次随机地分布在各个分支上，并且我们假设无限位点突变模型，则此关系给出在任何单个染色体上具有衍生等位基因的位点的预期比例 b：等于 m/n，即衍生等位基因的平均频率，如上文公式 6.13 所定义。这些关系表明，平均衍生等位基因的频率较低，对于中小样本，频率范围约为 0.20 到 0.50。从预期的衍生等位基因频率分布（图 6.5）一眼就可以猜到，分离位点的次要等位基因通常是衍生等位基因。事实上，Watterson 和 Guess (1977) 已经证明，等位基因代表祖先状态的概率等于其种群频率 p；样本中等位基因的频率也是如此（Griffiths 和 Tavaré 1998）。

对于参考基因组，上述结果表明，样本中相当一部分多态性位点（0.20–0.50）将具有参考基因组中的衍生等位基因。如上所述，在第 2 章中，这

图 6.8 样本中所有突变的比例发生在单个染色体的历史中。如果突变在树上均匀随机发生，则该比例等于单个染色体自 MRCA 以来的时间量除以树中的总时间。

对于用 * 表示的任意采样染色体，较粗的黑线表示其祖先谱系。该谱系在树中占总时间的比例由公式 6.14 给出。这一发现意味着，在个体样本中，碱基调用可能经常默认为衍生和祖先参考等位基因。如果碱基调用过于保守，则该过程将导致衍生等位基因频谱倾斜，因为多态性ms 与参考中的派生状态将被推向更高的表观采样频率。因此，当碱基调用默认为参考等位基因时，识别变异位点的保守标准可能会更糟。幸运的是，我们可以纠正这种偏差，方法是降低用于调用多态性位点的质量阈值，考虑参考基因组中固有的采样偏差，将低质量碱基分配为 Ns 而不是参考等位基因，或者迭代地识别变化的位点并根据位点与参考不同的更新概率重新调用基因型（例如，DePristo 等人，2011 年）。

选择 7

可以通过分别考虑选择对本身有利、有害或中性的突变的影响（直接选择）以及对与选择密切相关的突变的影响（连锁选择）来有效地划分自然选择对分子变异的影响。在这两种情况下，对多态性模式的期望通常不同，因此不同的方法最适合检测一种或另一种选择。

更深入地理解这些方法至少需要对其背后的理论和假设有一定的了解；因此，在本章中，我将考虑这些假设及其对推断直接选择的影响的影响。 这里讨论的方法已成功用于大量研究，以检查自然选择的影响，是理解 DNA 序列进化的关键工具。

序列分歧的积累

核苷酸替换率

约翰·梅纳德·史密斯 (John Maynard Smith) 称 k = m 是整个进化遗传学中最重要的方程之一 (Maynard Smith 2002)。我完全同意：它既是一个令人惊叹的优雅结果（由 Kimura 1968 年首次应用于分子数据，但也出现在 Wright 1938 年），也是中性理论的典范（在第 1 章中讨论过）。但它不仅是一个优雅的结果——它也是本书其余部分涵盖的许多自然选择测试的基石之一。下面我们开始探索它在分子进化中的作用。变量 k 定义为新等位基因的替换率——即等位基因在长期内固定的速率（有时也表示为 r）。替换率通常以每代或每百万年为单位进行测量。k 的值决定了两个序列预计会随着时间的推移而发散的速度。定义 d 为

两个直系同源序列之间的遗传距离，替换率对预期分歧量的贡献可以用以下公式看出：

E(d) = k2t + qAnc

其中 t 是物种分裂以来的时间（Gillespie 和 Langley 1979；另见 Li 1977）。为简单起见，假设我们计算每个位点的所有值，这样我们就可以忽略序列的长度。因此，在一个基因座内，k

表示跨位点的平均替换率。回想一下，我们使用 2t

是因为替换可以发生在系统发育树的两个分支上，并且我们只考虑序列之间的成对分歧。

我们添加了祖先中两个序列之间预期的核苷酸变异的平均量（qAnc

），因为在物种形成时，这些差异已经沿着两个谱系积累（图 7.1）。这个预期距离与用于统计量 dXY 的预期距离相同（第 5 章）；但是，我们只需从每个物种中抽取一个序列来计算 d。在发散水平远高于祖先物种的预期多态性水平的情况下，我们可以简单地写出：

E（d）≈k2t

这个方程表明了一种明显的关系：物种的不同是因为新的等位基因出现并固定，这个过程发生的速率和自谱系分裂以来经过的时间都对总发散量有贡献。然而，这个方程并没有告诉我们什么影响了替换率。有两个量决定了替换率。第一个是任何突变固定的概率，我们将其表示为 u。这个概率是突变频率 p 的函数，但现在我们只对新突变固定的概率感兴趣。

基因分裂时间

物种分裂时间

图 7.1 物种分裂以来的时间与被比较的基因分裂以来的时间之间的区别。单个祖先物种分裂成两个后代物种的历史由较暗的线条界定。在两个后代物种（A 和 B）中分别采样的单个直系同源序列的历史由较细的蓝线表示。在物种分裂之前，祖先种群中的两个序列预计平均具有 qAnc

d差异。如果祖先种群符合 Wright-Fisher 模型的假设，则基因分裂的时间预计在物种分裂时间 t 之前 2Ne 代。直接选择突变，对于常染色体基因座，其频率始终为 1/2N。当突变对适应度没有影响时（当它相对于替代等位基因为中性时），其固定的概率等于其当前频率。因此，对于 p = 0.4 的突变，u = 0.4，有 40% 的机会它最终会变成频率 1（相反，有 60% 的机会它会丢失）。对于频率为 1/2N 的新的中性突变，其固定概率因此仅为 u0

对于新的有利突变 (s > 0) 和较大的有效种群规模，固定概率为 (Haldane 1927; Fisher 1930; Wright 1931)：

（为了在这些论文中表示选择优势，Haldane 使用符号 k，Fisher 使用符号 a，Wright 使用现在标准的 s。）

是杂合子中新等位基因的选择优势，2sa

是纯合子中的优势。对于没有大影响的新的有害突变 (s < 0)，固定概率为 (Kimura 1957)：

是杂合子中新等位基因的选择劣势，

是纯合子中的劣势。图 7.2 显示了新选定突变相对于中性突变固定的概率，这些概率在 Ns = 0 时相等（即它们的比率为 1）。可以看出，略微有利的突变固定的可能性并不比中性突变高多少，而且即使是略微有害的突变也有一定的固定概率。即使这些略微有害的突变（通常称为近中性）必须具有极小的选择效应才能使 Ns 接近 0，但它们预计会对较小群体中的序列发散做出贡献（见下文）。影响替换率的第二个量是出现并可能被固定的突变总数。如果每一代核苷酸发生突变的概率为 n，那么在 N 个二倍体个体的群体中，在常染色体位点上，总群体中每代将有 2Nn 个新突变（在单个位点）。在这些突变中，2Nnf0（=2Nm）将是中性的，其中 f0 代表中性突变的比例。其余突变由

有利的部分 fa

有害的部分 fd

组成（因此

图 7.2 相对于中性等位基因的固定概率，新的选定突变。中性突变的固定概率为

1/2N，有利（Ns > 0）和有害（Ns < 0）突变分别具有增加和减少的固定概率。

值是根据公式 7.4 计算得出的，可用于正向或负向选择系数。

134 第 7 章

= 1）。因此，每代也有 2Nnfa

个新的有利突变

个新的有害突变。

我们现在可以将替代率的两个决定因素放在一起考虑。如果没有有利突变，并且有害突变

未在可察觉的水平上固定，则替换率仅取决于

发生的中性突变总数以及每个突变固定的概率。我们可以将此关系写为：

ν ν ( ) = 

  = k N f N f 2 1

2 0 0 (7.5)

或者，用符号 m 代替中性突变发生的总速率

nf0

，中性突变的替换率为：

k = m (7.6)

换句话说，此关系表示，当仅考虑中性突变时，

替换率等于突变率，无论种群大小如何。这种与种群规模无关的情况乍一看似乎不直观，但可以通过考虑以下事实来理解：虽然在大群体中会出现更多突变，但每个突变最终固定下来的可能性较小。同样，任何单个新突变在小群体中固定下来的可能性更大，但总体而言，突变较少。因此，在考虑中性突变时，种群规模完全抵消了这一影响。在我们探索自然选择发挥作用的案例的过程中，我们将一次又一次地使用这一结果。最后，让我们考虑非中性等位基因的替换率。在总突变率 n 中，我们已经将 fa 定义为所有有利和有害突变的比例，将 ua 定义为每种类型新突变的固定概率。现在，我们可以将有利突变的替换率写为：

k = (2Nn fa

) = 4Nn fa

为了获得有害突变的替换率，我们使用公式 7.4 来查找：

=

−()()−−

从公式 7.7 和 7.8 中，我们可以看到一个关键的区别中性突变和选择性突变之间的区别：N，即种群规模，在选择性突变的替换率中起着重要作用。在其他条件相同的情况下，更有利的突变将在较大的种群中固定，而不是在较小的种群中，而更有害的突变将在较小的种群中固定，而不是在较大的种群中。接下来，我们考虑如何确定每种类型的突变对物种间分化的贡献。

估计序列分化

量化自然选择影响的最直接方法之一根本不使用种群水平的数据。虽然这种直接选择特性似乎不适合用于种群遗传学文本，但它是每个种群遗传学家工具箱中的标准测试。

而且由于物种间的替换是物种内多态性的长期结果（Kimura 和 Ohta 1971），这些方法是理解使用多态性数据的更强大方法的良好起点。

之前我们将 d 定义为两个直系同源序列之间的遗传距离。具体来说，我们通常通过从每个物种中取出一个序列并计算它们之间不同的位置数，然后除以对齐的核苷酸总数来计算 d。与计算同一物种中个体之间的 p（第 3 章）一样，我们只计算没有间隙的位置上的核苷酸差异数。如果我们考虑多序列比对，这个过程通常会导致比对中任何序列中有间隙的位置在所有物种中都被忽略。因此，从为每对物种生成的比对计算出的成对发散值可能与从包括所有物种的单个多序列比对计算出的发散值略有不同。例如，图 7.3 显示了来自不同物种的 15 个核苷酸长的序列比对。使用整个多物种比对，我们将计算出序列 1 和 2 之间的发散为 0.273（3 个差异/11 个无间隙位点）；但是，如果我们只重新比对序列 1 和 2 并重新计算发散，我们将得到 d = 0.23。显然，如果在同一比对中同时包含近亲和远亲序列，则这种差异可能导致发散估计有偏差，因为只有那些不易插入和删除的区域才会被纳入分析。计算 d 时还有两个重要问题需要考虑。

首先，由于它是根据每个物种（或每个基因座，因为 d 也可以在旁系同源物之间计算）的单个序列计算得出的，因此该统计数据本身并不代表固定的差异。也就是说，在用于计算 d 的序列中存在的任何衍生的多态性等位基因都将被计入种间差异。之前我们还考虑了祖先多态性等位基因的贡献，但这些基因实际上在当前物种之间是固定的，因此不会混淆多态性和差异。这种差异度量包括多态性位点，因为从自然界采样的任何染色体预计平均包含 q 个在种群中仍然具有多态性的衍生等位基因。您可以通过考虑种群最近共同祖先的预计时间为 4N 代（公式 6.3）来了解这种预期来自何处。

这意味着在导致任何序列的谱系中发生了 4Nm 个突变（并且仍然是多态性的）（图 7.4）。假设被比较的两个物种的多样性水平相同，这意味着直系同源序列之间大约 2p 的差异实际上不是物种之间的固定差异。这些 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

T T ACA A T CCGA TCGT

图 7.3 四个序列的示例比对，每个序列来自不同的物种

物种间基因分裂时间

物种分裂时间

物种内基因分裂时间

图 7.4 来自两个不同物种的单个序列的比较包括固定和多态性差异。单个祖先物种分裂成两个后代物种的历史由较暗的线条界定。在两个后代物种中分别采样的单个直系同源序列的历史用较细的蓝线表示；假设序列 A 是从物种 1 中采样的，而序列 C 是从物种 2 中采样的。为清楚起见，显示了每个物种中两个最大发散样本的谱系（物种 1 中的 A 和 B；物种 2 中的 C 和 D）。每个物种内样本的最近共同祖先的时间不一定相同，但为简单起见假设相同。从物种间基因分裂时间到物种内基因分裂时间发生的所有变化都是固定的差异。自物种内基因分裂以来，导致 A 和 C 的单个分支上发生的突变是多态性的。单个序列中此类差异的数量预计为 A 和 C 的 p，因此 A 和 C 之间的差异总数 d 包括 2p ) 个多态性位点。计算意味着我们可以正确地将两个物种之间的固定差异量计算为 d* = d − 2p。这种方法类似于第 5 章中用于计算统计 da 的逻辑，但在这种情况下，我们试图消除当前的多态性而不是祖先的多态性。如果多样性水平不相等，我们可以从每个物种中减去不同的 p 值，这样 d* = d − (p1

大多数情况下，这种校正被忽略，因为分歧水平远大于多态性水平（即 d >> p）。

第二个问题出现在被比较的物种距离太远的情况下。早期的分子进化实践者认识到，当多个突变出现并固定在同一位点时，简单地计算两个物种之间不同位点的比例会低估真正的差异数量。在最极端的情况下，当每个位点都发生了多次变化时，我们仍然会有 25% 的碱基随机匹配（假设没有间隙）并且 d = 0.75。解释多次替换影响的最简单方法是进行 Jukes-Cantor 校正（Jukes 和 Cantor 1969）。将 a 定义为观察到的两个序列之间不同的位点比例（也称为 p 或 p 距离），Jukes

Cantor 校正的 d 值为：

   d a 3

Kimura (1980) 进行了调整，以解释转换突变（嘌呤之间或嘧啶之间）比颠换突变（嘌呤和嘧啶之间）更频繁发生的事实。 Kimura 双参数校正为：

( ) ( ) =− − − − − − d a b b 1

2 ln 1 2 1

4 ln 1 2 (7.10)

其中 a 现在是因转换而不同的位点比例，b 是因颠换而不同的位点比例。可以进行许多进一步的阐述，例如允许嘌呤与嘧啶的转换使用不同的参数（例如，Hasegawa、Kishino 和 Yano 1985；Tamura 和 Nei 1993）。当距离非常大时，这些方法都不能确保准确测量发散度，因为在这种情况下，即使是修正后的估计值也存在很大的方差。大于 0.75 的 a 值在公式 7.9 和 7.10 中均未定义，并且通常，大于 0.50 的 d 值将不如理想值准确。远大于 1 的 d 值基本上是对真实发散水平的猜测。利用分歧检测选择

除了测量两个序列之间的整体分歧之外，

在编码区域中，我们还可以分别测量由于

非同义和同义变化而导致的分歧。具体来说，我们定义 dN 为每个非同义位点的非同义差异数（也称为 dR ），每个同义位点的同义差异数（也称为 KS ）。计算两个序列之间非同义和同义差异数的方法在 BOX 7.1 和 BOX 7.2 中有详细描述。这里我们重点介绍如何使用量化来推断自然选择。计算 dn 的最简单方法首先找到比对中序列之间的非同义和同义差异数（分别表示为 Nd ），然后计算每个序列中的非同义和同义位点数（分别表示为 ）。最后，我们将每个序列的非同义差异数（继续）■ BOX 7.1（续）

非同义位点数 dn

每个同义位点的同义差异数 dS

。由于

的值对于比对中的每个序列都会略有不同，因此我们通常在计算 dn 时使用它们的平均值。有许多不同的方法可以计算非同义和同义位点，并且使用不同的方法通常会得到不同的答案。一种简单的方法假设转换和颠换的速率相等，并且没有密码子使用偏差（nei 和 Gojobori 1986）。例如，密码子 TTA 编码氨基酸亮氨酸。如果我们只考虑第三个位置的突变，TTG 亮氨酸 TTA 亮氨酸 → TTC 苯丙氨酸 TTT 苯丙氨酸 使用 nei 和 Gojobori 方法计算位点，第三个位置将是 1/3 同义和 2/3 非同义，因为三个可能的突变中只有一个是同义密码子，两个是非同义密码子；我们将此类密码子称为双重退化。密码子中的前两个位置都是 3/3 非同义的nd

0/3 同义，因为这两个位置上的所有突变都是非同义的。因此，该密码子的总计数为

= 8/3 和 Sc

其他方法在计算有效位点数时会考虑转换/颠换偏差（li 1993；Pamilo 和 Bianchi

1993；Comeron 1995；Ina 1995）。这种偏差可能对位点计数产生重要影响。例如，强烈的转换/颠换偏差将导致双重简并密码子在第三个位置上的同义性超过 1/3，因为遗传密码允许双重简并密码子在第三个位置上的所有转换都是同义的。从上面的例子可以看出，如果 A → G 的突变比 A → C 或 A → T 的突变更频繁，那么从某种意义上说，第三个位置的同义性超过 1/3。适当地加权这些位点（使用已知或估计的转换与颠换的比率）可以纠正这种偏差。还有其他方法考虑了序列内不平等的密码子使用情况，有效地纠正了一些同义密码子很少使用的事实（Yang and nielsen 2000）。还有许多不同的方法来计算序列之间的非同义和同义差异。当密码子仅在一个核苷酸处不同时，所有这些方法都会给出相同的结果。例如，在以下两个物种之间单个密码子的比对中：只有一个差异，并且它是同义的；因此 Nd

= 0 和 Sd

只有当两个序列之间单个密码子存在多个差异时，情况才会变得复杂；我们将这种情况称为复杂密码子。这些情况之所以更加复杂，是因为现在进化可能采取了两种可能的途径，具体取决于两种突变中的哪一种先出现。在某些情况下，核苷酸差异可能是同义的，也可能是非同义的，如下所示

两种可能的途径是：

途径 1 TTT (Phe) ↔ GTT (Val) ↔ GTA (Val)

途径 2 TTT (Phe) ↔ TTA (leu) ↔ GTA (Val)

在途径 1 中，第三位置的 T → A 变化是同义变化，总共有一个非同义差异和一个同义差异 (Nd

= 1)。在途径 2 中，T → A 变化是非同义的，总共有两个非同义替换 (Nd

= 0)。nei 和 Gojobori (1986) 方法

只是对不同的进化路径进行平均，并赋予每个路径等效权重；这样做会得到

= 1.5 和 Sd

= 0.5。其他方法会给予包含同义

直接选择 139 的路径更多权重 作为框 7.1 中概述的计算 dn 的“计数”方法的替代方法，已提出了几种可能性方法（并且通常用于此目的）（Goldman 和 Yang 1994；Muse

和 Gaut 1994）。这些方法基于密码子之间替换的马尔可夫链模型，指定 61 个非终止密码子之间独立转换的瞬时速率矩阵 Q。对于最流行的似然模型之一（Yang 1998），速率矩阵由每对密码子 i 和 j（i ≠ j）之间的单独转换速率 qij 组成，其中：这里，密码子 j 的平衡频率表示为 pj ，转换/颠换偏差表示为 k ，非同义替换相对于同义替换的相对概率表示为 w 。在这种情况下，我们可以将 w 视为非同义替换与同义替换的比率，或 dn ，该方法明确不是单独估计。事实上，虽然实现该模型的方法可以允许 w 在密码子之间变化，但它们假设同义替换率在整个基因中是恒定的 - 因此，推断密码子之间的差异仅仅是由于对非同义变化的选择差异。但是，可以放宽简单模型的许多假设：可以允许密码子之间有多个变化（例如，Whelan 和 Goldman 2004），可以允许密码子之间的相关变化率（例如，Siepel 和 Haussler 2004a），并且可以允许密码子之间的同义替换率变化（例如，Pond 和 Muse 2005）。当与多序列比对和预先指定的系统发育拓扑一起使用时，可以通过最大似然法估计上述模型的参数。

最常用的软件包是 PAML (Yang 2007) 和 HyPhy

          

         

i 0，如果 和 相差不止一个位置

，对于同义颠换

，对于同义转换

，对于非同义颠换

，对于非同义转换

替换 (li, Wu, and luo 1985) 或包含常见非同义变化的途径 (Miyata and Yasunaga 1980)，然后对途径取加权平均值

。除非序列高度分化，否则所有这些计算非同义和同义差异的方法应该会给出相似的结果。

最后，在对比对中的所有密码子的 Nc 求和后，建议对单个位点的多个可能的核苷酸替换的两个 dn 进行校正。

公式 7.9 和 7.10 描述了执行此类校正的两种常用方法。有关所有这些方法的更多详细信息，请参阅 nei 和 Kumar (2000) 和 Yang (2006)。

（续）

140 第 7 章

（Pond、Frost 和 Muse 2005）；在这里我描述了 PAML 中可用的方法。要计算所有密码子的平均自然选择强度 w，最简单的模型表示为 M0。

该模型对 w 的取值没有任何限制，如果需要，可以为树的每个分支估计单独的 w 参数。但是，除非正向选择作用于整个基因的许多位点，否则 M0 不太可能有 w > 1。相反，允许 w 在密码子之间变化的模型（称为位点模型）对于检测正向选择最有效。要检验比对中是否存在有证据显示存在正向选择的密码子，最直接的方法是进行两个嵌套位点模型之间的似然比检验，一个模型带有一个允许 w > 1 的参数，另一个模型不带有。模型 M1a 有两个 w 参数：一个用于 w < 1 的密码子 (w0

一个用于 w = 1 的密码子 (w1

；还有用于分配给每个 w 参数的位点比例的参数，p0

）。模型 M2a 有三个 w 参数：一个用于 w 的密码子

），一个用于 w = 1 的密码子 (w1

一个用于 w > 1 的密码子 (w2

)，以及用于分配给每个 w 的位点比例的相应参数 (p0

因为模型 M1a 有两个自由参数，而模型 M2a 有四个自由参数——并且模型 M1a 是嵌套的在 M2a 内——可以使用具有两个自由度的似然比检验来测试允许位点在正向选择下进化的模型 (M2a) 是否能更好地拟合数据。此检验只能用于两个分类单元（即两个序列），但添加更多分类单元时功效会提高（Anisimova、Bielawski 和 Yang 2001）。但请注意，该方法假设底层系统发育树来自不同物种的序列，而不是物种内的序列；PAML 要求比对中的所有位点都有一棵树，而从物种内采样重组序列会违反这一假设，并可能导致错误的正向选择推断（Anisimova、nielsen 和 Yang 2003）。对位点模型的进一步细化允许对系统发育树特定分支上的一组密码子进行正向选择测试（称为分支位点模型），这反过来又允许识别在一组位点上经历适应性进化的特定谱系（Zhang、nielsen 和 Yang 2005）。有关这些方法及其用途的更多详细信息，请参阅 Yang（2006；2007）。从替代率的讨论中应该可以明显看出，自然选择对固定的非同义突变的数量有着深远的影响。如果我们使用中性、有利和有害突变的替换率期望（分别在公式 7.6、7.7 和 7.8 中给出），并将所有替换率乘以 2t 以将其转换为距离，那么我们可以写下总预期非同义突变

ν ν ν ( ) ( ) ( )( ) ( ) = 

     ( ) −

E d N f N t N f s t N f s

2 2 2 2 2 2 2

2 Ns N 0 a a d

简化为：

ν ( )= + +

     ( ) −

E d t f f Ns f Ns

1 Ns N 0 a a d

因为一个区域的总非同义变异是由所有三种类型的突变引起的，所以我们对 dN

的表达式包括所有三个项。这个等式

■ 框 7.2（续）

直接选择

表明，更高的潜在突变率 n 和更长的变异时间 t 将增加变异量；固定的有利突变的比例将取决于它们出现的频率和它们的平均选择效果；如果选择足够弱，有害突变也会导致变异。

虽然并非所有同义突变都是选择性等价的，但我们现在将做出这一假设，稍后再讨论选择对同义突变的影响。假设所有同义变化

都是中性的，即这些位点的 f0

= 1 和 fa

= 0，那么两个序列之间预期的同义差异总量为：

该表达式反映了公式 7.2，因为对于中性突变，替换率等于突变率。

为了询问自然选择对序列进化的影响，

我们不能只比较基因之间的 dN

值，即使它们都来自

同一对物种。这是因为整个基因组中潜在的突变率各不相同e 和染色体上的差异导致在选择可能以完全相同的方式起作用的基因座上产生不同的值。Kimura (1977) 非正式地建议比较基因上非同义和同义分歧的比率 dN ，以控制基因座间突变率的差异（另见 Miyata 和 Yasunaga 1980）。因为 n 和 t 对于同一基因中的非同义和同义位点都大致相同，所以将公式 7.12 除以 7.13 可得出：

E d f f Ns f Ns

这个表达式向我们表明，相对于同义分歧，非同义分歧的水平再次归因于中性、有利和有害突变的比例。请注意，为了避免下标，我没有使用不同的符号来表示同义和非同义位点的中性突变比例；公式 7.14 中的 f0 仅代表非同义突变。

解释 dN

虽然 dN

由中性、有利和有害突变的综合影响决定，但它可以告诉我们很多有关自然选择对序列进化的一般影响的信息。图 7.5 给出了自然界中观察到的范围值的示例，其中显示了拟南芥和 A. lyrata 之间的 11,492 个直系同源物的值（Yang 和 Gaut 2011）。作为对许多分类群的一般观察，物种间计算的

平均值为 0.15 到 0.25，尽管单个基因的值从 0 到 >2 不等（大鼠基因组测序计划联盟

2004；Mikkelsen 等人 2005；果蝇 12 基因组联盟 2007；Yang

和 Gaut 2011）。该范围表明至少 75% 到 85% 的非同义突变是有害的并且不会固定。这代表了对强烈有害突变比例的最低估计，原因是

图 7.5 拟南芥和 A. lyrata 之间 11,492 个直系同源物的 dn 值分布。 （根据 Yang 和 Gaut 2011 年的研究。）即使少量的有利突变也会对分化产生不成比例的贡献（见公式 7.14）。即使有利突变的比例 fa 低至 1%，并且没有中性突变，足够强的选择和较大的种群规模也可能使 dN 远远超过 0.01。例如，如果 99% 的非同义突变是强烈有害的，而 1% 是有利的，平均 Ns = 10，则等于 0.4。对有害突变比例的简单估计 1 − dN 会低估该比例约 50%（0.6 对 0.99）。仅当不存在有利突变且有害突变不会以可观的速率固定时，dN​​ 才会固定（因此为 1 − dN ）。虽然基于此假设的选择性约束估计值被广泛使用，但它们几乎肯定低估了所有有害突变的比例。单个基因的 dN 值并不能告诉我们固定的非同义突变中有多少比例是中性的、有利的或略有害的。每种突变频率、其平均选择效应和种群大小的多种值组合可以给出该统计数据的相同值。一种用于估计轻微有害突变比例的方法是比较已知种群规模不同的分类群之间的 dN 平均值：预测种群规模较小的物种的 dN 将增加，与轻微有害突变的比例成正比。

对多个物种的比较表明，所有突变中有 15% 到 30% 是轻微有害的 (Eyre-Walker 等人，2002 年)。例如，小鼠和大鼠之间的平均值为 0.14，人类和黑猩猩之间的平均值为 0.20 (Mikkelsen 等人，2005 年)。虽然没有先验理由相信灵长类动物总体上存在更多和/或更强的有利突变，但众所周知，该进化枝中的种群规模要小得多，这将导致更高的 dN 值

直接选择

只有一种普遍接受的方法可以证明 - 使用 - 基因历史上存在有利突变（即正向选择）：该统计量的值必须大于 1（Hill 和 Hastie 1987；Hughes 和 Nei 1988）。要了解原因，请考虑一个没有有利突变的基因。在这种情况下，dN 的最大值为 1，只有当所有可能的非同义突变都是选择性等价时才会发生这种情况。获得 dN

> 1（具有统计意义）的唯一方法是存在一定比例的有利突变（公式 7.14）。我们还可以看到，这是一个非常严格的要求，因为需要越来越多的有利突变才能获得 dN

，有害突变的比例也在增加。事实上，只有那些存在重复的氨基酸替换适应性固定的基因才有机会提高取 dN

的值除以 1。这些基因通常是与免疫和生殖功能相关的蛋白质，因为它们经常参与宿主和病原体之间或雄性和雌性之间的“军备竞赛”动态。

dN

> 1 的显著性检验已通过多种不同的方式完成。框 7.2 中概述的似然法允许在有和没有正向选择的模型之间进行似然比检验。框 7.1 中概述的更简单的计数方法不基于任何模型，因此不允许进行如此简单的测试程序。

然而，后一种方法已经产生了许多不同的方法。Hughes 和 Nei (1988) 仅比较序列对，计算了 dN

的标准误差，并使用 t 检验来检验 dN

> 0 的假设是否显著。 Messier 和 Stewart (1997) 对灵长类动物溶菌酶基因进化树的各个分支进行了类似的分析，重建了每个节点的祖先状态，并计算了父节点和子节点之间非同义和同义差异以及非同义和同义位点的数量；同样，测试通过 t 检验进行。这两种方法都存在问题。首先，t 检验假设 dN 近似呈正态分布，除非两种类型的替换都存在大量替换，否则这种情况不太可能发生。其次，使用重建祖先状态的方法未能考虑这些重建中的不确定性，从而导致不确定性。第一个问题的一个解决方案是使用小样本测试来比较非同义和同义替换的数量。 Zhang、Kumar 和 Nei (1997) 通过使用 2 × 2 独立性检验得出了这一解决方案，该检验对非同义和同义替换的计数与非同义和同义位点的数量（即框 7.1 中的 Nd）进行了独立性检验；通过 Fisher 精确检验进行测试。这种方法优于大样本近似法，尽管计数非同义和同义位点的方法的选择可能会对结果产生很大影响（见框 7.1）。似然法（框 7.2）可用于对重建的祖先状态进行平均，从而有效地考虑到不确定性。另一个复杂因素是蛋白质的不同部分可能受到不同的选择力——一些可能受到强烈的正选择，而一些则受到强烈的负选择。虽然所有区域的平均值可能具有 dN < 1，但在 Hill 和 Hastie [1987] 和 Hughes 和 Nei [1988] 的原始论文中，功能域的子集可能具有 dN。因此，如果提前确定了这些域，则在不同的域中检查此统计数据会很有用；事后滑动窗口分析可能会产生很大的误导性（Schmid 和 Yang 2008）。只要满足测试假设，似然法也可用于以统计上合理的方式查找具有 dN > 1 的单个密码子（框 7.2）。许多或大多数具有 dN < 1 的蛋白质很可能具有一些正向选择的历史，但这个统计数据不允许我们确定这是否属实。任何数量的氨基酸变化都可能是由于有利突变引起的，因此任何 0 < dN ≤ 1 的基因都可能具有中性、有利和轻微有害突变的贡献。由于每种突变类型的贡献未知，我们对选择的推断必须同样谨慎。以下是解释 dN 的一些一般准则

<< 1 绝大多数非同义突变都是有害的，并且负向（净化）选择占主导地位。

< 1 大多数非同义突变是有害的，但可能存在一些未知的有利突变部分。

= 1 这种情况可能发生在两种情况下：首先，没有选择，所有非同义突变都是中性的。其次，存在大量中性和有利突变（以及有害突变）。

> 1 有许多有利的非同义突变，并且正向选择占主导地位，但仍有许多有害突变。

尽管第 1 章已经进行了区分，但应在此重申，观察到基因的 dN

= 1 并不表示该基因正在“中性”进化。如上所述，有多种选择力的组合会导致 dN

= 1，其中只有一种是没有选择。更重要的是，没有选择约束与中性进化不同。许多 dN

< 1 的基因可能正在中性进化 - 也就是说，它们可能没有固定任何适应性替代。

还必须提到有关解释 dN

的几个注意事项。

首先，较小的 dS

值可能会使 dN

的估计值向上偏差。虽然这种偏差绝不会导致 dN

的值明显大于大于 1，由于 dS 的差异，对具有非常不同分化时间的物种对进行比较可能会产生误导（Wolf 等人，2009 年）。其次，我们假设所有同义突变都是中性的。显然，并非每个物种都是如此：虽然非随机使用同义密码子的原因有很多，但自然选择是一个重要因素（Ikemura，1981 年；Akashi，1994 年；Chamary、Parmley 和 Hurst，2006 年；Lawrie 等人，2013 年）。如果所有中性的同义突变的比例不等于 1，则 dS 的值应由公式 7.12 而不是 7.13 中给出的表达式确定。目前尚不清楚这将如何影响：如果负选择占主导地位，那么 dS

将小于没有约束的情况下的预期，因此 dN

将更高。这是否会导致在没有有利突变的情况下 dN

> 1 的情况（这实际上不需要对非同义突变进行约束，但需要对同义突变进行严格约束）尚不清楚（但可能性不大）。另一方面，如果存在许多有利的同义替换（例如，Resch 等人，2007 年），那么 dN

的值将低于预期，因此将推断负选择的强度大于实际强度。一般而言，dN

是一个相对稳健的统计数据，几乎没有关于样本收集方式、人口人口历史或真正群体遗传方法引入的任何其他问题的假设。另一方面，只有反复的正向选择才能提供有利突变的证据，因此只有更强大的方法才能为适应性进化提供大量信息。

使用多态性检测选择

选择对多态性频率的影响

在上一节中，我们考虑了选择如何反对或支持所选等位基因的最终固定，以及我们如何从物种间分化数据推断进化过程。我们现在将讨论选择对种群内多态性频率的影响，最初不考虑固定差异。许多不同的方法可以利用这种数据，尽管它们通常不如结合多态性和分化的方法那么强大。在本章后面，我们将看到一种这样的方法，即 McDonald-Kreitman 检验。在比突变固定所需时间更短的时间尺度上，选择将改变种群中等位基因的平均频率。正如从上面讨论的不同固定概率所预期的那样，有利等位基因的平均频率将高于中性等位基因，而有害等位基因的频率则相对较低。对于新的突变，在频率 q 处发现的多态性密度由以下公式给出（Wright 1969，第 381 页）：

ν ( ) ( ) = −

其中 n 再次是突变的总速率。有利突变的 s > 0，而有害突变的 s < 0。

有利、中性和有害等位基因的频谱示例如图 7.6 所示。从分布中可以看出，有利的等位基因向更高的频率移动，而有害的等位基因向更低的频率移动。平衡选择将使等位基因保持在中间频率，平衡选择的确切形式极大地影响等位基因频率的特定概率密度。分布

多态性位点数

衍生等位基因的频率

图 7.6 有利（Ns > 0）、中性（Ns = 0）和有害（Ns < 0）突变的等位基因频谱。（Fay、Wyckoff 和 Wu 2001 年）

非同义

同义

图 7.7 弱有害非同义多态性存在的检验（P = 0.001；Fisher 精确检验）。多态性是指样本中出现不止一次的多态性；单态是仅出现一次的多态性。数据来自大肠杆菌中的 gnd 基因。（根据 Sawyer、Dykhuizen 和 Hartl 1987）。对于具有相对较强选择的变体（例如，|Ns| > 1），有害和有利等位基因的偏差将非常大。最强的有害突变会立即从种群中移除，而最强的有利突变会非常迅速地固定（或遭受快速的随机损失）；两者都很少处于多态状态。这些趋势意味着，在只有非常有害的突变和中性突变的位点，在中等大小的样本中观察到的等位基因频谱将更接近中性频谱。 Sawyer、

Dykhuizen 和 Hartl (1987) 提出了一种非同义多态性的中性简单测试方法，该方法基于

轻微有害和轻微有害之间的等位基因频谱差异。d 中性变异。

Sawyer 及其同事统计了单态非同义多态性的数量（即仅存在于一条染色体上的多态性）和存在于多条染色体上的多态性的数量（“多态性”）。

他们的样本（大肠杆菌中的基因 gnd）中的所有 12 个非同义多态性都是单态。然后，他们将这两个数字（12 个单态 vs. 0 个多态性）与同义多态性的配置进行了比较：34 个同义多态性是单态的，32 个是多态的（图 7.7）。在 2×2 列联表中比较这些值具有高度显著性，表明单态非同义变异过多。如果非同义多态性平均而言略有危害，这一发现正是人们所期望的。

除了简单地表明等位基因频率因选择作用于多态性而不同之外，人们还可以使用公式 7.15 并仅做一些假设（例如位点之间的自由重组）来找到与一组多态性相关的平均选择系数（例如，Sawyer 和 Hartl 1992）。也就是说，我们可以采用观察到的大量多态性的等位基因频率谱来估计公式 7.15（或类似方程）的选择参数。将选择参数本身视为分布也很有用，这样不仅可以估计选择的平均效应，还可以估计有利和有害突变的相对贡献。直接选择的等位基因频率谱的预期假设了一个平衡种群，因此任何非平衡种群人口统计甚至连锁选择都会混淆推论，因为两者都会影响频率谱（见第 8 章和第 10 章）。由于这些混杂因素，推断非同义位点选择参数的常用方法也使用同义位点或内含子位点来同时控制非平衡过程（例如，Loewe 等人 2006 年；Keightley 和 Eyre-Walker 2007 年；Boyko 等人 2008 年）。图 7.8 显示了 40 条人类染色体上 301 个基因的非编码、同义和非同义多态性的展开等位基因频谱。低频非同义多态性大量过剩。所有类型的高频多态性略有过剩，这归因于祖先状态的错误指定。（Williamson 等人，2005 年；版权所有 2005，美国国家科学院。）基因（Williamson 等人，2005 年）。在这种情况下，所有三类突变都与中性平衡预期不同，但可以通过对类别进行比较推断出过多的有害非同义多态性。

使用这些方法准确推断选择参数需要来自大量位点和较少染色体（大约 n = 10）的数据（Keightley 和 Eyre-Walker 2010）。有些方法利用固定差异（即 p = 1.0 的频率类别），有些则不利用，但并非所有使用固定差异的方法都利用了多态性的全频谱（例如，Bustamante 等人 2002）。使用固定差异可以提高检测出非常有利的等位基因的能力，否则这些等位基因不会被检测为多态性。然而，一些使用固定差异的方法容易得出高度偏斜的 Ns 值，这可能是因为必须对选择强度分布的形状做出假设。尽管如此，在介绍 McDonald-Kreitman 检验之后，我们将在本章末尾讨论几种此类方法。直接选择对多态性的平均频率有着深远的影响。虽然估计平均选择系数（或更准确地说，选择系数和种群大小的乘积）通常需要来自多个基因座的数据，但可以测试这些值是否不同于 0，或者选择系数在多态性类别之间是否不同（例如，非同义与同义）。虽然这些信息并没有告诉我们某个特定基因的历史，但它确实提供了选择影响的全局视图。在下一节中，我将介绍基于此处描述的期望的附加统计测试，这些测试将使我们能够了解特定基因座的历史。

通过类比本章前面列出的 dN 比较逻辑，我们可以在一个物种内比较每个非同义位点的非同义差异的平均数量（称为 pN）与每个同义位点的同义差异的平均数量(pS

pN

的计算很简单，结合了计算 p 的方法（第 3 章）和框 7.1 中描述的计算非同义和同义变化和位点的方法。只需要来自单个种群或物种内的序列。在大多数情况下，计算通常比物种间比较更容易，因为有限的多态性数量也会限制必须考虑的复杂密码子的数量。

pN

比率的解释也与 dN 大致相同，但有一个重要的例外。低于 1 的 pN

值再次证明了净化选择占主导地位，绝大多数编码位点都小于 1。可以对一些稍微有害的多态性进行采样，但由于它们的频率较低，因此对平均杂合性的贡献相对较小。当值大于 1 时，解释 pN

的主要区别是：因为正向选择会迅速固定有利突变，但在多态性研究中很少发现这些适应性变化。即使在直接选择固定的过程中采样了一两个这样的突变，pN 的增加也不太可能导致 pN 的值远高于 1。相反，我们发现在多等位基因平衡选择下的基因中 pN > 1——也就是说，在大量氨基酸多态性由某种形式的平衡选择维持的情况下，例如杂合子优势或负频率依赖选择。很少有基因似乎受到这种选择，主要组织相容性复合体 (MHC) 基因和其他免疫系统相关或性别和生殖相关基因是最突出的例外。Hughes 和 Nei (1988) 首次计算了人类和小鼠 MHC 基因座并发现值大于 1；令人困惑的是，他们将其统计数据称为 dN

正如要求 dN

大于 1 才能证明正向选择具有强有力的证据一样，pN

> 1 是检测平衡选择的一个非常严格的标准。在非常强的选择下，单个位点永远不会对 pN

的值做出足够的贡献，使 pN

大于 1，因此使用这种方法只能检测到极少数多等位基因选择的情况。即使在这些情况下，平衡选择的证据也可能仅限于已知的功能区域（如 Hughes 和 Nei 1988 年所述）。尽管如此，pN

的比较是种群遗传学家工具箱中的标准方法，接下来我们将看到如何将这些统计数据的细微修改与物种间差异相结合，以创建真正强大的自然选择测试。使用多态性和发散性检测选择

McDonald-Kreitman 检验

我们已经看到，通过比较 dN，我们可以从理论和统计上严谨的角度检测出正向选择的作用。但是，由于大约三分之二的可能编码突变会导致氨基酸发生变化，因此，只有当物种间所有差异中有三分之二以上是非同义的时，dN 才会大于 1。发生这种情况的选择性条件可能仅限于存在大量非同义差异的反复适应性固定的条件。在本节中，我将介绍由 McDonald 和 Kreitman (1991) 首次提出的一种检验方法，它类似于 dN，但可以原则性地绕过检测正向选择的非常严格的条件。该测试结合了多态性和发散性，以测试中性模型的预测。无论选择的强度如何，新的测试仍然无法仅通过一个或几个固定的差异来检测适应性进化——要做到这一点，我们将不得不等待下一章介绍的方法。但 McDonald-Kreitman (MK) 测试提供了我们拥有的最强大、最可靠的检测自然选择作用的方法之一。我首先解释测试是如何进行的，然后是其背后的推理和结果的解释。 McDonald-Kreitman 检验的实验设计

对单个基因进行 MK 检验所需的四个量是

非同义多态性的数量（PN

）、同义多态性的数量（PS

）、非同义固定差异的数量（DN

多态性

非同义

同义

图 7.9 McDonald-Kreitman 检验。

所需的四个计数是

非同义固定差异的数量（Dn

同义固定差异的数量（DS

非同义多态性的数量

）和同义多态性的数量（PS

和同义固定差异的数量

）。这些值放在 2 × 2 列联表中（图 7.9），以及标准独立性检验——当计数较低时通常是 Fisher 精确检验，但 c2 或

大样本近似的 G 检验——可用于计算显著性。

用于计算显著性的一般实验设计确定这四个量需要从单个物种中取样多个染色体，以及从第二个物种中取样单个染色体或多个染色体。当从多个物种中取样多个染色体时，可以合并多态性计数（如 McDonald 和 Kreitman 对果蝇、拟果蝇和果蝇的 Adh 基因的原始研究），或者可以对每个物种分别进行单独的测试，特别是当物种之间的 PN 比率存在显着差异时。一般来说，MK 测试对于从真正的“种群”中取样染色体相当稳健：从一个物种中的多个亚种群中取样一个染色体不会导致错误的正向选择推断，尽管它可能导致其他错误的结论（见下文）。 MK 测试中使用的所有值都是计数——无需计算非同义或同义位点的数量。虽然当同一密码子中同时发生非同义和同义变化时，仍可能存在一些涉及复杂密码子的问题，但可以使用框 7.1 中概述的相同方法解决这些问题（尽管 MK 测试中仅使用整数值）。使用标准 MK 测试时，与样本量有关的确定偏差问题很少，因为不使用等位基因频率，但显然必须确定所有突变，而不管它们是非同义的还是同义的（基因分型平台通常会丰富非同义变体）。事实上，只要样本中存在足够的多态性和固定差异，人们实际上只需要从一个物种中取两条染色体和一个外群中取一条染色体就可以进行测试。多少才足够？除非行计数总和和列计数总和均为 4 或更大（即，无论突变类型如何，至少存在四个多态性和四个固定差异），否则在 2 × 2 测试中不可能获得小于 0.05 的 P 值。因此，应继续进行采样——根据需要从种群中或更远的外群中收集更多染色体——以增加统计能力。否则，将无法获得具有统计意义的结果。 MK 检验中散度计算的一个重要区别是，在 McDonald-Kreitman 检验中，我们明确计算样本中 DN 值的直接选择固定差异。在这种情况下，固定差异意味着在种群样本中发现的等位基因的身份与在外群中发现的核苷酸不同。这与计算 dN 时使用的术语散度含义略有不同。如果种群样本中某个位置有两个不同的等位基因，并且它们都与外群不同，那么我们将它们算作一个多态性和一个固定差异。虽然在少数染色体中发现的固定差异可能不是物种范围内的固定差异，但只要使用相同的标准来调用非同义和同义固定差异，这个定义就没有问题。事实上，MK 检验中的固定差异不一定来自仅两个物种的比较。在从三个物种中取样染色体的情况下，无论是否从所有物种中收集了多态性数据，我们都可以极化单个谱系的变化以用作固定差异。最好将谱系上的极化变化与从同一谱系的尖端收集的多态性数据进行比较，因为我们预计中性突变率在这种比较中最为相似。此外，可以通过比较基因重复之间的固定差异来进行 MK 检验，即使它们属于同一物种。然而，这种类型的比较有额外的注意事项，应仔细考虑（见下文）。对多态性和发散程度的期望 McDonald-Kreitman 检验明确依赖于基于分子进化中性理论的期望。上面列出的 2 × 2 表测试了多态性和发散性中非同义和同义变化的比例是否相同。测试中使用的四个值的中性期望为：

N e nonsyn nonsyn

其中 mnonsyn

是非同义位点的中性突变率，msyn

是同义位点的中性突变率，a 是 Watterson 对样本大小的校正（公式 3.6），Lnonsyn

分别表示非同义和同义位点的数量。回想一下，我们需要两个不同的中性突变率，因为对于任何特定基因的非同义和同义位点，所有中性突变的比例都会有所不同。之前我们假设 mnonsyn

等于某个 f比例，f0

总突变率，n，并且 msyn

等于总突变率，因为所有同义突变都是中性的。不管这些

分数是多少，也不管所有同义突变是否真正具有选择性

等效性，如果我们假设两种类型的位点内随时间推移的中性突变率相同，则两个 PN

的比率等于 mnonsyn

。或者，我们可以将 MK 检验视为非同义和同义突变的多态性与发散性比率的比较，解释没有变化。在这种情况下，两个 PN

等于 4Ne

/2t，至少在单个基因座内。

在任一比较中，2×2 表中的内置期望是比率应该相等。这一预期是在假设所有观察到的多态性和固定差异都是中性的前提下得出的：没有有利突变的固定，样本中也没有发现平衡的多态性或有害的多态性。如果满足这一假设（通常不满足），McDonald-Kreitman 检验将不显著。因此，该方法明确测试是否存在正向选择的固定，或者是否存在平衡和/或轻微有害的多态性。当存在强烈的净化选择时，MK 检验将不显著——在强烈的净化选择下比率的期望与没有选择时完全相同。与 dN 的比较不同，检测净化选择不是该方法的目的。如果我们将非同义和同义位点的数量作为两个附加单元格，并将表格变成 3 × 2 表格，我们将能够推断是否存在任何净化选择的特征，但这种方法通常不使用。此外，当多态性和固定非同义差异过多时（即 pN > 1 和 dN > 1），MK 检验将不显著。在这些情况下，表中的比率将相同，而且，同样，由于我们在测试中不使用位点数，我们将无法检测到这种偏离中立性的情况。解释 McDonald-Kreitman 检验

图 7.10 显示了原始

非同义

多态性

同义

图 7.10 带有正向选择证据的 McDonald-Kreitman 表，其中有

P = 0.007；Fisher 精确检验。数据来自果蝇、拟果蝇和 D. yakuba 的 Adh 基因的多态性和固定差异。（根据 McDonald 和 Kreitman 1991 年的研究。）

McDonald 和 Kreitman (1991) 对 Adh 基因的研究。非同义多态性与同义多态性的比例为 2:42，而发散的比例为 7:17。这些数据的 McDonald-Kreitman 检验得出 P = 0.007，通常被解释为显示适应性自然选择的证据。这种解释来自比较中非同义固定差异的明显过量。当然，这是一个 2×2 独立性检验，因此显著结果可能是由于任何细胞中计数过量或不足造成的。然而，正向选择的解释是生物学动机的，因为同义变体数量不会发生大的变化（而非同义变体数量不会随之发生变化）。在下一节中，我们将讨论几种由于违反了方法的假设而导致这种解释可能不正确的情况。直接选择

当氨基酸替换反复固定时，基因中非同义固定差异会显著增加，例如 dN

> 1。但请注意，MK 测试所需的固定次数比 dN 比较要低得多（事实上，对于同三种果蝇之间的 Adh，为 0.114）。无论选择有多强，MK 测试仍然无法检测到仅对一个或几个替换的正向选择，但其效力大大增加。此外，与 dN

和 MK 测试一样，我们不知道选择了哪些特定的氨基酸。

MK 检验的真正威力在于，它使用多态性数据 (PN) 作为非同义中性突变与同义中性突变比率的代理，而没有正向选择变异的贡献。虽然很明显，正向选择将如何影响固定差异的数量，但有利突变在群体中传播时作为多态性所花的时间太短，因此预计它们对 PN 的贡献不大。Smith 和 Eyre-Walker (2002) 提供了一个显著的例子，说明适应性突变对多态性和发散的贡献存在二分法：如果平均优势为 Nes = 25 的突变以中性突变率的 1% 发生，它们将控制占所有固定差异的 50%，但仅占所有多态性的 2%。因此，假设非同义和同义位点的中性突变率随时间保持不变，任何由于正向选择而产生的过量非同义固定差异都将在 2×2 表中显示为 DN 比率增加，并且 MK 检验将显著（假设采样了足够多的变异位点）。除了由于非同义固定差异过多而显著之外，当非同义多态性过多时，MK 检验也可能显著。样本中所有非同义多态性都是中性的假设经常被违反。事实上，即使是 McDonald 和 Kreitman 的原始 Adh 数据也包括众所周知的平衡氨基酸多态性（Hudson、Kreitman 和 Aguadé 1987）。虽然单个位点上的平衡选择预计会提高链接位点的多态性水平（见下一章），但这种增加应该同样影响链接的中性非同义和同义多态性；只有单一的平衡多态性才会导致 PN 比率的过度变化。一般来说，平衡选择似乎是非同义多态性显著过剩的有限解释。只有当存在强大的多等位基因平衡选择时，才会出现多种氨基酸多态性的过剩，在这些情况下，也可能存在非同义固定差异的过剩，导致不显著的氨基酸多态性显著过剩的更常见解释是这些变体略有有害。图 7.11 显示了 McDonald-Kreitman 表的一个示例，其中非同义多态性显著过量（Nachman、Boyer 和 Aquadro 1994）。如果这些变异有轻微的有害性，它们将增加非同义多态性的水平，但不会被固定。请注意，这仅适用于弱有害变异——强有害突变不会增加

非同义

多态性

同义

图 7.11 有弱有害多态性证据的 McDonald-Kreitman 表（P = 0.004；Fisher 精确检验）。数据来自 ND3 基因，多态性数据来自 Mus domesticus，并且相对于 M. spretus 有固定差异。 （nachman、

Boyer 和 Aquadro 1994 年）

频率可观，因此不太可能被采样；在这些情况下，MK 检验将不显著（Akashi 1999）。

但是，任何关于样本中存在轻微有害多态性的推断都可能受到染色体采样方式的强烈影响。尽管我之前说过，物种内染色体的非随机采样不太可能导致 MK 框架中非同义固定差异明显过剩，但非同义多态性过剩的情况并非如此。异常高数量的氨基酸多态性是由弱有害多态性引起的结论基于样本来自单一种群的假设。例如，如果从不同的亚种群中抽取单个个体，则每个亚种群中的局部适应可能会导致许多氨基酸的差异。这些差异可能会被误解为弱有害非同义多态性过多，而不是局部适应性多态性过多。虽然这种采样方案似乎不寻常，但它在亚种群内变异程度较低的物种中很常见，例如自交植物拟南芥。这些轻微有害的多态性预计在种群中出现的频率较低（见公式 7.15），这一事实无助于区分有害和局部适应性变体：因为每个地区只抽取一个种质，在这两种情况下，非同义变体在样本中的频率可能较低，但在亚种群中的频率可能较高。 MK 检验的方法扩展

非中性多态性（尤其是轻微有害的多态性）的存在可以掩盖正向选择模式。要了解为什么会出现这种情况，请考虑一个假设的位点，该位点既具有过多的固定适应性非同义替换，又具有过多的轻微有害的非同义多态性（图 7.12A）。如果只比较多态性和固定位点的数量，MK 检验就不会很显著（图 7.12B；P = 0.093）。这是因为 MK 检验假设所有多态性都是中性的：过多的有害非同义多态性会增加非同义与同义分离位点的比例，使得看起来在中性进化的情况下，大量非同义固定差异是可以预料的。

我们该如何处理令人困惑的存在轻微有害等位基因？ Templeton (1996) 建议扩展 MK 检验，将

单一多态性（两种类型）与所有其他多态性区分开来。 由此产生的 3 × 2 表（图 7.12C）试图解开轻微有害等位基因的影响

直接选择

分离位点数

非同义

同义

非同义

同义

1 2 3 4 5 6 7 8 9固定

具有衍生等位基因的个体数

多态性

单一

非同义

非同义

同义

多态性

同义

图 7.12 一个假设的基因例子，其中低频和固定非同义差异过多。 (A) 假设某个基因的非同义和同义变化的等位基因频率谱，样本染色体数为 n = 10。固定等位基因被包含在频率区间中。

(B) 该基因的标准 McDonald-Kreitman 表 (P = 0.09；Fisher 精确检验)。(C) Templeton (1996) 提出的 3 × 2 表，将单一多态性与所有其他多态性区分开来 (P = 0.07；Fisher 精确检验)。(D) Fay、Wyckoff 和 Wu (2001) 提出的 2 × 2 表，仅包括在一条以上染色体上发现的多态性 (P = 0.04；Fisher 精确检验)。

有害多态性（预计在低频下过度表达）与中性多态性（在较高频率下存在）之间的差异。虽然低频下存在中性非同义多态性，但在中性下也应该有一定比例的同义多态性（频率相同）。该表本质上是 Sawyer、Dykhuizen 和 Hartl（1987；图 7.7）提出的检验和 McDonald-Kreitman 检验（图 7.9）的组合。Akashi（1999）提出了对这种方法的进一步扩展，比较了非同义和同义等位基因（而不仅仅是单例和多例）的全频谱，并将固定差异作为频率类别之一（如图 7.12A 所示）。他表明，这种方法比标准 MK 检验和 Templeton 对该检验的扩展都更有效，尽管可能需要大量的多态性才能充分填充每个频率类别。在 3×2 表中或使用全频谱进行显著性检验可以提供更大的能力来检测过多的适应性固定差异，但当不同频率类别中非同义差异与同义差异的比率不同时，结果也可能难以解释。Fay、Wyckoff 和 Wu (2001) 进行了一项更为简单的比较，他们只是简单地删除了样本中频率 <15% 的所有多态性。他们表明，通过删除低频多态性 - 并使用常见多态性和固定差异简单地进行 2×2 检验（图 7.12D；P = 0.045） - 他们能够更轻松地检测出适应性自然选择的模式。尽管 Fay、Wyckoff 和 Wu 使用的 15% 截止值似乎有些武断，但事实证明，这个值是解决有害氨基酸分离问题的近乎最优的解决方案，尽管仍会错过许多正向选择的情况（Charlesworth 和 Eyre-Walker 2008）。在大多数中等大小的样本中，只需从 MK 检验的多态性计数中去除单个样本，就会产生类似的效果。

MK 检验的假设

在考虑进一步使用 McDonald-Kreitman 检验之前，值得回顾一下从该检验中得出推论时所做的主要假设。总体而言，我们可以相当有信心地说，使用 MK 检验得出的正向选择推论对大多数假设都是稳健的，包括染色体的采样方式、样本的近期种群历史、非中性多态性的存在、位点之间的重组，甚至所有同义突变的选择性等效性。最重要的假设可能只是在检测非同义与同义差异时没有偏差。然而，有几种情况下，检验的假设可能会被违反，从而导致拒绝零假设，特别是对于下一节中讨论的 MK 检验的几种新用途。MK 检验最不切实际的假设可能是，无论是非同义变化还是同义变化，中性突变率都会随时间保持不变。对于非同义变化，这意味着所有中性突变的比例不会随时间而变化（即选择性约束是恒定的），即使该值可能对有效种群大小 Ne 的变化非常敏感。明确地说，这个假设来自于这样一个事实：在 2×2 独立性检验中，对 PN 的期望等于 mnonsyn；如果这些速率随时间而变化，那么这个假设就不再成立。在大多数情况下，我们预计，即使约束确实在单个基因上随时间而变化，这种变化也不具有方向性，因此对 MK 检验的影响很小（Fay 和 Wu 2001；图 7.13）。

然而，在整个种群层面，有效种群规模随时间持续变化可能会导致约束水平出现全面差异。例如，如果一个感兴趣的物种的种群规模在遥远的过去（大多数固定差异积累时）平均小于样本的近期历史（大多数多态性积累时），那么我们预计固定差异的中性突变率会高于多态性，因为目前有一定比例的轻微有害突变实际上是中性的。这种直接选择发散µ非同义/µ同义多态性图 7.13 非同义突变与同义突变比率的假设变化，这些突变会固定或多态。对于单拷贝基因（红线），选择强度可能随时间而变化，箭头表示在多态性数据中发现的比率。对于新复制的基因（黑线），选择强度最初是放松的，然后在获得新功能时发生变化（用星号表示）。在这些情况下，将多态性与发散进行比较可能会产生误导。（根据 Fay 和 Wu 2001 修改。）这种情况可能导致拒绝不是由于正向选择的零假设。相反，种群规模可能在最近崩溃，增加了中性突变的比例，并使其看起来好像存在大量有害的分离多态性。最近的种群规模较小也可能导致分离多态性的数量较少，从而降低测试的整体统计能力（例如，Parsch、Zhang 和 Baines 2009）。一般而言，种群规模变化的幅度必须非常大才能对大多数正向选择推论产生影响（Eyre-Walker 2002）。为了排除最近种群规模增加是导致他们在果蝇中观察到的适应性自然选择模式的原因的可能性，Fay、Wyckoff 和 Wu (2002) 测试了增加种群规模模型的几个辅助预测。他们比较了非洲果蝇和非非洲果蝇的 PN 比率，以查看非洲以外（该物种的本土范围）的比率是否较低，结果没有发现差异。他们还询问非非洲样本中是否存在更多罕见的非同义变体（如果对有害多态性的平均选择效应增加，则可能出现这种情况），结果没有发现差异。最后，他们推断，总体人口规模的增加将影响基因组中的所有基因，导致 DN 比率全面增加；相反，他们发现基因之间存在异质性影响，一组中性进化基因和一组适应性进化基因之间存在显着差异。总之，这些结果表明，人口规模变化不是产生误导性正向选择模式的原因。有一种情况是，MK 测试（以及恒定中性突变率假设）将不断被违反，以至于在没有发生任何正向选择时会推断出正向选择。这种情况是最近的基因重复（如 Jones 和 Begun 2005；Thornton 和 Long 2005；Arguello 等人 2006 中所述）。基因重复进化的一个常见模式是，对新复制的基因的选择较为宽松（详见 Conant 和 Wolfe 2008；Hahn 2009）。如果在重复历史的早期几乎没有约束，那么许多非同义固定差异将会积累，因为它们是中性的。如果选择在最近变得更强——例如当重复发现新功能时——那么非同义多态性与同义多态性的比率将提供与发散截然不同的图景。这一结果可能导致拒绝中性假设，并解释基因重复历史中的正向选择，仅仅是因为相对于非同义多态性而言，固定的非同义差异过多。但编码序列上的正向选择并不一定发生，因为环境变化——甚至是一个有利的调节变化——可能是新功能和随之而来的选择压力变化的原因。MK 检验的第二个主要假设是样本中非同义和同义差异的平均谱系历史是相同的。在非重组区域中，这显然是正确的：只有一个（可能是未知的）谱系与所有样本相关，因此与所有位点相关。然而，在重组中，序列的不同部分具有不同的谱系，因此它们最近的共同祖先的时间可能更长或更短（因此分离位点更多或更少）。如果距离最近的共同祖先的时间与某个区域中的非同义或同义变化之间存在关联，则每种类型的分离位点数量可能会有所不同，这仅仅是因为各个区域之间的进化差异（有关进化差异的讨论，请参阅第 6 章）。明确地说，这个假设来自于这样一个事实：在 2×2 独立性检验中，对两个 PN 的期望都等于 4Ne /2t。如果样本的 MRCA 的预期时间（此处用 Ne 的值表示）在不同区域之间不同，则此假设不成立。但由于基因座内的非同义和同义变化通常相互交错，因此预计平均谱系相同。对于单个基因，即使非同义和同义变化在空间上是分开的，相邻区域的谱系中也不太可能存在较大的差异。在这两种情况下，2×2 表中需要考虑的唯一方差来源是染色体和突变之间的抽样方差。然而，越来越多的情况是，被比较的突变类别（例如，非同义和同义）在基因组中是空间分离的（见下一节）。在这些情况下，应用 2×2 独立性检验不再合适，因为相对于非重组位点，区域间分离位点预期数量的方差被夸大了。这种夸大的方差导致中性模型的错误拒绝增加（Andolfatto 2008）。尽管此类应用确实违反了 MK 检验的假设，但仍有几种解决方案。Andolfatto (2008) 建议进行重组合并模拟，以建立更准确的假设检验显著性水平。此类模拟可用于模拟基因座之间的任何重组量，因此可用于各种情况。或者，测试可以明确考虑进化和抽样方差。HKA 检验（见第 8 章）在概念上与 MK 检验相似，但旨在用于多个不相关的基因座；因此，该检验的统计框架纳入了预期的方差增加。下一章将详细阐述该检验方法，但现在只需说明，如果不同类别的位点没有相似的谱系，则必须使用考虑基因座间进化差异的检验方法。

MK 检验方法的生物学扩展

MK 检验方法的关键见解是，通过比较不同类别的位点（受选择的位点与假定不受选择的位点），我们可以在很大程度上控制这样一个事实：种群不太可能满足 Wright-Fisher 模型的所有平衡假设。换句话说，与假定无功能的多态性类别进行比较，可以为由于人口统计学或连锁选择而导致的非平衡历史提供内部控制。鉴于这一见解，我们自然会问，是否可以在其他类型的多态性之间进行这种比较，特别是如果它们可以自然地分为选定类别和非选定类别。事实上，已经对 MK 测试进行了几种不同的扩展，以测试不同类型突变中的选择。Akashi (1995) 比较了果蝇中优选和非优选同义替换的多态性和差异。在这种生物和其他生物中，经常发现编码相同氨基酸的多个同义密码子中的一个比其他密码子更受青睐，无论是在翻译准确性方面还是在翻译效率方面（参见 Plotkin 和 Kudla 2011 的评论）。在这些情况下，优选密码子被称为优选密码子，而非优选密码子被称为非优选密码子。Akashi 没有比较 MK 表中的非同义和同义变化，而是使用了优选和非优选同义变化，并能够证明许多非优选多态性具有弱有害性。由于同义变化再次在单个基因中相互交叉，因此 MK 检验的这种扩展不会出现与进化方差增加相关的问题。MK 检验的进一步扩展检查了调控序列的多态性和分歧。该应用最初由 Ludwig 和 Kreitman (1995) 以及 Jenkins、Ortori 和 Brookfield (1995) 提出，但此后已应用于越来越多的数据集（例如，Crawford、Segal 和 Barnett 1999；Kohn、Fang 和 Wu 2004；Andolfatto 2005；MacDonald 和 Long 2005；Holloway 等人 2007；Jeong 等人 2008）；类似于 dN 的测试也可用于调控区域（例如 Hahn 等人，2004 年）。这些测试比较非编码区域中的结合位点以及散布在不结合转录因子的位点或附近区域的同义位点。在比较多态性和不相连（或仅松散相连）区域的固定差异的情况下，必须再次谨慎考虑基因座之间的不同历史（Andolfatto 2008）。与解释进化方差几乎同样重要的是编码区非同义变化和非编码区调控变化的进化之间的几个重要差异（Hahn 2007）。首先，因为功能性调控序列通常仅在单个物种中表征，所以在物种间比较中将功能同源性分配给一组同源核苷酸并不总是必要的。实验验证的结合位点存在于一个研究充分的焦点物种中，而其他物种中可能不存在；相反，在焦点物种中没有已知功能的序列实际上可能是用于比较的其他物种的结合位点。这两种类型的错误都会导致核苷酸分类错误、选择强度和方向估计错误以及违反中性突变率恒定的假设。第二个警告是结合位点的遗传密码目前未知。我们几乎没有关于结合位点变化对结合亲和力的影响的信息。尽管将结合位点内的任何变化归类为选择在很大程度上是一种假设，但它可能与将蛋白质中的任何氨基酸变化视为功能相关一样好。将 MK 测试应用于调控序列的第三个主要警告涉及正向选择的作用方式。虽然蛋白质中氨基酸的重复替换似乎是正向选择的良好证据，但很难想象这在调控区域中究竟如何起作用。这种怀疑源于调控序列的一些重要特征：结合位点通常不限于特定位置，结合位点经常通过点突变出现，结合位点的多次变化通常会导致结合亲和力的完全丧失（Wray 等人 2003 年综述）。这些原因都不排除自然选择以这种方式发挥作用；相反，它们只是表明，由于定向选择而可检测到重复替换的情况很少见。考虑到所有这些警告，实际上有一些通过使用非编码区域上的 MK 测试检测到正向选择的例子。Crawford、Segal 和 Barnett（1999）研究了两种鳉鱼亚种 Fundulus heteroclitus 之间 Ldh-B 基因座调控区域的变异。他们发现，与散布位点的替换相比，负责转录因子结合的核苷酸存在过多的固定差异，并且物种间 Ldh-B 调节序列驱动的表达存在差异。这些结果与导致调节突变固定的重复正向选择一致。

总结整个基因组的选择

鉴于 McDonald-Kreitman 测试对多个基因的结果，我们希望能够比较基因之间的结果或总结单个基因组的结果。从 2×2 列联表测试中获取的 P 值并不适用于基因之间的比较。这种情况的一个原因是，具有正向选择证据的基因和具有分离有害多态性证据的基因可以具有相同的 P 值。第二个原因是 P 值高度依赖于每个细胞中的计数 — 因此，无论选择强度如何，更长或进化更快的基因将具有更极端的 P 值。简单地将基因间的列联表相加以增加计数也不是一个好的解决方案，因为它会导致误导性的结果（Shapiro 等人，2007 年；Stoletzki 和 Eyre-Walker，2011 年）。为了提供一个简单的汇总统计数据，为 MK 检验的 2 × 2 列联表提供可解释的值，Rand 和 Kann（1996 年）提出了中立指数 (NI)：在中立的情况下，中立指数的预期值为 1，因为 的比率应该相等。大于 1 的值表示非同义多态性过量，小于 1 的值表示非同义固定差异过量。中性指数或其变体（见下文）已被广泛用于比较生物体和基因之间的分子进化模式。如原始论文所述，“该指数旨在提供偏离中性的方向和程度的定性指标”（Rand and Kann 1996，第 737 页）。仍然必须计算 P 值才能推断出与 NI = 1 的统计显着偏差，但该指数允许快速比较基因之间的进化力量。

除了一般的偏差问题和作为比率的比率的统计数据中的方差增加问题之外，在解释 NI 时还存在几个问题。当 DN 的值为 0 时会出现一个问题。在这些情况下，NI 未定义，无法为此类基因计算任何值。此问题的常见解决方案是在 2×2 表中的每个单元格中添加一个“伪计数”1，这确保 NI 被定义，尽管假设检验应该在没有伪计数的情况下进行。第二个问题是因为 NI 值不围绕 1 对称，因为 NI = 2 和 NI = 0.5 代表与中性预期的等效偏差，但幅度却有很大差异（Stoletzki 和 Eyre-Walker 2011）。解决这个问题的方法是取 NI 值的对数（或负对数），这会产生对称偏差（例如，Tachida 2000；Presgraves 2005）。负对数还具有令人愉悦的美学副作用，即对显示具有正选择系数的替换证据的基因赋予正值，对显示具有负选择系数的多态性证据的基因赋予负值（Li、Costello 等人，2008 年）。

图 7.14 显示了酵母、苍蝇和人类中数千个基因的负对数（NI）值，并立即证明了苍蝇中显示正选择证据的大量基因。

一个相关的统计数据是 a，旨在表示由正选择固定的非同义替换的比例（Smith 和 Eyre-Walker，2002 年）。

a 的值由以下公式给出：

H. sapiens

D. simulans

−log10[中性指数]

S. cerevisiae

图 7.14 人类、苍蝇和酵母的中性指数。每个点代表特定基因的中性指数 (nI) 以及与该基因的 McDonald-Kreitman 检验相关的 P 值。为了绘图目的，nI 定义为 + 1)，并显示这些值的负 log10。每个图中的下部水平线代表 P = 0.05；上部水平线代表 Bonferroni 校正的 P = 0.05/n，其中 n 是进行的测试次数。（来自 li, Costello, et al. 2008.）这与 1 − NI 相同（Smith and Eyre-Walker 2002）。事实上，a 不是一个比例，因为它可以取负值。最好将 a 视为一个汇总统计数据，因为不清楚当 a 为负值时（即当存在过多的非同义多态性时）如何解释 a。然而，正值确实代表了由于适应性自然选择而导致的氨基酸替代过量。与为基因组中的每个基因计算单独的值相比，a 最常用于汇总基因间的数据。这可以通过以无偏的方式对基因间的 DN 计数进行平均来实现。不幸的是，概述此汇总方法的原始论文描述了违反辛普森悖论（即组合两个数据集可能导致两个数据集中均不存在的结果）或詹森不等式（即比率的平均值不等于平均值​​的比率）的程序。虽然存在用于计算 a 的最大似然法（Bierne 和 Eyre-Walker 2004；Welch 2006），但最近的研究提供了一种简单的无偏方法来计算大量基因的平均 a（Stoletzki 和 Eyre-Walker 2011）：

其中下标表示总共 m 个基因集合中的第 i 个基因。

计算果蝇和果蝇基因组中非同义差异的 a 发现 a = 0.54，这是正向选择替换的很大一部分（54%；Begun 等人 2007）。

最后，通过将等位基因频率谱数据与固定差异数据相结合，有几种方法能够更好地估计由正向选择固定的非同义替换的比例（例如，Boyko 等人，2008 年；Eyre-Walker 和 Keightley，2009 年）。这些方法使用 Templeton（1996 年）和 Fay、Wyckoff 和 Wu（2001 年）提出的推理，首先考虑轻微有害的非同义多态性的存在，然后再估计这一比例。虽然方法略有不同——Eyre-Walker 和 Keightley 使用折叠等位基因频率谱，而 Boyko 等人使用展开谱——两者都试图首先通过与同义多态性谱进行比较来估计非同义多态性的比例和频率。通过将人口统计学模型拟合到同义数据集并假设所有非同义多态性都是中性的或有害的，这些方法可以估计由于适应性自然选择而导致的非同义固定差异的比例。这种方法和类似的方法也可用于计算联合参数 g = 2Ne s 的值（也容易混淆地表示为 a），该值表示选择对非同义变体的平均影响（Keightley 和 Eyre-Walker 2010）。Si模拟结果表明，如果收集大量基因的数据，即使样本量很小，这些方法也可以提供平均值的准确估计；准确估计单个基因座的参数需要更大的样本量（Keightley 和 Eyre-Walker 2010）。随着测序变得越来越便宜，甚至对有害和有利突变频率的单基因座研究也将成为可能。

选择 8

在上一章中，我们讨论了自然选择对影响生物体适应性的个体突变的影响。如果单一类型的突变足够多（例如有利的非同义突变），那么我们通常可以对主要的选择模式做出有力的推断。然而，前面描述的方法都不能用于检测对单个有利突变的选择，无论它有多大的影响，也不能用于检测对单个平衡多态性的选择。为了识别这种选择性变异，我们必须考虑它们对连锁中性变异的影响：多态性本身对适应性没有影响，但会受到附近发生的选择的影响。虽然使用连锁中性变异的方法通常仍然局限于检测强烈的近期选择，但它们或许是了解基因组近期选择历史的最佳窗口。我研究了连锁选择对变异的三个不同方面的影响，并针对每个方面讨论了旨在检测​​选择对这种变异特定方面的作用的测试。

使用多态性的数量检测选择

选择对连锁中性多样性的水平有相对直接的影响：它可以提高或降低它们。了解中性变异水平上升或下降的一般条件将有助于我们直观地了解选择测试如何工作以及可以检测到哪些选择模式。

正向选择对连锁中性水平的影响

... (A) 有利突变（以星号表示）出现在中性谱系中。(B) 有利等位基因的频率增加，但尚未在种群中固定。(C) 有利等位基因已固定，所有染色体在其出现后不久就合并。为清晰起见，显示了从样本根部延伸的分支。

所有完全连锁的变异都被扫除，唯一剩余的变异发生在过渡到固定期间。短语搭便车（或 Maynard Smith 和 Haigh 的原始拼写为“搭便车”）强调了这样一个事实：恰好位于有利突变发生的幸运染色体上的中性多态性的频率也会增加（图 8.1）。

这些搭便车者也将在种群中固定，除非距离选定基因座更远的重组将它们与有利等位基因分开。

为了更定量地了解选择性清除如何减少基因位点的多态性，请考虑突变在群体中固定所需的时间。这段时间决定了会有多少变异：在最极端的例子中，所有染色体都来自一个包含上一代发生的高度有利突变的单一谱系。在这样的群体中几乎没有变异，因为所有多态性都必须在一代中发生。突变固定所需的时间越长，相关多态性积累的时间就越长。例如，在固定的条件下，中性突变大约需要 4Ne 代才能固定（Kimura 和 Ohta 1969），尽管与这一预期相关的方差很大（Kimura 1970）。中性突变的固定时间与中性谱系中所有谱系融合到其最近的共同祖先所需的时间之间存在明显的联系，这也接近 4Ne（公式 6.3）。因此，中性基因座的变异量（以及谱系的形状和高度）与中性变体固定所需的时间密切相关。现在想象一下，一旦在二倍体种群。这种突变的修复时间是

（Nei 1973，第 383 页）：