整合性全转录组关联研究与孟德尔随机化分析鉴定出LYNX1和MS4A14作为骨髓炎的治疗靶点

摘要

本研究结合转录组范围关联分析与孟德尔随机化方法，系统性地鉴定骨髓炎的药物靶点。利用遗传学数据库，作者确定了LYNX1和MS4A14为关键候选基因，揭示了调控骨感染的关键分子通路，为精准治疗策略的制定提供了依据。

摘要

骨髓炎是一种严重的感染性疾病，其特征是骨组织和骨髓发生显著炎症，主要由金黄色葡萄球菌（Staphylococcus aureus，S. aureus）等细菌病原体引起。然而，目前的治疗方案因抗生素耐药性的出现以及由强效生物膜形成和胞内细菌持续存在所导致的反复感染而受到显著限制。因此，迫切需要寻找传统抗菌药物之外的新型治疗靶点。本研究采用整合性计算方法，结合转录组范围关联研究（TWAS）与基于汇总数据的孟德尔随机化（SMR）分析，系统性地识别骨髓炎背后的因果易感基因。利用大规模遗传变异作为工具变量，预测多种组织中的基因表达，从而探索这些靶基因与疾病风险之间的因果关联。该整合分析框架鉴定出LYNX1和MS4A14为关键候选基因及潜在治疗靶点。具体而言，LYNX1与骨髓炎易感性增加相关，而MS4A14则表现出潜在的保护作用。这些发现凸显了这些基因在骨感染过程中对宿主免疫应答和炎症调控的调节作用。本研究弥合了全基因组关联研究结果与生物学解释之间的鸿沟，深化了对骨髓炎遗传基础的理解，并支持发展靶向精准医疗策略，以增强宿主防御能力并克服治疗耐药性。

引言

骨髓炎是一种严重的感染性疾病，其特征是骨组织和骨髓出现显著的炎症反应。该病常由细菌侵入引发，其中金黄色葡萄球菌（Staphylococcus aureus，S. aureus）是急性和慢性血源性骨髓炎中常见的病原体^1,2。尽管S. aureus在其中占据重要地位，但骨髓炎常为多微生物感染，且由于不可培养的微生物持续存在以及生物膜动态相互作用，仅清除S. aureus并不总能根除感染^3,4。在易感宿主中，病原体采用复杂的逃逸策略，包括在巨噬细胞内长期存活以及形成致密的生物膜。这些生物膜由细胞外DNA、多糖和特化蛋白组成，形成强大的物理屏障，可抵抗宿主免疫防御并限制常规抗菌药物的渗透^5,6。因此，尽管临床治疗通常包括积极的外科清创和长期抗生素治疗，但顽固性生物膜和多微生物群落的存在显著降低了治疗效果，导致病程迁延和感染反复发作^7,8。阐明宿主-病原体相互作用背后的分子结构，对于发现新的治疗靶点仍至关重要⁹。

全基因组关联研究（GWAS）通过鉴定多个与骨髓炎易感性相关的遗传位点，推动了对该疾病的理解¹⁰。然而，由于多重检验的严格阈值以及统计效能和样本量的限制，许多影响骨髓炎易感性的关键位点可能仍未被发现⁷。此外，GWAS研究结果的解读仍具挑战性，因为大多数已鉴定的变异位于非编码区或基因间区，提示其可能通过调控基因表达而非直接改变蛋白质结构来发挥作用¹¹。同时，原始GWAS数据缺乏直接的组织特异性，因此需要结合后续的补充分析，以确定这些变异发挥作用的细胞环境。转录组范围关联研究（TWAS）通过将参考面板中基于遗传预测的基因表达数据整合到现有的GWAS数据集中，克服了上述局限性，能够在多种组织中系统性地识别出与疾病相关且具有功能意义的基因¹²。

为了加强因果推断，孟德尔随机化（Mendelian randomization, MR）采用遗传变异作为工具变量，利用受精时等位基因的随机分配来减少混杂因素和反向因果关系的影响^13,14。表达数量性状位点（expression quantitative trait loci, eQTL）分析进一步将全基因组关联研究（GWAS）中的变异与基因转录关联起来¹⁵，而基于汇总统计量的孟德尔随机化（summary-based Mendelian randomization, SMR）则整合GWAS与QTL数据集，通过异质性检验优先筛选出因果基因，并将其与由连锁不平衡驱动的关联区分开来¹⁶。

本研究首次应用了一种整合的TWAS与MR分析框架，专门针对骨髓炎，旨在填补全基因组关联研究（GWAS）发现与临床可操作性见解之间的空白。通过将实验性转录组数据与大规模人群遗传数据相结合，系统地绘制了疾病易感基因并推断其因果关系。该整合策略旨在揭示先前未被认识的、参与宿主对骨髓炎抵抗的分子通路，并识别出用于精准治疗开发的候选靶点。

方案

孟德尔随机化（Mendelian Randomization, MR）和全转录组关联研究（Transcriptome-Wide Association Study, TWAS）分析中使用的所有汇总统计量均严格来源于先前已发表的、去标识化的数据集。原始研究的伦理审批和个体知情同意情况已在各自发表的文献中注明。因此，浙江省中医院机构审查委员会（Zhe Tongde Lunshen 2024 [Yan] No. 028-JY）豁免了本项数据挖掘研究所需的额外伦理审批。本研究所用工具列于材料表中。

1. RNA-seq 数据获取与处理

为了初步验证哺乳动物物种间先天免疫通路的保守性，从基因表达综合数据库（Gene Expression Omnibus, GEO，编号：GSE272198）获取了转录组数据¹⁷。骨髓来源的巨噬细胞（BMDMs）以金黄色葡萄球菌（S. aureus）感染（感染复数，MOI = 10）1小时，随后用溶葡萄球菌酶（20 µg/mL）和庆大霉素（50 µg/mL）处理，以清除胞外细菌。用磷酸盐缓冲液（PBS）洗涤三次后，将BMDMs继续培养24小时，用总RNA提取试剂裂解细胞，并进行测序。

采用自动电泳系统评估 RNA 质量以确保其完整性。文库来自三次独立实验，并在高通量测序平台上进行测序。使用 STAR（v2.7.10a）将原始测序读段比对至小鼠基因组（GRCm38，mm10）。采用 DESeq2（v1.38.0）鉴定差异表达基因（DEGs）。为减少假阳性结果，统计学显著性定义为校正后 p 值（FDR）< 0.05 且 |log₂ 倍数变化| > 1。使用 clusterProfiler（v4.6.0）进行基因本体（GO）分析，使用 GseaVis（v0.0.5）进行基因集富集分析（GSEA）。在 R（v4.2.0）中使用 pheatmap 包（v1.0.12）生成热图。

TWAS 分析

全血RNA测序和全基因组测序（WGS）数据来自基因型-组织表达（GTEx）项目（V8）¹⁸。预训练的基因表达模型来自公共数据库（https://doi.org/10.5281/zenodo.3842289）。TWAS所需的骨髓炎汇总统计信息来自FinnGen联盟，包含2,336例病例和473,264例对照¹²。

TWAS 采用三种算法进行：联合组织填补（joint-tissue imputation, JTI）、PrediXcan¹⁹ 和 UTMOST^12,20。JTI 通过估计基因表达相似性和表观遗传染色质可及性，以优化预测准确性。PrediXcan 采用弹性网络回归并结合五折交叉验证，而 UTMOST 则利用稀疏组 LASSO 方法整合多组织表达数据以提高准确性。Zhou 等人¹² 描述的改进版 UTMOST 框架对超参数进行了标准化处理，以实现无偏估计。保留具有稳定交叉验证得分的基因作为可填补基因，其预定义标准为相关系数 r > 0.1 且预测显著性 p < 0.05²¹。全血转录组模型使用来自 1000 Genomes 参考数据集的 SNP 协方差矩阵构建。

随后分析了预测基因表达与骨髓炎风险之间的关联。为校正多重检验，TWAS 的统计学显著性主要采用错误发现率（False Discovery Rate, FDR）< 0.05 作为阈值进行定义。鉴于本多阶段研究具有假设生成的性质，对于达到提示性（名义）阈值（p < 0.05）的位点，也优先用于后续的孟德尔随机化（SMR）和共定位分析。该整合策略旨在尽可能全面地捕捉潜在的调控驱动因素，同时依赖多组学交叉验证（TWAS + SMR）以确保所筛选候选位点的稳健性。

SMR 分析

本研究遵循《加强流行病学中观察性研究报告规范》（Strengthening the Reporting of Observational Studies in Epidemiology, STROBE）指南²²。为了通过计算方法定义一种代表线粒体功能障碍遗传易感性的表型（出于分析目的，下文简称为“mitodys”），我们从 MitoCarta3.0 数据库²³ 中提取了所有已知线粒体相关基因对应的转录本。该基因集合被用作后续多基因风险预测的预定义、基于生物学知识的基础。所有与“mitodys”相关的下游功能解释均源于该计算推断，应被视为预测性结果并用于提出假设。

使用编码序列上下游1000 kb范围内的变异生成表达数量性状位点（eQTL）工具（cis-eQTLs）。汇总统计信息来源于eQTLGen联盟和GTEx V8²⁴。根据全基因组显著性阈值P < 5E-8，筛选出与1,013个线粒体功能障碍相关转录本相关的8,932,843个SNP。骨髓炎结局的基线全基因组关联研究（GWAS）统计数据来自FinnGen²⁰。

采用SMR（版本1.0.3）在默认参数下进行基于汇总数据的孟德尔随机化（Summary-data-based Mendelian Randomization, SMR）分析，以估计基因表达特征与骨髓炎结局之间的多效性关联。因果效应beta_mitodys–osteomyelitis表示线粒体功能障碍对骨髓炎的估计对数优势比效应大小，其计算公式如下：

比值比（ORs）表示标准化基因表达水平每增加一个单位的自然对数时的变化。此外，采用依赖工具异质性（HEIDI）检验进一步评估共定位情况。

结果

批量RNA测序及功能富集的表征

采用分离未感染（对照）与S. aureus感染的BMDMs的方法，在实验组中鉴定出3,775个差异表达基因（DEGs），其中包括1,686个上调基因和2,089个下调基因。GO分析显示，这些基因在先天免疫激活相关通路中显著富集，包括先天免疫应答的正向调控（GO:0045089）和激活（GO:0002218），反映了典型的由细胞因子驱动的病原体防御机制。白细胞增殖的调控（GO:0070661）突出了关键的免疫效应动态。值得注意的是，对病毒刺激的应答（GO:0009615）也呈现富集，提示可能存在交叉反应，即细菌感知激活了与病原体相关分子模式（PAMP）重叠的信号级联反应。与上述结果一致，GSEA分析进一步证实了TNF、NOD样受体和Toll样受体信号通路的激活，确立了骨髓炎中宿主-病原体相互作用的基础性炎症框架（图1A–F）。

全转录组关...

讨论

通过在跨组学数据集中整合TWAS和SMR分析，本研究鉴定出与骨髓炎表型相关的核心遗传位点网络。LYNX1和MS4A14被确定为候选基因靶点，这一发现是本研究的重要贡献，通过引入转录组学背景，提供了超越传统GWAS结果的机制性见解。这种整合性方法凸显了精准医学策略在推进严重骨感染临床管理方面的潜力。

从方法学角度来看，利用孟德尔随机化推断基因表达与骨髓炎之间的因果关系具有显著的分析优势。由于遗传变异在受精时随机分配，孟德尔随机化近似于随机对照试验的设计，能够减少观察性研究中固有的混杂因素和反向因果关系的影响³。

一个核心发现是鉴定出 MS4A14 和 LYNX1 分别可能是骨髓炎的潜在保护性因素和风险相关因素。然而，由于这些基因在骨感染中的直接作用此前尚未明确，因此功能解释仍属假设性，需谨慎解读。MS4A14 被鉴定为一个具有遗传保护作用的轴心¹⁶。尽管其在骨髓炎中的具体生理功能尚不清楚，但该基因属于四跨膜蛋白 MS4A 家族，已有文...

材料

本文使用的材料清单

姓名	公司	目录编号	评论
clusterProfiler	Guangchuang Yu's Team (Sun Yat-sen University)	NA	用于差异表达基因的基因本体（GO）功能富集分析。
DESeq2	Bioconductor Project (Open Source)	NA	用于识别差异表达基因（DEGs）。
GSEA / GseaVis	Broad Institute / Open Source Community	NA	用于基因集富集分析（GSEA），以识别共表达的基因集。GseaVis 可能用于可视化分析结果。
HEIDI（依赖工具异质性）检验	SMR 软件中的集成功能	NA	用于检验 SMR 中观察到的关联是否由连锁不平衡驱动，而非共享的因果变异。P_HEIDI > 0.05 表示无显著异质性。
JTI（联合组织填补）	Zhou 等人（开源模型）	NA	一种 TWAS 算法，通过估计基因表达相似性和表观遗传染色质可及性，以优化预测准确性。
pheatmap	Raivo Kolde（R 软件包）	NA	用于生成差异表达基因的热图，以实现数据可视化。
PrediXcan	Gamazon 等人（开源模型）	NA	一种 TWAS 算法，采用弹性网络回归结合五折交叉验证，用于填补基因表达并检验其与疾病风险的关联。
SMR（基于汇总数据的孟德尔随机化）	由复旦大学团队开发（开源软件）	NA	用于执行基于汇总数据的孟德尔随机化分析，以探索基因表达与骨髓炎风险之间的因果关联。本研究使用版本 1.0.3。
STAR	Alexander Dobin（冷泉港实验室）	NA	用于将 RNA-seq 测序读段比对至参考基因组（小鼠基因组 GRCm38/mm10）。
UTMOST	Hu 等人（开源模型）	NA	一种 TWAS 算法，利用稀疏组 LASSO 框架整合多组织表达数据，以提高预测准确性。