基于EMS群体的基因水平全基因组关联分析算法

本发明涉及生物信息学，具体涉及基于ems群体的基因水平全基因组关联分析算法。

背景技术：

1、ems诱变（ethyl methanesulfonate, 甲基磺酸乙酯）是一种常用的化学诱变方法，用于通过化学物质引起生物dna的碱基突变。ems突变体库是通过使用ems对生物体（例如植物、动物、细菌等）进行突变处理，生成的大量突变体集合。ems突变体库的特点是通量高，突变分布广，但是ems诱变导致的突变位置通常是随机的。这导致每个碱基位点的最小等位基因频率(minor allele frequency，maf)达不到进行全基因组关联研究（genome-wide association studies，gwas）的要求，给gwas分析带来困难。而gwas是在全基因组范围内找出snp,筛选出与性状相关的snp,研究遗传突变和表型之间的相关性，可用于研究作物的产量、抗病性等重要性状,为农业育种提供重要依据。

2、在人类疾病研究中gwas分析，主要是利用大样本量，寻找病例对照组中常见变异(maf>1%)的差异来解释复杂疾病，忽略了罕见变异(maf<1%)的作用与信息，而罕见变异，被认为能解释一部分疾病的遗传关系，发生机制。所以开发了基于区域/基因的关联分析方法负荷检验(burden analysis)。负荷检验简单的说就是比较两组表型差异的样本，在同一个区域/基因上，携带的罕见变异的“总数”是不是有显著差异。

3、但是负荷检验无法应用于ems突变体库中。首先人类的罕见变异数量远低于ems位点的数量，相差百万级，ems突变体库中的个体上千株也远大于人类队列的几百例，负荷检验无法负担大规模的分析。其次疾病的表型是质量性状，而作物的表型是数量性状，使用负荷检验分析数量性状会使分析效力大打折扣。再次，负荷检验统计的是统计区域的突变数量，没有合理的对不同功能的突变位点进行加权，导致分析效力降低。最后，负荷检验还需要大量的对照样本，增加了分析门槛。

技术实现思路

1、针对现有技术中的上述不足，本发明提供了基于ems群体的基因水平全基因组关联分析算法，本发明可以直接精准定位到单个基因，提高了候选基因功能验证实验的效率，有效解决了现有技术无法负担大规模的分析、分析效力低以及分析门槛高的问题。

2、为实现上述目的，本发明解决其技术问题所采用的技术方案是：提供基于ems群体的基因水平全基因组关联分析算法，包括以下步骤：

3、s1、获取待分析样本进行外显子测序；

4、s2、获取ems样本的表型数据；

5、s3、将ems群体的外显子测序数据与其物种的基因组进行比对；

6、s4、对得到的比对结果使用snp检测工具进行snp检测，获得基因编码区域的snp信息；

7、s5、使用snpeff软件对snp进行注释和突变效应预测，获得每个snp的变异对基因功能的影响程度；

8、s6、统计每个样本的每个基因区域的受到的总突变影响，获得基因水平的受影响程度信息；

9、s7、使用线性模型和卡方检验进行基因型和表型的关联分析。

10、进一步，步骤s3中，使用bwa、bowtie2等比对软件对步骤s1中的外显子测序数据与其物种全基因组进行比对。

11、进一步，步骤s4中，对得到的比对结果使用gatk、samtools、bcftools等snp检测工具进行snp检测。

12、进一步，步骤s4中，snp信息保存为vcf文件格式。

13、进一步，步骤s4中，具体方法为：从vcf文件中提取所有样本的snp位点的基因型信息，得到snp位点的分型矩阵，根据排列检测确定的突变频率阈值过滤突变频率过高的异常位点。

14、进一步，将vcf文件中每一个snp的基因型信息与样本信息打乱，得到随机排列的样本基因型信息，统计每个位点的突变频率，从低到高排列，选取第5低的频率值，作为一次检验的结果，共进行n次，将n次结果从低到高排列，选取第n*0.01低的结果作为突变频率阈值。

15、进一步，步骤s5中，具体方法为：使用注释工具snpeff注释每个位点所在的具体基因，设定距离阈值k，k设置为小于3000的正数，获取基因上以及距离基因上下游k bp以内的snp位点作为基因区域内的snp位点，并预测每个snp的影响效应，影响效应的等级分为high、moderate、low和modifier。

16、进一步，步骤s7中，进行基因型和表型的关联分析，具体包括以下步骤：

17、s71、过滤假阳性snp；

18、s72、过滤表型异常值；

19、s73、建立基因水平突变效应与表型的正向统计模型，进行基因型与表型的关联分析，获取表型与基因的关联p值；

20、s74、建立突变表型与基因水平突变效应的反向统计模型，进行表型与基因型的关联分析，获取基因与表型的关联p值；

21、s75、将正向与反向结果进行合并，获得最终表型与基因的关联分数；

22、s76、依据统计检验原理，使用排列检测得到可靠的显著性阈值；

23、s77、筛选显著关联的基因。

24、进一步，步骤s73中，进行基因型和表型的关联分析，具体包括以下步骤：

25、s731、依据位点突变的影响效应等级对每个位点进行赋分，其中，modifier=1.000001，low=2.00001，moderate=3.0001，high=4.001；

26、s732、依据样本的基因型信息和基因位置注释，统计每个样本的每个基因的总突变影响效应分数，例如一个样本的a基因上有3个low水平突变和2个high水平突变，那么a基因的总突变效应分数就为14.00203；

27、s733、建立基因水平的基因型与表型的线性回归统计模型，进行基因型与表型的关联分析，获取表型与各基因的关联p值，选取前0.5%作为显著关联基因；

28、s734、建立基因突变情况与表型情况的检验模型，进行基因突变与表型异常的关联分析，获取表型与各基因的关联p值，选取前0.5%作为显著关联基因；

29、s735、将步骤s33和步骤s34的显著关联基因互相验证，合并取交集作为候选关联基因；

30、s736、通过排列检验确定显著性阈值，过滤显著性小于阈值的候选关联基因，得到最终结果。

31、进一步，步骤s733中，对于单个基因a来说，计算每个样本a基因的总突变效应分数，作为线性模型x，表型值作为y，对整个ems群体进行y=ax+b的线性回归拟合，确定线性模型拟合系数的显著性p值，作为a基因和表型关联性的p值。对所有基因进行统计，得到最终结果。

32、进一步，步骤s734中，对于a基因来说，通过箱线图确定表型显著差异的样本，统计数量，同时统计在a基因上有突变的样本数量和没有突变的样本数量，通过卡方检验得到表型异常的样本数量和a基因上有突变的样本数量的显著性关系，作为a基因和表型关联性的p值。对所有基因进行统计，得到最终结果。

33、进一步，步骤s736中，将群体的表型信息与样本信息打乱后，再进行s734和s735的关联分析，记录显著性最高的第五个基因p值。这样重复n次，将n次结果从低到高排列，选取第n*0.01低的结果作为显著性阈值。

34、进一步，步骤s77中，筛选显著关联的基因，具体包括以下步骤：

35、s771、设定距离阈值k，k设置为小于3000的正数，获取基因上以及距离基因上下游k bp以内的snp位点；

36、s772、使用模拟数据预测snp在群体中的突变频率阈值，过滤超过阈值的snp位点；

37、s773、根据步骤s73和步骤s74计算出的各基因关联p值，进行综合计算，得到基因最终关联分数；

38、s774、确定关联分数的阈值p0，筛选出关联分数大于p0的基因。

39、本发明具有以下有益效果：

40、1、ems突变群体的基因型有稀疏性和随机性的特点，稀疏性使得单个位点的maf过低，以至于无法利用现有的关联分析技术进行基因型-表型关联分析。将关联分析的基因型基本单位，从snp提高至基因水平，极大的提升了检测的有效性。同时本发明不是简单的统计基因有无突变，而是根据基因的突变效应给基因的整体效应进行打分，以此构建线性回归模型，提高了检测的准确性。

41、2、本发明还利用突变富集检验反向验证结果，降低了结果的错误率。根据数据分布和模型的复杂度，没有使用传统的临界值作为阈值，而是使用排列检测结果作为阈值，防止结果的假阳性。从基因水平去进行关联分析，得到的结果也直接是具体的基因，而非传统关联分析的基因组区域，极大减少了人力物力和时间成本。

42、3、以基因为基本单位进行关联分析，maf指标相较于碱基水平有极大提高。根据每个突变位点的作用进行加权，统计单个样本的单个基因中所有突变的总加权值，作为关联分析的基础，极大的减小了结果的假阳性率，提高了分析的统计效力。同时使用多种统计方法进行综合评判，以找到和表型最相关的基因。同时，相对于gwas只能定位到一个模糊的区间，结果可能包含多个基因，本发明可以直接精准定位到单个基因，提高了候选基因功能验证实验的效率。