厦门免费医学论文发表-一种从时间序列遗传数据中推断一般二倍体选择的新型期望最大化方法
亚当·法恩，马蒂亚斯·施泰因吕肯
抽象
检测和量化选择强度是群体遗传学的一个主要目标。由于选择作用于多代，因此已经开发了许多方法，使用在多个时间点采样的遗传数据来检测和量化选择。这种时间序列遗传数据通常使用隐马尔可夫模型进行分析，但在大多数情况下，是在加性选择的假设下进行的。然而，存在许多表现出非加性机制的遗传变异的例子，这使得开发能够在更一般的情况下表征选择的方法变得至关重要。在这里，我们将先前介绍的用于推断加性选择系数的期望最大化算法扩展到一般二倍体选择的情况，其中杂合子和纯合子适应度是独立参数化的。我们还引入了一个框架来识别从给定数据中选择二倍体的定制模式，一个解释可变种群大小的启发式方法，以及一个跨链接基因座聚合数据以提高功效和稳健性的程序。通过广泛的仿真研究，我们发现我们的方法在各种场景中准确有效地估计了不同二倍体选择模式的选择系数;但是，除非选择非常强大，否则对选择模式进行分类的能力很低。我们将我们的方法应用于过去 4,450 年中来自英国的古代 DNA 样本，并在六个基因组区域发现了选择的证据，包括特征明确的 LCT 基因座。我们的工作是第一次全基因组扫描，表征了一般二倍体选择的信号。
作者总结
自然选择增加了有益遗传变异从父母传递给后代的可能性，从而构成了遗传适应新环境的基础。在多个时间点采样的基因组数据，例如从古代遗骸（古代 DNA）中提取的遗传物质或来自进化和重测序实验的数据，可以比单独的当代样本更精确地识别受选择压力影响的遗传变异。然而，大多数在选择下识别遗传变异的方法都集中在加性选择上，其中杂合子的适应性恰好在纯合子之间处于中间位置。利用多个时间点的遗传数据，我们开发了一种方法来检测加性和非加性选择，并推断最可能的显性机制。我们将我们的方法应用于距今不到 4,450 年的英国古代 DNA 数据集，并确定了六个具有近期选择信号的区域，其中包括一个位于 TFR2 基因座的区域，该位置以前未作为选择目标报道。我们的工作能够更准确地量化非加性选择动力学，并可用于测试更复杂的选择模型。
数字
Fig 15Table 2Fig 1Fig 2Fig 3Fig 4Fig 5Fig 6Fig 7Fig 8Fig 9Fig 10Fig 11Fig 12Fig 13Table 1Fig 14Fig 15Table 2Fig 1Fig 2Fig 3
引文： Fine AG， Steinrücken M （2025） 一种从时间序列遗传数据中推断一般二倍体选择的新型期望最大化方法。PLoS 基因 21（7）： e1011769 号。 https://doi.org/10.1371/journal.pgen.1011769
编辑 器： Parul Johri，亚利桑那州立大学坦佩校区：美国亚利桑那州立大学
收到： 2024 年 12 月 13 日;接受： 2025 年 6 月 11 日;发表： 7月 22， 2025
版权所有： © 2025 Fine， Steinrücken。这是一篇根据知识共享署名许可条款分发的开放获取文章，该许可允许在任何媒体上不受限制地使用、分发和复制，前提是注明原作者和来源。
数据可用性： 生成结果的所有代码都可以在 https://github.com/steinrue/EMSel 上获得。用于分析的古代 DNA 数据是从 https://doi.org/10.7910/DVN/FFIDCW 的 Allen Ancient DNA Resource （AADR）（7.0 版）下载的。
资金： AGF 得到了教育部、国家需求领域研究生援助 （GAANN） 和 Grant #P200A210054 的支持。AGF 和 MS 得到了美国国立卫生研究院国家普通医学研究所 （NIGMS） 的支持，R01GM146051 授予 MS。资助者在研究设计、数据收集和分析、发表决定或手稿准备方面没有任何作用。
利益争夺： 作者已声明不存在相互竞争的利益。
介绍
赋予生物体优于其同类适应性优势的遗传变异往往会随着时间的推移在种群中增加频率，直到最终固定，如果它没有因遗传漂变而丢失的话。这种随机的选择过程最终构成了适应的基础。因此，检测选择的证据并量化其强度是进化生物学中的一个基本问题，其应用范围从寻找对早期原始人进化至关重要的突变 [1] 到预测肿瘤生长 [2]。因此，在群体遗传学中，已经开发了许多方法来检测当代群体基因组数据中过去选择性事件的特征[3–5]。
然而，由于选择作用于多代，因此在整个过程中的多个时间点观察到的遗传数据比单独使用当今样本更准确地量化选择性过程。最近的技术进步使研究人员能够在全基因组范围内收集此类时间序列遗传数据。时间序列遗传数据的一个主要来源是古代 DNA （aDNA），即从人类或其他物种的已故个体中提取的遗传物质 [6]。新一代测序使得从大量古代样本中收集遗传数据成为可能，特别是通过开发杂交富集等技术[7]。时间遗传数据的另一个主要来源是实验进化研究[8]。当代实验进化研究在进化和重测序 （E&R） 实验中对多个生物重复使用下一代测序技术，以获得对整个基因组中许多位点的时间等位基因频率变化的高质量估计 [9]。这些数据集为检测和表征塑造基因组变异的适应性过程提供了前所未有的机会 [10， 11]。
观察真实的潜在群体等位基因频率轨迹随时间变化将有助于高度准确地表征潜在的选择性过程。然而，在实践中获得的数据中，遗传变异通常只针对在有限数量的时间点采样的一组个体进行评估。因此，量化选择强度涉及对选择、遗传漂变和其他群体遗传过程对群体等位基因频率的未观察到轨迹的作用进行建模，并将采样数据视为对该潜在轨迹的不精确观察。
一种常用的时间序列数据分析框架是隐马尔可夫模型 （HMM） [12]。在这些 HMM 中，潜在的群体等位基因频率根据马尔可夫过程、Wright-Fisher 模型演变，并且样本在给定潜在群体等位基因频率的情况下被建模为二项式观察。该 HMM 框架已被用于估计各种参数：加性选择系数 s [12]、有益突变出现的时间 [13]、有效种群大小 Ne [14] 或测序错误率 [15]。在这个 HMM 框架中，[16] 引入了一种期望最大化 （EM） 方法，可用于估计加性选择系数以及亚群之间的迁移率。有关基于 HMM 的估计选择系数的方法的综述以及 aDNA 和 E&R 分析中方法之间的比较，请参见 [17] 和 [18]。
迄今为止，文献中介绍的上述 HMM 方法的大多数实施都旨在仅检测加性选择，但存在许多表现出非加性机制的遗传变异的例子。在人类中，非加性选择靶标的范围从镰状细胞等位基因的一个拷贝赋予杂合子优势的经典案例 [19， 20] 到最近对英国生物样本库数据分析中普遍占主导地位的证据 [21]。此外，对复杂性状的稳定选择被认为在人类中很普遍[22]，表现为影响性状的基因座的选择性动力学不占主导地位[23–26]。
在这里，我们扩展了 [16] 中的 EM 方法，以估计一般二倍体模型下的选择系数，即当纯合子和杂合子基因型的适应度值独立参数化时。与基于网格搜索的方法相比，使用 EM 方法估计选择参数迭代最大化可能性，可以更好地扩展到两个以上的参数 [27]。例如，使用 EM 算法使我们能够同时估计二倍体选择系数和表征等位基因初始频率的参数，而使用网格搜索估计这些参数将具有挑战性。
我们进一步开发了一个新的框架，用于确定给定时间数据集的最佳选择模式。虽然文献中已经介绍了其他估计一般二倍体选择系数的方法[27–32]，但它们并没有明确解决区分不同选择模式的统计问题。此外，这些方法都没有应用于人类群体的全基因组数据。据我们所知，我们的分析是第一个从全基因组规模的古代 DNA 数据中表征人类近期一般二倍体选择的分析。
本文的其余部分组织如下。在方法中，我们概述了用于有效估计一般二倍体选择系数的迭代 EM 算法，以及用于推断最可能的二倍体选择模式的统计程序。在仿真研究中，我们将算法和推理框架应用于各种仿真场景，以评估其准确性。我们发现，我们的方法通常具有很好的检测选择和估计其强度的能力，但是，对选择模式进行分类的能力是有限的。此外，在推断有效种群规模中，我们介绍了一种从时间数据估计恒定种群大小的程序。在来自英国的古代 DNA 数据集中，我们然后使用过去 4,450 年居住在英国的个体的公开可用的古代 DNA 数据对人类基因组中最近的二倍体选择进行全基因组扫描 [33，版本 54.1]，引入了一种在链接基因座上聚合 p 值的程序，并讨论了显示最近选择信号的六个基因组区域。在驯养马的毛色基因座 ASIP 中，我们还将我们的方法应用于涉及马毛着色的基因座 [34]，以证明我们的方法在探索非加性场景时的实用性。最后，我们在 讨论 中讨论了未来的方向。我们的方法 EMSel（用于检测选择的 EM 算法）和生成本手稿中图形的脚本可在 https://github.com/steinrue/EMSel 在线获得。
方法
参数化一般二倍体选择
在本文中，我们考虑了在恒定大小 N 的二倍体种群中给定的双等位基因座起作用的选择e.该基因座的群体等位基因频率的动态可以用离散 Wright-Fisher 模型来描述，其中我们用 A 和 a 表示基因座上的两个等位基因，用 和 1-Y 表示t分别是 Generation 中 A 和 a 等位基因的随机种群水平频率。假设具有基因型 AA、Aa 和 aa 的个体的相对适应度为 1 + s2、1 + 秒1和 1 分别。给定这些适应度值，如果当前一代 A 等位基因的频率为 Yt= p，则下一代 Y 中的等位基因频率T + 1由乘以 2N 的二项分布随机变量给出e小 S 的平局和成功概率1和 s2.这可以通过正态分布进一步近似，其中 ，通常称为 Wright-Fisher 扩散 [35， Ch. 5.3]。
我们使用术语一般二倍体选择来表示任意 s 的情况1>−1 和 s2>−1，即纯合子和杂合子的适应度值是独立参数化的。许多定制的选择模式对应于约束二倍体选择，其中 s 的可能值1和 s2受到限制。我们考虑以下选择模式：
添加剂： .相对适应度与 A 等位基因的拷贝数成正比。也称为单倍体或基因选择。
主导： .任意数量的 A 等位基因拷贝都赋予了完整的适应效应。
隐性： s1= 0 的只有 A 等位基因纯合的个体才具有适应效应。
优势： .杂合子个体的适应性最高。
下 优势： 。杂合子个体的适应度最低。
加法、显性、隐性选择是单参数模式，分别具有 、 和 。over- 和 underdominant 是二维子空间，但我们经常将版本用作描述完全 over- dominant 或 underdominance 的单参数选择模式。我们也将后两者的组合称为杂合子差异。虽然完全过度和欠显性的单参数模式更便于模拟，但在推断参数时将它们组合成一维模式杂合子差异可以避免统计伪影。在图 1 中，我们显示了作为一般二倍体选择的完整二维参数空间的子空间的不同选择模式的图表。
thumbnail下载：
PPT的PowerPoint 幻灯片
PNG放大图片
国际电影节原始图像
图 1. 一般二倍体选择的二维空间。
我们考虑的五种定制选择模式是：加法、显性、隐性、过度优势和欠显性。颜色表示相应模式的子空间。绘制了在每种模式下模拟的 50 个重复及其平均值，以说明每种模式下轨迹的特征形状。
https://doi.org/10.1371/journal.pgen.1011769.g001
用于推断二倍体选择系数的隐马尔可夫模型
EM-HMM 算法的推导。
为了推导出我们从采样的时间序列遗传数据中推断一般二倍体选择系数的新方法，我们扩展了 [16] 中开发的方法用于加性选择。如果已知 T 代的焦点等位基因的确切频率轨迹，则可以使用 Wright-Fisher 过程的正常近似值来定义该轨迹的可能性：
(1)
在这里，我们表示感兴趣的参数，选择系数 s1和 s2，上标 （C） 表示每代群体等位基因频率采用连续范围 [0,1] 的值的模型。在不失去通用性的情况下，焦点等位基因是 A 等位基因。最大化此似然会生成选择系数的最大似然估计器 （MLE）。在加法选择的情况下，这个估计量已在 [36] 中介绍。在 S1 文本的 S.1.1 和 S.1.2 节中，我们重新审视了 Watterson 的加性 MLE 公式，并将结果扩展到一般二倍体选择，从而产生了估计量
跟。
然而，正如 引言 中所讨论的，数据通常由在某些时间点从人群中采样的给定数量的个体组成，因此群体等位基因频率的完整轨迹通常是未知的。为了有效地整合这种不确定性，[12] 引入了一个 HMM 框架，其中给定时间的群体等位基因频率是隐藏状态，根据 Wright-Fisher 模型进化，采样的基因型是观察结果。具体来说，假设我们有时对总体进行了抽样。在每个时间点，数据由此时采样的单倍型数和观察到的局灶等位基因数组成。为了方便起见，我们设置了 nt： = 0 和 at： = 0 时未观察到数据，我们表示与在 t 处用 O 对数据进行采样相关的随机变量t.此外，我们将群体等位基因频率离散为 M 个隐藏状态，以实现高效计算：、 和 。然后，第 t 代中的隐藏状态是 t 处的离散化群体等位基因频率，我们用 表示。
要应用标准 HMM 框架 [37,38 Ch. 13.2]，我们必须定义初始概率、转换概率和发射概率。现在，我们假设给出了初始概率。这些可以是固定的或估计的，我们在 估计初始等位基因频率分布 中提供了有关估计程序的详细信息。对于从隐藏状态 g 的转换概率我到隐藏状态 Gj中，我们使用 Wright-Fisher 过程的法态近似值并定义
其中，用于一般二倍体选择 和 表示具有均值和方差 的正态分布的密度。在这里，我们还设置了 和 ，它将所有概率质量分配给区间 [0,1] 之外的转换到相应的边界点。最后，观察到tT 代中的局灶等位基因由二项分布给出
与 nt平局和成功概率 Ft，相应的群体等位基因频率。图 2A 描绘了该 HMM 的示意图。
thumbnail下载：
PPT的PowerPoint 幻灯片
PNG放大图片
国际电影节原始图像
图 2. HMM 来推断选择。
A） HMM 示意图。每个“红绿灯”代表给定时间的一个单倍体样本，具有一定数量的焦点等位基因。还绘制了通过隐藏状态空间的三种可能的群体等位基因频率轨迹。在计算算法的 E 步中的预期值时，总体频率更接近样本中频率的轨迹被赋予更多权重。B） 在加法选择下模拟的给定重复的每种模式下 EM-HMM 优化的对数似然表面和路径。
https://doi.org/10.1371/journal.pgen.1011769.g002
然后，可以使用前向-后向算法 [37,38 Ch. 13.2] 来获得 t 代中隐藏状态的后验概率
以观测数据为条件并给出选择系数 S1和 s2和隐藏状态的关节后方
基于这些后验分布，后验期望
和
可以分别针对任意函数以及边际频率和联合边际频率进行计算。这里，上标 （D） 表示使用具有离散化等位基因频率的模型。
为了找到二倍体选择系数的最大似然估计 （MLE），在此 HMM 下，我们使用 EM 算法 [37,38，Ch. 13.2]，类似于加法情况下的 [16]。我们将此算法称为 EM-HMM 算法。该算法从 initial parameter estimate 开始。在迭代 k 中，算法然后计算 and（E 步），并通过使用以下方式最大化条件对数似然来更新参数估计值
(2)
（M 步长），直到估计值收敛。请注意，在二项式模型下，发射概率与 s 无关1和 s2，并且不需要在此更新步骤中明确考虑。与从方程 1 推导出给定轨迹的 MLE 类似，方程 2 得到
其中 ，我们用 表示离散化的条件期望
为简洁起见。有关派生的详细信息，请参见 S1 Text 中的 Sect S.2.2。
请注意，这种方法不是一种精确的 EM 算法，因为它将离散化模型下计算的条件期望与连续模型中的可能性最大化相结合，类似于 [16]。这种混合方法的原因是，虽然在 [36] 中使用 MLE 的二倍体泛化可以在连续设置中最大化条件对数似然，但不能轻易计算后验期望。另一方面，在离散化模型中，可以计算后验期望，但条件对数似然不能通过解析最大化。然而，混合方法在计算上是可处理的，并产生高度准确的估计。
针对自定义选择模式的约束优化。
除了估计无约束的二倍体选择系数 s1和 s2，我们还想估计单参数选择模式加法、隐性、显性和杂合子差异中的选择系数。为此，我们在 s1和 s2使用拉格朗日乘子框架 [39， Ch. 7.5]。用 表示似然，优化问题可以表述为在方程 2 中最大化任意函数的条件对数似然，这可以通过引入拉格朗日并求解来解决。
上述所有感兴趣的单参数选择模式都可以表示为线性约束，即 for suitable 。因此，在拉格朗日乘子形式中，我们可以设置并求解约束优化问题，以获得相应选择模式下的 MLE。在 S1 Text 的 S.2.3 节中，我们为 s 的更新规则推导出显式表达式1和 s2对于任意 a、b。例如，图 2B 显示了加法选择下模拟的数据集上加法、隐性、优势和无约束的一般二倍体选择 EM-HMM 算法的迭代。
初始等位基因频率分布的估计。
初始等位基因频率的离散分布可以固定用于分析，也可以估计。在估计它时，我们将 beta 分布拟合到初始频率，因为它是一个灵活的分布，只有两个参数，用 表示，这避免了潜在的过度参数化。如果我们假设初始分布不依赖于选择系数，那么方程 2 中给出的 EM 更新规则的离散化版本变为
(3)
对于参数 .这里， 表示离散化 beta 分布的可能性，条件期望现在由 参数化，并且是第 m 个离散化区间内参数化的 beta 分布的积分。在 EM-HMM 算法的每次迭代 k 中，我们以数值方式求解方程 3 以更新并选择系数。由于选择系数的 EM 更新已经需要计算，因此更新初始条件的额外步骤以最小的计算成本出现。我们观察到，估计初始分布确实会影响选择系数估计的准确性，参见 S1 文本中的图 J，同时提供了更大的灵活性。
其他实现详细信息。
使用我们的 EM-HMM 算法分析数据需要选择隐藏状态或离散化区间 M 的数量，以及离散化点的放置位置。离散化点的常见选择是将它们等距隔开，即 。然而，我们观察到，当选择很强时，这种离散化表现不佳，详见 S1 Text 中的 S.5 节。因此，我们决定将 Chebychev 节点用于 。这些节点通常用于数值积分，因为它们减轻了 Runge 现象 [40， Ch. 4.5]。直观地说，当使用等距间距时，M 步中使用的过渡的正态分布方差变得小于边界附近离散化区间的大小，因此，离散化点处的密度不能很好地近似区间中的概率质量。使用 Chebychev 节点可以有效地增加边界附近离散化点的密度，并降低内部离散化点的密度，从而缓解了这个问题。
我们在选择系数估计的稳健性中使用隐藏状态数量 M 的不同选择分析了模拟数据。我们观察到，当使用少于 250 个隐藏状态时，推理表现不佳，但对于较高的值是稳定的。为了平衡 EM-HMM 算法的准确性和计算效率，我们在所有分析中选择 M = 500 个隐藏状态。除非另有说明，否则我们使用 初始化 EM-HMM 算法的选择系数，并将初始分布的起始参数设置为 ，这与均匀分布相对应。对于所有分析，我们使用收敛标准，即迭代 k 和迭代之间的对数似然差必须小于 10−3.此外，由于 EM 算法可能需要多次迭代才能满足此收敛标准，因此我们使用 SQUAREM 程序测试了加速 EM [41]。我们发现，虽然这种方法减少了所需的迭代总数，但由于每次迭代的额外计算成本，算法的总运行时间增加了。由于我们没有观察到估计的准确性明显提高，因此我们继续采用常规的 EM 方法。
我们注意到，特别是对于中性下的模拟，许多重复的 EM-HMM 将返回为 MLE，但实际 MLE 会略微偏离 0，对数似然增加 <10−3.尽管如此，这将导致许多重复报告的对数似然比正好为 1，因此，1 附近的 p 值没有得到很好的校准。为了避免对 p 值分布造成潜在不必要的扭曲，我们要求 EM-HMM 在停止之前至少执行 5 次迭代，如果进一步迭代不增加对数似然，则取前 5 次迭代中获得的最大对数似然。在实践中，这对筛选下模拟的重复影响最小，因为在大多数情况下，EM-HMM 需要 5 个以上的步骤才能收敛。
区分选择模式
单个分项选择的 P 值。
在描述我们处理具有多个选择模式的选择模式的完整问题的方法之前，我们首先概述了在单参数选择替代模式的情况下拒绝中立性任务的解决方案。对于给定的仿行和单参数选择模式，我们使用 EM-HMM 算法来计算 MLEs 以及对数似然 lls对于这些参数。此外，我们计算了 MLE 和对数似然 ll0在中立性 （S = 0） 下。将初始分布的参数视为令人讨厌的参数，然后我们使用似然比统计量执行标准似然比检验。当样本量趋于无穷大时，Wilks 定理 [42] 意味着在原假设 s = 0 下，D 应该以一个自由度分布，我们用 表示。然而，每个重复只由一组时间样本组成，因此我们不是在 Wilks 定理的渐近状态下作，并且对 D 的分布没有理论保证。尽管如此，在单一备择 p 值的验证中，我们表明，如果我们使用 计算 p 值，假设存在分布，则得到的 p 值是经过良好校准的，因此渐近公式提供了很好的近似值。因此，在给定的单参数替代方案下，以这种方式计算的 P 值可用于接受或拒绝给定数据集的中立性零假设。
多个备选方案。
给定数据集背后的选择模式可能不是先验已知的，因此我们设计了以下策略，旨在从给定数据中推断多个备选方案中的选择模式。我们发现，在完整的二维参数空间中对 MLE 的无约束估计比定制约束模式下的推理具有更高的方差。因此，我们推断选择模式的方法完全依赖于受约束的选择模式。由于检验统计量 D 在所有选择模式的单一替代检验的情况下都经过了很好的校准，我们提出了以下程序来确定显着性并对显着重复进行分类：
对于目标数据集中的所有重复，获取所有受约束的选择模式的 MLE 参数：加法、显性、隐性和杂合子差异。
对于每个仿行，计算检验统计量，其中 ll
0是中立性下的可能性，并且 llm是众态 m 的最大似然。
对于每个仿行，将其 p 值计算为 ，并在所需的显著性水平上拒绝中立性。
如果中性被拒绝，则根据
杂合子差异模式是最可能的模式对选择模式进行分类，如果 s
1>0 和 as underdominant if s1<0.
我们发现，使用自由度为 1 的分布根据统计计算的 p 值不如单个备择检验中的 p 值校准得好，并且可能导致分布尾部出现略微反保守的 p 值（请参阅单个备择 p 值的验证）。这并不意外，因为来自经过良好校准的单一备选 D 统计量的一些 p 值被来自其他选择模式的更高似然估计值所取代。但是，使用分布会导致更差的拟合。
此外，我们还探索了将参数化分布拟合到统计量的参数化 bootstrap 模拟。在 S1 文本的图 N 中，我们将使用 和 分布计算的 p 值与使用参数化分布（例如广义伽马分布）计算的 p 值进行了比较，这些值适合于在与我们在来自英国的古代 DNA 数据集中的数据分析相关的场景中进行 10,000 次重复的 bootstrap 模拟。虽然 bootstrap 程序在这种情况下实现了更好的校准，但该程序需要针对给定场景进行广泛的模拟以确定必要的参数。因此，我们建议使用分布，因为它是一个快速而灵活的过程，即使在分布的尾部也能产生相当良好的校准 p 值，并且部分是由理论驱动的。在实践中，我们建议在感兴趣的特定场景中模拟数据，以确认 p 值校准良好。
结果
仿真研究
为了评估我们的 EM-HMM 算法在估计选择系数并推断各种选择制度中的正确选择模式的性能，我们在多种参数组合下模拟了数据集，并展示了推理的准确性。
仿真参数。
我们在离散 Wright-Fisher 模型下模拟了给定世代数的群体等位基因频率轨迹，并在给定时间从二项分布中采样了一定数量的单倍型，成功概率由相应的群体等位基因频率给出。具体来说，我们在中性 （） 下模拟数据集，每个单参数选择模式（加法、隐性、显性），其中我们设置 s = s2，并在 s = s 时完成上占或下占1和 s2= 0 的我们使用选择系数和模拟轨迹的长度（以世代为单位）进行模拟。我们考虑了两种不同类型的初始条件：1） 每个重复的群体等位基因频率初始化为固定值，以及 2） 每个重复的初始频率是从概率成正比的集合中提取的。后者是泊松随机场模型[43]下中性分离位点的稳态分布，因此对应于从常存变化中选择。根据选择模式，我们将模拟数据条件为焦点等位基因 A 在样本中没有丢失或固定。具体来说，对于加性、显性和隐性选择，我们的条件是 A 在第一代没有固定，在最后一代也没有丢失;对于overdominance，模拟以上一代的A分离为条件;对于劣势，我们以第一代的 A 隔离为条件。
对于所有模拟，种群大小都设置为 Ne= 10,000 个。采样方案由 K = 11 个等距时间点组成，其中第一个和最后一个点分别与模拟轨迹的起点和终点对齐。每次，我们都会对单倍体进行采样。对于每个参数组合，我们在中性下模拟了 10,000 次重复，在每种选择模式和选择强度下模拟了 1,000 次重复。对于每个重复，我们使用相同的 N 进行推理e这用于四种定制选择模式下的模拟：加法、隐性、显性和杂合子差异。
在估计选择系数、单一替代 p 值的验证和选择模式推断的验证中，我们提供了 251 代和初始等位基因频率固定为 p = 0.25 的模拟结果。世代数和初始等位基因频率分布的其他组合的结果可以在 S1 文本的 S.3 部分中找到。一般来说，结果是相似的，但极端情况除外：对于初始频率较低的隐性或欠显性，选定的等位基因通常会丢失，从而导致估计不准确。
估计选择系数。
在图 3A 中，我们展示了在加法、隐性、显性、显性和欠显性模式下模拟的重复中估计的选择系数的分布。选择系数的估计通常是准确且无偏的。对于中强选择，超优势在估计值中的方差最大。可能的解释是，优势重复很快达到平稳性 （p = 0.5），导致信息样本减少，估计值方差增加。在 S1 文本的图 I 中，我们绘制了使用无约束 EM 为所有选择模式估计的选择系数，即估计两个 s1和 s2.加法、隐性、显性和超显性选择由无约束 EM 很好地估计，而欠显性的估计显示出更多的可变性。因此，我们建议首先探索定制的一维模式。
thumbnail下载：
PPT的PowerPoint 幻灯片
PNG放大图片
国际电影节原始图像
图 3. 从模拟数据进行推理的准确性。
A） 在每种选择模式下模拟的 1,000 个重复的箱线图。胡须延伸到 2.5% 和 97.5% 的图块。图上方指示的绘图范围之外的估计值数。估计值通常是无偏的，并且方差较低，但超优势除外。B） 中性重复的单一备选检验的 （p） 与 （p）] 的 Q-Q 图。Inset 显示原始值的相同图。p 值在所有选择模式下都经过了很好的校准。C） 基于每种选择模式的似然比和模拟的选择强度的 AUC 值表。对于 s>0.01，AUC 值接近 1，表示中性和非中性重复之间完美区分。D） （p） 与 （p）] 的 Q-Q 图，用于 和 分布下的统计量。分布校准得更好，尽管它在尾部略微反保守。对于所有模拟，代数为 T = 251，初始条件固定在频率 p = 0.25。
https://doi.org/10.1371/journal.pgen.1011769.g003
验证单个替代 p 值。
接下来，我们报告使用分布获得的单一备选 p 值，如单个备选的 P 值中所述，适用于所有四种模式：加性、隐性、显性和杂合子差异。图 3B 显示了一个 Q-Q 图，其中我们绘制了在中性条件下模拟的重复的经验 p 值与其预期值的关系。对于所有选择模式，这些点都紧密地沿线 y = x，这表明 p 值已得到很好的校准，并且似然比统计量 D 的分布与分布非常近似。
此外，我们还根据受试者-操作员特征曲线的曲线下面积 （AUC） 报告了单一替代测试的性能。AUC 值为 0.5 表示无法区分中性重复和非中性重复，而 AUC 值为 1 表示完全区分。图 3C 报告了在选择参数 s 范围内每种模拟选择模式的 AUC 值。我们观察到，在这些情况下，该方法完全有能力拒绝 的原假设。在隐性选择的情况下，AUC 值较低，因为 251 代可能不足以在低频率下重复以逃脱漂移主导的状态。
选择模式推理的验证。
我们进一步按照 多个替代方案 中概述的程序计算了在多个替代方案的情况下模拟数据集的 p 值。p 值与其预期值的 Q-Q 图如图 3D 所示。在此方案中，分布与统计量拟合良好，尽管它在分布的尾部略微反保守。然而，在其他情况下（参见 S1 文本中的 S.3 教派），分布可能略微过于保守。由于不同的替代选择模式都嵌入在二维参数空间中，因此图 3D 还显示了使用分布从统计量计算的 p 值。但是，我们观察到这些 p 值的校准效果不佳，因此建议使用分布。
然后，我们测试了我们的方法正确识别选择模式的能力。图 4 显示了一个混淆矩阵，该矩阵是通过 p 值显著性阈值 0.05 推断选择模式而得出的。混淆矩阵的每一行表示在特定选择模式下模拟的所有重复。相应列中的数字表示推断特定模式的仿行分数。对角线上的值反映了对正确模式的识别，而非对角线上的值反映了未正确推断模式的仿行。对于 s = 0,6.3% 的重复中立性被拒绝，表明分布是反保守的。对于 s = 0.05，对于所有选择模式的大多数仿行，都会推断出正确的模型。对于较弱的选择 s = 0.01，对于大多数重复，仅推断显性选择和显性选择。然而，在所有选择模式的重复中，超过 50% 的重复中性被拒绝，这表明存在拒绝中立性的能力，但正确推断选择模式的准确性更加有限。
thumbnail下载：
PPT的PowerPoint 幻灯片
PNG放大图片
国际电影节原始图像
图 4. 推断选择模式时的混淆矩阵。
p 值阈值 0.05 用于拒绝中立性。带有橙色边框的单元格表示正确的分类。对于所有选择模式，当 s = 0.05 时，对于大多数仿行，都会推断出正确的模型。对于所有模拟，代数为 T = 251，初始条件固定在频率 p = 0.25。
https://doi.org/10.1371/journal.pgen.1011769.g004
最后，我们研究了选择参数的 MLE 的准确性，前提是模式被正确分类。图 5 显示了被正确分类为加性、显性或隐性选择的重复，以及被错误分类的中性重复的估计值分布。对于较大的 s，正确分类的仿行更紧密地聚集在真实值周围，与使用所有仿行的分布相比，异常值要少得多。对于 s 和中性值的较低值，分类为非中性的仿行的估计值存在偏差。这主要是由于 “赢家诅咒” 效应，其中中立性仅对具有极端 MLE 的重复被拒绝。
thumbnail下载：
PPT的PowerPoint 幻灯片
PNG放大图片
国际电影节原始图像
图 5. 去除推断的 与 true s 的图。
以在多个备选方案中选择正确的模型为条件。在中性情况下，该图显示了分类为特定模式的中性仿行。对于中立性，条带图是从 10,000 个模拟中随机选择的 1,000 个重复。每个条带上的数字表示从 1,000 个重复中正确分类的重复数。对于 s = 0.005 且中性，由于“赢家诅咒”效应，推断的值具有很强的偏差。对于所有模拟，代数为 T = 251，初始条件固定在频率 p = 0.25。
https://doi.org/10.1371/journal.pgen.1011769.g005
选择系数估计的稳健性。
接下来，我们研究了如果我们改变在之前的仿真研究中固定的参数，估计准确性会受到怎样的影响。我们改变了样本的采集次数、每个时间点的单倍体样本数量、有效种群大小以及 HMM 中隐藏状态的数量。我们一次改变一个参数，并将其他参数固定为 K = 11， ， Ne= 10,000，M = 500。我们将所有仿真的选择系数 s 固定为 0.025，将总代数设置为 T = 251，并将初始分布设置为 p = 0.25 的固定频率。图 6 中显示了每种场景的加性、隐性和显性选择的估计选择系数的箱线图。与加法情况下的 [16] 类似，我们发现准确性并不强烈依赖于采样方案，除了非常低的值，例如每个时间点采样 6 个单倍型或仅在轨迹的开始和结束时采样。但是，对于较大的样本量和更多的采样时间，箱线图会更紧密。这意味着由于有限采样引起的方差和真实等位基因频率的方差都会导致选择系数估计的不确定性。同样，对于 250、500、1,000 和 2,000 个隐藏状态，估计的准确性也相当。但是，对于 100 个隐藏状态，估计值明显向下偏斜。因此，我们建议使用 M = 500 来平衡精度和计算效率。
thumbnail下载：
PPT的PowerPoint 幻灯片
PNG放大图片
国际电影节原始图像
图 6. 不同参数范围的箱线图。
胡须延伸到 2.5% 和 97.5% 的图块。我们各不相同：A） 每个时间点采集的样本数，B） 采样的时间点数，C） 有效种群规模，D） HMM 中隐藏状态的数量。除最低参数外，选择系数的估计值是无偏的。箱线图的宽度随着样本数量和采样时间点的增加而减小。对于所有模拟，代数为 T = 251，初始条件固定在频率 p = 0.25。当不变化时，使用的参数为 K = 11 采样时间，每个时间点的样本，Ne= 10,000，M = 500 个隐藏状态。
https://doi.org/10.1371/journal.pgen.1011769.g006
为了调查选择模式不正确时估计值的准确性，我们在错误指定模式的情况下重新分析了数据的子集。我们考虑两种情况：（1） 在加法、显性、隐性、过度或欠显性选择下进行模拟并在加性选择下进行分析，以及 （2） 在加性选择下进行模拟并在加性、显性、隐性或杂合子差异下进行分析。图 7A 和 7B 描述了在这两种错误指定情景下 MLE 的箱线图。
thumbnail下载：
PPT的PowerPoint 幻灯片
PNG放大图片
国际电影节原始图像
图 7. 使用不正确的选择模式进行推理。
A） 在每种选择模式下模拟的 1,000 个重复的箱线图，并使用加性选择进行分析。这些估计大多不准确。B） 在加法选择下模拟的 1,000 个重复的箱线图，并在每种选择模式下进行分析。同样，如果众数不正确，则估计值不准确。对于这两组箱线图，须线都扩展到 2.5% 和 97.5% 的图块。C） 使用假设加性选择获得的似然比在每种选择模式下模拟的数据的 AUC 值表。与在正确的选择模式下进行分析相比，值较低，但仍显示出拒绝中立性的巨大力量。D） 在加法选择下使用假设每种可能选择模式获得的似然比模拟数据的 AUC 值表。这些值在选择模式中非常相似。对于所有模拟，代数为 T = 251，初始条件固定在频率 p = 0.25。
https://doi.org/10.1371/journal.pgen.1011769.g007
当在加性 EM-HMM 下模拟每种选择模式和分析时，所有模式的估计值都是偏倚的，最强烈地用于高 s 下的占主导地位选择。经历强显性选择的等位基因很快达到 p = 0.5 的稳定频率。由于等位基因的频率不会继续增加，因此加性 EM-HMM 低估了真正的选择强度。在加性选择下进行模拟并在每种选择模式下进行分析的情况下，估计值对隐性和显性没有强烈的偏差，而是在高选择强度下对杂合子差异的反向符号。在这种情况下，等位基因频率增加到 p = 0.5 以上，因此杂合子差异 EM-HMM 从显性选择变为次显性选择。总体而言，估计中的偏差表明，使用不正确的单参数 EM-HMM 可提供不准确的 s 估计。
我们进一步调查了在每种错误指定情况下，单个替代方案的 P 值中描述的单一替代方案测试框架中拒绝或接受中立性的能力。得到的 AUC 如图 7C 和 7D 所示。对于在加法选择下分析其他模式的情况，AUC 值略低于图 3C 中的值，但仍大大大于 0.5。对于使用不正确的选择模式分析的加法数据集，AUC 值在不同模式之间没有太大差异。在这两种错误指定的情况下，单参数加性 EM-HMM 具有足够的能力来准确拒绝或无法拒绝中高选择系数的中性。因此，如果目标只是识别非中性进化，那么仅使用加性 EM-HMM 分析给定数据就足够了，但为了准确表征选择系数，需要使用正确的模式。
推断有效种群规模。
接下来，我们探索了应用我们的方法来推断有效种群规模 Ne的基础人群。方程 2 和方程 3 中提供的 EM-HMM 算法的更新规则的推导方法并不容易产生 N 的更新规则e.因此，我们改用基于网格的方法来估计它。具体来说，我们计算了在中性 HMM 下观察到给定复制的可能性，并且给定的种群大小为 Ne.为了组合跨重复或基因座的功效，我们计算 N 个网格上所有重复的对数似然之和e值，然后在这些网格值之间进行插值，并使用 N 的值e，它将这个复合似然曲面最大化为我们对 N 的估计e.
我们注意到，当模拟具有较大 N 值的数据时e，则生成的似然表面通常非常平坦，这使得估计具有挑战性。为了抵消这一点，我们将可能性的条件定为在任何时间点观察到至少一个多态性样本，方法是将可能性除以 。这会对 N 的高值进行惩罚e并导致更多的峰值似然曲面。S1 Text 的图 Q 显示了几个示例复合似然面。即使在这个调节程序之后，表面仍然相当平坦，但它们确实允许确定一个明显的最大值。此外，我们发现 Ne当计算可能性的初始条件固定为均匀分布时，估计值最准确，而不是像初始等位基因频率分布的估计那样估计。一种可能的解释是，估计初始条件的过程受 N 选择的影响e，这反过来又可能会使估计值产生偏差。
我们在 25 代的情景下模拟了 25 批 10,000 个中性重复，初始频率固定为 p = 0.25，以及模拟研究中使用的抽样方案。我们估计了 Ne对于每个批次，请使用上述过程。图 8 显示了推断的 N 的箱线图e对于每个批次。Ne从基于格网的过程估计的值紧密聚集在 True 值周围，尽管在所有情况下都略微向上倾斜。更强大的 N 估计方法e例如，在当代数据中使用 IBD 片段 [44]。然而，我们认为，在分析时间序列遗传数据时，基于网格的 HMM 程序的性能是可以接受的，并产生最合适的 N 估计e用于下游分析。
thumbnail下载：
PPT的PowerPoint 幻灯片
PNG放大图片
国际电影节原始图像
图 8. 估计有效种群规模 .
N 的箱线图e对于在中立性下模拟的数据，使用基于网格的 HMM 过程进行估计。对于每个模拟值 Ne，显示的 25 个估计值中的每一个都基于一批 10,000 个仿行的复合似然。Ne估计值略微向上偏倚，方差较低。对于所有模拟，代数为 T = 251，初始条件固定在频率 p = 0.25。胡须延伸到 2.5% 和 97.5% 的图块。
https://doi.org/10.1371/journal.pgen.1011769.g008
来自英国的古代 DNA 数据集
GB aDNA 数据集的描述。
在证明了我们的方法能够准确表征模拟数据中的选择后，我们将我们的方法应用于从人类古代 DNA 获得的时间序列遗传数据。为此，我们从Allen Ancient DNA Resource（AADR）v54.1版中的个体子集中提取了基因型信息，该资源库是一个经常更新的存储库，旨在包含当前发布的大多数古代DNA数据集[33]。
我们的方法假设数据来自单个泛种群。因此，我们遵循类似于 [45] 的基本原理，并将我们的分析限制在过去 4,450 年中来自英国的样本。我们手动移除了非来自大陆的样本。受限于这个地理区域和时间窗口，我们分析的数据不太可能实质性地违反单一泛种群的假设，因为据估计，上一次进入英国的重大混合事件发生在 4,450 年前 [46， 47]，尽管已经提供了一些更近的基因流动的证据 [48， 49]. S1 文本中的图 R 显示了一张带有样本位置和时间的地图。此外，S1 文本中的图 S 描绘了两个 PCA 图，表明样本在全球范围内与现代欧洲个体聚集，但在局部尺度上没有表现出强结构的证据。
当分析使用 1240K 捕获测定法进行基因分型的样品和使用全基因组测序进行基因分型的样品时，我们注意到了虚假的选择信号，参见 S1 文本中的图 V。因此，作为一种保守的方法，我们只分析了使用 1240K 捕获进行基因分型的样本，并从我们的数据集中排除了该地理区域和时间范围内的 174 个全基因组样本。我们注意到，这种保守的方法也删除了当今的样本。因此，我们得到的数据集，以下简称GB aDNA数据集，包括504个使用1240K测定法进行基因分型的古老伪单倍体样本，跨越125代，世代时间为28.1[50]。同一代中的样本被分箱在一起，以形成用于分析的最终数据集。这些个体是 [46， 48， 49] 中发表的个体的子集。
我们还对数据集中的每个 SNP 应用了三个过滤器。首先，每个 SNP 必须在两个或多个时间点具有基因分型样本。其次，每个 SNP 必须有超过 50 个（占 504 个的 10%）样本的基因分型。第三，在一段时间内合并给定 SNP 的所有样本时，次要等位基因频率必须大于 0.05。我们期望只有通过这些过滤器的 SNP 才能产生可靠的选择信号。总的来说，在可用的 1,150,638 个 SNP 中，有 743,417 个 （64.6%） 通过了这三个过滤器并用于最终分析。
数据匹配的模拟。
为了评估我们的方法在 GB aDNA 数据集的特定背景下的准确性，我们模拟了两个与 GB aDNA 数据集的采样方案和时间框架匹配的数据集：在第一个数据集中，我们称为 IBDNe 数据集，我们使用随时间变化的有效种群规模的历史来模拟数据， 使用先前从 UK10K 数据集推断的英国人口的种群规模历史 [44]。在第二个数据集中，我们称之为 const-Ne 数据集，我们在单个常数 N 下进行模拟e根据 GB aDNA 数据集估计（参见 GB aDNA 数据集中重要的选择信号）。我们在本节中显示了 IBDNe 数据集的结果，并显示了常数为 N 的数据集的结果e在 S1 文本的 S.6 教派中。
我们使用选择系数模拟了 T = 125 代的特定选择模式下的等位基因频率轨迹，并根据 GB aDNA 数据集的采样时间和大小对给定每个轨迹的单倍型进行采样。对于 IBDNe 数据集，我们使用 graphreader 工具 （https://www.graphreader.com） 提取 N 的值e在 [44] 中给出的图 4A 中每一代的 GB aDNA 数据集中样本所跨越的时期。然后将这些值用作时变 Ne在 Wright-Fisher 模拟中。S1 文本的图 O 显示了提取的 Ne值。此外，我们随机省略了采样的单倍型，其概率等于 GB aDNA 数据集中随机选择的 SNP 的缺失分数，以模拟相同程度的缺失数据。我们在 S1 文本的图 T 中提供了 GB aDNA 数据集中每个 SNP 缺失的采样单倍型分数的直方图。为了进一步确保模拟的重复与 GB aDNA 数据集匹配，我们应用相同的三个基于 SNP 的筛选标准，并且只保留通过所有筛选的重复。我们在中性下模拟了 10,000 次重复，并在五种单参数选择模式中的每一种模式下模拟了 10,000 次重复。
为了生成每个模拟重复的初始频率，我们首先估计了参数以及通过过滤器的 GB aDNA 数据集中每个 SNP 在中性 （s = 0） 下初始频率的 beta 分布，如初始等位基因频率分布的估计中所述。S1 文本中的图 U 显示了为每个 SNP 估计的初始分布的平均值的直方图。对于每个模拟的重复，我们然后随机均匀地选择一个 SNP，并将重复的初始频率设置为等于所选 SNP 的平均值。此过程可确保模拟数据的初始频率分布与 GB aDNA 数据集紧密匹配，并捕获任何潜在的偏差，例如，由于确定 1240K SNP 集。
为了启发式地解释模拟数据中的可变种群大小历史，我们遵循类似于 [51] 的策略：我们使用中性重复来估计共享常数，使用推断有效种群大小中描述的过程，并在估计每个重复的选择系数时使用推断。然后，我们使用 EM-HMM 计算每个重复的选择系数的 MLE，分析中的选择模式与模拟模式匹配，并在图 9A 中显示 MLE 的分布。与 Simulation study 中的模拟数据集一样，选择系数的估计值在很大程度上是无偏的。但是，与模拟数据集不同的是，对于强选择，数据匹配的重复具有轻微的向下偏差。这可能是因为数据匹配的模拟包含的代数是图 3A 所示模拟研究的一半。T = 101 代的模拟数据集在高 s 处显示出类似的向下偏差，但不太明显——参见 S1 文本中的图 C。此外，尽管时间跨度较短，但对 Underdominant 的估计方差低于 Simulation 研究中模拟数据集的方差。
thumbnail下载：
PPT的PowerPoint 幻灯片
PNG放大图片
国际电影节原始图像
图 9. 数据匹配模拟的推理准确性。
A） 所有单参数选择模式的 MLE 箱线图。每个箱线图显示模拟的 10,000 个随机仿行的 1,000 个。胡须延伸到 2.5% 和 97.5% 的图块。对于小 s，估计值在很大程度上是无偏的，而对于大 s，估计值略微向下偏倚。B） （p 值） 与 （p 值） ] 的 Q-Q 图，用于使用分布获得的单一替代检验。插图显示原始 p 值与预期 p 值的对比。对于模拟数据集，p 值得到了很好的校准。C） 使用每种单参数选择模式的似然比进行数据匹配模拟的 AUC 值表。D） 使用 和 分布的多重替代似然比统计量的 （p 值） 与 （p 值）] 的 Q-Q 图。该分布提供了一个很好的拟合，但尾部略微反保守。对于所有模拟仿行，代数为 T = 125，即 N 的值e在每一代中都源自 [44]，初始频率与 GB aDNA 数据集相匹配。
https://doi.org/10.1371/journal.pgen.1011769.g009
接下来，我们应用我们的程序来测试单个替代方案，并应用我们的程序来推断 IBDNe 模拟的选择模式。对于单个备选方案的检验，图 9B 显示了根据似然比统计量 D 计算的 p 值的 Q-Q 图，假设分布与其预期值相反。我们再次观察到这些 p 值是经过良好校准的。与单一备选 p 值的验证一样，我们计算了 AUC 值以评估单一备选测试对 IBDNe 模拟的功效，如图 9C 所示。对于 ，AUC 值介于 0.8 和 1 之间，但主导选择除外。一般来说，AUC 值低于图 3C;例如，在优势情况下，当 s = 0.05 时，IBDNe 仿真的 AUC 为 0.81，而仿真研究的 AUC 为 0.99。这可能是由于我们用于 IBDNe 匹配模拟的初始分布，该分布在接近 p = 0 时具有更大的权重，并且代数减少。这两个特性都导致等位基因频率的累积变化较小，因此区分中性重复的能力降低。
我们还应用了 多个备选方案 中介绍的多个备选方案框架。图 9D 显示了 IBDNe 数据集中性重复的 p 值，使用统计量 under 和 distributions 计算，并根据预期值绘制。使用分布时，p 值通常校准良好，尽管与选择模式推理验证中的模拟一样，p 值在分布的尾部略微反保守。
图 10 总结了使用 Multiple alternatives中描述的过程对 IBDNe 模拟的所有重复的选择模式的推断;再次使用 p 值阈值 0.05 来拒绝中立性。对于 s = 0.05，我们成功地拒绝了大部分重复的中性，在加法选择下进行模拟时，高达 95% 的情况，但仅对 30% 到 55% 的重复推断出正确的模型。对于较低的选择强度 s = 0.01，仅在一小部分模拟中被拒绝中性。与选择模式推理验证中介绍的模拟研究一样，我们发现我们可以检测到非中性进化，但推断正确选择模式的能力是有限的。请注意，我们无法拒绝 94% 的中性重复，而不是 95%，因为在分布下 p 值略微反保守。
thumbnail下载：
PPT的PowerPoint 幻灯片
PNG放大图片
国际电影节原始图像
图 10. 用于推断应用于数据匹配模拟的选择模式的过程的混淆矩阵。
使用了 0.05 的 p 值阈值。每个单元格表示被分类为特定模式的仿行分数。性能比 sec：simulation-study 中的模拟差 - 仅对 s = 0.05 的多个重复推断出正确的模型。对于所有模拟，在给定的选择模式和强度下模拟 10,000 次重复，代数为 T = 125，即 N 的值e在每一代中都源自 [44]，初始频率与 GB aDNA 数据集相匹配。
https://doi.org/10.1371/journal.pgen.1011769.g010
最后，图 11 显示了推断的选择系数的分布，条件是推断多个选择中的正确选择模式，以及归类为非中性的中性重复的推断选择系数的分布。对于 s > 0.01，估计的选择系数与无条件估计相似，这表明我们的模型推理程序不会严重偏倚该参数范围内的估计。但是，对于 s = 0.005 和 s = 0.01，大多数推断系数都高于真实值。与选择模式推理验证中的模拟数据一样，我们观察到较低选择系数的 “赢家诅咒 ”现象，其中只有等位基因频率变化极端的重复被归类为非中性，因此也很大。
thumbnail下载：
PPT的PowerPoint 幻灯片
PNG放大图片
国际电影节原始图像
图 11. 对于数据匹配的模拟，去除推断的与真实 s 的图。
以推断中性重复的非中性为条件，以及推断非中性重复的多个备选方案中的正确选择模式为条件。每个条带图用于从 10,000 个模拟中随机选择的 1,000 个重复。每个条带上方的数字表示正确分类的重复数。对于 s = 0.005 和 s = 0.01，由于“赢家诅咒”，推断值会向上偏斜。对于所有模拟，代数为 T = 125，即 N 的值e在每一代中都源自 [44]，初始频率与 GB aDNA 数据集相匹配。
https://doi.org/10.1371/journal.pgen.1011769.g011
GB aDNA 数据集中的重要选择信号。
然后，我们应用我们的 EM-HMM 算法和程序来推断 GB aDNA 数据集中通过过滤器的所有 743,417 个 SNP 的选择模式。例如，在 [44] 中所示，英国有效种群规模的历史随时间而变化。为了启发式地解释这一点，我们再次镜像了 [51] 中的策略，并首先使用 推断有效种群大小 中描述的程序估计了 SNP 之间共享的单个恒定有效种群大小，结果为 。对于每个 SNP，我们然后使用 EM-HMM 估计所有单参数选择模式的 MLE，将推断固定为恒定的有效种群大小。如数据匹配模拟中所示，这种启发式方法用于考虑时变 Ne产生选择系数的准确估计值和经过充分校准的 p 值。在这里，我们主要描述加法模式下的结果以及推断选择模式的过程的结果。其他单参数模式的结果，以及所有模式的 p 值的 Q-Q 图，可以在 S1 Text 的 S.9 节中找到。
图 12 显示了使用加性 EM-HMM 计算的单一备择似然比检验的 p 值的曼哈顿图。我们还指出了通过应用 Benjamini-Hochberg （BH） 程序 [52] 以错误发现率 （FDR） 校正多重测试而获得的显著性阈值。在 S1 Text 的 S.10 节中，我们在排列数据集中的采样时间时计算相同的 p 值，这表明我们观察到的低 p 值的富集是数据中的可靠信号。此外，我们在 2 、 5 和 6 号染色体上观察到低 p 值的簇。这种 p 值的聚集是意料之中的，因为由于遗传搭便车 [53]，靠近实际选择靶标的 SNP，因此与靶标的连锁不平衡 （LD） 也将表现出非中性动力学。
thumbnail下载：
PPT的PowerPoint 幻灯片
PNG放大图片
国际电影节原始图像
图 12. 加法 p 值的曼哈顿图。
在 GB aDNA 数据集中所有 SNP 的加性选择模式下，从似然比检验获得的 p 值的曼哈顿图。显著性阈值是通过 Benjamini-Hochberg 程序获得的，FDR 为 。我们在 2 、 5 和 6 号染色体上观察到显著 p 值的簇，以及几个孤立的信号。
https://doi.org/10.1371/journal.pgen.1011769.g012
除了这些簇之外，我们还观察到几个 p 值超过 BH 阈值的孤立 SNP，但没有周围的 SNP 显示选择的证据，这种模式在遗传搭便车下是意料不到的。因此，我们认为这些孤立的 SNP 不对应于真正的选择信号，它们可能是数据集中的伪影。然而，其他几个区域的 SNP 的 p 值超过 BH 阈值，并且附近的多个 SNP 表现出低 p 值，从而在一定程度上支持这些 SNP 可能对应于真实的选择信号。图 13A 和 13B 分别显示了基因组区域中的 p 值，其中 5 号染色体上显著 SNP 周围的 p 值较低，而 SNP 周围的区域 7 上 p 值呈虚假的低。
thumbnail下载：
PPT的PowerPoint 幻灯片
PNG放大图片
国际电影节原始图像
图 13. 围绕重要 SNP 的 P 值和频率轨迹。
A） 5 号染色体基因组区域 p 值的曼哈顿图，以 p 值超过 BH 阈值的 SNP 为中心。周围的 SNP 表现出低 p 值，正如在遗传搭便车下所预期的那样。后处理的 p 值超过相应的 BH 阈值。B） 7 号染色体基因组区域 p 值的曼哈顿图，以 p 值超过 BH 阈值的 SNP 为中心。周围的 SNP 没有显示选择的证据。后处理的 p 值不显示任何显著信号。C） 图 13A 中 p 值最低的 SNP 的 20 个 SNP 的分箱等位基因频率轨迹。附近的 SNP 显示相关的等位基因频率变化，表明遗传搭便车和真实信号。点的大小表示样本的数量;色相表示基因组位置：红色较小，蓝色较大。D） 图 13B 中 p 值最低的 SNP 的 20 个 SNP 的分箱等位基因频率轨迹。附近的 SNP 不表现出相关的等位基因频率变化，表明先导 SNP 处存在虚假信号。
https://doi.org/10.1371/journal.pgen.1011769.g013
为了去除虚假的SNP，同时将显著的SNP保留在对周围SNP有额外支持的区域中，我们使用Brown方法[54]的修改版本对p值进行后处理，以组合非独立的p值。具体来说，我们考虑围绕每个分析的 SNP 的 50 个 SNP 的重叠窗口。对于这些重叠窗口中的每一个，我们计算 p 值的对数的负和，包括焦点 SNP。然后，我们将缩放分布的参数拟合到这些总和，并使用此拟合分布来计算每个 SNP 的后处理 p 值。我们将 BH 方法应用于这组新的 p 值，以获得第二个 BH 阈值。分离的 SNP，例如图 13B 中的 SNP，不再超过相应的 BH 阈值，而包含显著 SNP 和其他支持的区域，例如图 13A 中的区域，具有超过相应 BH 阈值的宽峰。
应用此后处理程序后，我们将显著 p 值分组到不同的显著性区域。对于被视为显著的区域，我们要求该区域内至少有一个 SNP 必须同时具有原始 p 值和后处理 p 值，每个值都超过其各自的 BH 阈值。后处理过程拓宽了 p 值峰值;因此，我们将超过相应 BH 阈值的每个连续后处理 p 值块作为单独的候选区域。根据此标准，对于一维加性 EM-HMM 算法，有 8 个不同的候选区域。表 1 列出了这些区域（在 hg19 坐标中）、与该区域重叠的任何基因、对应于该区域中最小 p 值的 SNP 的 rsID、该 SNP 的参考和替代等位基因、每个区域中重要 SNP 的数量（预处理和后处理）、该区域中最小 p 值的负对数以及加法 对于具有最低 p 值和置信区间的 SNP 处的衍生等位基因。通过使用 模拟 1,000 个重复来获得置信区间，这些重复匹配每个先导 SNP 的采样方案、估计的初始频率和估计的选择系数，然后使用偏差校正所有值并取模拟重复的 0.025 和 0.975 分位数。此外，我们还在 S1 文本的 S.11 节中展示了每个重要基因组区域中 p 值和分箱等位基因频率轨迹的曼哈顿图。
thumbnail下载：
PPT的PowerPoint 幻灯片
PNG放大图片
国际电影节原始图像
表 1. 在加性选择下被鉴定为显著的基因组区域。
https://doi.org/10.1371/journal.pgen.1011769.t001
重要基因组区域的功能性遗传变异。
除 TFR2 外，我们通过应用加性 EM-HMM 检测到的所有具有重要性的基因组区域之前都已被确定为选择靶标，尽管并非所有基因组区域都在英国的人群中被确定为如此。在以下段落中，我们将在相关文献的背景下讨论这些地区中的每一个。
在我们的分析中，最强的信号是特征明确的LCT基因座，它已被确定为英国[45,55,56]和更广泛的西欧人群[57–59]的选择目标。 我们分析中的先导 SNP rs4988235 以前已被确定为 LCT 基因座最强的选择信号 [45， 58， 60]，并且该 SNP 的衍生等位基因与消化乳糖酶进入成年期的能力有关 [61]。我们估计的加性选择系数为 0.080 （CI： [0.057， 0.099]），与文献中提供的其他估计值一致 [45， 57， 60]。尽管该 SNP 的衍生等位基因与乳糖酶持久性有关，但最近的研究认为，牛奶消费的引入早于该等位基因频率的增加，并且最近强烈的选择压力可能是由于后来的饥荒造成的，其中等位基因被证明是有益的 [62， 63]。我们的方法还在 LCT 基因座附近确定了一个具有重要意义的次级区域，这可能是遗传搭便车的结果。
我们的两个候选区域 TLR10/1/6 和 HLA 的基因组变异参与免疫调节。TLR10/1/6基因簇中的多态性与多种癌症以及肺结核和麻风病的发病率有关[64–67]。在先前使用古代 DNA 的研究中，TLR10/1/6 位点已被确定为选择靶标 [59]，尽管这项工作汇集了来自欧亚大陆西部的样本。同样，使用当今个体的数据集，[68] 发现 TLR 基因（TLR1 除外）经历了强烈的负选择，并且 TLR10/1/6 簇最近在非非洲人群中经历了选择。
HLA 基因座跨越 6 号染色体上的一个大区域，编码一组对先天免疫系统功能至关重要的高度多态性基因。以前的研究已经在该区域发现了多个选择信号 [45， 56]，我们发现的 SNP 与其中一个信号重叠。HLA 区域内的单个基因座被认为处于平衡或频率依赖性选择下，以增加等位基因多样性 [69， 70]。出于这个原因，这里和文献 [45， 56] 中对正选择的基因组扫描在这些基因座上发现信号有点令人惊讶。此外，我们对杂合子差异的扫描，包括显性，一种平衡选择模型，显示 HLA 区域没有显着信号。可能的解释是，几个频率变化的等位基因被检测为正选择，但这里考虑的短时间范围不足以检测长期平衡选择下的等位基因。
SLC45A2位点也被公认为欧洲人群的选择靶点[71,72]。 该位点的多态性与头发和皮肤色素沉着的差异有关[73,74]。使用来自西欧和英国的古代 DNA 数据的研究特别确定了 SLC45A2 位点作为选择的目标 [45， 59]。此外，我们在先导SNP上的估计选择系数与从现代和古代DNA分析中获得的值相匹配：我们的MLE为s = 0.049（CI：[0.029,0.066]），[72]估计s在0.04至0.05的范围内，[45]估计s = 0.043。
该基因SLC22A4包含在较大的 IBD5 位点中，该位点由一组与胃肠道疾病（如克罗恩病）相关的多态性基因组成 [75]。[76] 发现SLC22A4的遗传变异增加了抗氧化剂麦角硫因的吸收，并显示出正选择的信号，这可能是由于新石器时代早期农民饮食中麦角硫因的含量低。此外，他们认为，与克罗恩病相关的变异可能通过遗传搭便车增加了频率。SLC22A4/IBD5 位点也使用西欧人的古代 DNA 被确定为选择靶点 [59]。
TFR2 基因座的遗传变异尚未被确定为使用当代或古代基因组样本的选择靶标。S1 文本中的图 AJ 显示，先导 SNP 和周围 SNP 的等位基因频率协同移动，表明这可能是真正的选择目标，而不是错误信号。TFR2 位点突变可导致 3 型遗传性血色病，其特征是全身铁水平异常高 [77]。此外，我们确定为选择下的先导 SNP 突变的单倍型与帕金森病相关 [78]。
HYDIN 基因座包含数据集中 p 值最显著的 SNP。然而，在检查 S1 文本中图 AK 中提供的该基因组区域中的 p 值时，我们观察到该基因座确实包含三个 p 值非常低的 SNP，但也包含两个没有任何 SNP 的 100 kbp 区域。这些空区域是我们过滤程序的结果，该程序去除了大量次要等位基因频率低于 0.05 的 SNP。此外，除最后一个时间点外，HYDIN 位点的重要 SNP 在任何地方都具有极低的分档次要等位基因频率，参见 S1 文本中的图 AK。由于 16 号染色体上的基因 HYDIN 在 1 号染色体上有一个假基因 [79]，这些不寻常的模式可能是序列读数映射错误的结果，因此我们认为 HYDIN 基因座的选择信号是虚假的。
最后，我们更明确地将我们的分析结果与 [45] 中的结果进行比较，[45] 是最近的一项研究，该研究在与我们类似的时间数据集中确定了选择目标。作者分析了来自英格兰本土的 AADR 的现代和古代 DNA 样本，其历史可以追溯到 4,450 BP 以下，并在五个基因组区域中发现了选择信号。在我们的研究中还确定了三个区域 （LCT、SLC45A2、HLA），在我们的分析中仅确定了三个区域 （TLR10/1/6、SLC22A4、TFR2），两个区域 （DHCR7、HERC2） 仅在 [45] 中被识别。在两项研究中确定的基因座中，最显著 SNP 的选择系数估计值基本一致：LCT、SLC45A2 和 HLA 区域的系数在我们的研究中分别为 0.080 （CI： [0.057， 0.099]）、0.049 （CI： [0.029， 0.066]） 和 0.046 （CI： [0.025， 0.061]），在 [45] 中分别为 0.064、0.043 和 0.046。S1 文本中的图 AL 和 AM 显示了使用我们的加法 EM-HMM 和分箱等位基因频率轨迹计算的 p 值，这两个区域在 [45] 中被确定为显着，但在我们的分析中并不显着。在 11 号染色体上的 DHCR7 基因座处，后处理的 p 值接近超过 BH 阈值，但没有原始 p 值达到显著性。15 号染色体上的 HERC2 位点也表现出低 p 值，尽管程度低于 DHCR7。与我们仅使用 1240K 分析进行基因分型的样本的保守方法相比，[45] 还分析使用全基因组数据进行基因分型的样本（包括当今的样本），使用一种有可能检测加性系数随时间变化的选择的方法，并在组合相邻 SNP 的信号时使用不同的后处理。因此，不希望信号完美对齐。
推断 GB aDNA 数据集中的选择模式。
除了加性 EM-HMM 之外，我们还分析了其他三种单参数选择模式（显性、隐性和杂合子差异）下的 GB aDNA 数据集，参见 S1 文本中的 S.9 节，其中我们为每种单参数选择模式提供了完整的曼哈顿图，以及类似于表 1 的表格.此外，我们提供了一个表格，该表是 S1 文本图 AB 中前 1% SNP 的所有单参数模式的对数似然比统计之间的全基因组 Spearman 秩相关系数。显性和隐性 EM-HMM 都与加性 EM-HMM 具有高度相关性，并且标识与加性分析相同的区域，尽管 LCT 基因座在显性分析中没有分成两个不同的区域。
相比之下，杂合子差异 EM-HMM 与其他模式的相关性较低，并且仅将 LCT、两个 SLC 位点和 HYDIN 识别为显著性。这种模式无法检测到 HLA 基因座的选择起初有些令人惊讶，因为平衡选择或频率依赖性选择应该导致类似于过度选择的动态。然而，为这些类型的选择提供有力支持的特征是在中性下保持中频时间比预期更长的等位基因，因此这里考虑的短时间范围可能不够。
我们还使用了 多种选择中详述的程序来推断每个 SNP 最可能的选择模式。我们计算了 GB aDNA 数据集中所有 SNP 的统计数据，并使用分布来获得每个基因座的 p 值。我们应用基于 Brown 方法的程序来识别重要区域，使用 FDR 为 的 BH 阈值。原始 p 值超过 BH 阈值的重要区域中的 SNP 被归类为具有最高对数似然的单参数选择模式。
模式推理过程将与加性单一替代过程相同的基因座集标识为显著。图 14 显示了在 LCT 基因座处得到的原始 p 值，其中重要的 SNP 由其推断的选择模式着色。在该区域超过 BH 阈值的 68 个 SNP 中，有 30 个被归类为加性SNP，但其他 SNP 显示出不同的选择模式。该区域中的大多数基因座被归类为加性基因座的事实可能表明支持 LCT 在加性选择下进化的假设;然而，先导 SNP rs4988235 的推断模式是显性的。乳糖酶持久性作为显性性状发挥作用 [80]，进一步支持了先导 SNP 显性选择的推断。此外，[45] 中提出的方法可以模拟非常量选择系数，作者提供了证据，证明 LCT 基因座的选择随着时间的推移而减弱;一种类似于恒定显性选择的动态。我们确实警告不要过度解释这些结果，因为数据匹配模拟中的模拟研究表明，在该数据集中，即使是确定恒定的选择模式也具有挑战性，因此可能需要更大的样本数量和密度才能进行准确分类。最后，我们注意到，对于每个已识别区域的先导 SNP 处具有最高对数似然比的单参数模式如下：LCT – 显性，TLR10/1/6 – 加法，SLC45A2 – 杂合子差异，SLC22A4 – 显性，HLA – 隐性，和 TFR2 – 加法。
thumbnail下载：
PPT的PowerPoint 幻灯片
PNG放大图片
国际电影节原始图像
图 14. LCT 轨迹处原始 p 值的曼哈顿图。
使用多种选择中描述的选择模式的识别程序计算 P 值，并且显着的 SNP 通过推断模式着色。先导 SNP 由较大的星形标记表示。该区域中的大多数 SNP 被归类为加性 SNP，尽管先导 SNP 被归类为显性 SNP。
https://doi.org/10.1371/journal.pgen.1011769.g014
家养马的毛色位点 ASIP
数据集的描述。
在本节中，我们将 out 方法应用于 [34] 中提供的数据集，其中作者从欧亚大陆的一组样本中提取了影响驯养马毛色的六个基因座的古代 DNA，并在 ASIP 和 MC1R 基因座找到了选择的证据。具体来说，我们将我们的方法应用于 ASIP 基因座。图 15A 显示了该基因座上衍生等位基因的样本等位基因频率随时间的变化。样品表现出衍生等位基因的频率急剧增加，然后在大约 0.5 的频率处趋于平稳。潜在的采样等位基因计数可以在 S1 Text 的表 E 中找到。
thumbnail下载：
PPT的PowerPoint 幻灯片
PNG放大图片
国际电影节原始图像
图 15. ASIP 基因座非中性的频率轨迹和证据。
A） 在 ASIP 基因座处随时间推移的等位基因频率。点的大小表示样本的数量。B） & C） 与ASIP基因座匹配的1,000个模拟中性重复的统计直方图，使用B）在中性下估计的初始频率或C）使用杂合子差异模式估计的初始频率。在这两种情况下，原始数据集的统计数据（由红色虚线表示）都高于任何模拟重复，从而提供了反对中立性的有力证据。
https://doi.org/10.1371/journal.pgen.1011769.g015
该数据集已在几项研究中在不同选择模式下进行了重新分析，结果不同;[13] 分析了隐性选择下的 ASIP 基因座，但没有找到选择的证据，[28] 找到了压性选择的证据，而最近 [81] 分析了隐性选择下的基因座并得出结论，它不在选择下。已知 ASIP 位点通过隐性机制起作用——具有两个衍生等位基因拷贝的马是黑色的，否则它们是月桂色的 [82]。因此，ASIP 基因座显示显性选择的证据而不是隐性选择的证据有点出乎意料。
ASIP 基因座的单备选和多备选推理。
我们应用我们的 EM-HMM 来估计在每种单参数选择模式下 ASIP 基因座初始分布的选择系数和参数。在 [81] 之后，我们假设 Ne= 16,000 和生成时间为 8 年。对于单参数模式，我们使用基于 的对数似然比检验来计算 p 值。我们在表 2 中报告了选择系数、对数似然差异和 p 值的结果估计值。我们观察到，优势选择的证据最强，单个备选 p 值为 。与 [13] 和 [81] 相反，我们发现 ASIP 基因座显示了隐性选择的证据，尽管 p 值 （） 不是很强。
thumbnail下载：
PPT的PowerPoint 幻灯片
PNG放大图片
国际电影节原始图像
表 2. 不同选择模式下 ASIP 基因座的结果。
https://doi.org/10.1371/journal.pgen.1011769.t002
尽管所有选择模式的估计选择系数都较低（例如 对于杂合子差异），该数据集包含大约 2,500 代，是 GB aDNA 数据集的十倍多。这增加了检测较弱选择的能力;例如，S1 文本图 B 中 100 代模拟的 s = 0.005 的 AUC 值仅略高于 0.5，表明检测选择的能力最小，而 S1 文本图 G 中绘制的 1,000 代模拟中的 AUC 值均高于 0.9。
除了单个备选方案外，我们还计算了多个备选方案的检验统计量（参见多个备选方案），但使用参数化 bootstrap 模拟来评估其显着性。为此，我们模拟了两组 1,000 个中性重复，将采样方案和世代数与 ASIP 基因座的原始数据相匹配。一组使用在中性 EM-HMM 下估计的初始频率进行模拟，另一组使用来自杂合子差异 EM-HMM 的初始频率进行模拟，这是单参数模式中可能性最高的方式。我们模拟这两组以涵盖不同的合理场景。图 15B 和 15C 分别显示了使用中性初始频率和杂合子差异初始频率的 1,000 个模拟重复的统计量直方图，原始数据的统计量由垂直线表示。在这两种情况下，原始数据的统计量都比任何模拟的重复都大得多，这表明非中性动力学的重要证据。优势选择（具体来说，具有 s 的杂合子差异单参数模式1>0） 具有最高的对数似然。因此，ASIP 基因座的数据支持优势选择是最可能的模式。此外，隐性选择在所有单参数选择模式中具有最低的对数似然。由于数据显示频率急剧增加，然后是 0.5 附近的平台，因此预计优势地位得到了最有力的支持。
ASIP 基因座衍生等位基因的遗传机制是隐性的，但数据强烈表明，优势是最强的，而隐性选择是最不可能的选择模式。对于这种差异，我们提出了以下两种可能的解释。首先，衍生等位基因的选择系数可能在马驯化过程中降低，并且在初始增加后可能起作用某种形式的平衡选择。[81]中的点估计确实表明，在驯化之前有较强的选择，之后是较弱的选择;这也与 [83] 中的发现一致。为了进一步探索这一假设，我们在前三个时间点之后截断了 ASIP 基因座的数据，以仅分析频率增加的时间段，并报告表 2 中的结果。我们发现，事实上，隐性选择的证据比过度选择的证据更强，然而，加性和显性选择也很重要，并且具有更强的 p 值。对我们研究结果的第二种可能的解释可能是上位性，因为ASIP基因座的衍生等位基因与另一个已被证明处于选择的被选中的被毛着色基因座MC1R基因座[82]存在上位性相互作用，这可能会影响有效的选择模式。未来的工作可能会解决这些问题，但可能需要额外的样本。
讨论
在这项工作中，我们提出了一种从时间序列遗传数据中计算一般二倍体选择系数的最大似然估计 （MLE） 的新方法。为此，我们扩展了 [16] 中的加性框架，并推导出了 EM-HMM 算法来估计二倍体选择的参数。我们表明，二倍体 EM-HMM 框架也可以通过拉格朗日乘子方法限制为定制的单参数选择模型。我们进一步引入了一种新的基于可能性的程序，用于从加性、隐性、显性选择和过度或欠显性选择之间的时间数据中推断出最合适的二倍体选择模式。据我们所知，我们的研究是第一个解决从给定时间序列遗传数据中明确确定选择模式的统计问题的研究。此外，我们实现了一种方法来估计恒定的种群大小 Ne对于给定的数据集，可以更好地对 HMM 中的遗传漂变动态进行建模。为了进一步提高检测选择和去除虚假信号的能力，我们还引入了一种基于 Brown 方法的程序，用于组合链接基因座中的 p 值。
使用模拟研究，我们表明估计的选择系数在一系列选择参数、种群参数和抽样方案中是准确的。然而，我们发现从时间序列数据中确定选择模式是具有挑战性的，并且只有在选择很强时才能产生可靠的结果。我们还证明，在分析在不同选择模式下模拟的数据时，假设加性选择会产生与在加性选择下模拟数据时相当的拒绝中性能力，这意味着分析给定数据假设加性选择可能足以扫描定向选择。但是，如果模式指定错误，则估计的选择系数将不准确。此外，我们证明，我们考虑可变总体规模的程序导致了经过良好校准的估计值和 p 值。然而，这可能与时间跨度短以及种群规模从 到 稳步增加的事实有关（参见 S1 文本中的图 O）。更极端的历史，如指数增长或严重的瓶颈，可能会更具挑战性，在这种情况下，从业者将不得不重新评估该方法。
我们将我们的方法应用于从英国 AADR [33] 中 504 个古代个体获得的时间序列遗传数据，其年代可追溯到 4,450 BP 以下，并确定了具有选择信号的 6 个基因组区域。除 TFR2 外，这些区域在以前的研究中已被确定为选择目标，我们在相关文献的背景下讨论它们。这些区域在多个方向选择模式 （加性、隐性、显性） 下被识别为显著的。然而，当从数据中对选择模式进行分类时，结果是不确定的：例如，我们发现 LCT 位点的大多数 SNP 为加性选择提供了证据，但先导 SNP 被归类为显性。此外，我们重新分析了一个时间序列数据集，该数据集由 2,400 代 146 个样本组成，这些样本来自参与马毛着色的 ASIP 基因座，并展示了在不同非加性模式下选择的证据。
请注意，我们的 HMM 实现使用 Chebychev 节点来准确计算隐藏状态空间中的单代转换矩阵，而积分考虑了在边界中吸收的概率质量。在实现模型时，准确捕获这些特征很重要 [17， 84]，因此，我们相信我们的算法可以从给定的时间数据集中识别出 MLE，并高精度地计算可能性。因此，我们在模拟研究和数据匹配模拟中的模拟研究中展示的 MLE 的统计特性和似然比检验可能不是我们的方法所独有的，而是可能在给定的群体遗传模型下表征 MLE 和各自测试的一般能力，无论我们的方法还是不同的基于可能性的方法用于分析。
此外，我们在一系列场景中表征了统计特性，但如果这些场景没有涵盖特定经验数据集中遇到的确切方案，则我们的模拟框架可以很容易地修改以表征相应场景中的统计特性。当然，识别和表征选择的能力确实取决于确切的抽样方案：即使样本仅限于短时间，也可以很容易地识别出强选择。但是，较弱的选择需要在较长的时间段内对更多数据进行采样。例如，在我们对 GB aDNA 数据集的分析中，图 9C 和 10 表明，在相应情况下检测 s = 0.01 强的选择的能力有限。
我们的方法计算一般二倍体选择参数的 MLE，我们相信这至少在两个方面对研究人员有用：（1） 我们的方法可用于从给定的时间数据集中推断选择模式。虽然我们证明选择需要相当强大才能进行可靠分类，但我们的框架可用于描述给定场景中的统计功效，并确定是否需要在潜在的额外时间点进行额外样本。（2） 如果对给定时间数据集中遗传变异进行作的选择模式是先验已知的，例如，LCT 基因座的优势或由于对复杂性状的稳定选择而导致的劣势动态，我们的方法使研究人员能够在正确的模型下准确估计选择系数。我们证明了假设错误的选择模式会产生不准确的估计。
在实践中，我们推荐以下方法来分析给定的时间序列遗传数据，可能在大量基因座上：如果计算资源有限，研究人员应应用加法 EM-HMM 来获得每个基因座的加性选择系数的 MLE，并使用标准似然比检验来识别异常值。如选择系数估计的稳健性中所述，加性选择下的似然比检验可以识别非中性重复，即使选择模式指定错误，但估计的系数不准确。作为参考，在我们的 GB aDNA 数据集中，对 125 代 504 个样本的 743,417 个 SNP 进行加法分析大约需要 5,000 个 CPU 小时。如果有额外的计算资源可用，我们建议在每种定制的单参数选择模式以及无约束模式下分析数据，以表征数据中加法选择未正确描述的信号。然后，结果也可以用于识别从数据中选择的模式，但正如我们所证明的，准确性是有限的。此外，我们强烈建议执行数据匹配的模拟研究，如数据匹配的模拟或 S1 文本中的 S.6 部分所述。这种数据匹配的模拟可以准确描述特定场景中该方法的统计功效和准确性。
对于类似的方法，我们方法的基础 HMM 假设总体是泛性的，违反此假设可能会稀释信号或引入虚假信号。解决这一缺陷的未来工作至少可以沿着两种可能的途径进行：（1）通过在估计过程中直接使用主成分作为协变量来控制种群结构[85,86]，或（2）在基础种群遗传模型中明确包括种群结构和迁移者交换[87]。
此外，我们的方法仅从焦点位点的时间数据中估计选择系数，而不纳入链接位点的等位基因频率动态。在我们对 GB aDNA 数据集的分析中，我们确实使用一种新的后处理方法利用跨基因座的信号，将 p 值组合到基因组窗口中。这种后处理可以降低基因组扫描中的信噪比，以便进行一般选择。然而，在 Wright-Fisher 过程中，使用适当的多位点动力学似然模型来整合多个 SNP 的遗传变异，有可能更准确地解释染色体连锁，并产生更稳健的推断 [88， 89]。
虽然我们在这项研究中专注于分析单个基因座的时间序列遗传数据，但在研究复杂性状的多基因选择时，我们的方法更普遍地表征选择模式的能力也具有潜在的好处：在围绕最佳性状稳定多基因选择的模型中，影响性状的遗传变异经历低于主导选择动力学， 这可以使用我们的框架轻松解决。
支持信息
补充材料
跳至 fig分享导航
很抱歉，我们无法加载您的数据。
下载
无花果分享
S1 文本。 补充材料
包含方法的其他详细信息以及补充分析的图表的文档
https://doi.org/10.1371/journal.pgen.1011769.s001
（PDF格式）
确认
我们要感谢 Xinyi Li、Xiaoheng Cheng、Constanza de la Fuente 和 Maanasa Raghavan 对该方法和数据分析的有益评论。此外，我们感谢 Jeremy Berg、Maryn Carlson 和 Rowan Hart 对手稿的评论。此外，我们还感谢 Raghavan、Berg 和 Novembre 实验室的成员在整个项目中提供的宝贵反馈。
引用
1.Bustamante CD、Fledel-Alon A、Williamson S、Nielsen R、Hubisz MT、Glanowski S 等人。人类基因组中蛋白质编码基因的自然选择。自然界。2005;437(7062):1153–7.PMID：16237444
查看文章PubMed/NCBI谷歌学术
2.Bignell GR、Greenman CD、Davies H、Butler AP、Edkins S、Andrews JM 等人。癌症基因组中突变和选择的特征。自然界。2010;463(7283):893–8.PMID：20164919
查看文章PubMed/NCBI谷歌学术
3.Nielsen R. 自然选择的分子特征。Annu Rev Genet.2005;39:197–218.PMID：16285858
查看文章PubMed/NCBI谷歌学术
4.维蒂 JJ，格罗斯曼 SR，萨贝蒂 PC。检测基因组数据中的自然选择。Annu Rev Genet.2013;47:97–120.PMID：24274750
查看文章PubMed/NCBI谷歌学术
5.Lachance J， Tishkoff SA.人类适应的群体基因组学。Annu Rev Ecol Evol Syst. 2013 年;44:123–43.PMID：25383060
查看文章PubMed/NCBI谷歌学术
6.Orlando L、Allaby R、Skoglund P、Der Sarkissian C、Stockhammer PW Avila-Arcos MC 等人。古代 DNA 分析。Nat Rev 方法引物。2021;1(1):1–26.
查看文章谷歌学术
7.Hofreiter M、Paijmans JLA、Goodchild H、Speller CF、Barlow A、Fortes GG 等人。古代 DNA 的未来：技术进步和概念转变。生物论文。2015;37(3):284–93.
查看文章谷歌学术
8.Barghi N， Tobler R， Nolte V， JakˇsiÂ'c AM， Mallard F， Otte KA.遗传冗余促进了果蝇的多基因适应。PLOS Biol. 2019;17（2）：e3000128。
查看文章谷歌学术
9.Schlötterer C， Kofler R， Versace E， Tobler R， Franssen SU.将实验进化与下一代测序相结合：研究从站立遗传变异适应的强大工具。遗传 （Edinb）。2015;114(5):431–40.PMID：25269380
查看文章PubMed/NCBI谷歌学术
10.Malaspinas A-S.使用时间序列样本表征选择性扫描的方法：古老的 DNA 视角。Mol Ecol. 2016;25(1):24–41.PMID：26613371
查看文章PubMed/NCBI谷歌学术
11.Dehasque M、Ávila-Arcos MC、Díez-Del-Molino D、Fumagalli M、Guschanski K、Lorenzen ED 等人。从古代 DNA 推断自然选择。Evol Lett. 2020;4(2):94–108.PMID：32313686
查看文章PubMed/NCBI谷歌学术
12.Bollback JP、York TL、Nielsen R. 从时间等位基因频率数据估计 2Nes。遗传学。2008;179(1):497–502.PMID：18493066
查看文章PubMed/NCBI谷歌学术
13.Malaspinas A-S， Malaspinas O， Evans SN， Slatkin M. 从时间序列数据估计等位基因年龄和选择系数。遗传学。2012;192(2):599–607.PMID：22851647
查看文章PubMed/NCBI谷歌学术
14.王 J.一种伪似然方法，用于从时间间隔样本中估计有效总体大小。Genet Res. 2001 年;78(3):243–57.PMID：11865714
查看文章PubMed/NCBI谷歌学术
15.Ferrer-Admetlla A， Leuenberger C， Jensen JD， Wegmann D.Wright-Fisher 扩散的近似马尔可夫模型及其在时间序列数据中的应用。遗传学。2016;203(2):831–46.PMID：27038112
查看文章PubMed/NCBI谷歌学术
16.Mathieson I， McVean G. 从等位基因频率的时间序列数据估计空间结构总体中的选择系数。遗传学。2013;193(3):973–84.PMID：23307902
查看文章PubMed/NCBI谷歌学术
17.Tataru P， Simonsen M， Bataillon T， Hobolth A. 使用等位基因频率数据在 Wright-Fisher 模型中进行统计推断。系统生物学 2017;66（1）：e30-46。PMID：28173553
查看文章PubMed/NCBI谷歌学术
18.Vlachos C、Burny C、Pelizzola M、Borges R、Futschik A、Kofler R 等人。用于检测和定量进化和重测序研究中选择的基准测试软件工具。基因组生物学 2019;20(1):169.PMID：31416462
查看文章PubMed/NCBI谷歌学术
19.新泽西州 Gemmell J. 石板 J. 杂合子的繁殖力优势。公共科学图书馆一号。2006;1（1）：e125.PMID：17205129
查看文章PubMed/NCBI谷歌学术
20.海德里克 PW。杂合子优势选择的证据是什么？.Trends Ecol Evol.2012;27(12):698–704.PMID：22975220
查看文章PubMed/NCBI谷歌学术
21.Palmer DS， 周 W， Abbott L， Wigdor EM， Baya N， Churchhouse C， et al.英国生物样本库中的遗传显性分析。科学。2023;379(6639):1341–8.PMID：36996212
查看文章PubMed/NCBI谷歌学术
22.Sanjak JS， Sidorenko J， Robinson MR， Thornton KR， Visscher PM. 当代人类定向和稳定选择的证据。美国国家科学院院刊 2018 年;115(1):151–6.PMID：29255044
查看文章PubMed/NCBI谷歌学术
23.巴顿 N.通过突变和稳定选择之间的平衡来维持多基因变异。Genet Res. 1986 年;47(3):209–16.
查看文章谷歌学术
24.de Vladar HP， Barton N. 多基因性状的稳定性和对稳定选择和突变的响应。遗传学。2014;197(2):749–67.PMID：24709633
查看文章PubMed/NCBI谷歌学术
25.西蒙斯 YB、Bullaughey K、Hudson RR、Sella G。人类数量性状的 GWAS 发现的群体遗传解释。PLoS 生物学 2018;16（3）：e2002985。PMID：29547617
查看文章PubMed/NCBI谷歌学术
26.Koch E、Connally NJ、Baya N、Reeve MP、Daly M、Neale B 等人。遗传关联数据与塑造人类常见疾病和性状的稳定选择大体一致。冷泉港实验室;2024. https://doi.org/10.1101/2024.06.19.599789
27.Cheng X， Steinrücken M. diplo-locus：用于在一般二倍体选择下推理和模拟时间序列遗传数据的轻量级工具包。生物Rxiv。2025;2023.10.12.562101.PMID：37905072
查看文章PubMed/NCBI谷歌学术
28.Steinrücken M， Bhaskar A， Song YS.一种从时间序列遗传数据中推断一般二倍体选择的新型光谱方法。Ann Appl Stat. 2014 年;8(4):2203–22.PMID：25598858
查看文章PubMed/NCBI谷歌学术
29.Foll M， Shim H， Jensen JD.WFABC：一种基于 Wright-Fisher ABC 的方法，用于从时间采样数据中推断有效种群规模和选择系数。Mol Ecol 资源。2015;15(1):87–98.
查看文章谷歌学术
30.Schraiber JG、Evans SN、Slatkin M. 从等位基因频率时间序列中对自然选择的贝叶斯推断。遗传学。2016;203(1):493–511.PMID：27010022
查看文章PubMed/NCBI谷歌学术
31.Iranmehr A， Akbari A， Schlötterer C， Bafna V. Clear：进化和重序列实验的可能性组成。遗传学。2017;206(2):1011–23.PMID：28396506
查看文章PubMed/NCBI谷歌学术
32.Taus T， Futschik A， Schlötterer C. 使用 pool-seq 时间序列数据量化选择。分子生物学 Evol.2017;34(11):3023–34.
查看文章谷歌学术
33.Mallick S、Micco A、Mah M、Ringbauer H、Lazaridis I、Olalde I 等人。Allen Ancient DNA Resource （AADR） 是一本精心策划的古代人类基因组纲要。科学数据。2024;11(1):182.PMID：38341426
查看文章PubMed/NCBI谷歌学术
34.Ludwig A、Pruvost M、Reissmann M、Benecke N、Brockmann GA、Castanos P 等人。马驯化初期的毛色变化。科学。2009;324(5926):485.
查看文章谷歌学术
35.EWENS WJ。数学群体遗传学。2 ed. Springer;2004.
36.Watterson GA. 在单个基因座测试选择。生物测定学。1982;38(2):323–31.PMID：7115865
查看文章PubMed/NCBI谷歌学术
37.拉宾纳 LR.有关语音识别中隐藏马尔可夫模型和选定应用的教程。程序 IEEE。1989;77(2):257–86.
查看文章谷歌学术
38.模式识别和机器学习。斯普林格;2006.
39.霍夫曼 LD、布拉德利 GL、罗森 KH。商业、经济学、社会和生命科学的应用微积分。扩展第 10 版麦格劳-希尔;2010.
40.Mathews J， Fink K. 使用 matlab 的数值方法。第 4 版，皮尔逊;2003.
41.Varadhan R， Roland C. 用于加速任何 EM 算法收敛的简单全局收敛方法。Scand J Stat. 2008 年;35(2):335–53.
查看文章谷歌学术
42.威尔克斯 SS。用于检验复合假设的似然比的大样本分布。Ann Math Statist.1938;9(1):60–2.
查看文章谷歌学术
43.Sawyer SA，哈特尔 DL。多态性和分歧的群体遗传学。遗传学。1992;132(4):1161–76.PMID：1459433
查看文章PubMed/NCBI谷歌学术
44.Browning SR， Browning BL. 根据血统身份片段对近期有效人口规模的准确非参数估计。Am J Hum Genet.2015;97(3):404–18.PMID：26299365
查看文章PubMed/NCBI谷歌学术
45.Mathieson I， Terhorst J. 青铜时代英国自然选择的直接检测。基因组研究 2022;32(11–12):2057–67.PMID：36316157
查看文章PubMed/NCBI谷歌学术
46.Olalde I、Brace S、Allentoft ME、Armit I、Kristiansen K、Booth T 等人。Beaker 现象和西北欧的基因组转化。自然界。2018;555(7695):190–6.PMID：29466337
查看文章PubMed/NCBI谷歌学术
47.钦塔拉帕蒂 M， 帕特森 N， 穆尔贾尼 P.欧洲全新世期间主要人类混合事件的时空模式。eLife。2022;11：e77625.
查看文章谷歌学术
48.Patterson N、Isakov M、Booth T、Büster L、Fischer C-E、Olalde I 等人。青铜时代中晚期向英国的大规模移民。自然界。2022;601(7894):588–94.PMID：34937049
查看文章PubMed/NCBI谷歌学术
49.Gretzinger J、Sayer D、Justeau P、Altena E、Pala M、Dulias K 等人。盎格鲁-撒克逊人的迁移和早期英国基因库的形成。自然界。2022;610(7930):112–9.PMID：36131019
查看文章PubMed/NCBI谷歌学术
50.Moorjani P、Sankararaman S、Fu Q、Przeworski M、Patterson N、Reich D。一种用于确定古代基因组年代的遗传方法提供了过去 45,000 年人类世代间隔的直接估计。美国国家科学院院刊 2016 年;113(20):5652–7.PMID：27140627
查看文章PubMed/NCBI谷歌学术
51.Foll M、Poh Y-P、Renzette N、Ferrer-Admetlla A、Bank C、Shim H 等人。流感病毒耐药性：时间采样群体遗传学视角。PLoS 基因。2014;10（2）：e1004185。PMID：24586206
查看文章PubMed/NCBI谷歌学术
52.Benjamini Y， Hochberg Y. 控制错误发现率：一种实用而强大的多重测试方法。J Roy Statist Soc Ser B： 国家主义方法。1995;57(1):289–300.
查看文章谷歌学术
53.史密斯 JM，海格 J.有利基因的搭便车效应。Genet Res. 1974 年;23(1):23–35.
查看文章谷歌学术
54.布朗 MB。一种组合非独立、单侧显著性检验的方法。生物测定学。1975;31(4):987–92.
查看文章谷歌学术
55.Nait Saada J、Kalantzis G、Shyr D、Cooper F、Robinson M、Gusev A 等人。对 487,409 个英国样本的血统鉴定检测揭示了精细的种群结构和极其罕见的变异关联。Nat Commun.2020;11(1):6130.PMID：33257650
查看文章PubMed/NCBI谷歌学术
56.Field Y、Boyle EA、Telis N、Gao Z、Gaulton KJ、Golan D 等人。过去 2000 年人类适应的检测。科学。2016;354(6313):760–4.PMID：27738015
查看文章PubMed/NCBI谷歌学术
57.Bersaglieri T、Sabeti PC、Patterson N、Vanderploeg T、Schaffner SF、Drake JA 等人。乳糖酶基因近期强阳性选择的遗传特征。Am J Hum Genet.2004;74(6):1111–20.PMID：15114531
查看文章PubMed/NCBI谷歌学术
58.Itan Y， Powell A， Beaumont MA， Burger J， Thomas MG.欧洲乳糖酶持久性的起源。PLoS 计算生物学 2009;5（8）：e1000491。PMID：19714206
查看文章PubMed/NCBI谷歌学术
59.Mathieson I、Lazaridis I、Rohland N、Mallick S、Patterson N、Roodenberg SA 等人。230 名古代欧亚人的全基因组选择模式。自然界。2015;528(7583):499–503.PMID：26595274
查看文章PubMed/NCBI谷歌学术
60.Peter BM， Huerta-Sanchez E， Nielsen R. 区分选择性扫描与站立变异和从头突变。PLoS 基因。2012;8（10）：e1003011。PMID：23071458
查看文章PubMed/NCBI谷歌学术
61.Enattah NS， Sahi T， Savilahti E， Terwilliger JD， Peltonen L， Järvelä I. 鉴定与成人型低乳酸症相关的变异。Nat Genet.2002;30(2):233–7.PMID：11788828
查看文章PubMed/NCBI谷歌学术
62.Burger J、Link V、Blocher J、Schulz A、Sell C、Pochon Z 等人。青铜时代欧洲乳糖酶持久性的低患病率表明过去 3,000 年来持续的强烈选择。Curr Biol. 2020;30(21):4307–15.
查看文章谷歌学术
63.Evershed RP， Davey Smith G， Roffet-Salque M， Timpson A， Diekmann Y， Lyon MS. 欧洲乳品、疾病和乳糖酶持久性的演变。自然界。2022;608:336–45.
查看文章谷歌学术
64.普渡大学 MP、Lan Q、Wang SS、Kricker A、Menashe I、Zheng T-Z 等人。对 Toll 样受体基因变异和非霍奇金淋巴瘤风险的汇总调查。发生。2009;30(2):275–81.PMID：19029192
查看文章PubMed/NCBI谷歌学术
65.Sun J、Wiklund F、Hsu FC、Bälter K、Zheng SL、Johansson JE 等人。白细胞介素 1 受体相关激酶 4 和 toll 样受体 6-1-10 基因簇中序列变异的相互作用增加了前列腺癌的风险。癌症流行病生物标志物 Prev. 2006;15(3):480–5.PMID：16537705
查看文章PubMed/NCBI谷歌学术
66.马 X， 刘 Y， 高文 BB， Graviss EA， 克拉克 AG， 马瑟 JM.全外显子重测序显示 toll 样受体变异导致人类对结核病的易感性。公共科学图书馆一号。2007;2（12）：e1318。PMID：18091991
查看文章PubMed/NCBI谷歌学术
67.Wong SH、Gochhait S、Malhotra D、Pettersson FH、Teo YY、Khor CC 等人。麻风病和人类 toll 样受体 1 的适应。PLoS 病原体。2010;6（7）：e1000979。PMID：20617178
查看文章PubMed/NCBI谷歌学术
68.Barreiro LB、Ben-Ali M、Quach H、Laval G、Patin E、Pickrell JK 等人。人类 Toll 样受体的进化动力学及其对宿主防御的不同贡献。PLoS 基因。2009;5（7）：e1000562。PMID：19609346
查看文章PubMed/NCBI谷歌学术
69.Hedrick PW， Thomson G. HLA 平衡选择的证据。遗传学。1983;104(3):449–56.PMID：6884768
查看文章PubMed/NCBI谷歌学术
70.布朗森 PG、马克 SJ、埃利希 HA、斯拉特金 M。基于序列的方法表明，经典人类白细胞抗原 （HLA） 基因座的平衡选择是不对称的。Hum Mol Genet.2013;22(2):252–61.
查看文章谷歌学术
71.Lao O， de Gruijter JM， van Duijn K， Navarro A， Kayser M. 从单核苷酸多态性分析中揭示的与人类皮肤色素沉着相关的基因中正选择的特征。安·胡姆·热内特。2007;71:354–69.
查看文章谷歌学术
72.Beleza S、Santos AM、McEvoy B、Alves I、Martinho C、Cameron E 等人。欧洲人色素沉着变亮的时间。分子生物学 Evol.2013;30(1):24–35.PMID：22923467
查看文章PubMed/NCBI谷歌学术
73.Soejima M， Koda Y. 色素沉着相关基因 SLC24A5 和 SLC45A2 中两个编码 SNP 的种群差异。国际法律医学杂志，2007 年;121(1):36–9.PMID：16847698
查看文章PubMed/NCBI谷歌学术
74.Hysi PG、Valdes AM、Liu F、Furlotte NA、Evans DM、Bataille V 等。对欧洲血统的个体进行全基因组关联荟萃分析，确定了解释很大一部分头发颜色变异和遗传力的新基因座。Nat Genet.2018;50(5):652–6.PMID：29662168
查看文章PubMed/NCBI谷歌学术
75.Fisher SA、Hampe J、Onnie CM、Daly MJ、Curley C、Purcell S 等人。复杂疾病中的直接或间接关联：SLC22A4 和 SLC22A5 功能变异在克罗恩病中的作用。Hum Mutat. 2006 年;27(8):778–85.PMID：16835882
查看文章PubMed/NCBI谷歌学术
76.Huff CD、Witherspoon DJ、Zhang Y、Gatenbee C、Denson LA、Kugathasan S. Crohn 病和 IBD5 的遗传搭便车。分子生物学 Evol.2012;29(1):101–11.
查看文章谷歌学术
77.Girelli D、Bozzini C、Roetto A、Alberti F、Daraio F、Colombari R 等人。由于转铁蛋白受体 2 基因的新突变引起的 3 型血色病的临床和病理发现。胃肠。2002;122(5):1295–302.PMID：11984516
查看文章PubMed/NCBI谷歌学术
78.Rhodes SL、Buchanan DD、Ahmed I、Taylor KD、Loriot MA、Sinsheimer JS 等。帕金森病中铁相关基因的汇总分析：与转铁蛋白的关联。神经生物学 Dis. 2014;62:172–8.PMID：24121126
查看文章PubMed/NCBI谷歌学术
79.Dutcher SK， Brody SL. HY-DIN' 在纤毛中：在原发性纤毛运动障碍中发现中央对相关突变。Am J Respir Cell Mol Biol. 2020;62(3):281–2.PMID：31604022
查看文章PubMed/NCBI谷歌学术
80.吞咽 DM。乳糖酶持久性和乳糖不耐症的遗传学。Annu Rev Genet.2003;37:197–219.PMID：14616060
查看文章PubMed/NCBI谷歌学术
81.He Z， Dai X， Lyu W， Beaumont M， Yu F. 从古代 DNA 数据中估计时间可变的选择强度。分子生物学 Evol.2023;40（3）：msad008。PMID：36661852
查看文章PubMed/NCBI谷歌学术
82.Rieder S， Taourit S， Mariat D， Langlois B， Guérin G. 刺豚鼠 （ASIP）、延伸 （MC1R） 和棕色 （TYRP1） 基因座的突变及其与马 （Equus caballus） 毛色表型的关联。哺乳动物基因组。2001;12(6):450–5.PMID：11353392
查看文章PubMed/NCBI谷歌学术
83.Wutke S、Benecke N、Sandoval-Castellanos E、Döhle H-J、Friederich S、Gonzalez J 等人。马的斑点表型在中世纪失去了吸引力。科学代表 2016;6：38548。PMID：27924839
查看文章PubMed/NCBI谷歌学术
84.Paris C， Servin B， Boitard S. 使用对 Wright-Fisher 模型的各种参数近似从遗传时间序列中推断选择。G3 （贝塞斯达）。2019;9(12):4073–86.PMID：31597676
查看文章PubMed/NCBI谷歌学术
85.Luu K， Bazin E， Blum MGB.pcadapt：一个 R 包，用于执行基因组扫描，以便根据主成分分析进行选择。Mol Ecol 资源。2017;17(1):67–77.
查看文章谷歌学术
86.Ju D，马蒂森 I。欧亚大陆西部皮肤色素沉着相关变异的演变。美国国家科学院院刊 2021 年;118（1）：e2009227118。PMID：33443182
查看文章PubMed/NCBI谷歌学术
87.Joseph TA， Pe'er I. 从时间序列基因型数据推断种群结构。Am J Hum Genet.2019;105(2):317–33.PMID：31256878
查看文章PubMed/NCBI谷歌学术
88.Terhorst J， Schlötterer C， Song YS.来自实验进化的基因组时间序列数据的多位点分析。PLoS 基因。2015;11（4）：e1005069。PMID：25849855
查看文章PubMed/NCBI谷歌学术