如何使用QIIME2的ancom结果进行Lefse分析?
最近学习了差异丰度分析“Lefse在线网页”的使用方法,本分析使用的原始数据是QIIME2分析过程中产生的“l6-ancom-subject.qzv”文件。由于在学习过程中未发现相关的资料讲解,便根据自己的学习情况进行了整理,因此过程中涉及的一些问题也没有解决。如果才识渊博的您碰巧读到这篇博客,又碰巧知道问题答案的话,那请您碰巧指导下我这个“生信小菜鸟”吧!(捂脸哭(T﹏T)
Lefse分析
LEfSe(Linear discriminant analysis Effect Size)一种广泛用于宏基因组生物高维数据的biomaker鉴别和解释的工具,通过分类比较、生物学一致性测试和效应大小估计,描述和验证两个或多个微生物群落之间的差异物种、基因或途径(Segata et al, 2011)。
LEfSe通过将具有统计学意义的标准检验与表征生物学一致性和效应相关性的附加检验结合起来,确定最有可能解释类别之间差异的特征(如物种、进化枝、基因或功能,统称为 biomaker)。类别比较方法通常用于预测 biomaker,这些标记由违反类别之间无差异零假设的特征组成;并通过检测特征子集的丰度模式,估计显著变化的程度;效应大小提供了每个特征对类别差异贡献大小的估计。
LEfSe分析的三大步骤
Step1:利用 Kruskal-Wallis 秩和检验检测所有的特征物种,通过检测不同组间的物种丰度差异,获得显著性差异物种。
Step2:再利用 Wilcoxon 秩和检验检测上步获得的显著性差异物种的所有亚种是否都趋同于同一分类级别。
Step3:最后用线性判别分析(LDA),得到最终的差异物种(即 biomarker)。
Lefse在线网址:http://huttenhower.sph.harvard.edu/galaxy/
[ 问题1:我的电脑一直登不上这个网站,重启过电脑,换过其他浏览器,都不行;所以这次的分析是在手机上学的;怀疑是不是防火墙把它拦截了?]
原始数据的处理
- 将“l6-ancom-subject.qzv”文件导入“qiime2 view”中查看,并下载下方的原始数据“percent-abundances.tsv”文件(位于图片左下角)。
percent-abundances.tsv - 对下载后的数据进行整理,使其可以导入Lefse在线分析网站。
1)Lefse网站的输入文件类型为 txt 或 excel 格式。
2)必有 group/class 主分类,sub_class 和 subject_id 可缺省。
3)将第一行的 percentile 删去。
4)微生物类别层级划分只保留到属。
[ 问题2:这与网页版数据类型不同,可能是后续分支图无分支的原因。]
5)每个数据都加1。如何对数据进行+1操作
[ 问题3:accom检验是减了1的,所以要加回1;原理不了解?]
其他网页版:
LEfSe的输入文件要求将微生物丰度表与分组信息合并,整理为如下所示的样式。(网站搜索的样式,可与自己的数据表进行比对。)
每一列代表一个样本;
class(必须存在),各样本的主分组名称;
sub_class(可缺省),各样本的次分组名称;
subject_id(可缺省),代表编号;
各微生物类群的层级划分,首先是分类水平的最高级,依次往下推,不同的分类水平需要用“|”隔开。在各样本中对应的数值,为相应微生物类群在该样本中的总相对丰度。
Lefse在线运行
网页界面如下:
上传数据
将标准格式的文件上传至网站:Get Data->Upload File->Choose local file->Start
LEfSe 分析-步骤 A)Format Data for LEfSe 格式转化
输出文件格式(该格式文件只是中间文件,具体意义不需要详细理解)如下:
LEfSe 分析-步骤 B) LDA Effect Size (LEfSe) 差异计算与统计
less_strict 设为2;more_strict 设为4;还可以设置为3.5。
输出结果:第一列:物种信息,即每种微生物类群名称
第二列:每种微生物类群在各分组类别中丰度平均值中最大值的log10,如果平均丰度小于 10 按照 10 来计算。(不过这不是原始表格中的丰度,而是标准化后的丰度,LEfSe自己有套标准化方法)
第三列:差异物种富集的组名称
第四列:LDA 值,用以评估其对观测到的组间差异的效应大小,该值越高代表该微生物类群越重要。
第五列:Kruskal-Wallis 秩和检验的 P 值,若显著,则将这些微生物作为解释类别之间差异的biomaker来看待;若不显著,则该列值为“-”。
LEfSe 分析 C-F 步骤均为作图
步骤 C)Plot LEfSe Results
输出结果:
结果文件以柱形图的形式展示了显著的微生物类群,颜色代表了其富集的分组。
数值为线性判别得分(已经过log10转化),代表了该类群的效应大小,值越高代表了其对观测到的组间差异贡献越大,表明该微生物类群是更重要的biomaker。默认预设值为2.0(看横坐标,只有LDA值的绝对值大于2才会显示在图中)。
其他网页版:
[ 问题4:Lefse分析柱形图为什么有的在左右两个方向,而有的只在右侧方向?]
步骤 D)Plot Cladogram(参数设置与步骤 C 类似)
输出结果:
小圆圈: 图中由内至外辐射的圆圈代表了由门至属的分类级别(最里面的那个黄圈圈是界)。不同分类级别上的每一个小圆圈代表该水平下的一个分类,小圆圈的直径大小代表了相对丰度的大小。
颜色: 无显著差异的物种统一着色为黄色,差异显著的物种,Biomarker跟随组别进行着色,红色节点表示在红色组别中起到重要作用的微生物类群,蓝色节点表示在蓝色组别中起到重要作用的微生物类群。未能在图中显示的Biomarker对应的物种名会展示在右侧,字母编号与图中对应(为了美观,右侧默认只显示门到科的差异物种)
其他网页版:
如果觉得该图中,黄色(不显著)的点过于密集影响观测,则可以通过编辑模块B的输出结果,将其中的一部分不显著的微生物类群去除。之后再将编辑后的表读入Galaxy,作为LEfSe模块D的输入,但上传数据时注意选择文件格式为“lefse_internal_res”。
步骤 E)Plot One Feature 绘制单个变量的分布(画其中一个物种在不同组样品中的柱状图)
输出结果:
图中不同分组用黑实线隔开,每组柱状图中的实线表示该组样品表达量的平均值,虚线代表该组样品表达量的中位值。
步骤 F)Plot Differential Features (绘制差异物种或基因柱状图,输入设置与 E 基本类似)
F 与 E 的区别在于,E 每次只能做一张图,F 可以绘制 biomarker(或所有的物种)柱状图。一般建议跳过步骤 E,直接做步骤 F,步骤 F 的结果一般以压缩包的形式,下载到本地解压后即可查看每个 biomarker 在不同样品中的表达柱状图。