释放收费版本-把握LEfSe核心算法使用R实现整套算法

写在前面

曾几何时，我们再网站上运行LEfSe，面对的是网速限制和分析框架的限制，后来我们使用本地按安装版本解决了网络分析的风险。

现在可能更加方便了使用一条命令，即可实现整个过程：

conda install lefse

其实我们除了分析，很大程度上是因为想使用Graphlan的可视化进行出图。这一点我门不得不佩服LEfSe，使用Graphlan出图捕获了大量的用户，其实这个有多简单，只不过是将门类信息转化为树的结构而已。但再R语言中很长一段时间都没有成熟的工具，直到去年，开发了microbiomeviz，模仿了LEfSe分析的树图，但是模仿Graphlan的功能似乎还不够，当然这不是我们今天的主题，所以先说到这里。虽然树图是特色，但是作为一个分析最核心的分析过程LEfSe使用的确实R语言。也就是说目前的LEfSe离不开R语言，核心的两个分析：非参数检验，使用的coin包，LDA判别使用的MASS包。虽然LEfSe没有展示代码，但这两个包如何做这两种分析确是很容易学会。所以这样的问题就回来了，我们R语言都可以掌握计算的核心了，那为什么还要被python封装起来的工具LEfSe放到底层，而不为人所知呢？你们使用了这么多年的LEfSe了，知道这个核心功能是R语言在默默支持吗？

为R语言正名—实现LEfSe核心算法

下面这张图形想必大家很了解，因为这是LEfSe原理图，目前在网上能看到的解释就是这样：

• A. 首先在多组样本中采用的非参数因子Kruskal-Wallis秩和检验检测不同分组间丰度差异显著的物种；

• B. 再利用Wilcoxon秩和检验检查在显著差异物种类中的所有亚种比较是否都趋同于同一分类级别；

• C. 最后用线性判别分析（LDA）对数据进行降维和评估差异显著的物种的影响力（即LDA score）。

我并不是说这个不对，因为今天我是为R语言正名的，也就是说上面的原理我们根本无法知道R语言在这张原理图的地位。是吧，根本不可能知道的。

其实在一般情况下step2是很少使用的，虽然下面的3个步骤中就是步骤2是最难理解的，但是它却不是我们经常使用的，比如：你想要使用两组或者三组数据做lefse，根本不会用到步骤2。这是为具有subclass的分组准备的，是什么意思呢？比如我们有两个组A，B，各12个样本，要做LEfSe，一个分组的情况下就只需要指定A，B信息作为class就足够了，但是如果A中存在差异，比如A是癌症病人患者，那么这里有男人和女人患者，这个性别信息定义为subclass，同样B中也有性别分类，这样指定性别分类后LEfSe就开始运行step2了。使用wlx检验A，B中全部子分组的差异，并且全部子分类差异符合整体A，B之间差异的才会定义为差异，这样就排除了性别的影响。当然我们一般是没有的。用不到。

这么看来，LEfSe原始是相当简单了，也就是说对不同class之间做kruskal.test检验，在不考虑subclss的情况下做LDA判别即可。

文章后续，为了评估R语言实现的lefse和python版本是否相同，后来我用同一批数据进行了测试，十分吻合。那么就意味着你可以使用啦。

如果需要，还可以用R语言做树状图：

R 语言实现LEfSe分析

导入数据（数据在微生信生物任意一个群获取，代码就在下面了）

付费文档是不能复制粘贴的，有需要的请大家支持小编文档，如果不想码字，直接私信我，发给大家源代码。

otu = read.csv("./data_normalized.csv",row.names = 1)
head(otu)
map = read.delim("./map.txt",row.names = 1)
head(map)
map = map[colnames(otu),]
claslbl= map$SampleType
#转置otu表格

这份数据只有otu水平，大家可以换为所有层次的都囊括的数据，当然我已经测试过了，这里为了减少数据，就放了最简单的。

使用kruskal.test做非参数检验,记得的是不要忘记做p值的矫正。

otu = t(otu)
#KW rank sum test
dat3t = otu
rawpvalues <- apply(otut 2, function(x) kruskal.test(x, claslbl)$p.value);

#--得到计算后得到的p值
ord.inx <- order(rawpvalues);
rawpvalues <- rawpvalues[ord.inx];
clapvalues <- p.adjust(rawpvalues, method ="fdr");
# p.adjust
dat3t <- dat3t[,ord.inx];
dim(dat3t)

使用lda做判别

wil_datadf <- as.data.frame(dat3t);
head(wil_datadf)
#if no subclass within classes then no wilcoxon rank sum test
#Linear Discriminant analysis (LDA)
library( MASS)
ldares <- lda(claslbl~ .,data = wil_datadf);
ldares

很简单的两个过程我们就是先了lefse的核心过程。

那么如何得到统计量呢？lefse中对于不同分组的统计量如何得到呢？别急：

class_no <<- length(unique(claslbl));
ldamean$max <- apply(ldamean[,1:class_no],1,max);
ldamean$min <- apply(ldamean[,1:class_no],1,min);

#---计算LDA
ldamean$LDAscore <- signif(log10(1+abs(ldamean$max-ldamean$min)/2),digits=3);

ldamean$Pvalues <- signif(rawpvalues,digits=5);
ldamean$FDR <- signif(clapvalues,digits=5);
resTable <- ldamean;
head(ldamean)

结果展示

gut left palm right palm tongue max min LDAscore Pvalues FDR
New_CleanUp_ReferenceOTU667 68318.794 0.000 7018.882 0 68318.79 0.000 4.53 6.2526e-06 0.00021183
New_CleanUp_ReferenceOTU827 44729.247 0.000 8710.801 0 44729.25 0.000 4.35 6.7947e-06 0.00021183
`975306` 29568.891 0.000 4359.095 0 29568.89 0.000 4.17 7.7472e-06 0.00021183
`365628` 64100.199 957.121 6538.642 0 64100.20 0.000 4.51 7.8587e-06 0.00021183
New_CleanUp_ReferenceOTU793 19126.207 0.000 2903.600 0 19126.21 0.000 3.98 8.5894e-06 0.00021183
`968675` 2402.696 275675.863 109391.715 1532526 1532525.64 2402.696 5.88 1.1473e-05 0.00021183

但是这个结果你们有没有看到是有点问题的，尤其是LDA统计量，这里我门做的很粗鲁，求取最大最小值的差的绝对值，这样的信息只能得到LDA统计量。也就是otu的重要性特征，但是这个柱状图就不能得到了：

所以我们需要在这里区分一下统计量主要是由那个分组贡献的。

我们只需要进行这样一个定义即可：

a = rep("A",length(ldamean$max))
i = 1
for (i in 1:length(ldamean$max)) {
name =colnames(ldamean[,1:class_no])
a[i] = name[ldamean[,1:class_no][i,] %in% ldamean$max[i]]
}
ldamean$class = a
head(ldamean)

我们标记后就可以使用得到类似的出图了

> head(ldamean)
gut left palm right palm tongue max min LDAscore Pvalues FDR class
New_CleanUp_ReferenceOTU667 68318.794 0.000 7018.882 0 68318.79 0.000 4.53 6.2526e-06 0.00021183 gut
New_CleanUp_ReferenceOTU827 44729.247 0.000 8710.801 0 44729.25 0.000 4.35 6.7947e-06 0.00021183 gut
`975306` 29568.891 0.000 4359.095 0 29568.89 0.000 4.17 7.7472e-06 0.00021183 gut
`365628` 64100.199 957.121 6538.642 0 64100.20 0.000 4.51 7.8587e-06 0.00021183 gut
New_CleanUp_ReferenceOTU793 19126.207 0.000 2903.600 0 19126.21 0.000 3.98 8.5894e-06 0.00021183 gut
`968675` 2402.696 275675.863 109391.715 1532526 1532525.64 2402.696 5.88 1.1473e-05 0.00021183 tongue

到此我们已经完全实现了lefse的核心过程，这也是为R语言正名了大家踊跃尝试吧。

为了验证这个过程和lefse是否具有差异，我对同一组数据进行了两种途径的lefse，发现差异吻合度十分不错：左侧是R语言的结果，后侧是在terminal中使用python脚本lefse的结果，发现显著的和LDA值最大的都基本是一样的。当然一样了，因为本来就是调用的R语言的工具进行的分析，现在终于回归R语言了。

写在后面

完成核心功能后使用ggplot很容易实现柱状图的绘制。目前microbiomeviz可以实现树形图的可视化，这部分大家安装试试，看看是否可以出来呢？如果成功，记得留言告诉我。

如果本篇推文点击在看超过30个，我将柱状图绘制和树图绘制代码直接贴出来。反正小编已经快完成了。