LDA Effect Size分析 LEfSe详解

LEfSe的作用
LEfSe的原理

LEfSe的作用

在介绍LEfSe的作用前，我们先解释一个概念——biomarker，维基百科给出的定义是

A bio-marker, or biological marker is a measurable indicator of some biological state or condition. Biomarkers are often measured and evaluated to examine normal biological processes, pathogenic processes, or pharmacologic responses to a therapeutic intervention.

用我们搞数据的人能理解的话讲，biomarker就是非常强力的用来分类的特征，它可以是基因、细胞或者物种分类单元等。比如（瞎编的例子，不能当真）某个研究团队发现脚气会影响中央后回位于Sylvian Fissure附近的区域，从而影响舌头的知觉，于是这个团队打算进一步研究脚气怎么通过肠脑轴影响舌头的知觉，他们随机调查了一批志愿者，记录了志愿者的一些demographic information以及病理信息，并记录了他们肠道菌群的物种分类信息与物种丰度，于是他们有了一张数据表：

Group	有脚气	无脚气
界.门.纲.目.科.属1	丰度	丰度
界.门.纲.目.科.属2	丰度	丰度
。。。

他们想知道哪个属的细菌的丰度在有脚气与无脚气的志愿者之间是存在显著差异的，这个时候就需要LDA Effect Size分析了。

也就是说LDA Effect Size分析的作用是发现不同group之间存在显著差异的biomarker。下面我们介绍LDA Effect Size分析的原理。

LEfSe的原理

首先我们写出数据，用iii表示第iii个志愿者，i=1,2,⋯,ni=1,2,\cdots,ni=1,2,⋯,n，用yiy_iyi表示第iii个志愿者所在的group，yi∈{1,2,⋯,K}y_i \in \{1,2,\cdots,K\}yi∈{1,2,⋯,K} (比如讨论有无脚气时K=2K=2K=2，我们可以用yi=1y_i=1yi=1表示志愿者iii有脚气，用yi=2y_i=2yi=2表示志愿者iii无脚气)，用xix_ixi表示第iii个志愿者的肠道菌群物种分类信息，
xi=(xi1,⋯,xiM)Tx_i = (x_{i1},\cdots,x_{iM})^Txi=(xi1,⋯,xiM)T

比如xi1x_{i1}xi1可以表示是Bacteroidaceae(拟杆菌科)、xi2x_{i2}xi2可以表示Coriobacteriaceae（红蝽杆菌科）；用ziz_izi表示第iii个志愿者的demographic information以及病理信息，
zi=(zi1,⋯,ziL)Tz_i = (z_{i1},\cdots,z_{iL})^Tzi=(zi1,⋯,ziL)T

比如zi1z_{i1}zi1可以表示性别、zi2z_{i2}zi2表示种族、zi3z_{i3}zi3表示BMI等。

第一步 Kruskal-Wallis检验
因为MMM通常是一个非常大的数，也就是说肠道内包含的细菌属非常多，但我们只需要找到效果最好的那几个作为biomarker，所以在使用LDA分析前，我们最好先做一个预处理：选出在两个group间具有显著差异的属，筛掉那些不具有显著差异的属。

用μmj\mu_{mj}μmj表示第jjj个group的第KKK个细菌属的平均丰度，j=1,⋯,Kj=1,\cdots,Kj=1,⋯,K，m=1,⋯,Mm=1,\cdots,Mm=1,⋯,M，我们想要检验：
H0:μm1=μm2=⋯=μmKHa:至少有一个与其他的不等H_0:\mu_{m1} = \mu_{m2} = \cdots = \mu_{mK} \\ H_a:至少有一个与其他的不等H0:μm1=μm2=⋯=μmKHa:至少有一个与其他的不等

在介绍单因素试验设计的时候，我们介绍了检验多总体均值的方法——ANOVA，然而ANOVA需要总体服从正态分布，然而物种丰度通常会服从有偏的分布，所以我们并不能使用ANOVA来做这个检验，不过我们也讨论过，当总体的正态假设不满足时，我们可以使用非参的Kruskal-Wallis检验。

第二步 Wilcoxon秩和检验
接下来我们再做一个数据预处理，我们希望菌群丰度差异是group不同造成的，（比如菌群丰度主要是由有无脚气造成的），那么我就需要证明其他因素不会造成菌群丰度的显著差异。

比如在完成第一步以后，如果我们发现有无脚气的志愿者他们肠道红蝽杆菌科的丰度具有显著差异，但又有研究者提出，红蝽杆菌科丰度的差异可能与性别有关，于是我们按照性别把第j,1≤j≤Kj,1 \le j \le Kj,1≤j≤K个group每一个都分为2个subgroup，我们希望同一个group的两个subgroup之间的红蝽杆菌科丰度没有显著差异，但是不同group的subgroup之间红蝽杆菌科丰度存在显著差异。为了实现这个比较，我们需要给这2K2K2K个subgroup的红蝽杆菌科丰度做两两比较。在单因素试验设计中我们介绍过三种做两两比较的方法：Tukey检验、Fisher Least Significant Difference方法、Dunnett方法。遗憾的是，这三种方法都要求基于ANOVA的假设进行，因此在这个应用场景下，我们这三种方法同样一个也不能用。这时我们同样就需要非参数检验了，Wilcoxon秩和检验就可以在subgroup之间进行两两比较。

对第一步选出来的那些biomarker都使用Wilcoxon秩和检验，同一个group的两个subgroup之间没有显著差异，不同group的subgroup之间存在显著差异的biomarker我们可以保留下来，然后扔掉其他的biomarker。需要注意的是把group分成不同的subgroup是基于ziz_izi进行的，如果存在某个zilz_{il}zil它可能造成biomarker的显著差异，我们就按zilz_{il}zil的不同取值分subgroup。

比如我们按照性别把第j,1≤j≤Kj,1 \le j \le Kj,1≤j≤K个group每一个都分为2个subgroup，但是我们发现某些group的两个subgroup之间的红蝽杆菌科丰度存在显著差异，这就说明group之间红蝽杆菌科丰度的显著差异可能并不是由group(有无脚气)导致的，而是由性别决定的，于是我们排除掉红蝽杆菌科丰度。如果我们发现某些group的两个subgroup之间的红蝽杆菌科丰度没有显著差异，那么我们可以留下红蝽杆菌科丰度，并将其作为指示脚气的一个重要biomarker。

要找存在显著影响的zilz_{il}zil也很简单，对丰度使用adonis，用zi1,⋯,ziLz_{i1},\cdots,z_{iL}zi1,⋯,ziL以及group作为feature，zi1,⋯,ziLz_{i1},\cdots,z_{iL}zi1,⋯,ziL中显著的那些feature就应该用来做一下Wilcoxon秩和检验。

第三步 LDA
现在我们就有了质量还不错的数据了，这些数据就是显著的biomarker，并且它们的显著差异主要是由group不同造成的，而不会受其他confounding variable的显著影响。

那么最后我们要做的就是评估这些biomarker的重要性了，我们可以用LDA进行分析。简单解释一下LDA的原理，之前我没写完的统计学习的博客里推导了LDA的输赢比的公式：

ln⁡P(Y=1∣X=x)P(Y=−1∣X=x)=ln⁡π1exp⁡{−12(x−μ1)TΣ−1(x−μ1)}π0exp⁡{−12(x−μ0)TΣ−1(x−μ0)}=[−12(x−μ1)TΣ−1(x−μ1)+ln⁡π1]−[−12(x−μ0)TΣ−1(x−μ0)+ln⁡π0]\ln{\frac{P(Y=1|X=x)}{P(Y=-1|X=x)}} = \ln{\frac{\pi_1 \exp{\{-\frac{1}{2}(x-\mu_1)^T\Sigma^{-1}(x-\mu_1)\}}} {\pi_0 \exp{\{-\frac{1}{2}(x-\mu_0)^T\Sigma^{-1}(x-\mu_0)\}}}} \\ = \left[-\frac{1}{2}(x-\mu_1)^T\Sigma^{-1}(x-\mu_1) + \ln \pi_1 \right] - \left[-\frac{1}{2}(x-\mu_0)^T\Sigma^{-1}(x-\mu_0) + \ln \pi_0 \right] lnP(Y=−1∣X=x)P(Y=1∣X=x)=lnπ0exp{−21(x−μ0)TΣ−1(x−μ0)}π1exp{−21(x−μ1)TΣ−1(x−μ1)}=[−21(x−μ1)TΣ−1(x−μ1)+lnπ1]−[−21(x−μ0)TΣ−1(x−μ0)+lnπ0]

每个中括号里面的式子具有一样的构造，它代表的是数据点到centroid的距离，这个式子的值越大，那么对应的biomarker的效果就越强。

比如还是在脚气的例子中，假设我们通过第一步和第二步筛选出NNN个biomarker，那么nnn个志愿者就可以表示成nnn个在NNN维空间中的点，点的第aaa个坐标表示第aaa种biomarker的丰度(1≤a≤N1 \le a \le N1≤a≤N)。这nnn个点会集中在两个中心，这两个中心的坐标分别表示有脚气/无脚气的志愿者第aaa种biomarker的平均丰度(1≤a≤N1 \le a \le N1≤a≤N)，我们称这两个中心为centroid。我们可以把这nnn个点投影到第aaa个坐标轴上，就得到了nnn个位于直线上的点，有脚气的会落在一些区域，没有脚气的会落在另一些区域，如果这两个区域几乎不重叠，那么第aaa个biomarker就是一个非常重要的biomarker。我们可以用上面的对数输赢比来度量biomarker的重要性，比如设置对数输赢比大于2表示biomarker非常重要，那么我们就可以用对数输赢比的柱状图或者Cladogram等可视化方法来展示我们发现的biomarker的重要性以及它们的演化规律。