本节来介绍如何R来自定义构建函数来进行多样性指数计算及绘图
加载R包
library(tidyverse)
library(vegan)
library(rstatix)
library(ggpubr)
library(magrittr)
导入数据
# otu数据表转置，行为样本，列为OTU编号
alpha <- read.delim("otu_taxa_table-2.xls",sep="\t",row.names = 1) %>%
t() %>% as.data.frame()
# 分组文件
group <- read_tsv("group.xls") %>% set_colnames(c("sample","group"))
定义函数计算多样性指数
alpha_diversity <- function(x,y) {
Shannon <- diversity(x, index = 'shannon')
Simpson <- diversity(x, index = 'simpson')
observed_species <- specnumber(x)
Chao1 <- estimateR(x)[2,]
ACE <- estimateR(x)[4,]
pielou <- diversity(x,index = "shannon")/log(specnumber(x),exp(1))
result <- data.frame(Shannon,Simpson,observed_species,Chao1,ACE,pielou) %>%
rownames_to_column("sample") %>%
left_join(.,y,by="sample")
return(result)
}
导出指数计算结果
alpha_diversity(alpha,group) %>% write.table(file="alpah.xls",sep="\t",quote = F,row.names = F)
sample Shannon Simpson observed_species Chao1 ACE pielou group
1 A1 6.022632 0.9898521 2061 2676.254 2715.836 0.7892378 A
2 A2 5.197605 0.9761031 1516 2181.152 2212.486 0.7096840 A
3 A3 5.199001 0.9763631 1533 2129.035 2184.918 0.7087953 A
4 A4 5.347510 0.9810318 1506 1995.025 1997.526 0.7308124 A
5 A5 5.395820 0.9857777 1362 1734.167 1765.997 0.7476843 A
6 B1 5.546399 0.9888982 1445 1917.961 1887.453 0.7623010 B
7 B2 5.541260 0.9889996 1302 1612.939 1602.964 0.7726610 B
8 B3 5.892619 0.9913167 1844 2268.215 2309.436 0.7836250 B
9 B4 6.130304 0.9940395 1852 2143.778 2156.409 0.8147643 B
10 B5 5.598551 0.9902785 1393 1854.570 1786.834 0.7733644 B
11 C1 4.255467 0.9194131 1241 1574.346 1627.599 0.5973698 C
整理绘图文件
在此选取我们需要的数据来进行展示
df <- alpha_diversity(alpha,group) %>% select(-sample,-observed_species,-Simpson) %>%
pivot_longer(-group)
# 先给定一些颜色
col <- c("#1F78B4","#33A02C","#FB9A99","#E31A1C","#FDBF6F","#B2DF8A",
"#A6CEE3","#BA7A70","#9D4E3F","#829BAB")
定义绘图函数
make_plot <- function(data,x,y,z){
ggplot(data,aes(x={{x}},y={{y}},fill={{x}}))+
geom_violin(trim=F)+
stat_boxplot(geom="errorbar",position=position_dodge(width=0.2),width=0.1)+
geom_boxplot(position=position_dodge(width =0.8),width=0.1,fill="white")+
scale_fill_manual(values={{z}})+
facet_wrap(.~name,scales = "free")+
theme_bw()+
theme(panel.spacing.x = unit(0.2,"cm"),
panel.spacing.y = unit(0.1, "cm"),
axis.title = element_blank(),
strip.text.x = element_text(size=12,color="black"),
axis.text = element_text(color="black"),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
legend.position = "non",
plot.margin=unit(c(0.3,0.3,0.3,0.3),units=,"cm"))
}
数据可视化
make_plot(df,group,value,col)
数据获取
可以看到通过自定义两个函数后，续如果再有类似的需求可以直接调用非常的方便；当然此图还可以在上面添加方差分析显著性标记，这个以后在做介绍；需要数据的请评论区留言
