CCLasso: correlation inference for compositional data through Lasso

doi:10.1093/bioinformatics/btv349

Warning: The NCBI web site requires JavaScript to function. more...

登录NCBI

生物信息学. 2015年10月1日；31（19）：3172–3180。

2015年6月4日在线发布。内政部：10.1093/生物信息学/btv349

PMCID:PMC4693003型

PMID编号：26048598

CCLasso：通过套索对成分数据进行相关推断

华英坊,^1个，^2个，^三黄成成,⁴ 赵洪宇,⁵和邓明华^1个，^三，^6个，^*

作者信息文章注释版权和许可信息免责声明

这篇文章引用人PMC中的其他文章。

关联数据

补充资料: 补充资料

补充31_19_3172_uindex.html（863字节）
GUID:16EDF234-8E9B-4665-AAC1-5C29E0015F46

supp_btv349_supp_cclasso.pdf（458千）
GUID:582F6B88-192F-4E47-8D3C-FAC864480641

摘要

动机：随着16srrna基因图谱高通量测序技术的发展，直接分析环境和人体内微生物群落变得更加方便和可靠。推断微生物群落成员之间的相关关系是进行基因组调查研究的基础。传统的Pearson相关分析将观测数据视为微生物的绝对丰度，可能导致错误的结果，因为数据只代表相对丰度。在对这些成分数据进行相关分析之前，需要特别注意和适当的方法。

结果：在这篇文章中，我们首先讨论了组成数据潜在变量的相关定义。然后我们提出了一种新的基于最小二乘法的CCLasso方法ℓ₁惩罚从元基因组数据推断组分数据潜在变量的相关网络。利用增广拉格朗日方法中的一种有效的交替方向算法来求解优化问题。仿真结果表明，CCLasso算法在成分数据的边缘恢复方面优于现有的SparCC方法。该方法与SparCC在估算人类微生物群落项目中微生物种类的相关网络方面也有较好的比较。

可用性和实施：CCLasso是开源的，可以从https://github.com/huayingfang/CCLasso在GNU LGPL v3下。

联系人： nc.ude.ukp@hmgned

补充信息： 补充资料可在生物信息学在线。

1简介

微生物在环境和人类生活中扮演着重要的角色。细菌和古生菌是在极端条件下被发现的，比如深海高温喷口和地表下的岩石钻孔(皮库塔et公司 艾尔。2007年). 微生物影响它们存在的环境，反之亦然。据估计，居住在人体内的微生物细胞大约是人类细胞的10倍(萨维奇，1977年). 微生物影响着人类的饮食、健康和医学(腮et公司 艾尔。2006年). 微生物影响人类健康的方式在很大程度上仍然是未知的。对人类微生物群的分析可能有助于我们更好地了解我们自己的基因组。

测序技术不断提高的质量和降低的成本为通过测序分析微生物群落提供了很大的机会。与传统的微生物研究相比，这是一个很大的进步，传统的微生物研究受到几个限制因素的阻碍。首先，只有一小部分微生物可以在实验室条件下培养。第二，实验室只能研究单个微生物，但众所周知，大多数微生物需要其他微生物才能生存。相比之下，测序技术可以让研究人员直接从自然环境中收集群落中所有微生物的全基因组信息，从而促进混合基因组调查(汉德斯曼et公司 艾尔。1998年).

当许多社区都有数据时，微生物之间的依赖关系（可以用相关性来衡量）可能为微生物之间的相互作用提供重要线索。然而，基于测序的调查数据的一个独特之处在于，它们只提供了群落中不同微生物的相对丰度，因为测序结果是测序深度和生物样本量的函数(镍et公司 艾尔。2013年). 因此，从混合基因组调查研究中收集的元基因组数据属于统计学中所谓的组分数据。有人指出皮尔逊（1897）一个多世纪前，为绝对值设计的相关分析方法可能导致成分数据的虚假相关性。对成分数据进行适当的分析和解释需要高度重视和专门的方法。菲尔兹莫瑟和赫隆（2009）提出了一种基于天平的组分数据相关性度量方法，但天平所定义的组分并不总是明确定义和分离的。浮士德et公司 艾尔。（2012年）提出了基于置换和bootstrap的CCREPE来推断相关显著性，但很难解释置换和bootstrap样本之间的差异。弗里德曼和阿尔姆（2012）引入了基于组合数据对数比变换的潜在变量相关概念，提出了稀疏假设下关联矩阵的近似推导方法SparCC。但SparCC没有考虑合成数据中误差的影响，这可能会降低估计精度。此外，不能保证从SparCC推断出的协方差矩阵是正定的，甚至相关系数也可能超出[-1，1]。

在本文中，我们提出了一种新的基于最小二乘法的方法ℓ₁对原始成分数据进行对数比转换后的惩罚，通过潜在变量模型推断微生物之间的相关性，称为套索成分数据相关性推断（CCLasso）。与SparCC类似，CCLasso在相关分析中明确考虑了宏基因组数据的组成性质，并且其额外的好处是，组成数据的潜在变量的估计相关矩阵是正定的。我们还提出了一种有效的交替方向的增广拉格朗日方法来解决该方法所涉及的优化问题。通过交叉验证选择了平衡损失函数和稀疏假设的调谐参数。

通过仿真研究，比较了CCLasso和SparCC的性能，采用了几种相关网络结构和样本大小。仿真结果表明，CCLasso算法比SparCC算法能更准确地估计相关矩阵，并具有更好的边缘恢复能力。将CCLasso和SparCC应用于人类微生物组工程（HMP）中微生物的相关网络估计时，我们发现CCLasso与SparCC具有一致的准确性和可重复性。但是对于经过混洗的HMP数据集，假设任何物种都没有相关性，SparCC总是会产生一些小的相关性，而CCLasso则会将这些小值缩小为0。我们认为CCLasso可以应用于研究自然环境和人体中由宏基因组数据产生的组分数据之间的相关性，它也可以广泛地应用于许多其他有兴趣评估成分数据变量相关性的环境中。

2种方法

2.1成分数据潜在变量的相关性

假设有p微生物种类及其绝对丰度是随机向量是的= (是的₁, …, 是的_p)这在实践中是无法直接观察到的。相反，只有合成随机向量十= (十₁, …, 十_p),

十_{我} = \frac{{是的}_{我}}{\sum_{k = 1}^{p} {是的}_{k}},

（一）

可以从生物学实验中观察到。绝对丰度是的被称为潜在变量，因为它们不能被直接观察到。加性对数正态分布(艾奇森和沈，1980年)是一个特殊情况方程（1）什么时候是的来自多元对数正态分布。三者之间的关系是的比十在实践和理论上都是。微生物之间的相互作用描述如下是的而组成向量呈负相关趋势十即使在没有任何相关性的情况下是的,

\sum_{k = 1}^{p} 十_{k} = 1 \Rightarrow \sum_{k \neq 我} Cov公司 (x_{我}, 十_{k}) = - 变量 (十_{我}) .

让 $w = \sum_{k = 1}^{p} {是的}_{k}$ 是微生物种类的绝对总丰度。绝对协方差是的，以及它的组成代表十，它们可以通过基部联系起来方程（1）,

\begin{matrix} Cov公司 (自然对数 十_{我}, 自然对数 十_{j}) = Cov公司 (自然对数 {是的}_{我}, 自然对数 {是的}_{j}) - Cov公司 (自然对数 {是的}_{我}, 自然对数 w) \\ - Cov公司 (自然对数 w, 自然对数 {是的}_{j}) + 变量 (自然对数 w), \end{matrix}

从ln开始十_我=ln⁡是的_我−ln⁡w. 让Σ_ln⁡十=变量（ln⁡十),Σ_ln⁡是的=变量（ln⁡是的)以及一英寸（Cov）是的，ln⁡w)−Var（ln⁡w)1_p/2其中1_p是一个p×1向量的1′s，则矩阵形式之间的联系Σ_ln⁡十和Σ_ln⁡是的可以这样描述

Σ_ln⁡十= Σ_ln⁡是的 − 一1^T−1个一^T.

（二）

我们可以把重点放在是的我们也叫ln⁡是的潜在变量。什么时候？十由加性logistic正态分布得到ln⁡之间的独立性是的相当于是的。因为那里有信息丢失是的到x通过规范化程序(方程式1)，估计问题Σ_ln⁡是的从样本估计Σ_ln⁡十在没有任何假设的情况下是未定义的。这很容易从中看出方程式（2）在那里p(p+1）/2等式但是p(p+1）/2+p未知参数。

解决这个问题的一个方法是假设Σ_ln⁡是的是稀疏的，这意味着与完全连通的网络相比，微生物物种间的相互作用网络在所有可能的边缘中所占的比例很小。稀疏结构是线性回归模型等欠确定问题的一种常见假设(蒂比拉尼，1996年)，高斯图形模型(2007年，林和袁)以及压缩感知(2005年)其中未知参数的数量比数据点的数量更大，有时甚至更大。对于合成数据，可能存在多个对应的稀疏网络Σ_ln⁡十因为是的以及它的缩放形式C(是的)是的(C(是的)是任意正随机变量，它是一个比例因子）不能与基数区分方程式（1）如果两者都有Σ_ln⁡是的和Σ_{ln⁡(C(是的)是的)}是稀疏的。稀疏程度Σ_ln⁡是的是因为最多有一个稀疏网络Σ_ln⁡是的边缘密度不大于 $\frac{1}{2} - \frac{1}{p - 1}$ 对应于相同的Σ_ln⁡十. 而且这种稀疏密度条件不能放松（参见补充材料). 研究潜在变量ln⁡之间的相关性的统计方法很少是的除了最近引入的一些方法，例如SparCC(弗里德曼和阿尔姆，2012年).

删除一在方程式（2），我们可以选择秩（F）=p−1和F1_p=0且两边有多个F方程式（2）,

F Σ_{自然对数 十} F^{T} = F Σ_{自然对数 是的} F^{T} - F 一 1_{p}^{T} F^{T} - F 1 一_{p}^{T} F^{T} = F Σ_{自然对数 是的} F^{T} .

（三）

左手方程（3）是Fln的方差⁡十与右对应的Fln⁡是的. 他们的关系可以看作是

法兰十平方英尺是的−1个_pln⁡w)=Fln⁡是的.

上述关系可以解释F的两个约束条件。秩（F）=p−1确保十和Fln⁡十因为有常数和约束十。因此，从Fln⁡的统计推断中不存在信息损失十而不是十。一层楼_p=0有助于取消公分母w日志转换后。有许多这样的变换矩阵满足这两个约束条件，例如F=（e_p−1个，−1个_p−1个)线性变换的参考比是对数十_p和 $F = E_{p} - 1_{p} 1_{p}^{T} / p$ 对于组分数据的中心对数比，其中E_p是一个p × p单位矩阵(Chiaitson，1982年).

让Σ= Σ_ln⁡是的= [σ_我j]_p×p. 示例版本S属于Σ_ln⁡十可以通过贝叶斯伪计数方法从原始数据如宏基因组数据中进行分数估计后得到(阿格雷斯蒂和希区柯克，2005年). 从方程（3）以及样本估计S对于Σ_ln⁡十，我们可以得到下面的估计方程，

FΣF^T=FSF^T.

（四）

因为等级（F）=p−1和∑是p × p正定矩阵∑不能通过方程（3）. ∑的额外稀疏假设在许多应用环境中是合理的，例如在元基因组数据中，因为当组分数量很大时，大多数变量对不需要相关。因此，我们可以在没有其他先验信息的情况下，施加一些稀疏约束来帮助建模和推断∑。

2.2 SparCC及其局限性

弗里德曼和阿尔姆（2012）提出了一种迭代逼近方法SparCC来求解估计方程（4）对于变换矩阵F的几种特殊形式，SparCC首先得到了潜在变量ln⁡方差的粗略估计是的_我以及相应的相关矩阵。然后用一个阈值去除最相关的对，重复估计方差和相关性，直到满足某些终止条件。

根据上述符号，SparCC的算法可以概括如下。首先，SparCC从一个粗略的近似中得到∑对角线的估计

\sum_{j \neq 我} σ_{我 j} = 0, \forall 我 .

（五）

粗略近似(式5)附加用品p方程方程（4）. 这一假设意味着，平均而言，每个成分与其他成分没有相关性。让F₁=（−1_p−1个，E_p−1个)以及Σ¹²= (Σ²¹)^T=Cov（ln⁡是的₁，ln⁡是的_−1个), Σ²²=变量（ln⁡是的_−1个)其中ln⁡是的_−1个=（ln⁡是的₂，…，ln⁡是的_p)^T，那么方程（4）可以这样写，

\begin{matrix} (- 1_{p - 1}, E_{p - 1}) Σ {(- 1_{p - 1}, E_{p - 1})}^{T} = F_{1} S F_{1}^{T} \\ \Rightarrow 1_{p - 1} σ_{11} 1_{p - 1}^{T} - 1_{p - 1} Σ^{21} - Σ^{12} 1_{p - 1}^{T} + Σ^{22} = F_{1} S F_{1}^{T} . \end{matrix}

计算上述方程两边的轨迹，我们得到

(p - 1) σ_{11} + \sum_{我 = 2}^{p} σ_{我 我} - 2 \sum_{我 = 2}^{p} σ_{1 我} = tr公司 (F_{1} S F_{1}^{T}) .

如果 $\sum_{我 = 2}^{p} σ_{1 我} = 0$ ，那么 $(p - 1) σ_{11} + \sum_{我 = 2}^{p} σ_{我我} = tr公司 (F_{1} S F_{1}^{T}) .$ 让F_我, 我=2p是加法对数比变换矩阵，其中十_我是参考变量，例如F_p=（E，−1_p−1个). 然后类似于F₁，我们有 $(p - 1) σ_{我我} + \sum_{j \neq 我} σ_{j j} = tr公司 (F_{我} S F_{我}^{T})$ 对于我=2，…=2p从假设中(式5). 相应的解决方案是

σ_{我 我} = \frac{1}{p - 2} (tr公司 (F_{我} S F_{我}^{T}) - \frac{1}{2 (p - 1)} \sum_{我 = 1}^{p} tr公司 (F_{我} S F_{我}^{T})),

（六）

自从 ${(E_{p} + \frac{1}{c} 1_{p} 1_{p}^{T})}^{- 1} = (E_{p} - \frac{1}{c + p} 1_{p} 1_{p}^{T})$ . 然后代入得到基本系数方程式（6）进入之内方程（4）. 其实，以上程序只是一种解决的方法方程（4）和方程式（5）. 这在中称为基本SparCC弗里德曼和阿尔姆（2012）. SparCC的一个潜在问题是σ_二在里面方程式（6）可能是负数，所以是最小值五_最小需要更换负极σ_二.第二，SparCC采用迭代求精方案，在对应的幅值超过给定阈值时排除最强相关对α.这个σ_二根据另一个假设，通过删除最重要的相关对来更新方程式（5）,

\sum_{j \notin C_{我}} σ_{我 j} = 0,

（七）

哪里C_我表示ln⁡的一组指数是的_j与ln强相关是的_我. 最后，SparCC重复前两个步骤来更新方差σ_二以及ln⁡的相关矩阵是的穿过门槛α在给定的迭代时间内，或者直到没有新的强相关对被识别或只剩下三个分量为止。SparCC选择一个相关阈值来给出一个交互网络。

据我们所知，SparCC是推断潜在变量之间相关性的第一种方法是的对于成分数据。第二步是消除强假设的有效方法(式5)在第一近似步骤中方程式（7）在移除最强的对之后。尽管SparCC在分析成分数据方面取得了重大进展，但它在近似计算方面有一些局限性。首先，SparCC直接解决方程（4）在一系列近似假设下方程（4）受这些近似值产生的误差的影响。第二，没有考虑估计的相关矩阵的整体性质。SparCC不能保证推断的相关矩阵是正定的，甚至估计的相关性也可能超出[-1，1]。

2.3套索

我们注意到第一个版本的向量化方程（3）加上样本方差S是ε=（F⊗F）向量(Σ − S)，其中ε满足E(ε)Var=0和0(ε)=（F⊗F）Var（向量(S))（六）^T⊗F^T). 让五_S=Var（向量(S))，则逆方差加权损失函数可给出如下：，

\begin{matrix} {损失}_{1} (Σ) = \frac{1}{2} {(vec公司 (Σ - S))}^{T} (F^{T} \otimes F^{T}) {((F \otimes F) 五_{S} (F^{T} \otimes F^{T}))}^{- 1} \\ (F \otimes F) (vec公司 (Σ - S)), \end{matrix}

（八）

这里的逆符号米^−1个Moore Penrose的伪逆米(彭罗斯和托德，1955年). 减少损失的解决方案₁(Σ)在方程式（8）满足估计值方程（4）. 损失函数的一个重要性质(式8)它对线性变换矩阵F的任何选择都是不变的。这个性质是由于原始数据中∑的信息在对数比变换后保持不变。

功能丧失₁(Σ)在方程式（8）太复杂，无法处理高维协方差矩阵五_S(p² × p²). 灵感来自张和邹（2012）对于样本方差的方差近似，我们可以用下面的损失函数来代替方程式（8）,

损失₁′(Σ)=½tr（（F(Σ − S)F^T)（六）SF^T)^−1个（六）(Σ − S)F^T)).

（九）

变换矩阵F方程式（9）应该合理选择。考虑到组件的对称性，让 $F_{0} = E_{p} - \frac{1}{p} 1_{p} 1_{p}^{T}$ 具有对称投影性质的中心对数比变换矩阵 $F_{0}^{2} = F_{0}, F_{0}^{T} = F_{0}$ . 研究表明，在某些高维问题中，将损失函数中的加权协方差矩阵视为对角矩阵具有很好的效果(陈et公司 艾尔。2013年). 让 $五 = {(诊断 (F_{0} S F_{0}^{T}))}^{- 1}$ . 我们可以考虑另一种替代损失函数，

\begin{matrix} 损失 (Σ) = \frac{1}{2} tr公司 ((F_{0} (Σ - S) F_{0}^{T}) 五 (F_{0} (Σ - S) F_{0}^{T})) \\ = \frac{1}{2} | | F_{0} (Σ - S) F_{0} | |_{五}^{2} . \end{matrix}

（十）

对角矩阵五可以看作F的标准化矩阵₀(Σ − S)F₀. 该方法的核心思想是利用损失函数(式10)因为它简单。

合并∑稀疏假设的合理方法是最小化损失函数和适当的惩罚。理想的惩罚函数是Σ⁻它是∑的反对角线。但是当优化涉及到||Σ⁻||₀是一个指数复杂度的组合优化问题。一种常用的方法是替换ℓ₀-规范ℓ₁-标准(蒂比拉尼，1996年;袁和林，2007). 我们考虑以下目标函数结合损失函数和ℓ₁处罚，

f (Σ) = 损失 (Σ) + 笔 (Σ) = \frac{1}{2} | | F_{0} (Σ - S) F_{0} | |_{五}^{2} + λ_{n} | | Σ^{-} | |_{1},

（十一）

潘在哪里(Σ) = λ_n||Σ⁻||₁. 参数整定λ_n≥0英寸方程式（11）用于平衡模型的拟合度(三)以及∑的稀疏性假设。CCLasso的目标是找到一个正定矩阵 $\hat{Σ}$ 以便

\hat{Σ} = \underset{Σ ≻ 0}{精氨酸 最小} f (Σ) = \underset{Σ ≻ 0}{精氨酸 最小} \frac{1}{2} | | F_{0} (Σ - S) F_{0} | |_{五}^{2} + λ_{n} | | Σ^{-} | |_{1},

（十二）

哪里Σ≻0表示∑应为正定。相应的相关矩阵估计可以通过标准化 $\hat{Σ}$ . 优化所涉及的问题方程式（12）因为两个目标函数都是凸的f(Σ)以及约束区域{Σ|Σ≻0}是凸的。所以局部最小化方程式（12）是全球性的。

与SparCC相比，CCLasso明确考虑了估计背后的误差项方程（4）通过损失函数(式10). 稀疏假设直接通过一个附加的ℓ₁-与附加假设相反的类型惩罚函数(式7)对于SparCC。估计的相关矩阵方程式（12）是正定的，其元素位于正定限制的[-1，1]中。

2.4优化算法及选择λ_n

针对CCLasso中的约束优化问题，提出了一种基于交替方向法的高效算法(张、邹，2012). 轻松的版本方程式（12）可以在去掉正定约束后得到，

\tilde{Σ} = \underset{Σ = Σ^{T}}{精氨酸 最小} \frac{1}{2} | | F_{0} (Σ - S) F_{0} | |_{五}^{2} + λ_{n} | | Σ^{-} | |_{1} .

（十三）

如果解决方案 $\tilde{Σ}$ 在里面式（13）是正定的， $\hat{Σ} = \tilde{Σ}$ . 否则最近的正定矩阵 $\tilde{Σ}$ 用作 $\hat{Σ}$ .

推导一种交替方向法方程（13），我们引入了一个新的矩阵Σ₁重写方程（13）具体如下：，

(\tilde{Σ}, {\tilde{Σ}}_{1}) = \underset{Σ = Σ^{T}, Σ_{1} = Σ}{精氨酸 最小} \frac{1}{2} | | F_{0} (Σ - S) F_{0} | |_{五}^{2} + λ_{n} | | Σ_{1}^{-} | |_{1} .

我们考虑增广拉格朗日函数

\begin{matrix} 我 (Σ, Σ_{1}, Λ) = \frac{1}{2} | | F_{0} (Σ - S) F_{0} | |_{五}^{2} + λ_{n} | | Σ_{1}^{-} | |_{1} \\ + tr公司 (Λ (Σ - Σ_{1})) + (ρ / 2) | | Σ - Σ_{1} | |_{F}^{2}, \end{matrix}

其中| |··||_F是矩阵Frobenius范数。让 $(Σ^{k}, Σ_{1}^{k}, Λ^{k})$ 逐步解决问题k，我们更新(Σ, Σ₁, Λ)根据

Σ^{k + 1} = \underset{Σ = Σ^{T}}{精氨酸 最小} 我 (Σ, Σ_{1}^{k}, Λ^{k}),

（十四）

Σ_{1}^{k + 1} = \underset{_{Σ_{1}}}{精氨酸 最小} 我 (Σ^{k + 1}, Σ_{1}, Λ^{k}),

（十五）

和 $Λ^{k + 1} = Λ^{k} + ρ (Σ^{k + 1} - Σ_{1}^{k + 1})$ . 让Σ^k+1= S+ Δ^k+1为(14)，我们可以写

Δ^{k + 1} = \underset{Δ = Δ^{T}}{精氨酸 最小} \frac{1}{2} | | F_{0} Δ F_{0} | |_{五}^{2} + (ρ / 2) | | Δ | |_{F}^{2} + tr公司 (Δ (Λ^{k} + ρ (S - Σ_{1}^{k}))) .

上面右边的目标函数是Δ和的二次函数Δ^k+1是下式的解，

\frac{1}{2} (F_{0} 五 F_{0} Δ F_{0} + F_{0} Δ F_{0} 五 F_{0}) + ρ Δ = - (Λ^{k} + ρ (S - Σ_{1}^{k})) .

让 $F_{0} = 美国 D_{0} {美国}^{T}, {({美国}^{T} 五美国)}_{11} / 2 ρ + E_{p - 1} / 2 = {美国}_{0} D {美国}_{0}^{T}$ （下标₁₁表示删除最后一行和最后一列）作为相应矩阵的特征值分解 $米 = - {美国}^{T} (Λ^{k} + ρ (S - Σ_{1}^{k})) 美国 / ρ$ ，则上述方程的解为

Δ^{k + 1} = 美国 (\begin{matrix} {美国}_{0} {({美国}_{0}^{T} 米_{11} {美国}_{0}) \cdot C} {美国}_{0}^{T} & 米_{12} \\ 米_{21} & 米_{22} \end{matrix}) {美国}^{T},

（十六）

式中，°是矩阵和 $C_{我 j} = \frac{1}{D_{我我} + D_{j j}}$ .

更新 $Σ_{1}^{k + 1}$ ，我们定义了一个运算符G(A, λ)具体如下：，

G {(A, λ)}_{我 j} = {\begin{matrix} A_{我 j} & 我 = j, \\ A_{我 j} - λ & 我 \neq j, A_{我 j} > λ, \\ A_{我 j} + λ & 我 \neq j, A_{我 j} < - λ, \\ 0, & 我 \neq j, - λ \leq A_{我 j} \leq λ . \end{matrix}

从方程式（15），我们写道

Σ_{1}^{k + 1} = \underset{Σ_{1}}{精氨酸 最小} (ρ / 2) | | Σ_{1} | |_{F}^{2} - tr公司 (Σ_{1} (Λ^{k} + ρ Σ^{k + 1})),

那么上述问题的解决方案是 $Σ_{1}^{k + 1} = G (\frac{Λ^{k}}{ρ} + Σ^{k + 1}, \frac{λ_{n}}{ρ}) .$ 下面的算法详细总结了上述交替方向法求解优化问题的实现方法(式12)为了卡索。

初始化：k = 0，Λ⁰, $Σ_{1}^{0} = E_{p}$ .
重复（a）-（d）直到Σ^k和 $Σ_{1}^{k}$ 聚合：
1. Σ^k+1 ← S+ Δ^k+1哪里Δ^k+1是屈服的方程式（16）;
2. $Σ_{1}^{k + 1} \leftarrow G (\frac{Λ^{k}}{ρ} + Σ^{k + 1}, \frac{λ_{n}}{ρ})$ ;
3. $Λ^{k + 1} \leftarrow Λ^{k} + ρ (Σ^{k + 1} - Σ_{1}^{k + 1})$ ;
4. k ← k+1。
返回聚合的Σ^k作为解决方案 $\tilde{Σ}$ 定义于方程（13）.

调谐参数λ_n≥0英寸方程式（11）必须进行调整，因为它控制了模型适应度之间的平衡(方程式3)稀疏性假设。AK-损失函数的一般交叉验证(式10)用于选择λ_n在本文中。首先将所有样品分为K不相交的子群，如褶皱我_k对于k=1K. 这些折叠和测试将轮流使用。第二，每个人k=1K，计算S_k和S_−k对应于Var（ln⁡）的样本估计p)通过我_k和我₁, …, 我_k−1个, 我_k+1, …, 我_K. 下标-k意味着使用所有样本k-向左折叠。权重矩阵五因为训练数据和测试都是基于所有的数据。第三，让S= S_−k然后计算估计值 ${\hat{Σ}}_{- k}$ 通过方程式（12）每1≤k ≤ K. 然后计算K-调整参数的折叠交叉验证错误λ_n,

个人简历 (λ_{n}) = \frac{1}{K} \sum_{k = 1}^{K} \frac{1}{2} | | F_{0} ({\hat{Σ}}_{- k} - S_{k}) F_{0} | |_{五}^{2} .

最后，我们选择 $λ_{n}^{*} = 精氨酸 {最小}_{λ_{n}} 个人简历 (λ_{n})$ 作为最终的调谐参数。

3结果

3.1模拟研究

虽然基因组调查研究的目标是从丰度计数矩阵推断微生物群落成员之间的相关性，但可以比较使用CCLasso或SparCC的相关矩阵的估计精度，以评估它们的相对性能，因为它们都基于相同的潜在假设，如方程（1）. 这两种方法的本质区别在于得到分数估计后的估计过程。

成分数据由加性logistic正态分布模拟，具有给定的均值和协方差矩阵，

自然对数 是的 \sim N (μ, Σ), 十_{我} = \frac{{是的}_{我}}{\sum_{k = 1}^{p} {是的}_{k}} .

变化参数μ控制部件的不平衡。每一个元素μ由[−0.5，0.5]的均匀分布生成。在我们的仿真中，我们重点比较了SparCC和CCLasso在不同稀疏度的稀疏相关矩阵上的性能。考虑五种协方差结构：

随机模型：每对分量以给定概率0.3连接，相关强度为±0.15，等概率为0.5。
邻居模型：随机选择p[0,1]中的点²飞机。然后用相关强度0.5连接每个点的10个最近邻。
AR（4）型号：连接对(我,j)如果|我 − j|≤4，当距离分别为1、2、3或4时，将相关强度设为0.4、0.2、0.2或0.1。
枢纽模型：随机选取3个点作为枢纽，另一个作为枢纽p–3点作为共同点。然后以0.7的概率将每个轮毂连接到其他轮毂上，同时在公共点之间创建概率为0.2的边，所有边强度设置为0.2。
块模型：分割p平均分为5个区块。以概率0.6和相关强度0.4连接同一块中的每一对，以概率0.2和相关强度0.2连接不同块中的点。

为了使协方差矩阵正定，将∑的对角元素设置得足够大，然后将其归一化为1。随机模型是一种非常常见的图模型，其中每个可能的边都以相同的概率独立出现。通过等概率设定强度为±0.15，随机模型基本满足假设方程式（5）和（七）斯巴克的。邻域模型是一种二维地理模型，其中最近邻之间存在边。AR（4）模型可以看作是一个点沿着距离不大于4的节点之间存在边的直线进行线性排序的模型，并且相关性随着距离的增加而减小。集线器模型描述了一个图，其中一些特殊的节点（称为集线器）以比其他节点之间的连接概率更高的概率连接到其他节点。块模型定义了组内边缘概率高于组间边缘概率的网络聚类。所有模型都具有不同程度的稀疏性。邻域模型和AR（4）模型中的期望边数与p虽然p²对于随机、集线器和块模型。

对于所有型号，我们p = 50，并考虑不同的样本量n = 200、300和500。每种型号，以及(p,n)，我们重复模拟100次。调谐参数λ_n是通过3次交叉验证确定的，所有数据都用来估计∑和相关矩阵。我们使用R和默认的调优参数重新实现了SparCCα=0.1，k_最大值 = 10和五_米我n=10^−4个在最终相关性被-1和1截断作为下限和上限时使用。SparCC对于它的调优参数是健壮的，因为在每次迭代中只删除最强的一对(补充图S1).

为了比较CCLasso和SparCC在模型设置和样本量的每个组合下的性能，我们用平均绝对误差定义了相关推断的精度 $d_{1} (\hat{ρ}, ρ) = \frac{2}{p (p - 1)} \sum_{我 < j} | {\hat{ρ}}_{我 j} - ρ_{我 j} |$ 以及Frobenius范数距离 $d_{F} (\hat{ρ}, ρ) = | | \hat{ρ} - ρ | |_{F}$ 估计的相关矩阵之间 $\hat{ρ}$ 而真正的那个呢ρ.利用接收机工作特性曲线（AUC）下的面积来评估CCLasso和SparCC在恢复稀疏协方差矩阵∑中非零项的性能，以避免阈值参数的选择。

表1总结了CCLasso和SparCC在仿真研究中的性能d₁和d_F距离和AUC。当样本量从200增加到500时d₁和d_F在每个模拟设置中减少CCLasso和SparCC。CCLasso的估计误差小于SparCC。相应的结果表明CCLasso比SparCC具有更好的模拟性能。这可能是因为CCLasso考虑了随机错误，而SparCC没有。在边缘恢复方面，除了随机图模型在样本大小为200和300时，CCLasso的性能也优于SparCC。这可以解释为随机模型大致满足假设方程式（5）和（七）斯巴克的。准确度和AUC不一致，例如d₁和d_F在随机模型中，CCLasso比SparCC小，但SparCC的AUC大于CCLasso。造成这种现象的原因是，准确度衡量的是估计值与真实值之间的连续距离，而AUC则显示了非零和零之间的区别。

表1。

基于仿真结果的CCLasso和SparCC性能比较

n	方法	d₁	d_F	AUC
随机模型
200	套索	0.033（0.001）	2.954（0.049）	0.791（0.015）
	斯巴克	0.057（0.001）	3.528（0.080）	0.823（0.014）
300	套索	0.028（0.001）	2.409（0.057）	0.885（0.012）
	斯巴克	0.047（0.001）	2.901（0.059）	0.891（0.011）
500	套索	0.023（0.001）	1.994（0.053）	0.953（0.007）
	斯巴克	0.038（0.001）	2.332（0.056）	0.951（0.006）
邻域模型
200	套索	0.039（0.003）	3.355（0.206）	0.948（0.015）
	斯巴克	0.076（0.001）	4.606（0.081）	0.888（0.014）
300	套索	0.033（0.002）	2.675（0.151）	0.986（0.006）
	斯巴克	0.070（0.001）	4.176（0.060）	0.931（0.009）
500	套索	0.026（0.002）	2.064（0.121）	999.0.001（999.0.001）
	斯巴克	0.065（0.001）	3.800（0.041）	0.967（0.006）
AR（4）模型
200	套索	0.021（0.001）	2.444（0.134）	0.885（0.021）
	斯巴克	0.061（0.001）	3.766（0.087）	0.858（0.019）
300	套索	0.018（0.001）	1.994（0.133）	0.922（0.017）
	斯巴克	0.052（0.001）	3.210（0.078）	0.890（0.017）
500	套索	0.015（0.001）	1.549（0.087）	0.958（0.011）
	斯巴克	0.044（0.001）	2.693（0.059）	0.918（0.011）
轮毂型号
200	套索	0.037（0.001）	3.453（0.037）	0.749（0.021）
	斯巴克	0.067（0.001）	4.194（0.070）	0.690（0.014）
300	套索	0.036（0.001）	3.133（0.047）	0.768（0.021）
	斯巴克	0.059（0.001）	3.686（0.049）	0.735（0.012）
500	套索	0.032（0.001）	2.918（0.048）	0.828（0.018）
	斯巴克	0.051（0.001）	3.248（0.043）	0.788（0.010）
块体模型
200	套索	0.039（0.001）	3.307（0.113）	0.782（0.014）
	斯巴克	0.070（0.001）	4.268（0.072）	0.734（0.010）
300	套索	0.035（0.001）	2.773（0.079）	0.854（0.014）
	斯巴克	0.062（0.001）	3.788（0.052）	0.765（0.011）
500	套索	0.029（0.001）	2.258（0.076）	0.924（0.011）
	斯巴克	0.057（0.001）	3.374（0.038）	0.796（0.012）

d₁和d_F估计的相关矩阵与文本中定义的真实相关矩阵之间的两个距离。AUC是接收器工作特性曲线下的面积。结果是100次模拟运行的平均值，括号内有标准偏差。

有关ROC的更详细结果如所示图1. 随着样本量的增加，CCLasso和SparCC之间的差距也随之增大。对于0.1等较低的假阳性率，除随机图模型外，CCLasso的真阳性率均大于SparCC。一个有趣的现象是，CCLasso和SparCC对于hub模型的性能都很差，但是随着样本量的增加，估计效率提高。对于某些特殊的图形结构，例如hub模型和块模型，在给定的精度下，应该使用比其他模型大得多的样本量。我们还通过ROC对CCREPE和CCLasso进行了比较，发现CCREPE的性能与SparCC相似(补充图S2).

图1。

CCLasso和SparCC的ROC曲线。确定假阳性率后，平均100次重复的真阳性率，灰色线为基线参考

3.2 HMP数据

由于我们和我们体内的微生物有着密切的关系，人类微生物组项目联盟(2012年a,b)目的探讨微生物在人类健康和疾病中的基本作用。利用HMP健康个体16S可变区3-5（V35）的高质量测序读数，探索18个身体部位微生物之间的相互关系，并从HMPOC数据集中获得相应的操作分类学单位（OTU），可在http://www.hmpdacc.org/HMMCP/. 我们考虑第一阶段生产研究（2010年5月1日）和从同一个体的身体部位采集多个样本的第一个样本。通过移除读取次数少于500次或收集到超过60%0的样本，以及移除平均每个样本读取次数少于2次或超过60%0的OTU，进一步过滤数据。由于CCLasso和SparCC都假设0为OTU计数不是实的0分数。CCLasso将所有计数加上最大舍入误差0.5，然后将计数标准化以获得成分数据。弗里德曼和阿尔姆（2012）为SparCC提供了从计数中估计分数的贝叶斯框架。SparCC的最终估计值是来自分数后验分布的20个重复样本中估计值的中位数。

由于在实际数据中没有分类单元-分类单元相互作用的真实相关网络的先验信息，我们使用一致的准确性和再现性来比较CCLasso和SparCC的性能。首先，所有数据被用来构建CCLasso和SparCC的金标准参考相关矩阵。由于使用了所有数据，因此该步骤中估计的相关矩阵被视为“已知”。其次，我们随机选取半个样本，通过CCLasso和SparCC估计相关矩阵。通过第一步和第二步估计的相关矩阵之间的Frobenius范数距离来衡量一致精度。一致的再现性通过在第一个黄金参考网络（仅使用顶部1/4边缘）中为这两个步骤共享的相同边缘的分数来测量。此程序重复20次，以获得稳定的结果。

结果总结如下表2. CCLasso和SparCC在一致的准确性和再现性方面具有相似的性能。当样本量较小时，重复性较低。即使对于大样本量，如左前肘窝，CCLasso和SparCC的重复性也只有0.64。从模拟数据中我们可以发现一致的准确度和再现性不是很好的标准(补充表S1). 由于CCLasso的交叉验证有多个优化过程，SparCC比CCLasso更快(补充表S2). 再现性对于上边缘的选择是稳健的(补充表S3). 我们还比较了CCLasso和SparCC对所有身体部位的所有样本推断出的相关网络，发现它们的结果非常相似(补充图S3和表S4).

表2。

从HMP数据得出的CCLasso和SparCC在不同身体部位的一致的Frobenius准确性和再现性

身体	样品	弗罗贝尼乌斯准确度		再现性
现场	大小	套索	斯巴克	套索	斯巴克
安特纳	152	2.28（0.17）	第二十二条第二款	0.70（0.05）	0.68（0.05）
阿克林	193	1.71（0.13）	1.56（0.14）	0.75（0.04）	0.77（0.05）
鼻涕	196	2.47（0.17）	2.11（0.11）	0.71（0.03）	0.72（0.03）
哈帕	197	2.57（0.18）	2.24（0.13）	0.79（0.03）	0.80（0.03）
兰特福斯	51	4.35（0.31）	6.57（0.51）	0.64（0.05）	0.64（0.04）
LRetCre公司	123	2.24（0.21）	2.30（0.15）	0.69（0.04）	0.66（0.04）
迷走神经	45	2.69（0.48）	3.35（0.59）	0.64（0.08）	0.64（0.07）
帕顿	203	2.76（0.17）	2.31（0.17）	0.83（0.02）	0.83（0.02）
PosFor公司	22	3.34（0.68）	2.98（0.80）	0.67（0.14）	0.68（0.11）
兰特福斯	54	3.31（0.37）	6.32（0.32）	0.54（0.03）	0.60（0.04）
RRetCre公司	85	2.64（0.17）	3.50（0.20）	0.63（0.04）	0.61（0.05）
唾液	184	2.95（0.14）	2.75（0.13）	0.75（0.03）	0.77（0.03）
大便	190	1.81（0.13）	2.10（0.16）	0.72（0.04）	0.71（0.03）
次级	205	2.89（0.23）	2.45（0.18）	0.82（0.02）	0.84（0.03）
供应	207	2.74（0.22）	2.31（0.12）	0.82（0.02）	0.85（0.02）
喉部	197	2.79（0.15）	2.48（0.14）	0.79（0.03）	0.80（0.02）
通多	207	2.52（0.22）	2.00（0.17）	0.83（0.02）	0.84（0.02）
流浪汉	52	2.93（0.27）	2.94（0.21）	0.65（0.06）	0.63（0.05）

AntNar，前鼻孔；AKerGin，附着角质化牙龈；BucMuc，颊粘膜；HarPal，硬腭；LAntFos，左侧肘前窝；LRetCre，左耳后皱褶；中迷走神经，阴道中部；PalTon，腭扁桃体；PosFor，后穹窿；RAntFos，右侧肘前窝；RRetCre，右侧耳后皱褶；Subla，龈下斑块；供应，龈上菌斑；托多，舌背；迷走神经，阴道入口。

结果是20次复制运行的平均值，括号内有标准偏差。

我们还通过随机HMP数据比较了CCLasso和SparCC的性能。每个otu的个体计数被置换，因此它不应该发现物种之间的任何相关性。图2显示了通过CCLasso和SparCC对无序数据集估计的相关性的柱状图。CCLasso几乎没有检测到任何相关性，但是SparCC中总是有一些小的相关性。在这种情况下，CCLasso的表现优于SparCC。我们将CCLasso和SparCC用于另一个数据集，发现SparCC检测到太多强且无意义的边缘(补充图S4).

图2。

通过CCLasso和SparCC对随机HMP数据集估计相关性的直方图

4讨论

尽管成分数据在许多实际问题中自然出现，但研究人员通常对这些数据背后的潜在变量更感兴趣。例如，在基因组调查研究中，从观察到的细菌相对丰度（而不是绝对丰度）推断不同细菌之间的依赖性是非常有趣的。因此，有必要为成分数据推断潜在变量之间的相关矩阵。在这篇文章中，我们提出了一种新的方法来推断组分数据的潜在变量之间的相关性。通过求解常和约束问题，利用稀疏假设帮助估计潜在变量的相关矩阵。仿真结果表明，CCLasso比SparCC具有更好的性能，SparCC是我们所知道的从潜在变量的角度来解决这个问题的唯一可用方法。对于HMP数据，CCLasso与SparCC具有相似的一致精度和再现性。但是，从混合HMP数据集中，我们发现SparCC总是给出一些非零估计。

尽管CCLasso在仿真研究中的性能优于SparCC，但它也存在着与SparCC相似的困难，比如分量分数估计的可靠性和线性关系的解释。我们采用简单的伪计数0.5来避免HMP数据集的0个组件。还有其他规范化方法可以解释采样不足，例如保尔森et公司 艾尔。（2013年）介绍了一种评估稀疏高通量微生物标记基因调查数据中差异丰度的方法。最近，比斯瓦et公司 艾尔。（2014年）提出了一个Poisson多元正态分层模型来学习直接交互作用，消除了宏基因组测序实验中的混杂预测因子效应。假设比斯瓦et公司 艾尔。（2014年）从潜在模型的角度来看，与CCLasso相似，但两者的本质区别在于成分假设。未来的工作将集中在改进组分估计以解释欠采样和探索微生物之间的非线性关系。

基金

这项工作得到了国家自然科学基金的资助。31171262;,31428012;,31471246)国家重点基础研究项目（第1期）。2015CB910303). H、 F.的工作部分得到了北京大学研究生院的一项研究金的支持。H、 Z.的工作部分得到了NIH的支持GM59507标准.

利益冲突：未声明。

补充材料

补充资料：

单击此处查看。

工具书类

Agresti A.，Hitchcock D.B.（2005年）。分类数据分析的贝叶斯推理.统计方法应用。,14公元297-330年。[谷歌学者]
艾奇森J.（1982年）。成分数据的统计分析.J、 R.统计Soc。B,44第139-177页。[谷歌学者]
艾奇森.J.，沈S.M.（1980）。Logistic正态分布的一些性质和应用.生物计量学,67第261-272页。[谷歌学者]
Biswas S.等人（2014年）。从宏基因组数据学习微生物相互作用网络。arXiv:1412.0207v1[q-bio.QM]。[谷歌学者]
Candes E.J.，Tao T.（2005年）。线性规划译码.IEEE T.通知。理论,51，4203-4215年。[谷歌学者]
陈杰等（2013）。结构约束稀疏典型相关分析及其在微生物数据分析中的应用.生物统计学,14，244-258年。[PMC免费文章][公共医疗][谷歌学者]
Faust K.等人（2012年）。人类微生物群落中的微生物共生关系.公共科学图书馆计算机。生物。,8，e1002606。[PMC免费文章][公共医疗][谷歌学者]
Filzmoser P.，Hron K.（2009年）。成分数据的相关分析.数学。Geosci公司。,41公元905-919年。[谷歌学者]
Friedman J.，Alm E.J.（2012年）。从基因组调查数据推断相关网络.公共科学图书馆计算机。生物。,8，e1002687。[PMC免费文章][公共医疗][谷歌学者]
Gill S.R.等人（2006年）。人远端肠道微生物组的宏基因组分析.科学类,3121355年至1359年。[PMC免费文章][公共医疗][谷歌学者]
Handelsman J.等人（1998年）。未知土壤微生物化学的分子生物学研究：天然产物的新前沿.化学。生物。,5，R245–R249。[公共医疗][谷歌学者]
人类微生物群项目联盟。（2012年a）。人体微生物研究框架.自然,486，215–221页。[PMC免费文章][公共医疗][谷歌学者]
人类微生物群项目联盟。（2012年b）。健康人体微生物群的结构、功能和多样性.自然,486公元207-214年。[PMC免费文章][公共医疗][谷歌学者]
Ni J.等人（2013年）。有多少宏基因组测序足以达到既定目标？.科学。代表。,三，1968年，doi:10.1038/srep1968。[PMC免费文章][公共医疗][谷歌学者]
Paulson J.N.等人（2013年）。微生物标记基因调查的差异丰度分析.纳特。方法,101200年至1202年。[PMC免费文章][公共医疗][谷歌学者]
皮尔逊K.（1897）。关于在器官测量中使用指数时可能出现的一种伪相关形式.程序。R、 Soc。隆德。,60公元489-502年。[谷歌学者]
彭罗斯R.，托德J.A.（1955年）。矩阵的广义逆.数学。程序。剑桥,51.[谷歌学者]
Pikuta E.V.等人（2007年）。生命极限微生物.致命一击。版次。微生物学。,33183-209年。[公共医疗][谷歌学者]
萨维奇特区（1977年）。胃肠道微生物生态学.每年。版次。微生物学。,31第107-133页。[公共医疗][谷歌学者]
Tibshirani R.（1996年）。回归收缩和套索选择.J、 R.统计Soc。B,58公元267-288年。[谷歌学者]
袁明，林勇（2007）。高斯图形模型中的模型选择与估计.生物计量学,9419-35岁。[谷歌学者]
张涛，邹赫（2012）。基于lasso惩罚d-迹损失的稀疏精度矩阵估计.生物计量学,99，第1-18页。[谷歌学者]

文章来自生物信息学在这里提供牛津大学出版社

支援中心支援中心

项目管理咨询公司

CCLasso：通过套索对成分数据进行相关推断

华英坊

黄成成

赵洪宇

邓明华

关联数据

摘要

1简介