R语言一键批量完成差异统计和可视化
撰文：文涛 南京农大
责编：刘永鑫 中科院遗传发育所
R语言一键完成差异检测从数据到展示
单因素差异分析的完整方案
关键词：正态性检验；方差齐性；非参数检验；秩和检验；多重比较；带显著性字母柱状图或箱线图
由于作者水平有限，大家可以添加我的个人微信讨论细节、bug和可改进的地方（微信号：nanjingxuezi）
方案优点
这份方案有一下几个特征，在展示之前给大家mark一下
- 可视化图形标注方差齐性结果，正态性检验结果等
- 自动识别数据指标，自动保存分析结果，文件名称标记分析指标及其分析操作内容和方式。
引子
记得从2016年入学以来，老板叫做的第一个分析便是单因素方差分析,对R 来讲也就是一个aov函数。单因素方差分析确实在很多情况下确实是大多数人的需求，基于R语言的实现也很简单。但是做完之后我便是被吐槽，没有正态性检验，没有方差齐性检验，于是之后的一天，就做了一个简单流程，当时发布在我的个人公众号：微生信生物：《R语言绘制带有显著性字母标记的柱状图》；内容是：首先QQ图，方差齐性检验，后又做了aov和多重比较。并写了个简单的循环。大家可以看到明显不够完整。今天我来的目的就是完善单因素方差分析并且在不适合方差分析的情况下的非参数检验也加入方案中，其次可视化也做了一个完善。（大家有想法的留言，我将加在本框架内，完善差异检测方案）
那么为什么选在这样一个日子来完善这样一份代码呢？刘老师NBT上线不久，其中fig6让我很是怀念之前使用R语言出带显著性标记的柱状图。于是才有了今天可视化方案中的planB。当然刘老师亲自做的NBT的分析值得我逐行运行，尤其是其中的物种分类的GraPhlan，算是到目前为止，我见过的最为漂亮的物种分类树。相信刘老师在之后的文章解读部分会为大家详细解读。
我结合之前的工作同时结合刘老师的工作，设计了一个这样的思路：试图方便自己的工作，同时希望帮到大家。
单因素差异检测完整方案实现思路
试验中经常测定的指标共同的特征是：不同的处理，有重复，需要做差异检测。因此这里我首先对数据进行正态性检验和方差齐性检验，判断符合后进行多重比较并选择喜欢的可视化方案（这里我提供了两种可视化方法，分别是：柱状图，箱线图），判断不符合后，进行非参数检验，首先进行kruskal.test检验（同时对多组进行差异检验），如果有差异，我将继续进行Wilcoxon秩和检验，之后便选择两种可视化方案中的一种进行可视化。
我写了一个骨架：
主要函数解读
两种差异表示方案及其代码
下面是进行LSD多重比较及其添加表征差异的字母柱状图代码：
#值得注意的是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
字母标记箱线图代码
这部分代码来源自刘老师NBT中的fig6，也正是这份代码，让我产生了重新升级差异分析的念头。说句题外话，刘老师的代码书写思路我总有一种亲切的感觉。
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 + 箱线图 + 连线差异标注
由于两组之间的连线需要指定两组信息，这里我又想将所有组之间的差异展示出来，所以使用combn函数得到分组信息两两匹配的结果，并使用tapply结合函数将矩阵改变为列表。完成这一工作。
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
实战
测序数据和Rmd代码，生信宝典后台回复“anova”获取。
写在后面
值得注意的是，我们在方差分析之后可以选择t检验两两比对差异，并使用P值矫正。同样非参数检验也使用类似的方法进行两两比对，但是就两组之间的显著性结果转化为字母标记，我并没有实现，大家如果有想法，多多赐教。