NC教你：孟德尔随机化下一个新风口-单细胞MR

关注医豌豆 助力SCI

1. csMR概述

csMR是基于Snakemake工作流程管理系统开发的。图1展示了分析框架。csMR现已在https://github.com/rhhao/csMR上公开提供。

2. 数据集

2.1 eQTL数据

在这项研究中，我们通过整合单细胞eQTL数据将遗传变异分为不同的脑细胞类型。汇总水平的大脑单细胞eQTL数据来自Bryois等人进行的一项研究[7]。这是目前唯一一项通过进行单细胞RNA测序对成人人脑中所有主要细胞类型进行的eQTL研究。简而言之，单细胞eQTL数据集包含主要来自192个欧洲个体的大脑前额叶和颞叶皮质的8种细胞类型（星形胶质细胞、内皮细胞、兴奋性神经元、抑制性神经元、小胶质细胞、少突胶质细胞、少突胶质前体细胞和周细胞）的eQTL结果。根据原始研究中的描述，定量了约14,595个基因的表达，总共对530万个SNP进行了基因分型。

为了了解使用单细胞eQTL划分遗传变异在发现更多潜在因果关系方面是否优于使用批量/bulk-eQTL，我们另外获得了来自GTEx v8报告的bulk脑组织的eQTL结果。在这项研究中，我们使用了来自大脑皮层和额叶皮层的bulk-eQTL，因为它们是从与单细胞eQTL相似的大脑区域产生的。样本量也与单细胞eQTL研究相当。在每个组织中，测试了约1040万个SNP与超过20,000个基因的表达的相关性。本研究中使用的eQTL数据集的详细信息可在补充数据1中找到。

2.2 肥胖相关性状的GWAS数据

在这里，我们测试了csMR在肥胖相关特征的GWAS数据集上的性能。首先，我们利用了对GIANT联盟和UKBB的GWAS数据进行的meta分析得到的BMI的GWAS汇总统计数据 ，总样本量为681,275名欧洲血统的个体。我们还纳入了另外两个肥胖指标WHRadjBMI和体脂百分比(https://gwas.mrcieu.ac.uk/datasets/ukb-b-8909/)的汇总水平GWAS数据，分别包含694,649和454,633名欧洲个体。根据这些原始研究中样本纳入的描述，超过60%的受试者来自UKBB。因此，我们在后续分析中使用了来自UKBB的50,000名欧洲个体的基因组数据作为参考面板。我们通过对不相关的欧洲祖先个体进行采样来构建基于UKBB(https://biobank.ndph.ox.ac.uk/)获得的遗传亲属信息（数据字段22021）、种族背景信息（数据字段21000）和遗传种族分组信息（数据字段22006）的参考数据。补充数据1总结了肥胖GWAS数据资源的更多信息。

2.3 结局的GWAS数据

为了研究肥胖对复杂的疾病是否存在细胞分层因果效应，我们进行了孟德尔随机化分析，将18种疾病特征作为先前从观察性或孟德尔随机化研究中证明与肥胖密切相关的结果。根据医学主题词（MeSH）（https://meshb-prev.nlm.nih.gov/）定义的疾病类别术语，将这些结果进一步分为4类，包括神经系统疾病、精神障碍性疾病、肌肉骨骼疾病、心血管和代谢疾病。每个类别包含4-5种疾病特征，我们在补充数据4中总结了支持它们与肥胖关系的相关出版物。疾病特征的GWAS汇总统计数据从补充数据4中列出的不同来源获得。所有研究都声称没有来自UKBB和GIANT联盟的重叠个体，所有参与样本或>80%的样本都是欧洲人。这些研究的样本量从9954到446,696不等。

3. 共定位分析

共定位方法已经成功地解决了在分子性状(例如基因表达)和疾病性状之间发现共有因果变异的问题。在这项研究中，我们通过将共定位策略应用于GWAS和单细胞eQTL数据来对遗传效应进行分层，以确定在特定的细胞类型中，与GWAS性状和基因表达相关的因果变异。使用“coloc”包进行了共定位分析，该软件包采用了SuSiE方法在共定位前进行fine mapping。coloc的作者推荐SuSiE方法，因为它优于其他方法，尤其是在存在多个因果变异的情况下。基于包含如上所述的UKBB中50,000个不相关个体的参考面板生成LD矩阵。排除MHC区内的变异(chr6：2847779733448354，GRCh37组装)。通过设置“coverage=0.9”来识别可信集。我们使用默认参数来评估GWAS信号和100 kb区域内基因表达之间的共定位。最终估计后验概率以支持5个竞争假设：（1）与任一性状均无关联（PPH0）；（2）仅与性状1相关（PPH1）；（3）仅与性状2相关（PPH2）；（4）与两个性状相关，但具有不同的因果变异（PPH3）；（5）与两个性状的相关，具有共同的因果变异（PPH4）。然后使用PPH4评估GWAS和基因表达之间共定位的证据。使用普遍接受的阈值(PPH4>0.8）来寻找共定位的有力证据。“coloc”的输出还包括GWAS和eQTL关联的推定因果变异。具有与单细胞eQTL共定位的强有力证据的GWAS hit被进一步视为后续MR分析的候选工具变量。

4. 遗传工具的选择

接下来，我们从共定位分析采用的因果变异中选择合格的SNP用于MR分析。合格的遗传工具变量应满足三个核心假设：（1）它应与暴露密切相关（相关性假设）；（2）不应与暴露-结果关联的任何潜在混杂因素相关（独立性假设）；（3）仅通过暴露与结局相关（排除限制假设）。根据这些假设，我们首先选择了与暴露密切相关的SNP（P<5e-8）。通过使用“phenoscanner”v1.0 R包检查变异或其代理变异（r2>0.8）是否与混杂因素相关，进一步过滤掉变异。考虑了三个潜在的混杂因素，包括教育、饮酒和吸烟行为。接下来，我们使用PLINK提供的LD clumping程序，通过设置0.001的r2阈值、10,000 kb的窗口大小和5e−8的P值阈值，从剩余IVs中鉴定index SNP。然后将这些SNP与结局数据取交集。对于结局数据集中未显示的数据，我们用包含的代理（r2>0.8）替换了它们。接下来，使用“TwoSampleMR”v0.5.6 R包进行数据haormonise，以确保对暴露和结局的影响对应于相同的等位基因。同时，具有中等等位基因频率（>0.42）的回文SNP被去除。

5. 工具变量的质量控制

对于满足三个MR假设的IVs，我们执行质量控制流程以确保IVs质量。我们首先使用“RadialMR”1.0 R包65实施异质性测试以过滤掉异常多效性SNP。分别使用参数“ivw_radial(α=0.05，权重=1，tol=0.0001)”和“egger_radial(α=0.05，权重=1)”进行Cochran's Q检验和Rrucker's Q检验以评估水平和定向多效性效应。丢弃P阈值为0.05的异常值。此外，我们计算了F-statistic，以测量IVs的功效，公式为F=R2（N-K-1）/k（1-R2），其中N代表暴露的样本量；k表示IV的数量；R2代表由遗传变异解释的表型变异的比例，并且是每个snp解释的表型变异的总和。使用公式R2=β 2/(β 2+SE2 × N)根据效应大小(β)和标准误差(SE)计算R2。为避免工具偏倚，F统计量应至少为10。

6. 双样本MR

通过进行双样本MR分析来估计在每种脑细胞类型中鉴定的遗传工具（本研究中的肥胖相关变异）与疾病的相关性。基于不同假设的补充方法用于寻找稳健的关联，包括IVW、MR-Egger、weighted median、MR-robustadjustedprofile评分（MR-RAPS）和Weighted Mode。IVW方法假设平衡多效性，并实施乘法随机效应模型，以提供对因果效应的一致估计。MR-Egger方法通过Egger回归的斜率系数估计因果效应，该系数基于独立于直接效应（内部）假设的工具强度。MR-RAPS方法利用稳健的调整概况分数来完成统计推断，这也受到不平衡多效性的偏倚。weighted median法在假设工具总权重的至少50%来自有效变量的情况下估计因果效应。假设最大的SNP组是有效的工具，weighted mode提供了一致的估计。当只有一个可用的遗传工具时，Wald Ratio Method用于估计因果效应。我们使用“TwoSampleMR”v0.5.6 R包进行这些MR分析。在这项研究中，在对180项MR测试（18种疾病结果×（8种脑细胞类型+2种脑组织）进行Bonferroni校正后，确定了可靠的因果推断；P<2.78e-4（0.05/180））。

7. 多效性和敏感性分析

• （a）如果仅有一个遗传工具，则使用Wald ratio来估计因果效应。

• （b）对于多种遗传工具，如果没有检测到定向多效性，即Cochran’s Q、MR-Egger截距和MR-PRESSO的检验的P>0.05，我们选择IVW作为主要方法，因为它在检测因果效应方面比其他方法更有效。

• （c）在检测到定向多效性的情况下，违反了IVW和MR-RAPS的假设，则如果Rucker’s Q检验的P>0.05，则认为MR-Egger是主要的MR方法。

• （d）如果检测到定向多效性，并且Rucker’s Q检验的P<0.05，则weighted median和weighted mode可用于估计因果效应。weighted median更优选，因为其在检测因果效应方面具有更高的功效。

结果

1. 细胞分层MR（cell-stratified MR）分析的框架

csMR的概述如图1所示。基本上，csMR采用三个步骤来推断每种细胞类型的因果关系：（1）在通过fine-mapping识别因果变异后，找到GWAS和单细胞eQTL数据之间的遗传共定位；(2)过滤每种细胞类型中的共定位信号以选择合适的工具变异(IV)；(3)用选定的IVs进行双样本MR分析。脑单细胞eQTL数据来自Bryois等人，其中包括源自大脑皮层的8种细胞类型（兴奋性神经元（ExN）、抑制性神经元（InN）、星形胶质细胞、小胶质细胞、少突胶质细胞（ODC）、少突胶质前体细胞（OPC）、内皮细胞（EC）和周细胞）的汇总统计数据。在多个因果变异的假设下，我们采用SuSiE方法进行fine-mapping，采用“coloc”框架进行共定位分析。假设4（PPH4）的后验概率表明GWAS和QTL之间是否存在共同的因果关联，用于评估共定位。在每种细胞类型中，如果PPH4大于0.8，我们确定了显著的共定位。然后将共定位的因果SNP传递到下一步，以选择合适的IV用于MR分析，其包括SNP过滤、clump和质量控制（图1）。最后用细胞分层方法进行双样本MR分析，以评估特定细胞类型背景下暴露特征对结果的因果影响。我们在这里实施了六种MR方法来补偿每种方法的内在偏差，并根据多效性分析的结果对主要MR方法进行了优先排序。

2. BMI的共定位特征提示风险变异的细胞类型特异性调节

我们将BMI作为感兴趣的特征，以显示csMR如何从共定位分析中将遗传变异划分到不同的脑细胞。我们获得了Yengo等人报告的681,275名欧洲血统参与者的BMI的GWAS汇总数据。这是来自GIANT联盟和UKBB的meta分析。在最初的研究中，发现941个接近独立的SNP与BMI相关（在修订的全基因组显著性阈值P<1e-8）。在共定位分析后，如图2a所示，我们在脑细胞中鉴定了约235个共定位SNP，平均数量为171个。发现这些SNP影响约168个基因的表达（平均120个基因）。对于共定位的SNP，细胞特异性结果的比例范围为61.3%至82.2%，对于共定位的基因，细胞特异性结果的比例范围为74.2%至91.3%（图2a）。接下来，我们研究了其他细胞类型中共定位SNP的平均eQTL效应，并发现这些SNP在其他测试细胞中具有不太明显的效应大小（图2b），这表明BMI相关SNP倾向于以细胞类型特异性方式调节基因表达。此外，对于共定位基因，我们还观察到90.2%的这些基因在特定细胞类型中被识别（图2c）。接下来，我们有兴趣知道共定位基因如何在不同的脑细胞类型中表达。我们从Morabito等人获得了健康人脑的scRNA-seq数据。已经产生了星形胶质细胞、ExN、InN、小胶质细胞、ODC和OPC的基因表达矩阵。跨细胞的成对比较表明，所有共定位基因和在特定细胞类型中鉴定的基因都具有细胞类型特异性表达趋势（图2d）。总之，这些结果表明，BMI相关遗传变异可能通过影响特定细胞类型中特定基因的表达而充当脑eQTL。除了这里报道的BMI，在几种脑部疾病中也观察到了细胞特异性共定位。

在这里，我们展示了一个实例，其中BMI相关基因座与源自特定脑细胞类型的eQTL共定位。如图2e中所示的locuszoom图，在表达POMC的4种细胞类型中，我们仅在InN（PPH4=1）中检测到强烈的共定位证据。因果变异rs564667与BMI变异（P=6.99e-40）和InN中POMC表达（P=0.018）显著相关。POMC基因编码阿黑皮质素原，长期以来一直被研究其在能量稳态和体重控制中的作用。我们仅在抑制性神经元中观察到POMC的共定位信号，这与其已知的神经调节功能一致。

3. 孟德尔随机化推断BMI的细胞分层因果效应

鉴于共定位分析表明的风险变异的细胞类型特异性效应，csMR后续使用在每个细胞中识别为IVs的共定位SNP，并评估感兴趣性状/疾病的细胞分层因果关系。对于BMI，我们探索了属于4个不同疾病类别的18种疾病的因果效应，这些疾病先前已显示出与BMI的相关性。我们首先使用P<5e-8的所有相关变异进行MR分析，以推断BMI和18种结果之间的因果关系。结果显示，除了阿尔茨海默病和帕金森病，BMI与几乎所有疾病结局都有因果关系。先前的MR研究也得出结论，BMI和阿尔茨海默病或帕金森病之间没有因果关系，这与我们的结果一致。然后，我们专注于调查推断的因果关系是否按脑细胞分层。进行IV选择以过滤共定位的SNP。进一步对与BMI（P<5e-8）但与混杂因素（吸烟、饮酒、教育）无关的SNP进行连锁不平衡（LD）修剪和异常值去除。通过所有质量控制测试的IV（NIV）数量范围为4至34。经过MR和额外的多效性和敏感性分析，我们发现BMI对7种复杂疾病有显著的因果影响，P<2.78e-4（0.05/180）（图3），包括睡眠障碍（神经系统疾病）、注意力缺陷多动障碍（ADHD）（精神障碍）、痛风和骨质疏松症（肌肉骨骼疾病）、冠状动脉疾病（CAD）、心肌梗死和II型糖尿病（T2D）（心血管和代谢疾病）。这些因果关系主要在6种脑细胞类型中观察到，即星形胶质细胞、EC、ExN、小胶质细胞、ODC和OPC。

值得注意的是，我们在多种细胞类型中发现了BMI与睡眠障碍、多动症、痛风、CAD和T2D的因果关系。特别是，由EC、ExN、小胶质细胞、ODC和OPC中的细胞分层变异预测的较高BMI与睡眠障碍风险增加相关（BMI每增加1-SD的IVW β【95%CI】：EC，0.66【0.37-0.94】，P=7.61 × 10−6，NIV=19；ExN，0.84【0.62-1.06】，P=1.64 × 10−13，NIV=26；小胶质细胞，1.03【0.61-1.46】，P=2.23 × 10−6，NIV=14；ODC，0.58【0.32-0.83】，P=8.47 × 10−6，NIV=14；OPC，0.70【0.44-0.96】，P=8.16 × 10−7，NIV=21）。BMI和ADHD之间的正相关主要在星形胶质细胞和EC中观察到（IVW β【95%CI】：星形胶质细胞，0.66【0.36-0.97】，P=1.85 × 105，NIV=22；EC，0.66【0.34-0.98】，P=5.70 × 105，NIV=21）。此外，我们在ExN、小胶质细胞和ODC中发现了BMI和痛风之间的显著因果关系（ExN，1.09【0.65-1.52】，P=8.32 × 107，NIV=34；小胶质细胞，1.71【0.86-2.56】，P=8.18 × 105，NIV=13；ODC，1.33【0.77-1.89】，P=3.19 × 106，NIV=21）。此外，在星形胶质细胞和ExN中观察到与CAD的细胞分层正因果关系（IVW β【95%CI】：星形胶质细胞，0.21【0.12-0.30】，P=8.36 × 106，NIV=26；ExN，0.24【0.15-0.33】，P=1.56 × 107，NIV=26）。我们还发现，星形胶质细胞、EC、ExN和ODC中较高的遗传预测BMI与较高的T2D风险有因果关系（IVW β[95%CI]：星形胶质细胞，0.24[0.12-0.35]，P=4.76 × 105，NIV=20；EC，0.41[0.25-0.56]，P=2.14 × 107，NIV=16；ExN，0.26[0.16-0.35]，P=1.59 × 107，NIV=27；ODC，0.25[0.14-0.36]，P=7.74 × 106，NIV=21）。这些结果强调，BMI对一些复杂疾病的影响是由来自多个脑细胞的调节后果潜在驱动的。

此外，我们检测到骨质疏松症和心肌梗死的细胞特异性因果关系。心肌梗死的较高风险受EC分层IVs估计的较高BMI的因果影响（IVW β【95%CI】：0.27【0.15-0.40】，P=2.55 × 105，NIV=16）。虽然根据ExN分层IV估计，BMI增加1 SD与骨质疏松症风险降低相关（IVW β[95%CI]：-0.33[-0.49--0.17]，P=5.64 × 10-5，NIV=30）。这些观察结果强调了BMI对骨质疏松症和心肌梗死的细胞特异性影响。

4. 通过单细胞eQTL对因果关系进行分层比使用基于组织的数据发现更多的关联

我们有兴趣知道使用单细胞eQTL确定的分层因果关系是否与使用基于组织的eQTL数据确定的分层因果关系不同。由于csMR也与其他汇总水平的QTL结果兼容，我们接下来将csMR应用于来自大脑皮层和额叶皮层的GTEx eQTL数据。如图3所示，我们使用大脑额叶皮层的分层风险变异成功复制了BMI对睡眠障碍和T2D的因果效应（IVW β【95%CI】：睡眠障碍，0.59【0.34-0.83】，P=2.79 × 106，NIV=26；T2D，0.30【0.20-0.39】，P=3.19 × 109，NIV=33）。然而，使用单细胞eQTL数据指出的与其他5种疾病的关联在Bonferroni校正后未能达到显著水平（P<2.78e-4），但仅显示出ADHD、痛风、CAD和心肌梗死的名义显著性（P<0.05）。此外，我们没有使用基于组织的eQTL检测到新的关联。这些结果表明了整合单细胞数据以发现隐藏的指导性信息的优势。

5. 将csMR应用于其他肥胖相关特征突出了对不同结局的因果影响

接下来，我们试图调查csMR在其他性状上的表现。尽管BMI通常用于近似总体过量的身体脂肪，但也有其他替代物用于捕获与肥胖相关的异常代谢后果。例如，中心性肥胖的评估可以通过腰臀比(WHR)，而体脂百分比的测量可以提供脂肪组织的准确脂肪质量百分比。作为公认的肥胖指标，BMI和其他肥胖相关特征都被发现与大脑有关，而最强的遗传因素似乎在很大程度上是不同的。在这里，我们将csMR应用于另外两个与肥胖相关的特征，调整了BMI的WHR（WHRadjBMI）和体脂百分比，以观察它们是否表现出与BMI不同的细胞分层因果关系。

在18个结果中，我们最终识别到了WHRadjBMI对中风、CAD和T2D的因果影响（图4a）。未从BMI中发现对卒中的影响。确定了3种有效的细胞类型来支持WHRadjBMI和中风之间的因果关系（IVW β【95%CI】：星形胶质细胞，0.76【0.42-1.09】，P=9.71 × 106，NIV=16；ExN，0.66【0.33-0.98】，P=7.96 × 105，NIV=21；InN，0.61【0.30-0.92】，P=1.14 × 104，NIV=18）。值得注意的是，尽管发现CAD和T2D与BMI和WHRadjBMI都相关，但有效细胞类型是不同的。在InN中观察到WHRadjBMI对CAD的因果效应（IVW β【95%CI】：0.30【0.19-0.42】，P=3.84 × 107，NIV=12），而在星形胶质细胞和ExN中观察到BMI对CAD的影响。在InN和OPC中观察到WHRadjBMI对T2D的影响（IVW β【95%CI】：InN，0.24【0.12-0.35】，P=7.95 × 105，NIV=16；OPC，0.26【0.13-0.39】，P=7.61 × 105，NIV=12），然而，在星形胶质细胞、EC、ExN和ODC中发现了BMI和T2D之间的关联。

此外，我们认识到体脂百分比与癫痫、ADHD、精神分裂症以及膝关节和髋关节骨关节炎之间的因果关系（图4b）。除ADHD外，在关注BMI和WHRadjBMI的分析中未观察到这些受影响的结果。特别是，我们发现体脂百分比与OPC中的癫痫正相关（IVW β【95%CI】：0.97【0.46-1.48】，P=1.75 × 104，NIV=7），与ODC和OPC中的膝关节和髋关节骨关节炎正相关（IVW β【95%CI】：ODC，0.89【0.42-1.35】，P=1.80 × 104，NIV=17；OPC，1.39【0.68-2.10】，P=1.18 × 104，NIV=8），而与小胶质细胞中的精神分裂症负相关（IVW β【95%CI】：1.08【1.58-0.57】，P=3.14 × 105，NIV=9）。总之，这些结果表明，不同的肥胖测量可能通过脑细胞特异性调节途径在遗传上导致不同疾病结局的发生率。

结论

本研究开发了一个分析框架--csMR（cell-Stratified MR），通过整合GWAS和单细胞eQTL的汇总水平数据来发现细胞分层因果关系。在将其应用于BMI·GWAS数据后，我们展示了csMR的应用价值，它有助于对靶细胞进行优先排序，以推断BMI和18种疾病结果之间的因果关系，这通常在批量分析中是隐藏的。我们的方法具有时间和成本效益，将有助于更好地理解复杂疾病之间的细胞特异性遗传联系。

文章解读到此结束，看完之后有没有觉得研究课题这就来了！