差异分析完整解决方案升级版（更正版） / 四六文摘

写在前面

上次的推送差异分析完整解决方案升级版存在一个致命错误，这里我先给大家道歉，由于我的疏忽正态检验过程我没有进行分组正态检验，这将造成大量的数据原本正态分布，却误判。今天我来给大家带来忏悔版本的差异分析完整解决方案。

我本来想自己写一个多组正态检验的函数的，但是出于谨慎，我在网上搜了一下，得到了shapiro.test.multi函数。就不用重复造轮子了。我将输出做一个修改，同修改差异分析函数内容。这里我放上作者的链接https://blog.csdn.net/yanlingbin/article/details/81194134。有兴趣的朋友点击查看

文末提供这个原始数据和代码及其测试脚本。

shapiro.test.multi <- function( #定义函数名 data, #定义函数第一个参数 value, #第2参数 group) #第3参数 { #开始计算


require(magrittr) #按需要加载管道函数包
table(data[,group]) %>% #提取分组信息，此处即为统计group中值出现的次数，达到了去重的目的

data.frame(.) -> a1 #将提取信息从table格式转为数据库data.frame并存为a1，这样才能提取其中一列转为向量
a2 <- as.vector(a1[,1]) #将a1数据的第一列转为向量，这一列即为不重复的分组信息
data = data.frame(group = data[,group], #对数据集进行关键变量的提取，提取分组变量为data数据集的group变量

value = data[,value]) #提取计算值为data数据集的value
test.result <- data.frame(No=0, #行号

Name=0, #分组名

W=0, #W值

p.value=0, #p值

norm.test=0) #检测结果
for (i in (1:length(a2))){ #定义for循环计算，从1到a2的长度值这样一个区间，a2的长度即为分组的数量

subset(data, #指定取数据集 换行写代码使层次更清晰

group == a2[i], #定义去子集的条件，“==”为判断

select = value) %>% #定义需要取集的变量/列，“=”为定义

.[,1] %>% # "."定义计算结果放置的位置

shapiro.test(.) -> t.r #进行正态检验分布并存储为t.r
test.result[i,1] = i #存储组序号

test.result[i,2] = a2[i] #存储分组名

test.result[i,3] = t.r$statistic #存储W统计量

test.result[i,4] = t.r$p.value #存储计算的p值
if #if判断

(t.r$p.value > 0.05) #判断条件

test.result[i,5] = "Norm" #通过判断后的命令

else

test.result[i,5] = "Other_situation" #未通过判断执行的命令

} #结束循环计算
test.result[nrow( test.result)+1,1] = "Test Method:" #给数据框加上检验正态分布方法信息，在最后一行之后加上一行，在第1列放入次数据

test.result[nrow( test.result),2] = "Shapiro-Wilk" #同上行，存在第二列

test.result #显示用于存储计算结果的数据框 } #脚本结束

一行代码就跑过去了，这里我提取了判断指标，也就是全部分组都满足正态分布才被判断为真。

length(a$p.value[-dim(a)[1]][!a$p.value[-dim(a)[1]] >= 0.05]) == 0

我看见过最全的差异分析教程，但是好像还没有最智能的，机械的差异分析教程。因为这些标准都是死的，我们完全可以使用 if 来解放我们的大脑。

之前我有一篇推送：R语言一键批量完成差异统计和可视化提供了这个差异检测编码流程，但是做的比较粗糙，代码冗余过多，出现错误向上排查比较困难。所以这次我将本来1000行的代码修改到了300行，精简了许多。

一张神图，解决科研统计80%的问题。如果如作者所言：神图，解决科研统计80%的问题，那么我的这个流程把这张图90%的内容都进行了自动化数据识别和自动分析。顺便画张图，还有两种出图方式任你选择（建议箱线图，更准确）。

同时有几个地方我进行了升级：

aov添加多种多重比较的方法，并映射到图形上。
提供两种作图方式，两种差异展示方式。
提供表格形式的结果文件方便查看。

这次代码升级有几个地方需要认真编写：

1. 柱状图的差异分析line形式的差异自动添加：如图：

这里我还是本来想使用ggpubr进行做的，发现这里面根本就没有柱状图添加line形式的差异的方法，最后还是使用了ggsignif包，其实就这一个包就够了。

if (sig_show == "line") { zuhe = combn(aa$group,2) xxxx <- tapply(zuhe,rep(1:ncol(zuhe),each=nrow(zuhe)),function(i)i) xxxx sig_lis = rep("a",dim(zuhe)[2]) for (i in 1:dim(zuhe)[2]) { library(tidyverse) library("ggsignif")


if (filter(aa, group == xxxx[[i]][1])$groups == filter(aa, group == xxxx[[i]][2])$groups) {

sig_lis[i] = "no_sig"

}
if (filter(aa, group == xxxx[[i]][1])$groups != filter(aa, group == xxxx[[i]][2])$groups) {

sig_lis[i] = "*"

}
}

2. 两两非参数检验结果是成对的，之前我可能提高多如何将成对差异结果表示为字母，这是本次升级的又一个技巧：multcompLetters函数。

names(wilcox_levels) = paste(xx$group1,xx$group2,sep = "-")

wilcox.labels <- data.frame(multcompLetters(wilcox_levels, threshold = 0.05)['Letters']) colnames(wilcox.labels) = "groups" aa = wilcox.labels

3. 那如果将字母形式的显著性标记结果转化为两两比较的结果呢，这里有我自己编写了一个脚本，做了这个转化；

zuhe = combn(aa$group,2) xxxx <- tapply(zuhe,rep(1:ncol(zuhe),each=nrow(zuhe)),function(i)i) xxxx sig_lis = rep("a",dim(zuhe)[2]) for (i in 1:dim(zuhe)[2]) { library(tidyverse) library("ggsignif")


if (filter(aa, group == xxxx[[i]][1])$groups == filter(aa, group == xxxx[[i]][2])$groups) {

sig_lis[i] = "no_sig"

}
if (filter(aa, group == xxxx[[i]][1])$groups != filter(aa, group == xxxx[[i]][2])$groups) {

sig_lis[i] = "*"

}

}

到这里我们似乎就没什么难度了，将整个函数编写完毕，运行不同的参数进行调试即可：

01. 导入数据

# 读入实验设计 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)

02. 下面运行函数，首先测试五种多重比较方法

single_diff_and_plot (data_wt,plotname = "南京农业大学_",plot = "bar",method_Mc ="LSD") single_diff_and_plot (data_wt,plotname = "南京农业大学_",plot = "bar",method_Mc ="scheffe") single_diff_and_plot (data_wt,plotname = "资源与环境科学学院_",plot = "box",method_Mc ="LSD" ) single_diff_and_plot (data_wt,plotname = "资源与环境科学学院_",plot = "box",method_Mc = "scheffe" ) method_Mc == "Duncan" single_diff_and_plot (data_wt,plotname = "资源与环境科学学院_",plot = "box",method_Mc = "Duncan" )

03. 接着测试显著性标注方式

single_diff_and_plot (data_wt,plotname = "资源与环境科学学院_",plot = "box",method_Mc = "Duncan",sig_show = "abc" )

single_diff_and_plot (data_wt,plotname = "资源与环境科学学院_",plot = "box",method_Mc = "Duncan",sig_show = "line" ) single_diff_and_plot (data_wt,plotname = "资源与环境科学学院_",plot = "box",method_Mc = "Duncan",sig_show = "abc" )

函数默认参数如下：

data_wt = data_wt,##默认数据集 plotname = "wentao_",#出图同意带有这个前缀 plot = "bar",#出图形式为bar

#method_cv: 方差齐心方法选择,有两种方法：method_cv == "leveneTest",method_cv == "bartlett.test" method_cv = "bartlett.test" , method_Mc = "Tukey", sig_show = "abc"

04. 成图和差异表格列表欣赏

更正版函数

https://github.com/taowenmicro/easyMicrobiome

差异分析完整解决方案升级版（更正版）

相关推荐