CN114898809A - 适用复杂性状的基因-环境交互的分析方法及存储介质 - Google Patents

适用复杂性状的基因-环境交互的分析方法及存储介质 Download PDF

Info

Publication number
CN114898809A
CN114898809A CN202210373636.5A CN202210373636A CN114898809A CN 114898809 A CN114898809 A CN 114898809A CN 202210373636 A CN202210373636 A CN 202210373636A CN 114898809 A CN114898809 A CN 114898809A
Authority
CN
China
Prior art keywords
marginal
gene
environment interaction
normal distribution
approximation method
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202210373636.5A
Other languages
English (en)
Other versions
CN114898809B (zh
Inventor
张纪峰
毕文健
马雨茁
赵延龙
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Academy of Mathematics and Systems Science of CAS
Original Assignee
Academy of Mathematics and Systems Science of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Academy of Mathematics and Systems Science of CAS filed Critical Academy of Mathematics and Systems Science of CAS
Priority to CN202210373636.5A priority Critical patent/CN114898809B/zh
Publication of CN114898809A publication Critical patent/CN114898809A/zh
Application granted granted Critical
Publication of CN114898809B publication Critical patent/CN114898809B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/20—Allele or variant detection, e.g. single nucleotide polymorphism [SNP] detection
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Medical Informatics (AREA)
  • Biophysics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • General Health & Medical Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Biotechnology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Bioethics (AREA)
  • Public Health (AREA)
  • Evolutionary Computation (AREA)
  • Epidemiology (AREA)
  • Databases & Information Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Software Systems (AREA)
  • Artificial Intelligence (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Genetics & Genomics (AREA)
  • Molecular Biology (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Complex Calculations (AREA)

Abstract

本发明公开了一种适用复杂性状的基因‑环境交互的分析方法及存储介质，方法包括：获取表型数据、基因型数据、环境因素数据和混杂因素数据；选择分析所需的广义线性回归模型，在边际遗传效应与边际基因‑环境交互作用均为0的假设下拟合约束模型、估计模型参数并计算残差；对于待检验的位点，估计其次要等位基因频率，基于score统计量检验其边际遗传效应的显著性；选择并计算相应的用于检验边际基因‑环境交互作用显著性的检验统计量；使用正态分布近似方法与鞍点近似方法相结合的混合检验策略计算统计p值以检验边际基因‑环境交互作用的显著性，实现全基因组范围的关联分析。本发明具有适用范围广、分析速度快、准确度高等优点。

Description

适用复杂性状的基因-环境交互的分析方法及存储介质
技术领域
本发明涉及***生物学技术领域，涉及一种生物医学大数据处理方法，特别涉及一种适用复杂性状的基因-环境交互的分析方法。
背景技术
随着高通量基因测序技术、高精度成像技术以及电子健康记录***的发展，生物医学进入了健康医疗信息化的大数据新时代。海量的生物医学大数据为***生物学研究提供了丰富的研究资源。比如英国生物样本库(UK Biobank)收集了50万英国人的数据并向全球科研人员开放，其中包含基因测序、磁共振扫描成像、临床指标、生活方式等多维度、跨尺度数据。我国的大型生物样本库近年来也得到了快速发展，中国慢性病前瞻性研究、张江国际脑库等数据库收集了数以万计甚至数以十万计的样本信息并进行了全面的数据采集。利用快速、有效的算法从大数据中挖掘出可理解的模式，对于推进生物医学的基础研究和临床研究，特别是对于精准医疗、个体化预防、复杂疾病的智能诊疗等领域具有重要的研究意义。
生物医学研究中大量的复杂性状(如身高、体重等)和疾病(如糖尿病、高血压等)等表型呈现出明显的家族聚集性，比如患者直系后代的发病风险常高于普通人，这些现象表明遗传和环境因素在个体发育和患病过程中起到了重要作用。全基因组关联分析(genome-wide association study,GWAS)是一种通过检验全基因组遗传标记与表型变异关联的显著性来定位与性状相关的遗传位点,在群体水平上解析性状遗传基础的方法。简单的讲,就是从人类全基因组范围内的序列变异(单核苷酸多态,Single NucleotidePolymorphism,SNP),筛选出与疾病性状关联的SNP位点。全基因组关联分析可以在全基因水平上对复杂性状的遗传变异进行关联分析。
复杂性状的产生与发展不能完全由遗传变异来解释，而是遗传变异和环境因素共同作用的结果。基因-环境交互作用分析可以捕捉基因与环境因素之间的交互作用，对复杂表型具有重要的影响，识别基因-环境交互作用对于个性化和分层预防与治疗尤为重要，基于关联分析找出基因-环境交互作用(G×Eeffect)有助于我们了解疾病的发病机制和设计个体化治疗方案，值得深入研究。但是基因-环境交互作用的关联性较弱，因此常常需要较大的数据样本量才能达到足够的检验效能，这大大增加了发现关联性的难度。
最近几年，随着测序技术和电子健康记录(Electronic Health Record,EHR)的发展，很多大型的生物样本资源库(Biobank)为研究者提供了非常详尽的大样本信息，这使得我们可以在全基因组(genome-wide)乃至于全表型组(phenome-wide)尺度下进行基因-环境交互作用(G×E)的关联分析。然而，样本量的大量增长对全基因组关联分析算法的运算速度提出了更严苛的要求。比如，当分析40万样本数据时，拟合几百万次回归模型需要数月的运算时间。此外，表型具有不平衡的分布也可能会使得现有分析的方法会失效。
现有的大多数基因-环境交互作用的有效分析方法都侧重于数量性状或者二元性状，目前还没有一种统一的快速准确的统计方法可以应用于大型生物样本库中各种不同复杂性状(例如生存数据表型、多分类表型)的全基因组范围的基因-环境交互作用研究。现有的方法都存在一定的缺陷，比如：
(1)似然比检验方法
对于每一个待检验的位点，似然比检验需要拟合一次备择假设下的无约束模型和一次零假设下的约束模型，对于全基因组分析需要进行上百万次的模型拟合，运算时间长，需要几个月时间。
(2)Score检验方法
由于分析基因-环境交互作用所采用的约束模型仍然包含边际遗传效应，因此，对于每一个待检验的位点，Score检验仍然需要拟合一次约束模型，与似然比检验类似，对于全基因组分析需要进行上百万次的模型拟合，运算时间很长。
(3)两步骤方法
两步骤方法根据边际遗传效应的关联筛选出位点，然而，由于此方法在筛选步骤中排除掉了大多数的位点，因此可能会遗漏掉部分具有基因-环境交互作用的位点，并且无法生成全基因组范围的汇总统计数据，这些汇总统计数据对于全表型组关联分析和荟萃分析具有重要的作用。
(4)SPAGE方法
SPAGE方法是一种针对大样本量下的二元性状的全基因组基因-环境交互作用分析的快速准确的算法，该算法利用条件期望构造了一种新的Score统计量并利用鞍点近似方法(saddlepoint approximation,SPA)来计算统计p值。对于一次全基因组范围的基因-环境交互作用分析，SPAGE方法只需在零假设下进行一次逻辑回归模型，这极大地提升了运算速度。对于高度的病例-对照不平衡(病例对照比为1:100)情形下针对罕见变异(MAF＝0.001)的分析，SPAGE方法仍然可以非常好地控制第一类错误率，同时具有相当高的统计效力。然而，SPAGE方法仅适用于二元性状的分析，并不适用于针对其他复杂性状的基因-环境交互作用分析。
综上，目前迫切需要一种适用复杂性状的基因-环境交互的分析方法来避免大量的假阳性或假阴性结果。
发明内容
本发明的目的在于克服上述现有技术存在的缺陷，以提供一种准确的适用复杂性状的基因-环境交互的分析方法及存储介质。
本发明的目的可以通过以下技术方案来实现：
一种适用复杂性状的基因-环境交互的分析方法，包括以下步骤：
获取分析所需的表型数据、基因型数据、环境因素数据和混杂因素数据；
基于表型数据、基因型数据、环境因素数据和混杂因素数据，选择分析所需的广义线性回归模型，在边际遗传效应与边际基因-环境交互作用均为0的假设下拟合约束模型、估计模型参数并计算残差；
对于待检验的位点，估计其次要等位基因频率(MAF)，基于score统计量检验其边际遗传效应的显著性；
根据边际遗传效应的显著性来选择并计算相应的用于检验边际基因-环境交互作用显著性的检验统计量；
对于检验边际基因-环境交互作用的检验统计量，使用正态分布近似方法与鞍点近似方法相结合的混合检验策略计算统计p值以检验边际基因-环境交互作用的显著性，从而实现全基因组范围的关联分析。
进一步地，所述基于表型数据、基因型数据、环境因素数据和混杂因素数据，选择分析所需的广义线性回归模型，在边际遗传效应与边际基因-环境交互作用均为0的假设下拟合约束模型、估计模型参数并计算残差，包括：
假设参与研究的个体样本数量为n，对于个体i,1≤i≤n,令Xi表示一个k维的混杂因素向量(包括年龄、性别等混杂因素)，Gi表示基因型(Gi＝0,1,2)，Ei表示环境因素，Yi表示某一表型，线性预测项
Figure BDA0003589875320000051
Figure BDA0003589875320000052
其中βX表示混杂因素Xi的k维系数向量，βE、βG和βG×E分别表示环境因素Ei、基因型Gi和基因-环境交互项GiEi所对应的系数，假设表型Yi的期望可以表示成ηi的函数。此假设对于GWAS中大部分的模型是成立的，对数似然函数l(βX,βE,βG,βG×E)＝lnL(βX,βE,βG,βG×E)是ηi的函数，可以表示为
Figure BDA0003589875320000053
的形式，其中f(.)表示某个一元函数，对于不同的复杂性状，采用与表型相对应的广义线性回归模型进行分析；
为了检验边际基因-环境交互作用的显著性，在假设Hc:βG＝βG×E＝0，即边际遗传效应和边际基因-环境交互作用均为0的假设下拟合模型，得到参数(βX,βE)在假设Hc:βG＝βG×E＝0下的极大似然估计
Figure BDA0003589875320000061
并定义残差向量
Figure BDA0003589875320000062
其中，
Figure BDA0003589875320000063
f'(.)表示f(.)的一阶导数。
进一步地，所述对于待检验的位点，估计其次要等位基因频率(MAF)，基于score统计量检验其边际遗传效应的显著性，包括：
假设待检验位点的次要等位基因频率(MAF)为q，令
Figure BDA0003589875320000064
Figure BDA0003589875320000065
作为q的估计；
定义基因型向量G和基因-环境交互向量GE分别为：
Figure BDA0003589875320000066
计算用于检验边际遗传效应的score检验统计量
Figure BDA0003589875320000067
并认为检验统计量
Figure BDA0003589875320000068
在假设HG:βG＝0成立时服从期望为0、方差为
Figure BDA0003589875320000069
的正态分布，其中
Figure BDA00035898753200000610
表示
Figure BDA00035898753200000611
的方差
Figure BDA00035898753200000612
的估计，令
Figure BDA00035898753200000613
作为使用正态分布近似方法检验假设HG:βG＝0得到的双侧p值，其中
Figure BDA00035898753200000614
为
Figure BDA00035898753200000615
的观测值，Φ(.)表示标准正态分布的累积分布函数。
进一步地，所述根据边际遗传效应的显著性来选择并计算相应的用于检验边际基因-环境交互作用显著性的检验统计量，包括：
检验边际基因-环境交互作用的显著性即检验零假设H0:βG×E＝0是否成立；
令∈为某一给定的介于0和1之间的正数(如∈＝0.001)；
若对假设HG:βG＝0进行检验得到的p值小于或等于∈，则使用
Figure BDA0003589875320000071
作为检验零假设H0:βG×E＝0的检验统计量，其中
Figure BDA0003589875320000072
若对假设HG:βG＝0进行检验得到的p值大于∈，则使用
Figure BDA0003589875320000073
Figure BDA0003589875320000074
作为检验零假设H0:βG×E＝0的检验统计量，其中
Figure BDA0003589875320000075
Figure BDA0003589875320000076
为适应于基因型的残差向量，In为n阶单位矩阵，W＝(1n,G)为n×2矩阵，1n＝(1,1,…,1)T为元素全部为1的n维列向量。
进一步地，所述正态分布近似方法，包括：
若检验零假设H0:βG×E＝0的统计量为SG×E，则正态分布近似方法使用期望为0，方差为
Figure BDA0003589875320000077
的正态分布对零假设H0:βG×E＝0成立时统计量SG×E的分布进行近似估计，计算得到统计p值为：
Figure BDA0003589875320000078
其中sG×E为SG×E的观测值，Φ(.)是标准正态分布的累积分布函数；
若检验零假设H0:βG×E＝0的统计量为
Figure BDA0003589875320000081
则正态分布近似方法使用期望为0，方差为
Figure BDA0003589875320000082
的正态分布对零假设H0:βG×E＝0成立时统计量
Figure BDA0003589875320000083
的分布进行近似估计，其中，
Figure BDA0003589875320000084
为
Figure BDA0003589875320000085
的方差
Figure BDA0003589875320000086
的估计，使用正态分布近似方法计算得到统计p值为：
Figure BDA0003589875320000087
其中
Figure BDA0003589875320000088
为
Figure BDA0003589875320000089
的观测值，Φ(.)是标准正态分布的累积分布函数。
进一步地，所述鞍点近似方法，包括：
若检验零假设H0:βG×E＝0的统计量为SG×E，则在零假设H0:βG×E＝0成立时将Gi,i≤n视作服从二项分布B(2,q)的独立同分布的随机变量序列，将Ri,i≤n和Ei,i≤n视作常系数；为构造SG×E在零假设H0成立时的累积量生成函数(CGF)，首先估计Gi的矩生成函数(MGF)为：
Figure BDA00035898753200000810
Figure BDA00035898753200000811
的一阶导数和二阶导数分别为：
Figure BDA00035898753200000812
Figure BDA00035898753200000813
Gi的累积量生成函数的估计为
Figure BDA00035898753200000814
其一阶导数和二阶导数分别为：
Figure BDA00035898753200000815
因此，SG×E在H0成立时的累积量生成函数的估计为：
Figure BDA00035898753200000816
其一阶导数和二阶导数分别为：
Figure BDA0003589875320000091
Figure BDA0003589875320000092
给定检验统计量SG×E的观测值sG×E，首先计算使得
Figure BDA0003589875320000093
的ζ，然后计算得到
Figure BDA0003589875320000094
和
Figure BDA0003589875320000095
为了近似估计统计量SG×E在零假设H0下的分布，根据Barndorff-Nielsen鞍点近似公式，使用
Figure BDA0003589875320000096
来近似估计概率Pr(SG×E<sG×E)，其中Φ(.)是标准正态分布的累积分布函数，因此，使用鞍点近似方法计算得到的双侧p值为：
Figure BDA0003589875320000097
其中，
Figure BDA0003589875320000098
和
Figure BDA0003589875320000099
分别表示由上述鞍点近似方法得到的Pr(SG×E<-|sG×E|)和Pr(SG×E<|sG×E|)的估计。
若检验零假设H0:βG×E＝0的统计量为
Figure BDA00035898753200000910
则在零假设H0:βG×E＝0成立时将Gi,i≤n视作服从二项分布B(2,q)的独立同分布的随机变量序列，将
Figure BDA00035898753200000911
和Ei,i≤n视作常系数；为构造
Figure BDA00035898753200000912
在零假设H0成立时的累积量生成函数(CGF)，首先估计Gi的矩生成函数(MGF)为：
Figure BDA00035898753200000913
Figure BDA00035898753200000914
的一阶导数和二阶导数分别为：
Figure BDA00035898753200000915
Figure BDA00035898753200000916
Gi的累积量生成函数的估计为
Figure BDA00035898753200000917
其一阶导数和二阶导数分别为：
Figure BDA00035898753200000918
因此，
Figure BDA0003589875320000101
在H0成立时的累积量生成函数的估计为：
Figure BDA0003589875320000102
其一阶导数和二阶导数分别为：
Figure BDA0003589875320000103
Figure BDA0003589875320000104
给定检验统计量
Figure BDA0003589875320000105
的观测值
Figure BDA0003589875320000106
首先计算使得
Figure BDA0003589875320000107
的
Figure BDA0003589875320000108
然后计算得到
Figure BDA0003589875320000109
和
Figure BDA00035898753200001010
为了近似估计统计量
Figure BDA00035898753200001011
在零假设H0下的分布，根据Barndorff-Nielsen鞍点近似公式，使用
Figure BDA00035898753200001012
来近似估计概率
Figure BDA00035898753200001013
其中Φ(.)是标准正态分布的累积分布函数，因此，使用鞍点近似方法计算得到的双侧p值为：
Figure BDA00035898753200001014
其中，
Figure BDA00035898753200001015
和
Figure BDA00035898753200001016
分别表示由上述鞍点近似方法得到的
Figure BDA00035898753200001017
和
Figure BDA00035898753200001018
的估计。
进一步地，所述对于检验边际基因-环境交互作用的检验统计量，使用正态分布近似方法与鞍点近似方法相结合的混合检验策略计算统计p值以检验边际基因-环境交互作用的显著性，从而实现全基因组范围的关联分析，包括：
令r为一个预先选定的正数(如r＝2)；
若检验零假设H0:βG×E＝0的统计量为SG×E，sG×E为SG×E的观测值，如果
Figure BDA00035898753200001019
则使用所述正态分布近似方法来计算统计p值；否则，如果
Figure BDA0003589875320000111
则使用所述鞍点近似方法来计算统计p值；
若检验零假设H0:βG×E＝0的统计量为
Figure BDA0003589875320000112
为
Figure BDA0003589875320000113
的观测值，如果
Figure BDA0003589875320000114
则使用所述正态分布近似方法来计算统计p值；否则，如果
Figure BDA0003589875320000115
则使用所述鞍点近似方法来计算统计p值；
获得统计p值后，在给定的显著性水平下，进行全基因组范围的关联分析。
又一方面，本发明还公开一种计算机可读存储介质，存储有计算机程序，所述计算机程序被处理器执行时，使得所述处理器执行如上所述适用复杂性状的基因-环境交互的分析方法。
再一方面，本发明还公开一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行如上所述适用复杂性状的基因-环境交互的分析。
与现有技术相比，本发明存在以下技术效果：本发明提出了一种适用复杂性状的基因-环境交互的分析方法，此方法在应用时拟合与基因型无关的约束模型，即在边际遗传效应与边际基因-环境交互作用均为0的假设下拟合模型并计算残差；接下来估计位点的次要等位基因频率(MAF)，并基于score统计量检验其边际遗传效应的显著性；根据检验边际遗传效应的显著性来选择不同的检验统计量以检验边际基因-环境交互作用的显著性，兼顾了准确性与运算效率；提出了一个新的视角，在检验边际基因-环境交互作用的显著性时，在零假设下将基因型视作随机变量，将残差视作固定参数；此外，在计算统计p值时使用一种混合检验策略，将正态分布近似方法和鞍点近似方法相结合，这使得本发明在大规模的全基因组关联分析中快速且准确。
附图说明
图1是一种适用复杂性状的基因-环境交互的分析方法的流程图；
图2正态分布近似方法与鞍点近似方法相结合的混合检验策略的示意图；
图3是生存数据表型分析中针对位点的边际基因-环境交互作用和边际遗传效应均为0的情形的第一类错误率数值模拟结果图，其中横轴代表次要等位基因频率，纵轴代表经验第一类错误率，ER:Event Rate表示事件发生率，MAF:Minor Allele Frequency表示次要等位基因频率，EmpGxE表示本发明所提出的正态分布近似与鞍点近似相结合的混合检验方法，Normal表示正态分布近似方法；
图4是生存数据表型分析中针对位点的边际基因-环境交互作用和边际遗传效应均为0的情形的数值模拟的统计p值分位图，其中横轴代表期望的-log10(p值)，纵轴代表观测的-log10(p值)，ER:Event Rate表示事件发生率，MAF:Minor Allele Frequency表示次要等位基因频率，EmpGxE表示本发明所提出的正态分布近似与鞍点近似相结合的混合检验方法，Normal表示正态分布近似方法；
图5是生存数据表型分析中针对位点的边际基因-环境交互作用为0、边际遗传效应不为0的情形的数值模拟的统计p值分位图，其中横轴代表期望的-log10(p值)，纵轴代表观测的-log10(p值)，ER:Event Rate表示事件发生率，EmpGxE表示本发明所提出的正态分布近似与鞍点近似相结合的混合检验方法，Normal表示正态分布近似方法；
图6是生存数据表型分析中的统计效力的数值模拟结果图，其中横轴代表边际基因-环境交互作用，纵轴代表经验统计效力，ER:Event Rate表示事件发生率，MAF:MinorAllele Frequency表示次要等位基因频率，EmpGxE表示本发明所提出的正态分布近似与鞍点近似相结合的混合检验方法，Normal表示正态分布近似方法。
具体实施方式
为了进一步说明本发明的特征，下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。所附图仅供参考与说明之用，并非用来对本发明的保护范围加以限制。
如图1所示，本实施例公开了一种适用复杂性状的基因-环境交互的分析方法，并将其应用于生存数据表型的基因-环境交互作用分析，包括如下步骤S1至S5：
S1、获取分析所需的表型数据、基因型数据、环境因素数据和混杂因素数据；
S2、基于表型数据、基因型数据、环境因素数据和混杂因素数据，选择分析所需的广义线性回归模型，在边际遗传效应与边际基因-环境交互作用均为0的假设下拟合约束模型、估计模型参数并计算残差；
S3、对于待检验的位点，估计其次要等位基因频率(MAF)，基于score统计量检验其边际遗传效应的显著性；
S4、根据边际遗传效应的显著性来选择并计算相应的用于检验边际基因-环境交互作用显著性的检验统计量；
S5、对于检验边际基因-环境交互作用的检验统计量，使用正态分布近似方法与鞍点近似方法相结合的混合检验策略计算统计p值以检验边际基因-环境交互作用的显著性，从而实现全基因组范围的关联分析。
具体来说，在上述步骤S1中，假设参与研究的个体样本数量为n，对于个体i,1≤i≤n,令Xi表示一个k维的混杂因素向量(包括年龄、性别等混杂因素)，Gi表示基因型(Gi＝0,1,2)，Ei表示环境因素，Yi表示某一表型。
具体来说，在上述步骤S2中，假设表型Yi的期望可以表示成线性预测项
Figure BDA0003589875320000141
的函数,其中βX表示混杂因素Xi的k维系数向量，βE、βG和βG×E分别表示环境因素Ei、基因型Gi和基因-环境交互项GiEi所对应的系数，此假设对于GWAS中大部分的模型是成立的。对数似然函数l(βX,βE,βG,βG×E)＝lnL(βX,βE,βG,βG×E)是ηi的函数，可以表示为
Figure BDA0003589875320000142
的形式，其中f(.)表示某个一元函数；对于不同的复杂性状，采用与表型相对应的广义线性回归模型进行分析。关于边际遗传效应βG和边际基因-环境交互作用βG×E的得分函数分别为l(βX,βE,βG,βG×E)关于βG和关于βG×E的偏导数，具有以下形式：
Figure BDA0003589875320000151
Figure BDA0003589875320000152
其中，f'(.)表示f(.)的一阶导数。
为了检验边际基因-环境交互作用的显著性，传统的score检验在零假设H0:βG×E＝0下拟合模型并计算参数(βX,βE,βG)的极大似然估计
Figure BDA0003589875320000153
以此来计算用于检验边际基因-环境交互作用的score统计量
Figure BDA0003589875320000154
其中
Figure BDA0003589875320000155
Figure BDA0003589875320000156
然而此策略需要对每个位点拟合一次模型，对于全基因组范围的分析，其计算成本非常高昂。
为了提高计算效率，本发明并非在零假设H0:βG×E＝0下拟合模型，而在假设Hc:βG＝βG×E＝0，即在边际遗传效应和边际基因-环境交互作用均为0的假设下拟合模型，进而得到参数(βX,βE在假设Hc:βG＝βG×E＝0下的极大似然估计
Figure BDA0003589875320000157
并定义残差向量
Figure BDA0003589875320000158
其中
Figure BDA0003589875320000159
下面两个例子展示了在生存数据表型和多分类表型的基因-环境交互作用分析中基于表型数据、基因型数据、环境因素数据和混杂因素数据，选择分析所需的广义线性回归模型，在边际遗传效应与边际基因-环境交互作用均为0的假设下拟合模型，估计模型参数并计算得到残差的过程：
(1)生存数据表型：
假设表型是生存数据，参与研究的个体样本数量为n，使用Cox比例风险回归模型进行分析：
Figure BDA0003589875320000161
其中，λ(t；Xi)是风险函数，λ0(t)是基线风险函数。令Ci表示事件的删失时间，
Figure BDA0003589875320000162
表示事件的发生时间，
Figure BDA0003589875320000163
表示观测到的生存数据，
Figure BDA0003589875320000164
用于显示事件是否发生，其中I(.)表示示性函数。函数f(ηi)和f'(ηi)分别为
Figure BDA0003589875320000165
f'(ηi)＝δi-Λi
其中，
Figure BDA0003589875320000166
表示在时刻ti的风险集，
Figure BDA0003589875320000167
表示个体i直至时刻ti的累计风险。
为了检验边际基因-环境交互作用，传统的score检验在因-环境交互作用为0的零假设H0:βG×E＝0下拟合模型，得到参数(βX,βE,βG)的极大似然估计
Figure BDA0003589875320000168
并计算
Figure BDA0003589875320000169
检验零假设H0:βG×E＝0的score统计量为
Figure BDA0003589875320000171
Figure BDA0003589875320000172
其中
Figure BDA0003589875320000173
表示Λi在零假设H0:βG×E＝0下拟合模型后得到的估计。
为了提高计算效率，不在零假设H0:βG×E＝0下拟合模型，而在边际遗传效应与边际基因-环境交互作用均为0的假设Hc:βG＝βG×E＝0下拟合模型，得到参数(βX,βE)的极大似然估计
Figure BDA0003589875320000174
并计算在假设Hc下拟合模型后得到的线性预测项
Figure BDA0003589875320000175
以及Λi的估计
Figure BDA0003589875320000176
令
Figure BDA0003589875320000177
表示在假设Hc:βG＝βG×E＝0下拟合模型后计算得到的鞅残差，R＝(R1,…,Rn)T表示相应的残差向量。
(2)多分类数据表型：
假设表型是多分类数据，参与研究的个体样本数量为n，令J表示类别数量，令Yi＝1,2,…,J表示个体i的多分类表型，使用以下比例优势逻辑回归模型进行分析：
Figure BDA0003589875320000178
其中，vij＝Pr(Yi≤j|Xi,Ei,Gi)表示在给定Xi、Ei和Gi的条件下，Yi≤j的累计概率，分割点ε:ε1<…<εJ＝∞用于对数据进行分类。
对于个体i,1≤i≤n，定义一个J×1向量
Figure BDA0003589875320000179
作为多分类表型Yi的一个等价的表达形式；如果Yi＝j，那么Yij＝1并且
Figure BDA00035898753200001710
中其余的元素均为0。令μij表示Yij的期望，即μij＝E(Yij)＝Pr(Yij＝1)＝Pr(Yi＝j)＝Pr(Yi≤j)-Pr(Yi≤j-1)，其中E(.)表示某随机变量的期望；令
Figure BDA0003589875320000181
Figure BDA0003589875320000182
函数f(ηi)和f'(ηi)分别为
Figure BDA0003589875320000183
Figure BDA0003589875320000184
为了检验边际基因-环境交互作用，传统的score检验在边际基因-环境交互作用为0的零假设H0:βG×E＝0下拟合模型，得到参数βX,βE,βG,εj,1≤j≤J的极大似然估计
Figure BDA0003589875320000185
1≤j≤J，并计算
Figure BDA0003589875320000186
接下来分别计算得到vij、μij和Vij在H0:βG×E＝0下拟合模型后得到的估计量
Figure BDA0003589875320000187
Figure BDA0003589875320000188
以及
Figure BDA0003589875320000189
检验零假设H0:βG×E＝0的Score统计量为
Figure BDA00035898753200001810
Figure BDA00035898753200001811
为了提高计算效率，不在零假设H0:βG×E＝0下拟合模型，而在边际遗传效应与边际基因-环境交互作用均为0的假设Hc:βG＝βG×E＝0下拟合模型，得到参数βX,βE,εj,1≤j≤J的极大似然估计
Figure BDA00035898753200001812
1≤j≤J，并计算在假设Hc:βG＝βG×E＝0下拟合模型后得到的线性预测项
Figure BDA00035898753200001813
以及vij、μij和Vij的估计量
Figure BDA00035898753200001814
Figure BDA0003589875320000191
和
Figure BDA0003589875320000192
Figure BDA0003589875320000193
定义
Figure BDA0003589875320000194
为在假设Hc:βG＝βG×E＝0下拟合模型后计算得到的残差，令R＝(R1,…,Rn)T表示相应的残差向量。
具体来说，在上述步骤S3中，假设待检验位点的次要等位基因频率(MAF)为q，令
Figure BDA0003589875320000195
作为q的估计，
Figure BDA0003589875320000196
是准确且易于计算的估计量。
定义基因型向量G和基因-环境交互向量GE分别为：
Figure BDA0003589875320000197
计算用于检验边际遗传效应的score检验统计量
Figure BDA0003589875320000198
对于
Figure BDA0003589875320000199
在零假设H0:βG×E＝0成立时将残差Ri,i≤n视作常系数，将Gi,i≤n视作服从二项分布B(2,q)的独立同分布的随机变量序列，并认为检验统计量
Figure BDA00035898753200001910
在假设HG:βG＝0成立时服从期望为0、方差为
Figure BDA00035898753200001911
的正态分布，其中
Figure BDA00035898753200001912
表示
Figure BDA00035898753200001913
的方差
Figure BDA00035898753200001914
的估计，令
Figure BDA00035898753200001915
作为使用正态分布近似方法检验假设HG:βG＝0得到的双侧p值，其中
Figure BDA00035898753200001916
表示
Figure BDA00035898753200001917
的观测值，Φ(.)表示标准正态分布的累积分布函数。
具体来说，在上述步骤S4中，为了在假设Hc:βG＝βG×E＝0下拟合模型并对零假设H0:βG×E＝0进行检验，一种自然的方式是使用
Figure BDA0003589875320000201
作为检验统计量。然而，对于边际遗传效应βG≠0的情形，此策略是不可行的，因为它忽略了边际遗传效应对表型的作用。为了适应于边际遗传效应，需要提出新的检验统计量来对零假设H0:βG×E＝0进行检验。下面根据在步骤S3中对假设HG:βG＝0进行检验得到的p值，来选择并计算不同的检验统计量对零假设H0进行检验。
令∈为某一给定的介于0和1之间的正数(如∈＝0.001)。
若在步骤S3中对假设HG:βG＝0进行检验得到的统计p值小于或等于∈，则为了适应于边际遗传效应，提出一个新的检验统计量：
Figure BDA0003589875320000202
其中，
Figure BDA0003589875320000203
对于SG×E提出一个新的视角，在零假设H0:βG×E＝0成立时将残差Ri,i≤n、环境因素Ei,i≤n以及λ视作常系数，将Gi,i≤n视作服从二项分布B(2,q)的独立同分布的随机变量序列，其中q是次要等位基因频率(MAF)。这个假设在Hc:βG＝βG×E＝0成立时是合理的，因为在实际分析中获取某次表型向量的观测值后，会用全基因组的位点进行与所获取表型观测值之间的关联分析。在零假设H0:βG×E＝0成立时，将Ri,i≤n视作常系数是一种近似的方法，当边际遗传效应在0附近时是合理的。由于Ri,i≤n满足线性约束条件
Figure BDA0003589875320000211
和
Figure BDA0003589875320000212
因此在假设Hc:βG＝βG×E＝0成立时，
Figure BDA0003589875320000213
和
Figure BDA0003589875320000214
的期望和方差分别为：
Figure BDA0003589875320000215
Figure BDA0003589875320000216
在假设Hc:βG＝βG×E＝0成立时，S1和S2之间的协方差为：
Figure BDA0003589875320000217
其中，Ec(.)、Varc(.)和Covc(.,.)分别表示某统计量在假设Hc成立时的期望、方差以及某两个统计量之间的协方差。因此，λ可以表示为
Figure BDA0003589875320000218
注意到
Figure BDA0003589875320000219
包含了边际遗传效应的信息，若
Figure BDA00035898753200002110
其中，
Figure BDA00035898753200002111
表示S1和S2之间的协方差，
Figure BDA00035898753200002112
表示
Figure BDA00035898753200002113
的方差，则SG×E与
Figure BDA00035898753200002114
是不相关的，因为：
Figure BDA00035898753200002115
这表明在假设Hc成立时，通过减去
Figure BDA00035898753200002116
边际遗传效应的影响被从
Figure BDA00035898753200002117
中去除了。在边际遗传效应在0附近的情形
Figure BDA00035898753200002118
Figure BDA00035898753200002119
因此，在此情形下可以近似认为检验统计量SG×E仅包含了边际基因-环境交互作用的信息；对于边际遗传效应远离0的情形，λ与
Figure BDA0003589875320000221
相差较大，检验统计量SG×E仍然包含了较多的边际遗传效应的信息，在此情形下使用SG×E检验边际基因-环境交互作用会得到不可靠的统计推断。
若在步骤S2中对假设HG:βG＝0进行检验得到的统计p值大于∈，则使用
Figure BDA0003589875320000222
作为检验零假设H0:βG×E＝0的检验统计量，其中
Figure BDA0003589875320000223
为适应于基因型的残差向量，In为n阶单位矩阵，W＝(1n,G)为n×2矩阵，1n＝(1,1,…,1)T为元素全部为1的n维列向量。
Figure BDA0003589875320000224
可以被解释为残差向量R对基因型向量G进行线性回归并使用最小二乘法估计模型参数后得到的新的残差向量。因此，在统计量
Figure BDA0003589875320000225
中，边际遗传效应已通过对R的变换(In-W(WTW)-1WT)R而去除，
Figure BDA0003589875320000226
是仅包含边际基因-环境交互作用的信息的检验统计量。
注意到，计算统计量SG×E比计算统计量
Figure BDA0003589875320000227
的速度更快，因为在计算
Figure BDA0003589875320000228
的过程中需要更多的时间用于计算
Figure BDA0003589875320000229
但是相较于SG×E，使用统计量
Figure BDA00035898753200002210
可以获得更加准确的统计推断。因此，上述根据检验边际遗传效应得到的统计p值来选择并计算相应的用于检验边际基因-环境交互作用的检验统计量的策略兼顾了运算效率与准确性，在两者间达到了平衡。
具体来说，分别介绍在上述步骤S5中提及的正态分布近似方法、鞍点近似方法以及正态分布近似方法与鞍点近似方法相结合的混合检验策略。
(1)正态分布近似方法
若检验零假设H0:βG×E＝0的统计量为SG×E，则正态分布近似方法使用期望为0，方差为
Figure BDA0003589875320000231
的正态分布对零假设H0:βG×E＝0成立时统计量SG×E的分布进行近似估计，其中，
Figure BDA0003589875320000232
为SG×E的方差Var(SG×E)的估计，使用正态分布近似方法计算得到的统计p值为：
Figure BDA0003589875320000233
其中sG×E为SG×E的观测值，Φ(.)是标准正态分布的累积分布函数。
若检验零假设H0:βG×E＝0的统计量为
Figure BDA0003589875320000234
则正态分布近似方法使用期望为0，方差为
Figure BDA0003589875320000235
的正态分布对零假设H0:βG×E＝0成立时统计量
Figure BDA0003589875320000236
的分布进行近似估计，其中，
Figure BDA0003589875320000237
为
Figure BDA0003589875320000238
的方差
Figure BDA0003589875320000239
的估计，使用正态分布近似方法计算得到统计p值为：
Figure BDA00035898753200002310
其中
Figure BDA00035898753200002311
为
Figure BDA00035898753200002312
的观测值，Φ(.)是标准正态分布的累积分布函数。
(2)鞍点近似方法
在边际遗传效应在0附近，并且统计量的观测值在其分布期望附近的情况下，正态分布近似方法的效果很好；然而，在表型分布不平衡的情形，尤其是当所检测位点的次要等位基因频率(MAF)很低时，正态分布近似方法对尾部分布的近似并不准确，此时使用正太分布近似方法不能准确地控制第一类错误率，会造成大量假判真。
相较于仅仅使用了前两阶矩(期望和方差)的信息的正态分布近似方法，鞍点近似方法(SPA)利用累积量生成函数(CGF)中包含的全部高阶矩的信息来近似统计量的分布，因此鞍点近似通常非常准确，在近似尾部概率分布时仍然保持准确性；并且，鞍点近似不局限于某一特定的分布形状，如正态分布所限制的形状。总而言之，鞍点近似方法使用了更多关于分布的信息，因此近似估计的精度更高。为了在检验边际基因-环境交互作用的显著性时得到准确的统计推断，此处提出新的经验鞍点近似算法：
若检验零假设H0:βG×E＝0的统计量为SG×E，则在零假设H0:βG×E＝0成立时将Gi,i≤n视作服从二项分布B(2,q)的独立同分布的随机变量序列，将Ri,i≤n和Ei,i≤n视作常系数。为构造SG×E在零假设H0成立时的累积量生成函数(CGF)，首先估计Gi的矩生成函数(MGF)为：
Figure BDA0003589875320000241
Figure BDA0003589875320000242
的一阶导数和二阶导数分别为：
Figure BDA0003589875320000243
Figure BDA0003589875320000244
Gi的累积量生成函数的估计为
Figure BDA0003589875320000245
其一阶导数和二阶导数分别为：
Figure BDA0003589875320000246
因此，SG×E在H0成立时的累积量生成函数的估计为：
Figure BDA0003589875320000247
其一阶导数和二阶导数分别为：
Figure BDA0003589875320000248
Figure BDA0003589875320000249
给定检验统计量SG×E的观测值sG×E，首先计算使得
Figure BDA0003589875320000251
成立的ζ，然后计算得到
Figure BDA0003589875320000252
和
Figure BDA0003589875320000253
为了近似估计统计量SG×E在零假设H0下的分布，根据Barndorff-Nielsen的鞍点近似公式，使用
Figure BDA0003589875320000254
来近似估计概率Pr(SG×E<sG×E)；使用鞍点近似方法计算得到的双侧p值为：
Figure BDA0003589875320000255
其中，
Figure BDA0003589875320000256
和
Figure BDA0003589875320000257
分别表示由上述鞍点近似方法得到的Pr(SG×E<-|sG×E|)和Pr(SG×E<|sG×E|)的估计。
若检验零假设H0:βG×E＝0的统计量为
Figure BDA0003589875320000258
则在假设H0:βG×E＝0成立时将Gi,i≤n视作服从二项分布B(2,q)的独立同分布的随机变量序列，将
Figure BDA0003589875320000259
和Ei,i≤n视作常系数。为构造
Figure BDA00035898753200002510
在零假设H0成立时的累积量生成函数(CGF)，首先估计Gi的矩生成函数(MGF)为：
Figure BDA00035898753200002511
Figure BDA00035898753200002512
的一阶导数和二阶导数分别为：
Figure BDA00035898753200002513
Figure BDA00035898753200002514
Gi的累积量生成函数的估计为
Figure BDA00035898753200002515
其一阶导数和二阶导数分别为：
Figure BDA00035898753200002516
因此，
Figure BDA00035898753200002517
在H0成立时的累积量生成函数的估计为：
Figure BDA00035898753200002518
其一阶导数和二阶导数分别为：
Figure BDA0003589875320000261
Figure BDA0003589875320000262
给定检验统计量
Figure BDA0003589875320000263
的观测值
Figure BDA0003589875320000264
首先计算使得
Figure BDA0003589875320000265
成立的
Figure BDA0003589875320000266
然后计算得到
Figure BDA0003589875320000267
和
Figure BDA0003589875320000268
为了近似估计统计量
Figure BDA0003589875320000269
在零假设H0下的分布，根据Barndorff-Nielsen的鞍点近似公式，使用
Figure BDA00035898753200002610
来近似估计概率
Figure BDA00035898753200002611
因此，使用鞍点近似方法计算得到的双侧p值为：
Figure BDA00035898753200002612
其中，
Figure BDA00035898753200002613
和
Figure BDA00035898753200002614
分别表示由上述鞍点近似方法得到的
Figure BDA00035898753200002615
和
Figure BDA00035898753200002616
的估计。
(3)正态分布近似方法与鞍点近似方法相结合的混合检验策略
在计算统计p值时，鞍点近似方法的运算速度比正态分布近似方法更慢，为了提高计算效率，使用一种将正态分布近似方法与鞍点近似方法相结合的混合检验策略。当统计量的观测值在其分布的期望附近时，鞍点近似方法非常不稳定，而正态分布近似方法在分布的期望附近的近似很准确，因此，当检验统计量观测值在0附近时，使用正态分布近似方法。
令r为一个预先选定的正数(如r＝2)。
若检验零假设H0:βG×E＝0的统计量为SG×E，sG×E为SG×E的观测值，如果
Figure BDA00035898753200002617
则使用所述正态分布近似方法来计算统计p值；否则，如果
Figure BDA0003589875320000271
则使用所述鞍点近似方法来计算统计p值。
若检验零假设H0:βG×E＝0的统计量为
Figure BDA0003589875320000272
为
Figure BDA0003589875320000273
的观测值，如果
Figure BDA0003589875320000274
则使用所述正态分布近似方法来计算统计p值；否则，如果
Figure BDA0003589875320000275
则使用所述鞍点近似方法来计算统计p值。
获得统计p值后，在给定的显著性水平下，进行全基因组范围的关联分析。
由于正态分布近似方法的运算时间比鞍点近似方法要更少，因此，将正态分布近似方法与鞍点近似方法相结合的混合检验策略在保证统计推断准确的同时也兼顾了运算效率。
对于某一性状的基因-环境交互作用分析，选择分析所需的广义线性回归模型后，只需获取在假设Hc:βG＝βG×E＝0下拟合模型得到的残差，本实施例中的正态分布近似方法与鞍点近似方法相结合的混合检验策略就可以在全基因组的范围内应用并计算统计p值，从而进行全基因组范围的关联分析。因此，本发明的应用范围十分广泛，适用于多种复杂性状的基因-环境交互作用的分析。
本实施例还公开了一种计算机可读存储介质，包括供电子设备的一个或多个处理器执行的一个或多个程序，所述一个或多个程序被处理器执行时可实现如上所述适用复杂性状的基因-环境交互的分析方法。
可理解的是，本发明实施例提供的***、存储介质及设备与本发明实施例提供的方法相对应，相关内容的解释、举例和有益效果可以参考上述方法中的相应部分。
本领域普通技术人员可以理解：实现上述方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成，前述的程序可以存储于一计算机可读取存储介质中，该程序在执行时，执行包括上述方法实施例的步骤；而前述的存储介质包括：ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。
以下提供一个具体数值模拟示例：
本实施例提供了一个数值模拟示例以验证本文所提方法的有效性。本实施例将所提算法应用于针对生存数据表型的基因-环境交互作用分析，并进行数值模拟对所提方法的第一类错误率和统计效力进行评估。
具体实现步骤如下：
(一)生存数据表型的生成
在数值模拟中，将个体样本量设置为n＝10000，且个体之间相互独立。对于个体i，首先生成事件的发生时间
Figure BDA0003589875320000281
和删失时间Ci，然后计算生存数据表型
Figure BDA0003589875320000282
和指示量
Figure BDA0003589875320000283
删失时间Ci由尺度参数为0.15、形状参数为1的Weibull分布模拟生成。事件的发生时间
Figure BDA0003589875320000284
由一个具有Weibull基准风险函数的Cox比例风险模型生成：
Figure BDA0003589875320000285
其中，Ui由均匀分布U(0,1)生成，线性预测项ηi＝0.5X1i+0.5X2i+0.5Ei+GiβG+GiEiβG×E，混杂因素X1i由伯努利分布Bernoulli(0.5)生成，混杂因素X2i由标准正态分布生成，环境因素Ei由标准正态分布生成，基因型Gi服从哈迪-温伯格定律，即由二项分布B(2,q)生成，其中q是次要等位基因频率(MAF)，参数βG表示边际遗传效应，参数βG×E表示边际基因-环境交互作用，尺度参数α根据所给定的事件发生率(ER)进行选择。在本实施例中，第一类错误率和统计效力的数值模拟考虑1％、10％和50％三个事件发生率。
(二)生存数据表型分析的第一类错误率数值模拟：
为了评估第一类错误率，在边际基因-环境交互作用βG×E＝0的情况下模拟生成表型。
针对位点的边际基因-环境交互作用为0、边际遗传效应不为0的情形，首先生成1000个相互独立的位点的基因型，每个位点的MAF由相互独立的均匀分布U(0.05,0.5)生成，将线性预测项设置为
Figure BDA0003589875320000291
Figure BDA0003589875320000292
并使用其生成表型，其中Gki表示所生成的第k个位点的基因型取值，边际遗传效应
Figure BDA0003589875320000293
由相互独立的均匀分布U(-0.4,0.4)生成。对于每一个事件发生率(ER)，模拟生成了1000个包含表型、混杂因素及环境因素的数据集。因此，针对边际基因-环境交互作用为0、边际遗传效应不为0的情形，对于每个事件发生率，共进行了106次对于边际基因-环境交互作用的检验。
针对位点的边际基因-环境交互作用和边际遗传效应均为0的情形，在生成位点的基因型时，考虑0.3、0.05和0.005三个次要等位基因频率(MAF)，对于每一个MAF，使用二项分布B(2,MAF)生成104个相互独立的不具有边际基因-环境交互作用和边际遗传效应的位点的基因型，因此，针对边际基因-环境交互作用和边际遗传效应均为0的情形，对于每一对次要等位基因频率(MAF)和事件发生率(ER)，共进行了107次对于边际基因-环境交互作用的检验。
在显著性水平α＝5×10-5下，针对位点的边际基因-环境交互作用和边际遗传效应均为0的情形，对正态分布近似方法与鞍点近似方法相结合的混合检验方法(EmpGxE)和正态分布近似方法(Normal)的经验第一类错误率进行了评估比较。此外，为了进一步评估第一类错误率，针对位点的边际基因-环境交互作用和边际遗传效应均为0的情形以及位点的边际基因-环境交互作用为0、边际遗传效应不为0的情形，对正态分布近似与鞍点近似相结合的混合检验方法(EmpGxE)计算得到的统计p值以及正态分布近似方法(Normal)计算得到的统计p值进行了评估比较。
(三)生存数据表型分析的统计效力数值模拟：
为了评估统计效力，在边际基因-环境交互作用βG×E≠0的模型下生成表型。将边际遗传效应设置为βG＝0，将边际基因-环境交互作用βG×E取值范围设置为[0,3]。考虑0.3、0.05和0.005三个次要等位基因频率(MAF)，对于每一个MAF，使用二项分布B(2,MAF)生成1000个相互独立的具有边际基因-环境交互作用的位点的基因型。对于每一组次要等位基因频率(MAF)和事件发生率(ER)和βG×E，使用下述线性预测项来模拟生成表型
ηi＝0.5X1i+0.5X2i+0.5Ei+GiEiβG×E
并进行了1000次对于边际基因-环境交互作用的检验。本实施例在显著性水平α＝5×10-5对正态分布近似方法与鞍点近似方法相结合的混合检验方法(EmpGxE)和正态分布近似方法(Normal)的统计效力进行了评估比较。
(四)生存数据表型分析的第一类错误率数值模拟结果：
针对位点的边际基因-环境交互作用和边际遗传效应均为0的情形，上述基于107次检验的第一类错误率数值模拟仿真的结果如图3所示。在显著性水平α＝5×10-5下，当检验低频位点或者罕见位点(MAF＝0.05和0.005)时，尤其是在事件发生率很低(ER＝1％)的情形，正态分布近似方法(Normal)会导致偏大的第一类错误率；而正态分布近似与鞍点近似相结合的混合检验方法(EmpGxE)在任何次要等位基因频率(MAF)和事件发生率(ER)的设置下都可以相对较好地控制第一类错误率。
针对位点的边际基因-环境交互作用和边际遗传效应均为0的情形，上述基于107检验所获得的统计p值的分位图如图4所示。正如所料，在表型分布平衡，即事件发生率ER＝50％的情形下，无论对于常见位点、低频位点或者罕见位点，正态分布近似与鞍点近似相结合的混合检验方法(EmpGxE)以及正态分布近似方法(Normal)都可以得到校准良好的统计p值，即p值的经验分布近似服从均匀分布U(0,1)；在事件发生率(ER)很低的情形下，尤其当所检验的位点的次要等位基因频率(MAF)较低时，正态分布近似方法(Normal)会得到偏小的统计p值；而正态分布近似与鞍点近似相结合的混合检验方法(EmpGxE)在任何次要等位基因频率(MAF)和事件发生率(ER)的设置下都可以得到校准良好的统计p值，即p值的经验分布近似服从均匀分布U(0,1)。
针对位点的边际基因-环境交互作用为0、边际遗传效应不为0的情形，图5展示了上述基于106次检验所获得的p值的分位图，正态分布近似与鞍点近似相结合的混合检验方法(EmpGxE)以及正态分布近似方法(Normal)都可以得到校准良好的统计p值，即p值的经验分布近似服从均匀分布U(0,1)。
(五)生存数据表型分析的统计效力数值模拟结果：
图6展示了在显著性水平α＝5×10-5下，正态分布近似方法与鞍点近似方法相结合的混合检验方法(EmpGxE)和正态分布近似方法(Normal)的经验统计效力。在所有次要等位基因频率(MAF)和事件发生率(ER)的设置下，正态分布近似方法与鞍点近似方法相结合的混合检验方法(EmpGxE)和正态分布近似方法(Normal)的经验统计效力几乎完全相同。正如所料，事件发生率越高、位点的MAF越大，经验统计效力越高；事件发生率越低、位点的MAF越小，经验统计效力越低。
针对生存数据的基因-环境交互作用分析的数值模拟仿真结果表明：所提出的正态分布近似与鞍点近似相结合的混合检验的方法(EmpGxE)可以控制第一类错误率；在所检验存位点在边际遗传效应的情形也可以适应于边际遗传效应，得到良好校准的统计p值；并且，也具有较为理想的统计效力。
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种适用复杂性状的基因-环境交互的分析方法，其特征在于，包括以下步骤：
获取分析所需的表型数据、基因型数据、环境因素数据和混杂因素数据；
基于表型数据、基因型数据、环境因素数据和混杂因素数据，选择分析所需的广义线性回归模型，在边际遗传效应与边际基因-环境交互作用均为0的假设下拟合约束模型、估计模型参数并计算残差；
对于待检验的位点，估计其次要等位基因频率，基于score统计量检验其边际遗传效应的显著性；
根据边际遗传效应的显著性来选择并计算相应的用于检验边际基因-环境交互作用显著性的检验统计量；
对于检验边际基因-环境交互作用的检验统计量，使用正态分布近似方法与鞍点近似方法相结合的混合检验策略计算统计p值以检验边际基因-环境交互作用的显著性，从而实现全基因组范围的关联分析。
2.如权利要求1所述的适用复杂性状的基因-环境交互的分析方法，其特征在于，所述基于表型数据、基因型数据、环境因素数据和混杂因素数据，选择分析所需的广义线性回归模型，在边际遗传效应与边际基因-环境交互作用均为0的假设下拟合约束模型、估计模型参数并计算残差，包括：
假设参与研究的个体样本数量为n，对于个体i，1≤i≤n，令Xi表示一个k维的混杂因素向量，Gi表示基因型，Gi＝0，1，2，Ei表示环境因素，Yi表示某一表型，线性预测项
Figure FDA0003589875310000011
其中βX表示混杂因素Xi的k维系数向量，βE、βG和βG×E分别表示环境因素Ei、基因型Gi和基因-环境交互项GiEi所对应的系数，假设表型Yi的期望可以表示成ηi的函数，此假设对于GWAS中大部分的模型是成立的，对数似然函数l(βX，βE，βG，βG×E)＝lnL(βX，βE，βG，βG×E)是ηi的函数，表示为
Figure FDA0003589875310000021
的形式，其中f(.)表示某个一元函数；对于不同的复杂性状，采用与表型相对应的广义线性回归模型进行分析；
为了检验边际基因-环境交互作用的显著性，在假设Hc：βG＝βG×E＝0，即边际遗传效应和边际基因-环境交互作用均为0的假设下拟合模型，得到参数(βX，βE)的极大似然估计
Figure FDA0003589875310000022
并定义残差向量
Figure FDA0003589875310000023
其中，
Figure FDA0003589875310000024
f′(.)表示f(.)的一阶导数。
3.如权利要求2所述的适用复杂性状的基因-环境交互的分析方法，其特征在于，所述对于待检验的位点，估计其次要等位基因频率(MAF)，基于score统计量检验其边际遗传效应的显著性，包括：
假设待检验位点的次要等位基因频率为q，令
Figure FDA0003589875310000025
作为q的估计；
定义基因型向量G和基因-环境交互向量GE分别为：
Figure FDA0003589875310000026
计算用于检验边际遗传效应的score检验统计量
Figure FDA0003589875310000027
并认为检验统计量
Figure FDA0003589875310000031
在假设HG：βG＝0成立时服从期望为0、方差为
Figure FDA0003589875310000032
的正态分布，其中
Figure FDA0003589875310000033
表示
Figure FDA0003589875310000034
的方差
Figure FDA0003589875310000035
的估计，令
Figure FDA0003589875310000036
作为使用正态分布近似方法检验假设HG：βG＝0得到的双侧p值，其中
Figure FDA0003589875310000037
为
Figure FDA0003589875310000038
的观测值，Φ(.)表示标准正态分布的累积分布函数。
4.如权利要求3所述的适用复杂性状的基因-环境交互的分析方法，其特征在于，所述根据边际遗传效应的显著性来选择并计算相应的用于检验边际基因-环境交互作用显著性的检验统计量，包括：
检验边际基因-环境交互作用的显著性即检验零假设H0：βG×E＝0是否成立；
令∈为某一给定的介于0和1之间的正数；
若对假设HG：βG＝0进行检验得到的p值小于或等于∈，则使用
Figure FDA0003589875310000039
作为检验零假设H0：βG×E＝0的检验统计量，其中
Figure FDA00035898753100000310
若对假设HG：βG＝0进行检验得到的p值大于∈，则使用
Figure FDA00035898753100000311
Figure FDA00035898753100000312
作为检验零假设H0：βG×E＝0的检验统计量，其中
Figure FDA00035898753100000313
Figure FDA00035898753100000314
为适应于基因型的残差向量，In为n阶单位矩阵，W＝(1n，G)为n×2矩阵，1n＝(1，1，...，1)T为元素全部为1的n维列向量。
5.如权利要求4所述的适用复杂性状的基因-环境交互的分析方法，其特征在于，所述正态分布近似方法，包括：
若检验零假设H0：βG×E＝0的统计量为SG×E，则正态分布近似方法使用期望为0，方差为
Figure FDA0003589875310000041
的正态分布对零假设H0：βG×E＝0成立时统计量SG×E的分布进行近似估计，其中，
Figure FDA0003589875310000042
为SG×E的方差Var(SG×E)的估计，使用正态分布近似方法计算得到的统计p值为：
Figure FDA0003589875310000043
其中sG×E为SG×E的观测值，Φ(.)是标准正态分布的累积分布函数；
若检验零假设H0：βG×E＝0的统计量为
Figure FDA0003589875310000044
则正态分布近似方法使用期望为0，方差为
Figure FDA0003589875310000045
的正态分布对零假设H0：βG×E＝0成立时统计量
Figure FDA0003589875310000046
的分布进行近似估计，其中，
Figure FDA0003589875310000047
为
Figure FDA0003589875310000048
的方差
Figure FDA0003589875310000049
的估计，使用正态分布近似方法计算得到统计p值为：
Figure FDA00035898753100000410
其中
Figure FDA00035898753100000411
为
Figure FDA00035898753100000412
的观测值，Φ(.)是标准正态分布的累积分布函数。
6.如权利要求5所述的适用复杂性状的基因-环境交互的分析方法，其特征在于，所述鞍点近似方法，包括：
若检验零假设H0：βG×E＝0的统计量为SG×E，则在零假设H0：βG×E＝0成立时将Gi，i≤n视作服从二项分布B(2，q)的独立同分布的随机变量序列，将Ri，i≤n和Ei，i≤n视作常系数；为构造SG×E在零假设H0成立时的累积量生成函数(CGF)，首先估计Gi的矩生成函数(MGF)为：
Figure FDA00035898753100000413
Figure FDA00035898753100000414
的一阶导数和二阶导数分别为：
Figure FDA0003589875310000051
Figure FDA0003589875310000052
Gi的累积量生成函数的估计为
Figure FDA0003589875310000053
其一阶导数和二阶导数分别为：
Figure FDA0003589875310000054
因此，SG×E在H0成立时的累积量生成函数的估计为：
Figure FDA0003589875310000055
其一阶导数和二阶导数分别为：
Figure FDA0003589875310000056
Figure FDA0003589875310000057
给定检验统计量SG×F的观测值sG×E，首先计算使得
Figure FDA0003589875310000058
的ζ，然后计算得到
Figure FDA0003589875310000059
和
Figure FDA00035898753100000510
为了近似估计统计量SG×E在零假设H0下的分布，根据Barndorff-Nielsen鞍点近似公式，使用
Figure FDA00035898753100000511
来近似估计概率Pr(SG×E＜sG×E)，其中Φ(.)是标准正态分布的累积分布函数，因此，使用鞍点近似方法计算得到的双侧p值为：
Figure FDA00035898753100000512
其中，
Figure FDA00035898753100000513
和
Figure FDA00035898753100000514
分别表示由上述鞍点近似方法得到的Pr(SG×E＜-|sG×E|)和Pr(SG×E＜|sG×E|)的估计；
若检验零假设H0：βG×E＝0的统计量为
Figure FDA00035898753100000515
则在零假设H0：βG×E＝0成立时将Gi，i≤n视作服从二项分布B(2，q)的独立同分布的随机变量序列，将
Figure FDA0003589875310000061
和Ei，i≤n视作常系数；为构造
Figure FDA0003589875310000062
在零假设H0成立时的累积量生成函数，首先估计Gi的矩生成函数为：
Figure FDA0003589875310000063
Figure FDA0003589875310000064
的一阶导数和二阶导数分别为：
Figure FDA0003589875310000065
Figure FDA0003589875310000066
Gi的累积量生成函数的估计为
Figure FDA0003589875310000067
其一阶导数和二阶导数分别为：
Figure FDA0003589875310000068
因此，
Figure FDA0003589875310000069
在H0成立时的累积量生成函数的估计为：
Figure FDA00035898753100000610
其一阶导数和二阶导数分别为：
Figure FDA00035898753100000611
Figure FDA00035898753100000612
给定检验统计量
Figure FDA00035898753100000613
的观测值
Figure FDA00035898753100000614
首先计算使得
Figure FDA00035898753100000615
的
Figure FDA00035898753100000616
然后计算得到
Figure FDA00035898753100000617
和
Figure FDA00035898753100000618
为了近似估计统计量
Figure FDA00035898753100000619
在零假设H0下的分布，根据Barndorff-Nielsen鞍点近似公式，使用
Figure FDA00035898753100000620
来近似估计概率
Figure FDA00035898753100000621
其中Φ(.)是标准正态分布的累积分布函数，因此，使用鞍点近似方法计算得到的双侧p值为：
Figure FDA00035898753100000622
其中，
Figure FDA0003589875310000071
和
Figure FDA0003589875310000072
分别表示由上述鞍点近似方法得到的
Figure FDA0003589875310000073
和
Figure FDA0003589875310000074
的估计。
7.如权利要求5和6所述的适用复杂性状的基因-环境交互的分析方法，其特征在于，所述对于检验边际基因-环境交互作用的检验统计量，使用正态分布近似方法与鞍点近似方法相结合的混合检验策略计算统计p值以检验边际基因-环境交互作用的显著性，从而实现全基因组范围的关联分析，包括：
令r为一个预先选定的正数；
若检验零假设H0：βG×E＝0的统计量为SG×E，sG×E为SG×E的观测值，如果
Figure FDA0003589875310000075
则使用所述正态分布近似方法来计算统计p值；否则，如果
Figure FDA0003589875310000076
则使用所述鞍点近似方法来计算统计p值；
若检验零假设H0：βG×E＝0的统计量为
Figure FDA0003589875310000077
为
Figure FDA0003589875310000078
的观测值，如果
Figure FDA0003589875310000079
则使用所述正态分布近似方法来计算统计p值；否则，如果
Figure FDA00035898753100000710
则使用所述鞍点近似方法来计算统计p值；
获得统计p值后，在给定的显著性水平下，进行全基因组范围的关联分析。
8.一种计算机可读存储介质，存储有计算机程序，所述计算机程序被处理器执行时，使得所述处理器执行如权利要求1至7中任一项所述方法的步骤。
CN202210373636.5A 2022-04-11 2022-04-11 适用复杂性状的基因-环境交互的分析方法及存储介质 Active CN114898809B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210373636.5A CN114898809B (zh) 2022-04-11 2022-04-11 适用复杂性状的基因-环境交互的分析方法及存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210373636.5A CN114898809B (zh) 2022-04-11 2022-04-11 适用复杂性状的基因-环境交互的分析方法及存储介质

Publications (2)

Publication Number Publication Date
CN114898809A true CN114898809A (zh) 2022-08-12
CN114898809B CN114898809B (zh) 2022-12-23

Family

ID=82714844

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210373636.5A Active CN114898809B (zh) 2022-04-11 2022-04-11 适用复杂性状的基因-环境交互的分析方法及存储介质

Country Status (1)

Country Link
CN (1) CN114898809B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117831637A (zh) * 2024-03-05 2024-04-05 中国农业科学院作物科学研究所 一种基于机器学习的基因型和环境互作方法及其应用

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102369531A (zh) * 2009-02-06 2012-03-07 先正达参股股份有限公司 用于选择统计上确认的候选基因的方法
CN102758010A (zh) * 2012-06-07 2012-10-31 中国医学科学院阜外心血管病医院 与冠心病相关的多个基因单核苷酸多态性位点与环境因素组合及其应用
CN105740649A (zh) * 2016-01-22 2016-07-06 浙江大学 一种基于混合线性模型的多性状关联分析方法
CN107066781A (zh) * 2016-11-03 2017-08-18 西南大学 基于遗传和环境相关的结直肠癌数据模型的分析方法
CN110059749A (zh) * 2019-04-19 2019-07-26 成都四方伟业软件股份有限公司 重要特征的筛选方法、装置及电子设备
AU2018375721A1 (en) * 2017-11-29 2020-06-11 Intuit Inc. System and method for generating aggregated statistics over sets of user data while enforcing data governance policy
CN112185464A (zh) * 2020-09-29 2021-01-05 山东大学 一种基于孟德尔随机化的多性状全转录组关联分析方法、***及其应用
CN112582023A (zh) * 2020-12-17 2021-03-30 河南省农业科学院粮食作物研究所 一种基于全基因组关联分析和多环境预测模型的玉米分子育种方法
CN113223606A (zh) * 2021-05-13 2021-08-06 浙江大学 一种用于复杂性状遗传改良的基因组选择方法

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102369531A (zh) * 2009-02-06 2012-03-07 先正达参股股份有限公司 用于选择统计上确认的候选基因的方法
CN102758010A (zh) * 2012-06-07 2012-10-31 中国医学科学院阜外心血管病医院 与冠心病相关的多个基因单核苷酸多态性位点与环境因素组合及其应用
CN105740649A (zh) * 2016-01-22 2016-07-06 浙江大学 一种基于混合线性模型的多性状关联分析方法
CN107066781A (zh) * 2016-11-03 2017-08-18 西南大学 基于遗传和环境相关的结直肠癌数据模型的分析方法
AU2018375721A1 (en) * 2017-11-29 2020-06-11 Intuit Inc. System and method for generating aggregated statistics over sets of user data while enforcing data governance policy
CN110059749A (zh) * 2019-04-19 2019-07-26 成都四方伟业软件股份有限公司 重要特征的筛选方法、装置及电子设备
CN112185464A (zh) * 2020-09-29 2021-01-05 山东大学 一种基于孟德尔随机化的多性状全转录组关联分析方法、***及其应用
CN112582023A (zh) * 2020-12-17 2021-03-30 河南省农业科学院粮食作物研究所 一种基于全基因组关联分析和多环境预测模型的玉米分子育种方法
CN113223606A (zh) * 2021-05-13 2021-08-06 浙江大学 一种用于复杂性状遗传改良的基因组选择方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
WENJIAN BI ET AL.: "A Fast and Accurate Method for Genome-wide Scale Phenome-wide G×E Analysis and Its Application to UK Biobank", 《THE AMERICAN JOURNAL OF HUMAN GENETICS》 *
WENJIAN BI ET AL.: "A Fast and Accurate Method for Genome-Wide Time-to-Event Data Analysis and Its Application to UK Biobank", 《THE AMERICAN JOURNAL OF HUMAN GENETICS》 *
苏矿喜: "基因与基因交互效应检测高效统计方法研究", 《中国优秀硕士学位论文全文数据库 基础科学辑》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117831637A (zh) * 2024-03-05 2024-04-05 中国农业科学院作物科学研究所 一种基于机器学习的基因型和环境互作方法及其应用
CN117831637B (zh) * 2024-03-05 2024-05-28 中国农业科学院作物科学研究所 一种基于机器学习的基因型和环境互作方法及其应用

Also Published As

Publication number Publication date
CN114898809B (zh) 2022-12-23

Similar Documents

Publication Publication Date Title
Blangero et al. A kernel of truth: statistical advances in polygenic variance component models for complex human pedigrees
US20210375392A1 (en) Machine learning platform for generating risk models
Chen et al. Inference on the order of a normal mixture
Tan et al. Bayesian inference for multiple Gaussian graphical models with application to metabolic association networks
Ockuly et al. Response surface experiments: A meta-analysis
Overstall et al. conting: An R package for Bayesian analysis of complete and incomplete contingency tables
CN114898809B (zh) 适用复杂性状的基因-环境交互的分析方法及存储介质
Shin et al. Statistical power for identifying nucleotide markers associated with quantitative traits in genome-wide association analysis using a mixed model
Cartwright et al. A family-based probabilistic method for capturing de novo mutations from high-throughput short-read sequencing data
Montinaro et al. Revisiting the out of Africa event with a deep-learning approach
Hopkins et al. Phenotypic screening models for rapid diagnosis of genetic variants and discovery of personalized therapeutics
Sztepanacz et al. Heritable micro-environmental variance covaries with fitness in an outbred population of Drosophila serrata
Meyer et al. Predictive variable selection in generalized linear models
Sesia et al. Controlling the false discovery rate in GWAS with population structure
Ballard et al. Shared components of heritability across genetically correlated traits
Achlioptas et al. Two-locus association mapping in subquadratic time
Montazeri et al. Shrinkage estimation of effect sizes as an alternative to hypothesis testing followed by estimation in high-dimensional biology: Applications to differential gene expression
St-Pierre et al. Efficient penalized generalized linear mixed models for variable selection and genetic risk prediction in high-dimensional data
CN116756510A (zh) 可控制样本中亲缘相关性的分析方法、***及储存介质
Kirchler et al. Training normalizing flows from dependent data
CN114171110A (zh) 一种基于联合似然的孟德尔随机化分析方法
Cui et al. Improving neural networks for genotype-phenotype prediction using published summary statistics
Bhattacharya et al. Effects of gene–environment and gene–gene interactions in case-control studies: A novel Bayesian semiparametric approach
US20050177316A1 (en) Algorithm for estimating and testing association between a haplotype and quantitative phenotype
Tyagi et al. Application of Genetic Algorithms for Medical Diagnosis of Diabetes Mellitus

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant