CN105404793A - 基于概率框架和重测序技术快速发现表型相关基因的方法 - Google Patents

Abstract

本发明公开了一种基于概率框架和重测序技术快速发现表型相关基因的方法，建立了一种基于概率框架的方法，对基于基因组重测序技术的正向遗传学研究过程建模，通过计算四个指标评估研究过程中每一个步骤的设计对研究整体有效性造成的影响，从而指导对当前实验方案和分析方法的优化，达到使用尽可能少的样本，快速发现可能与特定表型相关的基因的目的。首次提出用四个指标来评估基于基因组重测序技术的正向遗传学研究过程的整体有效性，其中方法成功率和非孟德尔表型显著性是创新提出的指标，对于指导整个研究过程优化有重要价值。

Description

基于概率框架和重测序技术快速发现表型相关基因的方法

技术领域

本发明通过测序若干表型相同的独立个体的基因组，可快速发现表型相关的基因，具体地说，是一种基于概率框架和重测序技术快速发现表型相关基因的方法。

背景技术

测序技术飞速发展，在精度不断提高的同时成本不断降低。当前已普及的测序技术可以达到在一个月内完成一个人类个体的全基因组重测序，产生碱基准确度在99.9％以上测序数据对基因组平均30倍以上的覆盖；而其成本已经降到了一个普通的科研课题组可以经常使用的水平。

但是，基因组数据的大量产生并没有带来期望中的生物学、医学、药学、农学的爆炸式提升，主要原因是我们仍然缺乏基因与表型间关联关系的知识，无法应用基因组数据来对我们关注的诸多生物表型进行控制。为了获得基因与表型间的关联关系，我们可以使用正向遗传学方法和反向遗传学方法。反向遗传学通过对基因功能进行干预，观察产生的表型变化，从而获得特定基因与哪些表型相关的知识。正向遗传学通过研究一组具有特定表型特征的个体的基因组，获得特定表型与哪些基因相关的知识。本发明所阐述的方法，是一种基于基因组重测序技术的高效正向遗传学方法。

基因组重测序通常产生大量的短片段基因组读取数据。通过序列比对软件(如，BWA、Bowtie、Bowtie2、SOAP等)，可以把这些基因组片段序列定位到所研究物种的参考基因组上。然后，通过分析实测序列与参考基因组序列的异同，可以判断在基因组的特定位置上，实测个体是否存在相对于参考基因组的序列变异。序列变异分为两类，第一类是单核苷酸多态性(SNP)和小的插入删除(Smallindel)；第二类是结构变异，包括拷贝数变异(CNV)，大片段插入、删除、复制、倒向，和其它染色体结构变异。第一类变异可以由GATK、Samtools等软件分析获得。第二类变异可以由BreakDancer、DELLY、forestSV等软件分析获得。

正向遗传学方法测定具有特定表型特征的一组个体的基因组序列，把观察到的个体基因组变异注释到基因，然后判断哪些基因的功能受到影响与表型的产生相关。但是，由于个体基因组存在不同于参考基因组的遗传背景，个体基因组中存在较多的与所研究表型无关的基因组变异。这些变异中，一些对基因功能影响小，与所有表型无关；另一些虽然可以产生其它表型，但与所研究的表型无关。对于不能产生表型的变异，我们可以通过变异注释软件或变异对功能影响的预测软件来进行筛除，如SNPEff、ANNOVAR等，但是通过这些软件筛除存在一定的假阳性和假阴性率。对于能够产生其它与所研究表型无关表型的变异，可以依据对照样本，与表型样本遗传背景接近但不显示所研究表型的样本的基因组变异分析结果筛除。但是，对照样本的基因组变异分析也存在一定几率的假阳性和假阴性，也会导致筛除存在一定的假阳性和假阴性率。

因此，在进行了个体测序、参考基因组比对、变异分析、变异筛选这些步骤之后，对于每一个测序的个体，一方面与表型相关的基因组变异可能因为各种导致假阴性的原因而无法被观察到；另一方面与表型无关的基因组变异又可能因为各种导致假阳性的原因而被观察到。结果是，在最多的显示表型的个体中观察到变异的基因往往与表型无关。只有通过测序大量同一表型个体的基因组，才能有效发现表型相关的基因。但是，具有严格同一表型的样本往往获取困难。使用尽可能少的样本发现可能与特定表型相关的基因不但可以节约研究的时间和成本，包括验证基因与表型关联所需的时间和成本，还可以发现一些原先受制于样本量不足而无法发现的基因与表型间的关联关系。

发明内容

本发明正是针对现有技术的不足之处所作出的改进，建立了一种基于概率框架的方法，对基于基因组重测序技术的正向遗传学研究过程建模，通过计算四个指标评估研究过程中每一个步骤的设计对研究整体有效性造成的影响，从而指导对当前实验方案和分析方法的优化，达到使用尽可能少的样本，快速发现可能与特定表型相关的基因的目的。

本发明是通过以下技术方案来实现的：

本发明公开了一种基于概率框架和重测序技术快速发现表型相关基因的方法，对基于基因组重测序技术的正向遗传学研究过程建模，通过计算四个指标评估研究过程中每一个步骤的设计对研究整体有效性造成的影响，从而指导对当前实验方案和分析方法的优化，达到使用尽可能少的样本，快速发现可能与特定表型相关的基因。

作为进一步地改进，本发明所述的四个指标是：方法成功率、假阳性发现数、候选基因显著性和非孟德尔表型显著性，特定表型是孟德尔性状表型或具有明显主效基因的数量性状表型。

作为进一步地改进，本发明所述的方法成功率是真正与所研究表型相关的基因能够被发现的概率，假阳性发现数是在发现的全部基因中，假阳性发现的基因个数的期望值，候选基因显著性是每一个发现的基因与表型相关的统计显著性，非孟德尔表型显著性是当无法发现任何基因或所有已发现的基因都经由实验验证与表型无关时，纳入实验研究的个体实际上不是具有同一特定表型的个体，这一可能性的显著性。

作为进一步地改进，本发明所述的四个指标的概率计算框架式为：

$H\left(M,N\right)=P\left(\left|\mathrm{\Φ}\right|\≥M\right)=\underset{\∀\mathrm{\Φ}:\left|\mathrm{\Φ}\right|\≥M}{\Σ}P\left(\mathrm{\Φ}\right)=\underset{\∀\mathrm{\Φ}:\left|\mathrm{\Φ}\right|\≥M}{\Σ}\left(\underset{S_i\∈\mathrm{\Φ}}{\Π}{D}_{S_i}\×\underset{S_i\∉\mathrm{\Φ}}{\Π}\left(1-D_{S_i}\right)\right)$

$\begin{array}{c}{Q}_{G_j}=P\left(\left|{\mathrm{\Φ}}_{G_j}\right|\≥M\right)=\underset{\∀{\mathrm{\Φ}}_{G_j}:\left|{\mathrm{\Φ}}_{G_j}\right|\≥M}{\Σ}P\left({\mathrm{\Φ}}_{G_j}\right)\\ =\underset{\∀{\mathrm{\Φ}}_{G_j}:\left|{\mathrm{\Φ}}_{G_j}\right|\≥M}{\Σ}\left(\underset{S_i\∈{\mathrm{\Φ}}_{G_j}}{\Π}\left(1-{\left(1-B_{S_i}\right)}^{L_{G_j}}\right)\×\underset{S_i\∉{\mathrm{\Φ}}_{G_j}}{\Π}{\left(1-B_{S_i}\right)}^{L_{G_j}}\right)\end{array}$

$E(M,N)=\underset{\∀G_j}{\Σ}{Q}_{G_j}$

$Z(w^{\′},N)\≤Z(w,N)=P\left(\right|\mathrm{\Φ}|\≤w)={\Σ}_{\∀\mathrm{\Φ}:\left|\mathrm{\Φ}\right|\≤w}P\left(\mathrm{\Φ}\right)={\Σ}_{\∀\mathrm{\Φ}:\left|\mathrm{\Φ}\right|\≤w}\left({\Π}_{S_i\∈\mathrm{\Φ}}{D}_{S_i}\×{\Π}_{S_i\∉\mathrm{\Φ}}\left(1-D_{S_i}\right)\right)$

作为进一步地改进，本发明所述的四个指标基于递推思想的计算式为：

${\mathrm{PV}}_{G_j,S_i}=1-{(1-B_{S_i})}^{L_{G_j}}$

$P(A_u=v)=\left\{\begin{array}{cc}{D}_{S_u}\&times;P(A_{u-1}=v-1)+(1-D_{S_u})\&times;P(A_{u-1}=v)& (0<v<u)\\ {\&Pi;}_{i=1}^u(1-D_{S_u})& (v=0)\\ {\&Pi;}_{i=1}^u{D}_{S_u}& (v=u)\end{array}\right.$

$P(R_{G_j,u}=v)=\left\{\begin{array}{cc}{\mathrm{PV}}_{G_j,S_u}\&times;P(R_{G_j,u-1}=v-1)+(1-{{\mathrm{PV}}_{G_j,}}_{S_u})\&times;P(R_{G_j,u-1}=v)& (0<v<u)\\ {\&Pi;}_{i=1}^u(1-{\mathrm{PV}}_{G_j,S_u})& (v=0)\\ {\&Pi;}_{i=1}^u{\mathrm{PV}}_{G_j,S_u}& (v=u)\end{array}\right.$

$H(M,N)=P(A_N\&GreaterEqual;M)=\underset{k=M}{\overset{N}{\&Sigma;}}P(A_N=k)$

$Q_{G_j}=P(R_{G_j,N}\&GreaterEqual;M)=\underset{k=M}{\overset{N}{\&Sigma;}}P(R_{G_j,N}=k)$

$E(M,N)=\underset{\&ForAll;G_j}{\&Sigma;}{Q}_{G_j}$

$Z(w^{\&prime;},N)\&le;Z(w,N)=P(A_N\&le;w)={\&Sigma;}_{k=0}^{w}P(A_N=k).$

作为进一步地改进，本发明所述的计算式中，

S_i：一个样本；

G_j：一个基因；

基因组中基因G_j的区域长度；

一个分析方法在样本S_i中能够检测到与表型相关的那个或那些基因组变异的概率，表征单个样本的测序结果中观察到与所研究表型相关基因组变异的敏感度；

一个分析方法在经过基因组变异筛选的步骤以后在样本S_i上发现的基因组变异的频率，每个碱基有多大概率对应一个基因组变异，表征单个样本的测序结果中观察到与所研究表型相关基因组变异的特异性；

N：纳入研究的相同表型独立样本的总数；

M：研究人员指定的候选基因报告标准；候选基因必须在不少于M个相同表型独立样本中至少存在一个发现的基因组变异；

w′：当一个基于基因组重测序技术的正向遗传学研究过程无法发现可能与表型相关的基因或已发现的基因都已经被验证为假阳性发现时，真正与表型相关的基因在w′个相同表型独立样本中至少存在一个发现的基因组变异，w′<M；

w：当一个基于基因组重测序技术的正向遗传学研究过程无法发现可能与表型相关的基因或已发现的基因都已经被验证为假阳性发现时，不满足候选基因报告条件，未被证明与表型无关的基因，最多在w个相同表型独立样本中至少存在一个发现的基因组变异，w′≤w<M；

C_N：从N个样本中取出任意个样本的组合；

Φ：一种在样本中发现基因组变异的可能性。Φ代表一种样本组合：属于该组合的样本中，真正与表型相关的基因组变异被发现了，不属于该组合的样本中，真正与表型相关的基因组变异未被发现，Φ∈C_N；

一种在样本中发现基因组变异的可能性，代表一种样本组合：属于该组合的样本中，在G_j区域都发现了至少一个的基因组变异，不属于该组合的样本中，在G_j区域均没有发现的基因组变异，

||：计算组合(Φ或)中样本的个数；

H(M,N)：方法成功率指标；

候选基因显著性指标；

E(M,N)：假阳性发现数指标；

Z(w′,N)：非孟德尔表型显著性指标；

u，v：递推变量，可在递推公式中作为样本号替代S_i的下标i；

随机情况下，样本S_i中基因G_j至少包含一个发现的基因组变异的概率；

A_u：N个样本中前u个样本中，真正与表型相关的基因组变异被发现的样本数；

N个样本中前u个样本中，随机情况下，基因G_j至少包含一个(经过基因组变异筛选的步骤以后)发现的基因组变异的样本数。

作为进一步地改进，本发明所述的和通过以下计算式计算：

$D_{S_i}=H_{S_i}\&times;(1-{\mathrm{\&Phi;}}_{S_i});$

$B_{S_i}={\mathrm{Nv}}_{S_i}/l;$

其中：

样本S_i的基因组变异检测灵敏度VCS；

样本S_i的筛选器误删率FIR；

样本S_i经过筛选器筛选以后，剩下的可能与表型相关的基因组变异个数；

l：实验研究所关注的基因组区域。

作为进一步地改进，本发明所述的通过以下计算式计算：

${\mathrm{\&Phi;}}_{S_i}=1-(1-{\mathrm{\&Phi;}}_{S_i,F})\&times;(1-{\mathrm{\&Phi;}}_{S_i,S}),$ 其中：

第一类筛选器的整体误删率；

第二类筛选器的整体误删率。

第一类筛选器的整体误删率通过用同样的一组第一类筛选器筛选所研究物种的所有已知基因组变异，评估多少已知能够产生非特定孟德尔表型的变异被筛除，来估计；

第二类筛选器的整体误删率通过计算所有二类筛选方法所认定的不可能基因组区域占全部实验研究所关注的基因组区域的比例，来估计。

作为进一步地改进，本发明筛选器误删率计算方法中，有效区筛选和显著差异筛选属于第一类筛选器；稠密度筛选、同祖先筛选、对照组筛选属于第二类筛选器。

作为进一步地改进，本发明所述的通过计算四个指标评估研究过程中每一个步骤的设计对研究整体有效性造成的影响，从而指导对当前实验方案和分析方法的优化的具体方法是：

1)、当一个基于基因组重测序技术的正向遗传学研究过程只发现了少数几个候选基因，且候选基因显著性较强，在经过严格的多重统计矫正后仍具有显著性时，不论方法成功率高低，实验验证发现的候选基因；

2)，当一个基于基因组重测序技术的正向遗传学研究过程找到的与表型相关的基因很多，超过了可以直接进行实验验证的合理数量时，参考假阳性发现数，假阳性发现数较低，则从候选基因显著性最强的候选基因开始，逐个进行实验验证，假阳性发现数较高，则说明所采用的数据分析方法标准过松，分析中发现了大量假阳性基因组变异，从而导致发现了很多的假阳性候选基因，对数据分析流程包括完成各个分析步骤的软件选择和参数设定，以及最终的候选基因报告标准进行优化，使用更为严格的基因组变异分析方法、基因组变异筛选方法和/或更严格的候选基因报告标准，对分析方法的不同优化有不同的有效性，最有效的优化方法应该是在不显著降低方法成功率的条件下，极大地降低假阳性发现数，研究人员可以尝试多种优化方法，评估优化后的分析方法的四个指标，然后选出最有效的优化方法；

3)、当一个基于基因组重测序技术的正向遗传学研究过程无法发现可能与表型相关的基因或已发现的基因都已经被验证为假阳性发现时，参考非孟德尔表型显著性，非孟德尔表型显著性不强，可能有两个原因：一是数据分析所采用的方法标准过严，导致与表型相关的突变未被发现，从而使表型相关基因未能在足够多的样本中被发现含有突变；二是样本量不够，无法有效发现表型相关基因；首先对数据分析流程包括完成各个分析步骤的软件选择和参数设定，以及最终的候选基因报告标准进行优化，使用更为宽松的基因组变异分析方法、基因组变异筛选方法和/或更宽松的候选基因报告标准，最有效的优化方法应该是在不显著增加假阳性发现数的条件下，极大地增加方法成功率；

4)、方法优化后可以找到符合情况一或情况二的候选基因，则进行实验验证；

5)、不论怎样优化方法，都无法发现符合情况一或情况二的候选基因，且非孟德尔表型显著性仍然不强，则考虑是由于纳入研究的相同表型独立样本数过少，建议增加样本量；

6)、经过了方法优化和/或进一步研究，仍无法发现与表型相关的基因或已发现的基因都已经被验证为假阳性发现，且非孟德尔表型显著性较强，首先检查纳入研究的样本的表型是否存在观察错误，如果没有错误，则测量更多的表型特征，以确认和/或细分纳入研究的样本的表型，然后依据更严格的表型标准来进行研究，由于基因通常具有多效性，两个功能类似的基因的变异可能产生类似的主要表型，但在次要表型上总会存在不同，利用更严格的表型标准，获得同一基因突变的表型个体，从而使得基于基因组重测序技术的正向遗传学研究过程可以正确发现这一表型相关基因。

本发明的有益效果在于：

1)、本专利首次对整个基于基因组重测序技术的正向遗传学研究过程进行数学建模，利用精确的数学方法评估研究过程的有效性，并对研究过程进行最优化以提升研究效率；

2)、本专利首次提出用四个指标来评估基于基因组重测序技术的正向遗传学研究过程的整体有效性，其中方法成功率和非孟德尔表型显著性是创新提出的指标，对于指导整个研究过程优化有重要价值，假阳性发现数与候选基因显著性在此前的技术中有类似指标的提及，但它们没有一起，和/或与方法成功率和非孟德尔表型显著性一起被用于评估整个研究过程的有效性。

3)、本专利首次设计了基于四个指标优化基因组重测序技术的正向遗传学研究过程的方法；

4)、本专利首次提出了，基于单个样本的测序结果中观察到与所研究表型相关基因组变异的敏感度和特异性，来推导整个研究过程有效性的四个指标的概率计算框架；

5)、本专利首次提出了利用递推公式，基于单个样本的测序结果中观察到与所研究表型相关基因组变异的敏感度和特异性，来推导整个研究过程有效性的四个指标的概率计算框架，极大地简化了计算；

6)、本专利首次提出了在一个实测个体的基因组中引入假想的基因组变异，并根据实测数据的数据质量，变更实测数据以模拟个体基因组引入假想变异后的测序结果，并依据多大比例的假想变异能够正确地被特定的分析过程发现，来评估这一特定的分析过程依据这一组测序数据发现对应个体基因组中存在变异的敏感度的方法；

7)、本专利首次提出将筛除与所研究表型无关的基因组变异的筛选方法分为两类。第一类是依据基因组变异的预测功能筛除变异的筛选方法，第二类是依据测序质量或对照样本认定特定基因组区域不存在与所研究表型相关变异的筛选方法。所有第一类筛选方法对检出表型相关变异的敏感度的整体影响，可以通过用同样的一组方法，筛选一个已知能够产生表型的变异库，评估多大比例的已知能产生表型的变异会被筛除来估计，所有第二类筛选方法对检出表型相关变异的敏感度的整体影响，可以通过评估所有这些方法认定的不可能基因组区域占全部可能检出区域的比例来估计，这一基因组变异筛选方法的分类以及每类方法对检出表型相关变异的敏感度的整体影响的计算方法是创新的；

8)、基于第7点，本专利首次提出外显子组测序等针对特定区域的测序技术，以及在分析中仅关注特定基因组区域的分析方法，本质上是认定其它区域的基因组变异不具备产生表型的功能，因此是一种第一类筛选方法，其对检出表型相关变异的敏感度的影响可与其它第一类筛选方法一起整体估计；

9)、基于第7点，本专利首次提出依据变异注释软件或变异对功能影响的预测软件来筛除功能较弱的变异的方法，包括但不限于SNPEff、ANNOVAR、依据氨基酸残基相似性矩阵筛除相似的非同义编码突变等，均属于第一类筛选方法，其对检出表型相关变异的敏感度的影响可与其它第一类筛选方法一起整体估计；

10)、基于第7点，本专利首次提出依据特定表型的多个独立个体产生表型的基因组变异应独立不相同这一假设，认为多个独立个体中均出现的基因组变异为研究样本的背景基因组变异，在寻找表型相关基因时予以筛除的方法，属于第二类筛选方法，其对检出表型相关变异的敏感度的影响可与其它第二类筛选方法一起整体估计；

11)、基于第7点，本专利首次提出在含有家系或其它对照的分析过程中，依据对照样本出现基因组变异的区域筛除表型样本可能在对应位置出现的突变的方法，属于第二类筛选方法，其对检出表型相关变异的敏感度的影响可与其它第二类筛选方法一起整体估计；

12)、基于第7点，本专利首次提出依据均匀分布的基因组变异不太可能过分集中于较小的基因组区域，集中于较小基因组区域的多个变异往往是由于序列比对结果不精确造成的假阳性基因组变异这一假说，筛除一个特定较小基因组区域内出现的多个变异的方法，属于第二类筛选方法，其对检出表型相关变异的敏感度的影响可与其它第二类筛选方法一起整体估计。

13)、本发明的技术方案可应用于生物学、医学、药学、农学的相关研究和应用领域。包括但不限于，生物学基础研究，如发现表型调控的相关机制等，个性化医疗，如发现病人的致病相关基因并针对性地进行干预或设计长期疾病风险防范和/或生殖风险防范的措施等，精准治疗，如发现病人的致病相关基因或致病微生物的重要表型相关基因并针对性地选择药物或其它干预方法以达到最好的疗效和/或最低的风险等，分子育种，如发现农艺表型相关基因并进行相关遗传改造或转基因以提升种质等，精准农药和虫害防控如发现病害源生物的重要表型相关基因以利用对其杀灭或降低其致害或繁殖能力。

具体实施方式

本发明公开了一种基于概率框架和重测序技术快速发现表型相关基因的方法，对基于基因组重测序技术的正向遗传学研究过程建模，通过计算四个指标评估研究过程中每一个步骤的设计对研究整体有效性造成的影响，从而指导对当前实验方案和分析方法的优化，达到使用尽可能少的样本，快速发现可能与特定表型相关的基因。

特定表型指孟德尔性状表型或具有明显主效基因的数量性状表型。这类表型的发生与单一基因型，单个基因改变或几个基因同时改变相关，可应用本专利，没有明显主效基因的数量性状表型不由单一基因型控制，本专利不适用。

四个指标是指，在使用当前的实验设计和分析方法时，1)真正与所研究表型相关的基因能够被发现的概率，以下称为方法成功率；2)在发现的全部基因中，假阳性发现，与所研究表型无关的基因个数的期望值，以下称为假阳性发现数；3)每一个发现的基因与表型相关的统计显著性，以下称为候选基因显著性；4)当无法发现任何基因或所有已发现的基因都经由实验验证与表型无关时，纳入实验研究的个体实际上不是具有同一特定表型的个体，这一可能性的显著性，以下称为非孟德尔表型显著性。

基于基因组重测序技术的正向遗传学研究过程是指：

第一步，获得具有同一特定表型的多个独立个体。需要保证在这些待研究的表型个体中，基因组变异大致均匀分布在全基因组上。例如，使用化学诱变剂，如甲基磺酸乙酯等处理大批同一遗传背景的植物种子，获得的M2代植株群体中显示同一表型的多个植株；或者，患同一遗传病的多个独立家系，其中的患病者在排除了不患病的家系成员携带的与疾病无关的基因组变异后，剩下的可能与疾病相关的基因组变异也可大致认为均匀分布于患病者的基因组上。

第二步，对同一特定表型的多个独立个体进行基因组重测序。对于前述需要家系对照才能保证基因组变异均匀分布的情况，需要对家系对照个体也进行基因组重测序。

第三步，对测序结果进行分析，选择软件和参数，将测序结果比对到基因组上，分析获得每个个体的基因组变异，对这些变异进行筛选，获得可能与表型相关的变异。对于前述需要家系对照才能保证基因组变异均匀分布的情况，在显示表型的个体中发现的基因组变异，需要剔除在家系对照个体中发现的与表型无关的基因组变异后再进行后续分析。

第四步，选定一个报告标准，报告与表型相关的基因。例如，在N个具有特定表型的独立个体中，至少有M个个体的基因组在基因G_i的区域内包含至少一个经过筛选后仍可能与表型相关的变异，则报告基因Gi与表型相关。

通过计算四个指标评估研究过程中每一个步骤的设计对研究整体有效性造成的影响，从而指导对当前实验方案和分析方法的优化的具体方法是：

一，当一个基于基因组重测序技术的正向遗传学研究过程只发现了少数几个候选基因，且候选基因显著性较强，在经过严格的多重统计矫正后仍具有显著性时，不论方法成功率高低，都建议实验验证发现的候选基因。

二，当一个基于基因组重测序技术的正向遗传学研究过程找到的可能与表型相关的基因很多，超过了可以直接进行实验验证的合理数量时，可以参考假阳性发现数。如果假阳性发现数较低，则建议从候选基因显著性最强的候选基因开始，逐个进行实验验证。如假阳性发现数较高，则说明所采用的数据分析方法标准过松，分析中发现了大量假阳性基因组变异，从而导致发现了很多的假阳性候选基因。此时，建议对数据分析流程，包括完成各个分析步骤的软件选择和参数设定，以及最终的候选基因报告标准进行优化，使用更为严格的基因组变异分析方法、基因组变异筛选方法和/或更严格的候选基因报告标准。对分析方法的不同优化有不同的有效性。最有效的优化方法应该是在不显著降低方法成功率的条件下，极大地降低假阳性发现数。研究人员可以尝试多种优化方法，评估优化后的分析方法的四个指标，然后选出最有效的优化方法。

三，当一个基于基因组重测序技术的正向遗传学研究过程无法发现可能与表型相关的基因或已发现的基因都已经被验证为假阳性发现时，可参考非孟德尔表型显著性。如果非孟德尔表型显著性不强，可能有两个原因。一是数据分析所采用的方法标准过严，导致与表型相关的突变未被发现，从而使表型相关基因未能在足够多的样本中被发现含有突变。二是样本量不够，无法有效发现表型相关基因。因此，建议首先对数据分析流程，包括完成各个分析步骤的软件选择和参数设定，以及最终的候选基因报告标准进行优化，使用更为宽松的基因组变异分析方法、基因组变异筛选方法和/或更宽松的候选基因报告标准。最有效的优化方法应该是在不显著增加假阳性发现数的条件下，极大地增加方法成功率。

如果方法优化后可以找到符合情况一或情况二的候选基因，则进行实验验证。

如果不论怎样优化方法，都无法发现符合情况一或情况二的候选基因，且非孟德尔表型显著性仍然不强，则考虑是由于纳入研究的相同表型独立样本数过少，建议增加样本量。

如果经过了方法优化和/或进一步研究，仍无法发现可能与表型相关的基因或已发现的基因都已经被验证为假阳性发现，且非孟德尔表型显著性较强，例如P<0.05，则可能是由于纳入研究的样本并非严格意义上的相同表型独立样本。这时建议首先检查纳入研究的样本的表型是否存在观察错误。如果没有错误，则建议测量更多的表型特征，以确认和/或细分纳入研究的样本的表型。然后依据更严格的表型标准来进行研究。由于基因通常具有多效性，两个功能类似的基因的变异可能产生类似的主要表型，但在次要表型上总会存在不同。利用更严格的表型标准，可以获得同一基因突变的表型个体，从而使得基于基因组重测序技术的正向遗传学研究过程可以正确发现这一表型相关基因。

为了计算评价一个基于基因组重测序技术的正向遗传学研究过程的四个指标，本专利设计了一个概率计算框架，依据单个样本的测序结果中观察到与所研究表型相关基因组变异的敏感度和特异性，来推导整个研究过程有效性的四个指标，计算式为：

$H\left(M,N\right)=P\left(\left|\mathrm{\&Phi;}\right|\&GreaterEqual;M\right)=\underset{\&ForAll;\mathrm{\&Phi;}:\left|\mathrm{\&Phi;}\right|\&GreaterEqual;M}{\&Sigma;}P\left(\mathrm{\&Phi;}\right)=\underset{\&ForAll;\mathrm{\&Phi;}:\left|\mathrm{\&Phi;}\right|\&GreaterEqual;M}{\&Sigma;}\left(\underset{S_i\&Element;\mathrm{\&Phi;}}{\&Pi;}{D}_{S_i}\&times;\underset{S_i\&NotElement;\mathrm{\&Phi;}}{\&Pi;}\left(1-D_{S_i}\right)\right)$

$\begin{array}{c}{Q}_{G_j}=P\left(\left|{\mathrm{\&Phi;}}_{G_j}\right|\&GreaterEqual;M\right)=\underset{\&ForAll;{\mathrm{\&Phi;}}_{G_j}:\left|{\mathrm{\&Phi;}}_{G_j}\right|\&GreaterEqual;M}{\&Sigma;}P\left({\mathrm{\&Phi;}}_{G_j}\right)\\ =\underset{\&ForAll;{\mathrm{\&Phi;}}_{G_j}:\left|{\mathrm{\&Phi;}}_{G_j}\right|\&GreaterEqual;M}{\&Sigma;}\left(\underset{S_i\&Element;{\mathrm{\&Phi;}}_{G_j}}{\&Pi;}\left(1-{\left(1-B_{S_i}\right)}^{L_{G_j}}\right)\&times;\underset{S_i\&NotElement;{\mathrm{\&Phi;}}_{G_j}}{\&Pi;}{\left(1-B_{S_i}\right)}^{L_{G_j}}\right)\end{array}$

$E(M,N)=\underset{\&ForAll;G_j}{\&Sigma;}{Q}_{G_j}$

$Z(w^{\&prime;},N)\&le;Z(w,N)=P\left(\right|\mathrm{\&Phi;}|\&le;w)=\underset{\&ForAll;\mathrm{\&Phi;}:\left|\mathrm{\&Phi;}\right|\&le;w}{\&Sigma;}P\left(\mathrm{\&Phi;}\right)=\underset{\&ForAll;\mathrm{\&Phi;}:\left|\mathrm{\&Phi;}\right|\&le;w}{\&Sigma;}\left(\underset{S_i\&Element;\mathrm{\&Phi;}}{\&Pi;}{D}_{S_i}\&times;\underset{S_i\&NotElement;\mathrm{\&Phi;}}{\&Pi;}\left(1-D_{S_i}\right)\right)$

其中：

S_i：一个样本。

G_j：一个基因。

基因组中基因G_j的区域长度bp。

一个分析方法在样本S_i中能够检测到与表型相关的那个或那些基因组变异的概率。英文描述为Variantdetectionsensitivity，以下文章中简记为VDS，表征单个样本的测序结果中观察到与所研究表型相关基因组变异的敏感度。

一个分析方法在经过基因组变异筛选的步骤以后在样本S_i上发现的基因组变异的频率每个碱基有多大概率对应一个基因组变异。英文描述为Backgroundvariantfrequency，以下文章中简记为BVF，表征单个样本的测序结果中观察到与所研究表型相关基因组变异的特异性。

N：纳入研究的相同表型独立样本的总数。

M：研究人员指定的候选基因报告标准。候选基因必须在不少于M个相同表型独立样本中至少存在一个经过基因组变异筛选后发现的基因组变异。

w′：当一个基于基因组重测序技术的正向遗传学研究过程无法发现可能与表型相关的基因或已发现的基因都已经被验证为假阳性发现时，真正与表型相关的基因在w′个相同表型独立样本中至少存在一个经过基因组变异筛选后发现的基因组变异，w′<M。

w：当一个基于基因组重测序技术的正向遗传学研究过程无法发现可能与表型相关的基因或已发现的基因都已经被验证为假阳性发现时，不满足候选基因报告条件未被证明与表型无关的基因，最多在w个相同表型独立样本中至少存在一个经过基因组变异筛选后发现的基因组变异，w′≤w<M。

C_N：从N个样本中取出任意个样本的组合。

Φ：一种在样本中发现基因组变异的可能性。Φ代表一种样本组合：属于该组合的样本中，真正与表型相关的基因组变异，经过基因组变异筛选的步骤以后被发现了；不属于该组合的样本中，真正与表型相关的基因组变异经过基因组变异筛选的步骤以后未被发现，Φ∈C_N。

一种在样本中发现基因组变异的可能性。代表一种样本组合：属于该组合的样本中，在G_j区域都发现了至少一个，经过基因组变异筛选的步骤以后的基因组变异；不属于该组合的样本中，在G_j区域均没有发现，经过基因组变异筛选的步骤以后的基因组变异。

||：计算组合(Φ或)中样本的个数。

评价一个基于基因组重测序技术的正向遗传学研究过程的四个指标为：

H(M,N)：方法成功率指标。

候选基因Gj与表型关联的显著性，由与表型无关的基因组变异随机聚集于Gj基因区域使得Gj被作为候选基因发现的概率。

E(M,N)：假阳性发现数，存在由与表型无关的基因组变异随机聚集于某个基因区域使得该基因被作为候选基因发现，E(M,N)为全基因组内上述假阳性发现个数的期望值。

Z(w′,N)：非孟德尔表型显著性，如果在大量(N)相同表型独立样本中只有非常少(w′)的样本中同一个基因发生变异，且基因组变异的样本检测灵敏度不低，则应考虑纳入研究的样本并非严格意义上的相同表型独立样本的可能性，w′无法直接获得，但w可以直接观察，易得：Z(w′,N)≤Z(w,N)。Z(w,N)是对Z(w′,N)的保守估计。

以上计算方法实现了本专利的设计思想，但需要枚举从N个样本中取出任意个样本的组合情况，计算在每一种情况下的各种概率，为了简化计算，我们设计了一种基于递推思想的算法，可以迅速获得上述组合算法的精确解，极大地简化了计算：

${\mathrm{PV}}_{G_j,S_i}=1-{(1-B_{S_i})}^{L_{G_j}}$

$P(A_u=v)=\left\{\begin{array}{cc}{D}_{S_u}\&times;P(A_{u-1}=v-1)+(1-D_{S_u})\&times;P(A_{u-1}=v)& (0<v<u)\\ {\&Pi;}_{i=1}^u(1-D_{S_u})& (v=0)\\ {\&Pi;}_{i=1}^u{D}_{S_u}& (v=u)\end{array}\right.$

$P(R_{G_j,u}=v)=\left\{\begin{array}{cc}{\mathrm{PV}}_{G_j,S_u}\&times;P(R_{G_j,u-1}=v-1)+(1-{{\mathrm{PV}}_{G_j,}}_{S_u})\&times;P(R_{G_j,u-1}=v)& (0<v<u)\\ {\&Pi;}_{i=1}^u(1-{\mathrm{PV}}_{G_j,S_u})& (v=0)\\ {\&Pi;}_{i=1}^u{\mathrm{PV}}_{G_j,S_u}& (v=u)\end{array}\right.$

$H(M,N)=P(A_N\&GreaterEqual;M)=\underset{k=M}{\overset{N}{\&Sigma;}}P(A_N=k)$

$Q_{G_j}=P(R_{G_j,N}\&GreaterEqual;M)=\underset{k=M}{\overset{N}{\&Sigma;}}P(R_{G_j,N}=k)$

$E(M,N)=\underset{\&ForAll;G_j}{\&Sigma;}{Q}_{G_j}$

$Z(w^{\&prime;},N)\&le;Z(w,N)=P(A_N\&le;w)=\underset{k=0}{\overset{w}{\&Sigma;}}P(A_N=k)$

补充定义：

u,v：递推变量，可在递推公式中作为样本号替代S_i的下标i。

随机情况下，样本S_i中基因G_j至少包含一个，经过基因组变异筛选的步骤以后发现的基因组变异的概率。

A_u：N个样本中前u个样本中，真正与表型相关的基因组变异，经过基因组变异筛选的步骤以后被发现的样本数。

N个样本中前u个样本中，随机情况下，基因G_j至少包含一个经过基因组变异筛选的步骤以后发现的基因组变异的样本数。

上述两个计算公式中，仅有(单个样本的测序结果中观察到与所研究表型相关基因组变异的敏感度，VDSVariantcallingsensitivity)和(基因组变异的频率，BVFbackgroundvariantfrequency)需要通过实际测序数据估计，本专利设计了以下的估计方法：

$D_{S_i}=H_{S_i}\&times;(1-{\mathrm{\&Phi;}}_{S_i})$

$B_{S_i}={\mathrm{Nv}}_{S_i}/l$

样本S_i的基因组变异检测灵敏度VCS；

样本S_i的筛选器误删率FIR；

样本S_i经过筛选器筛选以后，剩下的可能与表型相关的基因组变异个数；

l：实验研究所关注的基因组区域，如外显子区域或全基因组的总长度，碱基数。

上述计算公式中，样本S_i的基因组变异检测灵敏度样本S_i的筛选器误删率需要估计，在基于基因组重测序技术的正向遗传学研究过程中，每个样本的和可以通过以下方式计算。可以由样本S_i的基因组变异检测灵敏度(VariantCallingSensitivity,VCS)，以及基因组变异检出后，各种变异筛选器错误滤除与表型相关变异的误删率(FalseIgnoranceRate,FIR)，两个指标计算。样本S_i经过筛选器筛选以后，剩下的可能与表型相关的基因组变异，相对于基因组区域的长度发生的频率，即为

其中，样本S_i的基因组变异检测灵敏度与样本S_i的测序质量、测序深度、分析其数据时使用的比对软件(bowtie2、BWA等)、分析其数据时所使用的基因组变异发现软件(GATK、samtools等)，以及这些软件的使用参数有关。

为估计本专利设计了一种基于模拟的方法：在样本S_i的基因组上引入假想的一组，例如10000个基因组变异，根据这个假想的基因组模拟一组和实测数据相同质量的测序数据，对于实测数据中每一条覆盖了假想基因组变异的序列，依据对应于假想变异位置的实测碱基的测序质量，模拟这个位置的碱基以一定的概率正确/错误地测定了这个假想的基因组变异。最后得到的数据集即是与实测数据质量与深度分布一致的，对假想基因组的测序模拟数据。然后，使用同样的基因组变异发现方法，使用同样的序列比对软件、变异发现软件或方法，选择同样的参数，分析假想基因组的测序模拟数据。所得结果能够正确发现之前引入的假想基因组变异的比例，即为一个对使用特定的分析方法分析样本S_i的一组实测数据能够发现真实的基因组变异的检测灵敏度

样本S_i的筛选器误删率仅与筛选器所采用的策略有关，与样本S_i的实测数据无关，但由于对不同的样本可能采用不同的筛选器，因此筛选器误删率仍与样本相关，记为为估计本专利设计了以下的方法：首先，将基因组变异的筛选策略分为两类，第一类通过预测基因组变异的功能，筛除不能产生任何表型的变异；第二类依据测序质量或对照样本认定特定基因组区域内不存在与所研究表型相关变异，可能存在产生其它表型的变异。

${\mathrm{\&Phi;}}_{S_i}=1-(1-{\mathrm{\&Phi;}}_{S_i,F})\&times;(1-{\mathrm{\&Phi;}}_{S_i,S})$

第一类筛选器的整体误删率。

第二类筛选器的整体误删率。

在上述公式中，第一类筛选器的整体误删率可以通过用同样的一组第一类筛选器筛选所研究物种的所有已知基因组变异，评估多少已知能够产生，非特定孟德尔表型的变异被筛除，来估计。这一方法要求所研究物种具有大量已知的基因型－表型关联关系。对于没有这样信息的物种，比如多数高等植物，我们也可以通过筛选人类的基因组库变异来估计误删率，在一定的精度下，高等生物相同基因组功能区域变异产生表型的概率可认为大致相等。ClinVAR数据库全面地收集了已知人类基因组变异与表型的关联关系。

第二类筛选器的整体误删率可以通过计算所有二类筛选方法所认定的不可能基因组区域占全部实验研究所关注的基因组区域的比例，来估计。这一方法的原理是估计当二类筛选方法的证据出现错误时，可能误删真正与表型相关变异的概率。

在之前的基因组重测序应用中，有一批常用的筛选与表型相关基因组变异的方法，主要分为五类，分别是，有效区筛选(Effectiveregionfilter)，显著差异筛选(Bigdifferencefilter)，稠密度筛选(Congestionfilter)，同祖先筛选(Ancestryfilter)，对照组筛选(Controlfilter)，以下分别说明：

在前述筛选器误删率计算方法中，有效区筛选和显著差异筛选属于第一类筛选器。

有效区筛选(Effectiveregionfilter)筛除不在用户所指定的“有效区”范围内的所有突变。有效区通常通过基因组功能区域指定，为一组功能区域的组合，包括，启动子区(promoting)、剪切位点区(splicing)、信使RNA的非翻译区(5-UTR、3-UTR)以及编码蛋白区(CDS)中的几个或全部。有效区筛选本质上是认定其它区域的基因组变异不具备产生表型的功能，因此是一种第一类筛选方法，其对检出表型相关变异的“敏感度”的影响可与其它第一类筛选方法一起整体估计。

显著性差异筛选(Bigdifferencefilter)筛除预测对于蛋白质功能(包括其活性调控)影响较小的变异，例如，筛除SNPEff、ANNOVAR等软件工具预测的低功能影响突变，和/或依据氨基酸残基相似性矩阵筛除相似的非同义编码突变等。这类筛选方法属于第一类筛选方法，其对检出表型相关变异的“敏感度”的影响可与其它第一类筛选方法一起整体估计。

稠密度筛选、同祖先筛选、对照组筛选(Controlfilter)属于第二类筛选器。

稠密度筛选(Congestionfilter)筛除集中于较小基因组区域的多个变异。随机产生的基因组变异不太可能集中分布于特定的基因组区域。集中于较小基因组区域(例如11个临近碱基内)的多个变异往往是由于序列比对结果不精确造成的假阳性基因组变异。基因组上存在序列高度相似的区域，给序列比对造成了困难。当比对位置出错时，比对错误的区域常发现假阳性的基因组变异。这一方法的本质是依据测序质量认定特定基因组区域内不存在与所研究表型相关变异，因此属于第二类筛选方法，其对检出表型相关变异的“敏感度”的影响可与其它第二类筛选方法一起整体估计。

同祖先筛选(Ancestryfilter)筛除在两个以上独立个体中发现的同一位点的突变。如果研究的样本间不存在亲缘关系，则两个独立样本中基因随机突变位于同一位点的概率极低。多个独立个体中均出现的基因组变异实质为研究样本的整体背景基因组变异，与所研究表型无关。这一方法的本质是依据样本对照认定特定基因组区域内不存在与所研究表型相关变异，因此属于第二类筛选方法，其对检出表型相关变异的“敏感度”的影响可与其它第二类筛选方法一起整体估计。

对照组筛选(Controlfilter)依据不表现表型的个体基因组，筛除表现表型的个体基因组中与表型无关的基因组变异。此外，在以家系为单位的研究中，家系内对照可以用于筛除与所研究表型无关的，家系相对于参考基因组的背景变异，使得留存的变异可被认为大致均匀分布(随机发生)于全基因组内。这一方法属于第二类筛选方法，其对检出表型相关变异的“敏感度”的影响可与其它第二类筛选方法一起整体估计。

以上描述的方法和公式可以完整地基于实际研究数据和公开数据计算一个基于基因组重测序技术的正向遗传学研究过程的四个研究有效性指标。在实际应用中，所有依据实际实验数据估计的指标都可以人工指定，用以模拟估算指标改变后对于整体分析流程有效性的影响，从而指导对于分析流程和实验设计的不断改进和优化。

为说明本专利的应用，通过具体实施例子来说明，这个例子中，我们通过EMS诱变具有纯合Pho2突变的水稻种子，筛选获得Pho2突变表型，正常磷供给条件下产生磷毒害表型发生抑制的M2代突变体，通过测序三个这样的突变体快速定位到了抑制Pho2突变表型的基因。

步骤1：获取样本数据。

材料：Tos17插入pho2基因的纯合M1代水稻种子，表型为在正常磷供给下表现磷中毒表型。

操作：EMS试剂诱导M1代水稻种子突变，筛选磷中毒表型减弱的M2代株系。通过筛选约13000个株系得到三个株系，记为：M28、M29、M249。对三个株系分别进行外显子组测序。测序结果可在SRA数据库中获取。

水稻基因组参考序列可以通过以下链接下载。

ftp://ftp.plantbiology.msu.edu/pub/data/Eukaryotic_Projects/o_sativa/annotation_dbs/pseudomolecules/version_7.0/all.dir/

步骤2：使用常规分析方法尝试寻找表型相关基因。

(1)使用bowtie2软件，对参考序列建立索引(buildindex)；

(2)使用bowtie2软件，对将测序原始文件(2个fastq文件)与参考基因组进行比对，生成bam文件；

(3)使用picard软件，对bam文件进行排序(sortandaddreadgroup)、删除重复序列(removeduplicates)等操作，并用samtools软件对去重复后的bam文件生成索引；

(4)使用samtools软件以默认参数(比对质量13以上，碱基质量0以上)，查找突变位点(callSNP)，生成相关的突变位点信息文件(VCF文件)；

(5)执行并评估表型相关基因发现流程。将3个水稻植株样本的VCF文件作为输入，根据VCF文件中突变位点所在的染色体、序列位置等信息，使用有效区筛选(Effectiveregionfilter，设定为筛除启动子区(promoting)、剪切位点区(splicing)和编码蛋白区(CDS)以外的突变)，显著差异筛选(Bigdifferencefilter，设定为筛除同义突变)，同祖先筛选(Ancestryfilter)3种筛选器对可能与表型无关的突变位点进行筛选。突变筛选后，若一个基因在全部三个样本中仍包含筛选剩下的(可能与表型相关的)突变，则认为该基因为一个可能与表型相关的候选基因。

使用这一方法一共发现了28个候选基因。这些基因及其在每个样本中发生的变异位点详见附录数据。

步骤3：用上文中所述的概率框架计算四个(类)指标，评估分析流程并进行优化。

上述常规分析方法发现了28个候选基因。这个数量规模的候选基因难以通过实验方法逐一验证与表型的关联。通过计算四类指标值，发现：一，这些候选基因与表型的关联显著性不高。最显著的基因与表型关联的显著性仅有1.9E-6，无法通过Bonferroni多重统计矫正。二，方法成功率为74.5％，有使用更严格的基因组变异筛选方法(将降低方法成功率但提高所发现的候选基因的准确性)进行候选基因发现的空间。三，每个样本的变异发现敏感度VCS都较高(>80％)，有使用更严格的基因组变异发现方法(将降低变异发现敏感度但提升所发现变异的可靠性)的空间。

基于以上分析，我们优化了常规分析方法，设计了针对这个数据集的专用方法。这一专用方法采用了更严格的基因组变异发现方法与基因组变异筛选方法。具体来说，我们采用相同的软件发现基因组变异，但要求碱基质量的分值从13调整到15，要求匹配质量的分值从0调整到20。在筛选基因组变异时，将有效区设定为剪切位点区和编码蛋白区(启动子区突变不再考虑)，显著差异筛选设定为筛除同义突变和PAM120打分矩阵中蛋白质残基相似性分值为正值的非同义突变(原先仅筛除同义突变)。使用这一数据集专用方法，仅发现了一个与表型有关的候选基因LOC_Os02g56510(OsPHO1；2)。其在三个样本中发现的具体突变位点如下表所示。

该候选基因与表型关联的显著性也从1.9E-6提升为9.5E-7，说明从概率角度，由随机分布于基因组上的与表型无关的变异随机聚合于某个基因而导致该基因被作为候选基因发现的概率显著降低。这一降低的原因是我们采用了更严格的基因组变异发现和筛选方法，使得最终在每个样本中发现的可能与表型相关的基因组变异数量显著减少。

当然，也并不是采用越严格的基因组变异发现和筛选方法越好。严格的方法可能导致真正与表型相关的基因组变异无法被发现。当方法无法发现可能与表型相关的基因或已发现的基因都已经被验证为假阳性发现时，或者当发现了过多的候选基因无法逐一进行验证希望挑选最显著的进行验证时，对分析流程的优化策略在本专利文本的前部已有详细描述(“通过计算四个指标评估研究过程中每一个步骤的设计对研究整体有效性造成的影响，从而指导对当前实验方案和分析方法的优化策略”部分)

步骤4：实验验证基因与表型的关联。

通过以下四个方面的证据，我们验证了通过数据集的专用方法发现的这个候选基因与

表型的关联关系。

(1)通过文献查阅，发现该基因是拟南芥中AtPHO1基因的同源基因。AtPHO1基因在拟南芥中负责磷的运输。研究显示，AtPHO1是AtPHO2基因的一个重要的下游基因。AtPHO2可以调控AtPHO1在内膜上的降解，以维持磷在拟南芥中的平衡。由拟南芥中同源基因的功能，推定该候选基因可能存在“抑制Pho2突变表型”这一功能。

(2)小规模重复了一次正向遗传学实验。通过EMS试剂重新诱导M1代水稻种子突变，筛选磷中毒表型减弱的M2代株系。通过筛选约5000个株系得到一个新的株系，记为：M358。直接针对该株系的该候选基因测序，结果与预期一致，发现了该候选基因的一个对功能影响较强的非同义突变。

(3)对所有经由高通量测序发现的该候选基因突变进行了传统的PCR－Sanger测序验证。结果确认这几个突变均为真实存在的突变。

(4)对所有突变体的表型进行了验证。通过PCR－Sanger测序确认了M28、M29、M249、M358株系的遗传背景(含有Tos17插入pho2基因的纯合突变)。通过测定植物根部磷含量确认了pho2基因突变导致植物相对于正常植物在根部积累磷元素，导致磷中毒；M28、M29、M249、M358突变株在pho2基因突变背景下，降低根部磷元素积累，抑制pho2基因突变表型。M28、M29、M249、M358降低根部磷元素积累的水平大致相同。

经过实验验证我们可以得到以下结果，相比未发生任何突变的野生型，未经EMS诱导的原始pho2基因存在缺陷的纯合突变株长势枯黄矮小。而经过EMS进行诱导后的M28、M29、M249、M358这几株是由pho2基因存在缺陷的纯合突变株经EMS试剂再次诱变筛选得到的pho2表型抑制株。经过EMS试剂再次诱变后，原始植株对磷浓度的耐受力显著增强，其株高虽不及没有任何突变的野生植株，但明显优于未经EMS诱导的原始pho2基因存在缺陷的纯合突变株，这说明pho2基因确实与植物对磷浓度的耐受力有关。

以上所述的仅是本发明的优选实施方式，应当指出，对于本技术领域中的普通技术人员来说，在不脱离本发明核心技术特征的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (10)

1.一种基于概率框架和重测序技术快速发现表型相关基因的方法，其特征在于，对基于基因组重测序技术的正向遗传学研究过程建模，通过计算四个指标评估研究过程中每一个步骤的设计对研究整体有效性造成的影响，从而指导对当前实验方案和分析方法的优化，达到使用尽可能少的样本，快速发现可能与特定表型相关的基因。

2.根据权利要求1所述的基于概率框架和重测序技术快速发现表型相关基因的方法，其特征在于，所述的四个指标是：方法成功率、假阳性发现数、候选基因显著性和非孟德尔表型显著性，所述的特定表型是孟德尔性状表型或具有明显主效基因的数量性状表型。

3.根据权利要求2所述的基于概率框架和重测序技术快速发现表型相关基因的方法，其特征在于，所述的方法成功率是真正与所研究表型相关的基因能够被发现的概率，所述的假阳性发现数是在发现的全部基因中，假阳性发现的基因个数的期望值，所述的候选基因显著性是每一个发现的基因与表型相关的统计显著性，所述的非孟德尔表型显著性是当无法发现任何基因或所有已发现的基因都经由实验验证与表型无关时，纳入实验研究的个体实际上不是具有同一特定表型的个体，这一可能性的显著性。

4.根据权利要求3所述的基于概率框架和重测序技术快速发现表型相关基因的方法，其特征在于，所述的四个指标的概率计算框架式为：

$H(M,N)=P\left(\right|\mathrm{\&Phi;}|\&GreaterEqual;M)=\underset{\&ForAll;\mathrm{\&Phi;}:\left|\mathrm{\&Phi;}\right|\&GreaterEqual;M}{\&Sigma;}P\left(\mathrm{\&Phi;}\right)=\underset{\&ForAll;\mathrm{\&Phi;}:\left|\mathrm{\&Phi;}\right|\&GreaterEqual;M}{\&Sigma;}(\underset{S_i\&Element;\mathrm{\&Phi;}}{\mathrm{\&Pi;}}{D}_{S_i}\&times;\underset{S_i\&NotElement;\mathrm{\&Phi;}}{\mathrm{\&Pi;}}(1-D_{S_i}\left)\right)$

$\begin{array}{c}{Q}_{G_j}=P\left(\right|{\mathrm{\&Phi;}}_{G_j}|\&GreaterEqual;M)=\underset{\&ForAll;{\mathrm{\&Phi;}}_{G_j}:\left|{\mathrm{\&Phi;}}_{G_j}\right|\&GreaterEqual;M}{\&Sigma;}P\left({\mathrm{\&Phi;}}_{G_j}\right|)\\ =\underset{\&ForAll;{\mathrm{\&Phi;}}_{G_j}:\left|{\mathrm{\&Phi;}}_{G_j}\right|\&GreaterEqual;M}{\&Sigma;}\left(\underset{S_i\&Element;{\mathrm{\&Phi;}}_{G_j}}{\mathrm{\&Pi;}}\right(1-{(1-B_{S_i})}^{L_{G_j}})\&times;\underset{S_i\&NotElement;{\mathrm{\&Phi;}}_{G_j}}{\mathrm{\&Pi;}}{(1-B_{S_i})}^{L_{G_j}})\end{array}$

$E(M,N)=\underset{\&ForAll;G_j}{\&Sigma;}{Q}_{G_j}$

$\begin{array}{c}Z(w^{\&prime;},N)\&le;Z(w,N)=P\left(\right|\mathrm{\&Phi;}|\&le;w)={\mathrm{\&Sigma;}}_{\&ForAll;\mathrm{\&Phi;}:\left|\mathrm{\&Phi;}\right|\&le;w}P\left(\mathrm{\&Phi;}\right)={\mathrm{\&Sigma;}}_{\&ForAll;\mathrm{\&Phi;}:\left|\mathrm{\&Phi;}\right|\&le;w}({\mathrm{\&Pi;}}_{S_i\&Element;\mathrm{\&Phi;}}{D}_{S_i}\&times;{\mathrm{\&Pi;}}_{S_i\&NotElement;\mathrm{\&Phi;}}(1-\\ D_{S_i}\left)\right).\end{array}$

5.根据权利要求3所述的基于概率框架和重测序技术快速发现表型相关基因的方法，其特征在于，所述的四个指标基于递推思想的计算式为：

${\mathrm{PV}}_{G_j,S_i}=1-{(1-B_{S_i})}^{L_{G_j}}$

$P(A_u=v)=\left\{\begin{array}{cc}{D}_{S_u}\&times;P(A_{u-1}=v-1)+(1-D_{S_u})\&times;P(A_{u-1}=v)& (0<v<u)\\ {\mathrm{\&Pi;}}_{i=1}^u(1-D_{S_u})& (v=0)\\ {\mathrm{\&Pi;}}_{i=1}^u{D}_{S_u}& (v=u)\end{array}\right.$

$\begin{array}{c}P(R_{G_j,u}=v)\\ =\left\{\begin{array}{cc}{\mathrm{PV}}_{G_j,S_u}\&times;P(R_{G_j,u-1}=v-1)+(1-{\mathrm{PV}}_{G_j,S_u})\&times;P(R_{G_j,u-1}=v)& (0<v<u)\\ {\mathrm{\&Pi;}}_{i=1}^u\&lsqb;1-{\mathrm{PV}}_{G_j,S_u}\&rsqb;& (v=0)\\ {\mathrm{\&Pi;}}_{i=1}^u{\mathrm{PV}}_{G_j,S_u}& (v=u)\end{array}\right.\end{array}$

$H(M,N)=P(A_N\&GreaterEqual;M)=\underset{k=M}{\overset{N}{\&Sigma;}}P(A_N=k)$

$Q_{G_j}=P(R_{G_j,N}\&GreaterEqual;M)=\underset{k=M}{\overset{N}{\&Sigma;}}P(R_{G_j,N}=k)$

$E(M,N)=\underset{\&ForAll;G_j}{\&Sigma;}{Q}_{G_j}$

$Z(w^{\&prime;},N)\&le;Z(w,N)=P(A_N\&le;w)={\mathrm{\&Sigma;}}_{k=0}^{w}P(A_N=k).$

6.根据权利要求4或5所述的基于概率框架和重测序技术快速发现表型相关基因的方法，其特征在于，所述的计算式中，

S_i：一个样本；

G_j：一个基因；

基因组中基因G_j的区域长度；

一个分析方法在样本S_i中能够检测到与表型相关的那个或那些基因组变异的概率，表征单个样本的测序结果中观察到与所研究表型相关基因组变异的敏感度；

一个分析方法在经过基因组变异筛选的步骤以后在样本S_i上发现的基因组变异的频率，每个碱基有多大概率对应一个基因组变异，表征单个样本的测序结果中观察到与所研究表型相关基因组变异的特异性；

N：纳入研究的相同表型独立样本的总数；

M：研究人员指定的候选基因报告标准；候选基因必须在不少于M个相同表型独立样本中至少存在一个发现的基因组变异；

w′：当一个“基于基因组重测序技术的正向遗传学研究过程”无法发现可能与表型相关的基因或已发现的基因都已经被验证为假阳性发现时，真正与表型相关的基因在w′个相同表型独立样本中至少存在一个发现的基因组变异，w′<M；

w：当一个“基于基因组重测序技术的正向遗传学研究过程”无法发现可能与表型相关的基因或已发现的基因都已经被验证为假阳性发现时，不满足候选基因报告条件，未被证明与表型无关的基因，最多在w个相同表型独立样本中至少存在一个发现的基因组变异，w′≤w<M；

C_N：从N个样本中取出任意个样本的组合；

Φ：一种在样本中发现基因组变异的可能性。Φ代表一种样本组合：属于该组合的样本中，真正与表型相关的基因组变异被发现了，不属于该组合的样本中，真正与表型相关的基因组变异未被发现，Φ∈C_N；

一种在样本中发现基因组变异的可能性，代表一种样本组合：属于该组合的样本中，在G_j区域都发现了至少一个的基因组变异，不属于该组合的样本中，在G_j区域均没有发现的基因组变异，

||：计算组合中样本的个数；

H(M,N)：方法成功率指标；

候选基因显著性指标；

E(M,N)：假阳性发现数指标；

Z(w′,N)：非孟德尔表型显著性指标；

u,υ：递推变量，可在递推公式中作为样本号替代S_i的下标i；

随机情况下，样本S_i中基因G_j至少包含一个发现的基因组变异的概率；

A_u：N个样本中前u个样本中，真正与表型相关的基因组变异被发现的样本

数；

N个样本中前u个样本中，随机情况下，基因G_j至少包含一个(经过基因组变异筛选的步骤以后)发现的基因组变异的样本数。

7.根据权利要求4或5所述的基于概率框架和重测序技术快速发现表型相关基因的方法，其特征在于，所述的和通过以下计算式计算：

$D_{S_i}=H_{S_i}\&times;(1-{\mathrm{\&Phi;}}_{S_i});$

$B_{S_i}={\mathrm{Nv}}_{S_i}/l;$

其中

样本S_i的基因组变异检测灵敏度VCS；

样本S_i的筛选器误删率FIR；

样本S_i经过筛选器筛选以后，剩下的可能与表型相关的基因组变异个数；

l：实验研究所关注的基因组区域。

8.根据权利要求7所述的基于概率框架和重测序技术快速发现表型相关基因的方法，其特征在于，所述的通过以下计算式计算：

${\mathrm{\&Phi;}}_{S_i}=1-(1-{\mathrm{\&Phi;}}_{S_i,F})\&times;(1-{\mathrm{\&Phi;}}_{S_i,S}),$ 其中：

第一类筛选器的整体误删率；

第二类筛选器的整体误删率。

第一类筛选器的整体误删率通过用同样的一组第一类筛选器筛选所研究物种的所有已知基因组变异，评估多少已知能够产生非特定孟德尔表型的变异被筛除，来估计；

第二类筛选器的整体误删率通过计算所有二类筛选方法所认定的不可能基因组区域占全部实验研究所关注的基因组区域的比例，来估计。

9.根据权利要求8所述的基于概率框架和重测序技术快速发现表型相关基因的方法，其特征在于，筛选器误删率计算方法中，有效区筛选和显著差异筛选属于第一类筛选器；稠密度筛选、同祖先筛选、对照组筛选属于第二类筛选器。

10.根据权利要求1或2或3或4或5或8或9所述的基于概率框架和重测序技术快速发现表型相关基因的方法，其特征在于，所述的通过计算四个指标评估研究过程中每一个步骤的设计对研究整体有效性造成的影响，从而指导对当前实验方案和分析方法的优化的具体方法是：

1)、当一个基于基因组重测序技术的正向遗传学研究过程只发现了少数几个候选基因，且候选基因显著性较强，在经过严格的多重统计矫正后仍具有显著性时，不论方法成功率高低，实验验证发现的候选基因；

2)，当一个基于基因组重测序技术的正向遗传学研究过程找到的与表型相关的基因很多，超过了可以直接进行实验验证的合理数量时，参考假阳性发现数，假阳性发现数较低，则从候选基因显著性最强的候选基因开始，逐个进行实验验证，假阳性发现数较高，则说明所采用的数据分析方法标准过松，分析中发现了大量假阳性基因组变异，从而导致发现了很多的假阳性候选基因，对数据分析流程包括完成各个分析步骤的软件选择和参数设定，以及最终的候选基因报告标准进行优化，使用更为严格的基因组变异分析方法、基因组变异筛选方法和/或更严格的候选基因报告标准，对分析方法的不同优化有不同的有效性，最有效的优化方法应该是在不显著降低方法成功率的条件下，极大地降低假阳性发现数，研究人员可以尝试多种优化方法，评估优化后的分析方法的四个指标，然后选出最有效的优化方法；

3)、当一个基于基因组重测序技术的正向遗传学研究过程无法发现可能与表型相关的基因或已发现的基因都已经被验证为假阳性发现时，参考非孟德尔表型显著性，非孟德尔表型显著性不强，可能有两个原因：一是数据分析所采用的方法标准过严，导致与表型相关的突变未被发现，从而使表型相关基因未能在足够多的样本中被发现含有突变；二是样本量不够，无法有效发现表型相关基因；首先对数据分析流程包括完成各个分析步骤的软件选择和参数设定，以及最终的候选基因报告标准进行优化，使用更为宽松的基因组变异分析方法、基因组变异筛选方法和/或更宽松的候选基因报告标准，最有效的优化方法应该是在不显著增加假阳性发现数的条件下，极大地增加方法成功率；

4)、方法优化后可以找到符合情况一或情况二的候选基因，则进行实验验证；

5)、不论怎样优化方法，都无法发现符合情况一或情况二的候选基因，且“非孟德尔表型显著性”仍然不强，则考虑是由于纳入研究的相同表型独立样本数过少，建议增加样本量；

6)、经过了方法优化和/或进一步研究，仍无法发现与表型相关的基因或已发现的基因都已经被验证为假阳性发现，且非孟德尔表型显著性较强，首先检查纳入研究的样本的表型是否存在观察错误，如果没有错误，则测量更多的表型特征，以确认和/或细分纳入研究的样本的表型，然后依据更严格的表型标准来进行研究，由于基因通常具有多效性，两个功能类似的基因的变异可能产生类似的主要表型，但在次要表型上总会存在不同，利用更严格的表型标准，获得同一基因突变的表型个体，从而使得基于基因组重测序技术的正向遗传学研究过程可以正确发现这一表型相关基因。