使用Easystat挖掘微生物组数据-看得到的就是方便和强大
写在前面
目前我们在Easystat用于数据差异分析及其可视化方面做了很多的工作，尤其是大大扩展了可视化的类型，相信大家总能找到一种适合自己的可视化方式。现在我带领大家使用此包进行微生物组数据的差异分析，并进行可视化。我在自己的分析流程中已经加入了这篇帖子的内容，这里也分享给大家做一个应用。由于默认Easystat的显著性字母标记和常规认识相反，所以这里我写了函数ord_sig，可以直接对显著性标记字母进行顺序重排。
实战
ps对象存储在ggClusterNet包中，该包在github上：taowenmicro/ggClusterNet，如果没有安装，可以在后台回复：数据，可以获取到ps_liu.rds文件。使用函数readRDS即可导入。
#---------
library(dplyr)
# library(EasyMicrobiome)
library(phyloseq)
library(EasyStat)
library(ggClusterNet)
library(ggplot2)
ps
写了一个小函数，用于按照注释文件tax合并otu。这里写这个函数有两个原因，基于phyloseq的函数得到的输出往往存在问题，对合并之后OTU的名字也命名的不够准确，并且速度很缓慢，所以这里这个小函数写出来还是很好用的。
tax_glom_wt <- function(ps = ps,ranks = "Phylum") {
otu <- as.data.frame(t(vegan_otu(ps)))
tax <- as.data.frame(vegan_tax(ps))
# building group
split <- split(otu,tax[[ranks]])
#calculate sum by group
apply <- lapply(split,function(x)colSums(x[]))
# result chack
otucon <- do.call(rbind,apply)
taxcon <- tax[1:match(ranks,colnames(tax))]
taxcon <- taxcon[!duplicated(tax[[ranks]]),]
row.names(taxcon) <- taxcon[[ranks]]
pscon <- phyloseq(
otu_table( as.matrix(otucon),taxa_are_rows = TRUE),
tax_table(as.matrix(taxcon)),
sample_data(ps)
)
return(pscon)
}
做差异分析之前，我们使用相对丰度转化，并且合并为门水平，减少微生物的数量。其次，使用vegan_otu函数提取otu表格等，这个函数存在于ggClusterNet中，可以安装使用。简单将map文件和otu的丰度文件合并后就可以继续进行分析了。
ps_rela = transform_sample_counts(ps, function(x) x / sum(x) );ps_rela
psdata <- tax_glom_wt(ps =ps_rela ,ranks = "Phylum")
otu = as.data.frame((vegan_otu(psdata)))
map = as.data.frame(sample_data(psdata))
map$ID = row.names(map)
mapp = map[,c(ncol(map),1)]
colnames(mapp) = c("ID","group")
data = cbind(mapp,otu)
开始使用MuiKwWlx函数进行非参数检验。
result = MuiKwWlx(data = data,num = c(3:ncol(data)))
result <- ord_sig(data = result,ID = NULL)
可视化方案一
res <- MuiPlotStackBar(data = data,i =c(3:ncol(data)) ,result = result,errbar = F)
#--提取图片
p3 = res[[1]]
p3
#--提取数据
plotdata = res[[2]]
#--可视化方案2
result1 = FacetMuiPlotresultBar(data = data,num = c(3:ncol(data)),result = result,sig_show ="abc",ncol = 3 )
p <- result1[[1]]
p
#---可视化方案3
result1 = FacetMuiPlotresultBox(data = data,num = c(3:ncol(data)),result = result,sig_show ="abc",ncol = 3 )
p <- result1[[1]]
p
#-可视化方案4
result1 = FacetMuiPlotReBoxBar(data = data,num = c(3:ncol(data)),result = result,sig_show ="abc",ncol = 3 )
p <- result1[[1]]
p
#--可视化方案5
result = GroupBar(data = data,num =c(3:ncol(data)),result = result,sig_show ="abc" )
# output result
p <- result[[1]]
p