CN112185464A - 一种基于孟德尔随机化的多性状全转录组关联分析方法、系统及其应用 - Google Patents

一种基于孟德尔随机化的多性状全转录组关联分析方法、系统及其应用 Download PDF

Info

Publication number
CN112185464A
CN112185464A CN202011049640.3A CN202011049640A CN112185464A CN 112185464 A CN112185464 A CN 112185464A CN 202011049640 A CN202011049640 A CN 202011049640A CN 112185464 A CN112185464 A CN 112185464A
Authority
CN
China
Prior art keywords
egger
mopmr
trait
data
pleiotropic
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
CN202011049640.3A
Other languages
English (en)
Other versions
CN112185464B (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.)
Shandong University
Original Assignee
Shandong University
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 Shandong University filed Critical Shandong University
Priority to CN202011049640.3A priority Critical patent/CN112185464B/zh
Publication of CN112185464A publication Critical patent/CN112185464A/zh
Application granted granted Critical
Publication of CN112185464B publication Critical patent/CN112185464B/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
    • G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
    • 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
    • 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
    • G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Evolutionary Biology (AREA)
  • Medical Informatics (AREA)
  • Theoretical Computer Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Biophysics (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Genetics & Genomics (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Physiology (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明提供一种基于孟德尔随机化的多性状全转录组关联分析方法、系统及其应用，属于生物技术和数据分析技术领域。该方法考虑了多个性状之间的相关性，适用于存在高度连锁不平衡的cis‑SNPs，借助广泛使用的Egger假设将水平基因多效性嵌入模型，并且基于极大似然理论进行推断，本发明将其命名为具有Egger假设的多结局概率孟德尔随机化(moPMR‑Egger)方法。模拟和实例分析表明，moPMR‑Egger在检验因果效应和水平多效性时可以控制一类错误，比单性状方法具有更高的检验效能，并且可以用于生物库规模数据的实证分析，计算效率高，因此具有良好的实际应用之价值。

Description

一种基于孟德尔随机化的多性状全转录组关联分析方法、系 统及其应用
技术领域
本发明属于生物技术和数据分析技术领域，具体涉及一种基于孟德尔随机化的多性状全转录组关联分析方法及其应用。
背景技术
公开该背景技术部分的信息仅仅旨在增加对本发明的总体背景的理解，而不必然被视为承认或以任何形式暗示该信息构成已经成为本领域一般技术人员所公知的现有技术。
全转录组关联研究(Transcriptome-wide association study,TWAS)在设计上旨在整合全基因组关联研究(Genome-wide association studies,GWASs)和基因表达谱研究，以探讨复杂疾病的内在分子遗传机制。TWAS将特定基因的顺式SNPs(cis-SNPs)作为工具变量，进而检验该特定基因与GWAS性状的因果效应，因此，TWAS的分析策略可以纳入孟德尔随机化(Mendelian randomization,MR)框架。目前已产生了许多TWAS统计分析方法，包括PrediXcan，TWAS，DPR，TIGAR，SMR，PMR-Egger， FOCUS和UTMOST等。不同的TWAS方法在以下诸多方面存在差异：纳入分析的cis-SNPs数目不同(有些方法使用一个cis-SNP，而有些方法使用所有cis-SNPs)；是否校正了水平基因多效性(有些方法考虑了水平多效性，而有些方法未考虑)；构建的cis-SNPs对基因表达的预测模型假设不同 (有些方法依赖于稀疏性假设，而有些方法依赖于基因多效性假设)；模型推断方式不同(有些方法基于似然进行推断，而有些方法基于两阶段回归策略进行推断)。
尽管当前TWAS分析方法众多，但在本质上都是针对单性状进行分析。然而，GWAS研究通常包含具有共同遗传基础的多个相关性状，许多位点已经被证实与多种表型之间存在多效关联，比如， CACNA1C与双相情感障碍和精神分裂症、SLC39A8与精神分裂症和帕金森病、RGS12与血脂和炎症性肠病。因此，对多个性状的联合分析势必能够提升TWAS的检验效能。事实上，多性状联合分析优于单性状分析已在GWAS中得到了证明。在识别同时影响多个性状的基因时，多性状联合TWAS 分析比单性状TWAS具有更高的检验效能。此外，将多性状间的相关性考虑在内，多性状联合分析也比单性状分析能够更有效地识别只影响其中一个性状的基因。因此，发展多性状联合的TWAS统计分析方法势在必行。
发明内容
针对上述现有技术的不足，本发明提供一种基于孟德尔随机化的多性状全转录组关联分析方法、系统及其应用。基于先前TWAS分析的似然推断框架，将其拓展到多个结局性状，提出了一种新型的多性状TWAS方法。该方法考虑了多个性状之间的相关性，适用于存在高度连锁不平衡(linkage disequilibrium,LD)的cis-SNPs，借助广泛使用的Egger假设将水平基因多效性嵌入模型，并且基于极大似然理论进行推断，本发明将其命名为具有Egger假设的多结局概率孟德尔随机化 (moPMR-Egger)方法。模拟和实例分析表明，moPMR-Egger在检验因果效应和水平多效性时可以控制一类错误，比单性状方法具有更高的检验效能，并且可以用于生物库规模数据的实证分析，计算效率高，因此具有良好的实际应用之价值。
本发明是通过如下技术方案实现的：
本发明的第一个方面，提供一种基于孟德尔随机化的多性状全转录组关联分析方法，所述分析方法包括依托两样本MR框架，在没有样本重叠的两项独立研究中测量暴露和结局，从而估计和检验暴露对结局的因果效应；具体的，一次检验一个基因，并将cis-SNPs作为工具变量，在似然框架下同时估计和检验基因对多个结局性状的因果效应，同时能将潜在的水平多效性嵌入模型。
所述分析方法包括对个体水平数据进行分析和汇总统计量数据进行分析。
其中，当对个体水平数据进行分析时，所述分析方法包括采用如下回归模型：
Figure BDA0002709151730000021
Figure BDA0002709151730000022
Figure BDA0002709151730000023
尽管上述三个公式针对的是两个独立的研究，但参数β将二者紧密的结合起来。式(2)-(3)可以合并为
Figure BDA0002709151730000024
其中，
Figure BDA0002709151730000025
遵循遗传学中常用的多基因建模假设，假设β的所有元素非零，且每个元素都服从正态分布
Figure BDA0002709151730000026
同时，遵循PMR-Egger假设，假设对于每个性状i,每个SNP的水平多效性相等。实际上，这种水平多效性假设在其他TWAS和稳健MR研究中得到了广泛认可和应用，它等价于
Figure BDA0002709151730000027
1p是一个元素全为1的p维向量,γ＝(γ1，…，γk)T。
在上述联合模型中，所关心的参数是基因在多个结局性状上的因果效应α。α的因果效应识别和可解释性可以基于因果推断中的决策理论加以证明。和传统的MR分析一样，α的因果解释至少需要两个假设成立：(1)工具变量与基因表达有关；(2)工具变量与基因表达和每个结局都有关的任何其他混杂因素无关。因为对水平多效性Γ做了明确的模型假设，moPMR-Egger不再要求满足MR的一般排除限制条件(即，工具变量只能通过基因表达途径影响每个结局)。然而，对Γ的模型假设需要满足工具变量强度与直接效应独立的假设(Instrument Strength Independent of Direct Effect,InSIDE), 这个假设要求工具变量-基因表达效应和工具变量-结局效应是相互独立的，有时称为弱排除限制条件。因为混杂因素往往是未知的，这些假设有许多在实践中是无法检验的，本发明所使用的术语“因果效应”，只是为了强调α的估计比Y对
Figure BDA0002709151730000033
的普通回归获得的效应值估计更可信。
在上述模型中，感兴趣的是在水平多效性γ存在时，估计因果效应α并检验原假设H0：α＝0,拒绝原假设α＝0意味着基因表达至少对一个性状有非零的因果效应。另外，还可对水平多效性γ进行估计和检验H0：γ＝0,拒绝原假设γ＝0说明cis-SNPs至少对一个性状有非零的水平多效性。在统计推断层面，发展了期望最大化(EM)算法的参数扩展版本(PX-EM)，通过最大化基于式(1)和(4) 定义的联合似然来进行参数推断。可以获得最大似然联合模型以及α和γ最大似然估计。此外，进一步将该EM算法应用于两个简化模型，一个模型中不包含α，另一个模型中不包含γ，分别获得对应的最大似然。最后，将从联合模型得到的最大似然函数与从两个简化模型得到的最大似然函数进行对比，实现对H0：α＝0或H0：γ＝0的似然比检验。
将此模型和算法一起称为具有Egger假设的多结局概率孟德尔随机化(moPMR-Egger)方法。其中，“mo”指的是对多个结局(multipleoutcome)的建模。“概率”是指moPMR-Egger模型完全模仿了TWAS背后的数据生成原理和似然推断程序。“Egger”代表水平多效性的假设，参数γ将Egger回归假设有效推广到到具备LD相关结构的工具变量和多个结局。
当所述分析方法对汇总统计量数据进行分析时，模型可表示如下：
Figure BDA0002709151730000031
Figure BDA0002709151730000032
其中，
Figure BDA0002709151730000034
是一个p维向量，表示SNPs对基因表达的边际效应估计；
Figure BDA0002709151730000035
是一个k×p维矩阵，表示 SNPs对k个结局的边际效应估计；ex是服从多元正态分布
Figure BDA0002709151730000036
的p维误差向量； Ey是一个k×p维误差矩阵，服从矩阵正态分布MN(0，Ω，∑2/(n2-1)),Ω是k×k维的行协方差矩阵，∑2/(n2-1)是p×p维的列协方差矩阵。调整了上述用于个体水平数据的PX-EM算法，以实现汇总统计量moPMR-Egger模型的极大似然推断。
本发明的第二个方面，提供一种基于孟德尔随机化的多性状全转录组关联分析系统，所述系统包括：
数据获取单元，用于获取相关数据信息并进行标准化处理；
模型构建单元，用于基于获取的数据进行多结局建模，在似然框架下同时估计和检验基因对多个结局性状的因果效应，同时能将潜在的水平多效性嵌入模型。
其中，当所述系统对个体水平数据进行分析时，
所述数据获取单元获取的相关数据信息包括但不限于基因表达的个体数据、表达数据对应的基因信息，全基因组关联分析的多维性状数据和全基因组关联分析对应的基因信息。
进一步的，上述相关数据信息均需要标准化处理。
所述模型构建单元中采用如下回归模型：
Figure BDA0002709151730000041
Figure BDA0002709151730000042
Figure BDA0002709151730000043
尽管上述三个公式针对的是两个独立的研究，但参数β将二者紧密的结合起来。式(2)-(3)可以合并为
Figure BDA0002709151730000044
其中，
Figure BDA0002709151730000045
遵循遗传学中常用的多基因建模假设，假设β的所有元素非零，且每个元素都服从正态分布
Figure BDA0002709151730000048
同时，遵循PMR-Egger假设，假设对于每个性状i,每个SNP的水平多效性相等。实际上，这种水平多效性假设在其他TWAS和稳健MR研究中得到了广泛认可和应用，它等价于
Figure BDA0002709151730000049
1p是一个元素全为1的p维向量,γ＝(γ1，…，γk)T。
当所述系统对汇总统计量数据进行分析时，
所述数据获取单元获取的相关数据信息包括但不限于基因表达的汇总数据，全基因组关联分析的多维性状汇总数据，基因表达遗传变异的连锁不平衡矩阵，全基因组关联分析遗传变异的连锁不平衡矩阵和纳入分析的性状间的相关矩阵。
所述模型构建单元中采用如下模型：
Figure BDA0002709151730000046
Figure BDA0002709151730000047
其中，
Figure BDA00027091517300000410
是一个p维向量，表示SNPs对基因表达的边际效应估计；
Figure BDA00027091517300000412
是一个k×p维矩阵，表示SNPs对k个结局的边际效应估计；ex是服从多元正态分布
Figure BDA00027091517300000411
的p维误差向量； Ey是一个k×p维误差矩阵，服从矩阵正态分布MN(0，Ω，∑2/(n2-1)),Ω是k×k维的行协方差矩阵，∑2/(n2-1)是p×p维的列协方差矩阵。调整了上述用于个体水平数据的PX-EM算法，以实现汇总统计量moPMR-Egger模型的极大似然推断。
本发明的第三个方面，提供一种电子设备，包括存储器和处理器以及存储在存储器上并在处理器上运行的计算机指令，所述计算机指令被处理器运行时，完成上述基于孟德尔随机化的多性状全转录组关联分析方法所进行的步骤。
本发明的第四个方面，提供一种计算机可读存储介质，用于存储计算机指令，所述计算机指令被处理器执行时，完成上述基于孟德尔随机化的多性状全转录组关联分析方法所进行的步骤。
本发明的第五个方面，提供上述方法、系统、电子设备和/或计算机可读存储介质在如下任意一种或多种中的应用：
1)复杂疾病内在分子遗传机制研究；
2)疾病易感性研究；
3)药物反应研究；
4)药物靶标研究；
5)药物筛选研究。
上述一个或多个技术方案的有益技术效果：
上述技术方案提供一个多性状TWAS联合分析方法，称之为moPMR-Egger。moPMR-Egger一次检验一个基因，将存在特定基因包含的连锁不平衡顺式SNPs作为工具变量，实现对多个性状因果效应的联合检验。同时，moPMR-Egger能够检验和控制工具变量潜在的水平多效性，从而能够最大限度提高检验效能且避免了TWAS中的错误关联。大量统计模拟证实moPMR-Egger对因果效应的检验具有稳定的一类错误，且比现有TWAS分析方法具有更高的检验效能。进一步，使用moPMR-Egger 分析了英国生物样本库中五个类别的11个性状，结果显示，moPMR-Egger比单变量TWAS方法多发现了约13.15％的基因关联，并清晰揭示了收缩压和舒张压的不同生物调节机制。因此具有良好的实际应用之价值。
附图说明
构成本发明的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。
图1为本发明moPMR-Egger模型示意图以及在不同情形下moPMR-Egger与单性状方法的检验效能对比。(A)moPMR-Egger可以应用于TWAS研究，在混杂因素(U，未展示)存在的情况下，以cis-SNPs (z)作为工具变量，来估计基因的表达(x)对多个性状(y1,…,yk)的因果效应。moPMR-Egger依赖于一个联合似然框架，可以矫正水平多效性(y1,…,yk)和多性状之间的相关性。(B-G)基于moPMR-Egger的多性状联合模型在不同场景下都优于单性状方法，比如每次只能分析一种性状的PMR-Egger方法。对两个性状(y1,y2)的简单示例性模拟可以概括为以下六种情形。(B)基因影响两个相关性状，其对这两个性状的因果效应与这两个性状的关联方向相反。(C)基因影响两个不相关的性状。(D)基因影响两个相关性状，其对这两个性状的因果效应与这两个性状的关联方向相同。(E)基因只影响两个相关性状中的一个。(F)基因只影响两个不相关性状中的一个。(G)基因不影响任何性状。B-F中y轴表示检验效能，G中y轴表示I类错误。
图2为本发明不同模拟情形下因果效应检验的I类错误和检验效能。(A,B)在零假设成立时，水平多效性存在或不存在时因果效应检验的-log10p值的分位数图。根据不同水平多效性效应大小构建不同零假设的模拟：(A)γ＝(0，0，0，0)T(B)γj从(0，1×10-4，5×10-4，1×10-3，2×10-3)中随机选取， j＝1，2，3，4.在一定的水平多效性范围内，moPMR-Egger具有较好的因果效应检验的I类错误。(C-F)在异质性因果效应设定下，基于Bonferroni多重校正阈值因果效应的检验效能为y轴，第一个性状的 PVEzy为x轴，其它四个性状的PVEzy是第一个性状PVEzy的15％，50％和85％。比较方法包括 moPMR-Egger、PMR-Egger、PrediXcan和TWAS。不同的线条符号代表四种性状是否相关(C)、因果效应的方向(D，E，F)。在不存在水平多效性效应的情况下，对受影响的性状数量从1到4进行模拟(C-F).
图3为本发明不同模拟情形下水平多效性检验的I类错误和检验效能。(A，B)在零假设下，因果效应存在或者不存在时检验水平多效性的-log10p值的分位数图。根据不同因果效应值构建不同零假设的模拟：(A)PVEzy＝(0，0，0，0)T；(B)PVEzy，j从(0，0.005，0.01，0.015，0.02)中随机选取，j＝1，2，3，4. (A)γ＝(0，0，0，0)T(B)γj从(0，1×10-4，5×10-4，1×10-3，2×10-3)中随机选取，j＝1，2，3，4.在一定的因果效应范围内，moPMR-Egger具有较好的水平多效性检验的I类错误。(C-F)在不同的模拟情形下，检验水平多效性效应的moPMR-Egger和PMR-Egger的检验效能。基于Bonferroni多重校正阈值的水平多效性检验效能为y轴，PVEzy为x轴。模拟在相关性状(D，F)，独立性状(C，E)，水平多效性为γ＝(0，1×10-4，5×10-4，1×10-3)T(C，D)，水平多效性为 (1×10-4，5×10-4，1×10-3，2×10-3)T(E，F)下进行。
图4为本发明moPMR-Egger和PMR-Egger方法对UK Biobank中五个性状类别的TWAS分析结果。检验不同性状因果效应值的-log10p值分位数图：(A)血压性状类；(B)身体测量指标类；(C)血细胞计数类；(D)白细胞指标类；(E)红细胞指标类。检验不同性状水平多效性效应的-log10p值分位数图： (F)血压性状类；(G)身体测量指标类；(H)血细胞计数类；(I)白细胞指标类；(J)红细胞指标类。比较方法有moPMR-Egger、单性状PMR-Egger(不同性状表示为不同颜色)。SBP，收缩压；DBP，舒张压； BMI，身体质量指数；FVC，最大肺活量；RBC，红细胞计数。
图5为本发明moPMR-Egger和PMR-Egger对血压性状类的基因进行GO功能和KEGG通路富集分析。点图展示了对SBP和DBP识别出具有相同因果效应方向的基因(A)和具有相反因果效应方向的基因(C)前十个富集GO BP、CC和MF条目的结果。点图展示了对SBP和DBP识别出具有相同因果效应方向的基因(B)和具有相反因果效应方向的基因(D)前十个富集KEGG通路的结果。点的颜色表示基于q值富集分析的统计显著性，点的大小表示每项包含基因的比例。SBP，收缩压；DBP，舒张压；GO,基因本体；BP,生化过程；CC,细胞成分；MF,分子功能；KEGG,京都基因和基因组百科全书。
具体实施方式
应该指出，以下详细说明都是例示性的，旨在对本发明提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本发明所属技术领域的普通技术人员通常理解的相同含义。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本发明的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。应理解，本发明的保护范围不局限于下述特定的具体实施方案；还应当理解，本发明实施例中使用的术语是为了描述特定的具体实施方案，而不是为了限制本发明的保护范围。下列具体实施方式中如果未注明具体条件的实验方法，通常按照本领域技术内的分子生物学的常规方法和条件，这种技术和条件在文献中有完整解释。
以下通过实施例对本发明做进一步解释说明，但不构成对本发明的限制。应理解这些实施例仅用于说明本发明而不用于限制本发明的范围。
实施例
材料和方法
moPMR-Egger概述
moPMR-Egger依托两样本MR框架，该框架在没有样本重叠的两项独立研究中测量暴露和结局，旨在估计和检验暴露对结局的因果效应。在TWAS框架下，暴露是在基因表达研究中测量的基因表达水平，而结局是在GWAS中测量的多个相关性状。moPMR-Egger一次检验一个基因，并将 cis-SNPs作为工具变量，能够在似然框架下同时估计和检验基因对多个结局性状的因果效应，同时能够将潜在的水平多效性嵌入模型。图1A给出了moPMR-Egger模型的示意图。
具体模型的数学表达式可以表述如下，记x为n1维向量，表示n1个个体在基因表达研究中的基因表达水平，Zx为n1×p维的cis-SNPs基因型矩阵，表示在基因表达研究中对应于该基因的p个工具变量，Y为一个k×n2维的矩阵，表示GWAS研究中n2个个体的k个相关性状的测量值，ZY为对应于p个相同工具变量的GWAS研究中n2×p维的基因型矩阵。不失一般性，假设x和YT,Zx,ZY的每一列都实行了标准化。考虑如下三个回归模型：
Figure BDA0002709151730000071
Figure BDA0002709151730000072
Figure BDA0002709151730000073
式(1)刻画了基因表达研究中特定基因的cis-SNPs与其表达水平的关系，式(2)-(3)描述了GWAS 研究中预测的基因表达和多维相关性状之间的关系。μx是一个标量，
Figure BDA0002709151730000082
是一个k维向量，二者都代表截距项；
Figure BDA0002709151730000083
是一个n2维向量，代表GWAS研究中不可观测的基因表达预测值；β为p维向量，刻画了cis-SNPs对基因表达的效应大小；Γ是一个k×p的水平多效性矩阵，代表p个工具变量对k个结局变量的水平多效性；εx是n1维的残差向量，其中的每个元素独立且服从正态分布
Figure BDA0002709151730000084
Figure BDA0002709151730000085
是n2维的残差向量，其中的每个元素独立且同服从正态分布
Figure BDA0002709151730000086
Figure BDA0002709151730000087
是一个k×n2的残差矩阵，它的每一列都服从一个多元正态分布MVNk(0，Ω),Ω是一个k×k的协方差矩阵。尽管上述三个公式针对的是两个独立的研究，但参数β将二者紧密的结合起来。式(2)-(3)可以合并为
Figure BDA0002709151730000081
其中
Figure BDA0002709151730000088
需要特别说明的是，组学数据研究中p往往远大于n1,因此需要对β做一定的建模假设以保证模型可识别性。类似的，也需要对Γ做额外的建模假设，否则，式(4)中定义的两个工具变量效应项：垂直多效性α(ZYβ)T和水平多效性
Figure BDA0002709151730000089
将存在无法识别的可能性。这里，遵循遗传学中常用的多基因建模假设，假设β的所有元素非零，且每个元素都服从正态分布
Figure BDA00027091517300000811
同时，遵循PMR-Egger 假设，假设对于每个性状i,每个SNP的水平多效性相等。实际上，这种水平多效性假设在其他TWAS 和稳健MR研究中得到了广泛认可和应用，它等价于
Figure BDA00027091517300000810
1p是一个元素全为1的p维向量, γ＝(γ1，…，γk)T。
在上述联合模型中，所关心的参数是基因在多个结局性状上的因果效应α。α的因果效应识别和可解释性可以基于因果推断中的决策理论加以证明。和传统的MR分析一样，α的因果解释至少需要两个假设成立：(1)工具变量与基因表达有关；(2)工具变量与基因表达和每个结局都有关的任何其他混杂因素无关。因为对水平多效性Γ做了明确的模型假设，moPMR-Egger不再要求满足MR的一般排除限制条件(即，工具变量只能通过基因表达途径影响每个结局)。然而，对Γ的模型假设需要满足工具变量强度与直接效应独立的假设(Instrument Strength Independent of Direct Effect,InSIDE), 这个假设要求工具变量-基因表达效应和工具变量-结局效应是相互独立的，有时称为弱排除限制条件。因为混杂因素往往是未知的，这些假设有许多在实践中是无法检验的，本研究所使用的术语“因果效应”，只是为了强调α的估计比Y对
Figure BDA00027091517300000812
的普通回归获得的效应值估计更可信。
在上述模型中，感兴趣的是在水平多效性γ存在时，估计因果效应α并检验原假设H0：α＝0,拒绝原假设α＝0意味着基因表达至少对一个性状有非零的因果效应。另外，还可对水平多效性γ进行估计和检验H0：γ＝0,拒绝原假设γ＝0说明cis-SNPs至少对一个性状有非零的水平多效性。在统计推断层面，发展了期望最大化(EM)算法的参数扩展版本(PX-EM),通过最大化基于式(1)和(4) 定义的联合似然来进行参数推断。可以获得最大似然联合模型以及α和γ最大似然估计。此外，进一步将该EM算法应用于两个简化模型，一个模型中不包含α，另一个模型中不包含γ，分别获得对应的最大似然。最后，将从联合模型得到的最大似然函数与从两个简化模型得到的最大似然函数进行对比，实现对H0：α＝0或H0：γ＝0的似然比检验。
将此模型和算法一起称为具有Egger假设的多结局概率孟德尔随机化(moPMR-Egger)方法。其中，“mo”指的是对多个结局(multipleoutcome)的建模。“概率”是指moPMR-Egger模型完全模仿了 TWAS背后的数据生成原理和和似然推断程序。“Egger”代表水平多效性的假设，参数γ将Egger回归假设有效推广到到具备LD相关结构的工具变量和多个结局。
基于汇总数据(Summarystatistics)的moPMR-Egger模型
考虑到对汇总数据的分析已经成为科学热点，进一步发展了基于汇总数据的moPMR-Egger模型。记∑1为基因表达研究中特定基因cis-SNPs间的SNP-SNP相关矩阵(即，LD矩阵)。∑2为GWAS 数据中对应的相关矩阵。这两个矩阵都是p×p维的半正定矩阵。通常，∑1，∑2都可以用LD参考面板数据中相应种族的个体数据进行估算(例如，千人基因组计划)。基于汇总数据的moPMR-Egger 模型可表示如下
Figure BDA0002709151730000091
Figure BDA0002709151730000092
其中，
Figure BDA0002709151730000093
是一个p维向量，表示SNPs对基因表达的边际效应估计；
Figure BDA0002709151730000094
是一个k×p维矩阵，表示 SNPs对k个结局的边际效应估计；ex是服从多元正态分布
Figure BDA0002709151730000095
的p维误差向量； Ey是一个k×p维误差矩阵，服从矩阵正态分布MN(0，Ω，∑2/(n2-1)),Ω是k×k维的行协方差矩阵，∑2/(n2-1)是p×p维的列协方差矩阵。调整了上述用于个体水平数据的PX-EM算法，以实现汇总统计量moPMR-Egger模型的极大似然推断。
模拟
进行了跨基因模拟设计来评估moPMR-Egger的性能，并将其与现有TWAS方法进行比较。先从GEUVADIS数据中随机选择10,000个基因，并提取这些基因对应的cis-SNPs，基因cis-SNPs的中位数为576个(最小值为11；最大值为7409)。对每个特定基因，分别利用基因表达研究和GWAS 研究获得的基因型数据，依次模拟其表达水平和结局性状。具体来讲，首先从GEUVADIS数据中获得了每个基因的p个cis-SNPs的基因型，得到了标准化的基因型矩阵Zx之后，从正态分布N(0，PVEzx/p)中模拟SNPs对基因表达的效应β，其中，PVEzx表示基因表达变异被遗传效应解释的比例。对所有cis-SNPs的遗传效应求和为Zxβ。进一步，从正态分布N(0，1-PVEzx)中模拟了残差εx，最后将遗传效应和残差相加，得到基因表达水平。
接下来，在Kaiser Permanente/UCSF成人健康和老龄化遗传流行病学(KaiserPermanente/UCSF Genetic Epidemiology Research Study on Adult Health andAging，GERA)数据中随机从对照人群中选择2000个个体相同的p个cis-SNPs，获得标准化的基因型矩阵Zy，针对四个相关结局性状使用基因表达数据中相同的SNP效应大小β，设置四个因果效应
Figure BDA0002709151730000101
其中，PVEzy是一个4维向量，它的每个元素代表由SNPs的因果效应所解释的对应表型变异的比例。然后，模拟了来自多元正态分布MVN(0，Ω)的每个个体的四种表型的残差。为了使模拟更加贴近实际，借助NFBC1966数据，计算的四个血脂性状间的相关矩阵(总胆固醇(TC)、低密度脂蛋白(LDL)、高密度脂蛋白(HDL)和甘油三酯(TG))作为Ω.最后，求出水平多效性效应、垂直多效性效应和残差的和，得到了四个模拟的相关结局性状，对四个性状水平多效性效应γ的设置，将在下一段中详述。
最基本的模拟情形设置为：n1＝465，n2＝2000，PVEzx＝10％，PVEzy＝(0，0，0，0)T，Γkj＝0(j＝1，…，p；k＝1，…，4).在该基本设置的基础上，每次改变一个参数，以探索各种参数对模型的影响。对于PVEzx，设置为1％，5％或10％，这三个值都接近基因表达遗传度估计的中位数。对于β，考察了偏离多基因假设的不同分布，随机选择1个SNP，1％或10％的SNPs对基因表达有非零效应，并且从正态分布中模拟非零效应来解释一个固定PVEzx.对于水平多效性效应Γ，为避免人为设定每个γj(j＝1，2，3，4)的取值对结果的影响，随机从0，1×10-4，5×10-4， 1×10-3，2×10-3选择一个取值，，这些取值是依托真实数据估计出来的水平多效性效应的0％、 50％、70％、90％、95％分位数。同时，从0％、10％、30％、50％，100％中随机选择四个比例值，作为四个性状上具有非零水平多效性的SNPs的比例。进一步，还研究了定向多效性设置(水平多效性效应负与正的SNPs比值为0∶10)、近似定向多效性设置(1∶9或3∶7)和平衡多效性设置(5∶5).对于Ω，除了使用从真实NFBC数据估计的Ω之外，还研究了Ω为单位矩阵及其Ω的每个元素有指数协方差结构Ωi，j＝ρ|i-j|的情形，ρ取0.5，0.7，0.9.此外，研究了β之间相关的情形，模拟SNPs对基因表达的效应服从协方差矩阵为w∑的多元正态分布，∑是SNPs的LD矩阵，w是一个标量，以确保PVEzx＝10％.对于PVEzy，除了上述零因果效应的基本设置外，还考察了同质因果效应设置和异质因果效应设置两种情形。在同质因果效应设置中，PVEzy的所有元素都具有相同的值v，v为0.5％/4，1％/4，1.5％/4，2％/4.在异质的因果效应的设置中，第一个性状因果效应为v，其他三个性状的因果效应为0.15v，0.85v和0.50v，v分别是0.5％，1％，1.5％或2％.在效应的方向上，分别在同质和异质因果效应设置中随机选取了基因表达所影响的四个相关性状中的一、二、三、四个性状。当基因表达影响两个性状时，对其中一个性状的影响与对另一个性状的影响效应符号是相反的。当基因表达影响三个性状时，其对随机选择的两个性状的影响与对第三个性状的影响是效应符号相反的。当基因表达影响四个性状时，其对随机选择的两种或三种性状的影响与对另外两种或一种性状的影响效应符号是相反的。
对于I型错误检验，对上述每个设置情形执行了10,000次模拟。对于检验效能计算，使用了两种不同的计算方法。在第一种方法中，进行了1,000次备择假设成立时的模拟，并基于Bonferroni 多重校正的p值阈值比较了检验效能(多变量方法为0.05/1000，单变量方法为0.05/4000)。在第二种方法中，为每个模拟情形执行了100个备择假设成立的模拟和900个原假设成立的模拟，并在错误发现率(false discovery rate，FDR)为0.05的标准下计算检验效能。同时，基于一种上述模拟情形，也对基于汇总数据的moPMR-Egger模型进行了评估(详情见讨论)。
为了更好地说明moPMR-Egger模型，考虑了水平多效性存在时两个性状的不同基因因果关联模式，设置n1＝465，n2＝2000PVEzx＝10％，γ＝(0.0001，0.0005)T，基因表达影响两种性状： PVEzy＝(0.01，0.01)T，基因表达只影响一种性状：PVEzy＝(0.01，0)T，基因表达对这两种性状都没有影响：PVEzy＝(0，0)T，性状间的相关性ρ＝0.66，性状独立时ρ＝0.需要指出的是，ρ＝0.66对应于英国生物样本库(UK Biobank)中估算的收缩压(systemicblood pressure，SBP)与舒张压(diastolic blood pressure，DBP)的相关系数。总体而言，研究了水平多效性存在时，考虑三种因果效应值设置和两个性状相关性设置的六种基因因果关联模式。考虑到单性状TWAS方法中，只有 PMR-Egger可以解释水平多效性，因此，对于每一种情况，比较了moPMR-Egger与PMR-Egger的性能。对于I类错误，进行了100次1，000个零假设下的模拟，计算了moPMR-Egger在100次重复中p值小于0.05/1000的频率，PMR-Egger100次重复中p值小于0.05/2/1000的频率。对于检验效能，进行了1,000次备择假设成立时的模拟，并基于Bonferroni多重校正的p值阈值比较了检验效能。
与moPMR-Egger相比较的方法
目前还没有多性状TWAS分析方法，因此，比较了三种常用的单性状TWAS方法，包括：(1) PMR-Egger，使用多个相关的SNP工具变量对水平多效性进行检验和控制。PMR-Egger模型使用了所有cis-SNPs作为工具变量，并使用R包PMR中的PMR_individual函数来检验因果效应和水平多效性效应。(2)PrediXcan，使用多个相关的SNP工具变量，但不控制水平多效性。PrediXcan使用了所有cis-SNPs，并使用R包glmnet执行ElasticNet来获得cis-SNPs对基因表达的系数估计。(3)TWAS，采用多个相关的SNP工具变量，但不控制水平多效性。TWAS使用了所有cis-SNPs，并使用GEMMA 软件中BSLMM方法来获得cis-SNPs对基因表达的系数估计。总体而言，这些方法都属于两样本设计，都能计算p值来检验因果效应αi，i＝1，…，4.然而，这三种方法的先验假设有所不同：PrediXcan 依赖于弹性网假设；TWAS依赖于贝叶斯稀疏线性混合模型假设；而PMR-Egger依赖于正态先验假设。在推断上，PrediXcan和TWAS依赖于两阶段回归，而PMR-Egger则是基于最大似然。在比较过程中，采用单性状TWAS方法一次分析一个性状，从而根据Bonferroni多重校正的p值阈值判定基因关联显著性，该阈值调整了检验的基因数量和性状数量。
此外，考虑到Bonferroni矫正或许会降低单性状TWAS方法的检验效能，进一步提出了单性状 TWAS分析的最小p值策略，对于每种单性状方法和一个特定基因，获得多个性状的最小p值，将其作为该基因的因果关联证据。因为获得最小p值的统计量分布并不容易，所以设定FDR水平为 0.05，比较了moPMR-Egger和单性状最小p值法的检验效能。对于实际数据分析，首先对个体进行 10次随机置换，同时保持性状之间的相关性。然后，将每种方法应用于置换后的数据，以获得最小 p值的统计量分布，据此计算FDR，并基于FDR阈值0.05计算因果关联基因的个数。
实际数据应用
通过整合GEUVADIS基因表达数据和UKBiobank的GWASs，应用moPMR-Egger进行多性状 TWAS分析。获得了GEUVADIS数据作为基因表达数据，并从UKBiobank中提取11个性状，这些性状在先前研究中均有报道，其遗传度均大于0.2。将这些结局性状分为以下五类：血压(SBP和 DBP)、身体测量指标(身高、体重指数、肺活量(FVC)和FEV1-FVC比率)、血细胞计数(血小板计数、红细胞计数、嗜酸性粒细胞数和白细胞计数)、白细胞指标(嗜酸性粒细胞计数和白细胞计数) 和红细胞指标(红细胞计数和红细胞分布宽度(RDW))。可以发现，白细胞指标类是血细胞计数类的一个子集。对于每一个类别，使用moPMR-Egger和其所比较的方法来分析这五个类别中的性状。下面详细介绍GEUVADIS和UKBiobank的数据处理流程和步骤。
GEUVADIS数据包含了465个个体的基因表达测量值，这些个体来自5个不同的群体，包括 CEPH(CEU)、芬兰人(FIN)、英国人(GBR)、托斯卡尼人(TSI)和约鲁巴人(YRI)。在表达数据中，只关注在GENCODE(release 12)中注释的蛋白编码基因和lincRNAs.在这些基因中，剔除了至少在一半的个体中没有计数的低表达基因，最终得到15,810个基因。根据之前的研究，进行了PEER归一化以消除混杂效应和多余的变异。随后，为了消除剩余的群体分层，对每个群体中个体的基因表达测量值进行分位数标准化，以达到标准正态分布，然后对来自所有五个群体的个体的基因表达测量值进行分位数标准化，以达到标准正态分布。除了表达数据外，GEUVADIS的所有个体也在千人基因组计划中进行了基因型测序。千人基因组项目3期中获得基因型数据。过滤掉了具有哈迪-温伯格平衡(HWE)p值<10-4，基因型召唤率<95％，或次要等位基因频率(MAF)<0.01的SNPs,总共保留了 7,072,917个SNPs进行分析。
UK Biobank数据包括487,298个个体和92,693,895个填补的SNPs。遵循Neale实验室 (https://github.com/nealelab/uk_biobank_gwas/tree/master/impuing-v2-gwas)相同的样本质量控制程序，保留了总计337,129个欧洲血统的个体。过滤了HWE p值<10-7、基因型召唤率<95％或MAF<0.001 的SNPs,总共得到13,876,958个SNPs。对每个性状依次进行性别和前10个主成分(PCs)的标准化表型回归得到残差，将残差标准化，使其均值为0，标准差为1，最后利用这些标准化的残差进行TWAS 分析。
将GEUVADIS数据与UK Biobank中的GWASs数据进行整合，进行TWAS分析。对于GEUVADIS数据中的每个特定基因，依次提取了位于转录起始位点(TSS)上游100kb或转录末端位点(TES)下游100kb区间内的cis-SNPs。将GEUVADIS中的这些SNPs与UK Biobank中获得的SNPs 取交集，得到共同的SNPs集。GEUVADIS和UKBiobank重叠的所有基因的cis-SNPs的平均个数为497(中位数为443,最小值为1,最大值为8085)。然后，对于每一对基因(来自GEUVADIS) 和性状类别(来自GWAS)，依次运行多性状moPMR-Egger方法和单性状PMR-Egger方法来检验基因表达和每个类别中多个性状之间的因果关系。为了比较moPMR-Egger和PMR-Egger，依据 Bonferroni多重校正阈值作为显著性的标准：moPMR-Egger为0.05/15,810,PMR-Egger为 0.05/(15,810*k)，k为每个类别中结局性状的数量。
在五个类别中，重点深入分析了血压这一类别。血压类包含SBP和DBP两个性状，两者之间的相关性估计为0.66。血压是一个复杂的性状，遗传度约在30-50％范围内。许多大规模的GWASs 研究了血压的遗传机制。重要的是，血压升高是导致中风和冠状动脉疾病风险的强有力且可改变的风险因素，而中风和冠状动脉疾病是全球死亡和疾病的主要原因。对于血压类别，除了moPMR-Egger 和PMR-Egger，还应用了PMR-Egger、TWAS和PrediXcan的最小p值方法。对于这些最小p值方法和每个特定基因，依次分析类别中的每个性状，并获得这些性状间的最小p值，作为该特定基因因果关联的证据。对个体进行了10次随机置换，借此获得最小p值的统计量分布，并据此计算出 FDR。在moPMR-Egger分析中，依据基因对两种性状的效应大小和方向，将具有统计显著性的基因分为两组：在两种性状上具有相同效应符号的基因和在两种性状上具有相反效应符号的基因，针对这两组不同的基因进一步通过R包clusterprofiler使用基因本体(GO)和京都基因和基因组百科全书 (KEGG)进行富集分析，并根据q值阈值0.05来声明基因集的显著性。
结果
moPMR-Egger模型概览
简单来讲，moPMR-Egger一次检验一个基因，并估计和检验某特定基因对多个结局性状的因果效应。与现有的TWAS方法不同，moPMR-Egger采用了具有LD相关关系的多个SNP工具变量，在最大似然框架下实现TWAS的统计推断，能够检验和控制TWAS中广泛存在的水平多效性效应，同时联合分析多个相关性状，借此提高了TWAS分析的能力。首先总结了图1B-G中两个示范性状的moPMR-Egger和PMR-Egger的比较结果。当基因影响两个性状，moPMR-Egger在基因效应与性状相关性方向相反时较高的检验效能(图1B)，紧随其后的是两个性状是独立的情形(图1C)、基因效应与性状相关性方向相同的情形(图1D)。而PMR-Egger的检验效能在这三种情况下几乎相同 (图1B-D)。出现这种结果的潜在可能的原因是，当二者方向相同时，moPMR-Egger很难区分因果效应和性状间相关性。当基因只影响一种性状时，通过利用两种性状之间的相关性，moPMR-Egger 的检验效能仍然高于PMR-Egger(图1E)，而当两种性状相互独立时，moPMR-Egger与PMR-Egger 的检验效能几乎相同(图1F)。总的来说，相对于PMR-Egger,moPMR-Egger的检验效能在大多数情况下要高得多。此外，当基因对这两种性状都没有因果效应时，moPMR-Egger得到了校准的I型误差控制，而PMR-Egger在0.05的期望值下略微膨胀(图1G)。
模拟：检验和估计因果效应
通过全面的模拟来检验moPMR-Egger的有效性，并将其与TWAS的现有方法进行了比较。材料和方法中提供了具体的模拟。简单地说，依靠真实的基因型数据，模拟基因表达以及四个结局性状，使用moPMR-Egger来分析这四个性状，一次一个基因。在总共287个模拟情形下(因果效果检验的25个零假设情形和152个备择假设成立的情形；水平多效性效果检验的22个零假设情形和 88个备择假设情形)考虑了moPMR-Egger因果效应检验和水平多效性检验的一类错误控制和检验效能。
第一组模拟针对因果效应检验，首先在零假设下检验了moPMR-Egger的I型错误控制。结果发现，无论是在存在水平多效性效应(图2A)还是不存在水平多效性效应(图2B)情况下，moPMR-Egger 都较好的控制了I型错误，而且无论基因表达的遗传结构是稀疏的还是多基因的，无论基因表达遗传度，无论SNP对基因表达的影响是否与LD相关，无论多个性状是否相关，moPMR-Egger始终能够较好的控制I型错误。需要注意的是，当基因表达遗传度较低时，moPMR-Egger的p值略微保守，这可能是由于联合似然提供的信息有限，渐近正态分布不能很好地近似。
moPMR-Egger对水平多效性做了一个相对强的Egger假设，该假设要求对于给定的性状，所有 SNP都具有相同的水平多效性效应。为了检验这种假设的稳健性，随机选择固定比例的SNP，而不是全部SNP，具有水平多效性效应。结果发现，无论水平多效性SNP的稀疏性如何，moPMR-Egger 都能很好的控制I类错误。除了正向和负向效应SNP的比例设置为0:10外，还通过随机分配相应比例的SNP，检验了两种近似方向多效性设置(1:9或3:7)和一种平衡设置(5:5)。结果表明，当每个性状的水平多效性效应较小或中等时，moPMR-Egger在近似方向性多效性设置或平衡设置中较好的控制I类错误。然而，当某一性状的水平多效性效应很大时，正如预期的那样，moPMR-Egger出现了 I类错误膨胀问题。总体而言，moPMR-Egger在检验因果效应零假设下表现良好。
接下来，评估了moPMR-Egger识别非零因果效应的检验效能，并将其与三种单变量方法进行比较。在大多数模拟情形下，moPMR-Egger比单变量TWAS方法更强大。在总共152个模拟情形中，与PMR-Egger、PrediXcan和TWAS相比，moPMR-Egger获得了平均53.12％、42.40％和36.79％的检验效能增益。不管因果效应在不同的性状间是同质的还是异质的(图2C-F)，无论基因影响一个性状或多个性状，不管性状是如何相互关联，存在或不存在水平多效性，moPMR-Egger始终具有大量的检验效能增益。唯一的例外是当四个性状是相互独立情形，且基因只与四个性状中的一个相关(图2C)。在这种情况下，moPMR-Egger的低检验效能可能是因为moPMR-Egger使用了额外的参数来对多个相关性状的非零随机效应进行建模，从而导致自由度的损失，进而降低了检验效能。然而，即使基因只与四个性状中的一个相关，当性状之间相互关联时，moPMR-Egger仍然比其他方法具有更大的检验效能(图2C)。值得注意的是，moPMR-Egger和PMR-Egger可以控制水平多效性，而PrediXcan和TWAS则不能控制水平多效性。因此，在水平多效性存在的情况下，PrediXcan 和TWAS无法控制I类错误，p值被夸大。如果将所有的p值进行校准，并依赖于一个名义p值阈值，而不是正确的I类错误阈值，那么PrediXcan和TWAS似乎比PMR-Egger具有更高的检验效能。尽管如此，基于名义p值阈值的PrediXcan和TWAS的检验效能始终低于moPMR-Egger，这进一步突出了多性状联合建模的重要性。
不同方法的检验效能和moPMR-Egger带来的检验效能增益取决于几个重要参数。首先，无论水平多效性是否存在，所有方法的检验效能都随着因果效应的减少而降低。然而，在所有方法中， moPMR-Egger的检验效能降低趋势是最不明显的，表明moPMR-Egger具有一定程度的稳健性能。第二，moPMR-Egger采用了多性状联合建模，并将性状间的相关关系考虑在内。因此，当性状之间的相关性增加时，moPMR-Egger的检验效能增加，而其他方法的检验效能基本保持不变。第三，对多性状联合建模所带来的检验效能增益取决于基因表达对性状的因果效应是否与性状相关性方向相同。当两种性状之间存在正相关关系，对两个性状的因果效应相反时的检验效能要高于对两个性状因果效应相同时候的检验效能；反之亦然(图2D)。最后，由于moPMR-Egger考虑了SNPs对多个性状的水平多效性效应。因此，无论受水平多效性影响的性状数量或水平多效性效应的大小如何变化，moPMR-Egger仍然有效，而其他方法则无效。
在上述模拟中，主要比较了moPMR-Egger与现有三种单性状的TWAS方法。进一步，还比较了moPMR-Egger和这些单性状TWAS方法最小p值策略。首先获得给定基因关联证据的不同性状间的最小p值，并在FDR 0.05的水准下，计算了不同方法的检验效能。与先前的模拟结果一致，无论因果效应对于不同性状是同质的还是异质的，无论水平多效性是否存在，moPMR-Egger仍然比 PMR-Egger,PrediXcan,TWAS的检验效能高。不同于上述Bonferroni调整阈值的结果，在FDR标准下，当存在水平多效性效应时，PrediXcan和TWAS均显示较低的检验效能。
最后，无论水平多效性是否存在，moPMR-Egger在零假设和各种备择假设下，均给出了因果效应的准确估计。
模拟：检验和估计水平多效性效应
第二组模拟主要是针对水平多效性效应的检验。与一般的TWAS/MR方法相比，moPMR-Egger 能够检验SNPs是否对任一结局性状表现出非零的水平多效性效应。
首先在零假设下检验了moPMR-Egger对水平多效性效应的一类错误控制情况，在零假设条件下，即，任何一个遗传变异对4个性状都不存在水平多效性效应。发现无论因果效应是否存在(图 3A,B)，无论多个性状之间的相关性大小，因果效应不存在时无论基因表达的基因结构是稀疏或多基因，moPMR-Egger对水平多效性的检验均能很好的控制I类错误。然而，当模型假设以多种方式不成立时，例如，当基因表达的基因结构是稀疏的，当基因影响两种以上性状时，moPMR-Egger 就会出现I类错误膨胀问题。总体来说，在零假设下，moPMR-Egger对水平多效性的检验均能较好的控制I类错误。
接下来，评估了moPMR-Egger识别非零水平多效性效应的效能。由于PMR-Egger是目前唯一可以为水平多效性效应提供检验的方法，为此，重点比较了moPMR-Egger和PMR-Egger的检验效能。模拟发现，PMR-Egger和moPMR-Egger的检验效能随着水平多效性的增加而增加，在一系列设置中，moPMR-Egger的表现优于PMR-Egger(图3C vs图3E,图3D vs图3F)。与上述因果效应的模拟一致，当四个性状彼此不相关时，moPMR-Egger与PMR-Egger具有类似的检验效能(图3C, E),然而，当性状间存在相关性时，moPMR-Egger优于PMR-Egger(图3D,F)。此外，无论是否存在因果效应，在一系列稀疏水平范围内，moPMR-Egger仍然比PMR-Egger更有效，这两种方法的检验效能均随着水平多效性稀疏程度的增加而减小。同样，与因果效应检验的结果一致，两种方法在多效性检验中，在没有定向水平多效性的情况下，效果都很差。进一步将moPMR-Egger与 PMR-Egger的最小p值方法进行了比较，得到了一致的结果。
最后，在存在定向水平多效性时，moPMR-Egger可以准确地估计水平多效性效应。然而，在没有定向多效性效应的情况下，多效性效应估计值是有偏的。
实际数据应用
将GEUVADIS的基因表达数据与来自UKBiobank五个类别的11个性状的GWAS数据进行整合(详见材料和方法部分)，应用moPMR-Egger进行了TWAS分析。GEUVADIS研究的基因表达数据共包括15,810个基因，UK Biobank的五个性状类别包括：血压类别、身体测量指标类别、血细胞计数类别、白细胞指标类别和红细胞指标类别。
展示了每个类别中检验每个基因对moPMR-Egger性状的因果效应的p值，以及每个类别中单性状PMR-Egger方法的p值(图4A-E)。PrediXcan和TWAS两种方法本质上不能控制水平多效性，且这两种方法往往出现p值膨胀问题，因此，本研究并未采用PrediXcan和TWAS方法分析，相反，将在下一段中采用这两种方法的最小p值策略进行深入分析。结果显示moPMR-Egger和PMR-Egger 的p值都表现良好，在全转录组Bonferroni多重校正阈值下，moPMR-Egger识别出的基因数目比 PMR-Egger多13.15％的基因(表1)，这突出了在TWAS中多性状联合分析的优势。PMR-Egger识别的基因(89.29％)大部分亦能被moPMR-Egger检测到，且moPMR-Egger方法的因果关联显著性增加。例如，TFRC基因被识别出与RDW相关(PMR-Egger p＝3.22×10-17)，但当所有红细胞指标类别的性状联合在一起时，关联证据明显变强(moPMR-Eggerp＝5.72×10-22)。TFRC编码经典转铁蛋白受体，参与细胞的铁吸收。TFRC中多个SNPs已被证实与GWAS研究中多种红细胞表型相关。这些相关的红细胞表型包括红细胞平均血红蛋白(mean corpuscular hemoglobin,MCH)和红细胞平均体积 (meancorpuscular volume,MCV),均与RDW密切相关。TFRC的变异可能正如在Tfrc缺乏的小鼠中观察到的那样，导致红细胞前体铁可用性降低，从而导致RDW测量的红细胞大小代偿性增加。更重要的是，moPMR-Egger识别了许多未能被PMR-Egger识别的基因。例如，在血压类中，只有 moPMR-Egger发现了EPHB4基因(moPMR-Egger p＝1.41×10-8,SBP PMR-Egger p＝0.98,DBP p＝5.33×10-4)。EPHB4编码Ephrin b型受体4，该受体与Ephrin-b2结合，在血管发育和血管生成中发挥重要作用。先前的研究表明，在小鼠中，Ephb4的缺失会导致低血压，这支持了Ephb4调节血压的因果作用。
随后，对血压类进行了深入的分析，它只包含两种性状(SBP和DBP)。除了与上述PMR-Egger 方法进行比较外，还将moPMR-Egger与PMR-Egger、PrediXcan、TWAS的最小p值方法进行了比较。计算了不同方法在FDR为0.05的标准下的结果。结果发现，moPMR-Egger比其他方法识别了更多的基因，再次验证了多性状联合建模的优势。moPMR-Egger、PMR-Egger、PrediXcan和TWAS 识别出的显著基因数量分别为765、691、90和83，这也说明基于似然的moPMR-Egger和PMR-Egger 方法比基于两阶段回归推断的PrediXcan和TWAS方法具有更高的效能。在moPMR-Egger鉴别出的显著基因中，约有40％(301/765)对这两个性状具有相反的因果效应方向。GO和KEGG通路富集分析(图5)显示，对两种性状具有相同因果效应方向的基因在抗原加工呈递(q＝1.96×10-4；图5A)、人类T细胞白血病病毒1感染(q＝0.03)、病毒性心肌炎(q＝0.03)和IgA产生的肠道免疫网络(q＝0.04；图 5B)通路显著富集。这种在炎症和诱导免疫反应途径中的富集支持了最近揭示的它们在调节血压中的作用。事实上，炎症激活先天和适应性免疫反应，导致血管系统、肾脏和交感神经系统(SNS)的改变，最终导致血压长期升高。相反，对两种性状具有相反因果效应方向的基因在COP9信号体(q＝0.02；图5C)和溶酶体(q＝5.14×10-4；图5D)通路中显著富集。COP9信号体和溶酶体途径的富集支持了它们在维持血管稳态和可塑性以及随后的血压调节中的关键功能。COP9信号体和溶酶体都控制着血管系统中的泛素化和随后的蛋白循环调节。血管中的蛋白质泛素化和循环决定了血管张力和硬度，这可能在细胞收缩和舒张的心脏周期中对SBP和DBP产生不同的影响。
接下来，关注水平多效性效应的检验。展示了5个类别两种方法(moPMR-Egger和PMR-Egger) 基因水平多效性检验的p值(图4F-J)。moPMR-Egger和PMR-Egger的p值均表现良好，moPMR-Egger 识别出的水平多效性基因总数比PMR-Egger多17.10％(表1)。具体来说，moPMR-Egger在血压类 中识别出56个基因，在身体测量类中识别出275个基因，在血细胞计数类中识别出383个基因，在 白细胞指标类中识别出146个基因，在红细胞指标类中识别出153个基因。相比之下，PMR-Egger 识别出的基因数分别是26、196、266、106和136。PMR-Egger发现的大部分基因(71.53％)也被 moPMR-Egger所包含。在moPMR-Egger识别的水平多效性基因中，36.16％具有显著的因果效应， 同时，具有显著的因果效应基因中，11.64％具有水平多效性效应。具有水平多效性效应的基因与具 有因果效应的基因之间明显的重叠，突出了在TWAS建模分析时同时考虑这两种效应的重要性。血 压类别的深入分析也显示，对DBP和SBP有显著因果效应的基因在抗原加工递呈和蛋白重折叠途 径中富集，而具有显著水平多效性效应的基因在各代谢过程和蛋白质丰富输出中富集。
最后，moPMR-Egger同时具备很高的计算效率，计算时间和物理内存需求与现有的TWAS 方法相似(表2)。
讨论
本发明提出了TWAS分析的moPMR-Egger方法，该方法扩展了单性状PMR-Egger到多性状联合分析。moPMR-Egger考虑了多个性状之间的相关结构，利用了所有LD相关的cis-SNPs，检验和控制水平多效性效应，并在似然框架下完成推断。在模拟和真实数据应用中，moPMR-Egger均能得到稳定的p值，且大大提高了现有单性状TWAS方法的检验效能。
moPMR-Egger的重要的建模假设是，SNPs在同一性状上表现出相同的水平多效性效应，该假设直接遵循了Egger回归假设，类似于通常用于罕见变异关联分析的负荷效应大小假设。与之前的研究一致，发现，在moPMR-Egger中采用的这种假设对于因果效应的估计和检验是相当稳健的。然而，必须充分认识到，这种等水平多效性假设在许多情况下可能过于严格。因此，尽管认为 moPMR-Egger是实现多性状联合分析的重要第一步，但在未来对水平多效性采用更现实的建模假设是必要的。
主要阐述了moPMR-Egger在分析个体水平数据方面的优势。moPMR-Egger可以很容易地扩展到对汇总统计量数据分析中(详见材料和方法部分)。具体来讲，基于汇总统计量的moPMR-Egger 需要基因表达和多结局性状的边际SNPs效应大小估计及其标准误，cis-SNPs间的LD相关矩阵，这个矩阵可以根据参考数据集构建。此外，还需要多个性状间相关矩阵的估计。利用每个性状的全基因组SNPs的边际Z得分统计量，可以很容易地估计性状间相关矩阵：因为两个性状之间的相关性约等于两个性状零假设下边际Z得分统计量之间相关性。因而，对于每个性状，可以选择p值大于预先设定的显著性阈值(例如10-5)的SNPs，使用这些边际Z得分统计量间的相关性作为性状间相关矩阵的估计。在模拟中验证了基于汇总统计量的moPMR-Egger基于个体水平数据的 moPMR-Egger的拟合结果基本一致，这大大增加了moPMR-Egger的适用范围。基于个体水平和汇总统计量水平的moPMR-Egger的详细算法见表3，4。
表1在Bonferroni多重校正阈值下，英国生物库中五类性状的显著基因数
Figure BDA0002709151730000191
表2在英国生物样本库五个类别中TWAS分析时每个基因的平均CPU时间和内存使用
Figure BDA0002709151730000201
计算是在Intel XeonE5-2697 v3的单线程上进行的。表中列出的计算时间和内存使用情况15,810个基因的平均值。
表3算法1：基于个体水平数据的moPMR-Egger,应用R包PMR中的moPMR_individual函数。
Figure BDA0002709151730000211
Figure BDA0002709151730000221
表4算法2：基于汇总统计量数据的moPMR-Egger,应用R包PMR中的moPMR_summary函数。
Figure BDA0002709151730000222
Figure BDA0002709151730000231
应注意的是，以上实例仅用于说明本发明的技术方案而非对其进行限制。尽管参照所给出的实例对本发明进行了详细说明，但是本领域的普通技术人员可根据需要对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围。

Claims (11)

1.一种基于孟德尔随机化的多性状全转录组关联分析方法，其特征在于，所述分析方法包括依托两样本MR框架，在没有样本重叠的两项独立研究中测量暴露和结局，从而估计和检验暴露对结局的因果效应；具体的，一次检验一个基因，并将cis-SNPs作为工具变量，在似然框架下同时估计和检验基因对多个结局性状的因果效应，同时能将潜在的水平多效性嵌入模型。
2.如权利要求1所述的多性状全转录组关联分析方法，其特征在于，所述分析方法包括对个体水平数据进行分析和汇总统计量数据进行分析。
3.如权利要求2所述的多性状全转录组关联分析方法，其特征在于，当对个体水平数据进行分析时，所述分析方法包括采用如下回归模型：
Figure RE-FDA0002750042100000011
Figure RE-FDA0002750042100000012
Figure RE-FDA0002750042100000013
式(2)-(3)合并为：
Figure RE-FDA0002750042100000014
其中，
Figure RE-FDA0002750042100000015
优选的，遵循遗传学中多基因建模假设，假设β的所有元素非零，且每个元素都服从正态分布
Figure RE-FDA0002750042100000016
同时，遵循PMR-Egger假设，假设对于每个性状i，每个SNP的水平多效性相等。
4.如如权利要求2所述的多性状全转录组关联分析方法，其特征在于，当所述分析方法对汇总统计量数据进行分析时，模型可表示如下：
Figure RE-FDA0002750042100000017
Figure RE-FDA0002750042100000018
5.一种基于孟德尔随机化的多性状全转录组关联分析系统，其特征在于，所述系统包括：
数据获取单元，用于获取相关数据信息并进行标准化处理；
模型构建单元，用于基于获取的数据进行多结局建模，在似然框架下同时估计和检验基因对多个结局性状的因果效应，同时能将潜在的水平多效性嵌入模型。
6.如权利要求5所述的多性状全转录组关联分析系统，其特征在于，当所述系统对个体水平数据进行分析时，
所述数据获取单元获取的相关数据信息包括基因表达的个体数据、表达数据对应的基因信息，全基因组关联分析的多维性状数据和全基因组关联分析对应的基因信息；
优选的，所述相关数据信息均标准化处理。
7.如权利要求5所述的多性状全转录组关联分析系统，其特征在于，所述模型构建单元中采用如下回归模型：
Figure RE-FDA0002750042100000021
Figure RE-FDA0002750042100000022
Figure RE-FDA0002750042100000023
式(2)-(3)合并为：
Figure RE-FDA0002750042100000024
其中，
Figure RE-FDA0002750042100000025
优选的，遵循遗传学中多基因建模假设，假设β的所有元素非零，且每个元素都服从正态分布
Figure RE-FDA0002750042100000026
同时，遵循PMR-Egger假设，假设对于每个性状i，每个SNP的水平多效性相等。
8.如权利要求5所述的多性状全转录组关联分析系统，其特征在于，当所述系统对汇总统计量数据进行分析时，
所述数据获取单元获取的相关数据信息包括基因表达的汇总数据，全基因组关联分析的多维性状汇总数据，基因表达遗传变异的连锁不平衡矩阵，全基因组关联分析遗传变异的连锁不平衡矩阵和纳入分析的性状间的相关矩阵；
所述模型构建单元中采用如下模型：
Figure RE-FDA0002750042100000027
Figure RE-FDA0002750042100000028
9.一种电子设备，包括存储器和处理器以及存储在存储器上并在处理器上运行的计算机指令，所述计算机指令被处理器运行时，完成权利要求1-4任一项所述基于孟德尔随机化的多性状全转录组关联分析方法所进行的步骤。
10.一种计算机可读存储介质，用于存储计算机指令，所述计算机指令被处理器执行时，完成权利要求1-3任一项所述基于孟德尔随机化的多性状全转录组关联分析方法所进行的步骤。
11.权利要求1-4任一项所述方法、权利要求5-8任一项所述系统、权利要求9所述电子设备和/或权利要求10所述计算机可读存储介质在如下任意一种或多种中的应用：
1)复杂疾病内在分子遗传机制研究；
2)疾病易感性研究；
3)药物反应研究；
4)药物靶标研究；
5)药物筛选研究。
CN202011049640.3A 2020-09-29 2020-09-29 一种基于孟德尔随机化的多性状全转录组关联分析方法、系统及其应用 Active CN112185464B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011049640.3A CN112185464B (zh) 2020-09-29 2020-09-29 一种基于孟德尔随机化的多性状全转录组关联分析方法、系统及其应用

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011049640.3A CN112185464B (zh) 2020-09-29 2020-09-29 一种基于孟德尔随机化的多性状全转录组关联分析方法、系统及其应用

Publications (2)

Publication Number Publication Date
CN112185464A true CN112185464A (zh) 2021-01-05
CN112185464B CN112185464B (zh) 2022-12-13

Family

ID=73946599

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011049640.3A Active CN112185464B (zh) 2020-09-29 2020-09-29 一种基于孟德尔随机化的多性状全转录组关联分析方法、系统及其应用

Country Status (1)

Country Link
CN (1) CN112185464B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113113141A (zh) * 2021-04-02 2021-07-13 北京果壳生物科技有限公司 一种基于孟德尔随机化评估微量营养素与精神类疾病因果关系的方法
CN113808665A (zh) * 2021-09-29 2021-12-17 山东大学 全基因组致病snp精细定位的因果关联分析方法
CN114171110A (zh) * 2021-10-20 2022-03-11 山东大学 一种基于联合似然的孟德尔随机化分析方法
CN114898809A (zh) * 2022-04-11 2022-08-12 中国科学院数学与系统科学研究院 适用复杂性状的基因-环境交互的分析方法及存储介质

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105205344A (zh) * 2015-05-18 2015-12-30 上海交通大学 基于多目标蚁群优化算法的基因位点挖掘方法
CN105740649A (zh) * 2016-01-22 2016-07-06 浙江大学 一种基于混合线性模型的多性状关联分析方法
CN111341448A (zh) * 2020-03-03 2020-06-26 西安交通大学 一种基于孟德尔随机化预测复杂疾病及表型相关代谢物的方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105205344A (zh) * 2015-05-18 2015-12-30 上海交通大学 基于多目标蚁群优化算法的基因位点挖掘方法
CN105740649A (zh) * 2016-01-22 2016-07-06 浙江大学 一种基于混合线性模型的多性状关联分析方法
CN111341448A (zh) * 2020-03-03 2020-06-26 西安交通大学 一种基于孟德尔随机化预测复杂疾病及表型相关代谢物的方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
ZHONGSHANG YUAN ET AL.: "Testing and controlling for horizontal pleiotropy with probabilistic Mendelian randomization in transcriptome-wide association studies", 《NATURE COMMUNICATIONS》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113113141A (zh) * 2021-04-02 2021-07-13 北京果壳生物科技有限公司 一种基于孟德尔随机化评估微量营养素与精神类疾病因果关系的方法
CN113808665A (zh) * 2021-09-29 2021-12-17 山东大学 全基因组致病snp精细定位的因果关联分析方法
CN113808665B (zh) * 2021-09-29 2024-03-08 山东大学 全基因组致病snp精细定位的因果关联分析方法
CN114171110A (zh) * 2021-10-20 2022-03-11 山东大学 一种基于联合似然的孟德尔随机化分析方法
CN114898809A (zh) * 2022-04-11 2022-08-12 中国科学院数学与系统科学研究院 适用复杂性状的基因-环境交互的分析方法及存储介质
CN114898809B (zh) * 2022-04-11 2022-12-23 中国科学院数学与系统科学研究院 适用复杂性状的基因-环境交互的分析方法及存储介质

Also Published As

Publication number Publication date
CN112185464B (zh) 2022-12-13

Similar Documents

Publication Publication Date Title
CN112185464B (zh) 一种基于孟德尔随机化的多性状全转录组关联分析方法、系统及其应用
Zeng et al. Signatures of negative selection in the genetic architecture of human complex traits
CN107229841A (zh) 一种基因变异评估方法及系统
IL258999A (en) Methods for detecting copy-number variations in next-generation sequencing
EP3871232A1 (en) Methods and apparatus for phenotype-driven clinical genomics using a likelihood ratio paradigm
Wu et al. Joint analysis of GWAS and multi-omics QTL summary statistics reveals a large fraction of GWAS signals shared with molecular phenotypes
EP4035163A1 (en) Single cell rna-seq data processing
Hernandez et al. Singleton variants dominate the genetic architecture of human gene expression
CN108268752B (zh) 一种染色体异常检测装置
CN116564410A (zh) 一种预测突变位点顺式调控基因的方法、设备和介质
CN109461473B (zh) 胎儿游离dna浓度获取方法和装置
US20240221954A1 (en) Disease prediction methods and devices, electronic devices, and computer readable storage media
Binder et al. Cluster-localized sparse logistic regression for SNP data
Zhang et al. MaLAdapt reveals novel targets of adaptive introgression from Neanderthals and Denisovans in worldwide human populations
WO2019242445A1 (zh) 病原体操作组的检测方法、装置、计算机设备和存储介质
CN114388062A (zh) 基于机器学习预测抗生素抗性表型的方法、设备及应用
CN116486913B (zh) 基于单细胞测序从头预测调控突变的系统、设备和介质
CA3190381A1 (en) Detecting cross-contamination in sequencing data
Liu et al. Analyzing association mapping in pedigree‐based GWAS using a penalized multitrait mixed model
CN107109324B (zh) 确定胎儿核酸含量的方法和装置
Wijsman et al. Joint oligogenic segregation and linkage analysis using Bayesian Markov chain Monte Carlo methods
Steuerman et al. Exploiting gene-expression deconvolution to probe the genetics of the immune system
KR20240046481A (ko) 지문 분석을 이용하여 화합물을 생리학적 조건과 연관시키는 시스템 및 방법
Zhang An improved nonparametric approach for detecting differentially expressed genes with replicated microarray data
Temple et al. Modeling recent positive selection in Americans of European ancestry

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