《Environmental Science & Technology》：Evaluating Methods for High-Dimensional Mediation in Metabolomics Data
引言
代谢组作为细胞、组织或生物体内小分子化学物质的集合，提供了基因-环境相互作用和基因转录的功能性生物学测量。尽管代谢组学已成为一种敏感的分析平台，有助于描述介导暴露-结局关系的生物学基础，但近期研究指出，现有方法往往未能充分处理中介物之间的相关性，当中介物数量远超过样本量时扩展性差，并可能产生保守或有偏的估计。改进这一理解将有助于阐明暴露、代谢组与疾病之间的复杂关系。高维技术的例子包括最初为表观遗传学等领域开发的中介分析方法，例如HIMA和HDMA，以及“中间相遇”（MITM）方法。
全氟和多氟烷基物质（PFAS）与许多生物扰动相关，包括氨基酸、脂肪酸、甘油磷脂、甘油脂、磷酸鞘脂、胆汁酸、神经酰胺、嘌呤和酰基肉碱。高维中介方法已被用于识别哪些代谢物和代谢通路可能 underlie 环境污染物（如PFAS）暴露与不良健康结局（包括胎儿生长受损）之间的联系。这些初步研究凸显了高维中介在代谢组学中的可行性和潜力。例如，产前PFAS暴露与支链氨基酸水平改变相关，这可能通过扰动氨基酸和能量代谢导致胎儿生长受损。中介分析提供了一个框架来正式检验此类通路，使研究人员能够评估支链氨基酸或其他代谢物类别是否在统计上介导了PFAS暴露与出生结局之间的关系。量化中介效应可以从统计上了解暴露-结局关系中由特定代谢通路解释的比例，从而更清晰地理解生物学机制和潜在的干预靶点。
中介变量是一个受暴露影响并沿着因果通路影响结局的变量。经典因果中介分析，如Baron和Kenny最初描述的，提供了一个识别介导暴露与结局之间关联的变量的框架。这可以通过暴露对结局的直接效应和间接效应来描述。间接效应指的是暴露通过某个第三方中间变量（即中介变量）对结局产生的效应，而直接效应指的是暴露对结局产生的、不依赖于中间变量的效应。当使用数据估计这些效应时，将直接和间接效应估计值解释为因果所需的条件相当强。一个关键假设是条件可交换性，这要求（1）暴露与中介变量之间、（2）暴露与结局之间以及（3）中介变量与结局之间没有未控制的混杂。使用经典因果中介方法时所需的第四个额外假设是，不存在受暴露影响的中介-结局混杂因素。如果违反这第四个假设，则必须使用基于一系列“g方法”的更复杂技术。否则会出现几个问题，包括潜在的碰撞分层偏倚。
不幸的是，由于生物过程的复杂相互影响以及代谢组数据中每位参与者固有的数千个已识别特征，经典中介方法无法应用于高维中介设置。例如，一项最近研究多氯联苯（PCBs）对代谢组影响的研究，在同一生化级联中识别出21种代谢物。许多已识别的代谢物表现出强烈的相互依赖性；例如，4-氨基丁酸（GABA）是谷氨酸的直接前体。如果假设使用现有的高维中介方法将这些代谢物与健康结局联系起来，那么这些方法所要求的因果假设将被违反。具体来说，这些方法要求中介变量独立作用且无混杂，然而实际上，几乎不可能假设某个给定的中介变量（例如GABA）没有同时也受暴露影响的中介-结局混杂因素（例如谷氨酸）。用于检验高维中介的适当统计方法才刚刚开始被探索。在研究人员可用的中介方法中，有些方法在处理复杂的代谢组学数据时可能表现优于其他方法。因此，识别和探索不同的中介选项，并评估它们在代谢组数据中的性能，将提高我们对使用何种方法的理解。
模拟研究在统计研究和方法学发展中起着至关重要的作用，它通过应用与真实世界数据高度相似的计算机生成数据来实现。这些研究为研究人员提供了一个宝贵的工具，用于评估统计方法在特定数据背景下解决特定研究问题的适用性。在复杂的高通量数据和先进机器学习算法应用的时代，模拟研究变得更加重要。通过创建具有已知特征的合成数据集，研究人员可以系统评估不同方法在各种条件下的性能。这使研究人员能够评估这些方法在揭示模拟“真相”或估计目标方面的能力，而这在真实世界数据中是无法实现的，因为真实值未知。对这些方法行为和性能的洞察使研究人员能够识别不同方法的优势和局限性，从而改进方法学并预防未来研究中潜在的可重复性问题。模拟研究为研究人员提供了一个受控的环境来探索、改进和验证统计方法，从而增强了实际数据应用中分析的可靠性和稳健性。
目前，对于非靶向代谢组学数据如何进行中介分析尚无标准方法。许多复杂因素包括高维性、中介集合通常大于样本量以及复杂的生物相互作用。本研究的目标是通过比较三种已在代谢组学数据和其他高维环境中应用的中介方法，并检验它们在模拟代谢组学数据中的性能，来填补这一关键空白。这些方法包括HIMA、HDMA和MITM方法。
方法
进行了一项模拟研究，以增进对高维中介分析方法在代谢组学数据中表现的理解。本研究使用基于亚特兰大非裔美国人母婴队列（Atlanta AA cohort）中观察到的代谢组学特征模拟的数据，以了解这种高维中介物的变化如何影响各种方法估计通过所有代谢特征的总间接效应（TIE）、通过单个特征的成分间接效应（CIE）以及所捕获代谢物的敏感度和特异度的能力。TIE代表由中间变量介导或解释的总效应部分。TIE是总效应的一部分，总效应包括直接效应（暴露对结局的、不经由中间变量介导的效应）和间接效应（暴露对结局的、通过中间变量运作的效应）。CIE允许理解单个特征对总间接效应的贡献程度。通过将总效应分解为其直接和间接组成部分，研究人员可以理解暴露、中间变量和结局之间的整体因果关系。
高维中介分析
本研究评估了三种高维中介分析方法：HIMA、HDMA和MITM。这些方法没有考虑受暴露影响的中介-结局混杂，这是应用于真实世界代谢组学数据时的一个关键考虑因素，因为真实的代谢组学数据中复杂的反馈回路和暴露诱导的混杂很常见，并且可能使中介效应估计产生偏倚。
- 1.
HIMA是一种基于惩罚回归的方法，采用自多中介模型的框架，并将此方法扩展到高维设置（即中介变量数量大于样本量的情况）。HIMA在表观遗传数据中已证明有效，但尚未在代谢组学数据中进行正式评估。HIMA可输出TIE和CIE。简而言之，HIMA遵循三个一般步骤：（1）将潜在中介因子池减少到小于样本中个体数量的规模；（2）使用最小最大凹惩罚（MCP）进行变量选择；（3）进行联合显著性检验以评估中介效应。在第1步中，为了将中介变量数量减少到更易管理的水平，对中介-结局模型应用确定独立筛选（SIS），以选择对结局变量影响最大的中介变量。简而言之，SIS是一种变量选择技术，通过根据协变量与结局的边际关联对它们进行排序来降低高维数据的维度。值得注意的是，中介变量将被标准化以确保系数代表相同的尺度。此步骤将中介变量数量减少到由研究人员控制但小于样本量的值。在第2步中，将选定的候选中介变量在中介-结局模型中建模，依赖MCP将中介-结局效应向零收缩。Zheng等人选择MCP作为首选的收缩方法，因为与其他方法（例如HDMA中的LASSO）相比，它能以计算高效、准确且无偏的方式进行变量选择。在第3步中，HIMA使用最大化p值来评估中介效应的显著性，并使用Bonferroni校正对多重检验的结果进行调整。每个选中中介变量的CIE计算为暴露到中介变量的系数与中介变量到结局的系数的乘积，TIE是所选中介变量的CIE之和。我们使用了R包hima中的默认HIMA实现。
- 2.
HDMA是一种基于去稀疏LASSO估计量的高维中介技术。HDMA可输出TIE和CIE。首先，应用SIS程序来降低中介集合的维度，使其小于参与者数量，使用标准：d = n / log(n)，其中d是中介集合，n是样本量。接下来，应用惩罚回归模型。该模型引入了一个惩罚项，有效地将中介集合中较不重要变量的系数收缩至0。惩罚线性结局模型的目标是找到与结局变量具有最强关联的预测变量集合，同时最小化不相关或弱相关变量的影响。为了抵消可能引入系数估计的任何偏差，应用去偏程序以试图提供更准确的估计。Gao等人确定HDMA方法比其他选项在统计上更稳健和高效，因为它在一个回归步骤中拟合多个中介变量，而不是像其他高维中介选项（例如HIMA）那样一次测试一个中介变量的中介效应。HDMA中的CIE计算方式与HIMA类似：作为每个选中中介变量的暴露到中介变量和中介变量到结局系数的乘积。TIE是所选CIE的总和。我们使用了R包hdmed中可用的默认HDMA实现。
- 3.
MITM方法结合非靶向代谢组学分析和靶向代谢组学分析，以揭示与感兴趣暴露和感兴趣结局相关的特征。首先，分别对暴露与代谢组之间以及代谢组与结局之间进行独立的代谢组全关联研究（MWAS）。识别重叠的特征和代谢通路，并将其视为与暴露和结局相关的交叉信号的代表。将使用传统系数乘积方法的一个简单扩展来确定每个已识别特征的TIE和CIE。首先，使用边际线性回归对暴露对每个特征的影响以及每个特征对结局的影响（控制暴露）进行建模。TIE将计算为CIE的总和。CIE将计算为两步回归过程系数的乘积。该方法在R中实现。
数据生成
使用有向无环图（DAG）模拟数据。最外生变量（即指向该变量的箭头最少的变量），即混杂因素集，首先被模拟，然后是暴露、代谢物，最后是结局。使用以下方程引入相关依赖关系和分布：
C = C0，其中C0是从标准正态分布N(0,1)生成的混杂因素列表。
E = E0+ 0.2Ci，其中E0是从标准正态分布N(0,1)生成的暴露，Ci是混杂因素数量1···i。选择0.2的beta值是为了在所有测试场景中提供一致、中等大小的效应量。
代谢物可以是独立的或相关的。如果是独立的，每个代谢物从标准正态分布N(0,1)生成。作为相关代谢物生成的基础，我们从多元正态分布开始抽样，均值为0，标准差从正定均匀相关矩阵生成，代谢物之间的相关性范围从-0.42到0.45。在所有生成的代谢物集合中，选择一个子集来介导模拟暴露与结局之间的关系。这是通过将用于生成代谢物的正态分布的均值从0修改为Mj= Mj+ βE来实现的，其中Mj是中介变量数量1···j，β可以在数据生成过程中由用户定义。
然后结局从线性回归模型生成，定义为O = O0+ 0.8E + 0.01Ci+ βMj，其中Ci是混杂因素数量1···i，Mj是中介变量数量1···j。如果确定某个代谢物是中介变量，则在数据生成过程中β定义为非零；否则，beta值设置为0。
模拟场景
模拟旨在理解高维中介技术在非靶向代谢组学环境中的行为，这已在先前研究中应用。下面描述的模拟场景考虑到了这种情况。我们使用了1000次蒙特卡洛模拟来评估每种方法的性能。
为了理解样本量对不同高维中介技术的影响，我们模拟了250个体的小样本量、500个体的中等样本量和1000个体的大样本量。选择这些样本量是因为它们反映或大于非靶向代谢组学分析中通常所见的情况。这可以展示随着该领域的成熟，当功效增加时这些技术的效率和准确性可能如何变化。中介集合的大小变化为200、400和600，以理解代谢物数量的变化如何影响技术的行为。接下来，研究表明真正中介变量的数量是实际识别出的代谢物的一小部分。为了反映这一点，将2%、5%和10%的代谢物设置为真正中介变量。接下来，将暴露到中介变量和中介变量到结局关系的β值设置为0.1或0.3，以反映代谢组学研究中观察到的相对较小的beta值，并观察这如何影响每种方法识别中介集合的能力。最后，所有场景都在独立和相关的中介集合下运行，以理解违反可交换性假设4如何影响HIMA、HDMA和MITM。为了进行模拟，我们使用了由戴尔Power Edge系统组成的高性能计算集群，配备英特尔至强处理器，总共32个核心和100GB内存。模拟相关的高维代谢组学数据集，并运行每个场景的1000次重复，使用标准笔记本电脑的内存和处理能力是不可行的。使用集群使我们能够高效完成计算，而不限制个人计算可用性。请注意，本研究中检查的中介分析工具在标准硬件上运行相当，当应用于真实数据集时并不固有地需要HPC环境。
为了提高图的清晰度，正文中仅呈现了最大中介集合（p = 600）的结果，因为它最接近非靶向代谢组学中常见的维度，并提供了方法性能的保守表示。p = 200和p = 400的结果包含在支持信息中，因为观察到跨中介集合大小的类似性能趋势。
分析
本研究利用了几个关键评估指标，包括CIE、TIE以及每种方法的敏感度和特异度，以全面评估不同分析方法的表现及其对估计量准确性的影响。
TIE可以使用潜在结果来定义，这是支持本项工作的因果理解。将X定义为暴露，Y定义为结局，M定义为中介变量。TIE可以定义如下：
E[YxMx – YxMx] = ΣmE[Y | x, m] {P(m | x) – P(m | x)}
方程左边，E[YxMx – YxMx]，代表当暴露（x）保持不变，而中介变量（M）处于不同暴露水平（x或x）时潜在结果（Y）的差异。在方程右边，表达式ΣmE[Y | x, m]是当暴露（x）和中介变量（m）都设定为特定水平时期望的结果，并对所有中介变量值求和。该求和由概率{P(m | x) – P(m | x)}加权，该概率代表中介变量（M）在两种暴露情景（x或x）下取特定值（m）的概率差异。
由于TIE捕获了由暴露和结局之间代谢物集合解释的所有信息，相关结构不影响该估计目标的因果解释。然而，CIE仅在中介变量独立的假设下才能被解释。如果数据中存在依赖性，隔离单个组分（例如M1）的效应将需要以下因果比较：
E[YxM1xM2x···Mjx – YxM1x*M2x···Mjx]
如果M1依赖于任何其他中介变量，则这种对比是不可能的，因为来自暴露的信息表明暴露状态是暴露的（x*），但来自所有其他中介变量（M2x···Mjx）的信息表明暴露状态是未暴露的（x）。这种不可能的悖论使得因果对比无法解释，这是由于“翻供证人”问题。为了处理这个问题，我们仅在假设中介变量之间独立的场景中评估CIE。
为了评估中介方法HDMA、HIMA和MITM的性能，我们评估了它们从整体集合中识别真正中介变量的敏感度和特异度。敏感度定义为每次模拟中正确识别的中介变量数量除以真正中介变量的总数。特异度定义为每次模拟中正确拒绝的非中介变量数量除以非中介变量的总数。这些指标使我们能够定量比较每种方法在模拟代谢组学数据中准确识别真正中介变量和正确拒绝非中介变量方面的有效性。分析在R（版本4.4.1）中进行。
为了评估每种中介方法的性能，我们为每个中介集合大小水平（p = 200, 400, 600）模拟了36个场景。这些场景代表了以下参数的所有组合：中介变量的相关结构（相关或独立）、样本量（n = 200, 500, 1000）、真正中介变量的比例（p的2%、5%或10%）以及暴露到中介变量和中介变量到结局关系的效应大小（β = 0.1或0.3）。场景1–18代表相关中介变量设置，而场景19–36对应独立中介变量设置。此外，场景1–9和19–27使用较小的中介变量beta值（β = 0.1），而其余场景使用较大的beta值（β = 0.3）。偏差计算为跨模拟重复的平均估计效应与数据生成过程中指定的相应真实效应之间的差异。这种结构使我们能够系统地评估每种方法在不同数据维度和效应大小条件下的表现。
结果
我们呈现了关于方法在估计CIE、TIE以及敏感度和特异度方面表现的模拟结果。在每个部分内，首先呈现具有独立中介结构的场景结果，然后是具有相关中介结构的场景结果（如果适用）。
成分间接效应（CIE）
CIE仅在独立代谢物场景中进行评估。对于所有测试的、假设代谢物独立的模拟场景，HDMA、HIMA和MITM之间的结果没有显着差异。总体而言，在独立代谢物设置中，HDMA和HIMA对CIE都是无偏的。相比之下，MITM在较小beta（β = 0.1）的设置中似乎系统地高估了CIE，而在较大中介变量beta（β = 0.3）的设置中系统地低估了CIE。中介变量beta指的是数据生成过程中用于表示暴露到中介变量和中介变量到结局关系强度的固定值。MITM方法估计CIE的偏差幅度在较小beta设置中范围为-0.001到0.058，在较大beta设置中范围为-0.075到-0.013。
总间接效应（TIE）
在代谢物独立和相关场景中都检查了总间接效应。在独立场景中，使用MITM方法通常低估了TIE。MITM方法的偏差幅度范围为-5.233到0.019。对于HIMA，在较低中介变量beta设置（β = 0.1）中TIE被低估，但在β = 0.3时通常有所改善。HIMA方法估计TIE的偏差幅度在较小beta设置中范围为-0.547到-0.027，在较大beta设置中范围为-3.43到0.171。HDMA在大多数独立场景中能够估计TIE，除了场景25和34在p = 400和600时，这两个都是高维设置。HDMA方法的偏差幅度范围为-3.109到0.046。在所有独立场景中，当模拟中介变量数量（p）小于样本量（n）时，即低维模拟设置中，方法表现更好。
在相关场景中结果相似，MITM低估了TIE，并且随着代谢物数量的增加，性能恶化。MITM方法的偏差幅度范围为-5.113到-0.005。对于HDMA和HIMA，与独立场景相比，在较低中介变量beta设置（β = 0.1）中，样本量（n）对TIE低估的影响较小。当β = 0.3时，HIMA表现最一致，尽管较小的样本量和高维场景（即当p > n时）导致该方法低估TIE。在这种设置下，与HIMA相比，HDMA似乎更受相关数据结构的影响，导致估计准确性较差。HIMA方法在相关设置中估计TIE的偏差幅度范围为-3.751到0.07，而HDMA在相关设置中的偏差幅度范围为-3.632到-0.014。
敏感度和特异度
敏感度范围很广，取决于测试的场景（捕获的真正中介变量<25%到100%）。所有方法的敏感度都随着样本量增加、中介变量beta值增大以及在低维设置中而提高。相比之下，所有方法在所有场景中的特异度始终很高（>90%）。MITM的敏感度往往低于HIMA和HDMA，尤其是在较低中介变量beta设置中（<60%），但通常具有最高的特异度（>97%）。HIMA和HDMA的敏感度往往相当，但HIMA的特异度始终高于HDMA，尽管两种方法都高于90%。
讨论
鉴于代谢组学在流行病学和环境健康研究中的应用迅速增长，识别高效准确的中介方法来研究介导暴露-结局关系的生物学基础至关重要。在我们的研究中，我们比较并评估了HIMA、HDMA和MITM方法在模拟代谢组学数据中的表现，以系统检验它们在不同场景下的性能，特别是当中介集合存在相关性以更好地模拟代谢组学数据固有的复杂生物相互作用时。在相关设置中，方法无法估计CIE，因为没有明确的因果对比，并且方法往往低估TIE。在较低中介变量beta设置（β = 0.1）中，样本量小和/或潜在中介代谢物数量多时，敏感度降至约50%或更低。然而，虽然许多真正的中介代谢物未被识别（反映在低敏感度上），但所有方法都稳健地识别出真正的中介代谢物，因为所有方法的特异度都>90%。
基于这些结果，有兴趣在代谢组学数据中使用HIMA、HDMA和MITM方法进行高维中介分析的研究人员，可能希望考虑将多种方法纳入其工作流程。尽管HIMA表现出最强的整体性能，可能是首选方法，但与MITM相比，它的实现和解释更为复杂。因此，将MITM等方法与HIMA结合使用可能是有利的，因为它使研究人员能够识别和交叉验证关键代谢物，突出显示两种方法都发现的代谢物
