CN103088120B - 基于SLAFseq技术的大规模样品基因分型方法 - Google Patents

Abstract

本发明提供一种基于SLAF-seq技术进行大规模样品基因分型的方法，其是利用SLAF-seq技术降低基因组的复杂度，并对大规模样品进行基因分型。利用该技术对基因组进行高通量测序，对样品进行标记开发、遗传图谱绘制及全基因组关联分析。该方法与传统的方法相比优点在于通量大幅度提高，成本大幅度降低。该方法主要应用于标记开发、遗传图谱绘制及全基因组关联分析。

Description

基于SLAFseq技术的大规模样品基因分型方法

技术领域

本发明提供一种基于SLAF-seq技术的大规模样品基因分型方法，其核心技术是利用SLAF-seq技术降低基因组复杂度并进行高通量测序，进行标记开发、遗传图谱绘制及全基因组关联分析。 

背景技术

随着具有高通量、低成本、测序错误率低、测序读长短特点的新一代测序技术和生物信息学的发展，使得该测序技术进行高通量标记开发成为可能。SLAF-seq（Specific Length Amplified Fragments sequencing）是一种简化基因组深度测序技术，在高通量测序技术的基础上，利用生物信息学方法对目标物种的参考基因组或已知BAC序列进行系统分析，根据基因组的GC含量、重复序列情况和基因组特点等信息，设计标记开发方案，以保证SLAF标签密度、均匀性、效率和分析的准确性。 

根据基因组特性，利用SLAF-seq技术对大规模样品进行研究，通过数学算法解决传统方法通量低和准确性不高的问题，提高生物学分析的准确性，降低成本，提高效率。目前还未见利用SLAF-seq技术对大规模样品进行基因分型的报道。 

发明内容

本发明的目的是利用SLAF-seq技术降低基因组复杂度并进行高通量测序的方法，对大规模样品进行研究，通过数学算法解决传统方法通量低和准确性不高的问题，提高生物学分析的准确性，降低成本，提高效率。 

为了实现本发明目的，本发明提供基于测序和SLAFseq技术的大规模样品基因分型方法，其基于SLAFseq技术降低基因组复杂度并进行高通量测序，对大规模样品进行标记开发、遗传图谱图谱、单体型图谱绘制及性状关联分析，包括如下步骤： 

1）对经性状鉴定的各样品的DNA进行检测； 

2）对样品基因组的复杂性进行降低处理，获得复杂性降低的DNA样品； 

3）利用二级引物库中的引物对复杂性降低后的样品DNA进行PCR扩增，使特异性长度片段扩增前后丰度一致； 

4）将扩增后的特异性长度片段连接标准测序接头，利用高通量测序技术进行测序； 

5）对于各样品的测序结果进行比较分析，获得SLAF标记，获得各样品的有效序列；进行遗传图谱绘制或全基因组关联分析。 

其中步骤1）通过琼脂糖电泳检测各样品DNA主带是否清晰，有无降解和污染，通过微量分光光度计如Nanodrop 2000检测DNA的浓度及纯度。 

其中步骤5）是对上述SLAF标签进行多态性分析，对存在多态性的标签中的所有样品进行基因分型，并使用自主研发打分体系进行打分，设置阈值，最终获得SLAF标记；通过上述步骤获得的SLAF标记，可用于遗传图谱绘制、全基因组关联分析。 

前述的基于测序和SLAFseq技术的大规模样品基因分型方法，其中所述大规模样品是植物、动物或微生物。 

前述的基于测序和SLAFseq技术的大规模样品基因分型方法，其中所述大规模样品为自然群体、遗传分离群体；所述遗传分离群体为经鉴定的性状分离群体，包括F2、BC1、DH等目标性状分离的群体。 

DNA复杂性降低方法是任何有选择性的减少基因组复杂度的方法，包括限制性内切酶酶切产物PCR扩增或酶切产物选择性吸附。 

其中最佳酶切方案需满足以下条件：a）保证序列标签在基因组上分布尽可能均匀；b）选择特定长度的酶切片段能够保证序列标签的数量；c）选择特定长度的酶切片段避免落在基因组高度重复区。 

其中最佳酶切方案确定需进行预实验，步骤及满足条件如下：1）通过生物信息学模拟，获得1~3套候选方案；2）对样品进行酶切连接、PCR扩 增及琼脂糖凝胶电泳；3）电泳结果各样品目标范围内无特异条带，且通过BMP图片解析软件如bmp2txt获得电泳胶图上各位点灰度值，利用灰度值模拟模板量进行扩增倍率分析，模板范围扩增倍率一致性高者即为最佳方案。 

前述的基于测序和SLAFseq技术的大规模样品基因分型方法，其中所述步骤3）二级引物库是为了利用高通量测序特点，节约测序成本，特为大样本量项目设计的扩增引物，后文中称二级引物，所述二级引物差异序列长度为3~7bp，通过生物信息学方法对ATCG组合进行相似性评估，剔除相似度高序列，保障组合内序列完全识别。所述的二级引物库中的引物序列如SEQ ID NO 4~112所示。 

前述的基于测序和SLAFseq技术的大规模样品基因分型方法，其中所述二级引物结合Solexa高通量测序标准引物，测序样品数量将以乘积形式增加，如：Solexa引物用12个，二级引物用96个，即可实现1次1152个样品的测序和数据分类需求； 

前述的基于测序和SLAFseq技术的大规模样品基因分型方法，其中所述步骤4）中Solexa技术是一种基于边合成边测序技术（Sequencing-By-Synthesis，SBS）的新型测序方法。通过利用单分子阵列实现在小型芯片（FlowCell）上进行桥式PCR反应。由于新的可逆阻断技术可以实现每次只合成一个碱基，并标记荧光基团，再利用相应的激光激发荧光基团，捕获激发光，从而读取碱基信息。 

前述的基于测序和SLAFseq技术的大规模样品基因分型方法，其中所述步骤4）中双端测序的优势如下：a）相对于单端测序，双端测序序列有效长度增倍，标记比对、定位更为准确；b）二级引物序列测序重复1次，提高碱基准确性，碱基错误率由0.01降低到0.0001。 

前述步骤4）的标准测序接头为： 

5'-GATCGGAAGAGCACACGTCT-3' 

序和SLAFseq技术的大规模样品基因分型方法，其中打分体系对SLAF标记分型结果给出基于SNP和深度观测的条件错误率，错误率的计算使用贝叶斯公式，具体方法如下：1）假设某一标记位点有m个等位基因，记为{A1,A2,…Am}，则对特定个体而言，在该位点所有可能的分型种类有m×(m+1)/2种，每种基因型的先验概率为该基因型出现的理论频率，因此纯合基因型先验概率为1/m^2,杂合基因型先验概率为2/m^2;对于二倍体物种，单次测序杂合基因型两个等位基因Ai,Aj间因测序错误产生交换的概率计算如下： 

$P_{\mathrm{ij}}^e=\frac{{(1-\mathrm{\ϵ})}^{{(\mathrm{rl}-s)}_{\mathrm{\ϵ}}s}}{3}$

其中：ε为测序的平均错误率，模型中取0.015，rl为测序序列长度，s为snp数目。基于该错误概率可以计算每种分型下单次测序测得等位基因Ai的概率P(A_i|A_uA_v)；2）对每一特定的等位基因Ai，假设其深度观测为di;则每种可能分型下产生该种深度观测的条件概率计算如下； 

$P\left((d_1,d_2,...,d_m)\right|A_i{A}_j)=\frac{\left({\mathrm{\Σ}}_{u=1}^{u=m}{d}_u\right)!}{{\mathrm{\Π}}_{u=1}^{u=m}{d}_u!}\×{\mathrm{\Π}}_{v=1}^{v=m}P{\left(A_v\right|A_i{A}_j)}^{d_v}$

3）最后由贝叶斯公式计算相应深度观测下，产生当前分型的概率： 

$P\left(A_i{A}_j\right|(d_1,d_2,...,d_m))$

$=P\left((d_1,d_2...,d_m)\right|A_i{A}_j)$

$\×\frac{P_{A_i{A}_j}}{{\mathrm{\Σ}}_{s=1}^{s=m}{\mathrm{\Σ}}_{t=1}^{t=m}(P\left((d_1,d_2,...,d_m)\right|A_s{A}_t)\×P_{A_s{A}_t})}$

 4）通过Phead score公式将该概率转换为最终的score，作为对该种分型错误率的度量： 

Score=-10×log₁₀(1-P(A_iA_j|(d₁，d₂，…，d_m)))

5）对杂合基因型一个等位基因深度为1，另一等位基因深度小于等于5的错误率直接使用公式： 

$P_{\mathrm{ij}}^e=\frac{{(1-\mathrm{\ϵ})}^{{(\mathrm{rl}-s)}_{\mathrm{\ϵ}}s}}{3}$ 计算。Score的上界限定为30。 

本发明所述的基于测序和SLAFseq技术的大规模样品基因分型方法，比传统方法的优点在于1）通过生物信息学方法可对酶切组合及特异长度片段数量进行预测，结合预实验选择最佳方案，实现SLAF标签数量可控、分布可控，并有效避开重复序列区，大大简化基因组，节省大量成本；2）采用第二代Solexa高通量双端测序方法，标签有效长度倍增，双端引物序列测序两次，单碱基错误率由0.01降至0.0001，准确性大幅提升；3）结合Solexa标准引物，设计足量二级引物可一次满足大规模样品测序及数据分类需求，大大减少建库及测序成本；4）大规模样品基因分型并使用自主研发打分体系保障分型准确及SLAF标签质量。 

附图说明

图1A为模拟AluI+MseI酶切片段在斑马鱼各染色体上的密度分布（a）和数量分布（b）。图1B为鲤鱼样品AluI+MseI酶切连接产物电泳图，其中1-20为群体，父本、母本：为亲本；SD：水稻为对照；M:100bp DNA ladder。图1C为扩增倍率分布图。 

图2A为鲤鱼样品AluI+MseI酶切连接PCR产物电泳图，其中，43-45,46-61为群体，父本、母本：为亲本；水稻：为对照；M:100bp DNAladder。图2B为鲤鱼样品DNA扩增产物纯化后混样电泳图，混样1为群体1-45的DNA混样，混样2为群体46-90的DNA混样，混样3为群体91-135的DNA混样，混样4为群体136-180的DNA混样，混样5为群体181-224的DNA混样。 

图3A为鲤鱼混合样品加入Solexa标准测序接头切胶前电泳图，混样1为群体1-45的DNA混样，混样2为群体46-90的DNA混样，混样3为群体91-135的DNA混样，混样4为群体136-180的DNA混样，混样5为群体181-224的DNA混样。图3B为鲤鱼混合样品加入Solexa标准测序接头切胶后电泳图。图3C为鲤鱼样品被切胶选择片段回收后电泳图。 

图4为SLAF标记分型数据格式图，其中MarkerID为SLAF标签编号，type为SLAF标签类型，为遗传图中的分型类型，aa-av为群体编号，对应的列表示群体所在的SLAF中的分型类型。 

图5为SLAF标记基因型类型数量分布图，横坐标表示SLAF标记在遗传图中的分型类型，其中1为abxcd，2为efxeg，3为hkxhk，4为lmxll，5为nnxnp，纵坐标表示SLAF标记的个数。 

图6为遗传图谱示意图，表示构建的遗传图中的四个连锁群的示意。 

图7为大豆样品SLAF标签及基因组重复序列的染色体分布图。 

图8为大豆样品SLAF标签在各染色体上的数量分布图。 

图9为大豆样品SNP位点在各染色体上的数量分布图。 

图10为大豆样品单体型连接示意图，R2表示染色体位置上完全连锁平衡（R2=0）到完全连锁不平衡（R2=1）的过程。 

图11为大豆样品下胚轴颜色性状全基因组关联分析结果示意图，横坐标1-20表示大豆的染色体编号，纵坐标表示关联强度。 

具体实施方式

以下实施例用于说明本发明，但不用来限制本发明的范围。 

实施例1基于测序和SLAFseq技术的大规模样品基因分型方法应用于鲤鱼224株F1群体遗传图谱构建 

鲤鱼遗传图谱构建，包括如下步骤： 

1、鲤鱼群体（购自黑龙江水产研究所）的选择：群体为F1群体，群体数量为224株，亲本2株，开发50,000标签。 

2、如表1所示为利用生物信息学方法，对相近物种斑马鱼的基因组序列进行酶切预测，基因组版本号为Ensembl公布的Zv9，所选最佳方案AluI+MseI的各参数评估结果，结果显示特异长度为330-380的酶切片段数量为64,918个，满足标签开发需求，同时低拷贝标签比率为94.97%，有效避免高度重复序列；如图1A所示为该酶切方案所得标签在基因组各染色体上数量分布评估图，结果显示酶切片段分布均匀；如图1B、1C所示为该酶切方 案预实验结果及扩增倍率分析结果，结果显示目标范围内特异长度DNA片段丰度一致性高。 

表1酶切方案AluI+MseI的各参数评估结果 

注：特异性标签比率：拷贝数为1的标签比率；低拷贝标签比率：拷贝数<=10的标签比率 

3、如图2A所示，为AluI+MseI双酶切对照样品、鲤鱼亲本和群体DNA后加接头并进行扩增后的电泳图；连接所用接头为：5'-GACGATGAGTCCTGAGTACTCAGGACTCAT-3'；扩增所用二级引物为二级引物库中的48个引物，分别为序列表中的第4，7，8，12，13，16，17，19，20-23，26，27，31-34，38，42-45，47-50，54，55，59，61，62，64，71，72，75，79，80，87，104-112。所选用的部分引物见表2，表2中混池列表示每一个混池所混合的群体编号，表2中下划线的碱基代表差异序列，表2中在五个混池中样品1，46，91，136，181使用的二级引物的编号为SEQID No.19。PCR扩增条件如下：98℃预变性3分钟；接着18个循环98℃变性10秒，65℃退火30秒，72℃延伸30秒；最后72℃延伸5分钟。随机选择二级引物，二级引物的数量由群体数和测序量决定，如本项目共224个群体，2个亲本，1个对照，根据测序量计算可以放在5个库进行测序，平均每个库中放45个群体，每个库中加入亲本和对照共使用48个二级引物以区分不同样品，每个库可使用不同测序引物以区分。 

表2扩增所用二级引物部分序列 

4、如图2B所示，为将上述对照、亲本、群体DNA扩增产物根据表2分 组，按3:4:5混合后DNA电泳图；M为100bp DNA ladder，混样1-5为5个DNA混池。 

5、如图3A、3B所示，分别为样品混1~混5加入Solexa标准测序接头（5'-GATCGGAAGAGCACACGTCT-3'和5'ACACTCTTTCCCTACACGACGCTCTTCCGATCT-3'）后在切胶前、切胶后的电泳图；图3C所示为样品混1~混5的被切胶产物的再次电泳图；切胶范围为450-500bp，FI909-1~FI909-5为450-500bp特异长度条带；A、B中M为50bpDNA ladder，C中M为100bp DNA ladder。 

6、将切胶得到的450-500bp的特异长度DNA片段进行Solexa双端测序。 

7、对原始测序数据进行数据识别处理，根据各样品连接的二级引物和测序引物序列的唯一性，得到各样品有效序列，经聚类分析及纠错，获得SLAF标签50,457个；如表3所示为SLAF标签多态性分类结果；如表4所示为各样品测到的标记数、深度及平均深度统计结果。 

表3总测序标签数量及测序深度统计 

表4每个样品测序标签数量及测序深度统计 

8、如表5所示为F1群体标记编码规则，选择上述多态性SLAF标记，依据表5规则进行基因型编码，完成226个样品在11,472个SLAF标记中的分型判断，共计2,592,672次；过滤掉父母本缺失及不适合F1群体作图的标记，得到有效SLAF标记9,017个，如图4所示为标记分型数据格式示意图，如图5所示为标记基因型类型分布图。 

表5多态性标记类型 

注：--表示子代基因型缺失 

9、对上述9,017个有效标记的分型结果，通过自主研发打分体系给出基于SNP和深度观测的条件错误率并进行打分，过滤掉分数低于阈值的完整度不高的标记和个体，最后选出高质量的SLAF标记4,051个，保留个体135个。 

10、对上述标记做两点连锁分析，计算重组率和LOD值，LOD值满足3的标记采取最短距离法进行聚类分析，将4,026个标记划分为50个连锁群；以连锁群为单位，通过Joinmap4.0获得连锁群内标记的线性排列，通过多点分析估算相邻标记间的遗传距离，最终绘制到图谱上的标记数量为4,026， 总图距为4,319厘摩；如图6所示为图谱结果示意图。 

通过对酶切方案及特异长度片段数量进行预测，保证特异性标签比率达到82%以上，有效避开重复序列区域，节省测序成本；二级引物库与Solexa测序引物配合，降低了大规模样品（本例为226个样品）基因分型的建库成本；自主研发的打分体系保障了分型准确，在成本降低的同时，依然完成了4,026个高质量标记的图谱绘制。 

实施例2基于测序和SLAFseq技术的大规模样品基因分型方法应用于大豆512个自然群体的单体型图谱构建和性状关联分析 

大豆单体型图谱构建和性状关联分析，包括如下步骤： 

1、大豆群体的选择：大豆自然群体共512份种质资源，设计50,000个SLAF标签。 

2、如图7所示，为利用生物信息学对大豆基因组进行酶切模拟并结合预实验确定的最佳方案HaeIII+MseI的避开高度重复区效果评估图，结果显示该方案酶切片段在基因组上均匀分布且有效降低基因组重复序列比例；如表6所示，为该方案简化基因组效果评估统计，结果显示基因组简化率为0.49%，大大减少测序量节约项目成本。 

表6简化基因组效果评估 

3、对于对照样品1个和大豆样品512个进行HaeIII+MseI双酶切，酶切温度37℃，酶切时间15小时，然后加接头，接头序列为：5'-GACGATGAGTCCTGAGTACTCAGGACTCAT-3'。连接产物使用二级引 物库中的引物进行PCR扩增、电泳，PCR扩增条件如下：98℃预变性3分钟；接着18个循环98℃变性10秒，65℃退火30秒，72℃延伸30秒；最后72℃延伸5分钟，此项目共使用42个二级引物，分别为序列表中的第4，7，8，11-13，16，17，19-23，26，27，31-34，38，39，42-45，48-50，54，55，59，61，62，64，71，72，75，78-80，87，104条引物。 

4、将对照样品和大豆样品以单株3：5混合后进行电泳。 

5、对上述混样连接Solexa标准测序接头后进行电泳，接头序列如下：5'- 

GATCGGAAGAGCACACGTCT-3'和 

5'-ACACTCTTTCCCTACACGACGCTCTTCCGATCT-3' 

6、对上述电泳结果进行切胶，切胶范围为500-550bp。 

7、将截取的500-550bp的特异长度DNA片段进行Solexa双端测序。 

8、对原始测序数据进行数据识别处理，根据各样品连接的二级引物和测序引物序列的唯一性，得到各样品有效序列，经比对和聚类纠错，最终定位得到SLAF标签59,494个；如图8所示，为深度最深的5万个SLAF标签在全基因组上的分布情况，结果显示实际所得SLAF标签在基因组上分布均匀。 

9、根据上述5万SLAF标签的分型结果，在全基因组上初步检测到SNP位点共76,001个。如表9所示为所有SNP位点完整度和杂合率等情况统计，结果显示SNP完整度达到了82%，整体杂合率15%，在基因组上SNP密度达到平均每12.50K含1个SNP位点；如图9所示为SNP位点在基因组上的分布图，整体上SNP在基因组上均有分布，在特异区域分布密集说明该区域突变比较活跃。 

10、利用上述得到的具有多态性SNP的Marker，通过单体型分析软件haploview，根据连锁不平衡的原理，然后选取频率高的相邻单体型进行动态连接，连接结果示意图见图10。 

11、利用SNP分型数据，使用MLM方法用Q+K的混合线性模型，实现全基因组关联分析。用Structure软件得到群体结构Q，用SPAGeDi软件得到 亲缘关系K，再结合基因型数据和表型数据，用TASSEL软件进行comparedMLM分析，得到经过Bonferroni校正的关联结果。如图11为下胚轴颜色性状关联结果，染色体12-13之间位置，关联强度超过域值，即为检测到的性状相关区域。 

通过对酶切方案及特异长度片段数量进行预测，重复序列区域SLAF标签占总标签数量的25.32%，低于大豆基因组42.96%的重复序列比例，有效避开重复序列区域，节省测序成本；二级引物库与Solexa测序引物配合，降低了大规模样品（本例为512个样品）基因分型的建库成本；在成本降低的同时，依然获得了76,001个高质量标记的SNP标记。 

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

Claims (7)

1.基于SLAFseq的大规模样品基因分型方法，其基于SLAFseq技术降低基因组复杂度并进行高通量测序，对大规模样品进行标记开发及遗传图谱绘制、全基因组关联分析，包括如下步骤：

1)对经性状鉴定的各样品的DNA进行检测，通过琼脂糖电泳检测各样品DNA主带是否清晰，有无降解和污染，通过微量分光光度计检测DNA的浓度及纯度；

2)对样品基因组的复杂性进行降低处理，获得复杂性降低的DNA样品；

3)利用二级引物库中的引物对复杂性降低后的样品DNA进行PCR扩增，使特异性长度片段扩增前后丰度一致；二级引物库序列如SEQ ID NO.4～112所示；

4)将扩增后的特异性长度片段连接标准测序接头，利用Solexa双端测序方式进行测序；标准测序接头为：

5′-GATCGGAAGAGCACACGTCT-3′

5′-ACACTCTTTCCCTACACGACGCTCTTCCGATCT-3′；

5)比较测序结果，对上述SLAF标签进行多态性分析，对存在多态性的标签中的所有样品进行基因分型，并进行打分，设置阈值，最终获得SLAF标记，获得SLAF标记后，进行遗传图谱绘制或全基因组关联分析；

通过上述步骤获得的SLAF标记，用于遗传图谱绘制、全基因组关联分析；

其中，所述大规模样品来自动物、植物或微生物；

所述SLAFseq是一种简化基因组深度测序技术，是在高通量测序技术的基础上，利用生物信息学方法对目标物种的参考基因组或已知BAC序列进行系统分析，根据基因组的GC含量、重复序列情况和基因组特点等信息，设计标记开发方案，以保证SLAF标签密度、均匀性、效率和分析的准确性；

所述样品为自然群体或遗传分离群体；当所述样品为遗传分离群体时，所述遗传分离群体为经鉴定的性状分离群体。

2.如权利要求1所述的方法，其特征在于，所述样品为遗传分离群体；所述遗传分离群体为选自F2、BC1、DH目标性状分离的经鉴定的性状分离群体。

3.如权利要求1所述的方法，其特征在于，步骤2)中DNA复杂性降低方法是任何有选择性的减少基因组复杂度的方法，包括限制性内切酶酶切产物PCR扩增或酶切产物选择性吸附。

4.如权利要求3所述的方法，其特征在于，酶切方法满足以下条件：

a)保证序列标签在基因组上分布均匀；

b)选择特定长度的酶切片段能够保证序列标签的数量；

c)选择特定长度的酶切片段避免落在基因组高度重复区。

5.如权利要求3所述的方法，其特征在于，酶切方法的确定需进行预实验，包括如下步骤：

1)通过生物信息学模拟，获得1～3套候选方案；

2)对样品进行酶切、二级引物连接、PCR扩增及琼脂糖凝胶电泳；

3)电泳结果各样品目标范围内无特异条带，且通过图片解析软件获得电泳胶图上各位点灰度值，利用灰度值模拟模板量进行扩增倍率分析，模板范围扩增倍率一致性高者即为可用的酶切方法。

6.如权利要求1中所述的方法，其特征在于，所述步骤5)的比较测序结果包括特异性长度片段多态性和单核苷酸多态性分子标记。

7.权利要求1-6任一所述的方法在遗传图谱绘制或全基因组关联分析中的应用。