孟德尔随机化研究溃疡性结肠炎与甲状腺功能减退发生风险的因果关系

花胤 王晓燕 王振 辛永宁 刘守胜

1南京医科大学青岛临床医学院（青岛市市立医院）（山东青岛 266011）；青岛市市立医院2国际门诊，3感染性疾病科，4临床研究中心（山东青岛 266011）

溃疡性结肠炎（ulcerative colitis，UC）是一种病程反复的炎症性肠病（inflammatory bowel diseases，IBD），常由复杂的环境、免疫以及遗传因素共同引起，其从结肠远端开始发病并向近端扩展直至全结肠［1］。作为一种免疫介导的炎症性肠病，虽然UC主要影响消化道，但患者常出现肠道以外症状，包括肌肉骨骼、皮肤、眼部、肝胆、血液系统及内分泌疾病等［2］。甲状腺功能减退（hypothyroidism，HT）作为一种内分泌疾病，常由一种自身免疫疾病桥本病引起，表现为免疫细胞攻击甲状腺，最终引起甲状腺的功能丧失导致HT 的发生［3-5］。有研究［6-8］发现UC 患者发生甲状腺疾病包括HT 的风险较正常个体更高。然而另外的研究［9］则发现IBD 患者与正常个体之间的HT 发病风险无明显差异。近期的一项回顾性研究［10］发现，IBD 会降低患者发生HT 的风险，并不推荐IBD 患者在缺少临床症状的情况下实施甲状腺超声检查。因此，目前虽已有较多研究报道了UC 与HT 之间的关系［11］，但结论仍然具有争议，并且既往研究大多是靠观察性研究得到的关联结论，可能受到某些混杂因素及反向因果关系的干扰。孟德尔随机化（mendelian randomization，MR）作为一种疾病之间因果关系的推断方法在遗传病学领域广泛应用，MR 通过遗传变异作为工具变量来推断暴露变量能否影响结局，可以避免观察性研究中可能出现的缺点。因此，本研究拟通过两样本孟德尔随机化（two sample mendelian randomization，TSMR）分析来确定UC 与HT 之间发病风险的相关性。

1 资料与方法

1.1 研究原理 孟德尔随机化的基本原理是通过工具变量来分析暴露与结局关联的因果关系。其分析的准确性需要满足下述3 个核心假设：（1）关联性假设：工具变量需要与暴露因素强相关；（2）排他性假设：工具变量只能通过暴露因素对结局产生影响；（3）独立性假设：工具变量需要和结局及混杂因素没有关联［12］（图1）。

图1 孟德尔随机化的3 个核心假设Fig.1 Key assumptions of Two-Sample MR Analysis

1.2 数据来源 通过搜索全基因组关联研究（genome wide association study，GWAS）。分别获得UC 和HT 的GWAS 汇总数据，其中UC 样本数据来自芬兰生物银行（https：//www.finngen.fi/fi），共计211 498个样本（观察组：1 198；对照组：210 300），共计16 380 453 个SNP 位点。HT 样本数据来自IEU open GWAS（https：//gwas.mrcieu.ac.uk/datasets/ebi-a-GCST90018862/），此数据库样本量为410 141（观察组：30 155；对照组：379 986），共计24 138 872 个SNP 位点（表1）。

1.3 工具变量的筛选 根据以下条件筛选暴露的工具变量：（1）选取与暴露的关联强度达到全基因组显著关联水平（P＜ 5 × 10-8）的SNPs；（2）排除SNPs 之间的连锁不平衡影响（linkage disequilibrium，LD），设置连锁不平衡参数（r2）阈值为0.01，遗传距离为1 000 kb，从而确保SNP的独立性；（3）排除位于回文序列中的工具变量；（4）使用PhenoScanner数据库剔除与混杂因素相关的遗传工具变量，重点排除表型与结局有相关意义的SNP。

1.4 统计分析 分别采用逆方差加权法（inverse variance weighted，IVW）、Egger 回归法（MR-Egger）、中位数加权法（weighted median）开展TSMR，并采用森林图（forest plot），散点图（scatter plot）可视化上述方法的OR值（odds ratio）及95%CI（confidence interval）。IVW 法通过Wald估计值获得加权平均值，是本研究的主要分析工具，其主要原理是将每个遗传变异看作独立的工具变量，并通过荟萃分析合并工具变量的因果效应估值，当所有的遗传变异都是有效的工具变量时，IVW 法是因果效应的最有效评价工具。MR-Egger 回归法与IVW 法类似，但可以通过回归截距评估工具变量的水平多效性。中位数加权法对因果效应的评估最少只需要50%的有效工具变量，因此能在存在无效工具变量的情况下，提供稳定的因果效应值。关于敏感性分析，本研究利用MR-PRESSO（MR Pleiotropy Residual Sum and Outlier，MR-PRESSO）法评价是否存在离群的SNP 位点［13］，采用MREgger 的截距值评估水平多效性，采用Cochran′s Q检测及漏斗图（funnel test）评估SNPs 中的异质性，采用留一法分析（leave one out）分析是否存在单个SNP 能够影响因果评估的结果，以上统计分析均通过开源软件R（版本4.3.1）及R包“TwoSampleMR”（版本0.5.7）进行。以P＜ 0.05 为差异有统计学意义。

1.5 工具变量评价 工具变量的强度通过统计量F计算，其公式为F=R2× (N-k- 1)/(1 -R2) ×k。其中R2代表工具变量解释的变异比例，N为暴露样本的样本量，k为工具变量的数量，而R2=(2 ×EAF×(1 -EAF) ×beta2)/(N×SE2)，EAF为次要等位基因频率，beta为等位基因效应值，SE为标准误差［14］。

2 结果

2.1 UC 对HT 的MR 分析

2.1.1 遗传工具变量确定 本研究选择与UC 相关具有全基因组显著性（P＜ 5 × 10-8）的遗传变量SNP 位点进行汇总，筛选出同时满足假设1、2、3 的SNP，得到4 个符合假设的SNP。利用Pheno Scanner 数据库去除与结局HT 具有相关意义的SNP（n= 0），通过GWAS 提取结局HT 的数据，获取以上4 个SNP 与结局的显著关系。合并暴露与结局数据集，并剔除了存在连锁不平衡和回文结构的SNP，最终共纳入了4 个SNP 作为工具变量进行MR 分析。

2.1.2 TSMR 研究结果 本研究的TSMR 结果提示UC 可降低HT 的发生风险。IVW 法提供了主要的因果证据（OR= 0.957，95%CI：0.924 ～ 0.990，SE= 0.017，P＜ 0.05），三种MR 统计方法的结果见表2，MR-egger 法的结果（OR= 0.936，95%CI：0.924 ～ 0.990，SE= 0.032，P＞ 0.05）虽然不显著但OR方向同样支持IVW 法的结果。Weighted median法的结果（OR= 0.944，SE= 0.021，P＜ 0.001）其OR方向也支持了IVW 法的因果判断。因此可知UC对HT 具有负向因果关系，TSMR 分析的散点图和森林图见图2、3。

表2 TSMR 分析结果Tab.2 Results of the TSMR analysis

图2 两样本孟德尔随机化散点图Fig.2 Scatter Plot of Two-Sample MR Analysis

图3 两样本孟德尔随机化森林图Fig.3 Forest Plots of Two-Sample MR Analysis

2.1.3 敏感性分析 Cochran′s Q 检验评估异质性Q= 2.566，P= 0.463，表明纳入的SNPs 无明显异质性，漏斗图也进一步验证了无异质性的结果（图4）。对于水平多效性，MR-Egger 法截距值为0.011 7（P= 0.523），表明纳入的SNPs无明显水平多效性，说明工具变量并不通过暴露以外的途径影响结局。MR-PRESSO 法未发现任何异常的SNP（P= 0.548）。留一法检验显示去除任意SNP 后结果仍然稳定（图5）。

图4 两样本孟德尔随机化漏斗图Fig.4 Funnel plot of Two-Sample MR Analysis

图5 “留一法”敏感性分析结果Fig.5 The results of sensitivity analysis of “leave-one-out”

2.1.4 工具变量评价 利用R2及F值的计算公式得出，其中4 个SNP 的统计量F值均＞ 10（31.07 ～65.43），表明本研究的工具变量均为强工具变量。

2.2 HT 对UC 的反向因果分析

2.2.1 遗传工具变量确定 同样选择与HT相关具有全基因组显著性（P＜ 5 × 10-8）的遗传变量SNP位点进行汇总，筛选出同时满足假设1、2、3 的SNP，得到70 个符合假设的SNP。利用PhenoScanner 数据库去除与结局UC 具有相关意义的SNP（n= 0），通过GWAS 提取结局UC 的数据，获取以上70 个SNP 与结局的显著关系。合并暴露与结局数据集，并剔除了2 个存在回文结构的SNP（rs2412976、rs2921053），并根据MR-PRESSO 结果剔除1 个离群的SNP（rs9271365），最终共纳入了67 个SNP 作为工具变量进行MR 分析。

2.2.2 TSMR 研究结果 本研究的TSMR 结果提示HT 对UC 不存在反向的因果关系。IVW 法（P=0.521），MR-Egger 法（P= 0.231），中位数加权法（P= 0.544）的结果均不提示存在显著因果关系（P＞ 0.05）。

2.2.3 敏感性分析 同样采用IVW 法评估异质性Q= 78.831，P= 0.133，表明纳入的SNPs 无明显异质性，漏斗图也进一步加强了无异质性的结果。MR-Egger 法截距值为-0.014（P= 0.306），表明纳入的SNPs 无明显水平多效性，说明工具变量并不通过暴露以外的途径影响结局。此外留一法检验示去除任意SNP 后结果仍然稳定。

2.2.4 工具变量评价 反向因果关系统计量F均＞ 10（30.03 ～ 0523.14），都为强工具变量。

3 讨论

在本研究中，我们通过TSMR 评价了UC 和HT之间潜在的双向因果关系。结果证实了UC 对HT的负向因果关系，表明并不推荐对无症状的UC 患者开展额外的HT 检查。此外我们识别了4 个联系UC 与HT 的SNP 位点。这些发现有助于加深对UC 与HT 之间潜在的病理机制理解。

在UC 对HT 的MR 分析结果中我们发现，MREgger 法得到了与另两种方法不一致的结果。有研究认为MR-Egger 法只能用作来评估工具变量的敏感性分析，而不能代替IVW 法的结果，并且本研究的质量评估中，工具变量并不存在异质性和水平多效性［15］。所以尽管MR-egger 法存在不一致的结果，我们仍然判断在UC 和HT 之间存在负向的因果关系。

既往的观察性研究［7-8，10，16-20］已经发现了一些UC 和HT 之间的关系。DORE 等［10］在一项回顾性病例对照研究指出，UC 患者与健康个体相比患有HT 的风险更低（OR=0.33；95%CI：0.17 ～ 0.66）。RICART 等［19］也发现IBD 患者患有自身免疫性甲状腺疾病的风险较低。此外POORAN 等［20］在一项横断面研究中发现，同属于IBD 的克罗恩病患者发生甲亢的风险与正常个体相似，但HT 的风险更低。这些研究的结果与我们的孟德尔分析结果一致，即UC 患者与HT 之间存在负向因果关系。虽然也有一些文献指出UC 会增加HT 的风险，例如KAPPELMAN 等［7］的研究中发现儿童UC 患者发生HT 的风险更高，MAHMOUD 等［8］发现IBD 患者中患有HT 的频率更高。但是BAR YEHUDA 等［21］2019 年发表的一项纳入12 625 例IBD 患者的IBD与自身免疫性疾病的研究表明，IBD 患者与大部分自身免疫性疾病相关如：风湿性关节炎、过敏性哮喘、银屑病等，但这其中并不包括甲状腺疾病。针对本研究UC 对HT 负向的因果关系解释，有研究指出针对IBD 的免疫抑制治疗可能抑制了自身免疫性HT 的产生，PASCHOU 等［22］发现IBD患者经过抗肿瘤坏死因子（anti-TNF）如阿达利姆单抗（adalimumab），英夫利昔单抗（infliximab）治疗后，与一般治疗法相比甲状腺自身抗体呈下降趋势。

本研究的优势主要在于首次运用了孟德尔随机化法研究了UC 与HT 之间的因果关系，并且我们还进行了反向TSMR 分析，以验证可能的反向因果关系。但是在我们的研究中仍然存在一些局限性。第一，本研究的选用的数据集是汇总类型的数据，无法通过年龄分层或者性别来进行亚组分析。第二，数据样本集包括了欧洲人群，这限制了对非欧洲人群的应用。需要进一步的研究验证对非欧洲人群的适用性。第三，本研究只纳入了4 个工具变量，虽然都为强工具变量，但存在工具偏倚的可能，因此需要谨慎地解释研究结果。在未来的研究中，可纳入更多的遗传标记作为工具变量，并进一步探索环境和生活方式因素对UC与HT 之间的因果关系，以及多种混杂因素对UC与HT 之间因果关系的影响。

【Author contributions】HUA Yin and WANG Zhen participated in the planning process and wrote the manuscript draft.WANG Xiaoyan performed the data analyses.XIN Yongning and LIU Shousheng interpreted the results of the data analyses and approved the final manuscript.All authors read and approved the final manuscript as submitted.

【Conflict of interest】The authors declare no conflict of interest.