这一节我们来根据OTU表计算样本的richness即丰富度指数，统计样本中不为0的数目累加即可，喜欢的小伙伴可以关注个人公众号R语言数据分析指南,持续分析更多优质资源，江山父老能容我，不使人间造孽钱，在此先行拜谢了！！！
点击原文链接获取数据
richness(物种丰富度指数计算)
加载所需R包
rm(list=ls())
library(tidyverse)
library(magrittr)
library(ggsci)
library(ggpubr)
计算richness
otu <- read.delim("otu_table.tsv",header = T,sep="\t",
check.names = F,row.names = 1)
richness <- NULL
for (i in seq_len(ncol(otu))) {
richness <- c(richness,which(otu[,i] !=0) %>% length())
}
names(richness) <- otu %>% colnames()
richness <- richness %<>% data.frame(sample= names(.),value = .)
richness
#p <- richness %<>% data.frame(names(.),.) %>%
# set_colnames(c("sample","value"))
group <- read.delim("group.xls",header = T,sep="\t",check.names = F)
colnames(richness)[1] <- colnames(group[1])
上述代码统计样本中read数不为0的个数进行累加
连接数据并可是化
将richness信息与样本分组信息整合传递给ggplot2进行数据可视化
full_join(richness,group) %>%
mutate(V4="richness") %>%
ggplot(aes(group,value,fill=group))+
scale_fill_nejm()+
geom_violin(trim=FALSE)+
geom_boxplot(width=0.1,fill="white")+
facet_grid(. ~V4,scales = "free",space="free_x")+
xlab(NULL)+ylab(NULL)+
theme_bw()+
theme(legend.position = "none",
axis.text.x=element_text(angle =0,hjust=0.5,vjust=0.5,
colour = "black",size=8),
strip.text.x = element_text(colour ="black",size=13))+
stat_compare_means(method = "anova",label.y =1000,label.x = 2,size=5)
