一种测试纯系大豆新品种的特异性、一致性与稳定性的方法

摘要

本发明公开了一种测试纯系大豆新品种的特异性、一致性与稳定性的方法。所述方法包括：获得变异位点；确定待测大豆品种的测试区域；构建数据库；确定抽样量后，随机抽样混合并提取混合样本的DNA；制备引物；利用引物对混合样本的DNA进行扩增，扩增产物用于构建高通量测序文库；对高通量测序文库进行高通量测序，得到测序片段组；分析测序片段组，获得待测大豆品种基因型和杂株基因型；比较获得近似品种、变异位点和变异位点率；将杂株基因型与数据库中的基因型比较，获得杂株品种后，计算杂株率；利用变异位点、变异位点率和杂株率，判断待测大豆品种特异性、一致性和稳定性。该方法能够准确、完整地判断待测大豆品种的特异性、稳定性与一致性。

说明书

技术领域

本发明涉及生物技术领域，特别涉及一种测试纯系大豆新品种的特异性、 一致性与稳定性的方法。

背景技术

作为一种特化的知识产权，植物新品种已经成为一个公司及至一个国家的 核心竞争力。植物新品种授权与相关法律问题的解决依赖于DUS测试，即对待 测大豆品种的特异性(Distinctness)、一致性(Uniformity)和稳定性(Stability) 的田间种植鉴定或室内分子标记鉴定。田间种植鉴定流程为：将待测大豆品种 与近似品种同时植于田间，在2年及以上的生长季节内，观察它们的多个性状， 根据性状表现判断待测大豆品种与近似品种的差异显著性，即特异性，同时判 断群体内杂株比例，即一致性和稳定性；室内分子标记鉴定的流程为：分单株 提取待测大豆品种与近似品种中每个样本的DNA，并分别对每个样本的每个测 试区域进行PCR(Polymerase Chain Reaction，聚合酶链反应)，并对每个PCR 产物进行电泳或一代测序检测，根据检测结果，获得待测大豆品种与近似品种 的差异位点比例，根据差异位点比例，判断待测大豆品种的特异性。

田间种植鉴定的缺点是：周期长、工作量大，环境影响性状，导致判断不 准确。室内分子标记鉴定的缺点是：需要分别处理每个样本的每个测试区域， 工作量大，不能对样本与测试区域大量抽样，无法计算杂株率，因而无法进行 稳定性与一致性的测试。田间种植鉴定与室内分子标记鉴定的共同缺点是：均 由于工作量的原因，无法从现有品种中客观选择近似品种，只能由品种权申请 人提供，而基于商业利益等动机，品种权申请人提供的近似品种可能不真实， 从而造成错误品种授权的法律后果。

发明内容

为了解决现有技术中的问题，本发明实施例提供了一种测试纯系大豆新品 种的特异性、一致性与稳定性的方法。所述技术方案如下：

本发明实施例提供了一种测试纯系大豆新品种的特异性、一致性与稳定性 的方法，所述方法包括：

获得待测大豆品种属于同种内不同品种间的变异位点；

通过所述变异位点确定待测大豆品种的测试区域，所述测试区域包括通用 测试区域，至少部分所述变异位点包含在所述通用测试区域内；

构建包含所述不同大豆品种在所有所述测试区域的基因型的数据库；

确定所述待测大豆品种的抽样量SN后，随机抽样混合并提取混合样本的 DNA；

制备扩增所述测试区域的引物，所引物包括所述通用测试区域引物；

利用所述引物对所述混合样本的DNA进行扩增，得到所述测试区域的扩增 产物，所述扩增产物作为高通量测序文库；

对所述高通量测序文库进行高通量测序，得到测序片段组；

分析所述测序片段组，获得待测大豆品种基因型和杂株基因型；

将所述待测大豆品种基因型与所述数据库中的所述不同品种的基因型比 较，获得所述待测大豆品种的近似品种、变异位点和变异位点率；

将所述杂株基因型与所述数据库中的所述不同品种的基因型比较，获得杂 株品种后，计算杂株率；

利用所述变异位点、所述变异位点率和所述杂株率，判断所述待测大豆品 种的特异性、一致性和稳定性。

具体地，所述抽样量SN满足如下条件：BINOM.INV(SN,M,0.95)/SN ≤1.15*M，其中BINOM.INV为excel 2010中的函数，M为判断所述一致性和稳 定性时所选用的阈值，所述抽样量SN满足的条件含义为：即使所述杂株率只超 出一致性和稳定性时的判断阈值M的15％，所述抽样量在95％的概率保证下， 可正确判断所述待测大豆品种的稳定性与一致性。

具体地，所述高通量测序的深度CF满足如下条件：BINOM.DIST(10,10, BINOM.DIST(8,20,BINOM.DIST(0,CF,0.1％,TRUE),TRUE),FALSE)≥99.9％，1-BI NOM.DIST(10000,10000,1-BINOM.DIST(8,20,1-BINOM.DIST(99.99％*CF,CF,99. 9989％,TRUE),TRUE),FALSE)≤0.1％且BINOM.DIST(10*(1-M)*CF,10*CF,1-11 0％*M,TRUE)≥95.0％，其中，CF为所述高通量测序的深度，M为判断所述一致 性和稳定性时所选用的阈值，BINOM.DIST为excel 2010中的函数，所述高通 量测序的深度CF满足的条件含义为：在所述杂株率低至0.1％、所述杂株品种 为10个且所述杂株品种与所述大豆品种间平均仅有20个差异位点的条件下， 由所述高通量测序的深度CF决定的检出全部所述杂株品种的概率≥99.9％；在所 述数据库的品种为10000个且所述杂株品种与所述大豆品种间平均仅有20个差 异位点的条件下，由所述高通量测序的深度CF决定的存在误判所述杂株品种的 概率≤0.1％；在所述杂株品种为10个且真实杂株率仅超过判断特异性时所选用 的阈值的10％时，由所述高通量测序的深度CF决定的对稳定性与一致性的判定 结论正确的概率≥95.0％。

具体地，所述测试区域还包括非通用测试区域，所述引物还包括非通用测 试区域引物。

进一步地，所述非通用测试区域引物包括第一引物和第二引物，所述第一 引物包括第一正向引物和第一反向引物，所述第二引物包括第二正向引物和第 二反向引物，所述第一引物和所述第二引物分别进行单独扩增得到两个所述非 通用测试区域的扩增产物，将两个所述非通用测试区域的扩增产物等量混合用 于构建单独扩增的高通量测序文库；

所述第一正向引物的5’端连接有如序列表中SEQ ID NO:1所示的序列1， 所述第一反向引物中的5’端连接有如序列表中SEQ ID NO:2所示的序列2；

所述第二正向引物的5’端连接有如序列表中SEQ ID NO:2所示的序列2， 所述第二反向引物的5’端连接有如序列表中SEQ ID NO:1所示的序列1。

具体地，利用所述变异位点、所述变异位点率和所述杂株率，判断所述待 测大豆品种特异性、一致性和稳定性的方法包括：

当所述变异位点率≥SD或所述非通用测试区域存在所述变异位点时，所述 待测大豆品种具有特异性，当所述变异位点率＜SD且所述变异位点不存在于所 述非通用测试区域中时，所述待测大豆品种不具有特异性，其中，SD为判断特 异性时所选用的阈值；

当所述待测大豆品种的所述杂株率≤M时，所述待测大豆品种具有一致性和 稳定性，当所述待测大豆品种的所述杂株率大于＞M时，所述待测大豆品种不 具有一致性和稳定性，M为判断所述一致性和稳定性时所选用的阈值；

所述杂株率R＝R1+R2-R3-R4，其中：

$R1={\Sigma }_{i1=1}^{n1}\frac{{\Sigma }_{j1=\mathrm{Int}(0.8\times t1)+1}^{t1-\mathrm{Int}(0.1\times t1)}2\times R1i1j1}{t1-\mathrm{Int}(0.8\times t1)-\mathrm{Int}(0.1\times t1)},$其中，n1为细胞核杂株品种的数目， t1为第i1个所述细胞核杂株品种的所有特异杂株核基因型的数目，i1j1为第i1 个所述细胞核杂株品种的所有所述特异杂株核基因型按频率由低到高排序后， 第j1个所述特异杂株核基因型，R1i1j1为第i1j1个所述特异杂株核基因型的频 率；R1为由杂株核基因型计算的所述细胞核杂株品种的所述杂株率的总和，所 述细胞核杂株品种的所述杂株率为去掉所述细胞核杂株品种中最低的80％和最 高的10％的所述特异杂株核基因型的频率后，剩余的所述特异杂株核基因型的 频率的平均值的2倍；

其中，t2为除所述细胞核杂株品种拥有的杂 株核基因型之外的且频率≥0.17％的所述杂株核基因型的数目，i2为除所述细胞 核杂株品种拥有的所述杂株核基因型之外的所有所述杂株核基因型按频率由低 到高排序后，第i2个所述杂株核基因型，R2i2为第i2个所述杂株核基因型的频 率；R2是利用除所述细胞核杂株品种拥有的所述杂株核基因型计算的所述杂株 率，R2为去掉除所述细胞核杂株品种拥有的所述杂株核基因型的频率中最低的 80％和最高的10％的值后，剩余值的平均值的2倍；

$R3={\Sigma }_{i3=1}^{n2}R3i3-2\times R3\mathrm{ic},$其中，$R3i3=\frac{{\Sigma }_{j3=\mathrm{Int}(0.8\times t3)+1}^{t3-\mathrm{Int}(0.3\times t3)}R3i3j3}{t3-\mathrm{Int}(0.8\times t3)-\mathrm{Int}(0.1\times t3)},$n2为细胞质杂株品种的数目，R3i3为第i3个所述细胞质杂株品种的所述杂株率， R3ic为i3＝ic时R3i3的值，ic为当所述待测大豆品种为核质互作型不育系或保持 系时，对应的所述保持系或所述不育系的所述细胞质杂株品种，t3为第i3个所 述细胞质杂株品种的所有特异杂株质基因型的数目，i3j3为第i3个所述细胞质 杂株品种的所有所述特异杂株质基因型按频率由低到高排序后，第j3个所述特 异杂株质基因型，R3i3j3为第i3j3个所述特异杂株质基因型的频率，R3ic指混入 所述不育系中的所述保持系的杂株率或混入所述保持系中的所述不育系的杂株 率；R3为由杂株质基因型计算的所述细胞质杂株品种的所述杂株率的总和，所 述细胞质杂株品种的杂株率为去掉所述细胞质杂株品种中最低的80％和最高的 10％的所述特异杂株质基因型的频率后，剩余的所述特异杂株质基因型的频率的 平均值；

其中，t4为除所述细胞质杂株品种拥有的所 述杂株质基因型之外的且频率≥0.17％的所述杂株质基因型的数目，i4为除所述 细胞质杂株品种拥有的所述杂株质基因型之外的所有所述杂株质基因型按频率 由低到高排序后，第i4个所述杂株质基因型，R4i4为第i4个所述杂株质基因型 的频率；R4是利用除所述细胞质杂株品种拥有的所述杂株质基因型计算的所述 杂株率，R4为去掉除所述细胞质杂株品种拥有的所述杂株质基因型的频率中最 低的80％和最高的10％的值后，剩余值的平均值；

Int()为取整函数；

所述细胞核杂株品种是指仅利用核基因型计算获得的所述杂株品种，所述 细胞质杂株品种是指仅利用质基因型计算获得的所述杂株品种；所述特异杂株 核基因型是指仅为一个所述细胞核杂株品种所有的所述杂株核基因型；所述特 异杂株质基因型是指仅为一个所述细胞质杂株品种所有的所述杂株质基因型； 所述杂株核基因型是指所述杂株基因型为所述核基因型，所述核基因型指所述 基因型位于细胞核基因组上；所述杂株质基因型是指所述杂株基因型为所述质 基因型，所述质基因型是指所述基因型位于细胞质基因组上；基因型的频率是 指所述测序片段组中，代表所述基因型的测序片段数占所述基因型所在所述测 试区域的测序片段总数的比例。

进一步地，所述方法还包括采用以下方式判定所述待测大豆品种的一致性 和稳定性的结论正确的概率为：当所述待测大豆品种具有一致性和稳定性时， 结论正确的概率≥BINOM.DIST(M*SN,SN,R,TRUE)*BINOM.DIST(ΣSeN*M, ΣSeN,R,TRUE)；当所述待测大豆品种不具有所述一致性和稳定性时，结论正确 的概率≥BINOM.DIST((1-M)*SN,SN,(1-R),TRUE)*BINOM.DIST(ΣSeN*(1-M), ΣSeN,1-R,TRUE)；其中，ΣSeN为所有用于计算所述杂株率R的所述基因型的 频率所在所述测试区域的测序片段的总和，M为判断所述一致性和稳定性时所 选用的阈值，BINOM.DIST(M*SN,SN,R,TRUE)为所述待测大豆品种进行了SN 次抽样，实际抽得的所述杂株率R小于所述阈值M的概率，BINOM.DIST(ΣSeN *M,ΣSeN,R,TRUE)的意义为：对所述待测大豆品种进行了ΣSeN次抽样，实际抽 得的所述杂株率R小于阈值M的概率。

进一步地，当所述非通用测试区域不存在所述变异位点时，若判断所述待 测大豆品种具有特异性，结论正确的概率≥BINOMDIST((1-SD)*TRN,TRN,1-OD, TRUE)；若判断所述待测大豆品种不具有特异性，结论正确的概率≥BINOMDIS T(SD*TRN,TRN,OD,TRUE)，其中，TRN为检测成功的测试区域的数目，OD为 所述变异位点率，SD为判断特异性时所选用的阈值，BINOMDIST为excel 20 10中的函数，所述结论正确的概率表示为当判断所述待测大豆品种具有特异性 时，所述变异位点率大于SD的概率，当判断所述待测大豆品种不具有特异性时， 所述变异位点率小于SD的概率，所述检测成功的测试区域通过分析所述测序片 段组后获得。

具体地，获得所述杂株品种的方法包括：所述杂株品种为存在于所述数据 库中的品种，且所述杂株品种的潜在杂株基因型与所述杂株基因型间有相同基 因型的所述测试区域的数目占所述杂株品种具有所述潜在杂株基因型的所述测 试区域的总数的比例≥60％；所述杂株基因型指频率≥0.02％的所述潜在杂株基 因型；

所述潜在杂株基因型与所述待测大豆品种的所有基因型间的差异碱基的数 量≥2个或所述差异碱基中有非连续碱基的插入或缺失。

具体地，通过所述变异位点确定所述通用测试区域的方法为：

通过区分度计算区分度的值，其中，a为变异窗口区域中被检测 到的品种总数，bi为所述变异窗口区域中第i种基因型的品种数，且bi>1，k为 包含大于1个品种的基因型的数目，所述变异窗口区域为以每个单核苷酸变异 位点为中心，向所述单核苷酸变异位点的两侧各延伸测序列长度的1/2作为检测 的窗口；

所述通用测试区域为细胞质基因组上区分度大的区域或细胞核基因组上所 述区分度大且均匀分布的区域。

本发明实施例提供的技术方案带来的有益效果是：本发明实施例提供的方 法通过高通量测序和多位点扩增，实现了待测大豆品种的大样本抽样与种间个 体测试区域的大样本抽样，再利用定义杂株基因型、定义细胞质杂株品种和定 义杂株率计算公式等综合手段，成功地实现了准确、完整地判断待测大豆品种 的特异性、稳定性与一致性的目标，且测试速度更快，可在10天以内完成。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将对本发明实施方式 作进一步地详细描述。

实施例.测定大豆品种‘北93-406’的特异性、一致性与稳定性

本发明实施例提供的待测大豆品种为大豆品种“北93-406”，大豆品种“北 93-406”为纯系大豆且为公开使用品种，其与“北丰8号”杂交，通过系统选育 培育了“垦鉴豆26号”大豆品种。则测定该大豆品种的特异性、一致性与稳定 性的方法包括以下步骤。

一、获得不同大豆品种间的变异位点。

不同大豆品种的变异位点可以从已公布的文献资料中获取，但该方法所获 得的结果比较零星，在本实施例中，通过比较不同大豆品种的基因组序列获得 了大量的不同大豆品种间的变异位点。

进一步地，获得不同大豆品种的基因组序列的方法如下：

本实施例的不同大豆品种的基因组序列有两种来源，第一种为Lam等对31 个大豆品种的基因组的高通量测序序列，相关文献信息如下：Lam HM et al. Resequencing of 31wild and cultivated soybean genomes identifies patterns of  genetic diversity and selection.Nat Genet 2010,42:1053–1059。该31个大豆品种的 基因组序列公布于NCBI Short Read Archive(http://www.ncbi.nlm.nih.gov/sra)，接 收号为SRA020131；第二种为按Lam等的上述发表的文章中提供的方法对“北 丰8号”和“垦鉴豆26号”进行了高通量测序。本实施例共获得了33个大豆品种 的基因组的高通量测序序列。

进一步地，利用不同品种的基因组序列获得变异位点。

具体地，由于这33个大豆品种的测序深度都不高，仅能鉴定单核苷酸变异 (SNP)位点，若大豆品种的测序深度足够高，则能够鉴定出其它变异类型如重 复数变异，由于可信度低，不进行鉴定。利用Frederick Sanger比对软件(版本 号为0.4)将这33个大豆品种的基因组的高通量测序序列比对到“Williams_82” 大豆细胞核参考基因组(版本：Release v1.01，下载地址：http://genome.jgi-psf.org/) 和细胞质参考基因组上，该细胞质参考基因组包括线粒体参考基因组与叶绿体 参考基因组，其在NCBI(National Center for Biotechnology Information，美国国 立生物技术信息中心)上的接收号分别为JX463295.1和NC_007942.1。对比时， 插入片段长度设为500bp，其他参数设定为默认值。采用的Ssaha Pileup软件包 (版本号为0.5)鉴定每个大豆品种的SNP位点。该SNP位点定义为差异确定 的碱基对、单碱基的插入或单碱基的缺失。该差异确定的碱基对是指不包括差 异不确定的碱基对，差异不确定的碱基对是指某些简并碱基间的碱基对，如R 代表A或G，因此，A与R之间可能存在差异，也可能不存在差异，因此，A 与R间差异不明确，互不为SNP。因此，本发明实施例中的SNP位点为不包括 上述差异不确定的碱基对。按以上SNP位点的定义，本发明实施例在所有33 个大豆品种间共获得6350046个SNP位点，其中31937个SNP位点位于细胞质 基因组上，其余的SNP位点位于细胞核基因组上。后文提及的基因型即是指测 试区域内多个SNP位点的组合，核基因型指基因型位于细胞核基因组上，质基 因型是指基因型位于细胞质基因组上。例如，表1中第1个测试区域位于细胞 核基因组上，为核基因型，该测试区域共有3个SNP位点，该测试区域的基因 型即为这3个SNP位点的组合。

二、通过变异位点确定待测大豆品种的测试区，测试区域包括通用测试区 域，至少部分变异位点包含在通用测试区域内，其方法包括：

确定通用测试区域

通用测试区域为细胞质基因组上区分度大的区域或细胞核基因组上区分度 大且均匀分布的区域，其中，区分度其中，a为变异窗口区域中被 检测到的品种总数，bi为变异窗口区域中第i种基因型的品种数，且bi>1，k为 包含大于1个品种的基因型的数目，变异窗口区域为以每个单核苷酸变异位点 (SNP位点)为中心，向单核苷酸变异位点的两侧各延伸测序列长度的1/2作为 检测的窗口；测试区域为细胞质基因组上区分度大的区域或细胞核基因组上区 分度大且均匀分布的区域。区分度的计算原理如下：所有品种间的组合数为其中，同一基因型内的不同品种间的组合是不可区分的，其数目为那么， 不可被区分的品种组合的比例为可被区分的品种组合的比例即区分度 由此可见，区分度越大，越能将不同品种区分开，区分度大的变 异窗口区域对DUS测试更有效。若细胞核基因组上的变异窗口区域分布不均匀， 会导致某些区域相邻，从而连锁遗传，信息容易重叠，因此，细胞核基因组上 选择通用测试区域的综合原则是：区分度大且SNP位点均匀分布。细胞质基因 组无连锁遗传问题，所以，细胞质基因组上只需要选择区分度大的区域即可。

本发明实施例中采用Proton高通量测序仪进行高通量测序，其测序检测的 测试区域长度可达到200bp，为了获得最大信息量，本实施例中的最长测试区 域也为200bp。因此，本实施例提到的变异位点是指整个测试区域，其内部可 能包含多个SNP位点。

首先，以获得的每个SNP位点为中心，向左右各延伸99bp和100bp，构 成200bp的变异窗口。根据获得的6350046个SNP位点，可以获得6350046个 变异窗口，计算这些变异窗口区域的区分度例如，第1个变异窗 口区域中，共检测到了a＝29个品种，共有k＝2种基因型GTT、ACC，它们的品 种数分别为b1＝22个、b2＝5，因此，其含义是：通过第1 个变异窗口区域，可以将29个品种中的41％的品种组合区分开，另外59％的品 种组合无法区分开，需要更多的变异窗口才能区分开。按照同样的方法，计算 获得全部6350046个变异窗口的区分度并从中选取位于细胞核基因组中区分度 最大的8000个变异窗口和位于细胞质基因组中区分度最大的100个变异窗口。 逐个检查位于细胞核基因组的8000个变异窗口中，每个变异窗口与下一个变异 窗口间的距离，若距离超过300K(1K＝1000个碱基)，则放弃其中区分度较小 的变异窗口之后再检查，直至相邻查变异窗口的距离均大于300K为止。选择 300K的距离标准是因为大豆基因组大小约为975M(1M＝100万个碱基)，按最 终入选2000个位于细胞核基因组的通用测试区域计，平均的通用测试区域间距 离约为500K，但由于一些特异区域如着丝粒等很少有变异位点，因此，平均距 离应该小于500K。按以上方法，选出了4987个位于细胞核基因组的变异窗口， 它们与获得的位于细胞质基因组中区分度最大的100个变异窗口一起共5087个 变异窗口作为入选的通过测试区域。其中，选择区分度最大的200个变异窗口， 为经验值，该数量可以根据具体情况进行修改。

该测试区域还可以包括非通用测试区域，具体方法如下：

确定非通用测试区域：

非通用测试区域是指特殊品种需要检测的非通用位点。DUS测试需要检测 定点改造的特殊位点，定点改造是现代育种中常用的技术手段，如回交育种、 转基因育种等，定点改造品种也可以因其具有特异性而成为新品种。基于新品 种保护特异性的判定原则，非通用测试区域应不包括在通用测试区域内且为已 知控制质量性状的位点。本实施例中，由于待测大豆品种不是通过定点改造而 来的，无非通用位点需要检测，因此，无非通用测试区域。

三、制备扩增测试区域的引物，测试区域引物包括通用测试区域引物，具 体如下：

制备通用测试区域引物，该通用测试区域引物针对所有品种，具体地：

通用测试区域采用多重PCR技术进行检测，多重PCR技术是指在同一个 PCR反应中加入多个PCR引物，同时扩增基因组上的多个位点。该技术的关键 是设计并合成多重PCR引物，本实施例采用美国LifeTechnology公司提供的多 重PCR技术，其能够设置多至12000重PCR引物。

引物获取过程如下：登录LifeTechnology公司多重PCR引物在线设计网页 https://ampliseq.com/protected/help/pipelineDetails.action，按其要求提交相关信息 即可。在本实施例中，“Application type”选项选择“DNA Hotspot designs (single-pool)”。若选择multi-pool，则多重PCR将分多管进行，成本会有所增加， 而single-pool的引物只需要一次多重PCR即可，节省成本，缺点是某些通用测 试区域引物设计可能失败，但基因组上的备选的通用测试区域较多，因此，放 弃一些备选的通用测试区域并不影响结果。将待测大豆品种的细胞核参考基因 组和细胞质参考基因组融合为一个文件，并在“Select the genome you wish to use” 选项中选择“Custom”后，上传融合的文件作为设计多重PCR引物时的参考基 因组。DNA Type选项选择“Standard DNA”，在Add Hotspot选项中，添加需要 设计的通用测试区域内的SNP位点的位置信息即可，包括染色体信息、SNP的 起始位点和SNP的结束位点，其部分实例见表1。最后点击“Submit targets”按 钮提交并得到设计的多重PCR引物。本实施例中，从所有5087个通用测试区域 中，设计并成功验证了2488对多重PCR引物，用于扩增相应的2488个通用测 试区域。验证了多重PCR引物的方法为该公司按本发明提供的方法，提取同一 株大豆上的叶片基因组DNA，并利用设计的多重PCR引物对获得的基因组DNA 进行扩增、建库、高通量测序并分析测序片段组，去掉以下测试区域相应的引 物：该测试区域的测序片段数不足1000或存在杂株基因型，保留下来的引物即 为验证成功的多重PCR引物。由于基因组DNA来源于同一株大豆叶片，不可 能存在杂株品种，因此，杂株基因型是由测试区域的特殊结构造成的PCR或测 序偏好性错误，去掉这些测试区域避免了此类系统错误。验证成功的多重PCR 引物也由该公司混合好后以液体的形式提供给客户使用。上述成功设计了多重 PCR引物的2488个通用测试区域即为最终用于待测大豆品种检测的通用测试区 域，同时，构建的数据库中的每个品种也包含了上述2488个通用测试区域，其 中，47个通用测试区域位于细胞质基因组上，剩余的2441个通用测试区域位于 细胞核基因组上。

需要说明的是：通用测试区域的数目要求≥900个，理由如下：若低于900 个，存在误判的杂株品种的概率将超过1％，该阈值的推算方法见表2。由于可 能存在检测失败的测试区域，因此，测试区域数目一般≥1000个。

测试区域引物还可以包括非通用测试区域引物，该非通用测试区域引物针 对待测大豆品种，具体如下：

制备非通用测试区域引物：

非通用测试区域的引物包括第一引物和第二引物，第一引物包括第一正向 引物和第一反向引物，第二引物包括第二正向引物和第二反向引物，第一引物 和第二引物分别进行单独扩增得到两个非通用测试区域的扩增产物，将两个非 通用测试区域的扩增产物等量混合用于构建单独扩增的高通量测序文库。第一 正向引物的5’端连接有如序列表中SEQ ID NO:1所示的序列1，第一反向引物 中的5’端连接有如序列表中SEQ ID NO:2所示的序列2；第二正向引物的5’端 连接有如序列表中SEQ ID NO:2所示的序列2，第二反向引物的5’端连接有如 序列表中SEQ ID NO:1所示的序列1。

非通用测试区域引物的设计过程如下：第一步，按扩增长度不超过200bp 和包含非通用测试区域内所有SNP位点的要求，按普通PCR引物设计方法，设 计扩增非通用测试区域的PCR的正向引物和反向引物；第二步，将设计好的正 向引物与反向引物的5’端分别连接序列表中SEQ ID NO:1和序列表中SEQ ID  NO:2，分别获得第一引物的正向引物和第一引物的反向引物；第三步，将设计 好的正向引物与反向引物的5’端分别连接序列表中SEQ ID NO:2和序列表中 SEQ ID NO:1，分别获得第二引物的正向引物和第二引物的反向引物。序列表中 SEQ ID NO:1和序列表中SEQ ID NO:2为高通量测序所用的接头序列，从而使 用PCR产物带有高通量测序的接头序列，可以直接与扩增的通用测序区域的产 物混合后建立测序文库后一同测序，而不必经过片段化、连接接头等繁琐的建 库步骤，提高了工作效率并降低了成本。做成两对仅接头不一样的引物是为了 同时从非通用测试区域的两端测序。

本实施例中的待测大豆品种由于没有非通用测试区域，因此，无需非通用 测试区域引物。

四、构建包含不同大豆品种在所有测试区域的基因型的数据库的方法如下：

本实例获得了2488个通用测试区域引物和0个非通用测试区域引物，它们 对应的扩增区域即为待测大豆品种的测试区域。构建包含33个品种的2488测 试区域的基因型及其SNP的位置信息的数据库，部分结果见表1。

表1为数据库品种基因型及其位置、待测大豆品种基因型、杂株基因型及 其频率的部分实例

表1中，‘-’代表该SNP位点的位置在参考基因组上缺失；除ATGC外，其 它字母代表简并碱基。若基因型全由简并碱基N组成，称相应测试区域基因型 与SNP数据缺失，缺失的基因型或SNP与任何基因型或SNP比较时，均作无 差异处理。可按本发明提供的检测待测大豆品种基因型的方法检测数据库品种 并补全缺失的基因型。

由于篇幅限制，本实施例没有完整列出全部数据库内容，只列出了其中5 个品种的10个测试区域的信息。同样基于篇幅限制，本实施例中还有部分地方 也仅列出部分相关实例，其余未列出的数据可根据本实施例的方法补全。

五、确定待测大豆品种的抽样量SN后，随机抽样混合并提取混合样本的 DNA，方法如下：

计算待测大豆品种抽样量

抽样量SN应满足如下条件：BINOM.INV(SN,M,0.95)/SN≤1.15*M，其中， M为判断一致性和稳定性时所选用的阈值，BINOM.INV为excel 2010中的函数， 其使用方法与excel 2010中的定义相同，其含义是使得累积二项分布的函数值 大于或等于临界值的最小整数。抽样量SN满足的条件含义是：即使杂株率只超 出一致性和稳定性时的判断阈值M的15％，该抽样量在95％的概率保证下，可 正确判断待测大豆品种的稳定性与一致性。M值是根据作物种类、标记类型、 具体要求等条件人为确定的。在农业部新品种保护办公室发布中《植物新品种 特异性、一致性和稳定性测试指南-大豆》中规定：对于常规种，一致性判定时， 采用0.5％的群体标准和至少95％的接受概率。因此，本实施例中，选用中间 值0.5％作为M值。逐步加大SN值后，计算上述公式发现，当SN≥25218时， BINOM.INV(SN,0.5％,0.95)/SN≤1.15*0.5％成立。因此，本实施例中的待测大豆 品种抽样量应≥25218。

随机抽样混合并提取混合样本的DNA

在本实施例中，选取了30000粒种子发芽，随机选取26000个大小大致相 等的芽混合后置于研钵中，向研钵中加入液氮后充分研磨成粉。采用北京天根 生化科技有限公司生产的货号为DP305的植物基因组DNA提取试剂盒提取并 获得待测大豆品种混合样本的DNA，DNA提取方法按该试剂盒的操作手册进 行。利用美国Invitrigen公司生产的dsDNA HS Assay Kit(货号为Q32852) 及其说明书对获得的DNA进行定量，将定量后的待测大豆品种DNA稀释为 10.00ng/μl。

六、利用引物对混合样本的DNA进行扩增，得到测试区域的扩增产物，扩 增产物作为高通量测序文库的方法如下：

高通量测序文库包括：通用测试区域的高通量测序文库和非通用测试区域 的高通量测序文库，在本实施例中，分别构建通用测试区域和非通用测试区域 的高通量测序文库，将二者混合，得到所有测试区域的高通量测序文库。

构建通用测试区域的高通量测序文库的方法如下：

利用文库构建试剂盒2.0(由美国LifeTechnology公司生产，货号为 4475345)多重PCR扩增通用测试区域后，利用扩增产物构建高通量测序文库。 该试剂盒包括以下试剂：5×Ion AmpliSeq^TM HiFi Mix、FuPa试剂、转换试剂、 测序接头溶液和DNA连接酶。文库构建的方法按该试剂盒的操作手册《Ion  AmpliSeq^TM Library Preparation》(出版号：MAN0006735，版本：A.0)进行。 通过多重PCR扩增2488个通用测试区域，多重PCR的扩增体系如下：5×Ion  AmpliSeq^TMHiFi Mix 4μl、制备的通用测试区域引物混合液4μl、待测大豆品种 的DNA 10ng和无酶水11μl。多重PCR的扩增程序如下：99℃，2分钟；(99℃， 15秒；60℃，4分钟)×25个循环；10℃保温。利用FuPa试剂消化掉多重PCR 扩增产物中多余的引物后，再进行磷酸化，具体方法为：向多重PCR的扩增产 物中加入2μL FuPa试剂，混匀后，在PCR仪上按如下程序反应：50℃，10分 钟；55℃，10分钟；60℃，10分钟；10℃保存，得到混合物a，混合物a为含 有经过磷酸化的扩增产物溶液。将磷酸化的扩增产物连接上测序接头，具体方 法为：向混合物a中加入转换试剂4μL、测序接头溶液2μL和DNA连接酶2 μL，混匀后，在PCR仪上按如下程序反应：22℃，30分钟；72℃，10分钟； 10℃保存，得到混合液b。利用标准的乙醇沉淀方法纯化混合液b后溶解于10μL 无酶水中。利用美国Invitrigen公司生产的dsDNA HS Assay Kit(货号为 Q32852)并按照其说明书进行测定，并获得混合液b的质量浓度后，将纯化后 混合液b稀释至15ng/ml，得到浓度约100pM的通用测试区域的高通量测序文 库。

构建非通用测试区域的高通量测序文库的方法如下：

以待测大豆品种的DNA为模板，利用上述制备的非通用测试区域的第一引 物和第二引物分别进行单独PCR扩增，等量混合扩增产物后得到非通用测试区 域的高通量测序文库。具体操作按《Ion Amplicon Library Preparation(Fusion  Method)》(出版号：4468326)进行，大致过程如下：将第一引物的正向引物和 反向引物用水溶解为10μM的浓度后，等体积混合，获得第一引物溶液。配制 如下PCR反应体系：第一引物溶液1μL、30ng待测大豆品种DNA和PCR高 保真混合物(美国invirtrigen公司生产，货号为12532016)45μL，混匀后，在 PCR仪上按如下程序反应：94℃，3分钟；(94℃，30秒；58℃，30秒；68℃， 1分钟)×40个循环；4℃保温。PCR扩增产物按标准的乙醇沉淀的方法纯化后 溶解于10μL水中，利用DNA 1000试剂盒(货号为5067-1504)在美国安捷伦 公司生产的生物分析仪(型号为2100)上，按该试剂盒说明书测定并获得扩增 产物的摩尔浓度后，稀释为200pM，即为第一引物的扩增产物。采用相同的方 法，获得浓度为200pM的第二引物的扩增产物。将第一引物的扩增产物与第二 引物的扩增产物等体积混合，获得浓度为100pM的非通用测试区域高通量测序 文库。本实施例中，由于无非通用测试区域，因此，也无需构建非通用测试区 域的高通量测序文库。

获得所有测试区域的高通量测序文库

按通用测试区域的数目与非通用测试区域的数目的比例混合等摩尔浓度的 通用测试区域的高通量测序文库和非通用测试区域的高通量测序文库，得到的 混合物即为所有测试区域的高通量测序文库。在本实施例中，因无非通用测试 区域的高通量测序文库，因此，构建的高通量测序文库即为浓度为100pM的通 用测试区域的高通量测序文库。

七、对高通量测序文库进行高通量测序，得到测序片段组，方法如下：

确定高通量测序深度的原理：高通量测序的深度满足如下条件：BINOM.DI ST(10,10,BINOM.DIST(8,20,BINOM.DIST(0,CF,0.1％,TRUE),TRUE),FALSE)≥99. 9％，1-BINOM.DIST(10000,10000,1-BINOM.DIST(8,20,1-BINOM.DIST(99.99％* CF,CF,99.9989％,TRUE),TRUE),FALSE)≤0.1％且BINOM.DIST(10*(1-M)*CF,10* CF,1-110％*M,TRUE)≥95.0％，其中，CF为高通量测序的深度，也即平均每个测 试区域被覆盖的倍数，M为判断一致性和稳定性时所选用的阈值，BINOM.DIS T为excel 2010中的函数，其使用方法与excel 2010中的定义相同，其返回的 是二项式分布的概率。该三个函数的意义为：在杂株率低至0.1％、杂株品种多 达10个且杂株品种与待测大豆品种间平均仅20个差异位点的条件下，由高通 量测序深度决定的检出全部杂株品种的概率≥99.9％；在数据库品种多至10000 个且杂株品种与待测大豆品种间平均仅20个差异位点的条件下，由高通量测序 深度决定的存在误判杂株品种的概率≤0.1％；在杂株品种多达10个且真实杂株 率仅超过判断特异性时所选用的阈值的10％时，由高通量测序深度决定的对稳 定性与一致性的判定结论正确的概率≥95.0％。以上条件十分严格，因此，真实 效果优于上述阈值。以上概率的推算方法见表2。

表2为本实施例相关概率的计算方法

表2为Excel 2010数据表，其函数、单元格等均与Excel 2010的定义相同。 其中，“判断一致性和稳定性时所选用的阈值(M)”为单元格B31，其它单元格 编号以B31为参照，按Excel 2010的规则定义，例如“杂株率(R)”所在的单元 格在B2的基础上增加了4行1列，因此编号为C35，其它单元格编号规则与此 相同。

本实施例高通量测序深度的确定方法为：将M＝0.5％代入上述三个公式后， 逐步加大测序深度CF至6001时，可使上述三个方程成立，因此，本实施例测 序深度确定为≥6001倍。

利用高通量测序文库进行高通量测序

利用获得的所有测试区域的高通量测序文库和试剂盒Ion PI Template OT2 200Kit v2(美国invirtrigen公司生产，货号为4485146)进行测序前的ePCR (Emulsion PCR，乳化聚合酶链反应)扩增，操作方法按该试剂盒的操作手册 进行。利用ePCR产物和试剂盒Ion PI Sequencing 200Kit v2(美国invirtrigen公 司生产，货号为4485149)在Proton二代高通量测序仪上进行高通量测序，操作 方法按该试剂盒的操作手册进行。在本实施例中，高通量测序通量设置为平均 覆盖测试区域10000倍。

对高通量测序结果进行预处理

首先判断高通量测序的数据质量是否≥Q20，若<Q20(此情况极少)，则按 上述方法重新进行高通量测序，直到质量要求达到Q20标准，Q20标准满足了 表2中“测序错误为特定碱基的概率”≤0.33％的要求。将达到质量要求的高通量 测序片段比对到所有2488个测试区域，去掉比对不成功和基因型检测不全的测 序片段后，剩余的所有测序片段称为测序片段组。基因型检测不完全的测序片 段是指没能将表1中“SNP在参考基因组上的位置”所示的该测序片段所在的测 序区域中的所有SNP位点检测到，基因型检测不全的原因是测序片段过短，比 对不成功的原因是测序片段多为非特异扩增产物。

八、分析测序片段组，获得待测大豆品种基因型和杂株基因型，方法如下：

将测序片段组比对到所有测试区域，并统计每个测试区域中的测序片段数， 去掉测序片段数≤1000条的测试区域，剩余的测试区域为检测成功的测试区域。 在本实施例中，共获得2406个检测成功的测试区域。比对到测试区域的片段称 为该测试区域的测序片段，从测序片段中提取表1中“SNP在参考基因组上的位 置”所示的位置的碱基组合称为该测序片段的基因型。基因型的频率是指测序片 段组中，代表该基因型的测序片段数占该基因型所在测试区域的测序片段总数 的比例。频率最大的基因型称为待测大豆品种基因型。杂株基因型指频率≥0.02％ 的潜在杂株基因型，其中，潜在杂株基因型与待测大豆品种的所有基因型间的 差异碱基的数量≥2个或差异碱基中有非连续碱基的插入或缺失。杂株基因型定 义的原理为：高通量测序中，插入或缺失错误极为罕见，而因测序错误造成2 个固定差异碱基的概率低至(1％/3)2＝0.0011％，且要求杂株基因型频率≥0.02％， 在这些条件限制下，即使是30000的测序深度，因测序错误产生某种杂株基因 型的概率仅为0.0001％(计算方法见表2)。0.02％的频率满足目前最严格DUS 测试标准，即从1万粒种子中检出的低至2粒杂种子。若差异碱基数量＝1个， 则全部测试区域都会产生错误的杂株基因型(计算方法见表2)，若差异碱基数 量≥3个时，杂株基因型数量急剧减少，难以准确计算杂株率R，因此，差异碱 基数量≥2个的阈值是最优的。

例如，在测序片段组中，第1个测序区域的测序片段总数为9987条，有TCA、 TCG、TCC、TCT……共25种基因型，代表这些基因型的测序片段数分别9612 条、218条、1条、2条……，这些基因型的频率为9612/9987＝96.25％、 218/9987＝2.18％、1/9987＝0.01％、2/9987＝0.02％……。按待测大豆品种基因型和 杂株基因型的定义，TCA应该是待测大豆品种在第1个测试区域的待测大豆品 种基因型，而TCG的频率超过0.02％，但与待测大豆品种基因型TCA比较有1 个<2个碱基的差异，因此TCG不为杂株基因型，其它基因型为测序错误产生的 基因型。杂株核基因型是指杂株基因型为核基因型，杂株质基因型是指杂株基 因型为质基因型。按此定义，第一个测试区域的也无杂株核基因型。按相同的 方法，判断并获得全部2406个检测成功的测试区域的待测大豆品种基因型、杂 株基因型及其频率，并判断获得的杂株基因型是杂株核基因型还是杂株质基因 型。结果表明：共获得291个杂株基因型，其中，286个为杂株核基因型，5个 为杂株质基因型。

下面简单介绍一下本实施例中的标准样品检测方法，从待测大豆品种中取1 粒种子，播种并长成幼苗后，利用幼苗的叶片按与待测大豆品种相同的方法提 取基因组DNA，该DNA称为待测大豆品种的标准样品。与待测大豆品种同时 且按相同方法平行构建标准样品的高通量测序文库并高通量测序。其中，频率 最大的基因型称为标准样品基因型，标准样品杂株基因型的频率≥0.02％且与标 准样品基因型间的差异碱基的数量≥2个或差异碱基中有非连续碱基的插入或缺 失。按与待测大豆品种相同的方法，获得每个检测成功的测试区域中的标准样 品基因型与标准样品杂株基因型。若标准样品基因型和待测大豆品种基因型相 同的测试区域占标准样品与待测大豆品种均检测成功的测试区域的比例超过 90％，则标准样品正确，否则，重新从待测大豆品种中取1粒种子，重复以上过 程，直至获得正确的标准样品。将正确的标准样品的杂株基因型与待测大豆品 种对应测试区域的杂株基因型比较，获得相同的杂株基因型，去掉待测大豆品 种中相同的杂株基因型，正确的待测大豆品种杂株基因型被保留下来并用于后 续分析。以上措施去掉了因系统选择性错误造成的杂株基因型，系统选择性错 误主要是基因序列的特殊结构所造成的PCR选择性错误扩增。需要说明的是： 当数据库品种多，能广泛代表不同品种基因型时，可以要求杂株基因型与数据 库品种的某个基因型相同，同样可以起到与标准样品相同的功能，在此情况下， 可以不检测标准样品，达到减轻工作量的目的。本实施例中结果为：从获得的 291个杂株基因型中，共去掉了2个杂株基因型，其中2个为杂株核基因型，0 个为杂株质基因型，所保留下来的289个杂株基因型用于后续分析，部分结果 见表1。

九、将待测大豆品种基因型与数据库中的不同品种的基因型比较，获得近 似品种、变异位点和变异位点率，方法如下：

若在测试区域中，待测大豆品种与数据库品种的基因型均无缺失，称该测 试区域为待测大豆品种与该数据库品种的共有测试区域。在共有测试区域中， 若待测大豆品种与数据库品种的基因型不完全相同，则称该不完全相同的基因 型所在的测试区域为待测大豆品种与该数据库品种的差异位点，相应的基因型 互为差异基因型，差异位点率＝差异位点的数目/共有测试区域的数目。从数据库 中获得差异位率最小的品种称为待测大豆品种的近似品种，相应的差异位点称 为变异位点，变异位点率＝变异位点的数目/共有测试区域的数目。

在本实施例中，待测大豆品种与数据库的第1个品种“北丰8号”的共有测试 区域数为2335个。在第1个共有测试区域中，待测大豆品种与“北丰8号”基因 型分别为GTT和GTT，二者相同，因此，第1个共有测试区域不为待测大豆品 种与“北丰8号”的差异位点，GTT也不为待测大豆品种与“北丰8号”的差异基 因型。按相同的方法，将所有共有测试区域中，待测大豆品种与“北丰8号”基因 型比较，发现共有190个差异位点，差异位点率＝190/2335＝8.14％。按相同的方 法，获得待测大豆品种与数据库中所有33个品种间的差异位点率，并获得差异 位点率最小的品种为“垦鉴豆26号”，差异位点率为3.24％。因此，“垦鉴豆26 号”为待测大豆品种的近似品种，待测大豆品种的变异位点率为3.24％。

十、将杂株基因型与数据库中的不同品种的基因型比较，获得杂株品种后， 计算杂株率，方法如下：

获得杂株品种：杂株品种存在于数据库中的品种，且杂株品种的潜在杂株 基因型与杂株基因型间有相同基因型的测试区域的数目占杂株品种具有潜在杂 株基因型的测试区域的总数的比例≥60％，其中，潜在杂株基因型与待测大豆品 种的所有基因型间的差异碱基的数量≥2个或差异碱基中有非连续碱基的插入 或缺失。杂株品种分为细胞核杂株品种和细胞质杂株品种，其中，细胞核杂株 品种是指仅利用核基因型计算获得的杂株品种，细胞质杂株品种是指仅利用质 基因型计算获得的杂株品种。例如，假设数据库中的品种的基因型分别为AA、 AA、AA/TT(“/”表示该测试区域为杂合基因型，存在“/”前后两种不同的基因型)、 AA/TT、AA/TT、AA/TT和AA时，待测大豆品种相应的基因型分别为AA、 AA/TT、TT、AA、TT/CC、GG/CC和-A时，相应的潜在杂株基因型为：无、 无、AA、TT、AA、AA/TT和AA。一般纯系品种中不存在杂合基因型，但极 少数位点可能存在，另外，杂株多为杂交种，杂合位点较常见，因此列举了各 种可能情况。参数60％可以保证全部杂株品种检出概率为100％且存在误判的杂 株品种的概率为0％，该参数值的确定方法见表2。

在本实施例中，第1个测试区域中，数据库中第一个品种“北丰8号”与待测 大豆品种的基因型分别为GTT和GTT，二者间无碱基的差异，因此，GTT不为 潜在杂株基因型，因而第1个测试区域中，“北丰8号”也不存在与杂株基因型相 同的潜在杂株基因型，按相同的方法，逐个判断所有核基因型的测试区域中， 数据库中第一个品种“北丰8号”的基因型是否为潜在杂株基因型，若为潜在杂株 基因型，再判断潜在杂株基因型与杂株基因型间是否有相同基因型，结果表明， “北丰8号”共有134个具有潜在杂株基因型的测试区域，它们与相同测试区域的 杂株基因型间有相同基因型的测试区域数为133个，其比例为133/134＝99.25％> 60％，因此，判断“北丰8号”为细胞核杂株品种。按类似的方法，利用所有质基 因型的测试区域，判断“北丰8号”不为细胞质杂株品种。按相同的方法，判断数 据库中所有其它品种是否为细胞核杂株品种或细胞质杂株品种，结果表明：仅 “北丰8号”为细胞核杂株品种，没有发现细胞质杂株品种。以上结果说明：“北 丰8号”可能是通过飞花传粉而不是机械混杂，将基因型混入了待测大豆品种。

获得特异杂株基因型：特异杂株基因型是指仅为一个杂株品种所有的杂株 基因型，其包括特异杂株核基因型和特异杂株质基因型；特异杂株核基因型是 指仅为一个细胞核杂株品种所有的杂株核基因型，特异杂株质基因型是指仅为 一个细胞质杂株品种所有的杂株质基因型。本实施例中，共获得291个杂株基 因型，其中，286个为杂株核基因型，5个为杂株质基因型。第一个杂株核基因 型ATGA仅为细胞核杂株品种“北丰8号”所有，所以，ATGA为“北丰8号”的 特异杂株核基因型。按相同的方法，逐一判断所有获得的286个杂株基因型中， 133个为“北丰8号”所拥有的特异杂株核基因型。按类似的方法，判断5个杂株 质基因型均不为特异杂株质基因型。

计算杂株率R原理，具体如下：

杂株率R＝R1+R2-R3-R4，其中：其中， n1为细胞核杂株品种的数目，t1为第i1个细胞核杂株品种的所有特异杂株核基 因型的数目，i1j1为第i1个细胞核杂株品种的所有特异杂株核基因型按其频率 由低到高排序后，第j1个特异杂株核基因型，R1i1j1为第i1j1个特异杂株核基 因型的频率；R1为由杂株核基因型计算的细胞核杂株品种的杂株率的总和，细 胞核杂株品种的杂株率为去掉细胞核杂株品种中最低的80％和最高的10％的特 异杂株核基因型的频率后，剩余的特异杂株核基因型的频率的平均值的2倍； 其中，t2为除细胞核杂株品种拥有的杂株核基因 型之外的且频率≥0.17％的杂株核基因型的数目，i2为除细胞核杂株品种拥有的 杂株核基因型之外的所有杂株核基因型按其频率由低到高排序后，第i2个杂株核 基因型，R2i2为第i2个杂株核基因型的频率；R2是利用除细胞核杂株品种拥有 的杂株核基因型计算的杂株率，其为去掉除细胞核杂株品种拥有的杂株核基因 型的频率中最低的80％和最高的10％的值后，剩余值的平均值的2倍； $R3={\Sigma }_{i3=1}^{n2}R3i3-2\times R3\mathrm{ic},$其中，$R3i3=\frac{{\Sigma }_{j3=\mathrm{Int}(0.8\times t3)+1}^{t3-\mathrm{Int}(0.1\times t3)}R3i3j3}{t3-\mathrm{Int}(0.8\times t3)-\mathrm{Int}(0.1\times t3)},$n2为 细胞质杂株品种的数目，R3i3为第i3个细胞质杂株品种的杂株率，R3ic为i3＝ic 时R3i3的值，ic为当待测大豆品种为核质互作型不育系或保持系时，对应的保 持系或不育系的细胞质杂株品种，t3为第i3个细胞质杂株品种的所有特异杂株 质基因型的数目，i3j3为第i3个细胞质杂株品种的所有特异杂株质基因型按其 频率由低到高排序后，第j3个特异杂株质基因型，R3i3j3为第i3j3个特异杂株 质基因型的频率，R3ic指混入不育系中的保持系的杂株率或混入保持系中的不育 系的杂株率；R3为由杂株质基因型计算的细胞质杂株品种的杂株率的总和，细 胞质杂株品种的杂株率为去掉细胞质杂株品种中最低的80％和最高的10％的特 异杂株质基因型的频率后，剩余的特异杂株质基因型的频率的平均值； 其中，t4为除细胞质杂株品种拥有的杂株质基因 型之外的且频率≥0.17％的杂株质基因型的数目，i4为除细胞质杂株品种拥有的 杂株质基因型之外的所有杂株质基因型按其频率由低到高排序后，第i4个杂株质 基因型，R4i4为第i4个杂株质基因型的频率；R4是利用除细胞质杂株品种拥有 的杂株质基因型计算的杂株率，其为去掉除细胞质杂株品种拥有的杂株质基因 型的频率中最低的80％和最高的10％的值后，剩余值的平均值；Int()为取整函 数，返回括号中的数的整数部分。

待测大豆品种中的杂株来自于繁殖过程中的飞花传粉混杂和机械混杂，其 中，飞花传粉混杂是杂株品种混杂的主要来源。飞花传粉混杂是指杂株品种的 花粉通过风力等传到待测大豆品种并授粉形成的杂交种子，飞花传粉不可能引 入细胞质，因此只会造成杂株核基因型，其杂株率为杂株核基因型频率的2倍。 机械混杂是指杂株品种种子直接混入待测大豆品种中，同时引入了细胞核与细 胞质，同时形成杂株核基因型和杂株质基因型，其杂株率应该为杂株质基因型 的频率。在杂株率R的计算公式中，R1+R2将机械混杂的杂株率高估了1倍， 需校正，校正后的R＝R1+R2-R3-R4。区分机械混杂与飞花传粉混杂是一个技术 难题，本发明解决了这一难题。

在杂株率R的计算公式中，细胞核杂株品种的杂株率都是2×杂株核基因型 频率，其理由如下：二倍体或异源多倍体大豆在细胞核基因组的测试区域是2 个拷贝，因此，杂株率是相应杂株核基因型频率的2倍。若一定要选择有N份 拷贝的细胞核基因组的测试区域，则系数应调整为N，若拷贝数不明确，作N＝2 处理，若有误，将会在计算R时，通过去掉80％的低极端值的方式将它们排除。

在杂株率R的计算公式中，仅利用了杂株基因型频率值处于中间的10％的 进行计算，其原理为：同一杂株品种的不同杂株基因型由该杂株品种的杂株率 决定，所以频率的期望值是相等的，为频率间的差异由PCR扩增、高通量测序 过程中的误差引起。通过杂株基因型的定义与待测大豆品种标准样品，已基本 将这些错误值去掉了，去掉10％的极端值足以去掉极少量偏离真实杂株率的测 试区域。之所以去掉最小的80％，而最大的则只去除10％，原理如下：(1)最 大误差来源为测序错误，而测序错误产生的杂株基因型频率很低；(2)在除杂 株品种之外的杂株基因型的频率中，高值更可能为不同杂株的共同杂株基因型， 代表着真实的杂株率。

当待测大豆品种为核质互作型不育系时，若其中混有该不育系对应的保持 系杂株品种，那么，由于该保持系杂株品种的细胞质与待测大豆品种不一样， 将被检测为细胞质杂株品种，但由于不育系与保持系的细胞核完全一样，将不 会被检测为细胞核杂株品种，因此，R3ic的值在R1+R2中没有被计算到，但在 R3i3中被计算到了，因此，需要在R3中减去2×R3ic进行效正。同样的道理， 当待测大豆品种为核质互作型保持系时，也需要在R3中减去对应的不育系杂株 品种的2×R3ic进行效正。显然，当待测大豆品种既不为核质互作型不育系也不 为核质互作型保持系时，R3ic＝0。

在R2与R4的计算公式中，要求杂株基因型的频率≥0.17％，其原理如下： 当数据库中的品种数与检测位点均达到10000个时，平均将产生149次杂株基 因型误判，当设置杂株基因型频率≥0.17％时，无误判的杂株基因型的概率 ≥99.98％(推算方法见表2)，才可准确计算到R2与R4的值。数据库中的品种 数与检测位点均达到10000个已是现实中的极限，因此，杂株基因型的频率 ≥0.17％的阈值可以适用于各种情况。R2与R4的引入，使得本发明能够在数据 库品种为0，即没有数据库支持的情况下，计算杂株率R。

特别地，若杂株品种A的所有杂株基因型为杂株品种B和其它杂株品种所 拥有，因而，杂株品种A无特异杂株基因型。此时，在计算杂株率R时，不计 算杂株品种A与杂株品种B的杂株率，而计算杂株品种AB的杂株率。杂株品 种AB的杂株基因型定义为：杂株品种A与杂株品种B所共有杂株基因型。

杂株率R的计算公式为通用公式，在现实中待测大豆品种一般只混杂1种 杂株品种。

计算杂株率R的假定实例

表3假定了一个杂株率计算实例，以便更清楚说明杂株率R的计算过程。

表3为计算杂株率R的一个假定实例

表3中，细胞核杂株品种共A和B两个，所以n1＝2，细胞质杂株品种数仅 C一个，所以n2＝1。按特异杂株核基因型的定义，获得杂株品种A的特异杂株 核基因型为编号为1-10号的杂株核基因型AA、TT、TCC、GG、AC、TTC、 TCCC、GGC、ACC和AG，所以，t1＝10，它们的频率分别为0.10％、1.20％、 0.10％、0.10％、0.02％、0.10％、0.10％、0.10％、0.10％和0.10％，对这10个特 异杂株核基因型频率由低到高排序后，为R11111＝0.02％、R11121＝0.02％、 R11131＝0.10％、R11141＝0.10％、R11151＝0.10％、R11161＝0.10％、R11171＝0.10％、 R11181＝0.10％、R11191＝0.10％和R111101＝1.20％。从j 1＝Int(0.8×t1)+1＝Int (0.8×10)+1＝9到j 1＝t1-Int(0.1×t1)＝10-Int(0.1×10)+1＝9的R111j1的 值为R11191＝0.10％，所以细胞核杂株品种A的杂株率为按同样的方法，获得细胞核杂株品种B的杂株率为 $R121=\frac{2\times 0.20\%+2\times 0.20\%}{2-0-0}=0.40\%.$由此，获得细胞核杂株品种$R1={\Sigma }_{i1=1}^{2}R1i1=$$R111+R121=0.60\%.$按类似的方法，获得R2＝0.02％，细胞质杂株品种的杂株率 R4＝0.04％。因此，该假定实例中杂株率 R＝R1+R2-R3-R4＝0.60％+0.02％-0.10％-0.04％＝0.48％。

参照上述假定实例，计算本实施例中的杂株率R：在本实施例中，杂株品种 仅为“北丰8号”且为细胞核杂株品种，R2、R3与R4均为0，因而，R＝R1＝R111。 “北丰8号”共有133个特异杂株核基因型，频率为：1.02％、1.03％......(部分实 例见表1)，按R的计算规则，去除最小的80％(106个)和最小的10％(13个) 的频率值后，剩余的14个频率的平均值即为杂株率R＝1.03％。

十一、利用变异位点、变异位点率和杂株率，判断待测大豆品种的特异性、 一致性和稳定性，方法如下：

其中，SD为判断特异性时所选用的阈值，M为判断一致性和稳定性时所选 用的阈值。判断待测大豆品种特异性、一致性和稳定性的方法为：当变异位点 率≥SD或非通用测试区域存在变异位点时，待测大豆品种具有特异性，当变异 位点率＜SD且变异位点不存在于非通用测试区域中时，待测大豆品种不具有特 异性；当待测大豆品种的杂株率≤M时，待测大豆品种具有一致性和稳定性，当 待测大豆品种的杂株率大于＞M时，待测大豆品种不具有一致性和稳定性。与 M值一样，SD值是根据育种水平、要求的严格程度，标记特性等诸多因素，人 为确定的。在本实施例中，SD选用1％的标准。

在本实施例中，变异位点率为3.24％>SD＝1％，因此，判断待测大豆品种具 有特异性；待测大豆品种的杂株率1.03％≤M＝0.5％，因此，判断待测大豆品种不 具有一致性和稳定性。

进一步地，在判断待测大豆品种特异性、一致性与稳定性后，对判断的正 确率进行估计，方法如下：

本发明中的纯系大豆新品种是指以纯系基因型为目标而选育的常规种、自 交系、恢复系、保持系、不育系等类型。

特异性正确率计算：当非通用测试区域不存在变异位点时，若判断待测大 豆品种具有特异性，结论正确的概率≥BINOM.DIST((1-SD)*TRN,TRN,1-OD,TR UE)；若判断待测大豆品种不具有特异性，结论正确的概率≥BINOM.DIST(SD* TRN,TRN,OD,TRUE)，其中，TRN为成功检测的测试区域的数目，OD为变异 位点率,BINOM.DIST为excel 2010中的函数，其使用方法与excel 2010中的 定义相同，其返回的是二项式分布的概率。上述概率实际上计算的是：当判断 具有特异性时，变异位点率大于SD的概率；当判断不具有特异性时，变异位点 率小于SD的概率，检测成功的测试区域通过分析测序片段组后获得。

本实施例中，采用变异位点率来判断待测大豆品种的具有特异性，因此， 特异性结论正确的概率≥BINOM.DIST((1-1％)*2406,2406,1-3.24％,TRUE)＝100.0 0％，可见本实施例对特异性判定结论的正确率是很高的。

一致性与稳定性正确率计算

判定待测大豆品种的一致性和稳定性的结论正确的概率为：当待测大豆品 种具有一致性和稳定性时，结论正确的概率≥BINOM.DIST(M*SN,SN,R,TRUE)* BINOM.DIST(ΣSeN*M,ΣSeN,R,TRUE)；当待测大豆品种不具有一致性和稳定 性时，结论正确的概率≥BINOM.DIST((1-M)*SN,SN,(1-R),TRUE)*BINOM.DIST (ΣSeN*(1-M),ΣSeN,1-R,TRUE)，其中，ΣSeN为所有用于计算杂株率R的基因 型频率所在测试区域的测序片段的总和，也即去掉80％的最小值和10％的最大 值后，保留下来用于计算杂株率的测试区域的测试片段的总和，M为判断一致 性和稳定性时所选用的阈值。判断一致性与稳定性的正确率完全取决于杂株率 的正确率，而杂株率的正率确取决于以下三个步骤的正确率：第一，待测大豆 品种抽样正确率，第二，从抽出样本中检出杂株品种的正确率，第三，利用检 出的杂株品种计算杂株率的正确率。因此，判断待测大豆品种一致性与稳定性 的正确率是以上三步正确率的积。由于本发明即使是在最严格的条件下，检出 杂株品种的正确率也控制99.9％以上，实际上绝大部分是接近100％的。例如， 在本实施例中，全部杂株品种检出概率在100.0000％以上，存在误判的杂株品种 的概率在0.0000％以下(具体计算方法见表2)。因此，判断待测大豆品种一致 性与稳定性的正确率可估算为第一步和第三步的正确率的积，其分别为上述公 式中前后两个函数所计算的值。例如，BINOM.DIST(M*SN,SN,R,TRUE)的意义 为：待测大豆品种进行了SN次抽样，实际抽得的杂株率R小于阈值M的概率； 用于计算待测大豆品种杂株率的每一个测序片段，实质上也相当对待测大豆品 种进行了一次抽样，因此，BINOM.DIST(ΣSeN*M,ΣSeN,R,TRUE)的意义为： 对待测大豆品种进行了ΣSeN次抽样，实际抽得的杂株率R小于阈值M的概率。

本实施例中，去掉最小的80％和最大10％的杂株基因型频率后，共有6个 杂株基因型频率被用于计算杂株率R，它们对应的测试区域的测序片段总数为1 38586条，所以ΣSeN＝138586，也即相当于对抽得的26000个样本再进行了138 586次抽样，如此大的抽样量的误差是相当小的。本实施例中，判定待测大豆品 种不具有一致性和稳定性，因此，该判定结论正确的概率≥BINOM.DIST((1-0.5％) *26000,26000,(1-1.03％),TRUE)*BINOM.DIST(138586*(1-0.5％),138586,1-1.03％, TRUE)＝100.0000％。可见，本实施对待测大豆品种的一致性与稳定性的判定也 是很准确的。

结果验证

按《植物新品种特异性、一致性和稳定性测试指南-大豆》中的方法种植并 观察待测大豆品种及其近似品种“垦鉴豆26号”，发现待测大豆品种在叶色等多 个性状上与近似品种存在明显差异。《植物新品种特异性、一致性和稳定性测试 指南-大豆》中规定：至少在一个性状上与近似品种具有明显且可重现的差异时， 即可判定申请的待测大豆品种具备特异性。因此，判定待测大豆品种具有特异 性。在实验过程中，共种植了300株待测大豆品种与近似品种(150株一个小区， 共2个重复)，发现7株异形株，《植物新品种特异性、一致性和稳定性测试指 南-大豆》中规定：当样本大小为300株时，最多可以允许有4株异型株，由 此判定待测大豆品种不具有一致性。由于不具备一致性的品种可认为该品种不 具备稳定性。由此判定，待测大豆品种也不具有稳定性。通过以上实验表明： 本实施例中对待测大豆品种的特异性、稳定性与一致性的判定是正确的。

本发明实施例通过高通量测序和多位点扩增，实现了待测大豆品种的大样 本抽样与种间个体测试区域的大样本抽样，再利用定义杂株基因型、定义细胞 质杂株品种和定义杂株率计算公式等综合手段，成功地实现了准确、快速、完 整地判断待测大豆品种的特异性、稳定性与一致性的目标，其技术效果是已有 DUS测试方法都达不到的。现有分子DUS检测技术如芯片只检测固定测试区域， 不能根据个案，灵活选择非通用测试区域。而本发明检测的是PCR产物，可以 很方便地根据个案灵活设计引物，检测非通用测试区域。此外，本发明实施例 针对26000个个体的抽样量对于传统的DUS测试技术来说，工作大，无法完成， 例如，田间DUS测试中，抽样26000株大豆需要种植2亩以上，且需种植2年， 且每年每株大豆需调查多个性状。在广泛采用的SSR分子DUS测试中，需要分 别做26000次DNA提取，26000*2488次PCR和26000*2488次PCR产物检测 (假设与本实施例一样，检测了2488个通用测试区域)。因此，由于工作量过 大，已有分子DUS测试都没测试稳定性与一致性，田间DUS测试虽然检测一 致性和稳定性，但抽样样本量都在1000株以下，而本实施例抽样了26000株大 豆，其准确性显然更高。本实施例之所以可以加大抽样量，是因为所有26000 个样本都混合后作为一个样本处理，与田间DUS测试比较，工作量相当于缩减 为1/26000；进一步地，所有2488个通用测试区域都只做一次混合扩增和一次 高通量测序检测，与SSR分子DUS测试比较，工作量相当于缩减为1/ (26000*2488)。因此，本发明在工作量显著减轻的情况下，实现了大样本和多 位点检测，使DUS测试不但准确而且简单。同时本发明实施例中数据库品种基 因型为碱基组成，十分标准，在不同实验条件下按本发明的方法检测相同品种， 可得到完全相同的基因型，因而，不需要在不同的条件下重复DUS测试，因此， 本发明实施例可直接与数据库品种基因型比较，客观地选择待测大豆品种的近 似品种。而已有的DUS测试技术不够标准，需要同时平行地对待测大豆品种与 近似品种进行DUS测试，才能得到可靠的结论，为了减轻工作量，不得不由由 品种权申请者提供近似品种，若近似品种错误，则可能产生错误授权的法律后 果。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的 精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的 保护范围之内。