背景

随着基因分型平台纳入更多标记，基因分型成本不断下降，越来越大、越来越复杂的数据集正在被分析。因此，用于分析基因型数据集的计算效率高的非参数方法越来越多地被用于揭示种群结构。种群结构的解析揭示了群体间的进化关系。此外，在全基因组关联研究中必须考虑种群结构，以减少由病例和对照之间的祖先差异导致的虚假关联[1].

主成分分析（PCA）是一种广泛用于人口结构分析的非参数方法，它使用协方差矩阵进行特征分析。个体之间的变异量和变异轴分别捕获在特征值和特征向量中。此前，我们开发了一个用于人口结构分析的PCA框架，该框架扩展了PCA的使用范围，通过使用迭代过程简化人口结构模式，使其超出了通常用于可视化人口结构趋势的应用范围。其他人使用的迭代方法，例如[2,3]与我们的客观方法不同，主观地将个人分组依赖于可用的民族地理人口标签。

我们的框架，我们称之为迭代修剪PCA（ipPCA），使用聚类算法将个体分配到子种群中，而不强加任何先验假设[4]. ipPCA解析人口数据集中的所有子种群，从而报告原始子种群的总数K（K）除了分配其中包含的个人之外。术语“人口”与ipPCA的数据集是同义词，ipPCA是可用于分析的整个个体集合。术语“亚群”定义了ipPCA分配的一组个体，其中不存在其他重要的子结构。ipPCA通过使用基于投影数据点和簇质心之间的欧几里德距离的聚类算法，将个体系统地分为两个簇。分离个体的决定需要测试数据集（或算法后续迭代的嵌套数据集）中是否存在重要结构。为了测试各组个体之间的同质性，我们之前建议使用在EIGENSTRAT/SmartPCA算法中实现的测试统计量，该算法根据Tracy-Widom（TW）分布报告结构的概率[5]. 如果不存在显著的结构，则被测个体属于一个子种群，从而终止了迭代聚类过程。图中总结了ipPCA框架1使用模拟数据和实际数据的数据集，我们展示了ipPCA如何正确地将个体分配到子群体并推断K（K）然而，ipPCA的准确性可能会受到停止准则的影响。不适当的终止标准会导致对亚种群数量的低估或高估。此外，早期迭代中的个别赋值错误将被合成并结转至后期迭代。

图1
图1

ipPCA框架概述该框架由三个主要部分组成。首先，对遗传数据进行编码，以零位为中心并进行标准化。然后，将个体投影到由输入数据矩阵的主成分所跨越的空间上。然后，计算结构度量以决定是前进到聚类步骤还是终止算法。当度量未超过阈值时，将解析同质子种群，然后算法终止。否则，这些个体被一分为二。该算法迭代，直到所有个体都被分配到终端子种群中。

全尺寸图像

将个体聚类为子群体的参数算法，例如STRUCTURE、，弗拉佩、ADMIXTURE和BAPS在一个关键方面与ipPCA不同，即将个体分配到子种群簇的方法。上述参数算法分别推断每个个体的祖先比例，并将具有类似推断祖先模式的个体分组。ipPCA和其他非参数方法无法推断祖先。这些技术试图将具有相似遗传特征的个体组合在一起。因此，参数方法仍然提供了非参数分析无法看到的重要信息。然而，大型和高度结构化的人口数据集对于参数分析来说是难以处理的，因为K（K）祖先集群是有限的。这是因为用于估计亚群等位基因频率的可用样本数量有限。为了更好地观察固有的种群结构，应该对重新取样的个体进行“监督”结构分析。这种监督分析的个人选择是任意的，通常以可用的民族/地理标签为指导。尽管如此，需要仔细选择，以确保被比较的个体具有相似的祖先，否则对区分某些群体的个体来说重要的祖先信号可能太弱。

在本文中，我们提出了对ipPCA的一种改进，为迭代聚类过程引入了一种新的停止准则，称为EigenDev，它对大数据集中的伪模式更具鲁棒性。这种新算法被称为EigenDev ipPCA。为了区分ipPCA框架中的两种算法，我们在后面的章节中将之前提出的使用TW统计量作为终止准则的算法称为TW-ipPCA。此外，我们还提出了一种新的协议，该协议使用来自特征Dev-ipPCA的信息来指导参数分析。使用真实数据集，我们演示了该方法如何揭示无监督结构分析无法检测到的新的和结构信息型祖先模式。

方法

新的ipPCA终止准则

Tracy-Widom（TW）测试统计，在EIGENSTRAT/SmartPCA算法中实现[5]，用作TW-ipPCA算法的停止准则。虽然这种停止标准在一些数据集上很有效，但我们发现，当分析包含大约1000个以上个体的更大数据集时，TW-ipPCA解决的亚群远远超过预期。因此，我们怀疑，在某些情况下，当采样量较大时，所解析的子群体可能是虚假的，即i型错误。的确，正如[5]，潜在亚群的相对样本量影响TW检验统计。

除了我们在使用TW结构统计测试时发现的I型错误外，还有其他缺陷促使我们开发一种替代的终止标准。第一个问题是计算困难。为了获得TW检验统计量的最终值，需要估计太多的未知参数。这些参数没有最佳估计值，因此估计值的选择会影响结果。而不是使用第页-TW检验统计量的值作为阈值，我们提出了一个新的终止准则来确定数据是否结构化。新准则基于数据矩阵的特征值，称为特征开发启发式。特征值开发启发式遵循与TW理论相同的假设，即如果数据矩阵的第一个特征值显著大于其余特征值，则存在子结构。然而，我们扩展了这一观察，而不仅仅是测试第一个特征值的重要性，以考虑数据的剩余方差。这使我们能够观察更高维度的结构。我们受到启发，根据信号处理领域中应用的特征值Grads启发式算法开发了特征值开发工具[6]. 这项工作表明，如果数据只包含噪声而不包含信号，即非结构化数据，则按降序排列的特征值具有良好的线性拟合。在种群遗传数据中，噪声代表一个（亚）种群内的自然遗传变异。

为了测试种群结构，根据基因型数据计算特征Dev统计量。此计算首先要求由编码的、零位的和标准化的基因型数据构建数据矩阵，如[5]. 该矩阵包含对应于个体的行和对应于等位基因的列。因此，双等位基因SNP标记由两列中的条目编码，每列一个等位基因，STR由数据集中该标记位点的等位基因总数编码。等位基因的存在编码为1，其缺失编码为0。对于缺失的数据，即没有基因型调用的标记，它们被编码为全部0。

给定零位，标准化基因型数据矩阵X（X）（根据[5])包含米样品n个每个样本的等位基因列，我们构建样本协方差矩阵

然后可以根据以下公式计算EigenDev值

(1)

哪里

(2)

和

（3）

哪里，是第一个第页的特征值C类按降序排列。等式中的数量(1)在某些情况下可能是负面的。为了防止这种可能性，将编码条目归一化为平均值为零。这一步对于从公共元素中删除信号很重要，只留下个体之间的差异（遗传方差）用于特征分析。在对模拟数据和实际数据进行的所有实证研究中，我们发现数据中90%的方差总是导致正值，并且所讨论的凸性约束从未被违反。为了说明遇到负值的罕见情况，我们在算法中包含了一个检查步骤，以检测和报告负值。如果发现负值，则参数第页可以调整以确保平方根中的正数。回想一下第页<最小值{米,n个}是用于计算EigenDev统计的特征值数量。我们还使用对数变换来稳定方差。如果特征值较大，则被分析的个体组将包含多个亚群，ipPCA会将该组一分为二；否则，当EigenDev值低于阈值时，EigenDev-ipPCA算法终止。

结果

测试

为了测试EigenDev概念，分析了几个数据集：

  1. 1

    一个由来自同一群体的10000个个体组成的模拟数据集，每个个体包含10000个SNP标记，用于测试TW分布的拟合度。它是使用GENOME工具生成的[7]使用以下参数和以下树文件：

-pop 1 10000-N树.txt-C 20-S 500

树.txt：

0 10000

1-1

1 10000

GENOME从10000个创始人群体开始，产生了与创始人规模相同的第一代人。每个个体有20条染色体，每条染色体包含500个SNP。

  1. 2

    使用与第一个数据集相同的GENOME参数模拟第二个数据集，但使用不同的树文件：

树.txt：

0 5000 5000

1-1 2-1

40 5000

1-1 1-2

80 5000 5000

1-1 2-1

100万

生成两个5000个个体的亚群。

  1. 三。

    第三个数据集是Bovine HapMap项目收集的497个个体，来自19个不同品种，27203个SNP的基因型。可从以下网址公开获取：[8].

  2. 4

    第四个数据集公开于[9]. 它包含3945个个体，包括185个不同的种族/地理标签，为1327个标记（包括848个微卫星、476个indels和3个SNP）进行了分类，这些标记来自[10].

模拟和真实复杂数据集的ipPCA编码输入矩阵也可从下载http://www4a.biotec.or.th/GI/tools/ippca.

人口结构的测试指标

为了测试TW是如何受到采样的影响的，在10到200个个体的20个不同样本量下随机采样了一个没有子结构的模拟数据集。用于测试TW分布拟合度的相应概率-概率（p-p）图如图所示2.可以观察到，大多数样本大小都违反了TW分布；仅对70名受试者进行了良好拟合。因此，与TW分布的偏差将给出错误检测（I型误差），特别是对于大样本量。另一方面，由于TW测试基于非线性相位变化，因此对于检测结构非常敏感。只要有足够的数据可用，它就不容易出现II型错误[5]. 然而，相位变化的非线性意味着存在一种全或无的情况，即类型I的可能性无法控制，即使在很大的p值阈值范围内也是如此。

图2
图2

测试TW分布的拟合度使用合并模型模拟了一个具有5000个标记的10000个个体的种群。针对10至200个个体的样本量生成了p-p图。

全尺寸图像

为了测试EigenDev启发式算法的性能，我们模拟了第二个模拟数据集中100、200和500个人的三种不同样本大小的接收机工作特性（ROC）曲线，如图所示3为了获得曲线，特征值阈值在0.077和0.387之间变化。观察到阈值随着样本量的增加而增加，当样本量较大时，EigenDev的性能更好。实际数据集的分析使用了0.21的EigenDev阈值。该值是三种样本大小达到10%假阳性率所需阈值的平均值。该值是检测和解析所有结构之间的良好折衷，在实际数据集的典型样本大小下，伪结构最小。

图3
图3

EigenDev启发式算法的经验接收机工作特性曲线使用合并模型模拟了包含两个亚群（每个亚群5000个）的5000个标记的10000个个体的结构化群体。ROC是针对100、200和500个人的样本量生成的。

全尺寸图像

用ipPCA指导参数分析

结构[11]可以使用祖先组件信息执行无监督聚类。然而，STRUCTURE的计算复杂度很高，尤其是在求K（K）祖先集群实际上限制为K（K）=20或更少。因此，必须将高度复杂的数据集划分为子数据集，然后由STRUCTURE单独分析这些子数据集。按照惯例，这是使用先前的信息，例如种族-地理人口标签，以任意方式完成的。然而，先验信息可能会使聚类结果产生偏差。为了解决这个问题，我们建议使用ipPCA的无监督聚类功能，以更有效的方式帮助缩小STRUCTURE的搜索空间。实际上，可以选择ipPCA分配的子种群进行后续的结构分析。我们称这种方法为ipPCA-guided STRUCTURE。我们将此方法应用于牛HapMap数据集[8]，这是我们之前分析的扩展数据集[4]. 结果与之前报告的结果类似，即EigenDev-ipPCA解析了18个亚群，每个亚群主要由同一品种的个体组成，但一个亚群包含安格斯（ANG）和红安格斯（RGU）个体。

STRUCTURE使用默认参数和10000次烧入和10000次运行迭代。从Gir（Gir）、Brahman（BRM）和Nelore（NEL）品种中选择个体进行结构分析，以确定这些品种之间是否存在推断祖先的差异。此外，选择这三个品种是因为它们印度双歧杆菌品种，因此彼此之间的关系比彼此更密切B.金牛在数据集中繁殖。结构分析K（K）=3，如图所示4，揭示了以前没有报道过的具有繁殖特征的祖先模式。

图4
图4

基于ipPCA的牛HapMap数据集个体结构分析。对来自瘤牛品种（GIR、BRM和NEL）。结果与K（K）=2和K（K）=3如图所示。

全尺寸图像

基于ipPCA的大型人体数据集分析

来自Tishko等人的数据集[10]包含大量个体（3945）。此外，这些个体包含185个种族-语言区分标签，表明存在大量不同的遗传群体。该数据集由EigenDev-ipPCA分析，该分析分配了49个子群（图5). 分配的亚群与之前报告的模式基本一致[10]在非洲，地理上不同的个体群体在基因上是不同的，而在非洲，主要的文化和语言群体也在基因上不同（参见附加文件1更多信息）。相反，ipPCA使用TW停止标准（TW-ipPCA）分配了109个子种群。比较两种方法不同的亚群表明，总体上，TW-ipPCA分配的亚群是EigenDev-ipPCA分配给较大亚群的子簇。例如，所有印度人（15个种族标签）通过EigenDev ipPCA被分配到两个亚群（SP2和SP7），而印度人通过TW ipPCA被分配到11个亚群（见附加文件1).

图5
图5

使用EigenDev-ipPCA方法对Tishkoff等人数据集的人口分配.49分配的子种群标记为SP1到SP49。条形图的高度与每个亚群中指定个体的数量成比例。指定个体的种群标签显示在每个栏的右侧，括号中有相同标签的个体数量。为了帮助可视化个人分配，185个人口标签被分为14个颜色组，以反映地理区域。颜色组中的颜色渐变表示不同的人口标签。有关完整的配色方案，请参见Additional文件中的图s31.

全尺寸图像

ipPCA-指导结构分析

非裔美国人是一个用来描述具有自我认同的非洲血统的美国国民的术语，其中大多数是通过奴隶贸易来到美国的西非个人的后裔。不过，“非裔美国人”这一术语非常广泛，因为它涵盖了来自广泛地理范围的非洲祖先的后代，有些人也有最近的非非洲血统。EigenDev-ipPCA将非裔美国人分为四个亚群，即SP4、SP5、SP15和SP16。SP4和SP5亚群包括大多数非洲裔美国人，以及主要讲西非和中非尼日尔-霍尔多芬语的少数民族。五名非裔美国人被分配到SP15，其中主要包括来自苏丹的讲非洲库什语的贝扬人。两名非裔美国人被分配到SP16，SP16主要由尼罗·萨哈兰（Nilo-Saharan）苏丹语和非洲库西语混合种族的东非人组成。

然后，我们使用来自EigenDev-ipPCA的信息来指导STRUCTURE。分配到SP4、SP5、SP15和SP16的所有个体，包括所有非洲裔美国人，都通过STRUCTURE从K（K）=2至K（K）=5（参见图中的A部分6). 在K（K）=3或更大，EigenDev-ipPCA分配的四个亚群中的每一个都显示出不同的祖先模式，尽管SP15和SP16个体之间似乎有一些重叠。当关注非裔美国人时，也可以观察到不同的血统模式，特别是在比较SP4和SP5指定个体时（参见图中的B部分6).

图6
图6

ipPCA-对Tishkoff等数据集中选定个体的结构分析.A）分配给SP4、SP5、SP15和SP16的所有个体（见图5），包括所有非洲裔美国人，均通过STRUCTURE从K（K）=2至K（K）= 5. 根据ipPCA分配对个人进行分类。图中还显示了每个亚群体中个体的主要民族语言标签（完整列表见图5）。B） 从A）扩展了对非洲裔美国人的看法。

全尺寸图像

讨论

TW和EigenDev停止标准

群体遗传结构分析首先需要一种方法来检测数据集（或嵌套数据集，用于ipPCA的进一步迭代）中是否存在显著结构。获取此信息的当前方法是测试从主成分分析计算出的最大特征值的Tracy-Widom分布是否存在偏差。A类第页-值小于10-12被认为是一个可以接受的阈值，用于拒绝数据属于同质（子）总体的零假设，从而构建数据[5]. 对无结构模拟数据集的第一次实验表明，与预期分布存在显著偏差，尤其是在大样本（>70个人）情况下。我们从这个结果推断，当样本量较大时，由于TW分布的这种偏差，TW方法会出现I型误差。只需使用较低的第页-值阈值可能不会给出更好的结果，因为第页-实用的价值[5]. 当应用于实际数据集时，高密度采样的同质（子）种群可能被错误地解释为具有结构。在TW-ipPCA中，这将导致一组个体被分配到单独的子种群中，而实际上这些个体应该被视为属于单个（子）种群。

为了缓解TW检验统计量的缺点，我们提出了一种称为EigenDev统计量的新终止准则，该准则计算更简单，没有隐藏参数，并且对I型误差更具鲁棒性。为了简单起见，可以选择单个特征值作为ipPCA的通用停止标准，这需要通过经验确定。我们通过数据模拟确定了0.21的阈值，这也适用于本文分析的实际数据集。

牛HapMap数据集分析

EigenDev-ipPCA的亚群分配支持公认的观点，即牛品种具有独特的遗传特征。ANG和RGU被分配在同一个亚群中的发现表明，这些品种在可用的标记上在遗传上无法区分，这也被其他方法报道[12]. 然而，发现GIR、BRM和NEL品种被EigenDev-ipPCA分解为单独的亚群是新的，因为早期的非监督结构分析[12]在整个数据集上无法区分这些品种。基于ipPCA的牛HapMap数据集结构分析表明，这些品种的祖先存在差异，这与EigenDev-ipPCA的分配一致。在这些indicine品种中，有证据（高杂合度和独特SNP）表明BRM与其他品种（包括GIR和NEL）在遗传上不同[12]. 这些结果提出了一个问题，为什么结构分析在特征提取-ipPCA指导下进行时，可以揭示这些品种之间的差异，而这些差异在无监督结构分析中并不明显？最可能的解释是，与其他品种相比，这些indicine品种中信息标记的总数较低（只有19%的位点的次要等位基因频率大于0.3）[12]. 换句话说，与牛磺酸品种相比，indicine品种的等位基因频率高度相关。与其他组相比，具有高度相关等位基因频率的个体组倾向于被STRUCTURE合并[11].

对大型人类数据集的分析

EigenDev-ipPCA分配的49个亚群中，每个亚群都包含基本上共享相同种族语言标签/从属关系的个体，根据[10,13]. 值得注意的是，EigenDev-ipPCA将426名印度个体分为两个亚群。该分组与这些个体在[13]显示出微弱的结构证据。因此，与EigenDev ipPCA相比，TW ipPCA解决的更大程度的分层可能是虚假的。因此，TW-ipPCA解析出的伪结构归因于大样本量（426），远高于模拟数据分析中I类错误的阈值。

在非洲个体中，EigenDev-ipPCA分配了亚群，揭示了之前未描述的分层模式。例如，来自西非和中非的非俾格米人奈杰尔·霍尔多凡语个体在基因上无法区分[10]，但被EigenDev-ipPCA分配到SP3、SP4和SP5亚群。EigenDev-ipPCA将大多数非裔美国人分配给SP4和SP5（图5)表明他们有西非和中非尼日尔-霍尔多凡人的祖先，这与[10]. 另一方面，EigenDev-ipPCA将非裔美国人分配到不同的亚群，这表明这些人之间存在显著的结构。在中执行的监督结构运行[10]阐明非裔美国人的血统只能揭示非裔美国人之间微妙的变异模式。然而，EigenDev-ipPCA引导的结构分析显示，SP4和SP5非裔美国人的血统存在明显差异。SP15和SP16分配的非裔美国人也显示出不同于SP4和SP5分配的个体的血统，尽管分配给SP15和SP 16的个体数量较少，但不可能观察到这两个组之间的显著血统差异。

一些非裔美国人被EigenDev-ipPCA分配到SP15和SP16是出乎意料的。这些亚群中的当代非洲个体主要来自撒哈拉和东非。最近一项对非裔美国人血统的研究得出结论，一些人的主要祖先成分既不是西非尼日尔-霍尔多凡人，也不是欧洲人[14]. 这种异常祖先可能是撒哈拉人或东非人，这也可能反映在线粒体DNA单倍型上，因为一些非裔美国人具有未知非洲起源的异常单倍型[15,16]. 特征Dev-ipPCA引导的结构与在中执行的监督结构之间的差异[10]是由于分析中个人的选择。当使用与其他人的等位基因频率差异不适当的个体时，将忽略关键的祖先差异，这与牛数据分析中显示的相同。

结论

我们描述了用于分析种群结构的特征Dev-ipPCA。这种方法将个体分配到子种群，并确定存在的子种群总数。该算法采用了一种新的启发式算法，称为EigenDev，用于检测子结构，并将其应用于迭代聚类过程。EigenDev对人口抽样具有鲁棒性，允许我们以更高的精度分析大型复杂数据集。由EigenDev-ipPCA分配的子种群揭示了个体群体之间的总体遗传关系，然后可以使用这些遗传关系来指导结构。其他参数算法，如混合和弗拉佩也可以以相同的方式使用。因此，特征Dev-ipPCA和STRUCTURE的结合是互补的，可以一起用于执行强大的人口分层分析。Matlab源代码（m-file）和Windows和Linux（64位）上的可执行版本的软件都可以从以下网站下载：http://www4a.biotec.or.th/GI/tools/ippca.