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

技术领域

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

背景技术

测序技术飞速发展，在精度不断提高的同时成本不断降低。当前已普及的测序技术 可以达到在一个月内完成一个人类个体的全基因组重测序，产生碱基准确度在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|\Phi \right|\ge M\right)=\underset{\forall \Phi :\left|\Phi \right|\ge M}{\Sigma }P\left(\Phi \right)=\underset{\forall \Phi :\left|\Phi \right|\ge M}{\Sigma }\left(\underset{S_i\in \Phi }{\Pi }{D}_{S_i}\times \underset{S_i\notin \Phi }{\Pi }\left(1-D_{S_i}\right)\right)$

$\left(\begin{array}{c}{Q}_{G_j}=P\left(\left|{\Phi }_{G_j}\right|\ge M\right)=\underset{\forall {\Phi }_{G_j}:\left|{\Phi }_{G_j}\right|\ge M}{\Sigma }P\left({\Phi }_{G_j}\right)\\ =\underset{\forall {\Phi }_{G_j}:\left|{\Phi }_{G_j}\right|\ge M}{\Sigma }\left(\underset{S_i\in {\Phi }_{G_j}}{\Pi }\left(1-{\left(1-B_{S_i}\right)}^{L_{G_j}}\right)\times \underset{S_i\notin {\Phi }_{G_j}}{\Pi }{\left(1-B_{S_i}\right)}^{L_{G_j}}\right)\end{array}\right)$

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

$Z(w^{\prime },N)\le Z(w,N)=P\left(\right|\Phi |\le w)={\Sigma }_{\forall \Phi :\left|\Phi \right|\le w}P\left(\Phi \right)={\Sigma }_{\forall \Phi :\left|\Phi \right|\le w}\left({\Pi }_{S_i\in \Phi }{D}_{S_i}\times {\Pi }_{S_i\notin \Phi }\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\ge M)=\underset{k=M}{\overset{N}{\Sigma }}P(A_N=k)$

$Q_{G_j}=P(R_{G_j,N}\ge 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-{\Phi }_{S_i});$

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

其中：

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

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

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

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

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

${\Phi }_{S_i}=1-(1-{\Phi }_{S_i,F})\times (1-{\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|\Phi \right|\ge M\right)=\underset{\forall \Phi :\left|\Phi \right|\ge M}{\Sigma }P\left(\Phi \right)=\underset{\forall \Phi :\left|\Phi \right|\ge M}{\Sigma }\left(\underset{S_i\in \Phi }{\Pi }{D}_{S_i}\times \underset{S_i\notin \Phi }{\Pi }\left(1-D_{S_i}\right)\right)$

$\left(\begin{array}{c}{Q}_{G_j}=P\left(\left|{\Phi }_{G_j}\right|\ge M\right)=\underset{\forall {\Phi }_{G_j}:\left|{\Phi }_{G_j}\right|\ge M}{\Sigma }P\left({\Phi }_{G_j}\right)\\ =\underset{\forall {\Phi }_{G_j}:\left|{\Phi }_{G_j}\right|\ge M}{\Sigma }\left(\underset{S_i\in {\Phi }_{G_j}}{\Pi }\left(1-{\left(1-B_{S_i}\right)}^{L_{G_j}}\right)\times \underset{S_i\notin {\Phi }_{G_j}}{\Pi }{\left(1-B_{S_i}\right)}^{L_{G_j}}\right)\end{array}\right)$

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

$Z(w^{\prime },N)\le Z(w,N)=P\left(\right|\Phi |\le w)=\underset{\forall \Phi :\left|\Phi \right|\le w}{\Sigma }P\left(\Phi \right)=\underset{\forall \Phi :\left|\Phi \right|\le w}{\Sigma }\left(\underset{S_i\in \Phi }{\Pi }{D}_{S_i}\times \underset{S_i\notin \Phi }{\Pi }\left(1-D_{S_i}\right)\right)$

其中：

S_i：一个样本。

G_j：一个基因。

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

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

一个分析方法在经过基因组变异筛选的步骤以后在样本S_i上发现的基因组变异的频 率每个碱基有多大概率对应一个基因组变异。英文描述为Backgroundvariant frequency，以下文章中简记为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\ge M)=\underset{k=M}{\overset{N}{\Sigma }}P(A_N=k)$

$Q_{G_j}=P(R_{G_j,N}\ge 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)和(基因组变异的频率，BVF backgroundvariantfrequency)需要通过实际测序数据估计，本专利设计了以下的估计 方法：

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

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

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

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

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

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

在之前的基因组重测序应用中，有一批常用的筛选与表型相关基因组变异的方法， 主要分为五类，分别是，有效区筛选(Effectiveregionfilter)，显著差异筛选(Big differencefilter)，稠密度筛选(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/pseudom olecules/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文件中突变位点所在的染色体、序列位置等信息，使用有效区筛选(Effective regionfilter，设定为筛除启动子区(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基因确实与植物对磷浓度的耐受力有关。

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