对于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方法。