MR文献导读（一）：通过孟德尔随机化分析确定多发性硬化症的潜在药物靶点

今天我们来学习一篇新鲜出炉的MR文章，看看其思路和方法能不能为己所用吧~

为什么选择多发性硬化症？

多发性硬化症是一种复杂的自身免疫性疾病，目前已开发出多种治疗多发性硬化症的疗法，并得到广泛应用。

然而，现有的多发性硬化症药物由于不能抑制复发和缓解疾病进展，效果远不能令人满意。

因此，目前仍需要新的药物靶点来预防多发性硬化症。

用了什么方法？

我们利用国际多发性硬化症遗传学联合会（nCase = 47 429，nControl = 68 374）的汇总数据，并在英国生物银行（nCase = 1356，nControl = 395 209）和芬兰基因队列（nCase = 1326，nControl = 359 815）中进一步验证，进行了孟德尔随机化，以探索多发性硬化症的潜在药物靶点。

734 种血浆蛋白和 154 种CSF蛋白的遗传变异来自最近发表的全基因组关联研究。

为了进一步巩固孟德尔随机分析的结果，利用双向孟德尔随机化分析和 Steiger 过滤、贝叶斯共定位和表型扫描 [phenotype scanning] 进行了反向因果关系推测，以寻找之前报道的遗传变异与性状的关联。

此外，还进行了蛋白互作网络分析，以揭示蛋白质和/或现有多发性硬化症药物之间的潜在关联。

目的是什么？

将血浆蛋白和CSF（脑脊液总蛋白）蛋白作为多发性硬化症的潜在治疗靶点。

具体如何实现？

首先，我们利用国际多发性硬化症遗传学联合会（IMSGC）的 GWAS 数据、Zheng 研究中的血浆 pQTL 数据和 Yang 研究中的 CSF pQTL 数据，使用 MR 来识别多发性硬化症的潜在的致病血浆和 CSF 蛋白质成分。

其次，利用反向因果关系检测、贝叶斯共定位分析和表型检测进一步验证了主要发现。

第三，我们绘制了已识别蛋白质之间、基于血浆和 CSF 的蛋白质之间以及已识别蛋白质与当前多发性硬化症药物靶点之间的相互作用网络。

最后，我们利用英国生物银行和芬兰基因队列中的 GWAS 数据以及两项新发表的研究中的血浆 pQTL 数据进行了重复分析，作为外部验证，以加强我们的结论。

这里简要介绍一下GWAS & QTL的区别吧~

• GWAS

全基因组关联研究（Genome-Wide Association Study，简称 GWAS），是在全基因组范围内，检测多个个体的遗传变异多样性，获得群体中每个个体的基因型；然后与性状（即我们常说的表型）进行统计学关联分析，根据统计量（主要指 P 值）筛选出候选变异位点和基因。

• QTL
  • eQTL
  • pQTL
  • sQTL
  • meQTL

  数量性状位点（Quantitative Trait Locus mapping，简称 QTL 定位），指的是控制数量性状的基因在基因组中的位置。QTL可分为cis-QTL和trans-QTL。其原理就是测定一群个体的某个数量性状（表型），以及它们的基因型（就是基因组上的一些遗传标记，例如SNP），然后寻找基因型和表型的对应关系。 eQTL eQTL（expression QTL，基因表达量），即控制数量性状基因的mRNA表达水平高低的那些位点。 【比如，cis-eQTL指的是某个基因的eQTL定位到该基因所在的基因组区域，表明可能是该基因本身的差别引起的mRNA水平变化； trans-eQTL指的是某个基因的eQTL定位到其他基因组区域，表明其他基因的差别控制了该基因mRNA水平的差异。】 多数情况下，我们关心的是“这一个SNP位点在这一个位点及附近mRNA表达量之间的关系”，所以一般做的是cis-eQTL分析。 pQTL pQTL（protein QTL，蛋白质表达量），即控制数量性状基因的蛋白质表达水平高低的那些位点。 sQTL sQTL（splicing QTL，选择性剪接），即影响剪接调控的剪接数量性状位点。 meQTL meQTL（methylation QTL，DNA甲基化），表观调控相关的DNA甲基化数量性状位点。

接下来就是具体的材料和方法——

流程图

CSF 和血浆蛋白定量性状位点的获取

CSF pQTL原始数据来源于：Yang et al.（ 274 pQTLs of 184 CSF proteins）。

哪些pQTL 数据能够被纳入：

(i) 全基因组强相关（P < 5 × 10-8）； (ii) 位于主要组织相容性复合体（MHC）区域（chr6，26-34 Mb）之外； (iii) 连锁不平衡（LD） r2 < 0.001； (iv) 是cis-pQTL。

最后，为 154 种蛋白质确定了 154 个顺式-pQTL。

血浆pQTL数据来源于Zheng et al，整合了五个GWAS的公开数据。

根据上述 CSF pQTL 数据集的筛选标准，共纳入了 734 个蛋白质的 738 个顺式作用 SNPs。

此外，从【Pietzner 等人（在 10 708 名参与者中测量了 4775 种血浆蛋白）和 Ferkingstad 等人（在 35 559 名参与者中测量了 4907 种血浆蛋白）】这两项最近发表的研究中获取了血浆 pQTL 数据，用于后续外部验证。

ps：对于 QTL GWAS 数据中的任何缺失信息，如效应等位基因频率（effect allele frequency），我们使用与之匹配的人类基因组构建作为参考，以完成数据处理。

多发性硬化症的GWAS数据获取

数据取自 IMSGC 最大的 GWAS 数据集，其中包括 115 803 名欧洲血统的个体（nCase = 47 429，nControl = 68 374）。

验证数据集则从英国生物银行（截至 2017 年，nCase = 1356，nControl = 395 209）13 和 FinnGen 研究（nCase = 1326，nControl = 359 815，R6 版本）中获得。

注意看：在公共数据集中获得的数据，需要标明数据获取的具体时间和数据库版本号，因为这两个数据库是在保持更新的。

统计分析

MR分析

我们以血浆和CSF蛋白作为暴露对象，以 多发性硬化症 为结果，使用 "TwoSampleMR"（https://github.com/MRCIEU/TwoSampleMR）进行 MR 分析。

如果某一蛋白只有一个 pQTL，则使用 Wald 比率。当这个蛋白有两种或两种以上遗传变异时，采用反方差加权磁共振（MR-IVW），然后进行异质性分析。

多发性硬化症患病风险增加的风险比 (OR) 以血浆蛋白质水平每增加一个标准差 (SD) 和CSF蛋白质水平每增加 10 倍表示。

使用 Bonferroni 校正来调整多重检验，并使用 0.05/888 的临界 P 值（P < 5.63 × 10-5）来确定进一步分析结果的优先次序。

提问一下：888是怎么来的？为什么要用0.05除以888呢？

只对初步鉴定出的蛋白质进行 MR 外部验证，并将 P 值阈值设定为 0.05。

我们采用了同变异策略和显着变异策略：前者使用与初步分析相同的 SNP 作为遗传变异工具，后者使用全基因组显著 SNP 作为遗传变异来验证初步发现。

反向因果推断

按照与 pQTL 相同的筛选标准，从 IMSGC 的 GWAS数据 中筛选出 134 个多发性硬化症的遗传工具进行双向 MR 分析，以检测潜在的反向因果关系。

蛋白质的完整汇总统计数据来自之前的三项研究:

Yang C, Farias FHG, Ibanez L, et al. Genomic atlas of the proteome from brain, CSF and plasma prioritizes proteins implicated in neurological disorders. Nat Neurosci. 2021;24:1302-1312. Ferkingstad E, Sulem P, Atlason BA, et al. Large-scale integration of the plasma proteome with genetics and disease. Nat Genet. 2021;53:1712-1721. Sun BB, Maranville JC, Peters JE, et al. Genomic atlas of the human plasma proteome. Nature. 2018;558:73-79.

使用 MR-IVW、MR-Egger、加权中位数（weighted median）、简单模式和加权模式（simple mode and weighted mode）估算效应。我们还进行了 Steiger 过滤法，以确保蛋白质与多发性硬化症之间关联的方向性。

当 P < 0.05 时，结果具有统计学意义。

贝叶斯共定位分析

贝叶斯共定位分析使用 "coloc "软件包（https://github.com/chr1swallace/coloc）和默认参数来评估两个性状共享相同因果变异体的概率。

贝叶斯共定位提供了关于两个性状是否共享一个变异体的五个假设的后验概率。

我们测试了假设 3（PPH3）和假设 4（PPH4）的后验概率：

假设 3 中，蛋白质和 多发性硬化症 都是通过不同变体与该区域相关联的；

假设 4 中，蛋白质和 多发性硬化症 都是通过共享变体与该区域相关联的。

我们同时使用了 coloc.abf 和 coloc.susie 算法，并根据至少一种算法确定的基于基因的 PPH4 > 80% 的共定位证据来定义基因。

Phenotype scanning（表型扫描）

我们还进行了表型扫描，搜索以前的 GWAS数据以揭示已确定的 pQTL 与其他性状的关联。

表型扫描通过phenoscanner 和 Ferkingstad 等人的血浆蛋白质组 GWAS 研究进行。

满足以下条件的 SNP 被认为是存在基因多效性的：

(i)关联具有全基因组显著性（P < 5 × 10-8）； (ii)GWAS 在欧洲血统人群中进行； (iii)SNP 与任何已知的多发性硬化症风险因素相关，包括代谢特征、蛋白质或临床特征。

这一步的目的是为了剔除混杂因素的干扰，保证工具变量与混杂因素的相关性。

此外，我们还计算了蛋白质的 pQTLs 之间的 LD r2，以揭示它们之间潜在的联系。

phenoscanner web

[PhenoScanner (cam.ac.uk http://www.phenoscanner.medschl.cam.ac.uk/)

根据你的需要来设置“门槛”~

phenoscanner package

# install.packages("devtools")
library(devtools)
# install_github("phenoscanner/phenoscanner")
library(phenoscanner)

#  For example, to query multiple SNPs:
res <- phenoscanner(snpquery=c("rs10840293","rs10"))##查询与我们的工具变量存在显著关联的traits,如果这些traits可能影响结局，那么对应SNPs就是pleiotropic 的，要去掉。
head(res$results)
res$snps

PPI网络比较分析

我们假设，由于血脑屏障的存在，血浆和 CSF 鉴定出的 pQTL 之间几乎不存在相关性。

因此，我们通过斯皮尔曼相关性分析研究了利用MR分析得到的效应估计值在 CSF 和血浆中鉴定出的共享 pQTL 之间的相关性，并设定了不同的 P 值阈值，以探讨相关性是否会随着显著性水平的提高而发生变化。

我们对脑脊液或血浆分析中提示与多发性硬化症风险相关的蛋白质（ MR 分析 P <0.05）的蛋白质-蛋白质相互作用（PPI）网络进行了探索。

我们的目的是研究优先蛋白质之间的相互作用，以及利用血浆数据确定的蛋白质是否会与利用脑脊液数据确定的蛋白质发生相互作用。

此外，为了探讨这些多发性硬化症相关基因与已上市药物靶点之间的相互作用，我们从最近的一篇综述中获得了 13 种治疗多发性硬化症的药物，并根据 Drugbank 数据库（https://www.drugbank.ca）获得了相应的药物靶点。

McGinley MP, Goldschmidt CH, Rae-Grant AD. Diagnosis and treatment of multiple sclerosis: A review. JAMA. 2021;325: 765-779.

我们还搜索了靶向已确定的潜在致病蛋白的现有药物。

所有 PPI 分析均使用检索相互作用基因的搜索工具（STRING）数据库 11.5 版（https://string-db.org/）进行，要求的最低相互作用分值为 0.4。

此外，我们还使用 Wald 比率法进行了 MR 分析，并使用 coloc.abf 算法进行了贝叶斯共定位分析，优先考虑的蛋白既是暴露蛋白也是结果蛋白。我们将 MR 的 P 值小于 0.05 视为潜在的相互作用，将 PPH4 >0.8 视为潜在的共定位。

遗憾的是，这篇文献没有提供代码，但是作为初学者，先学习思路也是一个很不错的选择。 孟德尔随机化的合辑刚刚起步，小编也在持续更新自己的学习进度，计划以后孟德尔随机化的内容以文献分享和代码实操为主，穿插一些MR相关的分析方法的分享，如果大家有其他建议，也欢迎后台留言哦~

参考

• Potential drug targets for multiple sclerosis identified through Mendelian randomization analysis - PubMed (nih.gov)
• GWAS与QTL的分析内容与原理的比较 (qq.com) https://mp.weixin.qq.com/s/Zwysrj4UenHjpnGx9h26Mw)
• eQTL就那样～ https://www.jianshu.com/p/d23be33f5b6e)
• 不会做QTL分析，也能轻松找sQTL、eQTL、meQTL - https://zhuanlan.zhihu.com/p/404930954)
• 很多领域的研究用到的分子生物学的QTL（数量性状位点）定位其基本原理和算法是什么呢？- https://www.zhihu.com/question/27695566)
• 使用coloc 进行 QTL 共定位 - https://zhuanlan.zhihu.com/p/594971999)
• [GitHub - phenoscanner/phenoscanner: https://github.com/phenoscanner/phenoscanner)