一种基于定向基因组测序获取全基因组SNP位点的方法

技术领域

本发明涉及遗传学领域，具体说是基于定向基因组测序技术获取全基因组SNP标记位点的方法。

背景技术

单核苷酸多态性(Single nucleotide polymorphisms，SNP)是目前最重要的分子遗传标记，因其在基因组中大量广泛分布，具有极高的检测丰度，且遗传稳定性较高，在分子遗传学中占据了重要地位。在人类和动植物遗传育种领域中，SNP被广泛用于复杂性状遗传机制的解析、基因组选择等研究方向。

不同的研究内容对遗传标记的数量需求有所不同，在遗传育种应用方面提高标记数目仍然是研究的重点之一。例如，基因组选择(genomic selection，GS)利用覆盖整个基因组的高密度SNP来构建亲缘关系矩阵，从而估计候选个体的基因组育种值并进行选种，因此需要通过高密度的遗传标记捕获到更多的加性遗传变异；而在全基因组关联分析中，想要更准确地鉴定到目标表型真正的致因突变，也需要提高标记密度。

当前主要的全基因组高通量分型技术有SNP芯片和基因组测序两种。SNP芯片的优势在于操作简单、分型速度快、准确性高，对分型流程的环境设施要求较低，可以在生产一线直接应用，因此芯片分型成为早期全基因组分型的主流方法。其不足之处在于，芯片设计完成后，能鉴定的SNP数据集是固定的，难以满足所有类型的研究需求，后期不易拓展升级，并且当面对多个育种群体时，部分SNP在有些群体中已经固定，因此同一款芯片难以覆盖所有群体的遗传特点，导致一些群体遗传信息损失。

测序成本不断降低使大规模的测序在动植物育种领域得到广泛使用。而定向基因组测序，意即将感兴趣的基因组区域通过探针杂交、酶切富集、多重PCR扩增等方法进行富集后进行测序的研究策略，根据不同的应用，利用较少的数据量就能得到超高的灵敏度和准确度，实现变异位点的快速筛选。目前定向基因组测序的方法有很多，例如简化基因组测序(Genotyping by Sequencing，GBS)、靶向基因组测序、多重PCR测序等。这些定向测序在不同程度上存在的共性问题是，基因组不同区域的富集效率可能存在较大差异。利用现有分析手段对该问题最常见的解决办法是，将低富集区域的测序数据进行删除，只保留高富集区域的测序数据用于基因组突变鉴定，这就直接损失了很多本应该富集到的基因组目标区段。为了在测序数据中能够获得更高比例的高富集区域，一般通过提升基因组整体测序深度来保证，然而这大幅提升了测序成本。

以GBS为例，这类方法通过富集并测序基因组中很小比例的片段得到分布全基因组的SNP标记，同时实验流程简便且成本较低。Tassel软件是基于简化基因组开发的程序包，专门用于GBS数据的SNP检测，但是对于双端测序而言，Tassel无法直接对双端reads识别，需要对每端的reads单独处理，处理后，要对两端的reads进行合并。并且Tassel在发现变异位点的过程中存在将插入缺失变异(Indel)误判成SNP的可能，因此在用传统的Tassel程序对GBS数据进行变异位点发现时存在假阳性。此外，传统GBS分析流程，对于每个位点靶向测序的测序深度存在最低要求(一般均要求10×/位点以上)，对于部分未达到测序深度要求的位点，则无法准确判定其杂合子/纯合子的状态，如果删除该位点则损失了标记密度，如果保留该位点则很可能降低分型的准确性，准确性较低的分型结果会进一步影响到填充准确性，不利于育种实践。

发明内容

针对现有技术中存在的缺陷，本发明的目的在于提供一种基于定向基因组测序获取全基因组SNP位点的方法。利用本方法产出的定向基因组测序数据，目的在于提供一种新的定向基因组测序基因分型方法，在维持低成本的同时，能够获得更高密度的SNP标记，并提高分型的准确率。

为达到以上目的，本发明采取的技术方案是：

一种基于定向基因组测序获取全基因组SNP位点的方法，其特征在于，包括如下步骤：

步骤1，提取育种群体中待检测个体样本的基因组DNA，对上述基因组DNA进行定向基因组测序得到定向基因组原始测序数据；

步骤2，将步骤1得到的定向基因组原始测序数据按照每个个体特定的barcode进行拆分，对拆分后个体的测序数据进行过滤，剔除不符合过滤标准的个体；

步骤3，将通过步骤2过滤筛选的个体的测序结果比对到参考基因组，对比对后的全部bam文件进行sort排序和构建索引；

步骤4，考察经步骤3对比后的全部bam文件所记载的基因组中每个位点是否为多态位点，得到各多态位点的位置信息；

步骤5，对步骤4得到的各多态位点进行进一步填充分型：通过EM迭代算法处理所述育种群体的突变位点信息，输出该突变位点的基因分型结果。

进一步，步骤2的具体步骤为：

步骤2-1，对步骤1得到的定向基因组原始测序数据进行质控修剪，具体为：以4bp碱基为一个窗口，滑动窗口计算碱基平均质量，碱基质量阈值为15，如窗口内平均质量值低于阈值，则删除该窗口及之后的所有碱基；控制序列长度，剔除修剪后长度低于75bp的测序Reads；

步骤2-2，将步骤1得到的定向基因组原始测序数据进行拆分，拆分时识别酶切位点序列的同时匹配每个个体特有的barcode；提取每条测序Reads的barcode信息后，删除Reads中的barcode序列；

步骤2-3，将经过步骤2-1质控修剪后的定向基因组原始测序数据，根据构建文库时的barcode序列记录进行个体拆分；

进一步，步骤4所述考察经步骤3对比后的产生的bam文件所记载的基因组中每个位点是否为多态位点，其考察过滤条件为：

剔除复等位基因、剔除EAF＜0.01的位点、保留pcr≥0.5的位点。

上述EAF为估计的群体最小等位基因频率。

上述保留pcr≥0.5的位点指：对每一个候选多态位点的群体覆盖度(populationcoverage rate，pcr)进行评估过滤，具体为在该位点有reads覆盖(即该位点的测序乘数≥1)的样本数目占总样本数目的比例须大于等于0.5。

进一步，步骤5中所述每个多态位点的信息包括Info Score值，分型结果包括对上述Info Score值进行过滤处理的结果，对上述Info Score值进行过滤的标准为：次等位基因频率大于0.05、单个SNP在群体中的call rate大于0.6、Info_Score大于0.4。

本发明所述的一种基于定向基因组测序获取全基因组SNP位点的方法，其有益效果为：

本发明适用于定向基因组测序，建立了一种分析定向基因组测序数据多态位点的新流程，该方法避免了原有定向基因组测序流程中对测序深度敏感、分型准确度有限等不利因素，在低成本的基础上，可以更准确地得到分布于全基因组的、更高密度的SNP信息。

此外，本过程得到的分型结果相较于原始定向基因组测序数据分型流程的结果，可以进一步向更高密度乃至全基因组水平填充，同时保持更高的填充准确性。

附图说明

本发明有如下附图：

图1本发明提供的基于简化基因组测序检测SNP标记位点的方法的流程图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明中的附图，对本发明中的技术方案进行详细说明，下面所描述的实施例是本发明一部分实施例，只为起到说明的目的，而不是为了限制全部的适用范围。基于本发明中的实施例，本领域技术人员在没有背离本发明的前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例1

实施例1以GBS为例对本发明所述基于定向测序获取全基因组SNP位点的方法进行具体说明。

1、实验群体：

本实施例采用以快大速长型岭南黄鸡(High Quality Line A，HQLA)和广东惠阳胡须鸡(Huiyang Beard chicken，HB)为亲本的正反交深度杂交系为实验群体，F2代之前采用正反交建系方案，F2代后的建系，采取公母配比为1∶2的伪随机交配，每代总群体数目保持在1000只以上(F3代除外)，用于繁育后代配种的公鸡保持在80只以上，母鸡保持在160只以上。维持此配种方案将该群体传至F16代。

2、实验数据：

本实施示例选用上述群体F0代、F9代的简化基因组原始测序数据(GBS)。

3、原始数据质控

控制序列的质量值，以4bp碱基为一个窗口，滑动窗口计算碱基平均质量，碱基质量阈值为15，如窗口内平均质量值低于阈值，则删除该窗口及之后的所有碱基；控制序列长度，剔除修剪后长度低于75bp的测序reads。

4、数据拆分

对F0、F9代共807个个体的GBS原始测序数据进行拆分，拆分时识别酶切位点序列的同时匹配每个个体特有的barcode，其中单端测序的barcode位于EcoR I酶切位点处，长度为6-9bp不等。在提取每条测序Reads的barcode信息后，删除Reads中的barcode序列，将上一步质控修剪完成的GBS数据，根据构建文库时的barcode序列记录进行个体拆分，拆分后得到每个个体的FASTQ格式数据，当同时满足特定的酶切位点及barcode与个体号对应时可拆分成功。

5、比对参考基因组

本实施例中采用的比对软件为BWA，在比对过程总采用的参考基因组版本为鸡Gallus_gallus-6.0版本，对比对后的bam文件进行sort排序和构建索引。

6、变异位点的鉴定以及过滤

将比对至参考基因组后的GBS测序数据后的807个bam文件，使用Basevar(0.0.1.3)软件进行判定，判定内容为基因组每一个位点是否存在多态，共检测出1801326个SNP多态位点，然后根据以下过滤条件：①剔除复等位基因；②剔除EAF＜0.01的位点；③保留pcr≥0.5的位点。

按照以上条件过滤后，剩余411K位点。

7、基因分型以及质控

利用STITCH软件进行进一步的填充分型，分型的结果还包含的每个位点Info_Score值表明了该位点的的SNP分型质量，在进一步对SNP进行质控过滤时也需要将此参数考虑进去，进一步对Info_Score进行过滤，过滤标准如下：①次等位基因频率大于0.05；②单个SNP在群体中的call rate大于0.6；③Info_Score大于0.4。

经过以上过滤标准后，剩余273K位点，即本方法得到的最终结果。

对比例：

本实施例选用F0代全基因组重测序数据、F0、F9代GBS测序原始数据和Tassel处理后的基因型数据。

其中全基因组重测序数据主要采用HiseqX Ten测序平台双端151bp测序获得，平均测序深度为10×，比对软件为BWA，主要用到的变异鉴定软件为GATK，采用Hard Filter方法过滤后去除MAF＝0的位点后剩余12M位点的基因型数据(下文简称“GATK数据”)。该高深度重测序数据的分型结果作为本实施示例的参考数据集，即真值数据。

GBS原始数据以EcoRI-MseI双酶切建库的方法，通过单端91bp获得，测序平台为Nextseq500。使用Tassel处理GBS原始数据，作为本实施示例的对照数据集。Tassel处理后的GBS原始SNP数据(以下简称GBS-Tassel数据)按照如下参数对其进行过滤，具体的过滤参数为：去除次等位基因频率MAF小于0.05的SNP位点，去除测序质量值小于98的位点，去除测序深度小于5的位点，然后采用Beagle软件进行基因型填充，共得到189K位点(下文简称“Tassel数据”)。

实施例1和对比例所得多态位点的正确率及分型准确性的对比：

以高深度GATK数据作为金标准，其中分型准确性用F0代30个个体的三种基因型数据判断，位点正确率用与高深度位点的重合个数与全部位点数的比值表示，分型准确性用同一个个体的不同分型方式得到的基因型数据与高深度数据的基因型一致性表示。

Tassel数据的189K位点与高深度GATK数据比较分型准确性，与高深度GATK数据的基因型一致性平均为93.6％。

本发明得到的新方法数据得到的273K位点与高深度GATK数据比较分型准确性，与高深度GATK数据的基因型一致性平均达到98.5％。

本实施例证明，采用本发明中的流程对GBS数据进行变异位点的发掘，可以有效提高SNP位点的检出数目，并且提高分型的准确性。

虽然，上文中已经用一般性说明及具体实施方案对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。本说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。