专利名称：基于酶切建库测序与贝叶斯统计的单核苷酸多态性位点鉴定的方法
技术领域：
本发明涉及一种基于酶切建库测序与贝叶斯统计的单核苷酸多态性位点鉴定的方法。具体为基于酶切建库单末端(single-end)测序或者双末端(pair_end)测序获得的单核苷酸多态性(SNP)位点进行一种特殊贝叶斯统计检验，从而精确鉴定SNP基因型的方法；可以在缺少参考基因组序列的情况下，为SNP检验提供可靠统计学意义。该方法属于生物信息学技术领域。这对于缺少参考序列的非模式生物的研究以及基因型鉴定的准确性研究具有重要的意义。
背景技术：
SNP(Single Nuc·leotide Polymorphisms)单核苷酸多态性标记，是指在基因组上单个核苷酸的变异，其数量很多，在人类基因组中大概每1000个碱基就有一个SNP，多态性丰富。SNP是进行各类分子生物学研究的理想标记，如构建遗传图谱，基因分型，分子标记辅助育种，疾病预测，用药指导等领域。如今，第二代DNA测序技术是一种高通量低成本的测序技术，基本原理是边合成边测序。以solexa测序方法为例，先用物理方法将DNA链随机打断，然后在片段两端加上特定接头，接头上有扩增引物序列。测序时，DNA聚合酶合成待测片段的互补链，通过检测新合成碱基所携带的荧光信号读取碱基序列，从而获得待测片段的序列。第二代测序技术已经广泛应用于生物科学的许多领域，特别是研究一个物种不同个体之间的多态性。传统寻找SNP标记的方法是将个体进行测序，得到短reads，然后通过短序列比对软件将这些短reads比对回参考序列，从而得到测序个体的SNP信息。常见的流程有(大体过程如图1所示):使用BWA软件将reads比对回参考序列，使用SAMtools软件处理比对结果寻找SNP位点1，2;使用SOAP软件将reads比对回参考序列，使用SOAPsnp软件处理比对结果寻找SNP位点3，4。对于有参考序列的物种都可以很方便的进行SNP标记的查找，但是对于那些非模式生物而言，基本上是不存在参考序列的。而在没有参考序列的情况下，传统寻找SNP标记的方法存在着技术上的瓶颈。RAD测序技术采用了新的建库方式(酶切建库)，其测序具体过程如图2所示，用限制性内切酶切断DNA特定的位点，再用物理方法将酶切之后的DNA分子随机打断，通过琼脂糖胶DNA分离技术挑选特定长度的DNA分子，然后在挑选出来的DNA末端添加特定的扩增接头与测序接头，从而构建上机文库进行高通量测序。其中RAD测序方法为本领域公知的方法，例如可参考文献5，6，7，8。利用RAD测序技术鉴定SNP位点已经在很多领域取得成功，但是截止到本发明出现为止，所用的方法一般都是利用经验阈值进行筛选与过滤。例如文献6中将两种碱基深度比例在
0.25-0.75 (低深度碱基深度:高深度碱基深度)之间的位点判断成杂合位点，比例在0.1以下的判断成纯合位点。这种方法没有统计学意义，同时受到其他外界因素的影响比较大，如测序总量，鉴定得到的SNP基因型准确度无法保证。文献9在经验阈值方法的基础上改进鉴定方法，使用最大似然法进行位点基因型的修正，但是其最大问题在于无法确定统计方法中的错误率。文献:1.Li, H.and R.Durbin, Fast and accurate short read alignment withBurrows-WheeIer transform.Bioinformatics, 2009.25(14):p.1754-60.
2.Li, H., et al., The Sequence Alignment/Map format and SAMtools.Bioinformatics，2009.25 (16):p.2078-9.
3.Li，R.，et al.，SNP detection for massively parallel whole-genomeresequencing.Genome Res, 2009.19(6):p.1124-32.
4.Li, R., et al., SOAP: short oligonucleotide alignment program.Bioinformatics, 2008.24(5):p.713—4.
5.Houston, R.D., et al., Characterisation of QTL-1inked and genome-widerestriction site-associated DNA(RAD)markers in farmed Atlantic salmon.BMCGenomics, 2012.13(1):p.244.
6.Scaglionej D., et al., RAD tag sequencing as a source of SNP markers inCynara cardunculus L.BMC Genomics, 2012.13(1):p.3.
7.Daveyj J.W.，et a 1.，Special features of RAD Sequencingdata:1mplications for genotyping.Mol Ecolj 2012.
8.Dasmahapatraj K.K.，et al.，Butterfly genome reveals promiscuousexchange of mimicry adaptations among species.Nature, 2012.
9.Hohenlohej P.A., et al., Population genomics of parallel adaptation inthreespine stickleback using sequenced RAD tags.PLoS Genet.6 (2):p.e1000862.
发明内容
本发明的目的是提供一种基于酶切建库测序与贝叶斯统计的单核苷酸多态性位点鉴定的方法；它是一种通过处理基于酶切建库测序(RAD测序技术)得到的测序数据，在个体内或个体之间寻找单核苷酸多态性位点，并予以统计学检验的技术方案。本发明的目的通过以下`技术方案来实现:一种基于酶切建库测序与贝叶斯统计的单核苷酸多态性位点鉴定的方法，其步骤如下:I)在获得RAD高通量测序技术的测序结果后，对RAD酶切端测序序列进行过滤以去除不合格的测序序列。其中，RAD高通量测序技术可以为Illumina GA测序技术，也可以为现有的其他高
通量测序技术。所述的不合格的测序序列为测序质量低于预定的低质量阈值的碱基个数超过整条序列碱基个数的50%的序列，以及没有酶切特征的序列。2)根据测序个体基因组酶切一端的测序序列，利用序列的全同性生成每个个体堆的信息，并计算每个序列堆测序深度信息。例如，将每个个体过滤后的酶切一端的测序序列信息作为哈希的键，哈希的值指向一个链表，用于存放另一端的序列信息，并计算测序深度信息。可用任何一种编程语言实现该过程。3)将一个个体内所有序列堆进行两两比对，对堆进行聚类以确定个体内的候选杂合SNP位点。对于3)的聚类结果，只有一个堆的聚类结果表明在酶切一端测序片段上不存在杂合位点，只有两个堆的聚类结果表明在酶切一端测序片段上存在杂合位点。4)将不同个体内所有序列堆进行两两比对，对堆进行聚类以确定不同个体间的候选SNP位点。对于4)的聚类结果,只有一个堆的聚类结果表明两个个体不存在SNP位点，有两个堆的聚类结果表明个体间存在SNP位点。由于后续使用贝叶斯统计方法，所以此处不对每个堆或者每个聚类结果进行深度过滤，是该方法的优势之一，尽可能保留更多的SNP位点。5)利用贝叶斯统计方法对每个候选SNP位点上各种基因型的深度信息进行分析，鉴定候选SNP的准确性。由于缺少参考序列与每个位点上各种碱基出现的先验概率，所以无法获得实际测得的碱基在各位点的错误率。因此，该贝叶斯统计方法中使用步移法穷举所有可能出现的错误率，然后挑选使得基因型存在概率最高的情况作为该SNP位点的基因型。具体公式与计算步骤如下:对于每个候选SNP位点，会存在4中可能的碱基，即“ATCG”中的任一一种或者多种，定义出现频率最高(深度最大)的碱基型为G1，对应的深度为NI，依次递减定义其余碱基，即基因型Gi (i=l，2，3，4)，深度Ni (i=l，2，3，4)。生物学上，一般的物种在一个SNP位点上面只会出现两种碱基型，例如测序数据显示该SNP位点出现A或者T (Na彡Ντ)，那么该位点一定是 纯合或者杂合两种基因型。因此该贝叶斯方法只检测以上两种基因型在错误率为ε的条件下概率为:
■/ja” Λ\." Ji r；-,s) = -二-一 0.754、-—((}.25)n
ΙΛ^υοι」I，」3 I ^M ! M f M I I
i 穿 I，i r m f 3./ 爹 4 */"(,V1, —V,” ;V, (LiS) = -^-(0.5 — (I 25ε) Λ/1 ν:'(0.25)Λ.s/ Λ:/
Lwozi」ρ:， ρ 4."y f ΛΓ ! y ! Λ； !
a Γ -1 2' iV；r其中Ν=Ν1+Ν2+Ν3+Ν4,/，./.Ρ,2,3,4)。各种基因型的后验概率为:
P(miV4|G/)
興.‘)丁上上
■ P(N M2,N3,M4IGY)
itl'/Π.
因为没有序列与前期研究数据，错误率ε无法确定值，但是文献10，11报道Illumina GA测序错误率在1%左右。此处，设定ε从0.01%_5%，步距0.01%。最终使用的后验概率为:P (Nij) Final = max (P (Nij, E))i, j e {1,2, 3,4}, ε e
。如果P (Nij)nnal不小于0.95，表明该SNP位点的基因型为ij，否则定义为无法判断的数据(缺失数据)。本发明的技术方案采用了生物信息学分析方法，处理RAD (restriction-siteassociated DNA)测序数据，寻找RAD测序片段上的SNP位点信息，利用贝叶斯方法鉴定SNP基因型，以突破非模式生物缺少参考序列的瓶颈，在降低成本的同时获得准确的结果。在鉴定SNP位点基因型时首次引入贝叶斯统计方法，与之前经验阈值的方法相比，统计学意义显著提高，准确性也相应提升。文献:10.Li, Y., et al., State of the art de novo assembly of human genomes frommassively parallel sequencing data.Hum Genomics, 2010.4(4):p.271-7.
11.Xie, ff., et al., Parent-1ndependent genotyping for constructing anultrahigh-density linkage map based on population sequencing.Proc Natl Acad SciUSA.107(23):p.10578-83.
图1为传统寻找SNP位点方法的过程原理图；图2为RAD测序技术的测序具体过程原理图；图中，(A)限制性内切酶酶切基因组DNA，并加上Pl接头，每一个Pl接头含有不同的序列标签；(B)带有不同Pl接头的样品混合在一起，打断；(C)加上接头P2 ；⑶扩增富集RAD tags ；图3为RAD测序的例图；图4为生成堆的聚类过程图，图例中使用EcorI限制性内切酶；图5为堆中酶切一端序列信息示意图；图6为堆中个体内与个体间SNP位点寻找流程示意图；图7为候选SNP碱基型与深度信息的例图，图中有20个候选SNP位点分别在15个个体中的碱基型与深度 信息，“CI 9”表示该位点C测到9次，“CI 9: T I 3”表示该位点C测到9次、T测到3次。图8为贝叶斯统计后SNP位点的基因型结果的例图，图中“a”和“b”分别表示两种不同的纯合基因型，“h”表示杂合基因型，例如“a”表示AA，“b”表示CC，“h”表示AC，“X1:X2:X3:X4”中Xl表示纯合的后验概率，x2表示纯合后验概率最大时的错误率值，x3表示杂合的后验概率，x4表示杂合后验概率最大时的错误率值。图9为利用经验阈值方法与贝叶斯统计方法后数据的缺失情况。图10为经验阈值方法与贝叶斯统计方法不同位点的状态统计结果。图11为随机挑选经验阈值方法与贝叶斯统计方法不同的位点进行利用sanger测序验证的结果。
具体实施例方式下面结合附图与具体实施例进一步阐述本发明的技术特点。如图2所示，RAD测序与常规高通量测序不同的是在加接头之前需要利用限制性内切酶完全酶切基因组，然后加上RAD特异接头，后续建库过程与常规建库完全相同。图3为RAD酶切末端测序的例图。在图3中显示了用限制性内切酶EcorI，识别DNA分子上“G'AATTC”的回文序列，并在G与A之间将DNA分子切断，将酶切后的DNA分子用物理方法打断成短的序列片段，并在含有酶切末端序列的DNA片段两端加上接头，PCR富集并进行单末端(single-end)测序,测序读长一般为IOOnt,也可以为50nt。基于酶切建库测序与贝叶斯统计的单核苷酸多态性位点鉴定的方法，其步骤如下:I)在获得RAD高通量测序技术的测序结果后，对RAD双末端测序序列进行过滤以去除不合格的测序序列。其中，RAD高通量测序技术可以为Illumina GA测序技术，也可以为现有的其他高
通量测序技术。所述的不合格的测序序列为测序质量低于预定的低质量阈值的碱基个数超过整条序列碱基个数的50%的序列，以及没有酶切特征的序列。低质量阈值由具体测序技术及测序环境而定，例如设定为单碱基测序质量低于20 ;测序序列中测序结果不确定的碱基(如Illumina GA测序结果中的N)个数超过整条测序序列碱基个数的10%则认为是不合格序列；除样本接头序列外，与其它实验引入的外源序列比对，如各种接头序列。若序列中存在外源序列则认为是不合格序列；在酶切一端测序序列中，若起始的几个 碱基不是酶切末端序列则过滤掉(如限制性内切酶Ecorl，测序序列开头若不是“AATTC”则过滤掉整个测序序列)。2)根据测序个体基因 组酶切一端的测序序列，利用序列的全同性生成每个个体堆的信息。具体过程如图4所示。堆(Stack)中酶切一端的序列信息可以以图5的方式保存，在图5中，第一列表不的是酶切一端的序列信息；第二列表不的是酶切一端序列被测序的次数，即测序深度信息；第三列是该堆的ID，用于唯一确定堆。3)个体内堆进行比较，如果聚类没有出现图6中的情况，表明个体内部这些序列上面没有杂合SNP ;如果存在图6中的情况，表明个体内部存在杂合SNP。4)个体间堆进行两两比较时，如果出现图6 Ca)情况，表明两个个体间存在纯合SNP;如果出现图6 (b)情况，表明两个个体间存在杂合SNP。5)堆聚类与比较结束后获得每个候选SNP位点的碱基型与相应深度信息，如图7表示上述信息汇总。6)利用贝叶斯统计方法对碱基型与深度进行分析，如图8，得到最终的SNP基因型结果。实施例数据:瓠瓜F2群体构建遗传图谱项目，包括139个F2植株和两个亲本，对这141个个体进行RAD-PE测序。(说明:父本和母本杂交生成的后代叫Fl，Fl自交生成的后代叫F2 ;虽然使用RAD-PE测序，但是分析仅使用酶切端序列)材料来源:浙江省农业科学院。实施例具体操作流程:将两个亲本RAD-PE测序得到的测序数据，根据测序质量值，N的含量，以及是否含有酶切末端序列进行过滤，去除不合格的测序序列，得到的有效数据统计如表I所示。 表1:瓠瓜RAD测序有效数据统计
权利要求
1.基于酶切建库测序与贝叶斯统计的单核苷酸多态性位点鉴定的方法，其步骤如下: 1)在获得RAD高通量测序技术的测序结果后，对RAD测序序列进行过滤以去除不合格的测序序列； 2)利用序列的全同性生成每个个体序列堆的信息，并计算每个序列堆测序深度信息； 3)将一个个体内所有序列堆进行两两比对，对堆进行聚类以确定个体内的候选杂合SNP位点； 4)将不同个体内所有序列堆进行两两比对，对堆进行聚类以确定不同个体间的候选SNP位点； 5)利用贝叶斯统计方法对每个候选SNP位点上各种基因型的深度信息进行分析，鉴定候选SNP的准确性，用于后续群体分析或者实验等工作。
2.根据权利要求1所述的基于酶切建库测序与贝叶斯统计的单核苷酸多态性位点鉴定的方法，其特征在于:步骤I)中，RAD高通量测序技术为Illumina GA测序技术；测序类型默认是单末端测序，如果是双末端测序，只对酶切端进行单核苷酸多态性位点分析。
3.根据权利要求1所述的基于酶切建库测序与贝叶斯统计的单核苷酸多态性位点鉴定的方法，其特征在于:步骤I)中，所述的不合格的测序序列为测序质量低于预定的低质量阈值的碱基个数超过整条序列碱基个数的50%的序列，以及起始处没有酶切特征序列的序列。
4.根据权利要求1所述的基于酶切建库测序与贝叶斯统计的单核苷酸多态性位点鉴定的方法，其特征在于:步骤3)中，各序列堆相互比对时，仅允许出现两类不同的序列堆，过滤掉超过两类的序列堆。
5.根据权利要求1所述的基于酶切建库测序与贝叶斯统计的单核苷酸多态性位点鉴定的方法，其特征在于:步骤4)中， 保留所有比对获得的序列堆，不对超过两类的序列堆进行处理。
全文摘要
本发明公开了基于酶切建库测序与贝叶斯统计的单核苷酸多态性位点鉴定的新的算法。它用于处理RAD测序数据，寻找RAD测序片段上的候选SNP，采用贝叶斯统计的生物信息学分析方法鉴定SNP可靠性。该方法一方面可以用于模式与非模式生物，摆脱大量物种缺少参考序列的限制，减少了测序成本；另一方面解决了目前利用RAD数据进行SNP鉴定时无可靠统计方法的瓶颈，使获得的SNP位点准确性极大提升。
文档编号C12Q1/68GK103114150SQ201310077509
公开日2013年5月22日 申请日期2013年3月11日 优先权日2013年3月11日
发明者陶晔, 钱刚, 郑泽群, 胡秋萍 申请人:上海美吉生物医药科技有限公司