[toc]

R语言一键完成差异检测从数据到展示

单因素差异分析的完整方案

方案优点

1. 完整的差异分析思路及其R语言实现
2. 提供两种可视化方案：柱状图和箱线图；差异结果使用两种表示方法：字母进行排序表示，或者两组之间连线。

3. 可视化图形标注方差齐性结果，正态性检验结果等

4. 自动识别数据指标，自动保存分析结果，文件名称标记分析指标及其分析操作内容和方式。

主要函数解读

两种差异表示方案及其代码

#值得注意的是LSD多重比较输出的就是字母形式的结果，如果我们选择其他多重比较方法，注意提取差异显著字母格式的结果
out <- LSD.test(model,"group", p.adj="none")#进行多重比较
aa = out$group#结果显示：标记字母法 aa$group = row.names(aa)
wen1 = as.data.frame(tapply(dd,data_wt$group,mean,na.rm=TRUE)) wen2 = as.data.frame(tapply(dd,data_wt$group,sd,na.rm=TRUE))
went = cbind(wen1,wen2)
wentao = merge(aa,went, by="row.names",all=F)
colnames(wentao) = c(colnames(wentao[1:4]),"mean" ,"SD")
aa = mutate(wentao, ymin = mean - SD, ymax =  mean + SD)
a = max(aa$mean)*1.2 # ss <- round(wtx3$Pr(>F)[1],3)
p = ggplot(aa , aes(x = group, y = mean,colour= group)) +
geom_bar(aes(colour= group,fill = group),stat = "identity", width = 0.4,position = "dodge") +
geom_text(aes(label = groups,y=ymax, x = group,vjust = -0.3,size = 6))+
geom_errorbar(aes(ymin=ymin,
ymax=ymax),
colour="black",width=0.1,size = 1)+
scale_y_continuous(expand = c(0,0),limits = c(0,a))+
labs(x=paste(name_i,"of all group", sep = "_"),
y="group",
title = paste("Normality test",p1,"Homogeneity of variance",p2,"kruskal.test",sumkrusk[3],sep = ":"))
p

字母标记箱线图代码

out <- LSD.test(model,"group", p.adj="none")
aa = out$group aa$group = row.names(aa)
a = max(aa$dd)*1.2 data_box = data_wt[c(1,2,i)] colnames(data_box) = c("ID" , "group","dd" ) stat = out$groups
data_box$stat=stat[as.character(data_box$group),]$groups max=max(data_box[,c("dd")]) min=min(data_box[,c("dd")]) x = data_box[,c("group","dd")] y = x %>% group_by(group) %>% summarise_(Max=paste('max(',"dd",')',sep="")) y=as.data.frame(y) rownames(y)=y$group
data_box$y=y[as.character(data_box$group),]$Max + (max-min)*0.05 p = ggplot(data_box, aes(x=group, y=data_box[["dd"]], color=group)) + geom_boxplot(alpha=1, outlier.size=0, size=0.7, width=0.5, fill="transparent") + labs(x=paste(name_i," group", sep = "_"), y="group", title = paste("Normality test",p1,"Homogeneity of variance",p2,sep = ":"))+ geom_text(data=data_box, aes(x=group, y=y, color=group, label= stat)) + geom_jitter( position=position_jitter(0.17), size=1, alpha=0.7)+theme(legend.position="none") p if (length(unique(data_box$group))>3){    p=p+theme(axis.text.x=element_text(angle=45,vjust=1, hjust=1))}
FileName <- paste(name_i,"_aov_LSD_box", ".pdf", sep = "")
ggsave(FileName, p, width = 8, height = 8)

ggpubr + 箱线图 + 连线差异标注

wtq = levels(data_wt$group) lis = combn(levels(data_wt$group), 2)
x <-lis
my_comparisons <- tapply(x,rep(1:ncol(x),each=nrow(x)),function(i)i)
p = ggplot(data_box, aes(x=group, y=data_box[["dd"]], color=group)) +
geom_boxplot(alpha=1, outlier.size=0, size=0.7, width=0.5, fill="transparent") +
labs(x=paste(name_i,"of all group", sep = "_"),
y="group",
title = paste("Normality test",p1,"Homogeneity of variance",p2,sep = ":"))+
# geom_text(data=data_box, aes(x=group, y=y, color=group, label= stat)) +
geom_jitter( position=position_jitter(0.17), size=1, alpha=0.7)+theme(legend.position="none")+
stat_compare_means()+
stat_compare_means(comparisons=my_comparisons,label = "p.signif",hide.ns = F) # Add pairwise

p

实战

导入需要的包

{r pre analyse, echo=TRUE, message=FALSE, warning=FALSE, paged.print=FALSE}

install.packages(“tidyverse”, repo=site)

library(tidyverse)
library(agricolae)

library(car)

library(reshape2)
library(“ggplot2”)
library(“ggpubr”)

### 导入数据

{r imput data, echo=TRUE, message=FALSE, warning=FALSE, paged.print=FALSE}
# 读入实验设计
data_wt = read.table("./cs.txt", header=T, sep="\t");head(data_wt)
##数据由长变宽
data_wt = dcast(data_wt,ID +group ~ grou, value.var = "count")
#这里备注所需的数据格式
#前量列从第一列开始是ID，第二列是分组信息，剩下的列均为数据列
head(data_wt)

选择可视化方案

plot = "bar"
plot = "box"

运行函数

{r function for main, echo=TRUE, message=FALSE, warning=FALSE, paged.print=FALSE}
for (i in 3:ncol(data_wt)) {
ss <- data_wt[i]
colnames(ss) <- c(“count”)
ssgroup
xx <-  shapiro.test(ss$count) p1 <- xx[[2]] 方差齐性检验 xc <- bartlett.test(count~group,data=ss) p2 <- xc[[3]] if ( plot == “bar”) { if (p1 >.05& p2 >.05) { p1 <- round(p1,3) p2 <- round(p2,3) data_i = data_wt[i] ee <-as.matrix(data_i) dd <- as.vector(ee) name_i = colnames(data_wt[i]) model<-aov(dd ~ group, data=data_wt)#方差分析 wtx1 = summary(model) wtx2 = wtx1[[1]] wtx3 = wtx2[5]  if ( wtx3$Pr(>F)[1]< 0.05) {
out <- LSD.test(model,"group", p.adj="none")#进行多重比较，不矫正P值
aa = out$group#结果显示：标记字母法 aa$group = row.names(aa)
a = max(aa$dd)*1.2 aa wen1 = as.data.frame(tapply(dd,data_wt$group,mean,na.rm=TRUE))
wen2 = as.data.frame(tapply(dd,data_wt$group,sd,na.rm=TRUE)) went = cbind(wen1,wen2) wentao = merge(aa,went, by="row.names",all=F) colnames(wentao) = c(colnames(wentao[1:4]),"mean" ,"SD") aa = mutate(wentao, ymin = mean - SD, ymax = mean + SD) mi=c("#1B9E77" ,"#D95F02", "#7570B3","#E7298A") p = ggplot(aa , aes(x = group, y = dd,colour= group)) + geom_bar(aes(colour= group,fill = group),stat = "identity", width = 0.4,position = "dodge") + geom_text(aes(label = groups,y=ymax, x = group,vjust = -0.3,size = 6))+ geom_errorbar(aes(ymin=ymin, ymax=ymax), colour="black",width=0.1,size = 1)+ scale_y_continuous(expand = c(0,0),limits = c(0,a))+ labs(x=paste(name_i,"of all group", sep = "_"), y="group", title = paste("Normality test",p1,"Homogeneity of variance",p2,sep = ":")) p p=p+theme_bw()+ geom_hline(aes(yintercept=mean(dd)), colour="black", linetype=2) + geom_vline(aes(xintercept=0), colour="black", linetype="dashed") + # scale_fill_manual(values = mi, guide = guide_legend(title = NULL))+ theme( panel.grid.major=element_blank(), panel.grid.minor=element_blank(), plot.title = element_text(vjust = -8.5,hjust = 0.1), axis.title.y =element_text(size = 20,face = "bold",colour = "black"), axis.title.x =element_text(size = 24,face = "bold",colour = "black"), axis.text = element_text(size = 20,face = "bold"), axis.text.x = element_text(colour = "black",size = 14), axis.text.y = element_text(colour = "black",size = 14), legend.text = element_text(size = 15,face = "bold"), legend.position = "none"#是否删除图例 ) p FileName <- paste(name_i,"_aov_LSD_bar", ".pdf", sep = "_") ggsave(FileName, p, width = 8, height = 8) }else if ( wtx3$Pr(>F)[1]>= 0.05)  {
out <- LSD.test(model,"group", p.adj="none")#进行多重比较，不矫正P值
aa = out$group#结果显示：标记字母法 aa$group = row.names(aa)
a = max(aa$dd)*1.2 wen1 = as.data.frame(tapply(dd,data_wt$group,mean,na.rm=TRUE))
wen2 = as.data.frame(tapply(dd,data_wt$group,sd,na.rm=TRUE)) went = cbind(wen1,wen2) wentao = merge(aa,went, by="row.names",all=F) colnames(wentao) = c(colnames(wentao[1:4]),"mean" ,"SD") aa = mutate(wentao, ymin = mean - SD, ymax = mean + SD) ss <- round(wtx3$Pr(>F)[1],3)
mi=c("#1B9E77" ,"#D95F02", "#7570B3","#E7298A")
p = ggplot(aa , aes(x = group, y = dd,colour= group)) +
geom_bar(aes(colour= group,fill = group),stat = "identity", width = 0.4,position = "dodge") +
# geom_text(aes(label = groups,y=ymax, x = group,vjust = -0.3,size = 6))+
geom_errorbar(aes(ymin=ymin,
ymax=ymax),
colour="black",width=0.1,size = 1)+
scale_y_continuous(expand = c(0,0),limits = c(0,a))+
labs(x=paste(name_i,"of all group", sep = "_"),
y="group",
title = paste("Normality test",p1,"Homogeneity of variance",p2,"aov",ss,sep = ":"))
p
p=p+theme_bw()+
geom_hline(aes(yintercept=mean(dd)), colour="black", linetype=2) +
geom_vline(aes(xintercept=0), colour="black", linetype="dashed") +
# scale_fill_manual(values = mi, guide = guide_legend(title = NULL))+
theme(

panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),

plot.title = element_text(vjust = -8.5,hjust = 0.1),
axis.title.y =element_text(size = 20,face = "bold",colour = "black"),
axis.title.x =element_text(size = 24,face = "bold",colour = "black"),
axis.text = element_text(size = 20,face = "bold"),
axis.text.x = element_text(colour = "black",size = 14),
axis.text.y = element_text(colour = "black",size = 14),
legend.text = element_text(size = 15,face = "bold"),
legend.position = "none"#是否删除图例

)
p
FileName <- paste(name_i,"aov_nosig"," bar", ".pdf", sep = "_")
ggsave(FileName, p, width = 8, height = 8)
}

}else if (p1 <.05| p2 <.05){
p1 <- round(p1,3)
p2 <- round(p2,3)
data_i = data_wt[i]
ee    <-as.matrix(data_i)
dd <- as.vector(ee)
name_i = colnames(data_wt[i])
xax = data.frame(dd = dd,group = data_wt$group) krusk=compare_means(dd ~ group, data=xax, method = "kruskal.test") sumkrusk=as.data.frame(krusk) if ( sumkrusk[3]< 0.05) { out <- LSD.test(model,"group", p.adj="none")#进行多重比较，不矫正P值 aa = out$group#结果显示：标记字母法
aa$group = row.names(aa) out <- LSD.test(model,"group", p.adj="none")#进行多重比较，不矫正P值 aa = out$group#结果显示：标记字母法
aa$group = row.names(aa) aa wen1 = as.data.frame(tapply(dd,data_wt$group,mean,na.rm=TRUE))
wen2 = as.data.frame(tapply(dd,data_wt$group,sd,na.rm=TRUE)) went = cbind(wen1,wen2) wentao = merge(aa,went, by="row.names",all=F) colnames(wentao) = c(colnames(wentao[1:4]),"mean" ,"SD") aa = mutate(wentao, ymin = mean - SD, ymax = mean + SD) a = max(aa$mean)*1.2
# ss <- round(wtx3$Pr(>F)[1],3) mi=c("#1B9E77" ,"#D95F02", "#7570B3","#E7298A") p = ggplot(aa , aes(x = group, y = mean,colour= group)) + geom_bar(aes(colour= group,fill = group),stat = "identity", width = 0.4,position = "dodge") + # geom_text(aes(label = groups,y=ymax, x = group,vjust = -0.3,size = 6))+ geom_errorbar(aes(ymin=ymin, ymax=ymax), colour="black",width=0.1,size = 1)+ scale_y_continuous(expand = c(0,0),limits = c(0,a))+ labs(x=paste(name_i,"of all group", sep = "_"), y="group", title = paste("Normality test",p1,"Homogeneity of variance",p2,"kruskal.test",sumkrusk[3],sep = ":")) p p=p+theme_bw()+ geom_hline(aes(yintercept=mean(dd)), colour="black", linetype=2) + geom_vline(aes(xintercept=0), colour="black", linetype="dashed") + # scale_fill_manual(values = mi, guide = guide_legend(title = NULL))+ theme( panel.grid.major=element_blank(), panel.grid.minor=element_blank(), plot.title = element_text(vjust = -8.5,hjust = 0.1), axis.title.y =element_text(size = 20,face = "bold",colour = "black"), axis.title.x =element_text(size = 24,face = "bold",colour = "black"), axis.text = element_text(size = 20,face = "bold"), axis.text.x = element_text(colour = "black",size = 14), axis.text.y = element_text(colour = "black",size = 14), legend.text = element_text(size = 15,face = "bold"), legend.position = "none"#是否删除图例 ) p FileName <- paste(name_i,"_kruskal.test_YES_bar", ".pdf", sep = "_") ggsave(FileName, p, width = 8, height = 8) }else if ( sumkrusk[3] >= 0.05) { out <- LSD.test(model,"group", p.adj="none")#进行多重比较，不矫正P值 aa = out$group#结果显示：标记字母法
aa$group = row.names(aa) aa wen1 = as.data.frame(tapply(dd,data_wt$group,mean,na.rm=TRUE))
wen2 = as.data.frame(tapply(dd,data_wt$group,sd,na.rm=TRUE)) went = cbind(wen1,wen2) wentao = merge(aa,went, by="row.names",all=F) colnames(wentao) = c(colnames(wentao[1:4]),"mean" ,"SD") aa = mutate(wentao, ymin = mean - SD, ymax = mean + SD) a = max(aa$mean)*1.2
# ss <- round(wtx3$Pr(>F)[1],3) mi=c("#1B9E77" ,"#D95F02", "#7570B3","#E7298A") p = ggplot(aa , aes(x = group, y = mean,colour= group)) + geom_bar(aes(colour= group,fill = group),stat = "identity", width = 0.4,position = "dodge") + # geom_text(aes(label = groups,y=ymax, x = group,vjust = -0.3,size = 6))+ geom_errorbar(aes(ymin=ymin, ymax=ymax), colour="black",width=0.1,size = 1)+ scale_y_continuous(expand = c(0,0),limits = c(0,a))+ labs(x=paste(name_i,"of all group", sep = "_"), y="group", title = paste("Normality test",p1,"Homogeneity of variance",p2,"kruskal.test",sumkrusk[3],sep = ":")) p p=p+theme_bw()+ geom_hline(aes(yintercept=mean(mean)), colour="black", linetype=2) + geom_vline(aes(xintercept=0), colour="black", linetype="dashed") + # scale_fill_manual(values = mi, guide = guide_legend(title = NULL))+ theme( panel.grid.major=element_blank(), panel.grid.minor=element_blank(), plot.title = element_text(vjust = -8.5,hjust = 0.1), axis.title.y =element_text(size = 20,face = "bold",colour = "black"), axis.title.x =element_text(size = 24,face = "bold",colour = "black"), axis.text = element_text(size = 20,face = "bold"), axis.text.x = element_text(colour = "black",size = 14), axis.text.y = element_text(colour = "black",size = 14), legend.text = element_text(size = 15,face = "bold"), legend.position = "none"#是否删除图例 ) p FileName <- paste(name_i,"_kruskal.test_nosig_bar", ".pdf", sep = "") ggsave(FileName, p, width = 8, height = 8) } } }else if( plot == “box”) { if (p1 >.05& p2 >.05) { p1 <- round(p1,3) p2 <- round(p2,3) data_i = data_wt[i] ee <-as.matrix(data_i) dd <- as.vector(ee) name_i = colnames(data_wt[i]) model<-aov(dd ~ group, data=data_wt)#方差分析 wtx1 = summary(model) wtx2 = wtx1[[1]] wtx3 = wtx2[5] if ( wtx3$Pr(>F)[1]< 0.05) {
out <- LSD.test(model,"group", p.adj="none")#进行多重比较，不矫正P值
aa = out$group#结果显示：标记字母法 aa$group = row.names(aa)
a = max(aa$dd)*1.2 data_box = data_wt[c(1,2,i)] colnames(data_box) = c("ID" , "group","dd" ) out = LSD.test(model,"group", p.adj="none") # alternative fdr stat = out$groups
data_box$stat=stat[as.character(data_box$group),]$groups max=max(data_box[,c("dd")]) min=min(data_box[,c("dd")]) x = data_box[,c("group","dd")] y = x %>% group_by(group) %>% summarise_(Max=paste('max(',"dd",')',sep="")) y=as.data.frame(y) rownames(y)=y$group
data_box$y=y[as.character(data_box$group),]$Max + (max-min)*0.05 # mi=c("#1B9E77" ,"#D95F02", "#7570B3","#E7298A") p = ggplot(data_box, aes(x=group, y=data_box[["dd"]], color=group)) + geom_boxplot(alpha=1, outlier.size=0, size=0.7, width=0.5, fill="transparent") + labs(x=paste(name_i," group", sep = "_"), y="group", title = paste("Normality test",p1,"Homogeneity of variance",p2,sep = ":"))+ geom_text(data=data_box, aes(x=group, y=y, color=group, label= stat)) + geom_jitter( position=position_jitter(0.17), size=1, alpha=0.7)+theme(legend.position="none") p p=p+theme_bw()+ geom_hline(aes(yintercept=mean(dd)), colour="black", linetype=2) + geom_vline(aes(xintercept=0), colour="black", linetype="dashed") + # scale_fill_manual(values = mi, guide = guide_legend(title = NULL))+ theme( panel.grid.major=element_blank(), panel.grid.minor=element_blank(), plot.title = element_text(vjust = -8.5,hjust = 0.1), axis.title.y =element_text(size = 20,face = "bold",colour = "black"), axis.title.x =element_text(size = 24,face = "bold",colour = "black"), axis.text = element_text(size = 20,face = "bold"), axis.text.x = element_text(colour = "black",size = 14), axis.text.y = element_text(colour = "black",size = 14), legend.text = element_text(size = 15,face = "bold"), legend.position = "none"#是否删除图例 ) p if (length(unique(data_box$group))>3){    p=p+theme(axis.text.x=element_text(angle=45,vjust=1, hjust=1))}
FileName <- paste(name_i,"_aov_LSD_box", ".pdf", sep = "")
ggsave(FileName, p, width = 8, height = 8)
}else if ( wtx3$Pr(>F)[1]>= 0.05) { out <- LSD.test(model,"group", p.adj="none")#进行多重比较，不矫正P值 aa = out$group#结果显示：标记字母法
aa$group = row.names(aa) a = max(aa$dd)*1.2

data_box = data_wt[c(1,2,i)]
colnames(data_box) = c("ID" , "group","dd" )

out = LSD.test(model,"group", p.adj="none") # alternative fdr
stat = out$groups data_box$stat=stat[as.character(data_box$group),]$groups
max=max(data_box[,c("dd")])
min=min(data_box[,c("dd")])
x = data_box[,c("group","dd")]
y = x %>% group_by(group) %>% summarise_(Max=paste('max(',"dd",')',sep=""))
y=as.data.frame(y)
rownames(y)=y$group data_box$y=y[as.character(data_box$group),]$Max + (max-min)*0.05

ss <- round(wtx3$Pr(>F)[1],3) # mi=c("#1B9E77" ,"#D95F02", "#7570B3","#E7298A") p = ggplot(data_box, aes(x=group, y=data_box[["dd"]], color=group)) + geom_boxplot(alpha=1, outlier.size=0, size=0.7, width=0.5, fill="transparent") + labs(x=paste(name_i,"box", sep = "_"), y="group", title = paste("Normality test",p1,"Homogeneity of variance",p2,"aov",ss,sep = ":"))+ # geom_text(data=data_box, aes(x=group, y=y, color=group, label= stat)) + geom_jitter( position=position_jitter(0.17), size=1, alpha=0.7)+theme(legend.position="none") p p=p+theme_bw()+ geom_hline(aes(yintercept=mean(dd)), colour="black", linetype=2) + geom_vline(aes(xintercept=0), colour="black", linetype="dashed") + # scale_fill_manual(values = mi, guide = guide_legend(title = NULL))+ theme( panel.grid.major=element_blank(), panel.grid.minor=element_blank(), plot.title = element_text(vjust = -8.5,hjust = 0.1), axis.title.y =element_text(size = 20,face = "bold",colour = "black"), axis.title.x =element_text(size = 24,face = "bold",colour = "black"), axis.text = element_text(size = 20,face = "bold"), axis.text.x = element_text(colour = "black",size = 14), axis.text.y = element_text(colour = "black",size = 14), legend.text = element_text(size = 15,face = "bold"), legend.position = "none"#是否删除图例 ) p if (length(unique(data_box$group))>3){    p=p+theme(axis.text.x=element_text(angle=45,vjust=1, hjust=1))}
FileName <- paste(name_i,"_aov_nosig_box", ".pdf", sep = "")
ggsave(FileName, p, width = 8, height = 8)
}

}else if (p1 <.05| p2 <.05){
p1 <- round(p1,3)
p2 <- round(p2,3)
data_i = data_wt[i]
ee    <-as.matrix(data_i)
dd <- as.vector(ee)
name_i = colnames(data_wt[i])
xax = data.frame(dd = dd,group = data_wt$group) krusk=compare_means(dd ~ group, data=xax, method = "kruskal.test") sumkrusk=as.data.frame(krusk) if ( sumkrusk[3]< 0.05) { out <- LSD.test(model,"group", p.adj="none")#进行多重比较，不矫正P值 aa = out$group#结果显示：标记字母法
aa$group = row.names(aa) a = max(aa$dd)*1.2

data_box = data_wt[c(1,2,i)]
colnames(data_box) = c("ID" , "group","dd" )

out = LSD.test(model,"group", p.adj="none") # alternative fdr
stat = out$groups data_box$stat=stat[as.character(data_box$group),]$groups
max=max(data_box[,c("dd")])
min=min(data_box[,c("dd")])
x = data_box[,c("group","dd")]
y = x %>% group_by(group) %>% summarise_(Max=paste('max(',"dd",')',sep=""))
y=as.data.frame(y)
rownames(y)=y$group data_box$y=y[as.character(data_box$group),]$Max + (max-min)*0.05

# mi=c("#1B9E77" ,"#D95F02", "#7570B3","#E7298A")
wtq = levels(data_wt$group) lis = combn(levels(data_wt$group), 2)
x <-lis
my_comparisons <- tapply(x,rep(1:ncol(x),each=nrow(x)),function(i)i)

p = ggplot(data_box, aes(x=group, y=data_box[["dd"]], color=group)) +
geom_boxplot(alpha=1, outlier.size=0, size=0.7, width=0.5, fill="transparent") +
labs(x=paste(name_i,"of all group", sep = "_"),
y="group",
title = paste("Normality test",p1,"Homogeneity of variance",p2,sep = ":"))+
# geom_text(data=data_box, aes(x=group, y=y, color=group, label= stat)) +
geom_jitter( position=position_jitter(0.17), size=1, alpha=0.7)+theme(legend.position="none")+
stat_compare_means()+
stat_compare_means(comparisons=my_comparisons,label = "p.signif",hide.ns = F) # Add pairwise

p
p=p+theme_bw()+
geom_hline(aes(yintercept=mean(dd)), colour="black", linetype=2) +
geom_vline(aes(xintercept=0), colour="black", linetype="dashed") +
# scale_fill_manual(values = mi, guide = guide_legend(title = NULL))+
theme(

panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),

plot.title = element_text(vjust = -8.5,hjust = 0.1),
axis.title.y =element_text(size = 20,face = "bold",colour = "black"),
axis.title.x =element_text(size = 24,face = "bold",colour = "black"),
axis.text = element_text(size = 20,face = "bold"),
axis.text.x = element_text(colour = "black",size = 14),
axis.text.y = element_text(colour = "black",size = 14),
legend.text = element_text(size = 15,face = "bold"),
legend.position = "none"#是否删除图例

)
p
if (length(unique(data_box$group))>3){ p=p+theme(axis.text.x=element_text(angle=45,vjust=1, hjust=1))} FileName <- paste(name_i,"_kruskal.test_wlc_box_", ".pdf", sep = "") ggsave(FileName, p, width = 8, height = 8) }else if ( sumkrusk[3] >= 0.05) { out <- LSD.test(model,"group", p.adj="none")#进行多重比较，不矫正P值 aa = out$group#结果显示：标记字母法
aa$group = row.names(aa) a = max(aa$dd)*1.2

data_box = data_wt[c(1,2,i)]
colnames(data_box) = c("ID" , "group","dd" )

out = LSD.test(model,"group", p.adj="none") # alternative fdr
stat = out$groups data_box$stat=stat[as.character(data_box$group),]$groups
max=max(data_box[,c("dd")])
min=min(data_box[,c("dd")])
x = data_box[,c("group","dd")]
y = x %>% group_by(group) %>% summarise_(Max=paste('max(',"dd",')',sep=""))
y=as.data.frame(y)
rownames(y)=y$group data_box$y=y[as.character(data_box$group),]$Max + (max-min)*0.05

ss <- round(wtx3$Pr(>F)[1],3) # mi=c("#1B9E77" ,"#D95F02", "#7570B3","#E7298A") p = ggplot(data_box, aes(x=group, y=data_box[["dd"]], color=group)) + geom_boxplot(alpha=1, outlier.size=0, size=0.7, width=0.5, fill="transparent") + labs(x=paste(name_i,"box", sep = "_"), y="group", title = paste("Normality test",p1,"Homogeneity of variance",p2,"aov",ss,sep = ":"))+ # geom_text(data=data_box, aes(x=group, y=y, color=group, label= stat)) + geom_jitter( position=position_jitter(0.17), size=1, alpha=0.7)+theme(legend.position="none")+ stat_compare_means() p p=p+theme_bw()+ geom_hline(aes(yintercept=mean(dd)), colour="black", linetype=2) + geom_vline(aes(xintercept=0), colour="black", linetype="dashed") + # scale_fill_manual(values = mi, guide = guide_legend(title = NULL))+ theme( panel.grid.major=element_blank(), panel.grid.minor=element_blank(), plot.title = element_text(vjust = -8.5,hjust = 0.1), axis.title.y =element_text(size = 20,face = "bold",colour = "black"), axis.title.x =element_text(size = 24,face = "bold",colour = "black"), axis.text = element_text(size = 20,face = "bold"), axis.text.x = element_text(colour = "black",size = 14), axis.text.y = element_text(colour = "black",size = 14), legend.text = element_text(size = 15,face = "bold"), legend.position = "none"#是否删除图例 ) p if (length(unique(data_box$group))>3){    p=p+theme(axis.text.x=element_text(angle=45,vjust=1, hjust=1))}
FileName <- paste(name_i,"_kruskal.test_nosig_box", ".pdf", sep = "")
ggsave(FileName, p, width = 8, height = 8)
}

}

}
}

### 结果文件展示
![fig10](http://bailab.genetics.ac.cn/Note/WenTao/190521R_DA/fig10.png)

### 选择Tukey多重比较方法

Tukey多重计较方法，用来调用字母表示多重比较结果不同，下面是我整理好的代码：
{r Tukey, echo=TRUE, message=FALSE, warning=FALSE, paged.print=FALSE}
library(multcomp)

model<-aov(dd ~ group, data=data_wt)#方差分析
# model<-aov(total.ASA.mg.g.1FW. ~ gruop, data=data_wt)#方差分析
summary(model)

litter.mc <- glht(model, linfct = mcp(group = 'Tukey'))
summary(litter.mc)

insx = cld(litter.mc)
TUK_a <- insx$mcletters$monospacedLetters
TUK_a = as.data.frame(TUK_a)
colnames(TUK_a) = c("value_aov")
head(TUK_a)

写在后面

https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA

http://blog.sciencenet.cn/blog-3334560-1184225.html

全部精选博文导读

GMT+8, 2019-7-18 05:21

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社