歡迎關(guān)注R語言數(shù)據(jù)分析指南
本節(jié)來介紹如何在計(jì)算多樣性指數(shù)的基礎(chǔ)上來進(jìn)行顯著性標(biāo)記，可在文末找到獲取數(shù)據(jù)的方式
加載R包
library(tidyverse)
library(vegan)
library(magrittr)
library(multcompView)
導(dǎo)入數(shù)據(jù)
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"))
定義函數(shù)計(jì)算多樣性指數(shù)
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)
}
數(shù)據(jù)整理
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")
上面這些基本是上一篇文檔的內(nèi)容為了文檔結(jié)構(gòu)的完整疾棵，將其放置于此理郑；那么接下來就是本文的重點(diǎn)內(nèi)容多組之間進(jìn)行方差分析添加顯著性標(biāo)記
方差分析
p <- split(df,list(df$name))
aov_data <- data.frame()
str(p)
for(i in 1:4) {
anova <- aov(value ~ group,data=p[i] %>% as.data.frame() %>%
set_colnames(c("group","name","value")))
Tukey <- TukeyHSD(anova)
cld <- multcompLetters4(anova,Tukey)
dt <- p[i] %>% as.data.frame() %>%
set_colnames(c("group","name","value")) %>%
group_by(group,name) %>%
summarise(value_mean=mean(value),sd=sd(value)) %>%
ungroup() %>%
arrange(desc(value_mean)) %>%
as.data.frame()
cld <- as.data.frame.list(cld$`group`)
dt$Tukey <- cld$Letters
aov_data <- rbind(aov_data,dt)
}
構(gòu)建顯著性標(biāo)記數(shù)據(jù)集
df2 <- df %>% arrange(name) %>% left_join(.,aov_data,by=c("group","name"))
text <- df2 %>% group_by(group,name) %>% summarise(max(value)) %>% arrange(name) %>% ungroup() %>%
set_colnames(c("group","name","value")) %>%
left_join(.,df2 %>% select(1,2,6),by=c("group","name")) %>% distinct() %>%
mutate(value=case_when(name =="ACE" ~ value+90,
name =="Chao1" ~ value+90,
name =="pielou" ~ value +0.008,
name =="Shannon" ~ value+0.065))
由于循環(huán)構(gòu)建的為條形圖的數(shù)據(jù)坡疼，但顯著性標(biāo)記是不區(qū)分圖形的因此在此通過上面的代碼構(gòu)建箱線圖的數(shù)據(jù)控轿，由于還存在離群值因此做了過多的處理缠诅，各位觀眾老爺細(xì)細(xì)品味
定義繪圖函數(shù)
make_plot <- function(data,x,y,z){
ggplot(data,aes(x={{x}},y={{y}},fill={{x}}))+
stat_boxplot(geom="errorbar",position=position_dodge(width=0.2),width=0.2)+
geom_boxplot(position=position_dodge(width =0.2),width=0.5,outlier.shape = NA)+
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"))
}
數(shù)據(jù)可視化
make_plot(df,group,value,col)+
geom_text(data=text,aes(label=Tukey,y=value))
數(shù)據(jù)獲取
可以看到過程還是比較繁瑣的沒有直接調(diào)用R包來的方便哑蔫，但是通過此文你一定有所收獲钉寝，需要獲取數(shù)據(jù)的可以關(guān)注我的公眾號(hào)R語言數(shù)據(jù)分析指南會(huì)找到數(shù)據(jù)獲取方式
