遗留数据扰乱了基因组学研究

摘要

最近的报告已经确定了不同人群突变谱的差异。尽管这些报告中的一些已经在其他队列中复制，但大多数仅在1000基因组项目（1kGP）数据中报告。在调查日本人口中一个有趣的假定人口分层时，我们发现了一个之前未报告的批处理效应，导致1kGP数据中出现虚假突变呼叫，并导致明显的人口分层。由于1kGP数据被广泛使用，我们发现批量效应还导致了主要插补服务器和少数可疑GWAS关联的错误插补。因此，1kGP早期阶段的低质量数据继续以隐藏的方式污染现代研究。现在可能是退休或升级此类遗留测序数据的时候了。

介绍

老化参考队列数据中的批处理效应

最后5个 近年来，人类基因组序列数据的数量和质量急剧增加。参考群体，如国际HapMap项目(2005年国际HapMap联盟)1000基因组项目（1kGP）(2010年1000基因组项目联盟,2012;联合体等。2015)和西蒙斯基因组多样性项目（SGDP）(Mallick等人。2016)例如，已经公开了数千个基因组序列，用于人口和医学遗传分析。通过提供插补服务的服务器，可以间接获得更多的基因组(McCarthy等人。2016)或用于变量频率估计的汇总统计(Lek等人。2016)。

对1kGP中的第一个基因组进行了测序10 几年前(van Dijk等人。2014). 自那时以来，测序平台得到了迅速改进。与早期相比，1kGP的第二阶段实施了多项技术和分析改进(1000基因组项目联合会2012;联合体等。2015)导致在项目过程中进行异质样本制备和数据质量。

然而，由于免费获得的数据具有非凡的价值，来自1kGP的早期数据仍被广泛用于插补非类型变体、估计等位基因频率以及回答广泛的医学和进化问题。这就提出了一个问题，即是否以及如何将这些遗留数据与最近的队列一起纳入当代分析。Mafessoni等人。(2018)最近通过寻找远处变异中LD过量的个体，确定了1kGP中的批量效应。在这里，我们指出了1kGP早期阶段的这些和其他未报告的批效应是如何通过群体遗传分析和虚假GWAS关联导致错误的遗传结论的，这是使用1kGPs作为参考进行插补的结果。

突变特征

不同的诱变过程可能优先影响不同的DNA基序。例如，烟草烟雾中的某些诱变剂已被证明优先与某些基因组基序结合，导致过量的G到T颠换(Pfeifer等人。2002;Pleasance等人。2010). 因此，通过考虑多态性的DNA背景来寻找签名不同突变过程(Alexandrov等人。2013;Shiraishi等人。2015). 这种全基因组突变特征已被用作癌症的诊断工具(Alexandrov等人。2013;Shiraishi等人。2015)。

除了体细胞突变特征外，最近人们对种系突变特征的群体变异也产生了兴趣，这种变异特征可以在大型测序小组中发现。2015年，哈里斯报告TCC增加50%→ 与非洲人群相比，欧洲人群中的TTC突变，这在2017年的不同队列中重复出现(哈里斯2015;哈里斯和普里查德2017;Mathieson and Reich 2017年). 突变特征的强烈群体富集表明每个群体的历史中存在重要的遗传或环境差异(哈里斯2015;哈里斯和普里查德2017). Harris和Pritchard进一步确定了一系列种群的不同突变谱，Aikens等人在最近的一份出版物中对此进行了进一步研究(哈里斯和普里查德2017;Aikens等人。2019)。

另一个签名，*AC→ *与1kGP和SGDP中的非洲人相比，东亚人观察到CC的频率更高(哈里斯和普里查德2017;Aikens等人。2019). 这两项研究还发现，在1kGP的日本个体中，这种特征的频率是异质的。这种异质性令人感兴趣，因为种系特征的差异在许多代人中积累。日本人口内的系统差异表明，日本各亚群体之间存在持续的环境或遗传差异，几乎没有基因流动。由于日本样本数量较少，该观察结果无法在SGDP中重现(Aikens等人。2019). 因此，我们决定使用来自长滨的日本个体的新序列数据集来跟踪这一观察结果。虽然我们无法在日本人群中再现突变异质性，但我们可以追溯到1kGP数据中的技术伪影。除了在突变签名中产生偏见外，这种人为因素还导致虚假插补结果，这些结果在最近的出版物和在线资源中都有发现。

结果部分组织如下。我们首先尝试再现原始信号，并从1kGP中识别JPT队列中的有问题的变体。接下来，我们将分析扩展到1kGP中的其他人群，并确定显示技术偏见证据的变体列表。最后，我们研究了这些变异是如何影响现代基因组学分析的。

结果

日本特有的突变特征

哈里斯和普里查德（2017）报告了过量的三元取代模式*AC→*在1kGP的部分日本个体中进行CC。我们最初的目标是确定这种特征是否可以用人口结构或技术错误来解释。当我们试图在一个来自长滨的更大和更近的日本队列中对这一观察结果进行跟踪时，我们没有发现这种特殊的特征。当比较1kGP和这个更大的数据集中的日本个体的等位基因频率时，我们观察到两组中的一组个体存在一些单核苷酸多态性（SNP）(图1). 鉴于这两个种群的相似性，这强烈表明存在技术差异，而不是种群结构效应。尽管在这两个队列中只考虑满足严格质量掩码和Hardy–Weinberg平衡的站点，但这些不匹配仍然存在。

当从1kGP数据中删除不匹配站点时，*AC→*CC信号消失(图1). 为了确定差异的可能技术原因，我们对*AC的流行率进行了回归分析→*针对1kGP提供的不同个体级质量指标的CC突变签名（参见补充图S1,补充材料在线）。映射基的平均质量问每个个体都表现出强烈的相关性：问显示签名速率提高。因此，从低位调用的序列-问数据中包含的变异在研究中重现性较差，并显示出特定的突变特征。

为了确定在没有进入第二个队列的情况下，可能在队列中繁殖不良的SNP，我们在JPT中进行了一项与低水平问(图1A). 传统上，全基因组关联研究使用基因型作为自变量。在这里，我们进行了基因型条件关联测试（GCAT），其中基因型现在是我们使用连续变量预测的因变量问作为自变量(Song等人。2015). 我们使用基因型的logistic回归问并识别587个SNPP（P） < 10⁻⁸和1034个SNPP（P） < 10⁻⁶。在确定要排除的假定低质量SNPs时，使用更高的P（P）值阈值增加了过滤的严格性（即，排除具有以下特征的SNPP（P） < 10⁻⁶比排除SNPS更严格P（P） < 10⁻⁸). 与关联的变量问在*AC中有显著富集→*CC突变，GA*→GG*和GC*→GG*突变(图1C). 这三种富集可以概括为G**的过剩→低GG*患者问使用χ计算这些富集的统计显著性²测试跟踪哈里斯和普里查德（2017）.

因此，这种突变信号在问-相关SNP，但在非显著SNP中仍存在残余信号，可能是因为在低水平个体中发现了许多罕见的等位基因问由于缺乏动力，使用关联技术仍无法识别(补充图S2,补充材料在线）。将个人与问<30成功删除*AC→*CC信号在日本富集，但哈里斯和普里查德发现的其他信号似乎没有变化，包括*AC的大陆富集→*Harris和Pritchard报告的东亚人与非洲人的CC信号比较，并在SGDP中复制(补充图S3和S4,补充材料在线）。对于对稀有变异积累敏感的群体遗传分析问似乎优于过滤特定的低质量SNP。对于插补质量重要的其他分析，确定问-相关变体可能更可取。

识别1000基因组项目中的可疑变异

分布问在1kGP种群中，许多种群的分布问得分与JPT相当，尤其是在项目第一阶段测序的人群：在1kGP早期进行的测序变化较大，总体上倾向于包括质量较低的测序数据(图2和补充图S5,补充材料在线）。这种变异性可能是由不断发展的测序平台和协议或测序中心之间的变异造成的。到2011年，旧的测序技术被淘汰，方法变得更加一致，从而产生更高、更统一的质量。

因此，我们在所有人群中独立地执行了相同的反向GWAS方法，并进行了类似的鉴定问-在1kGP的26个种群中，23个种群中存在相关SNP，其中第1阶段种群受影响最大，与第3阶段种群相比，显著相关位点的数量平均为4倍。超过1165个变异与低问具有P（P） < 10⁻⁶在每个(补充图S6,补充材料在线）。

为了建立一个测试统计量来同时表示所有人群的相关性，我们进行了一个简单的logistic回归，根据问用logistic factor analysis（LFA）作为人口结构或基因型条件关联测试（GCAT）的抵消(Song等人。2015). 我们还考虑了两种解释混淆因素的替代方法，即使用前五个主成分和使用人口成员身份作为协变量。这些模型大体一致（参见补充图S7,补充材料在线）。

该方法识别了与问分布在基因组中，15270个通过1kGP严格屏蔽过滤器（补充图S8–S11，补充材料在线）。以下大多数分析侧重于满足严格掩码的15270个变体，因为这些变体不太可能被标准管道过滤。为了说明大量测试的原因，我们使用了两阶段Benjamini和Hochberg升压FDR控制程序来调整P（P）使用标称I类错误率的值α = 0.01 (Benjamini等人。2006). 我们分别测试了SNP、INDEL和重复区域，因为它们可能有不同的错误率(表1). 的列表问-相关变异和个体问在中提供补充材料在线。

问-相关变异体分布在整个基因组中，1号染色体显示出过多的此类变异体，而其他染色体则相对统一(补充图S12A类,补充材料在线）。染色体1在*AC中表现出强烈的富集→*CC信号与其他染色体相比，尽管每个染色体检测的变异数正常化(补充图S12B类,补充材料在线）。在a 1 Mb标度，我们还发现分布相当均匀，少数区域显示此类变体富集(补充图S12C类,补充材料在线）。三个外围1 染色体1、2和17中的Mb区域超过30个问-相关变体。这些地区的协会统计数据分布情况见补充图S13,补充材料在线。相比之下，未通过1kGP严格屏蔽的变异在基因组中分布更不均匀(补充图S12D类,补充材料在线）。

GCAT模型中此变体列表的突变3聚体替换模式与1kGP JPT单种群测试中确定的特征相似，因为*AC→*还有从TAC、TCT和TGT到同源核苷酸TTT的突变富集(补充图S14,补充材料在线）。

最近的出版物Mafessoni等人。(2018)使用一种使用连锁不平衡而非质量指标的方法来识别19196个可疑变异体，其中67%通过了1kGP严格屏蔽，从而在1kGP-中发现了批量效应(补充图S15A类,补充材料在线）。他们确定了17917个与异常LD模式显著相关的变异，这些变异与问。我们发现，1279（3%）个他们发现的变异体也在我们的可疑变异体列表中，并且具有相关性P（P）使用GCAT方法识别的值。我们还发现23111问-与异常LD模式无关的相关变体。有趣的是，我的Mafessoni等人鉴定的变体在前面描述的突变光谱中没有得到丰富(补充图S15B类,补充材料在线）。这些结果表明，在1kGP中可能存在大量的批量效应，而这些效应只能通过一系列关联测试来识别。

假变异体、偏倚基因型或细胞系伪制品？

评估是否问-相关变异是由测序伪影或表现出偏向性基因分型的真实变异引起的假变异，我们将原始的1kGP序列数据与最近使用相同细胞系的测序工作进行比较。问-在测序实验中不重复的相关变体很可能是1kGP中测序工件的结果。相比之下，问-在实验中重复出现的相关变异可能是由（反复发生的）细胞系突变引起的，或者更可能是由基因分型依赖于问最后，考虑到我们的名义错误发现率α = 0.01，我们预计为问-关联变体是错误发现，即与问因此，这应该在实验中重现。

2017年，Lan等人。(2017)从1kGP中重新测序83名汉族人。在296个这样的变体中问-在CHB或CHS的单种群测试中，重新排序的数据中只有6个(补充图S16,补充材料在线）。这稍微多了一点，但与根据α = 0.01标称错误发现率。因此，大多数问-CHB或CHS中的相关变异似乎是虚假变异。

我们使用GCAT模型中确定的所有变体（而不仅仅是与问CHB和CHS内部）。15270年问-在全球范围内发现的相关变异体中，有6307个在83个重测序个体的1kGP中具有多态性（参见补充图S17,补充材料在线）。绝大多数多态性与问在重新排序数据集中根本不存在，这表明它们是伪变量。五种变异在两个数据集之间显示出不同的频率，可能是由有偏见的基因型解释的。最后，1139（或18%）出现在可比频率的重测序数据中。由于这高于错误发现率，我们得出结论，队列GCAT测试主要识别假变异体，但也识别在CHB/CHD以外的人群中具有偏见基因型的变异体。

最后，在15270人中问-Illumina的Omni 2.5芯片上存在相关变体613（参见补充图S18,补充材料在线）。这些很可能是数据中出现的少数变异，但在1kGP中表现出有偏见的基因分型。

可疑变异影响现代基因组学分析

最先进的插补服务器使用许多数据库的组合，包括一些无法免费获得的数据库。从研究人员的角度来看，它们充当黑盒插补机，将观察到的基因型作为输入并返回插补基因型。

为了调查来自1kGP的可疑呼叫是否被输入到基因分型研究中，我们向密歇根插补服务器提交了1kGP基因型数据的前两条染色体的基因型数据。我们发现，所有与问在样本中重新插补。这表明插补参考小组仍然包括低水平的个体问而这些可疑的变异将在与低质量个体最接近的个体中进行插补。这些问-相关变异也可能影响邻近真实变异的插补，然而，当考虑来自Nagahama的日本个体基因型数据的插补分数时，似乎对邻近变异插补分数没有任何影响(补充图S19,补充材料在线）。

我们在文献中搜索了可能报告这些情况的GWAS问-相关变异体与某些生物性状显著相关，尽管这些变异体没有特定的原因与表型相关。NHGRI-EBI公布的全基因组关联研究目录确定了17份最近的出版物，其中报告这些变异接近或高于全基因组显著阈值(表2)。

其中11项研究将1kGP纳入其插补参考小组(Xu等人。2012;Lutz等人。2015;Park等人。2015;Astle等人。2016;Herold等人。2016;López-Mejías等人。2017;Nagy等人。2017;Spracklen等人。2017;Suhre等人。2017;Tian等人。2017;Gao等人。2018)另一个直接使用1kGP序列数据和细胞系(Mandage等人。2017). 一项研究使用内部参考小组进行插补(Nishida等人。2018)，两项研究对个体进行了基因分型，并使用HapMap II作为插补参考数据库对数据进行了插补(Kraja等人。2011;Ebejer等人。2013)，两项研究使用基因分型芯片数据(Yucesoy等人。2015;Ellinghaus等人。2016)。

所有这些文章都使用了各种严格的质量过滤器，包括Hardy–Weinberg平衡测试、预期等位基因频率偏差和测序数据质量阈值。他们还删除了罕见等位基因和高度缺失的等位基因。事实上，我们预计会有大量问-通过哈迪-温伯格平衡测试等质量控制筛选出相关变体。尽管这些研究使用了最先进的质量控制，但这些变异被输入到基因型数据中，并在全基因组范围内与生物性状相关。然而，事实上，其他研究中的一些变体没有被删除，并且这些变体中的绝大多数在高质量数据集中缺失，这意味着问-相关变体应标记为删除，以避免虚假关联。

这些关联并不一定是错误的——插补中的微弱但显著的偏差可能仍然会导致正确的关联。区分与问从具有强烈偏见的变异中，我们区分了低、高等位基因频率差异的变异-问大于2的因数（这自然会将两组变体分开补充图S20,补充材料在线）。大多数（92.7%）问-相关变异具有强烈的偏见，因为它们在低水平个体中的发生率是低水平个体的两倍以上-问与高-问数据。相比之下，大多数问-GWAS目录中报告的相关变体具有微弱的偏差（参见补充图S21,补充材料在线），但有三个例外。一项研究报告了七种高度偏见的关联问-将1kGP中爱泼斯坦-巴尔病毒（EBV）序列的拷贝数作为表型的相关变异(Mandage等人。2017). 这七个变体没有出现LD峰值，并已从下游分析中正确删除，但它们仍包含在NHGRI-EBI目录中(Mandage等人。2017). EBV拷贝数表型可能对绘图敏感，并被导致偏倚SNP调用的相同技术伪影所混淆。

讨论

本研究中确定的变体可以通过遗留技术的技术工件进行解释。不同的测序技术会有不同的错误模式。一份比较基因组分析仪II（GAII）和Illumina HiSeq的报告发现，GAII在质量分数低于30分时的读取率要高得多(米诺什等人。2011)例如，不同的质量模式会随着读取而降低。读取质量和错误配置文件的差异反过来需要不同的调用管道。

一条染色体突变谱的富集(补充图12B,补充材料在线）很难通过生物机制或测序技术进行解释，因为预计两者都不会在染色体上产生系统性偏差。由于染色体通常是分开处理和分析的，所以这些数据在整合或处理方式上的细微差异可以解释这种偏差。然而，要查明差异的准确技术来源，还需要进一步法医调查异质样品制备、测序技术组合和整个1kGP使用的数据处理管道的细节。鉴于自1kGP早期以来在测序和调用方面取得的进展(图2)最近的序列数据中可能不再积极引入这些偏差的来源。

然而，由于1kGP数据被广泛用作参考数据库，这些变体仍被输入到新的基因型数据中，从而影响各种表型的关联研究。尽管变体与质量指标的显著关联本身并不表明该变体是虚假的，但我们建议仔细检查GWAS与此类变体的关联，例如，通过重复分析而不将1kGP作为插补小组的一部分。

对于无法单独检测个体变异的分析（突变谱、等位基因频率分布、多基因风险评分），我们建议只需将具有问 < 30和筛选出由Mafessoni等人。(2018)和问-我们识别的相关变体（此类变体和样本ID的列表在补充材料在线）。我们还建议插补服务器丢弃低水平的个体问（或至少提供不进行插补的选项）。考虑到自由访问数据的价值问对社区来说也是一项有价值的投资。

结论

在技术方面，我们感到惊讶的是，在1kGP项目中，变异体和技术协变量之间的强大关联之前尚未被确定。基因型在质量指标上的全基因组logistic回归分析是简单易行的，可能是各种组学研究的标准。LFA的计算要求更高，但产生了更稳健的结果(Song等人。2015). 我们测试的所有三种人口结构方法都产生了可比较的结果。

更一般地说，为了提高基因组参考数据集的质量，我们可以通过添加新的更好的数据和更好地管理现有数据来进行。鉴于技术进步迅速，基因组研究的重点自然是数据生成方面。然而，清理现有数据库对于避免生成虚假结果也很重要。目前的研究结果表明，1kGP项目最终发布的数据中有相当一部分已经过期，无法退休或重新排序。

材料和方法

代码和数据可用性

由于这些数据主要使用公开可用的数据执行，我们在GitHub上提供了完全可复制和公开可用的脚本(https://github.com/LukeAndersonTrocme/LegacyData网站，上次访问时间为2019年9月15日）。该存储库包含用于数据下载、处理、分析和绘图的脚本。

元数据

此分析中使用的元数据是从1kGP文件系统中的每个索引文件编译而来的。映射基础的平均质量问每个样本都是从与每个比对文件关联的BAS文件中获得的。每个BAS文件都有关于每个样本的每个测序事件的元数据。如果一个样本被多次测序，我们取每个样本的平均值问每个测序实例的得分。序列索引文件中提供了分析中每个样本的提交日期和测序中心。

质量控制

对于突变谱分析，我们复制了Harris等人使用的质量控制和数据过滤管道，因为他们应用最新的质量阈值来删除可疑序列，以检测人群水平差异。应用了几个掩模文件来删除基因组中可能质量较低的区域，或者与基因组的其余部分相比，可能具有非常不同的突变率或碱基对复杂度的区域。使用1kGP严格掩码删除基因组的低质量区域，使用phastCons100way掩码文件删除高度保守区域，使用RepeatMasker的NestedRepeats掩码文件删除高重复区域。此外，只考虑缺失度<0.01、MAF<0.1和MAF>0.9的位点。总的来说，7786023个双列常染色体变体通过了突变谱分析的质量控制。我们使用类似的方法计算了JPT群体显著变异列表的碱基对三联体的突变谱，如(哈里斯和普里查德2017)。

在GCAT测试中，唯一使用的过滤是应用0.000599的次要和主要等位基因频率截止值（去除单倍体、双倍体和三倍体），结果总计S公司 = 测试中包括28516063个变体。我们还使用NestedRepeats掩码文件标记重复区域内的变量，因为这些变量是单独分析的，用于错误发现率估计。由1kGP严格掩码标记的变量包含在关联测试中，并包含在FDR调整中。这些变异只有在FDR之后才被删除，并且被排除在错误模式的下游讨论之外，因为大多数群体遗传学分析都使用严格的掩码作为过滤器，我们希望在过滤区域中发现有问题的变异。

检测质量与基因型的相关性

在零假设下，斜率系数β_秒为零，模型2减为模型1。β_秒表示与映射基的平均质量的关联问_我到基因型年_硅为了检验零假设，我们使用广义似然比检验统计量，其偏差是相加边际重要性的度量问_我在模型中。零模型下的偏差检验统计量近似为一个自由度的叉方分布。

我们总共运行了S公司回归，其中S公司是基因组位点的总数。考虑到大量的测试、大量的预期零假设以及基因组之间的正相关性，我们使用两阶段Benjamini和Hochberg的加速FDR控制程序来调整P（P）值(Benjamini等人。2006). 通过使用标称I类错误率α = 0.01，共发现15270个变异具有统计学意义。请参见补充材料在线查看变体列表并进行调整P（P）值。

个体特异性等位基因频率

推测

使用密歇根插补服务器，我们从1kGP中插补了染色体1和2的基因型数据。我们使用了来自1kGP Omni 2.5的基因型数据 M芯片基因型数据。然后下载从服务器返回的VCF文件，并用于搜索成功插补的重要变异数。

致谢

我们要感谢K.Harris分享她的突变谱脚本，F.Mafessoni、G.Santpere和A.Navarro对文章早期版本的反馈。我们还要感谢Gravel实验室的成员在编码和有用讨论方面的帮助。

工具书类