一种针对混杂人群的通用全基因组关联分析方法

本发明涉及系统生物学和生物统计学，涉及一种生物医学大数据分析处理方法，特别涉及一种针对混杂人群的通用全基因组关联分析方法。背景技术：1、全基因组关联分析(genome-wide association studies,gwas)旨在在人类全基因组范围内寻找与特定性状显著相关的序列变异。近年来，gwas开始出现一些新特点。随着基因测序和电子健康记录技术的进步，生物样本库积累了大量健康数据，为gwas提供了重要数据源。比如英国生物样本库收集了50万英国人的健康数据，包括基因测序、磁共振扫描成像、临床指标、生活方式等多维度、跨尺度数据，将广泛的遗传信息电子健康记录(ehr)与健康信息相结合，使gwas研究不再局限于传统的连续型性状或二分类性状，具有更加复杂结构的性状开始应用于gwas。与传统的简单性状(连续型性状或二分类性状)相比，复杂性状能够提供更丰富的信息，从而产生更有价值的结果，因此在大规模的全基因组关联分析研究中愈发受到关注。例如，使用生存数据除了可以刻画事件是否发生、发生的时间之外，还可以刻画数据的删失和截尾；多分类数据可以将疾病的严重程度进一步加以区分。此外，影像学数据、时间序列数据、树状数据也常被用于描述健康医疗相关的个体化特征。尽管新的数据类型含有更多信息，有助于gwas进一步挖掘出重要的遗传位点，但有效的数据分析方法仍然较少。许多研究者先将复杂的表型数据转化为连续性指标或二值型数据后再使用传统方法进行分析。这种处理方式会丢失部分信息，降低分析的效力，无法找出真正具有功能性的遗传位点。此外，全基因组测序数据正逐步替代微阵列基因分型数据用于检测高通量遗传信息。基因分型数据只能检测到预设的遗传位点，一般为常见变异位点，而测序数据可以检测到罕见变异位点，因此有潜力挖掘出更多的新发现。在当前的gwas研究中，表型数据和遗传数据都呈现出新的特点。根据这些新特点，通过设计并使用高效且准确的方法挖掘复杂性状中的有用信息，有助于发现致病基因、探究复杂疾病机制、实现精准医疗和个体化预防等领域的发展，具有重要的实际意义。2、虽然分析复杂性状很重要，但在大数据量的gwas中分析复杂性状在技术上具有挑战性，亟需快速、准确、普适的方法去分析。大数据量的含义既包括样本个体的数量大，也包括遗传位点的数量很多，这给关联分析方法带来了巨大的运算负担。常见的传统检验方法如wald检验和似然比检验需要对每一个待检验的位点拟合一次备择假设下的无约束模型，在全基因组范围内需要拟合上百万次模型，运算时间很长。对于这种情况score检验是有效的，因为score检验在一次全基因组范围的分析中只需要拟合一次原假设下的约束模型，大大提升了运算效率。为了找出与表型具有显著相关性的遗传位点，gwas研究的核心问题之一在于如何估计检验统计量在原假设为真时的理论分布。传统的score检验基于正态分布的假设计算p值，然而对于不平衡的表型分布和具有较低遗传变异率的位点，score统计量的正态分布假设不准确，会造成偏大的第一类错误率，进而导致统计推断结果的假阳性率过高。为了解决检验统计量的正态分布近似不准确的问题，鞍点近似方法(saddlepointapproximation,spa)在gwas中越来越受到关注。相较于仅利用了前两节矩的正态分布近似方法，鞍点近似方法使用了更多高阶矩信息，因而比传统的正态分布近似方法更准确，可以避免过高的假阳性率或假阴性率。随着鞍点近似方法理论研究的不断深入，已有多个基于鞍点近似方法的快速且准确的gwas分析算法被设计开发。3、群体分层或人群结构是gwas中的一个重要混杂因素，大型生物库可能包括大量来自异质人群的个体，例如混血人群中的个体。比如uk biobank中约40万人为白人族群，约10万人为其他族群，在白人族群中也存在着多个亚群，不同族群或亚群中遗传位点的遗传变异率以及表型的分布常有差异。混杂人群可能由于各种原因而产生，如迁徙、殖民化和来自不同种群的个体之间的通婚。将遗传位点衍生的遗传主成分(pc)作为协变量纳入分析所用的回归模型可以帮助控制人群混杂带来的影响，但仍可能导致偏高的假阳性率或假阴性率。例如，混杂人群中个体的祖先人群可能具有表型异质性(如连续型性状的方差异质性和生存表型的基线风险函数异质性)的情况。由于人群结构的带来的混杂，混杂人群中的个体通常被排除在分析之外。4、最近提出的spacox算法利用经验鞍点近似方法分析生存数据表型，通过分析ukbiobank的英国白种人数据，该算法找出了多个与复杂疾病相关的遗传位点。spacox算法理论上也可以作为一种通用的分析框架，用于分析其他的复杂性状。然而，spacox算法被设计开发用于分析同质人群，当spacox算法应用于混杂人群或多个族群的分析时，会出现大量的假阳性或假阴性结果，造成不可靠的分析结果。5、综上所述，目前迫切需要一种快速、准确、普遍适用于复杂性状且可以适用于混杂人群的全基因组关联分析方法，以避免混杂人群的分析结果中出现过高的假阳性率或假阴性率。技术实现思路1、基于背景技术存在的技术问题与缺陷，本发明提出了一种针对混杂人群的通用全基因组关联分析方法，适用于混杂人群的分析，可以控制第一类错误率不超过给定的显著性水平，且具有良好的统计功效，避免分析结果中出现过高的假阳性率或假阴性率。2、本发明提出的一种针对混杂人群的通用全基因组关联分析方法，包括如下步骤：3、s1：获取表型数据、基因型数据、混杂因素(年龄、性别以及遗传主成分等)数据；4、s2：基于样本数据构建回归模型，拟合位点遗传效应为0的约束模型并计算残差；5、s3：对于待检验的位点，基于遗传主成分估计个体水平的位点遗传变异率；6、s4：对于待检验的位点，计算用于检验位点遗传效应显著性的score检验统计量，基于个体水平的位点遗传变异率，使用正态分布近似方法与鞍点近似方法相结合的混合检验策略计算统计p值，以实现全基因组范围的关联分析。7、进一步地，所述基于样本数据构建回归模型，拟合位点遗传效应为0的约束模型并计算残差，具体包括：8、假定参与研究的个体样本数量为n，n个个体可以来自不同的族群或亚群，个体之间不具有亲缘相关性，对于个体i,1≤i≤n，令xi表示一个k维的混杂因素向量，gi表示位点的基因型，取值0,1或2，yi表示某一表型，对于不同的复杂性状，采用相应的回归模型进行分析，所用回归模型中的线性预测项为其中，βx表示混杂因素xi的k维系数向量，系数βg表示位点基因型gi的遗传效应；9、为了检验位点遗传效应的显著性，在原假设h0:βg＝0，即位点的遗传效应为0的约束条件下拟合模型，获得模型参数在约束条件下的极大似然估计，并计算约束模型下的n维残差向量r＝(r1,…,rn)t。10、进一步地，所述对于待检验的位点，基于遗传主成分估计个体水平的位点遗传变异率，具体包括：11、对于待检验的位点，假定n个个体的个体水平遗传变异率向量为q＝(q1,…,qn)t，其中，qi表示个体i在该位点的遗传变异率；12、令g＝(g1,…,gn)t表示位点的基因型向量，令表示包含元素全部为1的列向量和包含全部人群结构信息的前d个遗传主成分向量的n×(d+1)维遗传主成分矩阵；13、基于遗传主成分，可以通过拟合响应变量和解释变量分别为g和遗传主成分xpc的线性回归模型计算得到q的一个估计并令qc＝0，可以得到q的一个估计14、其中，15、基于遗传主成分，可以通过拟合响应变量和解释变量分别为和遗传主成分xpc的逻辑回归模型来得到q的一个估计其中，函数σ(x)＝1/(1+exp(-x))，表示该逻辑回归模型中的线性预测项的极大似然估计；16、为了兼顾估计的准确性与运算效率，若中的元素在区间[0,1]中的比例超过事先给定的介于0和1之间的正数(默认值0.9)，则使用作为q的估计；若中的元素在区间[0,1]中的比例没有超过事先给定的介于0和1之间的正数(默认值0.9)，则使用作为q的估计。17、进一步地，所述对于待检验的位点，计算用于检验位点遗传效应显著性的score检验统计量，基于个体水平的位点遗传变异率，使用正态分布近似方法与鞍点近似方法相结合的混合检验策略计算统计p值，以实现全基因组范围的关联分析，具体包括：18、对于待检验的位点，检验位点遗传效应显著性即检验原假设h0:βg＝0是否成立；19、假定该位点的样本基因型向量为g＝(g1,…,gn)t，拟合βg＝0的约束模型得到的残差向量为r＝(r1,…,rn)t，则用于检验位点遗传效应显著性的score检验统计量为20、在哈迪-温伯格平衡的假设下，将每个个体的位点基因型gi,i≤n视作相互独立的服从二项分布binom(2,qi)的随机变量；21、对于一次全基因组范围的分析，对于每个待检验的位点，使用相同的残差r的观测值来构造score检验统计量，基于个体水平的位点遗传变异率，在给定残差r的条件下，通过近似score检验统计量s在原假设h0:βg＝0为真时的条件分布来计算统计p值；22、(a)使用正态分布近似方法计算统计p值时，具体如下：23、在给定残差r的条件下，score检验统计量s在原假设h0:βg＝0为真时的期望和方差分别为：24、25、26、正态分布近似方法使用期望为方差为的正态分布对原假设为真时score检验统计量s的条件分布进行近似估计，其中，和分别表示原假设为真时score检验统计量s在给定r的条件下的期望和方差的估计；27、使用正态分布近似方法计算得到的双侧p值为：其中，s为s的观测值，φ(.)表示标准正态分布的累积分布函数；28、(b)使用鞍点近似方法计算统计p值时，具体如下：29、在给定残差r的条件下，为构造score检验统计量s在原假设h0:βg＝0为真时的累积量生成函数，首先估计位点基因型gi,1≤i≤n的矩生成函数：30、31、的一阶导数和二阶导数分别为：32、33、34、位点基因型gi的累积量生成函数的估计为其一阶导数和二阶导数分别为：35、36、37、因此，在给定残差r的条件下，score检验统计量s在原假设h0:βg＝0为真时的累积量生成函数的估计为：38、39、其一阶导数和二阶导数分别为：40、41、42、对于给定的实数s0和残差r，计算使得等式成立的并计算得到以及则根据barndorff-nielsen鞍点近似公式，在原假设h0:βg＝0为真时，概率pr(s