利用MicrobiotaProcess完成LEFse分析（近似）
2021-12-15 本文已影响0人
leoxiaobei
对于Lefse来说，首先测试所有特征在不同类中的值是否存在差异分布。其次，对显著不同的特征进行事后检验，不同类中子类之间的所有两两比较应都明显符合类趋势。最后，利用线性判别分析(LDA)或随机森林(RF)对显著判别特征进行评估。
与Lefse不同，在MicrobiotaProcess中，
①省去Lefse分析前整理数据的步骤，可以直接使用physeq对象进行操作，节省大量时间
②每一步中使用的检验方法都可以自定义：
第一步中可以是oneway.test，也可以是kruskal_test
第二步中可以是wilcox.test，也可以是lm
第三步可以是lda，也可以是rf
③Lefse中由于LDA效应大小是通过随机重采样计算的，因此不能保证结果完全可以重现，MicrobiotaProcess中引入随机种子设置以实现结果的重复性
④第一步检验中提供p值校正策略，过滤方法也可以选择pvalue或是padjust。
library(coin) # for the kruskal_test and wilcox_test
library(MicrobiotaProcess)
library(phyloseq)
load("1.rdata")#加载physeq对象
set.seed(1024)#由于LDA效应大小是通过随机重采样计算的，因此应设置随机种子以实现结果的重复性
deres <- diff_analysis(obj = physeq,
classgroup = "Category",#分组
firstcomfun = "kruskal_test",
padjust = "fdr",#p值校正方法
filtermod = "pvalue",#以pvalue列过滤
firstalpha = 0.05,
strictmod = TRUE,#是否进行一对一事后检验
secondcomfun = "wilcox_test",
secondalpha = 0.01,
mlfun = "lda",#线性判别分析，可选随机森林
ldascore=3#线性判别分数
)
deres
# The original data: 253 features and 14 samples
# The sample data: 1 variables and 14 samples
# The taxda contained 664 by 7 rank
# after first test (kruskal_test) number of feature (pvalue<=0.05):61
# after second test (wilcox_test and generalizedFC) number of significantly discriminative feature:31
# after lda, Number of discriminative features: 19 (certain taxonomy classification:16; uncertain taxonomy classication: 3)
ggdiffclade画优势物种进化分支图
#visualization of different results by ggdiffclade
ggdiffclade(obj=deres,
layout="circular",#布局类型
alpha=0.3, #树分支的背景透明度
linewd=0.2, #树中连线粗细
skpointsize=0.8, #树骨架中国点的大小
taxlevel=2, #要展示的树分支水平2（门）及3以下（纲目科属种）
cladetext=4,#文本大小
setColors=F,#自定义颜色
removeUnknown=T,#不删分支，但移除分类中有un_的物种注释
reduce=T)+ # 移除分类不明的物种分支和其注释，为T时removeUnknown参数无效
scale_fill_manual(values=c("#00AED7", "#FD9347")) +
guides(color = guide_legend(keywidth = 0.1, keyheight = 0.6,order = 3,ncol=1)) +
theme(panel.background=element_rect(fill=NA),
legend.position="right",
plot.margin=margin(0,0,0,0),
legend.spacing.y=unit(0.02, "cm"),
legend.title=element_text(size=12),
legend.text=element_text(size=10),
legend.box.spacing=unit(0.02,"cm"))
ggdiffbox画LDA效应图带物种丰度箱线图
# visualization of different results by ggdiffbox
diffbox <- ggdiffbox(obj=deres, box_notch=FALSE,
colorlist=c("#00AED7", "#FD9347"),
l_xlabtext="relative abundance")
diffbox
ggeffectsize放大LDA效应图
effectsize <- ggeffectsize(obj=deres, lineheight=0.1,linewidth=0.3,pointsize=3) +
scale_color_manual(values=c("#00AED7", "#FD9347"))
effectsize
ggdifftaxbar可视化每个LDA显著性物种在样本的分布情况
ggdifftaxbar(obj=deres, xtextsize=1.5, output="each_biomarkder_barplot",coloslist=c("#00AED7", "#FD9347"))
示例
PS：MicrobiotaProcess新版本已全面支持MPSE流操作，但mp_diff_analysis函数同样参数计算出来有OTU而不只是物种等级，尽管具有物种等级数目两个函数都是19个，示例如下：
sss <- as.MPSE(physeq)
sss %<>% mp_rrarefy()
sss %<>% mp_cal_abundance(.abundance = RareAbundance,force=T) %>% mp_cal_abundance(.abundance=RareAbundance,.group=Category,force = T)
sss %<>% mp_diff_analysis(
.abundance = RelRareAbundanceBySample,
.group = Category,
p.adjust="fdr",
filter.p="pvalue",
first.test.alpha = 0.05,
second.test.alpha = 0.01,
ldascore = 3
)
taxa.tree <- sss %>% mp_extract_tree(type="taxatree")
taxa.tree %>% select(label, nodeClass, LDAupper, LDAmean, LDAlower, Sign_Category, pvalue, fdr) %>% filter(!is.na(fdr))
#从第八行开始
# # A tibble: 26 x 8
# label nodeClass LDAupper LDAmean LDAlower Sign_Category pvalue fdr
# <chr> <chr> <dbl> <dbl> <dbl> <chr> <dbl> <dbl>
# 1 OTU_80 OTU 3.32 3.28 3.23 M 0.00671 0.0975
# 2 OTU_17 OTU 4.14 4.12 4.10 N 0.00195 0.0759
# 3 OTU_184 OTU 3.16 3.12 3.06 N 0.00298 0.0759
# 4 OTU_4 OTU 4.06 4.03 3.99 M 0.00671 0.0975
# 5 OTU_24 OTU 3.93 3.90 3.88 M 0.00671 0.0975
# 6 OTU_289 OTU 3.76 3.70 3.62 M 0.00192 0.0759
# 7 OTU_71 OTU 3.68 3.62 3.54 N 0.00671 0.0975
# 8 p__ Bacteroidetes Phylum 5.15 5.13 5.11 M 0.00451 0.0856
# 9 p__ Firmicutes Phylum 4.99 4.97 4.94 N 0.00451 0.0856
# 10 c__ Bacteroidia Class 5.15 5.13 5.11 M 0.00451 0.0856
# # ... with 16 more rows
其次随后的mp_plot_diff_res画图功能还有所欠缺，尤其是cladogram中无法有效去除冗余物种以及显示LDA差异物种具体等级信息，如下：
#mp_plot_diff_res函数以分步呈现
#plot treeskeleton
p1 <- ggtree(taxa.tree,layout="radial",size = 0.3) +
geom_point(
data = td_filter(!isTip),
fill="white",
size=1,
shape=21)
# display the high light of phylum clade.
p2 <- p1 +
geom_hilight(mapping = aes(subset = nodeClass == "Phylum",
node = node,
fill = label))
# display the relative abundance of features(OTU)
#按样本
p3 <- p2+
ggnewscale::new_scale("fill") +
geom_fruit(
data = td_unnest(RareAbundanceBySample),
geom = geom_star,
mapping = aes(x = fct_reorder(Sample, Category, .fun=min),
size = RelRareAbundanceBySample,
fill = Category,
subset = RelRareAbundanceBySample > 0.05),
starshape = 13,
starstroke = 0.25,
offset = 0.04,
pwidth = 0.8,
grid.params = list(linetype=2)) +
scale_size_continuous(name="Relative Abundance (%)",range = c(1, 3)) +
scale_fill_manual(values=c("#1B9E77", "#D95F02"))
p3
# display the tip labels of taxa tree
p4 <- p3 +
geom_tiplab(size=2, offset=7.5)
p4
# display the LDA of significant OTU.
p5 <- p4 +
ggnewscale::new_scale("fill") +
geom_fruit(
geom = geom_col,
mapping = aes(x = LDAmean,
fill = Sign_Category),
orientation = "y",
offset = 0.4,
pwidth = 0.5,
axis.params = list(axis = "x",
title = "Log10(LDA)",
title.height = 0.01,
title.size = 2,
text.size = 1.8,
vjust = 1),
grid.params = list(linetype = 2)
)
p5
# display the significant (FDR) taxonomy after kruskal.test (default)
p6 <- p5 +
ggnewscale::new_scale("size") +
geom_point(
data=td_filter(!is.na(fdr)),
mapping = aes(size = -log10(pvalue),
fill = Sign_Category),
shape = 21) +
scale_size_continuous(range=c(1, 3)) +
scale_fill_manual(values=c("#1B9E77", "#D95F02"))
p6 + theme(
legend.key.height = unit(0.3, "cm"),
legend.key.width = unit(0.3, "cm"),
legend.spacing.y = unit(0.02, "cm"),
legend.text = element_text(size = 7),
legend.title = element_text(size = 9),
)
按样本
#修改p3和p4，其余不变
# display the relative abundance of features(OTU)
# 按分组
p3 <- p2+
ggnewscale::new_scale("fill") +
geom_fruit(
data = td_unnest(RareAbundanceByCategory),
geom = geom_star,
mapping = aes(x=Category,
size = RelRareAbundanceByCategory,
fill = Category,
subset = RelRareAbundanceByCategory > 0.05),
starshape = 13,
starstroke = 0.25,
offset = 0.04,
pwidth = 0.1,
grid.params = list(linetype=2)) +
scale_size_continuous(name="Relative Abundance (%)",range = c(1, 3)) +
scale_fill_manual(values=c("#1B9E77", "#D95F02"))
p3
# display the tip labels of taxa tree
p4 <- p3 +
geom_tiplab(size=2, offset=2,mapping = aes(subset=!is.na(fdr))
p4
按分组
目前来看，后一个mp_diff_analysis函数所提供的的画图功能有所不足，其进化分支图主要展示的为OTU，这也解释了其结果中为什么会出现OTU，而不是仅有物种等级，推荐使用传统diff_analysis方法。
