Sanger法测序中确定单核苷酸多态性的方法

摘要

本发明提供一种Sanger法测序中确定单核苷酸多态性的方法，所述方法包括将来源于同一个待测序列的所述高质量的碱基序列文件通过构建kmer哈希表，以kmer序列为哈希表的键，以kmer序列在整条序列中出现的次数为该键所对应的值，使用kmer值进行所述高质量碱基序列拼接组装。本发明的Sanger法测序中确定单核苷酸多态性的方法可以实现全分析流程自动化，只需设置好相关参数，中间无需人工参与，直接输出结果文件。

说明书

技术领域

本发明涉及生物技术领域，具体涉及Sanger法测序中确定单核苷酸多态性的方法。

背景技术

Sanger测序是DNA测序技术的“金标准”，曾在人类基因组计划中发挥了关键推动作用，并且在现在仍被用来获得高度准确且可信赖的测序数据。

在Sanger测序期间，DNA聚合酶通过向生长链(延伸产物)中加入核苷酸来复制单链DNA模板。链延长发生在引物的3′端，通过与模板的碱基对互补来选择加入到扩增产物中的脱氧核苷酸。

Sanger测序利用一种DNA聚合酶来延伸结合在待定序列模板上的引物。直到掺入一种链终止核苷酸为止。每一次序列测定由一套四个单独的反应构成，每个反应含有所有四种脱氧核苷酸三磷酸(dNTP)，并混入限量的一种不同的双脱氧核苷三磷酸(ddNTP)。由于ddNTP缺乏延伸所需要的3-OH基团，使延长的寡聚核苷酸选择性地在G、A、T或C处终止。终止点由反应中相应的双脱氧而定。每一种dNTPs和ddNTPs的相对浓度可以调整，使反应得到一组长几百至几千碱基的链终止产物。

目前，Sanger测序数据的碱基质量需要根据abi文件中每个碱基对应的峰图值进行人工辨别以及手工去除低质量的碱基，具有一定的人为操作误差，同时现有的方法进行此步骤处理时，只能一条序列进行处理，结束以后输出结果再进行下一条序列的操作，严重影响了分析进度。而进行序列之间的拼接时，同样存在上述问题，一次只能进行一个样本的序列拼接，分析结果完成后才能进行下一个样本的分析。而针对一代Sanger测序数据进行SNP检测时，需要根据abi文件中同一位置的多个峰图值进行辨别，只有当次峰值与主峰值达到一定比值范围，才能判定此碱基位置存在SNP变异，由于目前方法中峰值无法进行数字量化，只能人工根据峰图高度进行估算，存在严重的人为操作误差，如果序列较长时，进行人工SNP判别需要花费大量时间，无法提高工作效率。而且基于一代Sanger测序数据的分析方法，无法在已有方法基础上进行改进以实现质控、拼接以及SNP检测的全自动化分析。

发明内容

在一种实施方式中，本发明提供一种Sanger法测序中确定单核苷酸多态性的方法，所述方法包括以下步骤：步骤1：根据每个待测序列的长度，设计N对引物进行PCR扩增，N为不小于2的整数，N对引物能够完全覆盖待测序列；步骤2：对每个待测序列的扩增产物进行Sanger法测序，每个待测序列生成2N个Sanger测序abi文件，每个待测序列的abi文件进行命名，以便于根据该命名来识别来源于同一个待测序列；步骤3：将所述Sanger测序abi文件转化成文本格式文件，并对碱基信号值做归一化处理；步骤4：通过滑窗方法删除低于预设碱基质量值的序列，同步删除低于预设碱基质量值的序列区域所对应的碱基质量值，获得高质量碱基序列及相应的碱基质量值；步骤5：将来源于同一个待测序列的所述高质量的碱基序列文件通过构建kmer哈希表，以kmer序列为哈希表的键，以kmer序列在整条序列中出现的次数为该键所对应的值，使用kmer值进行所述高质量碱基序列拼接组装；和步骤6：基于步骤4得到的高质量碱基序列及相应的碱基质量值，获得每个碱基位点的次最大碱基信号值和最大碱基信号值的比值，当该值大于预设值时，评估每条拼接组装后的序列的多态性位点，然后储存并输出每条待测序列的多态性位点信息。

在一种实施方式中，所述滑窗方法中使用的滑窗范围为5-20bp，优选地为5-10bp。

在一种实施方式中，所述预设碱基质量值为30-60，优选地为质量值范围为50-60。

在一种实施方式中，在步骤2中，所述待测序列的每个abi文件以待测序列名称+引物名称方式进行命名。

在一种实施方式中，将来源于同一个待测序列的高质量序列文件按照待测序列名称中的引物名称进行排序，分别构建相邻的两条待拼接序列的kmer哈希表，以kmer序列为哈希表的键，以kmer序列在整条序列中出现的次数为该键所对应的值，分别从相邻两条序列所对应的哈希表的键中检索是否存在相同的代表kmer序列的键，并且该键仅对应唯一的值，当两条序列中不存在相同的键时，对kmer值减1，继续进行搜索，直至找到最大kmer值为止，基于相邻两条序列中同时存在且唯一的所有kmer序列所对应的位置信息，定位出两序列之间最大的重叠区间，获得该区间在两条序列中的位置信息。

在一种实施方式中，所述kmer值范围为90-150bp。

在一种实施方式中，所述预设值不小于0.25。

在一种实施方式中，步骤6中，单核苷酸多态性位点识别结果以Excel文件格式输出，记录单核苷酸多态性位点在拼接序列中的坐标位置以及相对应的碱基。

采用本发明的Sanger法测序中确定单核苷酸多态性的方法，可以在普通windows系统电脑上实现自动读取数据输入文件夹中的多个abi格式序列文件，根据文件名称将来源于同一样本的多条序列进行滑窗法自动过滤两端低质量碱基，完成序列拼接后进行SNP检测，输出SNP位点简并碱基，拼接序列及对应峰图结果到数据结果文件夹中，然后依次处理输入数据文件夹中的数据文件，直至完成所有数据文件。

本发明的Sanger法测序中确定单核苷酸多态性的方法可以实现全分析流程自动化，只需设置好相关参数，中间无需人工参与，直接输出结果文件。此方法可以避免在序列质控过滤以及SNP检测过程中人工判读产生的误差，得到的结果更精确。同时全自动流程可以大量节约分析时间，如果采用商业化方法进行序列质控和拼接，每次只能输入一个样本的序列进行分析，然后人工判读SNP位点信息后再输出相应的结果，一个样本大约需要3-5分钟的时间，而采用本发明方法只需5-6秒的时间即可完成同样的工作，极大的提升了工作效率。

附图说明

为了更清楚地说明本申请实施例中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施例，对于本领域普通技术人员来说，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。

图1是样本LSYZ2017020_SeqF3原始下机序列及滑窗长度(win)分别为5bp、10bp、15bp和20bp时进行低质量碱基去除后序列的5’和3’端峰图；

图2是样本LSYZ2017020_SeqF3原始下机序列及固定滑窗长度(win)为5bp，碱基质量的平均值低于预设的碱基质量值(QC)30，40，50和60时被去除的序列的5’和3’端峰图；

图3是样本LSYZ2017020的六条过滤后高质量序列文件，按照排序SeqR3-SeqR4-SeqRT-R2-SeqPRT-F2-SeqF3-SeqF4进行序列拼接组装，Kmer分别设置为20，25，30，35，40，45，50，60，70，80，90，100，150，160，180和200的组装结果序列展示图；

图4是样本LSYZ2017020的六条过滤后高质量序列文件拼接组装示意图及实际数据组装展示图；

图5是样本LSYZ2017020的不同SNP cutoff值检测SNP位点峰图结果展示图；

图6是样本LSYZ2017020的数据拼接及单核苷酸多态性位点识别结果展示图；

图7是样本LSYZ2017032_SeqF3原始下机序列及滑窗长度(win)分别为5bp、10bp、15bp和20bp时进行低质量碱基去除后序列的5’和3’端峰图；

图8是样本LSYZ2017032_SeqF3原始下机序列及固定滑窗长度(win)为5bp，碱基质量的平均值低于预设的碱基质量值(QC)30，40，50和60时被去除的序列的5’和3’端峰图；

图9是LSYZ2017032的六条过滤后高质量序列文件，按照排序SeqR3-SeqR4-SeqRT-R2-SeqPRT-F2-SeqF3-SeqF4进行序列拼接组装，Kmer分别设置为20，25，30，35，40，45，50，60，70，80，90，100，150，160，180和200的组装结果序列展示图；

图10是LSYZ2017041的六条过滤后高质量序列文件，按照排序SeqR3-SeqR4-SeqRT-R2-SeqPRT-F2-SeqF3-SeqF4进行序列拼接组装，Kmer分别设置为20，25，30，35，40，45，50，60，70，80，90，100，150，160，180和200的组装结果序列展示图；

图11是样本LSYZ2017032的六条过滤后高质量序列文件拼接组装示意图及实际数据组装展示图；

图12是样本LSYZ2017041的六条过滤后高质量序列文件拼接组装示意图及实际数据组装展示图；

图13是样本LSYZ2017032的不同SNP cutoff值检测SNP位点峰图结果展示图；

图14是样本LSYZ2017041的不同SNP cutoff值检测SNP位点峰图结果展示图；和

图15是样本LSYZ2017032和样本LSYZ2017041的数据拼接及单核苷酸多态性位点识别结果展示图。

具体实施方式

为了使本领域技术领域人员更好地理解本申请中的技术方案，下面将结合实施例对本发明作进一步说明，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都应当属于本申请保护的范围。下面结合附图及实施例对本发明作进一步描述。

实施例一：本发明方法测单序列中单核苷酸多态性的方法

一、待测目标序列长度：～1200bp，设计3对测序引物，具体如表1所示。

表1扩增引物信息表

引物名称碱基构成(5′-3′)位置(HXB2)PRT-F2CTTTARCTTCCCTCARATCACTCT2243-2266RT-R2CTTCTGTATGTCATTGACAGTCC3326-3304SeqF3AGTCCTATTGARACTGTRCCAG2556-2577SeqR3TTTYTCTTCTGTCAATGGCCA2639-2619SeqF4CAGTACTGGATGTGGGRGAYG2869-2889SeqR4TACTAGGTATGGTAAATGCAGT2952-2931

二、对待测序列的扩增产物进行Sanger法测序，获得六条Sanger测序abi文件，文件名以样本名称+引物名称组成，序列名称LSYZ2017020-(SeqPRT-F2)，LSYZ2017020-(SeqF3)，LSYZ2017020-(SeqF4)，LSYZ2017020-(SeqRT-R2)，LSYZ2017020-(SeqR3)和LSYZ2017020-(SeqR4)，将上述序列文件放置于同一数据输入文件夹中，计算机根据序列文件名称自动识别属于待测序列的测序文件，并进行后续分析。

三、将所述Sanger测序abi文件转化成文本格式文件，并对碱基质量信号值做归一化处理。

四、通过滑窗方法删除低于预设碱基质量值的序列，同步删除低于预设碱基质量值的序列区域所对应的碱基质量值，获得高质量碱基序列及相应的碱基质量值，具体过程如下所示：

1.设置滑窗长度分别为5bp，10bp，15bp和20bp检测序列碱基质量平均值，当碱基质量的平均值大于预设的碱基质量值(默认值为50)时，停止滑动，删除低于预设碱基质量值的序列，同步删除低碱基质量序列区域所对应的相关信息，获得高质量的碱基序列及其相关数据信息；采用此方法依次处理数据输入文件夹中的六条Sanger测序abi文件：LSYZ2017020-(SeqPRT-F2)，LSYZ2017020-(SeqF3)，LSYZ2017020-(SeqF4)，LSYZ2017020-(SeqRT-R2)，LSYZ2017020-(SeqR3)和LSYZ2017020-(SeqR4)。

以样本LSYZ2017020_SeqF3为例，将原始测序序列分别设置长度为5bp，10bp，15bp和20bp滑窗，检测序列碱基质量平均值，当碱基质量的平均值大于预设的碱基质量值(默认值为50)时，停止滑动，删除低于预设碱基质量值的序列，同步删除低碱基质量序列区域所对应的相关信息，获得高质量的碱基序列及其相关数据信息，并输出峰图文件。由于测序序列长度较长，不方便完全展示，为了更好的比较不同长度滑窗的质量评估效果，分别截取原始测序序列和5bp，10bp，15bp和20bp长度滑窗去除低质量碱基后序列的5’和3’端峰图进行平行比较，如图1所示。下述相关序列峰图比较均采用此种截取5’和3’端峰图方式，其余序列LSYZ2017020-SeqF4、LSYZ2017020-SeqPRT-F2、LSYZ2017020-SeqR3、LSYZ2017020-SeqR4和LSYZ2017020-SeqRT-R2依次按照上述方法进行低质量碱基去除及相关对应峰图比较。

根据原始序列和四种不同长度滑窗进行碱基质量评估以及去除低质量碱基以后的5’和3’端峰图结果进行比较分析，依据去除低质量碱基的数量以及保留高质量碱基的峰值信息，数据质量过滤时建议使用滑窗范围为5-10bp。

2.设置滑窗长度为5bp检测序列碱基质量平均值，当碱基质量的平均值大于预设的碱基质量值30，40，50和60时，停止滑动，删除低于预设碱基质量值的序列，同步删除低碱基质量序列区域所对应的相关信息，获得高质量的碱基序列及其相关数据信息；采用此方法依次处理数据输入文件夹中的六条Sanger测序abi文件：LSYZ2017020-(SeqPRT-F2)，LSYZ2017020-(SeqF3)，LSYZ2017020-(SeqF4)，LSYZ2017020-(SeqRT-R2)，LSYZ2017020-(SeqR3)和LSYZ2017020-(SeqR4)。

以样本LSYZ2017020SeqF3为例，设置滑窗长度为5bp检测序列碱基质量平均值，当碱基质量的平均值大于预设的碱基质量值30，40，50和60时，停止滑动，删除低于预设碱基质量值的序列，同步删除低碱基质量序列区域所对应的相关信息，获得高质量的碱基序列及其相关数据信息，并输出峰图文件。分别截取原始测序序列和碱基质量的平均值低于预设的碱基质量值30，40，50和60时被去除的序列的5’和3’端峰图进行平行比较，如图2所示。其余序列LSYZ2017020-SeqF4、LSYZ2017020-SeqPRT-F2、LSYZ2017020-SeqR3、LSYZ2017020-SeqR4和LSYZ2017020-SeqRT-R2依次按照上述方法进行低质量碱基去除及相关对应峰图比较。根据原始序列和四种不同碱基质量值进行碱基质量评估以及去除低质量碱基以后的5’和3’端峰图结果进行比较分析，依据去除低质量碱基的数量以及保留高质量碱基的峰值信息，数据质量过滤时建议使用质量值范围为50-60。

五、将来源于同一个待测序列的所述高质量的碱基序列文件通过构建kmer哈希表，以kmer序列为哈希表的键，以kmer序列在整条序列中出现的次数为该键所对应的值，使用kmer值进行所述高质量碱基序列拼接组装。

1.将来源于一个样本的六条过滤后高质量序列文件按照样本名称中引物名称进行排序(seq1-seq2-seq3-seq4-seq5-seq6)，分别构建两条待拼接序列(seql-seq2，seq2-seq3，seq3-seq4等)的kmer哈希表，以kmer序列为哈希表的键，以kmer序列在整条序列中出现的次数为该键所对应的值，分别从两条序列所对应的哈希表的键中检索是否存在相同的代表kmer序列的键，并且该键仅对应唯一的值，当两条序列中不存在相同的键时，对kmer值减1，继续进行搜索，直至找到最大kmer值为止，基于两条序列中同时存在且唯一的所有kmer序列所对应的位置信息，定位出两序列之间最大的重叠区间，获得该区间在两条序列中的位置信息；

以LSYZ2017020样本为例，通过识别来源于此样本的六条过滤后高质量序列文件，按照样本名称中引物名称进行排序SeqR3-SeqR4-SeqPRT-F2-SeqRT-R2-SeqF3-SeqF4，如图3所示；然后按照上述方法进行序列拼接组装，Kmer分别设置为20，25，30，35，40，45，50，60，70，80，90，100，150，160，180和200，评估序列组装结果。

将16种Kmer组装结果序列采用MEGA软件进行比对排序，结果如图3显示所有Kmer条件下，序列组装结果均为长度1061bp。图片左侧栏展示不同Kmer组装结果的序列名称，后侧栏展示组装的序列碱基比对排序后结果，ATCG四种碱基分别标识为绿色、红色、淡蓝色以及紫色等四种颜色，简并碱基无颜色标识，如果所有比对序列在同一位置的碱基一致，则在图片的上侧栏标识*，否则为空白。比对拼接组装结果发现，16种Kmer组装拼接结果长度完全一致，无差异，但是组装序列中SNP识别有差异。

表2 LSYZ2017020样本不同Kmer组装结果SNP统计表

将16种Kmer组装拼接结果中识别的SNP位点及对应碱基进行统计，如表2所示，左侧栏为SNP在序列中的位置坐标，上侧栏为16种Kmer值，后侧栏为SNP对应碱基，Ref表示参考序列碱基，Alt表示组装序列的突变碱基。从上表分析发现，Kmer值大于150bp以后，存在一些位点SNP未被检出的情况，结合多个样本的测试结果，建议使用本方法进行序列拼接组装时使用kmer取值范围为90-150bp。

2.六条高质量序列分别拼接组装结果展示

如图4，来源于同一样本的六条高质量序列按照样本名称中引物名称进行排序SeqR3-SeqR4-SeqPRT-F2-SeqRT-R2-SeqF3-SeqF4(图4-A)，然后按照上述Kmer方法进行序列拼接组装，图4-B展示SeqR3、SeqR4和SeqPRT-F2序列5’端比对一致性情况，依据拼接方案，此区域比对序列在组装后生成一条一致性序列，即为拼接组装后序列，图4-C展示SeqR3序列3’端、SeqR4序列3’端、SeqPRT-F2序列5’端和SeqRT-R2序列5’端比对一致性情况，图4-D展示SeqR4序列3’端、SeqPRT-F2序列3’端、SeqRT-R2序列3’端、SeqF3序列3’端和SeqF4序列5’端比对一致性情况。

六、基于步骤4得到的高质量碱基序列及相应的碱基质量值，获得每个碱基位点的次最大碱基信号值和最大碱基信号值的比值，当该值大于预设值时，评估每条拼接组装后的序列的多态性位点，然后储存并输出每条待测序列的多态性位点信息。

1.基于步骤4中过滤后得到的高质量序列及相关数据文件，获得每个碱基位点的次最大碱基信号值和最大碱基信号值的比值，当该值大于预设值(分别设置为0.2，0.25，0.33，0.5)时，评估多态性位点；

根据四种不同的SNP评估预设值(cut off)，识别组装序列中的SNP位点及对应的碱基，统计结果如表3所示，左侧栏为SNP所在序列中的坐标位置，上侧栏为预设值，后侧栏为检测SNP位点对应的碱基，Ref表示参考序列碱基，Alt表示组装序列的突变碱基，表示未检测到SNP。分析表3中同一位置在不同预设值0.2、0.25、0.33和0.5时，SNP检测结果有差异，选取坐标位置为5、401、463、554和812的碱基峰图(图5)进行详细分析，图中左侧栏为预设值，后侧栏为5个不同坐标位置对应碱基的峰图结果，根据实际峰图结果判定，分析进行多态性位点分析时建议使用cut off值最低为0.25。

表3不同SNP cutoff值检测结果

2.Sanger测序数据拼接及单核苷酸多态性位点识别结果

结果输出文件以输入文件名命名(图6-A)，将同一样本来源的六条序列拼接结果以fasta格式输出(图6-B)，其中单核苷酸多态性位点以简并碱基的形式记录，拼接结果中碱基质量值转换成图形信号以pdf格式文件输出(图6-C)，单核苷酸多态性位点以蓝色条柱标识，同时单核苷酸多态性位点识别结果以Excel文件格式输出(图6-D)，记录SNP位点在拼接序列中的坐标位置以及相对应的碱基。

实施例二：本发明方法测二个样本序列中单核苷酸多态性的方法。

一、以样本LSYZ2017032和LSYZ2017041为例，每个样本测序获得六条Sanger测序abi文件，文件名以样本名称+引物名称组成，序列名称分别为LSYZ2017032-(SeqPRT-F2)，LSYZ2017032-(SeqF3)，LSYZ2017032-(SeqF4)，LSYZ2017032-(SeqRT-R2)，LSYZ2017032-(SeqR3)和LSYZ2017032-(SeqR4)；LSYZ2017041-(SeqPRT-F2)，LSYZ2017041-(SeqF3)，LSYZ2017041-(SeqF4)，LSYZ2017041-(SeqRT-R2)，LSYZ2017041-(SeqR3)和LSYZ2017041-(SeqR4)。将上述序列文件放置于同一数据输入文件夹中，计算机根据序列文件名称自动识别属于待测序列的测序文件，并进行后续分析。

二、将所述Sanger测序abi文件转化成文本格式文件，并对碱基质量信号值做归一化处理。

三、通过滑窗方法删除低于预设碱基质量值的序列，同步删除低于预设碱基质量值的序列区域所对应的碱基质量值，获得高质量碱基序列及相应的碱基质量值：

1.设置滑窗长度分别为5bp，10bp，15bp和20bp检测序列碱基质量平均值，当碱基质量的平均值大于预设的碱基质量值(默认值为50)时，停止滑动，删除低于预设碱基质量值的序列，同步删除低碱基质量序列区域所对应的相关信息，获得高质量的碱基序列及其相关数据信息；采用此方法依次处理数据输入文件夹中的12条Sanger测序abi文件：LSYZ2017032-(SeqPRT-F2)，LSYZ2017032-(SeqF3)，LSYZ2017032-(SeqF4)，LSYZ2017032-(SeqRT-R2)，LSYZ2017032-(SeqR3)和LSYZ2017032-(SeqR4)；LSYZ2017041-(SeqPRT-F2)，LSYZ2017041-(SeqF3)，LSYZ2017041-(SeqF4)，LSYZ2017041-(SeqRT-R2)，LSYZ2017041-(SeqR3)和LSYZ2017041-(SeqR4)。

以样本LSYZ2017032_SeqF3为例，将原始测序序列分别设置长度为5bp，10bp，15bp和20bp滑窗，检测序列碱基质量平均值，当碱基质量的平均值大于预设的碱基质量值(默认值为50)时，停止滑动，删除低于预设碱基质量值的序列，同步删除低碱基质量序列区域所对应的相关信息，获得高质量的碱基序列及其相关数据信息，并输出峰图文件。由于测序序列长度较长，不方便完全展示，为了更好的比较不同长度滑窗的质量评估效果，分别截取原始测序序列和5bp，10bp，15bp和20bp长度滑窗去除低质量碱基后序列的5’和3’端峰图进行平行比较(图7)，下述相关序列峰图比较均采用此种截取5’和3’端峰图方式，其余序列LSYZ2017032-SeqF4、LSYZ2017032-SeqPRT-F2、LSYZ2017032-SeqR3、LSYZ2017032-SeqR4和LSYZ2017032-SeqRT-R依次按照上述方法进行低质量碱基去除及相关对应峰图比较。样本LSYZ2017041的六条序列也均采用上述方法进行低质量碱基去除及相关对应峰图比较。

根据原始序列和四种不同长度滑窗进行碱基质量评估以及去除低质量碱基以后的5’和3’端峰图结果进行比较分析，依据去除低质量碱基的数量以及保留高质量碱基的峰值信息，数据质量过滤时建议使用滑窗范围为5-10bp。

2.设置滑窗长度为5bp检测序列碱基质量平均值，当碱基质量的平均值大于预设的碱基质量值30，40，50和60时，停止滑动，删除低于预设碱基质量值的序列，同步删除低碱基质量序列区域所对应的相关信息，获得高质量的碱基序列及其相关数据信息；采用此方法依次处理数据输入文件夹中的12条Sanger测序abi文件，LSYZ2017032-(SeqPRT-F2)，LSYZ2017032-(SeqF3)，LSYZ2017032-(SeqF4)，LSYZ2017032-(SeqRT-R2)，LSYZ2017032-(SeqR3)和LSYZ2017032-(SeqR4)；LSYZ2017041-(SeqPRT-F2)，LSYZ2017041-(SeqF3)，LSYZ2017041-(SeqF4)，LSYZ2017041-(SeqRT-R2)，LSYZ2017041-(SeqR3)和LSYZ2017041-(SeqR4)。

以样本LSYZ2017032_SeqF3为例，设置滑窗长度为5bp检测序列碱基质量平均值，当碱基质量的平均值大于预设的碱基质量值30，40，50和60时，停止滑动，删除低于预设碱基质量值的序列，同步删除低碱基质量序列区域所对应的相关信息，获得高质量的碱基序列及其相关数据信息，并输出峰图文件。分别截取原始测序序列和碱基质量的平均值低于预设的碱基质量值30，40，50和60时被去除的序列的5’和3’端峰图进行平行比较(图8)，其余序列LSYZ2017032-SeqF4、LSYZ2017032-SeqPRT-F2、LSYZ2017032-SeqR3、LSYZ2017032-SeqR4和LSYZ2017032-SeqRT-R2依次按照上述方法进行低质量碱基去除及相关对应峰图比较。样本LSYZ2017041的六条序列也均采用上述方法进行低质量碱基去除及相关对应峰图比较。

根据原始序列和四种不同碱基质量值进行碱基质量评估以及去除低质量碱基以后的5’和3’端峰图结果进行比较分析，依据去除低质量碱基的数量以及保留高质量碱基的峰值信息，数据质量过滤时建议使用质量值范围为50-60。

四、将来源于同一个待测序列的所述高质量的碱基序列文件通过构建kmer哈希表，以kmer序列为哈希表的键，以kmer序列在整条序列中出现的次数为该键所对应的值，使用kmer值进行所述高质量碱基序列拼接组装：

1.将来源于一个样本的六条过滤后高质量序列文件按照样本名称中引物名称进行排序(seq1-seq2-seq3-seq4-seq5-seq6)，分别构建两条待拼接序列(seq1-seq2，seq2-seq3，seq3-seq4等)的kmer哈希表，以kmer序列为哈希表的键，以kmer序列在整条序列中出现的次数为该键所对应的值，分别从两条序列所对应的哈希表的键中检索是否存在相同的代表kmer序列的键，并且该键仅对应唯一的值，当两条序列中不存在相同的键时，对kmer值减1，继续进行搜索，直至找到最大kmer值为止，基于两条序列中同时存在且唯一的所有kmer序列所对应的位置信息，定位出两序列之间最大的重叠区间，获得该区间在两条序列中的位置信息；

以LSYZ2017032和LSYZ2017041样本为例，通过分别识别来源于样本的六条过滤后高质量序列文件，按照样本名称中引物名称进行排序SeqR3-SeqR4-SeqPRT-F2-SeqRT-R2-SeqF3-SeqF4，然后按照上述方法进行序列拼接组装，Kmer分别设置为20，25，30，35，40，45，50，60，70，80，90，100，150，160，180和200，评估序列组装结果。

将LSYZ2017032和LSYZ2017041样本的各16种Kmer组装结果序列采用MEGA软件进行比对排序，结果显示(图9和10)所有Kmer条件下，序列组装结果长度分别为1056bp和1057bp。图片左侧栏展示不同Kmer组装结果的序列名称，后侧栏展示组装的序列碱基比对排序后结果，ATCG四种碱基分别标识为绿色、红色、淡蓝色以及紫色等四种颜色，简并碱基无颜色标识，如果所有比对序列在同一位置的碱基一致，则在图片的上侧栏标识*，否则为空白。比对拼接组装结果发现，每个样本的16种Kmer组装拼接结果长度完全一致，无差异，但是组装序列中SNP识别有差异。

表4 LSYZ2017032样本不同Kmer组装结果SNP统计表

表5 LSYZ2017041样本不同Kmer组装结果SNP统计表

分别将两个样本的16种Kmer组装拼接结果中识别的SNP位点及对应碱基进行统计，如表4和5所示，左侧栏为SNP在序列中的位置坐标，上侧栏为16种Kmer值，后侧栏为SNP对应碱基，Ref表示参考序列碱基，Alt表示组装序列的突变碱基。从上表分析发现，Kmer值大于150bp以后，存在一些位点SNP未被检出的情况，结合多个样本的测试结果，建议使用本发明方法进行序列拼接组装时使用kmer取值范围为90-150bp。

来源于同一样本的六条高质量序列按照样本名称中引物名称进行排序SeqR3-SeqR4-SeqPRT-F2-SeqRT-R2-SeqF3-SeqF4(图11-A和12A)，然后按照上述Kmer方法进行序列拼接组装，图11-B和12-B展示SeqR3、SeqR4和SeqPRT-F2序列5’端比对一致性情况，依据拼接方案，此区域比对序列在组装后生成一条一致性序列，即为拼接组装后序列，图11-C和12-C展示SeqR3序列3’端、SeqR4序列3’端、SeqPRT-F2序列5’端和SeqRT-R2序列5’端比对一致性情况，图11-D和12-D展示SeqR4序列3’端、SeqPRT-F2序列3’端、SeqRT-R2序列3’端、SeqF3序列3’端和SeqF4序列5’端比对一致性情况。

五、基于步骤4中过滤后得到的高质量序列及相关数据文件，获得每个碱基位点的次最大碱基信号值和最大碱基信号值的比值，当该值大于预设值(分别设置为0.2，0.25，0.33，0.5)时，评估多态性位点。

根据四种不同的SNP评估预设值(cut off)，识别组装序列中的SNP位点及对应的碱基，LSYZ2017032和LSYZ2017041样本统计结果如表6和7所示，左侧栏为SNP所在序列中的坐标位置，上侧栏为预设值，后侧栏为检测SNP位点对应的碱基，Ref表示参考序列碱基，Alt表示组装序列的突变碱基，表示未检测到SNP。分析表6样本LSYZ2017032中同一位置在不同预设值0.2、0.25、0.33和0.5时，SNP检测结果有差异，选取坐标位置为28、100、254和1045的碱基峰图(图13)进行详细分析，分析表7样本LSYZ2017041中同一位置在不同预设值0.2、0.25、0.33和0.5时，SNP检测结果有差异，选取坐标位置为459、653、656、683、696和736的碱基峰图(图14)进行详细分析，图中左侧栏为预设值，后侧栏为5个不同坐标位置对应碱基的峰图结果，根据实际峰图结果判定，分析进行多态性位点分析时建议使用cut off值最低为0.25。

表6 LSYZ2017032样本不同SNP cutoff值检测结果

表7 LSYZ2017041样本不同SNP cutoff值检测结果

结果输出文件以输入文件名命名(图15-A)，将LSYZ2017032样本来源的六条序列拼接结果以fasta格式输出(图15-B)，其中单核苷酸多态性位点以简并碱基的形式记录，拼接结果中碱基质量值转换成图形信号以pdf格式文件输出(图15-D)，单核苷酸多态性位点以蓝色条柱标识，同时单核苷酸多态性位点识别结果以Excel文件格式输出(图16-F)，记录SNP位点在拼接序列中的坐标位置以及相对应的碱基；将LSYZ2017041样本来源的六条序列拼接结果以fasta格式输出(图15-C)，其中单核苷酸多态性位点以简并碱基的形式记录，拼接结果中碱基质量值转换成图形信号以pdf格式文件输出(图15-E)，单核苷酸多态性位点以蓝色条柱标识，同时单核苷酸多态性位点识别结果以Excel文件格式输出(图15-G)，记录SNP位点在拼接序列中的坐标位置以及相对应的碱基。

应该理解到披露的本发明不仅仅限于描述的特定的方法、方案和物质，因为这些均可变化。还应理解这里所用的术语仅仅是为了描述特定的实施方式方案的目的，而不是意欲限制本发明的范围，本发明的范围仅受限于所附的权利要求。

本领域的技术人员还将认识到，或者能够确认使用不超过常规实验，在本文中所述的本发明的具体的实施方案的许多等价物。这些等价物也包含在所附的权利要求中。