lyao222lll的个人博客分享 http://blog.sciencenet.cn/u/lyao222lll

博文

LEfSe分析的在线+本地运行教程

已有 5272 次阅读 2020-1-15 20:59 |个人分类:Software/Pipeline/Script|系统分类:科研笔记| _生信工具, _差异分析

LEfSe分析的在线+本地运行的详细教程

 

尽管网上的教程多的是,但还是想自己整理一份。

LEfSeLinear discriminant analysis Effect Size)一种广泛用于宏基因组生物高维数据的biomaker鉴别和解释的工具,通过分类比较、生物学一致性测试和效应大小估计,描述和验证两个或多个微生物群落之间的差异物种、基因或途径(Segata et al, 2011)。

本篇以某16S高通量测序数据所得的微生物群落数据为例,展示使用LEfSe寻找组间显著差异的微生物类群,以及确定它们的重要程度。

下文中涉及的示例数据、脚本代码、示例结果等,已存放至百度盘,可能会用得到。

https://pan.baidu.com/s/1ctwyoTqT1ofygMUnTj7CPA

 


LEfSe概述


LEfSe通过将具有统计学意义的标准检验与表征生物学一致性和效应相关性的附加检验结合起来,确定最有可能解释类别之间差异的特征(如物种、进化枝、基因或功能,统称为biomaker)。类别比较方法通常用于预测biomaker,这些标记由违反类别之间无差异零假设的特征组成;并通过检测特征子集的丰度模式,估计显著变化的程度;效应大小提供了每个特征对类别差异贡献大小的估计。

以下是对LEfSe工作原理的简单概括,细节部分可参Segata2011)的论文。

1.png

Data输入数据包含m个样本(列)的集合,每个样本由n个变量(行,通常进行行归一化,红色代表高值,绿色代表低值)组成;每个样本都对应了一些类别(即样本分组),代表正在研究的主要生物学假设,可以为两组或多组,也可以包含反映类内分组的子类别标签。

Step1通过Kruskall-Wallis检验分析所有变量,检验不同类别中的值是否存在差异分布。

Step2对于Step1中违反零假设的变量,即Kruskall-Wallis检验显著的结果,继续使用成对的Wilcoxon检验检验不同类别内的子类别之间所有成对比较是否与类别水平趋势显著一致。

Step3:所得的向量子集用于构建线性判别分析LDA)模型,利用类间的相对差异对变量进行降维和分类,输出包含一列与类别有关的特征列表,这些特征与类别中的子类别分组保持一致,并根据它们区分类别的效应大小进行排名。

 

因此LEfSe是用于对不同生物学方面的相关性进行排名以及进行进一步研究和分析的有价值的工具。它强调统计意义和生物相关性,能够在组与组之间寻找具有统计学差异的biomaker。在该方法中引入先验知识作为监督,解决了传统上与高维数据挖掘有关的挑战。

接下来,展示LEfSe工具的使用,以及对结果部分的说明。

 


LefSe输入文件格式


LEfSe的输入文件要求将微生物丰度表与分组信息合并,整理为如下所示的样式。

2.png

每一列代表一个样本,对于文件的前几行:

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

别问我为什么没用上面自己转化的数据,因为下文是先写的,上文代码是后写的……

3.png

由于我没找到关于该数据集的文字描述,所以下述几行文字是我根据它的分组名称猜的……

oxygen_availability,测量了样本来源的微生物群落的氧利用能力,并据此将样本归为3个类别:高氧利用(High_O2)、中氧利用(Mid_O2)以及低氧利用(Low_O2)。

body_site,根据样本来源,将样本归为6个类别:earoralgutnasalvaginaskin

subject_id,代表编号。

Bacteria等,为各类微生物类群的名称,及其在各样本中的相对丰度信息。

 


LEfSe在线运行


LefSe提供了在线分析方法,作为Galaxy平台的一个分析模块,方便了我们使用。

Galaxy平台:http://huttenhower.sph.harvard.edu/galaxy

下文的在线运行结果,已存放在了网盘附件中,见 LEfSe在线运行示例的结果文件”。

 

上传待分析的数据


首先,把已经生成好的LEfSe配置文件上传到云平台上。

4.png

 

上传成功后,进入LEfSe分析界面,设置参数并运行LEfSe分析,以及可视化。

在线版的LEfSe由执行以下步骤的六个模块组成,接下来我们一步步分析。

5.png

 

A)设置LEfSe的数据格式


A模块中的可选项主要包括指定样本的分组类别行、子类别行等,以及是否需要标准化数据。

执行后,LEfSe将预处理输入文件的格式,获得LEfSe内置识别文件格式。

6.png

 

BLEfSe分析


这一步是LEfSe的计算环节。

使用A模块中提供的预处理后的文件,设置程序运行参数执行分析。计算过程概要可见上文“LEfSe原理概述”中的描述。

主要参数中,包括了Kruskall-Wallis检验的p值阈值、Wilcoxon检验的p值阈值、以及线性判别得分(LDA得分)的阈值,用作识别显著的biomaker的标准。

7.png

该步结束后,就获得LEfSe的计算结果了。

结果下载后可用excel打开,共5列信息,包含了每种微生物类群名称及其对类间差异贡献效应的大小。

1列,每种微生物类群名称;

2列,每种微生物类群在各分组类别中丰度平均值中最大值的log10(不过这不是原始表格中的丰度,而是标准化后的丰度,LEfSe自己有套标准化方法);

3列,若该微生物类群存在组间显著差异,且线性判别得分高于设定阈值,则该列展示其所富集的分组类别名称;若不存在显著差异,则该列为空;

4列,若该微生物类群存在组间显著差异,且线性判别得分高于设定阈值,则该列展示其线性判别得分值,用以评估其对观测到的组间差异的效应大小,该值越高代表该微生物类群越重要;若不存在显著差异,则该列为空;

5列,统计检验的p值,若显著,则将这些微生物作为解释类别之间差异的biomaker来看待;若不显著,则该列值为“-”。

8.png

该步的结果可进一步为可视化模块(CDEF)提供输入。

 

C)绘制LEfSe得分值


以图形方式展示了发现的生物标志物(模块B的输出)及其效应大小。

9.png

查看结果文件,以柱形图的形式展示了显著的微生物类群,颜色代表了其富集的分组。

数值为线性判别得分(已经过log10转化),代表了该类群的效应大小,值越高代表了其对观测到的组间差异贡献越大,表明该微生物类群是更重要的biomaker

10.png

 

D)绘制进化分支图


以图形方式展示了由分层要素名称指定的分类树中发现的biomaker(模块B的输出)。

参数部分详见网页中的帮助说明。

11.png

查看结果文件,该结果就是常见的那种“物种进化树”样式。

图中每一个节点代表一种给定的微生物类群,对于不显著的微生物使用黄色节点标注,显著的则赋予其它颜色,颜色代表了其富集的分组,即表明这些分支在某些分组中具有更高的丰度。

12.png

 

如果觉得该图中,黄色(不显著)的点过于密集影响观测,则可以通过编辑模块B的输出结果,将其中的一部分不显著的微生物类群去除。之后再将编辑后的表读入Galaxy,作为LEfSe模块D的输入,如此出图中黄色的点就少了,图片就会“整齐”很多。

例如,我删掉了模块B输出结果中的一些内容,然后将修改过后的文件重新上传至Galaxy,同上述,通过左侧界面的“Get Data”的“Upload File”上传,上传时注意选择文件格式为“lefse_internal_res”。

13.png

然后在“LEfSe”中选择“D Plot Cladogram”,读取新上传后的过滤后的文件,重新作图。示例结果如下,节点变少了,图片“整齐”很多。(不过我好像删的数据有点多……大家视情况来吧)

14.png

 

E)绘制单个变量的分布


该步将模块AB的输出作为输入,将某特定变量(手动指定)的行值绘制为一个包含类和子类结构的丰度直方图。

选择展示的变量可以为已识别的biomaker或非biomaker,即可以为显著的也可以为非显著的,通过在绘图选项中指定。

15.png

例如选择了其中某微生物类群进行观测,结果图中以柱形图展示了该微生物在各样本中的丰度信息,颜色代表了样本的分组状况。

16.png

 

F)绘制差异特征


与模块E相比,该步可一次将所有变量的行值统一绘制为包含类和子类结构的丰度直方图,结果以zip压缩包形式保存。

可以只输出已识别为biomaker的变量,也可以输出所有变量,通过在绘图选项中指定。

17.png

例如上述只选择了识别为biomaker(即显著)的变量进行展示,可在结果中下载打包好的压缩包。下图是其中之一的展示图。

18.png

 


LEfSe本地运行


然后是本地版LEfSe工具的使用,可选项更多,更加灵活。

LEfSe目前也已打包在conda中,方便了我们安装使用。如下继续使用上述示例数据,展示本地版LEfSe工具的使用。

本地运行的输出结果示例,也放在了网盘附件中,见“LEfSe本地运行示例的结果文件”。由于我在本地运行和在线版运行时的设定参数不同,所以两个示例输出是不一致的,所以不要疑问为啥二者不一致。

各结果文件的说明,参考上文在线版的输出示例即可。

#conda 安装 LEfSe
conda 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

 


LEfSe分析的包装器Koeken


上文“LefSe输入文件格式”中所展示的转换代码,过程中借助QIIME工具实现各分类水平的统计,需要事先将txt格式的OTU丰度表转换为biom格式以方便使QIIME识别。

既然已经将OTU表转换为biom格式,并也已安装了QIIMELEfSe程序,那就继续延伸一个工具,Koeken,是为LEfSe分析提供的包装器。KoekenQIIMELEfSe结合在一起,为从QIIME OTU表到微生物类群生物标志物的鉴定提供了更便捷的工作流程,输入biom格式的带注释OTU丰度表和带有分组信息的mapping文件,即可直接运行获得LEfSe分析结果。

Koekengithub链接:https://github.com/twbattaglia/koeken

下文的示例数据也可在网盘附件中获取,见“Koeken示例数据和结果”。

 

安装Koeken

#使用 pip 安装 Koeken
pip install https://github.com/twbattaglia/koeken/zipball/master
koeken.py --help
 
#或者编译安装
wget https://codeload.github.com/twbattaglia/koeken/zip/master
unzip master
cd koeken-master
python setup.py install
koeken.py --help
 
#或者 conda 安装
conda install Koeken
koeken.py --help

 

Koeken的输入文件包括OTU丰度表和样本分组文件。

同样地,OTU表转化为biom格式;样本分组信息包含在mapping文件里,mapping文件也是QIIME的标准配置文件。mapping文件的格式可参考QIIME文档:

http://qiime.org/documentation/file_formats.html#mapping-file-overview

#某试验对微生物群落添加某种处理,分为对照组和处理组
#并在处理后的第2、3、4天取样,16S高通量测序获得细菌群落数据
 
#otu_table.txt,OTU 丰度表
#otu_tax.txt,OTU 注释文件,数据库随意
#mapping.txt,QIIME 的 mapping 配置文件,含样本分组信息
 
#常见的 OTU 丰度表转为 biom 格式,并添加 OTU 注释信息
cp otu_table.txt otu_table.tsv
sed -i '1i\# Constructed from biom file' otu_table.tsv
biom convert -i otu_table.tsv -o otu_table.biom --table-type="OTU table" --to-json
biom add-metadata -i otu_table.biom --observation-metadata-fp otu_tax.txt -o otu_table.tax.biom \
    --sc-separated taxonomy --observation-header OTUID,taxonomy

 

然后以下是一个简单示例,展示Koeken运行。

分别代表时间和处理的变量,DayTreatment,已经添加在mapping文件中,分析时指定该列即可。

#我们期望在每个时间点,比较处理和对照的细菌群落差异
#并从中找到最具代表性的差异类群
 
#运行 Koeken 的一个简单示例,细节参数详见 koeken.py --help
#例如,设定统计检验显著性 p 值 0.05,LDA 得分值 2.0
koeken.py --input otu_table.tax.biom --map mapping.txt --output koeken_out \
    --class Treatment --split Day --level 6 --pval 0.05 --effect 2 --clade --image pdf

因为Koeken只是作为LEfSe的包装器使用,所以对于统计结果文件,内容结构和LEfSe标准输出一致。详情参考上文即可。

此外,对于每个时间点的比较分析,各自输出一个统计结果文件。这也是Koeken的一个优势,允许一次同时并行分析多批次的数据,提高效率;同时各批次结果也能分开给出,便于我们各自比较。



参考资料


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).

 


更多精彩,可关注个人公众号“生信小白鱼”,感谢大家支持。




http://blog.sciencenet.cn/blog-3406804-1214106.html

上一篇:R语言中的判别分析方法
下一篇:关于二代测序中duplication占比问题的探讨

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-8-8 17:42

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部