将组学数据集与OmicsPLS包集成

摘要

背景

随着可用生物医学数据的指数级增长，需要能够提取数据集之间关系信息的数据集成方法。然而，这些数据集可能具有非常不同的特征。对于可解释的结果，需要量化特定数据的变化。为此，提出了双向正交偏最小二乘法（O2PLS）。为了促进该方法的应用和发展，需要免费的开源软件。然而，O2PLS并非如此。

结果

我们介绍OmicsPLS公司是R中O2PLS方法的开源实现。它可以有效地处理低维和高维数据集。实现了检查和可视化结果的通用方法。可以使用标准和更快的替代交叉验证方法来确定组件的数量。一项模拟研究表明，与替代方案相比，OmicsPLS在准确性和CPU运行时间方面具有良好的性能。我们通过整合遗传和血糖数据来证明OmicsPLS。

结论

我们提出了OmicsPLS R包：用于统计数据集成的O2PLS的免费开源实现。OmicsPLS可在https://cran.r-project.org/package=OmicsPLS并且可以通过安装在R中install.packages（“OmicsPLS”）.

背景

随着技术的迅速进步，可以测量出生物变化的几个水平。因此，同一组受试者可获得多个组学数据集。为了更好地理解潜在的生物系统，应该同时分析这些数据[1].

已经开发了几种数据集成方法，可以在忽略某些特定数据特征的情况下估计接头零件。例如，偏最小二乘法（PLS）[2]最大化协方差以计算联合主成分。典型相关分析（CCA）[三]考虑相关性而不是协方差。其他几种方法对数据集的级联版本执行分析，例如同步组件分析（SCA）[4]. 对于许多数据集成方法，都可以使用开源软件包[5]. 特别是，mixOmics R包实现了PCA、PLS和CCA的几种变体[6].

Omics数据集可能是异构的，因为它们在数据特定的特征上通常不同，例如大小、规模、分布和实验误差。这妨碍了对这些数据之间联合部分的估计。为了正确解释数据集成结果，应建模特定于数据的变化[7]. 此变体捕获特定于每个数据集的信息。此外，它可能会扭曲对估计接头部分的解释[8]. 因此，我们考虑评估联合成分和特定成分的方法。这些方法包括双向正交PLS（O2PLS）[8]、JIVE[9]和DISCO-SCA[10]. O2PLS考虑两个数据集，并在联合、特定和残差部分进行分解。通过最大化两个数据集之间的协变来计算联合部分，同时校正数据特定的变化。JIVE在多个数据集的连接上使用迭代PCA来交替查找关节和数据特定部分。DISCO-SCA执行SCA并旋转解决方案，以获得每个数据集的联合组件和特定组件。

在JIVE和DISCO-SCA方法中，关节和特定组件被约束为相互正交。此外，他们假设数据集共享完全相同的联合潜在变量。O2PLS仅强制每个部分内的分量的正交性，并假设每个数据集的相关联合潜在变量。因此，我们期望O2PLS在复杂情况下具有更好的性能。

O2PLS在软件包SIMCA中实现[11]，这是封闭来源和商业。源代码的不可用阻碍了方法的开发和扩展。据我们所知，没有免费的开源替代方案可以实现O2PLS。因此，我们建议使用OmicsPLS，这是一个免费的开源R软件包，用于将两个数据集分解为联合的特定部分。关于其他方法，DISCO-SCA[12]只能从商业计算环境MATLAB中获得，而JIVE在r.JIVE包中免费提供[13]. 因此，我们将OmicsPLS与r.jive进行了比较。

我们的目标是提供方便的方法和可视化工具，并促进更先进方法的开发。文章的其余部分组织如下。首先，我们详细讨论了OmicsPLS的实现。其次，OmicsPLS软件包使用克罗地亚人群队列的遗传和聚糖数据进行说明。我们还将JIVE应用于这些数据。出于数据分析的动机，我们进行了一项模拟研究，以比较OmicsPLS和r.jive在估计精度、执行时间和针对特定数据特征的鲁棒性方面的差异。最后，我们讨论了OmicsPLS的未来扩展。

实施

O2PLS模型

将观测数据收集到矩阵中X（X）=[x个₁,…,x个_第页] (N个×第页)和矩阵Y（Y）=[年₁,…,年_q个] (N个×q个). 在这里，N个表示受试者人数，以及第页和q个表示中的变量数X（X）和Y（Y）分别是。O2PLS方法分解X（X）和Y（Y）在两个接头、特定部件和残余部件中。接头部分的尺寸由n个，每个特定零件的尺寸由n个_X（X）和n个_Y（Y）分别是。连接部分由矩阵组成T型,单位（两者都是N个×n个),W公司(第页×n个)和C类(q个×n个). 矩阵T型和单位被称为联合得分或联合潜在成分，以及矩阵W公司和C类被称为联合载荷或联合主成分。这些连接部分表示X（X）和Y（Y）.具体部分由矩阵组成T型_Y（Y）⊥(N个×n个_X（X）),单位_X（X）⊥(N个×n个_Y（Y）),P（P）_Y（Y）⊥(第页×n个_X（X）)和P（P）_X（X）⊥(q个×n个_Y（Y）). 这些矩阵分别称为特定分数和载荷。剩余部分表示为E类(N个×第页)和如果(N个×q个). 然后，O2PLS分解为

$$\开始{对齐}{X}&=TW^{\top}+T_{Y}\bot}P_{Y}\fot}^{top}+{E}，\\{underset{Data}}{\下大括号{Y}}}&={undersset{Joint}}{\undersbrace{UC^{top{}}}+{underset{Specific}{{underbrace{U_{X}\bot}P_{X}\bot}^{\top}}}+{\underset{{Residuals}}{\underbrace{{F}}}}。\结束{对齐}$$

(1)

每行X（X）和Y（Y）包含上的测量值相同的主题。在本文中，假设X（X）和Y（Y）以零为中心。之间的关系T型和单位由线性模型给出单位=T型B类_T型+H（H）或T型=单位B类_单位+H（H）^′.给，B类_T型和B类_单位是大小方阵n个，表示两个模型的回归系数。由于O2PLS算法在X（X）和Y（Y）.

注意，在PLS中，每个数据集只考虑一个接头和一个残余零件。任何特定于数据的变化都会被这两部分吸收。这使得对PLS结果的解释更加困难，因为估计的载荷可能有偏差，并且联合评分之间的相关性通常看起来较弱。O2PLS限制接头载荷W公司和C类和具体得分T型_Y（Y）⊥和单位_X（X）⊥具有正交柱。JIVE和DISCO-SCA还限制了矩阵的列[W公司P（P）_Y（Y）⊥]和[C类P（P）_X（X）⊥]为正交。此外，这两种方法都假设单位=T型，而O2PLS仅假设单位和T型.

估算O2PLS组分的O2PLS算法见[8]. 简而言之，协方差矩阵的奇异向量X（X）^⊤Y（Y）已计算。根据这些向量，估计包含联合变异和特定变异的负荷和分数。然后，使用SVD估计特定变化并从数据中减去。最后，利用修正后的数据，对接头零件进行了重新估算。

解释在每个部分中，组件具有与PCA类似的解释。尤其是荷载值周_jk公司表明了变量 x个_j个对于组件k个.如果周_jk公司和\（幻影{\dot{i}\！}w{j^{prime}k}\）具有相同的符号，对应的变量x个_j个和\（幻影{\dot{i}\！}x{j^{prime}}\）组分内正相关k个其他部分也有相同的解释。分数可用于定义学科在每个组件中：例如，如果\（幻影{\dot{i}\！}t{ik}\近似于t{i^{prime}k}\），然后是受试者我和我^′组件相似k个.连接部件之间，在k个'第个接头组件，荷载值周_jk公司和\（幻影{\dot{i}\！}c{j^{prime}k}\）指示之间的相关性x个_j个和\（幻影{\dot{i}\！}y{j^{prime}}\）较高的正或负荷载值分别表示这些变量之间的该分量具有较高的正相关性或负相关性。因此，接头荷载值周_k个和c（c）_k个可以对中的变量排序X（X）和Y（Y）基于高协变量。

实施

OmicsPLS中的功能可以组织如下

• 交叉验证：用于确定O2PLS组件数量的函数。

• 拟合：用于拟合O2PLS模型的功能。

• 总结和可视化：总结和可视化结果的功能。

交叉验证。交叉验证是一种众所周知的技术，用于选择模型的调整参数，同时限制过盈风险。所有样本都划分为k个块（表示为褶皱），模型拟合在k个−1倍。左侧折叠用于评估模型拟合。对于O2PLS，确定分量数量的一种方法是在可能整数的三维网格上最大化预测误差，并选择三元组(n个,n个_X（X）,n个_Y（Y）)从而使该误差最小化。由于O2PLS在中是对称的X（X）和Y（Y），两个预测误差之和\（||Y-\hat{Y}||^{2}+||X-\hat}X}|||^}2}\）作为误差度量。在这里，\（||A||^{2}:=\sum_{ij}A_{ij}^{2{）。此方法在交叉阀o2m功能：

交叉值_o2m（X，Y，a，ax，ay，nr_folds）

在这里，一,斧头和是的是要考虑分量数的整数向量n个,n个_X（X）和n个_Y（Y）.矢量一必须有积极因素，而两者斧头和是的可以包含零。折叠次数由数量_文件夹并且应该介于2和之间N个. The交叉阀o2m函数返回带有预测错误的三维数组。

三维网格上的交叉验证可能需要大量计算，尤其是对于许多网格点。因此，我们提出了另一种交叉验证程序[14]. 这种方法的基本原理在于对特定部分的解释：数据中的特定变化将影响联合得分，从而减少T型和单位。校正特定变化将增加此协方差。另一方面，过度纠正会再次降低联合得分之间的协方差。候选人n个_X（X）和n个_Y（Y），给定n个，是联合得分协方差最大化的那些整数。这种方法被称为：

交叉值_o2m_adjR2（X，Y，a，ax，ay，nr_folds）

它在一维网格上执行交叉验证一，同时最大化联合得分之间的协方差T型和单位在给定的二维网格上斧头和天最后一个最大化不涉及交叉验证。因此，循环数量_文件夹二维中省略了折叠。这可以大大减少计算时间，同时通常会产生与完全交叉验证方法类似的最小化结果。输出是一个包含预测误差和组件数量的矩阵(n个,n个_X（X）,n个_Y（Y）).

请注意，这两种交叉验证策略可以结合使用：替代交叉验证用于找到预测误差的候选最小值。基于这些最小化器，构建了一个三维网格，在该网格上执行全面交叉验证。两种交叉验证实现都支持并行计算。

适合的在最简单的形式中，用于拟合O2PLS模型的函数调用是o2m（X、Y、n、nx、ny）输入参数是两个数据矩阵X（X）和Y（Y）以及接头中组件的数量，X（X）-具体和Y（Y）-特定部分。输出是一个包含分数和加载的列表，其符号为[8]以及解释方差和残差矩阵的比例；这些比例定义如下。用户可以通过添加striped=真作为参数，丢弃剩余矩阵（并减少内存使用）。

默认情况下，协方差矩阵的奇异值分解（SVD）X（X）和Y（Y）用于计算接头和特定零部件。如果两者都有X（X）和Y（Y）是高维的，协方差矩阵X（X）^⊤Y（Y）将使用大量内存。因此，OmicsPLS包中实现了另一种算法，名为NIPALS[2]. NIALS算法是一种迭代算法，避免了协方差矩阵的构建和存储。此外，如果NIALS迭代次数足够大，则基于NIALS的联合分量在数值上等于基于SVD的PLS分量（直到符号）。在这种情况下第页或q个不是太大，NIPALS方法可能比SVD方法慢。因此，对数据维度进行检查以确定适当的方法。默认情况下，阈值为第页=q个=3000，可以调整。

总结和可视化建模变化总结如下摘要（对象）.给，对象包含由200万呼叫。输出包括以下比例：

• 变化X（X）和Y（Y）由接头、特定零件和残余零件解释，例如||T型W公司^⊤||/||X（X）||.

• 变化单位和T型这是可以预测的T型相应的。单位例如||T型B类_T型||²/||单位||².

请注意Y（Y）通过X（X）就是那个时候||T型B类_T型||²/||单位||²×||单位C类^⊤||²/||Y（Y）||²=||T型B类_T型||²/||Y（Y）||².

OmicsPLS包提供了一个灵活的框架来绘制每个组件中的加载。由于该框架是在ggplot2包上构建的，因此可以添加几个绘图层，以增强可视化并帮助解释结果。构建绘图的命令是

plot（x，loading_name）。

在这里x个是O2PLS配合，是唯一需要的对象。参数加载_名称表示应绘制四个部分中的哪个部分（X接头、Y接头、X特定或Y特定）。plot命令调用地理文本来自ggplot2包。其文档包含有关编辑的信息，例如文本颜色、透明度和大小。这些属性可以在OmicsPLS中更改情节功能。

工作流程和教程图中提供了OmicsPLS分析的工作流1工作流程中的步骤基于下一节所示的遗传和糖组学数据分析。此外，还有一个在线补充教程，用免费提供的转录组学和代谢组学数据说明了OmicsPLS（附加文件1).

图1

OmicsPLS包的工作流。首先，对每个数据集进行预处理。其次，O2PLS用于将每个数据集分解为联合、特定和残差部分。最后，对输出进行可视化和解释

结果

遗传和血糖数据分析

我们认为第页=333858基因型单核苷酸多态性（SNP）和q个=20定量IgG1聚糖（糖肽）丰度，用纳米LC-ESI-MS测量N个=885名CROATIA_Korcula队列参与者[15]. 这两个数据集都包含高度相关的测量值，并且是异质的（因为它们在规模、分布和测量误差方面不同）。

我们的目的是通过测定IgG1聚糖数据和遗传数据之间的统计重叠来研究IgG1-聚糖如何与遗传变异相关，如等式(1). 为此，我们使用OmicsPLS软件包获得联合变异量的估计值，并估计遗传和IgG1聚糖测量值对该联合变异的贡献。

通过提取每个基因（在UCSC hg18数据库中）50千碱基范围内的所有SNP，并应用主成分分析，对SNP进行总结。对于每个基因，相应的SNP组被尽可能多的主成分所取代，以解释至少80%的SNP。这为每个基因提供了一个包含一个或多个变量的新数据集，我们称之为遗传PC。这个“遗传PC”数据集包含37819个变量，称为X（X）.聚糖测量值是对数转换、分批修正的[16]和分位数规范化[17]. 生成的数据矩阵称为Y（Y）.

屏幕图X（X）X（X）^⊤,Y（Y）^⊤Y（Y）和X（X）^⊤Y（Y）如图所示2通过在这些碎石图中识别弯头，可以确定接头和特定组件的数量。根据这些图，保留了5个联合和5个遗传特异性成分。请注意，未检测到聚糖特异性部分。O2PLS安装耗时约5秒。

图2

遗传和聚糖数据的协方差矩阵的特征值。每个特征值对所有特征值之和的相对贡献如遗传PC所示（面板一)和IgG1聚糖数据（面板b条)，及其协方差（面板c（c）)分别为

关于五种IgG1聚糖联合成分，它们占总IgG1-聚糖变异的96%。用遗传PC可以预测的IgG1变异量为70%。每个IgG1聚糖变量的负载值如图所示三第一个关节成分与“平均”IgG1聚糖成比例，因为所有聚糖的负载值大致相同。第二个关节成分区分岩藻糖基化（负负荷值）和非岩藻糖基化（大多为正负荷值）IgG1聚糖。该组分被称为“岩藻糖基化”组分。第三种关节成分尤其涉及非半乳糖基化（负负荷值）和二-半乳糖基化（正负荷值）IgG1聚糖，而单-半乳糖基化聚糖的估计负荷值约为零。该成分被称为“半乳糖基化”成分。在第四个关节成分中，G1NS和G2NS聚糖具有高负载值。除G1NS和G2NS外，第五个关节成分区分存在二等分GlcNAc（负荷载值）或不存在（正荷载值）的聚糖。

图3

用OmicsPLS R包获得的Genetic-Glycan联合主成分。每个IgG1聚糖变量的加载值按组件（面板）描述一-e（电子）). 颜色和形状代表聚糖的生物分组。在最后一行和最后一列中，显示了特定聚糖结构的图形表示（面板（f）)

关于遗传PC数据集中的五个联合成分，它们占总变异的0.8%。对于特定部件，该百分比为1.9%。每个遗传PC联合成分中的前五个基因如表所示1在第一个遗传PCs联合成分中，负载值最高的基因为DNAJC10公司相应的蛋白质参与识别和降解错误折叠的糖蛋白。第一个关节成分对应于第一个聚糖关节成分中的“平均”聚糖模式。与“岩藻糖基化”成分相对应的第二个关节成分中的顶部基因是FUT8（FUT8）编码岩藻糖基转移酶，对岩藻糖转化为糖肽进行催化。在第三个关节成分中，对应于“半乳糖基化”成分，基因AKAP9公司具有第二高加载值。它编码一种A激酶锚定蛋白，参与维持高尔基体的完整性。注意，高尔基体中发生糖基化（尤其是半乳糖基化）。在第四和第五组分中，没有发现直接相关的基因。需要更多的研究来进一步阐明这些关系。

表1遗传甘氨酸联合主要成分的前5个基因和负载值

为了进行比较，r.jive也应用于数据。然而，该算法在500次迭代（和3000秒）后没有收敛。我们将在模拟研究中调查可能的原因。

模拟研究

进行了仿真研究，比较了r.jive和OmicsPLS的精度和速度。为了深入了解r.jive的鲁棒性，研究了r.jive缺乏收敛性的可能原因。模拟数据遵循满足O2PLS和JIVE假设的模型：

$$\开始{对齐}X&=T W^{\top}+T_{Y}\bot}P_{Y\perp}^{top}+E，\\Y&=U C^{\top}+U_{X}\bot}P_{X\perp{^{top{+F，\end对齐}$$

(2)

哪里单位=T型注意，在O2PLS配方中，B类_T型=我_第页和H（H）=0.在第一个场景中，我们取N个=500,第页=q个=100,n个=2,n个_X（X）=3和n个_Y（Y）=1.在第二个场景中，我们考虑第页=q个=10⁴.的元素W公司,C类,P（P）_Y（Y）⊥和P（P）_X（X）⊥独立于标准正态分布绘制。JIVE约束通过使关节和特定零件中的每个柱相互正交来应用。的元素T型,T型_Y（Y）⊥和单位_X（X）⊥独立于标准正态分布绘制。噪音，表示为E类和如果，已添加到X（X）和Y（Y）约占总变化的10%。对于r.jive和OmicsPLS，都提取了加载矩阵。为了评估估计精度，计算相应列之间内积的绝对值。这里，较高的值表示较低的估计误差。对于每个场景，我们都生成了1000个副本。

为了研究r.jive在数据分析中缺乏收敛性，我们考虑了两个额外的场景。在第一个附加场景中单位标准偏差为10，即。，单位=10T型在第二个附加场景中，特定部分中的元素将以10的标准偏差进行正态分布。尺寸和样品尺寸如上所示。注意，这两种情况都代表了每个零件变化量的“不平衡”。这里，如果r.jive需要少于500步，则认为它是收敛的。在这些额外的场景中，我们生成了100个副本。

在表中2中位内积值以及中位绝对偏差（MAD）第页=q个=100。可以看出，对于平衡场景设置，OmicsPLS在中值内积方面的表现与r.jive一样好。的结果第页=q个=10⁴与这些结果非常相似（未显示）。

表2 OmicsPLS和r.jive的模拟结果：内积

在表中三，显示了经过的时间和收敛率。OmicsPLS在第一种情况下运行速度快约3500倍(第页=q个=100），第二秒速度快7倍(第页=q个=10⁴)场景。在关节和特定部件之间的变化量不平衡的两种附加场景中，r.jive在大多数运行中都没有收敛。万一单位=10T型，r.jive在90%以上的跑步中没有收敛。如果特定部分包含更多变化，则r.jive在100次运行中的74次和63次未收敛第页=q个=100和第页=q个=10⁴分别是。

表3 OmicsPLS和r.jive w.r.t.median（MAD）在1000次重复中总运行时间（秒）的性能比较，以及在100次运行中的收敛性

数据分析和模拟研究的R代码可作为在线补充（附加文件2).

讨论

在本文中，我们介绍了用于集成两个（组学）数据集的OmicsPLS包。我们通过模拟研究评估了其性能，并使用遗传和IgG1糖组学数据进行了验证。在数据分析方面，遗传PC数据集中联合变异的比例为0.8%。这一比例预计很小，因为不太可能有很大一部分遗传变异（尤其是SNP）与IgG1糖基化有关。在关节成分中，发现了几个可能在IgG1糖基化的遗传调控中起作用的基因。已知这些基因中的一些直接相关（例如。，FUT8（FUT8）)，而其他(DNAJC10公司和AKAP9公司)局限于大多数糖基化发生的细胞隔室（内质网和高尔基体）。然而，关于（IgG）糖基化的遗传调控仍有很多未知之处。

此外，我们考虑对这种类型的数据使用JIVE，但没有成功：算法没有收敛。这种收敛不足的一个潜在原因是两个数据集的不同数据特定特性。特别是，维度和变化量不同。因此，JIVE假设单位=T型可能不合理。我们的仿真证实了这一点：r.jive算法对两个关节部分之间或关节与特定部分之间的变化量的“不平衡”没有鲁棒性。特别是，当单位=10T型在90%以上的重复中，r.jive没有收敛。这表明r.jive可能不适合分析异构数据集（其中数据特定的特征在数据集之间不同）。注意，在DISCO-SCA中，相同的假设(单位=T型)因此，在分析异构数据时，我们预计此方法的性能也会次优。

作为OmicsPLS软件包未来更新的一部分，我们打算处理丢失的数据。为了插补缺失值并同时估计O2PLS成分，OmicsPLS算法可以扩展[2]. 插补步骤也可以在分析之前执行。对于多组学数据，集成回归插补[18]和多因素分析插补[19]已被提议。请注意，与所有插补方法一样，应评估缺失数据造成的不确定性，并将其提交给用户。O2PLS的概率框架将有助于插补，同时解决因数据缺失而产生的额外不确定性。

OmicsPLS的一个重要扩展涉及获得估计值的标准误差。为此，可以应用类似于PLS文献中发现的引导方法[20]. 使用重采样方法的一个缺点是计算负担，特别是对于高维数据集。为了避免此类程序，可以使用O2PLS的概率框架直接计算渐近标准误差。

通过扩展算法以生成稀疏估计，可以提高OmicsPLS输出的可解释性。此扩展可以通过考虑稀疏PLS来实现[21]或者通过考虑O2PLS的概率框架并获得惩罚最大似然估计。

我们目前正在研究概率O2PLS用于数据整合的可能性，这有助于多重插补和统计推断，例如计算渐近标准误差。通过惩罚可能性，可以获得稀疏估计。

由于OmicsPLS是开源的，所以很容易扩展当前的实现来处理更复杂的情况。例如，在GitHub存储库中，可以初始化几个“分支”，在其中可以开发新功能。

结论

我们提出OmicsPLS，这是一个开源且免费的R包，用于将异构数据与O2PLS进行稳健集成。它包括确定零部件数量、装配和检查结果的功能。对于高维数据，使用了节省内存的实现。

可用性和要求

• 项目名称：OmicsPLS

• 项目主页：https://github.com/selbouhaddani/OmicsPLS网站

• 操作系统：Linux、Mac OS、Windows

• 编程语言：R

• 许可证：GPL-3

• 非学者使用的任何限制：无。