微生物组数据再也别用相对丰度了，结合RastStat完美展示OTU

写在前面

题目错误，更正为：EasyStat包。我们在分析群落数据的时候，往往相对丰度转化后就展示了，在差异分析中往往使用相对丰度做差异分析不是上上策，因为目前edger和Desep2包中的标准化方式可能更加准确的评估序列的真实丰度。下面我们以edger包中的TMM标准化为例，使用我的EasyStat，共同展示部分微生物的丰度。相当完美。

实战

导入R包和功能

#-差异OTU分析作图
# library(EasyMicrobiome)
library(edgeR)
library(dplyr)
# #使用案例
library(EasyStat)
library(ggplot2)
# 导入功能
vegan_otu <- function(physeq){
OTU <- otu_table(physeq)
if(taxa_are_rows(OTU)){
OTU <- t(OTU)
}
return(as(OTU,"matrix"))
}

导入数据 指定OTU

OTU <- c("ASV_1700", "ASV_1574", "ASV_479", "ASV_1953")

ps = readRDS("../EasyMicrobiome/ori_data/ps_liu.rds")
ps

edger 包 的TMM方法标准化

otu_table = as.data.frame(vegan_otu(ps))
count = as.matrix(otu_table)
count <- t(count)
head(count)
#数据整理形式一######
sub_design <- as.data.frame(sample_data(ps))
dim(sub_design)
sub_design$SampleType = as.character(sub_design$Group)
sub_design$SampleType <- as.factor(sub_design$SampleType)

# create DGE list
d = DGEList(counts=count, group=sub_design$SampleType)
d$samples
d$counts
d = calcNormFactors(d)#默认为TMM标准化
# 生成实验设计矩阵
design.mat = model.matrix(~ 0 + d$samples$group)
colnames(design.mat)=levels(sub_design$SampleType)
d2 = estimateGLMCommonDisp(d, design.mat)
d2 = estimateGLMTagwiseDisp(d2, design.mat)
fit = glmFit(d2, design.mat)
otu = as.matrix(fit$fitted.values)
dim(otu)

重构phyloseq对象 生成EasyStat包输入表格

ps_sub1 <- phyloseq(otu_table(otu, taxa_are_rows=TRUE),
sample_data(ps),
tax_table(ps)

)
ps_sub1
ps_sub <- ps_sub1 %>%
subset_taxa(
row.names(tax_table(ps))%in% OTU
)
ps_sub
# ps_sub = transform_sample_counts(ps_sub, function(x) x / sum(x) );ps_sub

otu = as.data.frame(t(otu_table(ps_sub)))
head(otu)
map = as.data.frame(sample_data(ps_sub))
map = map[,1]
otuGroup = merge(otu,map,by= "row.names",all = F)
head(otuGroup)
colnames(otuGroup)[1] = "ID"
#-------
otuGroup <- select(otuGroup, ID, Group,everything())
colnames(otuGroup)[2] = "group"

然后两条函数即可完成完美展示

result = MuiKwWlx(data = otuGroup,num = c(3:6))
result
#
result1 = FacetMuiPlotresultBox(data = otuGroup,num = c(3:6),result = result,sig_show ="abc",ncol = 2 )
p <- result1[[1]]
p

使用柱状图展示

result1 = FacetMuiPlotresultBar(data = otuGroup,num = c(3:6),result = result,sig_show ="abc",ncol = 2 )
p <- result1[[1]]
p
ggsave("./1.png",p,width = 6,height = 6)

改良柱状图

result1 = FacetMuiPlotReBoxBar(data = otuGroup,num = c(3:6),result = result,sig_show ="abc",ncol = 2 )
p <- result1[[1]]
p
ggsave("./2.png",p,width = 6,height = 6)