❝本节来介绍如何R来自定义构建函数来进行多样性指数计算及绘图；
加载R包
代码语言：javascript复制library(tidyverse)
library(vegan)
library(rstatix)
library(ggpubr)
library(magrittr)
导入数据
代码语言：javascript复制# 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"))
定义函数计算多样性指数
代码语言：javascript复制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)
}
导出指数计算结果
代码语言：javascript复制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
整理绘图文件
代码语言：javascript复制❝在此选取我们需要的数据来进行展示 ❞
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")
定义绘图函数
代码语言：javascript复制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"))
}
数据可视化
代码语言：javascript复制make_plot(df,group,value,col)
