CN104805183A - 一种测试纯系植物新品种的特异性、一致性与稳定性的方法 - Google Patents

Abstract

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

Description

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

技术领域

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

背景技术

作为一种特化的知识产权，植物新品种已经成为一个公司及至一个国家的核心竞争力。植物新品种授权与相关法律问题的解决依赖于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满足如下条件：BINOMDIST(10,10,BINOMDIST(8,20,BINOMDIST(0,CF,0.1％,TRUE),TRUE),FALSE)≥99.9％，1-BINOMDIST(10000,10000,1-BINOMDIST(8,20,1-BINOMDIST(99.99％*CF,CF,99.9989％,TRUE),TRUE),FALSE)≤0.1％且BINOMDIST(10*(1-M)*CF,10*CF,1-110％*M,TRUE)≥95.0％，其中，CF为所述高通量测序的深度，M为判断所述一致性和稳定性时所选用的阈值，BINOMDIST为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时，所述待测品种不具有一致性和稳定性；

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

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

$R2=\frac{{\mathrm{\Σ}}_{i2=\mathrm{Int}(0.8\×t2)+1}^{t2-\mathrm{Int}(0.1\×t2)}2\×R2i2}{t2-\mathrm{Int}(0.8\×t2)-\mathrm{Int}(0.1\×t2)},$ 其中，t2为除所述细胞核杂株品种拥有的杂株核基因型之外的且频率≥0.17％的所述杂株核基因型的数目，i2为除所述细胞核杂株品种拥有的所述杂株核基因型之外的所有所述杂株核基因型按频率由低到高排序后，第i2个所述杂株核基因型，R2i2为第i2个所述杂株核基因型的频率；R2是利用除所述细胞核杂株品种拥有的所述杂株核基因型计算的所述杂株率，R2为去掉除所述细胞核杂株品种拥有的所述杂株核基因型的频率中最低的80％和最高的10％的值后，剩余值的平均值的2倍；

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

$R4=\frac{{\mathrm{\Σ}}_{i4=\mathrm{Int}(0.8\×t4)+1}^{t4-\mathrm{Int}(0.1\×t4)}R4i4}{t4-\mathrm{Int}(0.8\×t4)-\mathrm{Int}(0.1\×t4)},$ 其中，t4为除所述细胞质杂株品种拥有的所述杂株质基因型之外的且频率≥0.17％的所述杂株质基因型的数目，i4为除所述细胞质杂株品种拥有的所述杂株质基因型之外的所有所述杂株质基因型按其频率由低到高排序后，第i4个所述杂株质基因型，R4i4为第i4个所述杂株质基因型的频率；R4是利用除所述细胞质杂株品种拥有的所述杂株质基因型计算的杂株率，R4为去掉除所述细胞质杂株品种拥有的所述杂株质基因型的频率中最低的80％和最高的10％的值后，剩余值的平均值；

Int()为取整函数；

所述细胞核杂株品种是指仅利用核基因型计算获得的所述杂株品种，所述细胞质杂株品种是指仅利用质基因型计算获得的所述杂株品种；所述特异杂株核基因型是指仅为一个所述细胞核杂株品种所有的所述杂株核基因型；所述特异杂株质基因型是指仅为一个所述细胞质杂株品种所有的所述杂株质基因型；所述杂株核基因型是指所述杂株基因型为所述核基因型，所述核基因型指所述基因型位于细胞核基因组上；所述杂株质基因型是指所述杂株基因型为所述质基因型，所述质基因型是指所述基因型位于细胞质基因组上。

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

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

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

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

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

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

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

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

具体实施方式

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

实施例.测定水稻新品种‘R7723’的特异性、一致性与稳定性

本发明实施例提供的待测品种为水稻品种“R7723”，水稻品种“R7723”为纯系水稻且为公开授权品种，授权号为CNA20100474.1。水稻品种“R7723”的来源为水稻品种“R8377”与“IRBB23”杂交后，以“R8377”为轮回亲本，以抗白叶枯病性状为主要目标性状，按回交育种的方法培育而成。测定该水稻品种的特异性、一致性与稳定性的方法包括以下步骤。

一、获得待测品种所属种内不同品种间的变异位点。

待测品种的种为水稻，因此，获得待测品种所属种内不同品种间的变异位点则应为获取不同水稻品种间的变异位点，不同水稻品种间的变异位点可以从已公布的文献资料中获取，但该方法所获得的结果比较零星，在本实施例中，通过将不同水稻的基因组序列与参考水稻品种的基因组序列进行比对，获得了大量的不同水稻品种间的变异位点，其中参考水稻品种可以为“日本睛”水稻，该“日本睛”水稻可以替换为其他已知的参考水稻品种。

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

本实施例的不同水稻品种的基因组序列有三种来源，第一种为韩斌对1082个水稻品种的基因组的高通量测序序列，相关文献信息如下：Huang XH et al.Amap of rice genome variation reveals the origin of cultivated rice.Nature.2012；7:497–503。该1082个水稻品种的基因组序列公布于EuropeanNucleotideArchive(http://www.ebi.ac.uk/ena/)，接收号为ERP001143、ERP000729和ERP000106；第二种为徐讯对50个水稻品种的基因组的高通量测序序列，相关文献信息如下：Xun X et al.Resequencing 50accessions of cultivated and wildrice yields markers for identifying agronomically important genes.Nat Biotechnol.2011,30(1):105-11，该50个水稻品种的基因组序列公布于NCBI Short ReadArchive(http://www.ncbi.nlm.nih.gov/sra)，接收号为SRA023116；第三种为按韩斌上述发表的文章中提供的方法对“R8377”、“金科1A”、“IRBB23”、杂交种“金科1A/R7723”和杂株种“D优527”进行了高通量测序。本实施例共获得了1137个水稻品种的基因组的高通量测序序列。

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

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

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

确定通用测试区域

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

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

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

确定非通用测试区域

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

本实施例中，通过回交育种的方式，将高抗白叶枯的基因Xa23从亲本IRBB23中转入了亲本R8377后培育了本实施例提供的待测品种。因此，待测品种与亲本R8377遗传背景相同，仅白叶枯抗性不同。Xa23基因控制的白叶枯病抗性为质量性状，且Xa23来源于野生稻，不包括在通用测试区域中。基于以上理由，将Xa23基因作为非通用测试区域加以检测，Xa23基因已被克隆，其抗性由7个碱基的缺失引起，因此，待测品种的非通用测试区域即为这7个缺失的碱基，其位于日本睛参考基因组上第11染色体的24046820到24046825位，关于Xa23基因的更详细信息见：Wang,C.,X.Zhang,et al.(2014)."XA23is anexecutor R protein and confers broad-spectrum disease resistance in rice."Molecularplant:ssu132.

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

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

通用测试区域采用多重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引物。本实施例中，从上述获得的4261个通用测试区域中，设计并成功验证了2231对多重PCR引物，用于扩增相应的2231个通用测试区域。验证多重PCR引物的方法为按本发明提供的方法，提取同一株水稻上的叶片基因组DNA，并利用设计的多重PCR引物对获得的基因组DNA进行扩增、建库、高通量测序并分析测序片段组，去掉以下测试区域相应的引物：该测试区域的测序片段数不足1000或存在杂株基因型，保留下来的引物即为验证成功的多重PCR引物。由于基因组DNA来源于同一株水稻叶片，不可能存在杂株品种，因此，杂株基因型是由测试区域的特殊结构造成的PCR或测序偏好性错误，去掉这些测试区域避免了此类系统错误。验证成功的多重PCR引物也由该公司混合好后以液体的形式提供给客户使用。上述成功设计了多重PCR引物的2231个通用测试区域即为最终用于待测品种检测的通用测试区域，同时，构建的数据库中的每个品种也包含了上述2231个通用测试区域，其中，100个通用测试区域位于细胞质基因组上，剩余的2131个通用测试区域位于细胞核基因组上。

需要说明的是：通用测试区域的数目要求≥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 IDNO:2，分别获得第一引物的正向引物和第一引物的反向引物；第三步，将设计好的正向引物与反向引物的5’端分别连接序列表中SEQ ID NO:2和序列表中SEQ ID NO:1，分别获得第二引物的正向引物和第二引物的反向引物。序列表中SEQ ID NO:1和序列表中SEQ ID NO:2为高通量测序所用的接头序列，从而使用PCR产物带有高通量测序的接头序列，可以直接与扩增的通用测序区域的产物混合后建立测序文库后一同测序，而不必经过片段化、连接接头等繁琐的建库步骤，提高了工作效率并降低了成本。做成两对仅接头不一样的引物是为了同时从非通用测试区域的两端测序。

具体地，在本实施例中，所设计用于扩增待测品种非通用测试区域(Xa23基因)的普通PCR引物的正向引物序列为：TGCGGCATCACTAACATCAG，反向引物序列为：TGTTAGTGATGCGGGAGGAA。向其两端分别加上序列表中SEQ ID NO:1和序列表中SEQ ID NO:2后形成的第一引物的正向引物为：5’-CCATCTCATCCCTGCGTGTCTCCGACTCAGTGCGGCATCACTAACATCAG如序列表中SEQ ID NO:3；第一引物的反向引物为：5’-CCTCTCTATGGGCAGTCGGTGATTGTTAGTGATGCGGGAGGAA如序列表中SEQ ID NO:4；第二引物的正向引物为：5’-CCTCTCTATGGGCAGTCGGTGATTGCGGCATCACTAACATCAG如序列表中SEQ ID NO:5；第二引物的反向引物为：5’-CCATCTCATCCCTGCGTGTCTCCGACTCAGTGTTAGTGATGCGGGAGGAA如序列表中SEQID NO:6。所设计的非通用测试区域引物由美国LifeTechnology公司合成。

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

构建包含不同品种的在所有测试区域的基因型的数据库，具体地，在待测品种的测试区域上，获得不同品种在对应该测试区域上的基因型并组成数据库。本实例获得了2231个通用测试区域引物和1个非通用测试区域引物，它们对应的扩增区域即为待测品种的测试区域。构建包含1137个品种的2232测试区域的基因型及其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，其中，BINOM.INV为excel 2010中的函数，其使用方法与excel 2010中的定义相同，其含义是使得累积二项分布的函数值大于或等于临界值的最小整数。抽样量SN满足的条件含义为：即使杂株率只超出阈值M的15％，该抽样量在95％的概率保证下，可正确判断待测品种的稳定性与一致性。M值是根据作物种类、标记类型、具体要求等条件人为确定的。在农业部新品种保护办公室发布中《植物新品种特异性、一致性和稳定性测试指南-水稻》中规定：自交种样本大小为356～818株时，最多可以允许有2株异型株，相当于M值为0.24％～0.56％，因此，本实施例中，选用中间值0.40％作为M值。逐步加大SN值后，计算上述公式发现，当SN≥29783时，BINOM.INV(SN,0.40％,0.95)/SN≤1.15*0.40％成立。因此，本实施例中的待测样本抽样量应≥29783。

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

在本实施例中，选取了50000粒种子发芽，随机选取30000个大小大致相等的芽混合后置于研钵中，向研钵中加入液氮后充分研磨成粉。采用北京天根生化科技有限公司生产的货号为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连接酶。文库构建的方法按该试剂盒的操作手册《IonAmpliSeq^TM Library Preparation》(出版号：MAN0006735，版本：A.0)进行。通过多重PCR扩增2231个通用测试区域，多重PCR的扩增体系如下：5×IonAmpliSeq^TM HiFi Mix 4μl、制备的通用测试区域引物混合液4μl、待测品种的DNA10ng和无酶水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水中，利用DNA1000试剂盒(货号为5067-1504)在美国安捷伦公司生产的生物分析仪(型号为2100)上，按该试剂盒说明书测定并获得扩增产物的摩尔浓度后，稀释为200pM，即为第一引物的扩增产物。采用相同的方法，获得浓度为200pM的第二引物的扩增产物。将第一引物的扩增产物与第二引物的扩增产物等体积混合，获得浓度为100pM的非通用测试区域高通量测序文库。

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

按通用测试区域的数目与非通用测试区域的数目的比例混合等摩尔浓度的通用测试区域的高通量测序文库和非通用测试区域的高通量测序文库，得到的混合物即为所有测试区域的高通量测序文库。在本实施例中，取获得的通用测试区域的高通量测序文库2231μL和1μL非通用测试区域的高通量测序文库混匀后，得到浓度为100pM的所有测试区域的高通量测序文库。

七、对高通量测序文库进行高通量测序，得到测序片段组。

确定高通量测序深度的原理：高通量测序的深度满足如下条件：BINOMDIST(10,10,BINOMDIST(8,20,BINOMDIST(0,CF,0.1％,TRUE),TRUE),FALSE)≥99.9％，1-BINOMDIST(10000,10000,1-BINOMDIST(8,20,1-BINOMDIST(99.99％*CF,CF,99.9989％,TRUE),TRUE),FALSE)≤0.1％且BINOMDIST(10*(1-M)*CF,10*CF,1-110％*M,TRUE)≥95.0％，其中，CF为高通量测序的深度，M为判断一致性和稳定性时所选用的阈值，也即平均每个测试区域被覆盖的倍数，BINOMDIST为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)”为单元格B2，其它单元格编号以B2为参照，按Excel 2010的规则定义，例如“杂株率(R)”所在的单元格在B2的基础上增加了4行1列，因此编号为C6，其它单元格编号规则与此相同。

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

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

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

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

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

八、分析测序片段组，获得检测成功的测试区域、待测品种基因型和杂株基因型，方法如下；

将测序片段组比对到所有测试区域，并统计每个测试区域中的测序片段数，去掉测序片段数≤1000条的测试区域，剩余的测试区域为检测成功的测试区域。在本实施例中，共获得2030个检测成功的测试区域。比对到测试区域的片段称为该测试区域的测序片段，从测序片段中提取表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个测序区域的测序片段总数为33320条，有ACCC、CGTT、CCCC、GCCC……共42种基因型，代表这些基因型的测序片段数分别33001条、16条、1条、2条……，这些基因型的频率为33301/33320＝99.04％、16/33320＝0.05％、1/33320＝0.003％、2/33320＝0.006％……。按待测品种基因型和杂株基因型的定义，ACCC应该是待测品种在第1个测试区域的待测品种基因型，而CGTT的频率超过0.02％且与待测品种基因型ACCC比较有4个≥2个碱基的差异，因此CGTT为杂株基因型，其它基因型为测序错误产生的基因型。杂株核基因型是指杂株基因型为核基因型，杂株质基因型是指杂株基因型为质基因型。按此定义，第一个测试区域的杂株基因型CGTT也为杂株核基因型。按相同的方法，判断并获得全部2030个检测成功的测试区域的待测品种基因型、杂株基因型及其频率，并判断获得的杂株基因型是杂株核基因型还是杂株质基因型。结果表明：共获得109个杂株基因型，其中，108个为杂株核基因型，1个为杂株质基因型。

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

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

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

在本实施例中，待测品种与数据库的第1个品种“金科1A”的共有测试区域数为2025个。在第1个共有测试区域中，待测品种与“金科1A”基因型分别为ACCC和CGTT，二者不完全相同，因此，第1个共有测试区域为待测品种与“金科1A”的差异位点，CGTT与ACCC为待测品种与“金科1A”的差异基因型。按相同的方法，将所有共有测试区域中，待测品种与“金科1A”基因型比较，发现共有153个差异位点，差异位点率＝153/2025＝7.56％。按相同的方法，获得待测品种与数据库中所有1137个品种间的差异位点率，并获得差异位点率最小的品种为“R8377”，仅1个差异位点，其为编号10号的非通用测试区域(见表1)，差异位点率为0.05％。因此，“R8377”为待测品种的近似品种，待测品种的变异位点率为0.05％。

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

获得杂株品种：杂株品种存在于数据库中的品种，且杂株品种的潜在杂株基因型与杂株基因型间有相同基因型的测试区域的数目占杂株品种具有潜在杂株基因型的测试区域的总数的比例≥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个测试区域中，数据库中第一个品种“金科1A”与待测品种的基因型分别为GGTT和ACCC，二者间存在2个以上的碱基的差异，因此，为GGTT为潜在杂株基因型，而该潜在杂株基因型与第1个测试区域中的杂株基因型GGTT相同，按相同的方法，逐个判断所有核基因型的测试区域中，数据库中第一个品种“金科1A”的基因型是否为潜在杂株基因型，若为潜在杂株基因型，再判断潜在杂株基因型与杂株基因型间是否有相同基因型，结果表明，“金科1A”共有97个具有潜在杂株基因型的测试区域，它们全部与相同测试区域的杂株基因型间有相同基因型，其比例为97/97＝100％>60％，因此，判断“金科1A”为细胞核杂株品种，按类似的方法，利用所有质基因型的测试区域，判断“金科1A”不为细胞质杂株品种。按相同的方法，判断数据库中所有其它品种是否为细胞核杂株品种或细胞质杂株品种，结果表明：仅“金科1A”为细胞核杂株品种，没有发现细胞质杂株品种。以上结果说明：“金科1A”是通过飞花传粉而不是机械混杂，将基因型混入了待测品种，由于“金科1A”可育花粉很少，而“金科1A”与“金科1B”核基因型完全相同，因此，混入的花粉是与“金科1A”核基因型完全相同的“金科1B”的花粉。

获得特异杂株基因型：特异杂株基因型是指仅为一个杂株品种所有的杂株基因型，其包括特异杂株核基因型和特异杂株质基因型；特异杂株核基因型是指仅为一个细胞核杂株品种所有的杂株核基因型，特异杂株质基因型是指仅为一个细胞质杂株品种所有的杂株质基因型。本实施例中，共获得97个杂株基因型，且全部为杂株核基因型，第一个杂株核基因型CGTT仅为细胞核杂株品种“金科1A”所有，所以，CGTT为“金科1A”的特异杂株核基因型。按相同的方法，逐一判断所有获得的97个杂株基因型中，全部为“金科1A”所拥有的特异杂株核基因型。本实施例中，由于97个杂株基因型中，无质基因型，所以，也无特异杂株质基因型。

计算杂株率R原理

杂株率R＝R1+R2-R3-R4，其中： $R1={\mathrm{\&Sigma;}}_{i1=1}^{n1}\frac{{\mathrm{\&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倍； $R2=\frac{{\mathrm{\&Sigma;}}_{i2=\mathrm{Int}(0.8\&times;t2)+1}^{t2-\mathrm{Int}(0.1\&times;t2)}2\&times;R2i2}{t2-\mathrm{Int}(0.8\&times;t2)-\mathrm{Int}(0.1\&times;t2)},$ 其中，t2为除细胞核杂株品种拥有的杂株核基因型之外的且频率≥0.17％的杂株核基因型的数目，i2为除细胞核杂株品种拥有的杂株核基因型之外的所有杂株核基因型按其频率由低到高排序后，第i2个杂株核基因型，R2i2为第i2个杂株核基因型的频率；R2是利用除细胞核杂株品种拥有的杂株核基因型计算的杂株率，其为去掉除细胞核杂株品种拥有的杂株核基因型的频率中最低的80％和最高的10％的值后，剩余值的平均值的2倍； $R3={\mathrm{\&Sigma;}}_{i3=1}^{n2}R3i3-2\&times;R3\mathrm{ic},$ 其中， $R3i3=\frac{{\mathrm{\&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％的特异杂株质基因型的频率后，剩余的特异杂株质基因型的频率的平均值； $R4=\frac{{\mathrm{\&Sigma;}}_{i4=\mathrm{Int}(0.8\&times;t4)+1}^{t4-\mathrm{Int}(0.1\&times;t4)}R4i4}{t4-\mathrm{Int}(0.8\&times;t4)-\mathrm{Int}(0.1\&times;t4)},$ 其中，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={\mathrm{\&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：在本实施例中，杂株品种仅为“金科1A”且为细胞核杂株品种，R2、R3与R4均为0，因而，R＝R1＝R111。“金科1A”共有97个特异杂株核基因型，频率为：0.05％、0.05％、0.06％、0.05％......(部分实例见表1)，按R的计算规则，去除最小的80％(77个)和最小的10％(9个)的频率值后，剩余的11个频率的平均值即为杂株率R＝0.05％。

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

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

在本实施例中，变异位点率为0.05％<SD＝1％，但非通用测试区域(编号为10号的测试区域)存在变异位点(见表1)，因此，判断待测品种具有特异性；待测品种的杂株率0.05％≤M＝0.40％，因此，判断待测品种具有一致性和稳定性。

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

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

本实施未采用变异位点率来判断待测品种的特异性，因此，不计算特异性结论正确的概率。

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

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

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

结果验证

按《植物新品种特异性、一致性和稳定性测试指南-水稻》中的方法种植并观察待测水稻品种及其近似品种“R8377”，发现待测水稻品种的高抗白叶枯病，近似品种则高感。《植物新品种特异性、一致性和稳定性测试指南-水稻》中规定：至少在一个性状上与近似品种具有明显且可重现的差异时，即可判定申请的待测品种具备特异性。因此，判定待测水稻品种具有特异性。在实验过程中，共种植了400株待测品种与近似品种(200株一个小区，共2个重复)，发现1株异形株，待测水稻品种为自交种，《植物新品种特异性、一致性和稳定性测试指南-水稻》中规定：自交种样本大小为356～818株时，最多可以允许有2株异型株，由此判定待测品种具有一致性。《植物新品种特异性、一致性和稳定性测试指南-水稻》中规定：如果一个品种具备一致性，则可认为该品种具备稳定性。由此判定，待测品种也具有稳定性。待测品种目前已经通过中华人民共和国农业部新品种保护办公室的DUS测试，并获得了授权，授权号为CNA20100474.1，待测品种可以获得授权表明本实施例中对待测品种的特异性、稳定性与一致性的判定是正确的。

实施例二.测定油菜品种‘苏2051’的特异性、一致性与稳定性

本发明实施例提供的待测油菜品种为油菜品种“苏2051”，油菜品种“苏2051”为纯系油菜且为公开使用品种。则测定该油菜品种的特异性、一致性与稳定性的方法包括以下步骤。由于本实施例中的方法与实施例一中的方法类似，因此，相同的地方就不再重述，只列出不同的地方。

一、获得不同油菜品种间的变异位点。

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

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

本实施例的不同油菜品种的基因组序列有两种来源，第一种为Huang等对10个油菜品种的基因组的高通量测序序列，相关文献信息如下：Huang et al.:Identification of genome-wide single nucleotide polymorphisms in allopolyploid cropBrassica napus.BMC Genomics 2013 14:717。该10个油菜品种的基因组序列公布于NCBI Short Read Archive(http://www.ncbi.nlm.nih.gov/sra)，接收号为SRA057227；第二种为按Huang等的上述发表的文章中提供的方法对“430AB”、“P65”和杂交种“宁杂9号”进行了高通量测序。本实施例共获得了13个油菜品种的基因组的高通量测序序列。

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

按与实施例一相同的方法将这13个油菜品种的基因组的高通量测序序列分别比对到油菜细胞核参考基因组(版本：Release v1.01，下载地址：http://www.ncbi.nlm.nih.gov)和细胞质参考基因组上，该细胞质参考基因组包括线粒体参考基因组与叶绿体参考基因组，其在NCBI(National Center forBiotechnology Information，美国国立生物技术信息中心)上的接收号分别为NC_016734.1和AP006444.1。对比的参数与方法与实施例一相同，本实施例在所有13个油菜品种间共获得911346个SNP位点，其中18543个SNP位点位于细胞质基因组上，其余的SNP位点位于细胞核基因组上。

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

确定通用测试区域

按与实施例一相同的方法，获得的全部911346个变异窗口的区分度，并从中选取位于细胞核基因组中区分度最大的6000个变异窗口和位于细胞质基因组中区分度最大的100个变异窗口。逐个检查位于细胞核基因组的6000个变异窗口中，每个变异窗口与下一个变异窗口间的距离，若距离超过200K(1K＝1000个碱基)，则放弃其中区分度较小的变异窗口之后再检查，直至相邻查变异窗口的距离均大于200K为止。选择200K的距离标准是因为油菜基因组大小约为930M(1M＝100万个碱基)，按最终入选2000个位于细胞核基因组的通用测试区域计，平均的通用测试区域间距离约为500K，但由于一些特异区域如着丝粒等很少有变异位点，因此，平均距离应该小于500K。按以上方法，选出了4367个位于细胞核基因组的变异窗口，它们与获得的位于细胞质基因组中区分度最大的100个变异窗口一起共4467个变异窗口作为入选的通过测试区域。

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

按与实施例一相同的方法，在本实施例中，从所有4467个通用测试区域中，设计并成功验证了2302对多重PCR引物，用于扩增相应的2302个通用测试区域。

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

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

按与实施例一相同的方法，构建包含13个品种的2302测试区域的基因型及其SNP的位置信息的数据库，部分结果见表4。

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

表4中的符号含义与表1中的符号含义相同。

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

按与实施例一相同的方法，计算待测油菜品种抽样量应≥5957。在本实施例中，选取了10000粒种子发芽，随机选取8000个大小大致相等的芽混合后按与实施例一相同的方法提取待测油菜品种的DNA后，同样稀释至10.00ng/μl。

六、利用引物对混合样本的DNA进行扩增，得到测试区域的扩增产物，扩增产物用于构建高通量测序文库，其中：

按与实施例一相同的方法，利用引物对混合样本的DNA进行扩增，得到测试区域的扩增产物，扩增产物用于构建高通量测序文库。

获得高通量测序文库：在本实施例中，因无非通用测试区域的高通量测序文库，因此，构建的高通量测序文库即为浓度为100pM的通用测试区域的高通量测序文库。

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

按与实施例一相同的方法，确实本实施例测序深度为≥1935倍。按与实施例一相同的方法，利用高通量测序文库进行高通量测序，高通量测序通量设置为平均覆盖测试区域10000倍。按与实施例一相同的方法，对高通量测序结果进行预处理。

八、按与实施例一相同的方法，分析测序片段组，获得待测油菜品种基因型和杂株基因型，方法如下：

在本实施例中，共获得2117个检测成功的测试区域。判断并获得全部2117个检测成功的测试区域的待测油菜品种基因型、杂株基因型及其频率，并判断获得的杂株基因型是杂株核基因型还是杂株质基因型。结果表明：共获得177个杂株基因型，其中，174个为杂株核基因型，3个为杂株质基因型。

按与实施例一相同的方法，对本实施例中的标准样品进行检测方法，结果为：从获得的177个杂株基因型中，共去掉了2个杂株基因型，其中2个为杂株核基因型，0个为杂株质基因型，所保留下来的175个杂株基因型用于后续分析，部分结果见表4。

九、按与实施例一相同的方法，将待测油菜品种基因型与数据库中的不同品种的基因型比较，获得近似品种、变异位点和变异位点率，其中，获得差异位点率最小的品种为“P65”，差异位点率为3.68％。因此，“P65”为待测油菜品种的近似品种，待测油菜品种的变异位点率为3.68％。

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

按与实施例一相同的方法，获得仅“430AB”为细胞核杂株品种，没有发现细胞质杂株品种。按与实施例一相同的方法，共获得177个杂株基因型，其中，174个为杂株核基因型，3个为杂株质基因型。逐一判断所有获得的177个杂株基因型中，82个为“430AB”所拥有的特异杂株核基因型。判断10个杂株质基因型均不为特异杂株质基因型。

按与实施例一相同的方法，计算本实施例中的杂株率R＝2.52％。

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

在本实施例中，变异位点率为3.68％>SD＝1％，因此，判断待测油菜品种具有特异性；待测油菜品种的杂株率2.52％＞M＝2％，因此，判断待测油菜品种不具有一致性和稳定性。

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

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

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

按与实施例一相同的方法，计算一致性与稳定性的正确率。本实施例中，判定待测油菜品种不具有一致性和稳定性，因此，该判定结论正确的概率≥BINOM.DIST((1-M)*SN,SN,(1-R),TRUE)*BINOM.DIST(ΣSeN*(1-M),ΣSeN,1-R,TRUE)＝BINOM.DIST((1-2％)*8000,8000,(1-2.52％),TRUE)*BINOM.DIST(89091*(1-2％),89091,1-2.52％,TRUE)＝99.90％。可见，本实施对待测油菜品种的一致性与稳定性的判定也是很准确的。

结果验证

按《植物新品种特异性、一致性和稳定性测试指南-甘蓝型油菜》中的方法种植并观察待测油菜品种及其近似品种“430AB”，发现待测油菜品种在株高等多个性状上与近似品种存在明显差异。《植物新品种特异性、一致性和稳定性测试指南-甘蓝型油菜》中规定：至少在一个性状上与近似品种具有明显且可重现的差异时，即可判定申请的待测油菜品种具备特异性。因此，判定待测油菜品种具有特异性。在实验过程中，共种植了200株待测油菜品种与近似品种(100株一个小区，共2个重复)，发现12株异形株，《植物新品种特异性、一致性和稳定性测试指南-甘蓝型油菜》中规定：观测样本为200株时，最多可以允许有7株异型株，由此判定待测油菜品种不具有一致性。由于待测油菜品种已经不具有一致性，因此，也不可能具有稳定性。由此判定，待测油菜品种也不具有稳定性。通过以上实验表明：本实施例中对待测油菜品种的特异性、稳定性与一致性的判定是正确的。

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

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

Claims (10)

1.一种测试纯系植物新品种的特异性、一致性与稳定性的方法，其特征在于，所述方法包括：

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

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

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

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

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

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

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

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

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

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

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

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

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

4.根据权利要求1所述的方法，其特征在于，所述测试区域还包括非通用测试区域，所述引物还包括非通用测试区域引物。

5.根据权利要求4所述的方法，其特征在于，所述非通用测试区域引物包括第一引物和第二引物，所述第一引物包括第一正向引物和第一反向引物，所述第二引物包括第二正向引物和第二反向引物，所述第一引物和所述第二引物分别进行单独扩增得到两个所述非通用测试区域的扩增产物，将两个所述非通用测试区域的扩增产物等量混合用于构建单独扩增的高通量测序文库；

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

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

6.根据权利要求4所述的方法，其特征在于，利用所述变异位点、所述变异位点率和所述杂株率，判断所述待测品种特异性、一致性和稳定性的方法包括：

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

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

所述杂株率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是利用除所述细胞核杂株品种拥有的所述杂株核基因型计算的所述杂株率，R2为去掉除所述细胞核杂株品种拥有的所述杂株核基因型的频率中最低的80％和最高的10％的值后，剩余值的平均值的2倍；

$R3={\mathrm{\&Sigma;}}_{i3=1}^{n2}R3i3-2\&times;R3\mathrm{ic},$ 其中， $R3i3=\frac{{\mathrm{\&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是利用除所述细胞质杂株品种拥有的所述杂株质基因型计算的所述杂株率，R4为去掉除所述细胞质杂株品种拥有的所述杂株质基因型的频率中最低的80％和最高的10％的值后，剩余值的平均值；

Int()为取整函数；

所述细胞核杂株品种是指仅利用核基因型计算获得的所述杂株品种，所述细胞质杂株品种是指仅利用质基因型计算获得的所述杂株品种；所述特异杂株核基因型是指仅为一个所述细胞核杂株品种所有的所述杂株核基因型；所述特异杂株质基因型是指仅为一个所述细胞质杂株品种所有的所述杂株质基因型；所述杂株核基因型是指所述杂株基因型为所述核基因型，所述核基因型指所述基因型位于细胞核基因组上；所述杂株质基因型是指所述杂株基因型为所述质基因型，所述质基因型是指所述基因型位于细胞质基因组上。

7.根据权利要求6所述的方法，其特征在于，所述方法还包括采用以下方式判定所述待测品种的一致性和稳定性的结论正确的概率为：当所述待测品种具有一致性和稳定性时，结论正确的概率≥BINOMDIST(M*SN,SN,R,TRUE)*BINOMDIST(∑SeN*M,∑SeN,R,TRUE)；当所述待测品种不具有所述一致性和稳定性时，结论正确的概率≥BINOMDIST((1-M)*SN,SN,(1-R),TRUE)*BINOMDIST(∑SeN*(1-M),∑SeN,1-R,TRUE)；其中，M为判断所述一致性和稳定性时所选用的阈值，∑SeN为所有用于计算所述杂株率R的基因型的频率所在所述测试区域的测序片段的总和，BINOMDIST(∑SeN*M,∑SeN,R,TRUE)的意义为：对所述待测品种进行了∑SeN次抽样，实际抽得的所述杂株率R小于阈值M的概率；BINOMDIST((1-M)*SN,SN,(1-R),TRUE)为所述待测品种进行了SN次抽样，实际抽得的所述杂株率R大于所述阈值M的概率，所述基因型的频率是指所述测序片段组中，代表所述基因型的测序片段数占所述基因型所在所述测试区域的测序片段总数的比例。

8.根据权利要求6所述的方法，其特征在于，当所述非通用测试区域不存在所述变异位点时，若判断所述待测品种具有特异性，结论正确的概率≥BINOMDIST((1-SD)*TRN,TRN,1-OD,TRUE)；若判断所述待测品种不具有特异性，结论正确的概率≥BINOMDIST(SD*TRN,TRN,OD,TRUE)，其中，TRN为检测成功的测试区域的数目，OD为所述变异位点率，SD为判断特异性时所选用的阈值，BINOMDIST为excel 2010中的函数，所述结论正确的概率表示为当判断所述待测品种具有特异性时，所述变异位点率大于SD的概率，当判断所述待测品种不具有特异性时，所述变异位点率小于SD的概率，所述检测成功的测试区域通过分析所述测序片段组后获得。

9.根据权利要求1所述的方法，其特征在于，获得所述杂株品种的方法包括：所述杂株品种为存在于所述数据库中的品种，且所述杂株品种的潜在杂株基因型与所述杂株基因型间有相同基因型的所述测试区域的数目占所述杂株品种具有所述潜在杂株基因型的所述测试区域的总数的比例≥60％；所述杂株基因型指频率≥0.02％的所述潜在杂株基因型；

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

10.根据权利要求1所述的方法，其特征在于，通过所述变异位点确定所述通用测试区域的方法为：

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

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