一种基于大规模并行全基因组测序的单核苷酸多态性检测方法

检测方法

技术领域

本发明涉及一种基于新一代大规模并行全基因组测序的单核苷酸多态性的检测方法，属基因工程技术领域。

背景技术

遗传突变影响着生物体表型的变异，罹患疾病的风险，对药物和环境刺激的反应。全基因组连锁分析和定点克隆技术对受单基因影响的孟德尔遗传病的研究已经取得了巨大的成功。然而，绝大多数常见疾病，如糖尿病，心血管疾病和癌症等重要的数量性状具有复杂的遗传学基础，由多个基因以及基因与环境因素之间的相互作用共同决定。连锁分析在检测对疾病的影响近于中性的遗传变异时具有极大的局限性。相对于连锁分析和定点克隆来说，全基因组关联性分析能够为疾病关联基因位点定位提供一个更好的方法。人类基因组序列的测序完成，使得数以百万计的单核苷酸多态性位点得以鉴定，构建了高密度的单倍体型基因组变异图谱。这些研究进展开创了使用大规模全基因组的单核苷酸多态性检测技术来寻找引起各种人类疾病或与之相关的基因变异的时代。二十多年以来，Sanger法测序和荧光电泳技术一直在DNA测序领域占据着主导地位，使用随机鸟法策略或目标区域的PCR扩增法进行了大量的SNP检测。现有dbSNP数据库中的大多数SNP位点都是通过这些方法鉴定的。鸟法测序进行SNP检测的标准方法是将测序片段与基因组参考序列进行比对，并根据碱基的质量打分来过滤掉低质量的测序错误，得到较为可信的SNP结果。使用来自二倍体样本的PCR扩增序列进行直接测序，再通过对色谱图进行分析，检测出杂合的多态性位点，也是常见的方法，主要软件有SNPdetector，novoSNP，PolyPhred以及PolyScan。

与传统的毛细管电泳法测序相比，新一代测序技术如Illumina Genome Analyzer (GA), AB SOLiD，以及Roche 4 FLX系统显着地提高了测序通量，极大地降低了成本。Illumina GA一次运行可以产生约四千万条长度为50bp的测序片段。这种超高的测序通量使得新一代的测序技术特别适合于在已知参考基因组序列的基础上进行大规模个体的重测序从而进行基因变异的研究。截至当前，使用新的测序技术已经完成了两个人的基因组测序：James Watson的个人基因组测序（Roche 4FLX）和第一个亚洲人的基因组测序（Illumina GA）。此外，国际千人基因组计划执行委员会也决定使用这种测序技术对来自全世界的1000个个体基因组进行测序，得到最详细的人类基因组变异图谱。

随着新的测序技术的发展，相应的SNP检测方法也有了很好的发展；然而，由于新测序技术产生的测序片段与以往相比有显著的差异，新的精确的SNP检测方法也亟待开发。

发明内容

本发明的目的是提出了一种适应于大规模平行测序法Illumina GA技术特点的构建待测基因组一致序列和检测SNP的新方法。在分析了Illumina GA测序数据的错误特征以后，解决了由于测序过程和数据处理过程中实际问题所造成的一致序列不准确及SNP不准确等问题，从而为采用新一代测序技术进行高效、快速、准确的全基因组测序分析提供了可靠的手段。

本发明提出的一致序列构建方法和SNP检出方法，包括以下步骤： （1）将测序数据比对到参考基因组上。

使用短序列比对程序SOAP(http://soap.genomics.org.cn) 将Illumina GA的测序片段比对到参考基因组序列上。

（2）对所有唯一比对上的测序片段进行统计

通过对SOAP程序结果的文本处理，统计特定测序质量值和特定测序序列坐标下，每两种碱基之间的错配比例，将此比例作为对错配概率的估计，记录在一个四维概率矩阵里面，作为统计学模型中各项参数的基础。

（3）判别基因组上每个碱基的基因型

对于基因组上每一个位点，将比对在此位点上所有的测序片段碱基收集起来，记录其碱基类型、测序质量和在测序片段上的序列坐标，从四维概率矩阵中查出四种碱基观察到测序碱基的概率。对于二倍体基因组而言，其真实基因型的可能性共有10种（纯合基因型4种：AA、CC、GG、TT；杂合基因型6种：AC、AG、AT、CG、CT、GT）。从每一种真实基因型观察到覆盖该位点的所有碱基的概率，为观察到每一个单独碱基的概率之积，而后者是可以从步骤（2）中建立的概率矩阵中查到的。这样，我们就得到了每一种潜在可能的基因型得到此位点的观察碱基的似然概率（likelihood）。考虑参考基因组的碱基类型和已知多态性位点的信息，我们可以为每一个潜在的基因型赋予一个先验概率，与likelihood结合，得到后验概率。后验概率最高的基因型则为此基因组位点最有可能正确的基因型，即一致序列基因型。该基因型正确的概率，为其后验概率在所有10个基因型后验概率之和中所占的比例。我们将对基因型估计的正确率转换为一个质量分数。

（4）鉴定单核甘酸多态性位点

将步骤（3）所构建的待测基因组的一致序列中与参考基因组序列不一致的位点从结果中抽出，作为潜在的多态性位点。通过对质量分数设定阈值，将错误率降低至1%以下。

附图说明

图1. 大规模平行测序仪一致序列构建的方法框架 图2. 测序片段比对的准确性和唯一性

图3. 测序质量分数的不准确和错配的偏向性。 图4. 校正前后的起点数目分布。

图 5. 在不同质量阈值下全基因组覆盖度，基因分型位点的覆盖度和错误率。 图 6. 在不同测序深度下(a) 基因分型位点覆盖度 (b) 欠识别和过识别的几率。

具体实施方式

（一）计算最优基因型

我们用贝叶斯公式依据给定的等位基因类型和在染色体上的对应位置的质量值来推断基因型，图1描述了这个方法的具体实现步骤。本方法所需要的输入数据是由Illumina测序技术给出的测序片段(reads)。该方法使用的所有的reads都是都是在参考基因组序列唯一比对上的（例如，将一个亚洲个体基因组序列数据以NCBI人类基因组作为参考序列所做的比对结果），并用这些序列构建待测基因组的一致序列。首先，将序列的质量分数，也就是每个碱基测序错误率的估计值，通过唯一比对上reads进行统计，然后校正，使之更加准确。然后计算基因组每个位置上各种基因型的似然值，然后通过贝叶斯公式计算每种基因型的后验概率，从而确定最有可能的基因型，最后用秩和检验进一步消除错误的杂合位点。

1.基因型的先验概率估计。给定一个可用的参考基因组，可以估计待分析的基因组序列相对于参考基因组的突变率。例如，人类基因组上绝大部分位置每代发生突变的概率在10-8数量级，而两个个体和其共同祖先的差异大约要经过10000代的时间形成。据此推算出来的两个人的单倍染色体SNP概率大约为0.001。 每个单倍体跟参考序列大概有1000个核甘酸的差异。假设人的参考基因组序列有0.00001的错误率，则此错误率相对于真实多态性位点的发生率可

以忽略。综合这些数据，对于双倍体我们设发生纯合SNP的概率为0.0005，杂合子为0.001。

根据以往对NCBI dbSNPs的研究，转换突变发生的频率是颠换的4倍，但是在这两种变异中各类型之间却几乎是等频率出现的。由此，在我们的SNP检测模型中用到了这些比率。例如，假如参考序列上某位点基因型为G，单倍体碱基类型可能是A，C，和T的先验概率都是6.67E-4，G的概率是0.999，而对于双倍体出现GG组合的概率是0.9985，AA是3.33E-4，TT是8.33E-5，AC和AT是1.11E-7，GC和GT是1.67E-4，AG是6.67E-4，CT是2.78E-8（表1）。 表 1. 二倍体基因组的先验概率。假定参考序列基因组碱基型为脱氧鸟甘酸（G）

A C G T A C G T 3.33E-4 1.11E-7 6.67E-4 1.11E-7 8.33E-5 1.67E-4 2.78E-8 0.9985 1.67E-4 8.33E-5

2.用质量分数进行似然估计。每个位置上的候选等位基因型D可以从与参考序列完全比对上的reads观察出来。假设每种碱基类型Ti的可能性为

P(D|Ti)， 影响 P(D|Ti)主要因素有四个，包括a)碱基类型，b)测序质量，c)read上的位置，d)出现的次数。所有这些因素都被我们考虑在了本方法的似然估计模型中。

测序过程中的错误和短片段的比对错误都可能引起待测序列与参考序列之间差异的出现。既然下一代的测序片断长度很短，那么会有更多的reads出现比对错误，因此就出现了错误的SNP。通过设置reads支持最低数目阈值能够滤掉绝大部分的错配的碱基型。测序中，reads靠近3’端低质量的区域，错误通常都不是随机的：AC和ＧT错配在统计上显著地（ｐ<0.0001）高于其它错配。为解决这个问题，我们把ｒｅａｄ坐标,质量分数以及错误的偏向性均考虑在内，用数组校正质量分数。然后用这些经过校正的质量分数进行基因型的似然估计。为了避免错误之间的相关性，本方法对ＰＣＲ扩增所得的大量拷贝进行了质量罚分，以减少其对基因型估计的影响。

（二）测序片段比对的准确性

为了评价read覆盖的准确度和唯一性，我们以人类参考基因组１２号染色体作为参考序列，模拟了不同长度的测序片段（包含了0.001比例的ＳＮＰ和测序错误），然后将这些reads重新比对到人类全基因组参考序列中去。当reads定位到参考序列上具有最少碱基错配的位置时，我们认为是最佳匹配。如果一个read在参考序列上只有一个最佳匹配位置，则称此匹配称为唯一匹

配。有多个位置是最佳匹配（同一个read在参考序列上能够和多个位置以同样数目的错配数的比对上）时，我们认为它是重复比对。

从这些数据中我们可以计算出具有唯一比对reads的百分比。对于单向测序唯一比对上的reads，其唯一比对的比例在测序长度从15bp增长到25bp时变化较大，之后再增加测序长度时，唯一比对的比例仅有很少的变化。78.6%的25bp长度的reads和91.5%的50bp的reads能够唯一地比对到基因组上。Illumina测序技术典型的测序长度是35bp，在比对结果中发现，没有错误的reads有85.7%能够唯一地比对到基因组上，1个错配的reads有86.3%，两个错配的reads有85.9%。模拟数据的比对结果与炎黄一号基因组实测比对结果相似。从分析中，我们发现，用双向测序技术能够很大的提高唯一比对的

reads比例。双向测序技术插入片段长度在100bp至10kbp（±10%）范围内,能够唯一比对的reads比例随着插入片段长度的增加而增加；其中，在插入片段为200bp的时候有95.4%的成对reads能够唯一的比对上。

用短序列测序片段去比对的话，reads里面包含的SNP和一些测序错误都有可能导致reads比对到不正确的位置上。对于模拟的reads，由于原始的位置事先知道，所以我们能够估计不能比对正确的reads的含量。长度为25bp的单向测序reads，2.3%的具有1个测序错误的reads和3.5%的具有2个测序错误的reads没能正确的匹配上（图2C）。如果测序长度为50bp，这两个比例分别降到了0.6%和0.8%。同理模拟了含有一个测序错误的双向测序reads，插入片段从100bp-10kbp，分别有0.4%和0.06%的reads比对错误；具有2个测序错误的reads在同样的插入片段下有0.3%和0.06%的不能比上。参见（图2D）

我们从模拟数据中计算了碱基的错误率以检测错误识别的SNPs。至少95%的错误碱基在单向测序reads（35bp）和双向测序reads（插入片段200bp）中都只出现了一次。基于这些数据，我们设置了一个阈值滤除所有低频率的错误碱基（通常设为4）。用这个方法在炎黄一号基因组上查找SNPs，结果只有大约0.036%的错误等位基因没有被检出。在随机DNA片段的模拟中，36倍测序数据大下约有0.008%的测序错误未被滤除。为了区分具有高频碱基型的错误和真实的杂合多态性位点，我们利用二项分布（P=0.0001）来检测这两种碱基型的差别，发现剩余的错误碱基型中87.3%都被滤除了。总共，99.93%的由于reads测序错误导致的错误碱基型都被滤除了。

除此之外还有一种导致SNP错误的因素是reads中含有插入删除，这种错误跟reads的比对方式有关。由于SNPs数量大概是小片段的插入删除的5~10倍，比对中优先考虑了不加gap的比对情况，这类含有插入删除的reads可能比对失败。我们允许最长不超过3bp的插入或删除，通过这种方式比对出确实含有插入删除的reads。我们模拟了10000个小片段插入删除来评价SNPs检测的潜在影响，结果发现0.6%的包含插入删除的reads没有检测出来。如果我们要求至少4个reads支持才可信的话，只有3（0.03%）个错误SNP等位基因被检出。

（三）Illumina测序质量分数的校正

测序错误的累积导致测序片段3’端错误率要远高于5’端。Illumina测序也能给出一个质量打分，但是这个质量分数是通过信号强度来计算的，并不准确代表错误率的发生。为了纠正这个问题，我们对炎黄一号基因组测序的结果进行了评估，设计了一种方法根据比对错配率来校准质量分数。在这些reads中，真正的SNPs(发生率0.001)和测序错误(发生率0.01)混在一起。在校正中为了尽可能的避免SNPs造成的干扰，我们把已知的SNP位点排除在外。

Illumina测序技术通过每个测序循环的数据来校准质量分数。经过校准后的质量分数跟实际的错配率有一定的差异，而且这个差异是随着read坐标的不同上下波动的(图3A)。目前我们对质量分数的校正仍然是通过比对信息和原始序列信息或者经过Illumina校正后的分数得来。

除了测序误差随着测序循环的不断增加对质量分数造成影响之外，测序仪器本身对碱基的检测也影响了质量分数。Illumina测序技术利用两种不同频率的激光照射四种被标记过的碱基，A，C用一种激光表示，G，T用另一种激光表示。所以AC, GT测序错误出现的频率就高于其它的错配。我们发现AC, GT错配的概率要比我们模拟比对的情况分别高出58%和72%，同时CG错配概率大约要降低36%（图3B）。例如，质量分数为10（理论错误率0.1）的碱基中，reads与参考序列之间能够观察到的错配AC, CA, GT, TG分别为4.62%，5.27%，5.29%，4.62%，而其它的错配类型仅有1.62%~2.48%。鉴于此结果，我们也按错配类型校正了质量分数。

（四）重复序列的罚分策略

由于DNA建立文库过程中使用了PCR扩增，给定很少数量的DNA起始片断，能够产生大量的拷贝，因此我们能够获取大量一致长度的ＤＮＡ片断。然而，这些大量复制的片断对测序过程中的随机性有很大影响，有些区域的测序深度不是很理想。这也可能导致各种易与杂合等位基因位点混淆的错误的出现。特别是ＤＮＡ的损伤被ＰＣＲ扩增以后可能带来的冗余序列覆盖度，这些可重复的错误很难和多态性位点区分开。因此我们设置了针对扩增重复的罚分规则。如果DNA文库和测序过程都是随机的，那么序列的起点的分布也应该是服从泊松分布。炎黄一号基因组测序的深度是３６倍，使用reads长度是35bp，大约0.39％的染色体上的位置有超过６个测序片段起始；然而，理论上这个比例应该是仅有0.07％。因此我们根据一个经验公式来减少有着共同起始位置的reads的影响。我们使用Illumina 1M BeadChip对炎黄一号样品做了基因分型，并检验了纯合的位点。理论上，经过调整后测序错误出现的频率符合泊松分布。（图４）

（五）人类基因组深度测序的结果评估

我们用炎黄一号个人基因组３６Ｘ测序数据测试了以上方法，在Illumina 1M BeadChip上使用同一ＤＮＡ样本对一致序列进行基因分型。假设所有的基因分型都是正确的，我们可以通过过识别和欠识别来分类基因分型和测序不一样的位点。所谓欠识别就是在等位基因查找过程中少找了一个杂合位点，过识别就是指多找了一个不正确的碱基类型。

图5显示的是一致序列在参考序列上的覆盖度，Illumina 1M BeadChip分型位点的覆盖度和一致序列的错误率。如果没有质量过滤，整个基因组的覆盖度要低于基因分型位点的覆盖度。因为基因分型位点与基因组的特征是有差别的。增加Q0~Q40的质量，基因分型位点的覆盖度有微小的减少：从

98.98%降到98.75%；但是当临界值设的更高时，基因分型位点的覆盖度减少得就比较明显。这个可以用SNP位点的先验概率来解释，结果是相对于其它位点质量分数低了很多。过识别和欠识别的概率随着质量阈值的增加(欠识别：

0.046%Q0~0.024%Q20，过识别：0.096%Q0~0.067%Q20)而连续的减少。根据这些数据得到一个质量阈值为Q20，以平衡覆盖度和欠识别过识别几率。

我们也使用了五个其它的过滤步骤用来去除一些共有序列中不可靠的部分：1）要求单倍体最少两个reads，双倍体场染色体最少4个reads覆盖；2）整体深度，包括重复比对上的reads的覆盖，必须小于100；3）局部序列在基因组上的重复数少于2；4）至少一个双向测序片段支持；5）SNPs间隔至少5bp以上。

1.单倍体染色体X 我们对只含有一个X染色体的雄性基因组进行了测序，因此X染色体的一致序列查找与单倍体基因组是相同的。每个位点有四种不同的可能的基因型。在X染色体的37933个Illumina 1M BeadChip 基因分型位点中，99.59%的位点能够被很好的覆盖，而且基因分型和测序有99.96%的共同性（表2）。一致序列在参考基因组中的X染色体上有88.07%的覆盖度。未有效覆盖的染色体区域是高度重复的，并且几乎没有能够唯一比对的read。Y染色体主要由重复序列组成并且参考序列组装得不是很好，因此我们不对其进行讨论。

2.双倍体常染色体 为了评价一致序列和SNP检测的准确性，我们将炎黄一号所有常染色体组装好的一致序列和NCBI参考序列作了比较，发现一致序列在参考序列的常染色体有92.25%的覆盖，在Illumina 1M BeadChip 基因分型位点上有99.22%的覆盖，其中大约有99.92%位点是一致的（表2）。纯合基因分型位点，有0.062%被测序认为是杂合体。

表 2. 本方法所构建的一致序列在基因分型位点的覆盖度和准确性。

基因分型的碱基型 分型位点数 覆盖度 相同 过识别 欠识别 X染色体 纯合非突变 纯合突变 总数 常染色体 纯合非突变 纯合突变 杂合 总数l 27,196 10,737 37,933 0,878 208,436 250,667 999,981 99.75% 99.14% 99.58% 99.69% 99.26% 98.18% 99.22% 99.99% 99.87% 99.96% 99.94% 99.78% 99.99% 99.92% 0.007% 0.132% 0.042% 0.062% 0.222% 0.013% 0.083% - - - - - 0.103% 0.025% 为评估本方法对SNP的识别准确性，针对部分过识别的SNPs采用PCR扩增后再运用传统的Sanger测序技术进行测定。在57个测试的中，49（86.0%）个碱基型和芯片得出的结果是一致的。

综上所述，我们的SNP检测方法和基因分型芯片的结果有很高的一致性，并且 对于过识别的位点，GA测序技术精确度更高。

附图

图 1.

图 2a.

图 2b.

图 2c.

图 2d.

图 3a.

图 3b.

图 4.

图 5.

图 6a.

图 6b.