LEfSe分析的在线+本地运行的超详细教程
LEfSe概述
LEfSe通过将具有统计学意义的标准检验与表征生物学一致性和效应相关性的附加检验结合起来，确定最有可能解释类别之间差异的特征（如物种、进化枝、基因或功能，统称为biomaker）。类别比较方法通常用于预测biomaker，这些标记由违反类别之间无差异零假设的特征组成；并通过检测特征子集的丰度模式，估计显著变化的程度；效应大小提供了每个特征对类别差异贡献大小的估计。
以下是对LEfSe工作原理的简单概括，细节部分可参考Segata等（2011）的论文。
Data：输入数据包含m个样本（列）的集合，每个样本由n个变量（行，通常进行行归一化，红色代表高值，绿色代表低值）组成；每个样本都对应了一些类别（即样本分组），代表正在研究的主要生物学假设，可以为两组或多组，也可以包含反映类内分组的子类别标签。
Step1：通过Kruskall-Wallis检验分析所有变量，检验不同类别中的值是否存在差异分布。
Step2：对于Step1中违反零假设的变量，即Kruskall-Wallis检验显著的结果，继续使用成对的Wilcoxon检验检验不同类别内的子类别之间所有成对比较是否与类别水平趋势显著一致。
Step3：所得的向量子集用于构建线性判别分析（LDA）模型，利用类间的相对差异对变量进行降维和分类，输出包含一列与类别有关的特征列表，这些特征与类别中的子类别分组保持一致，并根据它们区分类别的效应大小进行排名。
因此LEfSe是用于对不同生物学方面的相关性进行排名以及进行进一步研究和分析的有价值的工具。它强调统计意义和生物相关性，能够在组与组之间寻找具有统计学差异的biomaker。在该方法中引入先验知识作为监督，解决了传统上与高维数据挖掘有关的挑战。
接下来，展示LEfSe工具的使用，以及对结果部分的说明。
LefSe输入文件格式
LEfSe的输入文件要求将微生物丰度表与分组信息合并，整理为如下所示的样式。
每一列代表一个样本，对于文件的前几行：
class（必须存在），各样本的主分组名称；
sub_class（可缺省），各样本的次分组名称；
subject_id（可缺省），代表编号。
之后的各行，为各微生物类群的层级划分，首先是分类水平的最高级，依次往下推，不断的增加分类水平，不同的分类水平需要用“|”隔开。在各样本中对应的数值，为相应微生物类群在该样本中的总相对丰度。
对于LEfSe识别的输入格式，可以手动使用Excel整理微生物丰度表（就是有点麻烦），也可自写代码完成。如下命令行展示了一个从OTU表到LEfSe输入格式的转换过程，代码可见网盘附件“OTU表转化为LEfSe输入格式的脚本”。
#OTU 表到 LEfSe 输入格式的转换过程##测试文件说明
#otu_table.txt，OTU 丰度表，绝对丰度相对丰度都可以
#otu_table_tax_assignments.txt，OTU 注释文件，我用 silva 注释的（大家随意）
#group.txt，样本分组文件
#1-summarize_trans.py、2-lefse_trans.py，两个自备的脚本
##首先借助 qiime 工具，从 OTU 水平的丰度表获得“界门纲目科属”水平的微生物类群丰度表
#可使用 conda 安装 qiime 环境
#conda create -n qiime1 qiime
#激活 qiime 环境
source activate qiime1
#在 otu_table.tsv 开头添加一行“# Constructed from biom file”，以便将 otu_table.tsv 转为 qiime 可识别样式
cp otu_table.txt otu_table.tsv
sed -i '1i\# Constructed from biom file' otu_table.tsv
#otu_table.tsv 转换为 biom 格式
biom convert -i otu_table.tsv -o otu_table.biom --table-type="OTU table" --to-json
#添加 otu 注释信息至 biom 格式的 otu 表（otu_table.biom ）的最后一列，并将列名称命名为 taxonomy
biom add-metadata -i otu_table.biom --observation-metadata-fp otu_table_tax_assignments.txt -o otu_table.silva.biom --sc-separated taxonomy --observation-header OTUID,taxonomy
#使用 qiime 实现对 otu 表的“界门纲目科属”统计
summarize_taxa.py -i otu_table.silva.biom -o taxonomy
#退出 qiime 环境
source deactivate qiime1
##接下来，使用自写脚本合并“界门纲目科属”统计的类群，并转换为 LEfSe 的输入格式
#合并“界门纲目科属”统计的类群，并排序，结束后 taxonomy 路径下获得“otu_table.silva_all.txt”
python2 1-summarize_trans.py -i taxonomy --prefix otu_table.silva
#过滤去除注释为“Other”的微生物类群
grep -v 'Other' taxonomy/otu_table.silva_all.txt > otu_table.silva_filt.txt
#将“otu_table.silva_filt.txt”转化为 LEfSe 的输入格式，并合并分组（group.txt 文件）
python2 2-lefse_trans.py otu_table.silva_filt.txt group.txt otu_table_for_lefse.txt
#最后的“otu_table_for_lefse.txt”就是了
#该脚本只支持添加 class 组，若存在 sub_class 组可后续手动添加即可
#就是过程有点繁琐......主要懒得自写代码统计“界门纲目科属”，就用 qiime 来做，就烦 qiime 必须识别 biom 文件，还得来回转格式
#后续哪天闲了重新整个不借助 qiime 的
下文开始展示LefSe在线运行过程，使用的示例数据集下载自（我也放网盘附件里了，lefse_test.txt）：
https://bitbucket.org/biobakery/biobakery/raw/tip/demos/biobakery_demos/data/lefse/input/hmp_small_aerobiosis.txt
别问我为什么没用上面自己转化的数据，因为下文是先写的，上文代码是后写的……
由于我没找到关于该数据集的文字描述，所以下述几行文字是我根据它的分组名称猜的……
oxygen_availability，测量了样本来源的微生物群落的氧利用能力，并据此将样本归为3个类别：高氧利用（High_O2）、中氧利用（Mid_O2）以及低氧利用（Low_O2）。
body_site，根据样本来源，将样本归为6个类别：ear、oral、gut、nasal、vagina、skin。
subject_id，代表编号。
Bacteria等，为各类微生物类群的名称，及其在各样本中的相对丰度信息。
LEfSe在线运行
LefSe提供了在线分析方法，作为Galaxy平台的一个分析模块，方便了我们使用。
Galaxy平台：http://huttenhower.sph.harvard.edu/galaxy
下文的在线运行结果，已存放在了网盘附件中，见 “LEfSe在线运行示例的结果文件”。
上传待分析的数据
首先，把已经生成好的LEfSe配置文件上传到云平台上。
上传成功后，进入LEfSe分析界面，设置参数并运行LEfSe分析，以及可视化。
在线版的LEfSe由执行以下步骤的六个模块组成，接下来我们一步步分析。
（A）设置LEfSe的数据格式
A模块中的可选项主要包括指定样本的分组类别行、子类别行等，以及是否需要标准化数据。
执行后，LEfSe将预处理输入文件的格式，获得LEfSe内置识别文件格式。
（B）LEfSe分析
这一步是LEfSe的计算环节。
使用A模块中提供的预处理后的文件，设置程序运行参数执行分析。计算过程概要可见上文“LEfSe原理概述”中的描述。
主要参数中，包括了Kruskall-Wallis检验的p值阈值、Wilcoxon检验的p值阈值、以及线性判别得分（LDA得分）的阈值，用作识别显著的biomaker的标准。
该步结束后，就获得LEfSe的计算结果了。
结果下载后可用excel打开，共5列信息，包含了每种微生物类群名称及其对类间差异贡献效应的大小。
第1列，每种微生物类群名称；
第2列，每种微生物类群在各分组类别中丰度平均值中最大值的log10，如果平均丰度小于10的按照10来计算（即该列的最小值为1）；
第3列，若该微生物类群存在组间显著差异，且线性判别得分高于设定阈值，则该列展示其所富集的分组类别名称；若不存在显著差异，则该列为空；
第4列，若该微生物类群存在组间显著差异，且线性判别得分高于设定阈值，则该列展示其线性判别得分值，用以评估其对观测到的组间差异的效应大小，该值越高代表该微生物类群越重要；若不存在显著差异，则该列为空；
第5列，统计检验的p值，若显著，则将这些微生物作为解释类别之间差异的biomaker来看待；若不显著，则该列值为“-”；
该步的结果可进一步为可视化模块（C、D、E、F）提供输入。
（C）绘制LEfSe得分值
以图形方式展示了发现的生物标志物（模块B的输出）及其效应大小。
查看结果文件，以柱形图的形式展示了显著的微生物类群，颜色代表了其富集的分组。
数值为线性判别得分（已经过log10转化），代表了该类群的效应大小，值越高代表了其对观测到的组间差异贡献越大，表明该微生物类群是更重要的biomaker。
（D）绘制进化分支图
以图形方式展示了由分层要素名称指定的分类树中发现的biomaker（模块B的输出）。
查看结果文件，该结果就是常见的那种“物种进化树”样式。
图中每一个节点代表一种给定的微生物类群，对于不显著的微生物使用黄色节点标注，显著的则赋予其它颜色，颜色代表了其富集的分组，即表明这些分支在某些分组中具有更高的丰度。
如果觉得该图中，黄色（不显著）的点过于密集影响观测，即觉得重叠严重，则可以通过编辑模块B的输出结果，将其中的一部分不显著的微生物类群去除。之后再将编辑后的表读入Galaxy，作为LEfSe模块D的输入，如此出图中黄色的点就少了，图片就会“整齐”很多。
例如，我删掉了模块B输出结果中的一些内容，然后将修改过后的文件重新上传至Galaxy，同上述，通过左侧界面的“Get Data”的“Upload File”上传，上传时注意选择文件格式为“lefse_internal_res”。
然后在“LEfSe”中选择“D） Plot Cladogram”，读取新上传后的过滤后的文件，重新作图。示例结果如下，节点变少了，图片“整齐”很多。（不过我好像删的数据有点多……大家视情况来吧）
（E）绘制单个变量的分布
该步将模块A、B的输出作为输入，将某特定变量（手动指定）的行值绘制为一个包含类和子类结构的丰度直方图。
选择展示的变量可以为已识别的biomaker或非biomaker，即可以为显著的也可以为非显著的，通过在绘图选项中指定。
例如选择了其中某微生物类群进行观测，结果图中以柱形图展示了该微生物在各样本中的丰度信息，颜色代表了样本的分组状况。
（F）绘制差异特征
与模块E相比，该步可一次将所有变量的行值统一绘制为包含类和子类结构的丰度直方图，结果以zip压缩包形式保存。
可以只输出已识别为biomaker的变量，也可以输出所有变量，通过在绘图选项中指定。
例如上述只选择了识别为biomaker（即显著）的变量进行展示，可在结果中下载打包好的压缩包。下图是其中之一的展示图。
LEfSe本地运行
然后是本地版LEfSe工具的使用，可选项更多，更加灵活。
LEfSe目前也已打包在conda中，方便了我们安装使用。如下继续使用上述示例数据，展示本地版LEfSe工具的使用，shell命令行。
本地运行的输出结果示例，也放在了网盘附件中，见 “LEfSe本地运行示例的结果文件”。由于我在本地运行和在线版运行时的设定参数不同，所以两个示例输出是不一致的，所以不要疑问为啥二者不一致。
各结果文件的说明，参考上文在线版的输出示例即可。
#conda 安装 LEfSeconda install -c biobakery lefse
#或者
conda install -c bioconda lefse
##对应于上述 A-F 6 个模块，本地版的命令行操作示例如下
#A，设置 LEfSe 的数据格式，详情 format_input.py -h
#-c，指定 class 的行（必须指定）；-s，指定 sub_class 的行（可缺省）；-u，指定 subject_id 的行（可缺省）；-o，设置归一化值，默认 -1 即不执行标准化
#注：版本问题，有时 format_input.py 命令找不到，可能为 lefse-format_input.py
format_input.py lefse_test.txt A_lefse_test.in -c 1 -s 2 -u 3 -o 1000000
#B，LEfSe 分析，详情 run_lefse.py -h
#-l 2.0，设定 LDA 得分的对数值的最低阈值为 2
run_lefse.py A_lefse_test.in B_lefse_test.res -l 2.0
#备注：对于 B 的输出，我们可以选择从中删一些不必要（不显著）的数据，以增强 C、D、E 作图时的美感
#C，绘制 LEfSe 得分值，详情 plot_res.py -h
#注：版本问题，有时 plot_res.py 命令找不到，可能为 lefse-plot_res.py
plot_res.py B_lefse_test.res C_lefse_test.lda.pdf --format pdf --dpi 150 --width 16
#D，绘制进化分支图，详情 plot_cladogram.py -h
#注：版本问题，有时 plot_cladogram.py 命令找不到，可能为 lefse-plot_cladogram.py
plot_cladogram.py B_lefse_test.res D_lefse_test.cladogram.pdf --format pdf --dpi 150
#E，单张图的展示略，直接使用 F 绘制所有的图
#F，绘制差异特征，详情 plot_features.py -h
#注：版本问题，有时 plot_features.py 命令找不到，可能为 lefse-plot_features.py
mkdir -p F_out
plot_features.py A_lefse_test.in B_lefse_test.res F_out/lefse_test --format pdf --dpi 200
参考资料
https://bitbucket.org/biobakery/biobakery/wiki/lefse
Segata N, Izard J, Waldron L, et al. Metagenomic biomarker discovery and explanation. Genome biology, 2011, 12(6).
R语言的多组间非参数差异分析Kruskal-Wallis检验