CN114898809A - 适用复杂性状的基因-环境交互的分析方法及存储介质 - Google Patents
适用复杂性状的基因-环境交互的分析方法及存储介质 Download PDFInfo
- 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
Links
- 230000003993 interaction Effects 0.000 title claims abstract description 106
- 238000004458 analytical method Methods 0.000 title claims abstract description 43
- 238000000034 method Methods 0.000 claims abstract description 144
- 238000009826 distribution Methods 0.000 claims abstract description 130
- 238000012360 testing method Methods 0.000 claims abstract description 102
- 230000002068 genetic effect Effects 0.000 claims abstract description 68
- 108700028369 Alleles Proteins 0.000 claims abstract description 27
- 230000007613 environmental effect Effects 0.000 claims abstract description 23
- 238000012098 association analyses Methods 0.000 claims abstract description 14
- 238000012417 linear regression Methods 0.000 claims abstract description 12
- 230000001186 cumulative effect Effects 0.000 claims description 27
- 238000005315 distribution function Methods 0.000 claims description 13
- 238000001604 Rao's score test Methods 0.000 claims description 8
- 238000007476 Maximum Likelihood Methods 0.000 claims description 6
- 239000011159 matrix material Substances 0.000 claims description 6
- 238000010219 correlation analysis Methods 0.000 claims description 5
- 108010014173 Factor X Proteins 0.000 claims description 4
- 230000002146 bilateral effect Effects 0.000 claims description 4
- 238000004590 computer program Methods 0.000 claims description 3
- 239000002184 metal Substances 0.000 claims description 3
- 101150043283 ccdA gene Proteins 0.000 claims 1
- 230000006870 function Effects 0.000 description 25
- 238000004088 simulation Methods 0.000 description 18
- 230000004083 survival effect Effects 0.000 description 18
- 238000010998 test method Methods 0.000 description 13
- 238000011160 research Methods 0.000 description 9
- 238000012482 interaction analysis Methods 0.000 description 8
- 239000000523 sample Substances 0.000 description 8
- 238000004422 calculation algorithm Methods 0.000 description 6
- 238000004364 calculation method Methods 0.000 description 6
- 238000009827 uniform distribution Methods 0.000 description 6
- 239000012472 biological sample Substances 0.000 description 4
- 201000010099 disease Diseases 0.000 description 4
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 4
- 230000008569 process Effects 0.000 description 4
- 108090000623 proteins and genes Proteins 0.000 description 4
- 238000012217 deletion Methods 0.000 description 3
- 230000037430 deletion Effects 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 230000018109 developmental process Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 230000007614 genetic variation Effects 0.000 description 3
- 230000036541 health Effects 0.000 description 3
- 238000012163 sequencing technique Methods 0.000 description 3
- 238000011282 treatment Methods 0.000 description 3
- 238000003657 Likelihood-ratio test Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 238000007477 logistic regression Methods 0.000 description 2
- 230000002265 prevention Effects 0.000 description 2
- 208000017667 Chronic Disease Diseases 0.000 description 1
- 206010020772 Hypertension Diseases 0.000 description 1
- 241000764238 Isis Species 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000002776 aggregation Effects 0.000 description 1
- 238000004220 aggregation Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 210000004556 brain Anatomy 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 206010012601 diabetes mellitus Diseases 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000011337 individualized treatment Methods 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000010197 meta-analysis Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000002773 nucleotide Substances 0.000 description 1
- 125000003729 nucleotide group Chemical group 0.000 description 1
- 230000005868 ontogenesis Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 239000003960 organic solvent Substances 0.000 description 1
- 230000008506 pathogenesis Effects 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
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表示某一表型,线性预测项 其中β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的函数,可以表示为的形式,其中f(.)表示某个一元函数,对于不同的复杂性状,采用与表型相对应的广义线性回归模型进行分析;
为了检验边际基因-环境交互作用的显著性,在假设Hc:βG=βG×E=0,即边际遗传效应和边际基因-环境交互作用均为0的假设下拟合模型,得到参数(βX,βE)在假设Hc:βG=βG×E=0下的极大似然估计并定义残差向量
进一步地,所述对于待检验的位点,估计其次要等位基因频率(MAF),基于score统计量检验其边际遗传效应的显著性,包括:
定义基因型向量G和基因-环境交互向量GE分别为:
计算用于检验边际遗传效应的score检验统计量
并认为检验统计量在假设HG:βG=0成立时服从期望为0、方差为的正态分布,其中表示的方差的估计,令作为使用正态分布近似方法检验假设HG:βG=0得到的双侧p值,其中为的观测值,Φ(.)表示标准正态分布的累积分布函数。
进一步地,所述根据边际遗传效应的显著性来选择并计算相应的用于检验边际基因-环境交互作用显著性的检验统计量,包括:
检验边际基因-环境交互作用的显著性即检验零假设H0:βG×E=0是否成立;
令∈为某一给定的介于0和1之间的正数(如∈=0.001);
若对假设HG:βG=0进行检验得到的p值大于∈,则使用 作为检验零假设H0:βG×E=0的检验统计量,其中 为适应于基因型的残差向量,In为n阶单位矩阵,W=(1n,G)为n×2矩阵,1n=(1,1,…,1)T为元素全部为1的n维列向量。
进一步地,所述正态分布近似方法,包括:
若检验零假设H0:βG×E=0的统计量为SG×E,则正态分布近似方法使用期望为0,方差为的正态分布对零假设H0:βG×E=0成立时统计量SG×E的分布进行近似估计,计算得到统计p值为:其中sG×E为SG×E的观测值,Φ(.)是标准正态分布的累积分布函数;
若检验零假设H0:βG×E=0的统计量为则正态分布近似方法使用期望为0,方差为的正态分布对零假设H0:βG×E=0成立时统计量的分布进行近似估计,其中,为的方差的估计,使用正态分布近似方法计算得到统计p值为:其中为的观测值,Φ(.)是标准正态分布的累积分布函数。
进一步地,所述鞍点近似方法,包括:
若检验零假设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)为:
因此,SG×E在H0成立时的累积量生成函数的估计为:
其一阶导数和二阶导数分别为:
给定检验统计量SG×E的观测值sG×E,首先计算使得的ζ,然后计算得到和为了近似估计统计量SG×E在零假设H0下的分布,根据Barndorff-Nielsen鞍点近似公式,使用来近似估计概率Pr(SG×E<sG×E),其中Φ(.)是标准正态分布的累积分布函数,因此,使用鞍点近似方法计算得到的双侧p值为:
若检验零假设H0:βG×E=0的统计量为则在零假设H0:βG×E=0成立时将Gi,i≤n视作服从二项分布B(2,q)的独立同分布的随机变量序列,将和Ei,i≤n视作常系数;为构造在零假设H0成立时的累积量生成函数(CGF),首先估计Gi的矩生成函数(MGF)为:
其一阶导数和二阶导数分别为:
给定检验统计量的观测值首先计算使得的然后计算得到和为了近似估计统计量在零假设H0下的分布,根据Barndorff-Nielsen鞍点近似公式,使用来近似估计概率其中Φ(.)是标准正态分布的累积分布函数,因此,使用鞍点近似方法计算得到的双侧p值为:
进一步地,所述对于检验边际基因-环境交互作用的检验统计量,使用正态分布近似方法与鞍点近似方法相结合的混合检验策略计算统计p值以检验边际基因-环境交互作用的显著性,从而实现全基因组范围的关联分析,包括:
令r为一个预先选定的正数(如r=2);
获得统计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的期望可以表示成线性预测项的函数,其中βX表示混杂因素Xi的k维系数向量,βE、βG和βG×E分别表示环境因素Ei、基因型Gi和基因-环境交互项GiEi所对应的系数,此假设对于GWAS中大部分的模型是成立的。对数似然函数l(βX,βE,βG,βG×E)=lnL(βX,βE,βG,βG×E)是ηi的函数,可以表示为的形式,其中f(.)表示某个一元函数;对于不同的复杂性状,采用与表型相对应的广义线性回归模型进行分析。关于边际遗传效应βG和边际基因-环境交互作用βG×E的得分函数分别为l(βX,βE,βG,βG×E)关于βG和关于βG×E的偏导数,具有以下形式:
其中,f'(.)表示f(.)的一阶导数。
为了检验边际基因-环境交互作用的显著性,传统的score检验在零假设H0:βG×E=0下拟合模型并计算参数(βX,βE,βG)的极大似然估计以此来计算用于检验边际基因-环境交互作用的score统计量其中 然而此策略需要对每个位点拟合一次模型,对于全基因组范围的分析,其计算成本非常高昂。
为了提高计算效率,本发明并非在零假设H0:βG×E=0下拟合模型,而在假设Hc:βG=βG×E=0,即在边际遗传效应和边际基因-环境交互作用均为0的假设下拟合模型,进而得到参数(βX,βE在假设Hc:βG=βG×E=0下的极大似然估计并定义残差向量
下面两个例子展示了在生存数据表型和多分类表型的基因-环境交互作用分析中基于表型数据、基因型数据、环境因素数据和混杂因素数据,选择分析所需的广义线性回归模型,在边际遗传效应与边际基因-环境交互作用均为0的假设下拟合模型,估计模型参数并计算得到残差的过程:
(1)生存数据表型:
假设表型是生存数据,参与研究的个体样本数量为n,使用Cox比例风险回归模型进行分析:
其中,λ(t;Xi)是风险函数,λ0(t)是基线风险函数。令Ci表示事件的删失时间,表示事件的发生时间,表示观测到的生存数据,用于显示事件是否发生,其中I(.)表示示性函数。函数f(ηi)和f'(ηi)分别为
f'(ηi)=δi-Λi
为了提高计算效率,不在零假设H0:βG×E=0下拟合模型,而在边际遗传效应与边际基因-环境交互作用均为0的假设Hc:βG=βG×E=0下拟合模型,得到参数(βX,βE)的极大似然估计并计算在假设Hc下拟合模型后得到的线性预测项以及Λi的估计
(2)多分类数据表型:
假设表型是多分类数据,参与研究的个体样本数量为n,令J表示类别数量,令Yi=1,2,…,J表示个体i的多分类表型,使用以下比例优势逻辑回归模型进行分析:
其中,vij=Pr(Yi≤j|Xi,Ei,Gi)表示在给定Xi、Ei和Gi的条件下,Yi≤j的累计概率,分割点ε:ε1<…<εJ=∞用于对数据进行分类。
对于个体i,1≤i≤n,定义一个J×1向量作为多分类表型Yi的一个等价的表达形式;如果Yi=j,那么Yij=1并且中其余的元素均为0。令μij表示Yij的期望,即μij=E(Yij)=Pr(Yij=1)=Pr(Yi=j)=Pr(Yi≤j)-Pr(Yi≤j-1),其中E(.)表示某随机变量的期望;令 函数f(ηi)和f'(ηi)分别为
为了检验边际基因-环境交互作用,传统的score检验在边际基因-环境交互作用为0的零假设H0:βG×E=0下拟合模型,得到参数βX,βE,βG,εj,1≤j≤J的极大似然估计1≤j≤J,并计算接下来分别计算得到vij、μij和Vij在H0:βG×E=0下拟合模型后得到的估计量 以及
为了提高计算效率,不在零假设H0:βG×E=0下拟合模型,而在边际遗传效应与边际基因-环境交互作用均为0的假设Hc:βG=βG×E=0下拟合模型,得到参数βX,βE,εj,1≤j≤J的极大似然估计1≤j≤J,并计算在假设Hc:βG=βG×E=0下拟合模型后得到的线性预测项以及vij、μij和Vij的估计量 和 定义为在假设Hc:βG=βG×E=0下拟合模型后计算得到的残差,令R=(R1,…,Rn)T表示相应的残差向量。
定义基因型向量G和基因-环境交互向量GE分别为:
计算用于检验边际遗传效应的score检验统计量
对于在零假设H0:βG×E=0成立时将残差Ri,i≤n视作常系数,将Gi,i≤n视作服从二项分布B(2,q)的独立同分布的随机变量序列,并认为检验统计量在假设HG:βG=0成立时服从期望为0、方差为的正态分布,其中表示的方差的估计,令作为使用正态分布近似方法检验假设HG:βG=0得到的双侧p值,其中表示的观测值,Φ(.)表示标准正态分布的累积分布函数。
具体来说,在上述步骤S4中,为了在假设Hc:βG=βG×E=0下拟合模型并对零假设H0:βG×E=0进行检验,一种自然的方式是使用
作为检验统计量。然而,对于边际遗传效应βG≠0的情形,此策略是不可行的,因为它忽略了边际遗传效应对表型的作用。为了适应于边际遗传效应,需要提出新的检验统计量来对零假设H0:βG×E=0进行检验。下面根据在步骤S3中对假设HG:βG=0进行检验得到的p值,来选择并计算不同的检验统计量对零假设H0进行检验。
令∈为某一给定的介于0和1之间的正数(如∈=0.001)。
若在步骤S3中对假设HG:βG=0进行检验得到的统计p值小于或等于∈,则为了适应于边际遗传效应,提出一个新的检验统计量:
对于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满足线性约束条件和因此在假设Hc:βG=βG×E=0成立时,和的期望和方差分别为:
在假设Hc:βG=βG×E=0成立时,S1和S2之间的协方差为:
其中,Ec(.)、Varc(.)和Covc(.,.)分别表示某统计量在假设Hc成立时的期望、方差以及某两个统计量之间的协方差。因此,λ可以表示为
这表明在假设Hc成立时,通过减去边际遗传效应的影响被从中去除了。在边际遗传效应在0附近的情形 因此,在此情形下可以近似认为检验统计量SG×E仅包含了边际基因-环境交互作用的信息;对于边际遗传效应远离0的情形,λ与相差较大,检验统计量SG×E仍然包含了较多的边际遗传效应的信息,在此情形下使用SG×E检验边际基因-环境交互作用会得到不可靠的统计推断。
若在步骤S2中对假设HG:βG=0进行检验得到的统计p值大于∈,则使用作为检验零假设H0:βG×E=0的检验统计量,其中为适应于基因型的残差向量,In为n阶单位矩阵,W=(1n,G)为n×2矩阵,1n=(1,1,…,1)T为元素全部为1的n维列向量。可以被解释为残差向量R对基因型向量G进行线性回归并使用最小二乘法估计模型参数后得到的新的残差向量。因此,在统计量中,边际遗传效应已通过对R的变换(In-W(WTW)-1WT)R而去除,是仅包含边际基因-环境交互作用的信息的检验统计量。
注意到,计算统计量SG×E比计算统计量的速度更快,因为在计算的过程中需要更多的时间用于计算但是相较于SG×E,使用统计量可以获得更加准确的统计推断。因此,上述根据检验边际遗传效应得到的统计p值来选择并计算相应的用于检验边际基因-环境交互作用的检验统计量的策略兼顾了运算效率与准确性,在两者间达到了平衡。
具体来说,分别介绍在上述步骤S5中提及的正态分布近似方法、鞍点近似方法以及正态分布近似方法与鞍点近似方法相结合的混合检验策略。
(1)正态分布近似方法
若检验零假设H0:βG×E=0的统计量为SG×E,则正态分布近似方法使用期望为0,方差为的正态分布对零假设H0:βG×E=0成立时统计量SG×E的分布进行近似估计,其中,为SG×E的方差Var(SG×E)的估计,使用正态分布近似方法计算得到的统计p值为:其中sG×E为SG×E的观测值,Φ(.)是标准正态分布的累积分布函数。
若检验零假设H0:βG×E=0的统计量为则正态分布近似方法使用期望为0,方差为的正态分布对零假设H0:βG×E=0成立时统计量的分布进行近似估计,其中,为的方差的估计,使用正态分布近似方法计算得到统计p值为:其中为的观测值,Φ(.)是标准正态分布的累积分布函数。
(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)为:
因此,SG×E在H0成立时的累积量生成函数的估计为:
其一阶导数和二阶导数分别为:
给定检验统计量SG×E的观测值sG×E,首先计算使得成立的ζ,然后计算得到和为了近似估计统计量SG×E在零假设H0下的分布,根据Barndorff-Nielsen的鞍点近似公式,使用来近似估计概率Pr(SG×E<sG×E);使用鞍点近似方法计算得到的双侧p值为:
若检验零假设H0:βG×E=0的统计量为则在假设H0:βG×E=0成立时将Gi,i≤n视作服从二项分布B(2,q)的独立同分布的随机变量序列,将和Ei,i≤n视作常系数。为构造在零假设H0成立时的累积量生成函数(CGF),首先估计Gi的矩生成函数(MGF)为:
其一阶导数和二阶导数分别为:
给定检验统计量的观测值首先计算使得成立的然后计算得到和为了近似估计统计量在零假设H0下的分布,根据Barndorff-Nielsen的鞍点近似公式,使用来近似估计概率因此,使用鞍点近似方法计算得到的双侧p值为:
(3)正态分布近似方法与鞍点近似方法相结合的混合检验策略
在计算统计p值时,鞍点近似方法的运算速度比正态分布近似方法更慢,为了提高计算效率,使用一种将正态分布近似方法与鞍点近似方法相结合的混合检验策略。当统计量的观测值在其分布的期望附近时,鞍点近似方法非常不稳定,而正态分布近似方法在分布的期望附近的近似很准确,因此,当检验统计量观测值在0附近时,使用正态分布近似方法。
令r为一个预先选定的正数(如r=2)。
获得统计p值后,在给定的显著性水平下,进行全基因组范围的关联分析。
由于正态分布近似方法的运算时间比鞍点近似方法要更少,因此,将正态分布近似方法与鞍点近似方法相结合的混合检验策略在保证统计推断准确的同时也兼顾了运算效率。
对于某一性状的基因-环境交互作用分析,选择分析所需的广义线性回归模型后,只需获取在假设Hc:βG=βG×E=0下拟合模型得到的残差,本实施例中的正态分布近似方法与鞍点近似方法相结合的混合检验策略就可以在全基因组的范围内应用并计算统计p值,从而进行全基因组范围的关联分析。因此,本发明的应用范围十分广泛,适用于多种复杂性状的基因-环境交互作用的分析。
本实施例还公开了一种计算机可读存储介质,包括供电子设备的一个或多个处理器执行的一个或多个程序,所述一个或多个程序被处理器执行时可实现如上所述适用复杂性状的基因-环境交互的分析方法。
可理解的是,本发明实施例提供的***、存储介质及设备与本发明实施例提供的方法相对应,相关内容的解释、举例和有益效果可以参考上述方法中的相应部分。
本领域普通技术人员可以理解:实现上述方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成,前述的程序可以存储于一计算机可读取存储介质中,该程序在执行时,执行包括上述方法实施例的步骤;而前述的存储介质包括:ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。
以下提供一个具体数值模拟示例:
本实施例提供了一个数值模拟示例以验证本文所提方法的有效性。本实施例将所提算法应用于针对生存数据表型的基因-环境交互作用分析,并进行数值模拟对所提方法的第一类错误率和统计效力进行评估。
具体实现步骤如下:
(一)生存数据表型的生成
在数值模拟中,将个体样本量设置为n=10000,且个体之间相互独立。对于个体i,首先生成事件的发生时间和删失时间Ci,然后计算生存数据表型和指示量删失时间Ci由尺度参数为0.15、形状参数为1的Weibull分布模拟生成。事件的发生时间由一个具有Weibull基准风险函数的Cox比例风险模型生成:
其中,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)生成,将线性预测项设置为 并使用其生成表型,其中Gki表示所生成的第k个位点的基因型取值,边际遗传效应由相互独立的均匀分布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表示某一表型,线性预测项其中β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的函数,表示为的形式,其中f(.)表示某个一元函数;对于不同的复杂性状,采用与表型相对应的广义线性回归模型进行分析;
4.如权利要求3所述的适用复杂性状的基因-环境交互的分析方法,其特征在于,所述根据边际遗传效应的显著性来选择并计算相应的用于检验边际基因-环境交互作用显著性的检验统计量,包括:
检验边际基因-环境交互作用的显著性即检验零假设H0:βG×E=0是否成立;
令∈为某一给定的介于0和1之间的正数;
5.如权利要求4所述的适用复杂性状的基因-环境交互的分析方法,其特征在于,所述正态分布近似方法,包括:
若检验零假设H0:βG×E=0的统计量为SG×E,则正态分布近似方法使用期望为0,方差为的正态分布对零假设H0:βG×E=0成立时统计量SG×E的分布进行近似估计,其中,为SG×E的方差Var(SG×E)的估计,使用正态分布近似方法计算得到的统计p值为:其中sG×E为SG×E的观测值,Φ(.)是标准正态分布的累积分布函数;
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)为:
因此,SG×E在H0成立时的累积量生成函数的估计为:
其一阶导数和二阶导数分别为:
给定检验统计量SG×F的观测值sG×E,首先计算使得的ζ,然后计算得到和为了近似估计统计量SG×E在零假设H0下的分布,根据Barndorff-Nielsen鞍点近似公式,使用来近似估计概率Pr(SG×E<sG×E),其中Φ(.)是标准正态分布的累积分布函数,因此,使用鞍点近似方法计算得到的双侧p值为:
若检验零假设H0:βG×E=0的统计量为则在零假设H0:βG×E=0成立时将Gi,i≤n视作服从二项分布B(2,q)的独立同分布的随机变量序列,将和Ei,i≤n视作常系数;为构造在零假设H0成立时的累积量生成函数,首先估计Gi的矩生成函数为:
其一阶导数和二阶导数分别为:
给定检验统计量的观测值首先计算使得的然后计算得到和为了近似估计统计量在零假设H0下的分布,根据Barndorff-Nielsen鞍点近似公式,使用来近似估计概率其中Φ(.)是标准正态分布的累积分布函数,因此,使用鞍点近似方法计算得到的双侧p值为:
7.如权利要求5和6所述的适用复杂性状的基因-环境交互的分析方法,其特征在于,所述对于检验边际基因-环境交互作用的检验统计量,使用正态分布近似方法与鞍点近似方法相结合的混合检验策略计算统计p值以检验边际基因-环境交互作用的显著性,从而实现全基因组范围的关联分析,包括:
令r为一个预先选定的正数;
获得统计p值后,在给定的显著性水平下,进行全基因组范围的关联分析。
8.一种计算机可读存储介质,存储有计算机程序,所述计算机程序被处理器执行时,使得所述处理器执行如权利要求1至7中任一项所述方法的步骤。
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)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117831637A (zh) * | 2024-03-05 | 2024-04-05 | 中国农业科学院作物科学研究所 | 一种基于机器学习的基因型和环境互作方法及其应用 |
Citations (9)
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 | 浙江大学 | 一种用于复杂性状遗传改良的基因组选择方法 |
-
2022
- 2022-04-11 CN CN202210373636.5A patent/CN114898809B/zh active Active
Patent Citations (9)
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)
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)
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 |