CN104651517A - 一种基于snpldb标记的限制性二阶段全基因组关联分析方法 - Google Patents

Abstract

本发明公开了一种基于SNPLDB标记的限制性二阶段全基因组关联分析方法，以解决传统方法无法估计复等位基因信息、假阳性率高以及在近交作物中检测功效低的问题。本发明结合基于单倍型区块构建的SNPLDB标记、近交群体关联分析模型偏差的矫正和多位点模型下二阶段关联分析策略，建立了适合于近交作物常规育种的GWAS方法。该方法将SNPLDB标记用于GWAS，为复等位基因估计提供了方法，第一阶段基于单位点模型来筛选候选位点，第二阶段基于多位点模型下的逐步回归分析方法作进一步筛选以平衡缺失遗传率和遗传率估计过高的问题，从而将最终遗传模型的解释率控制到性状遗传率。GWAS使用由SNPLDB标记估计的相似系数矩阵的特征向量和合适的显著水平来提高定位的准确性和功效。

Description

一种基于SNPLDB标记的限制性二阶段全基因组关联分析方法

技术领域

本发明属于分子数量遗传和分子育种技术领域，具体涉及一种基于SNPLDB标记的限制性二阶段全基因组关联分析方法。

背景技术

绿色革命通过群体株型改造育成高产品种是常规育种技术发展的成功。常规育种主要是创造超亲重组型，但以往育种家只能通过表型评估来设计组合并间接追踪超亲重组型。分子技术的快速发展为直接鉴别超亲重组基因型提供了可能。

基于种质资源群体的全基因组关联分析(GWAS)为全面解析农艺性状的遗传基础提供了方法。以往植物中的GWAS研究主要目的是发掘主效基因，但是通过提高显著水平来尽可能降低假阳性的同时却导致了检测功效大大降低。植物育种家为了准确检测全基因组数量性状基因座(QTL)，他们的GWAS策略必须符合种质资源中广泛存在复等位基因的情况，并需要在缺失遗传率和遗传率过高估计之间进行平衡，以及矫正由近交和迁移导致的群体偏差。

现有GWAS广泛使用的单核苷酸多态性(SNP)分子标记仅有两个等位变异，无法估计资源群体中大量存在的复等位基因效应，这进一步限制了其在常规育种中的应用。另外，GWAS的精度依赖于连锁不平衡(LD)的衰减距离。随机交配群体的LD衰减距离通常较短，GWAS的精度也较高。但自花授粉作物自然群体往往严重偏离随机交配群体，高度自交导致了群体较长的LD衰减距离，GWAS的精度也随之降低。

GWAS一直饱受高假阳性的困扰，现有GWAS方法研究也主要针对如何通过控制群体结构来降低假阳性。群体结构推断和主成分分析是GWAS中两种广泛使用的降低假阳性的方法，这两种方法均通过将推断的群体结构特征作为协变量引入GWAS统计模型以降低群体结构的影响。但已有报道显示群体内个体间的亲缘关系也会导致GWAS中的假阳性，GWAS应同时考虑群体结构和亲缘关系。相应地，研究者提出了一系列基于混合线性模型(LMM)的GWAS方法将个体间两两亲属关系考虑进来。LMM方法假定每个个体均抽样自不同的群体，并将群体背景作为随机效应，群体结构作为固定效应拟合到GWAS统计模型，并将群体的亲本系数(kinship)矩阵作为随机效应的协方差结构。目前，LMM方法被认为比基于群体结构和PCA的GWAS假阳性更低，已被广泛应用于动植物遗传研究。但是，大量基于LMM的GWAS研究结果仅能检测到少数几个位点，并且仅解释表型变异的很少部分，而实际上数量性状通常由许多效应大小不等的位点控制。因此，较高的假阴性率导致了GWAS的失踪遗传率问题，即关联位点的总遗传贡献率远低于性状遗传率。LMM方法主要依赖的kinship矩阵通常是由全基因组分子标记估计得到的，该矩阵实际上是个体间状态同样的估计。实际研究通常使用所有的分子标记来估计这种kinship矩阵，但是最新研究表明kinship估计中所用分子标记如果包含遗传位点SNP标记时，GWAS的功效将会降低，反之，假阳性会升高。研究者也相应的提出了几种方法来解决基于LMM方法的GWAS中kinship敏感的问题，但是对于解决失踪遗传率问题帮助甚微。

发明内容

本发明目的在于提供一种基于SNPLDB标记的限制性二阶段全基因组关联分析方法，该方法结合基于单倍型区块构建的SNPLDB标记、近交群体关联分析模型偏差矫正和多位点模型下二阶段关联分析策略，建立了适合于近交作物常规育种应用的GWAS方法。

本发明解决其技术问题所采取的技术方案是：一种基于SNPLDB标记的限制性二阶段全基因组关联分析方法，包括如下步骤：

a)构建全基因组SNPLDB标记：首先对已有种质资源群体的全基因组SNP基因型数据进行连锁不平衡分析，然后利用Haploview软件定义全基因组SNP分子标记的单倍型区块，阈值为D'>0.7，窗口设为估计的LD衰减距离；最后将单倍型区块内的SNP分子标记合并为新的标记，即SNPLDB标记，也就是将单倍型区块内的每一个单倍型视为位点的一个等位变异并进行编码；对于频率较低的单倍型(频率小于1％)，通过单倍型的聚类分析使用最为相似的单倍型替换低频率的单倍型。

近交群体较长的LD衰减距离说明基因型呈现出区块结构，即基因组可以分割为长度不等的区块，区块间重组频繁，区块内重组较少，一般的解释是染色体上存在着大量的重组热点。本发明发现由于重组热点区域频繁的重组导致了重组不频繁的区域呈现出区块的模式。这种区块内丰富的单倍型变异提供了类似复等位基因的变异特征，相比SNP分子标记更符合具有广泛遗传变异的种质资源群体的遗传特征，可利用此特性来估计自然群体中的复等位基因信息。

单倍型区块定义所用的Haploview软件参数为：-minMAF 0.01、-hwcutoff 0、-maxDistance200、-blockoutput GAB、-blockMAFThresh 0.01，其中-maxDistance 200指窗口大小为200kb；

b)近交群体关联分析模型偏差的矫正：直接基于构建的全基因组SNPLDB标记，计算其遗传相似系数矩阵作为亲属关系的估计，这种方法一定程度上反映了群体内个体间的遗传关系以及自然群体偏离随机交配群体的程度，可以用于矫正由近交导致的GWAS模型偏差。

假定群体包含n个个体，关联分析中对单个标记位点的假设测验的线性模型一般可以表示为

$y_i=\mathrm{\μ}+\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{w}_{\mathrm{ij}}{\mathrm{\α}}_j+\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{x}_{\mathrm{il}}{\mathrm{\β}}_l+{\mathrm{\ϵ}}_i---\left(1\right)$

其中y_i为第i个个体的表型观测值，μ为群体平均数，w_ij为反应第i个个体与第j个亚群体的关系的系数，α_j为第j个亚群体效应，x_il取值为0或1，表示第i个个体在标记位点上的基因型，如果个体在标记位点的基因型为第l个等位基因，则x_il为1，反之为0，β_l为标记位点第l个等位基因的效应并假定对于单个位点有Σβ_l＝0，ε_i为残差效应并服从N(0,σ²)，σ²为误差方差。

在近交群体中，准确估计模型(1)的群体结构效应(α_j)是提高关联分析功效和降低假阳性的关键。已有的方法不适合本发明构建的SNPLDB标记，因此我们提出利用全基因组SNPLDB标记计算的个体间遗传相似系数作为亲本系数的估计用于关联分析群体偏差的矫正，并记该方法为EigenIBS。

二倍体群体中，基于SNPLDB标记的个体间遗传相似系数可以简单定义为状态同样位点的比例，即Σn_k/2m，取值范围为[0,1]，其中n_k为在第k个标记上两个体共有的等位基因数目，m为总标记数目。这种遗传相似系数非常容易计算，并且适用于任何分子标记，在群体遗传研究中已有广泛应用。对于包含n个个体的群体，该相似系数矩阵为一个n×n的对称矩阵，使用该遗传相似矩阵的部分特征向量作为群体结构的估计用于关联分析群体偏差的矫正。

c)多位点模型下二阶段关联分析：第一阶段基于单位点模型使用0.05的显著水平进行筛选候选位点以排除大量无用干扰信息；第二阶段以控制遗传贡献总量为目的，基于多元逐步回归分析方法构建包含多个位点的最终遗传模型以控制总的表型变异解释率，两个阶段中均使用相似系数矩阵特征向量矫正由近交导致的模型偏差。。

在模型(1)的基础上，将关联分析模型拓展为多位点模型如下

$y_i=\mathrm{\μ}+\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{w}_{\mathrm{ij}}{\mathrm{\α}}_j+\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{l=1}{\overset{L_k}{\mathrm{\Σ}}}{x}_{\mathrm{ikl}}{\mathrm{\β}}_{\mathrm{kl}}+{\mathrm{\ϵ}}_i---\left(2\right)$

其中K为控制数量性状总位点数，其余符号与模型(1)相同。

本发明使用如下二阶段策略的多位点关联分析方法从所有标记中筛选显著关联的标记：第一阶段使用p＝0.05的显著水平基于模型(1)的EigenIBS方法对所有标记进行初步筛选，筛选到的标记将作为候选位点纳入第二阶段分析；第二阶段使用标准的多元逐步回归方法对候选位点进行第二轮筛选，逐步回归中使用EigenIBS方法来矫正群体偏差。大致过程为：首先在已入选位点条件计算所有剩余位点显著性测验的p值，然后选择低于显著水平p值最小的位点作为新位点加入模型，并重复该过程直到没有显著的位点可供选择。最后根据模型(2)拟合包含所有入选位点的遗传模型，删除大于显著水平p值最大的位点，并重复该过程直到模型中所有位点均显著。

由于涉及非常多的假设测验，上述二阶段策略中多元逐步回归方法通常具有严重的多重测验问题，进而可能导致模型过拟合，因此提出3种矫正多元逐步回归多重测验的方法。第一种方法是使用Bonferroni矫正后的显著水平(α/m)作为多元逐步回归中每一步的显著水平，其中α和m分别为显著水平和候选标记个数。第二种方法基于Holm矫正方法将多元逐步回归中第i步的显著水平调整为α/(m-i+1)。第三种方法基于假发现率(FDR)控制过程将多元逐步回归中第i步的显著水平调整为αi/m。

第二阶段分析中基于多重测验矫正的逐步回归较为严格，最终模型的表型解释率可能远低于性状遗传率，导致缺失遗传率问题，因此为将GWAS应用于作物常规育种，可以将性状遗传率作为最终模型表型解释率的上限，通过禁用多重测验矫正和调整显著水平来进行模型选择，使得GWAS能够解析全基因组的遗传位点。

本发明所述的SNPLDB指单核苷酸多态性分子标记连锁不平衡区块。

有益效果：

1、全基因组SNPLDB标记的构建，缩短了LD衰减距离，提高了近交作物GWAS的检测精度，并且SNPLDB标记具有复等位变异特性，因此提供了检测自然群体中广泛存在的复等位基因以及估计其效应的方法。

2、直接利用群体的遗传相似系数进行近交群体关联分析模型的矫正，遗传相似系数计算简单，且不受限于标记类型，适用于SNPLDB标记，该方法降低了由高度近交导致的近交作物GWAS的假阳性。

3、限制性二阶段分析策略大幅减少了运算量，通过第一阶段分析排除了大量的无用标记，减少了第二阶段分析的多重测验问题，基于多位点模型的分析方法大幅提高了GWAS的检测功效，使得GWAS能够有效应用于常规育种的标记辅助选择，为全面解析近交作物数量性状遗传基础及其常规育种应用提供了理论依据和技术方法，建立了适合于近交作物常规育种应用的GWAS方法。

附图说明

图1为本发明的方法流程图。

具体实施方式

下面结合说明书附图对本发明创造作进一步说明。

如图1所示，一种基于SNPLDB标记的限制性二阶段全基因组关联分析方法的流程图，设计了该方法的实施过程，该方法包括全基因组SNPLDB标记的构建、近交群体关联分析模型的矫正及多位点模型下的二阶段关联分析这三个步骤。

a)全基因组SNPLDB标记的构建

分布于全基因组的145558个SNP分子标记基因型数据为简化基因组测序(RAD-seq)后经组装及一系列质量控制程序后获得，每个SNP的最小等位基因频率大于1％。利用Haploview软件定义了SNP分子标记的单倍型区块，阈值设为D'>0.7，窗口设为200kb，基于包含1024份材料的中国大豆种质资源群体分布全基因组的145558个SNP分子标记，构建了分布全基因组的36952个SNPLDB标记，每个标记呈现2到14个变异类型，可见SNPLDB提供了复等位基因检测方法。比较了SNP和SNPLDB两种不同标记下群体LD的变化，结果(表1)显示基于SNP标记估计的平均LD衰减至0.6的距离约为3Mb，而对于SNPLDB标记的距离则缩短至500kb左右。由此可见，SNPLDB标记通过合并单倍型区块内的SNP位点不仅缩短了近交群体LD的衰减距离，而且提供了比SNP标记更丰富的适合育种应用的复等位变异信息，更适合应用于近交作物GWAS。

表1SNP和SNPLDB标记关于LD的比较

注：SNPLDB表示SNP连锁不平衡区块，LD采用D′度量，距离区间为物理位置区间，Mean为区间内D′平均数，Q1、Q2和Q3分别区间内D’的第1、2、3分位数。

b)近交群体关联分析模型的矫正

基于全基因组的36952个SNPLDB标记，计算了群体内个体间两两遗传相似系数，并相应构建群体遗传相似矩阵。利用线性代数标准的特征值分解方法，提取遗传相似矩阵的前10个特征向量并用于GWAS分析中作为协变量以控制群体结构。

本实施例使用了大豆种质资源群体SNP分子标记数据的100位点的模拟表型数据比较分析了EigenIBS与ADMIXTURE、EIGENSTRAT和LMM方法的在关联分析中的假阳性及检测功效。ADMIXTURE与STRUCTURE方法的基本原理相同，但计算更快速。LMM方法采用TASSEL软件的MLM功能进行计算。从模拟结果(表2)可以看出，相对与未经群体偏差矫正的方法，四种群体偏差矫正方法均有效降低了关联分析的假阳性，但是检测功效也大大降低，这说明群体偏差矫正对检测功效有非常大的影响，特别是LMM方法的检测功效仅为6.7。另外除了LMM方法的假阳性率最低外，其余三种方法均具有较高的假阳性。EigenIBS虽然假阳性略高于ADMIXTURE和EIGENSTRAT方法，但检测功效也略高，因此该方法能够有效降低由于群体偏差引起的假阳性。

表2五种矫正GWAS模型偏差方法的功效比较

注：括号内数值为100次模拟结果的标准差，方法不进行任何偏差矫正，ADMIXTURE分析设定了14个亚群，EIGENSTRAT及EigenIBS中均使用10个特征向量，LMM方法中未使用TASSEL软件MLM的压缩选项。使用Bonferroni校正的0.05作为显著水平。

以上模拟中EigenIBS使用了10个特征向量进行群体偏差矫正，为进一步了解特征向量个数对关联分析结果的敏感程度，我们进一步分别计算了1、5、10、20和40个特征向量时的关联分析的假阳性和检测功效。从模拟结果(表3)可以看出，当仅使用一个特征向量时，检测功效较高，随着特征向量个数的增加，假阳性和检测功效均逐渐下降，最终趋于稳定，这说明遗传相似矩阵信息可有少量的特征向量代替，因此当基本解析了大部分遗传相似矩阵的信息后，特征向量个数对关联分析的检测功效变得并不敏感，也就是说使用较多的特征向量一般不会明显降低检测功效。

表3不同特征向量数目时GWAS功效

注：括号内数值为100次模拟结果的标准差，使用Bonferroni校正的0.05作为显著水平。

c)多位点模型下的二阶段关联分析

从国家大豆改良中心获得的代表中国大豆种质资源的1024份大豆材料的百粒重表型数据，供试大豆材料包含野生大豆、地方大豆和育成品种三大类型。百粒重表型数据经2010、2011、2012这3年田间试验鉴定获得，采用随机区组试验设计。利用SAS/STAT软件的GLM过程对以上百粒重表型数据进行了方差分析，并估计了该性状的遗传率。

全基因组关联分析的第一阶段基于单位点模型对全基因组SNPLDB标记筛选与百粒重关联的位点，使用0.05作为筛选的阈值；第二阶段基于多位点模型使用逐步回归方法对SNPLDB候选标记位点作进一步筛选，第二阶段根据不同的目的设置合适的显著水平。如果以基因发掘为目的，则可使用本发明提供的三种多重测验方法。如果以检测全基因组QTL和常规育种应用为目的，则应设置一系列显著水平，并选择不同显著水平下模型的表型解释率接近性状遗传率的模型为最终模型。

使用基于大豆种质资源群体SNPLDB分子标记模拟的100位点表型数据对方法进行了评估。表4列出了基于不同多重测验过程的新方法的功效，可以看出假阳性从第一阶段的7852.63分别降低至了9.49和21.76，假发现率与LMM方法相当，而新方法的功效明显高于LMM方法的7.24，分别为17.45和24.53。因此新方法更适合发现更多的位点，而LMM方法适合发掘个别基因。此外，对于不同的显著水平矫正方法，新结果差异不大，基于Bonferroni和Holm矫正方法的新结果一样，基于FDR方法的新功效较高，假阳性也略高。

表4新方法与LMM方法的功效比较

注：所有方法的显著水平均设为0.05，LMM和NEW(Bonf)方法中显著水平使用Bonferroni矫正，NEW(Holm)和NEW(FDR)方法中显著水平分别使用Holm和FDR方法矫正。

本发明不限于上述实施例，一切采用等同替换或等效替换形成的技术方案均属于本发明要求保护的范围。

Claims (9)

1.一种基于SNPLDB标记的限制性二阶段全基因组关联分析方法，其特征在于，所述方法包括如下步骤：

a)构建全基因组SNPLDB标记：首先对全基因组分子标记进行单倍型区块分析，根据连锁不平衡分析结果对全基因组标记进行单倍型区块的定义，然后根据单倍型区块通过合并区块内的标记将原始分子标记数据转换为SNPLDB标记；

b)近交群体关联分析模型偏差的矫正：直接基于构建的全基因组SNPLDB标记，计算其遗传相似系数矩阵作为亲属关系的估计，用于矫正由近交导致的GWAS模型偏差；

c)多位点模型下二阶段关联分析：第一阶段基于单位点模型使用0.05的显著水平进行筛选候选位点；第二阶段基于多元逐步回归分析方法构建包含多个位点的最终遗传模型。

2.如权利要求1所述的一种基于SNPLDB标记的限制性二阶段全基因组关联分析方法，其特征在于，所述构建全基因组SNPLDB标记的具体步骤为：首先对已有种质资源群体的全基因组SNP基因型数据进行连锁不平衡分析，然后利用Haploview软件定义全基因组SNP分子标记的单倍型区块，阈值为D'>0.7，窗口设为估计的LD衰减距离；最后将单倍型区块内的SNP分子标记合并为新的标记SNPLDB，就是将单倍型区块内的每一个单倍型视为位点的一个等位变异并进行编码，对于频率小于1％的单倍型，通过单倍型的聚类分析使用最为相似的单倍型替换低频率的单倍型。

3.如权利要求1所述的一种基于SNPLDB标记的限制性二阶段全基因组关联分析方法，其特征在于，所述近交群体关联分析模型偏差的矫正的具体步骤为：直接基于构建的全基因组SNPLDB标记，计算其遗传相似系数矩阵作为亲属关系的估计，该方法为EigenIBS；二倍体群体中，基于SNPLDB标记的个体间遗传相似系数简单定义为状态同样位点的比例，即Σn_k/2m，取值范围为[0,1]，其中n_k为在第k个标记上两个体共有的等位基因数目，m为总标记数目；对于包含n个个体的群体，该相似系数矩阵为一个n×n的对称矩阵，使用该遗传相似矩阵的部分特征向量作为群体结构的估计用于关联分析群体偏差的矫正。

4.如权利要求1所述的一种基于SNPLDB标记的限制性二阶段全基因组关联分析方法，其特征在于，所述多位点模型下二阶段关联分析的具体步骤为：第一阶段基于模型的EigenIBS方法，使用p＝0.05的显著水平对所有标记进行初步筛选，筛选到的标记作为候选位点纳入第二阶段分析；第二阶段使用标准的多元逐步回归方法对候选位点进行第二轮筛选，逐步回归中使用EigenIBS方法来矫正群体偏差：首先在已入选位点条件计算所有剩余位点显著性测验的p值；然后选择低于显著水平p值最小的位点作为新位点加入模型，并重复该过程直到没有显著的位点可供选择；最后根据模型拟合包含所有入选位点的遗传模型，删除大于显著水平p值最大的位点，并重复该过程直到模型中所有位点均显著。

5.如权利要求1所述的一种基于SNPLDB标记的限制性二阶段全基因组关联分析方法，其特征在于：所述的SNPLDB标记是指单核苷酸多态性分子标记连锁不平衡区块。

6.如权利要求1所述的一种基于SNPLDB标记的限制性二阶段全基因组关联分析方法，其特征在于：所述单倍型区块定义所用的Haploview软件参数为：-minMAF 0.01、-hwcutoff 0、-maxDistance 200、-blockoutput GAB、-blockMAFThresh 0.01，其中-maxDistance 200指窗口大小为200kb。

7.如权利要求1所述的一种基于SNPLDB标记的限制性二阶段全基因组关联分析方法，其特征在于，所述多位点模型下二阶段关联分析的两个阶段中均使用相似系数矩阵特征向量矫正由近交导致的模型偏差。

8.如权利要求1所述的一种基于SNPLDB标记的限制性二阶段全基因组关联分析方法，其特征在于，所述多位点模型下二阶段关联分析过程中有3种矫正多元逐步回归多重测验的方法：使用Bonferroni矫正后的显著水平α/m作为多元逐步回归中每一步的显著水平，其中α和m分别为显著水平和候选标记个数或基于Holm矫正方法将多元逐步回归中第i步的显著水平调整为α/(m-i+1)或基于假发现率控制过程将多元逐步回归中第i步的显著水平调整为αi/m。

9.如权利要求1所述的一种基于SNPLDB标记的限制性二阶段全基因组关联分析方法，其特征在于：所述二阶段关联分析的第二阶段将性状遗传率作为最终模型表型解释率的上限，通过禁用多重测验矫正和调整显著水平来进行模型选择，使得GWAS能够解析全基因组的遗传位点。