浅析R语言多组定量资料非参数检验的多组比较及簇状柱形图显著性字母标记之分面与分组的图形艺术

R语言多组定量资料非参数检验的多组比较

非参数检验的应用

本流程是在刘永鑫老师提供的代码资料指导下完成

测试数据及代码链接：https://pan.baidu.com/s/1s0IX2eVe34jL2DViWCvZfQ
提取码：dunw

链接失效，后台回复”多组比较”，获得更新链接。

先简单说明一下非参数检验

参数检验是在已知总体分布的条件下(一般要求总体服从正态分布)对一些主要的参数(如均值、百分数、方差、相关系数等)进行的检验，有时还要求某些总体参数满足一定的条件。如独立样本的t检验和方差分析不仅要求总体符合正态分布，还要求各总体方差齐性。

而当我们的数据不满足参数检验时我们可以选择非参数检验-秩和检验(Wilcoxon/Kruskal-Wallis H)

非参数检验方法简便，不依赖于总体分布的具体形式，因而适用性强，但灵敏度和精确度不如参数检验。一般而言，非参数检验适用于以下3种情况

(1)顺序类型的数据资料，这类数据的分布形态一般是未知的。

(2)虽然是连续数据，但总体分布形态未知或者非正态，这和卡方检验一样，称为自由分布检验。

(3)总体分布虽然正态，数据也是连续类型，但样本容量极小，如10以下(虽然T检验被
称为小样本统计方法，但样本容量太小时，代表性毕竟很差，最好不要用要求较严格的参数检验法)。

非参数检验的缺点与不足

因为有这些特点，加上非参数检验法一般的原理和计算比较简单，因此常用于一些为正式研究进行探路的预备性研究的数据统计中。当然，由于非参数检验许多牵涉不到参数计算，对数据中的信息利用不够，因而其统计检验力相对参数检验则差得多。

关于使用非参数检验的纠结

在很多人看来，参数检验似乎是统计方法的“主流”，而非参数检验往往被当成“非主流”。大家似乎更喜欢用t检验、方差分析这样的参数检验。即使在数据不满足正态分布的条件下，仍然有使用参数检验的执念；总觉得非参数检验不够正式，或者就是不喜欢用。

而现实中的很多数据都不服从正态分布，这时候非参数检验的统计学效能要高于参数检验。并且，非参数检验更加稳健。参数统计方法是建立在严格假设条件基础上，一旦假设条件不符合，其推断的正确性就不存在了。非参数检验带有最弱的假设，对模型的限制很少，因而天然的具有稳健性。

多组定量资料非参数检验及多组比较的实现

kruskal.test()函数可以确定总体差异是否有统计学意义，但是不知道具体那些组之间存在差异。可以使用wilcox.test()函数每次进行两组数据比较。一种更为有效的方法是在控制犯第一类错误的概率前提下,执行可以同步进行的多组比较,这样可以直接完成所有组之间的成对比较。nparcomp包提供了所需要的非参数多组比较程序。此包中的nparcomp()函数接受的输入为一个两列的数据框,其中一列为因变量,另一列为分组变量。==本推文的第三种方法采用此方法==。

本文分别采用三个统计R包的三种不同函数方法进行非参数多组比较：

• “pgirmess”统计包的kruskalmc()函数

• “PMCMRplus”统计包的posthoc.kruskal.nemenyi.test()函数

• “nparcomp”统计包的nparcomp()函数

自己挑选喜欢的用

library(tidyverse)#数据整理与数据转换用了一些更好用更易懂的函数
library(ggprism)
dat <- read.delim('./vegan.txt')#读入α多样性数据
head(dat, n = 3)
design <- read.delim('./metadata.txt')#读入试验设计文件
head(design, n = 3)
dat <- design %>%dplyr::select(SampleID,Group) %>%inner_join(dat,by='SampleID')#按照SampleID内连接文件
head(dat, n = 3)
dat <- gather(dat,alpha,v,-(SampleID:Group))#宽数据变长数据
head(dat, n = 3)

函数和参数简介

函数参数设置：

• data就是上面整好的数据,

• group是你的分组信息列，比如α多样性的种类(或不同的基因),

• compare是每个α多样性要比较的不同处理(或每个gene要比较的不同处理),

• value 值就是要比较的α多样性/gene拷贝数的数值。

整体思想如下(例如本数据)：
首先给输入数据dat，根据alpha列分成不同的小子集，每个小子集比较不同Group下v值的差异情况，最后汇总输出。

“pgirmess”统计包的kruskalmc()函数

# 1 -----------------------------------------------------------------------
kruskalmc_compare1 <- function(data,group,compare,value){library(multcompView)library(pgirmess)##多组两两比较函数用到的R包library(multcomp)a <- data.frame(stringsAsFactors = F)type <- unique(data[,group])for (i in type){g1=comparesub_dat <- data[data[,group]==i,]names(sub_dat)[names(sub_dat)==compare] <- 'g1'names(sub_dat)[names(sub_dat)==value] <- 'value'sub_dat$g1 <- factor(sub_dat$g1)options(warn = -1)k <- kruskalmc(value ~ g1, data=sub_dat, probs=0.05)dif <- k$dif.com[['difference']]names(dif) <- rownames(k$dif.com)difL <- multcompLetters(dif)label <- data.frame(difL['Letters'], stringsAsFactors = FALSE)label$compare = rownames(label)label$type <- imean_sd <- merge(aggregate(sub_dat[['value']],by=list(sub_dat[,'g1']),FUN=sd),aggregate(sub_dat[['value']],by=list(sub_dat[,'g1']),FUN=mean),by="Group.1")names(mean_sd) <- c('compare','std','mean')a <- rbind(a,merge(mean_sd,label,by='compare'))}names(a) <- c(compare,'std','mean','Letters',group)return(a)
}

“PMCMRplus”统计包的posthoc.kruskal.nemenyi.test()函数

# 2 -----------------------------------------------------------------------
PMCMR_compare1 <- function(data,group,compare,value){library(multcompView)library(multcomp)library(PMCMRplus)library(PMCMR)##多组两两比较函数用到的R包a <- data.frame(stringsAsFactors = F)type <- unique(data[,group])for (i in type){g1=comparesub_dat <- data[data[,group]==i,]names(sub_dat)[names(sub_dat)==compare] <- 'g1'names(sub_dat)[names(sub_dat)==value] <- 'value'sub_dat$g1 <- factor(sub_dat$g1)options(warn = -1)k <- posthoc.kruskal.nemenyi.test(value ~ g1,data=sub_dat)n <- as.data.frame(k$p.value)h <- n %>%mutate(compare=rownames(n)) %>%gather(group,p,-compare,na.rm = TRUE) %>%unite(compare,group,col="G",sep="-")dif <- h$pnames(dif) <- h$GdifL <- multcompLetters(dif)K.labels <- data.frame(difL['Letters'], stringsAsFactors = FALSE)K.labels$compare = rownames(K.labels)K.labels$type <- imean_sd <- merge(aggregate(sub_dat[['value']],by=list(sub_dat[,'g1']),FUN=sd),aggregate(sub_dat[['value']],by=list(sub_dat[,'g1']),FUN=mean),by="Group.1")names(mean_sd) <- c('compare','std','mean')a <- rbind(a,merge(mean_sd,K.labels,by='compare'))}names(a) <- c(compare,'std','mean','Letters',group)return(a)
}

“nparcomp”统计包的nparcomp()函数

# 3 -----------------------------------------------------------------------
nparcomp_compare1 <- function(data,group,compare,value){library(nparcomp)##多组两两比较函数用到的R包library(multcompView)library(do)a <- data.frame(stringsAsFactors = F)type <- unique(data[,group])for (i in type){g1=comparesub_dat <- data[data[,group]==i,]names(sub_dat)[names(sub_dat)==compare] <- 'g1' names(sub_dat)[names(sub_dat)==value] <- 'value' sub_dat$g1 <- factor(sub_dat$g1)k <- nparcomp::nparcomp(value ~ g1,data=sub_dat, alternative = "two.sided")k$Analysisb <- k$Analysisdif <- b$p.Valuenames(dif) <- b$Comparisondifname <- names(dif) %>% Replace(from = ' , ',to='-') %>% Replace(from=c('p\\(','\\)'),to='')#正则表达式替换temp_name <- data.frame(tn=difname) %>% separate(col = 'tn',sep = '-',into = c('n1','n2'))difname <- paste(temp_name$n2,temp_name$n1,sep = '-') %>% Replace(from = ' - ',to='-')#正则表达式替换names(dif) <- difnamedifL <- multcompLetters(dif)labels <- data.frame(difL['Letters'], stringsAsFactors = FALSE)labels$compare = rownames(labels)labels$type <- imean_sd <- merge(aggregate(sub_dat[['value']],by=list(sub_dat[,'g1']),FUN=sd),aggregate(sub_dat[['value']],by=list(sub_dat[,'g1']),FUN=mean),by="Group.1")names(mean_sd) <- c('compare','std','mean')a <- rbind(a,merge(mean_sd,labels,by='compare'))}names(a) <- c(compare,'std','mean','Letters',group)return(a)
}

alpha多样性在不同处理下的差别

运行，这里拿alpha多样性测试，看不同alpha多样性在不同处理下的差别。

#df1 <- kruskalmc_compare1(data=dat,group='alpha',compare='Group',value='v')
df1 <- kruskalmc_compare1(dat,'alpha','Group','v')
df1######字母是反着标的(a<b<c)(测试时的字母顺序，下同，可能在跑你自己的数据时会出现反着标的可能，不过不影响使用，可在AI或PDF编辑器中调)

在此可以查看各个α多样性下不同处理间的多组比较字母标注结果，这也是本脚本的亮点之一

数据量很大的情况下，可以直接查看差异情况，不用一个个的做出图再点开看，很是方便。

p1 = ggplot(dat)+geom_boxplot(aes(x=Group,y=v,fill=Group))+geom_text(data=df1,aes(x=Group,y=mean+1.3*std,label=Letters))+facet_wrap(~alpha,scales = "free_y")+ labs(x='Group',y='AlphaDiv')+ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p1

本图一张即可包含所有数据情况，方便查看

#df2 <- PMCMR_compare1(data=dat,group='alpha',compare='Group',value='v')
df2 <- PMCMR_compare1(dat,'alpha','Group','v')
df2########字母标正着标(a>b>c)

p2 = ggplot(dat)+geom_boxplot(aes(x=Group,y=v,fill=Group))+geom_text(data=df2,aes(x=Group,y=mean+1.3*std,label=Letters))+facet_wrap(~alpha,scales = "free_y")+ labs(x='Group',y='AlphaDiv')+ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p2

#df3 <- nparcomp_compare1(data=dat,group='alpha',compare='Group',value='v')
df3 <- nparcomp_compare1 (dat,'alpha','Group','v')
df3########字母标正着标(a>b>c)

p3 = ggplot(dat)+geom_boxplot(aes(x=Group,y=v,fill=Group))+geom_text(data=df3,aes(x=Group,y=mean+1.3*std,label=Letters))+facet_wrap(~alpha,scales = "free_y")+ labs(x='Group',y='AlphaDiv')+ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p3

簇状柱形图显著性字母标记之分面与分组的图形艺术

只是展示方法(多样性数据可能不适合做这个图，但是纵轴标度一致的数据且有需求要做，适合做这个图)。

颜色设置

library(RColorBrewer)
library("scales")
#自己选颜色
display.brewer.all()

#选取颜色
display.brewer.pal(9,"Set1")

#可以直接选取
yanse=brewer.pal(9,"Set1")
#也可以显示出色号后
show_col(brewer.pal(9,"Set1"))

#再挑选喜欢的色号
yanse=c("#E41A1C","#377EB8","#4DAF4A","#984EA3","#FF7F00","#FFFF33","#A65628","#F781BF","#999999")

分面的图形艺术

数据用的df2

df1$alpha <- factor(df2$alpha,levels = c("richness","chao1","ACE","shannon","simpson","invsimpson"))
H = max(df2$mean)*1.2
#可以根据需求调整分面及分组顺序
p4 = ggplot(df2,aes(x =Group,y = mean,fill = Group))+ facet_wrap(~alpha,nrow=1)+geom_bar(stat = "identity" ,width=1, position = position_dodge(0.9)) + geom_text(aes(label = Letters),size=7,vjust=-1.2,position = position_dodge(0.9))+geom_errorbar(aes(ymin = mean - std/3^0.5,ymax = mean+std/3^0.5),position = position_dodge(0.9), width = .1,color="red",lwd=1)+labs(x='Group',y='AlphaDiv')+scale_y_continuous(expand = c(0,0),limits = c(0,H))+scale_fill_manual(values = yanse, guide = guide_legend(title = NULL))+ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p4

分组的图形艺术

数据用的df2

p5 = ggplot(df2,aes(x =alpha,y = mean,fill = Group))+ geom_bar(stat = "identity", width=0.9, position = position_dodge(0.9)) + geom_text(aes(label = Letters),size=7,vjust=-1.2,position = position_dodge(0.9))+geom_errorbar(aes(ymin = mean - std/3^0.5,ymax = mean+std/3^0.5),position = position_dodge(0.9), width = .1,color="red",lwd=1)+labs(x='Group',y='AlphaDiv')+scale_y_continuous(expand = c(0,0),limits = c(0,H))+scale_fill_manual(values = yanse, guide = guide_legend(title = NULL))+ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p5

Output figure width and height
Letter纸图片尺寸为单栏89 mm，双栏183 mm，页面最宽为247 mm
推荐比例16：10，即半版89 mm x 56 mm; 183 mm x 114 mm

# ggsave("./alpha1.pdf", p1, width = 350, height = 200, units = "mm")
# ggsave("./alpha2.pdf", p2, width = 350, height = 200, units = "mm")
# ggsave("./alpha3.pdf", p3, width = 350, height = 200, units = "mm")
# ggsave("./alpha4.pdf", p4, width = 600, height = 200, units = "mm")
# ggsave("./alpha5.pdf", p5, width = 600, height = 200, units = "mm")

猜你想学

R统计-正态性分布检验[Translation]

R统计-方差齐性检验[Translation]

参考资料

《R语言统计分析与应用》

浅析R语言单因素方差分析中的多重比较

在R语言上做Kruskal–Wallis秩和检验

Wilcoxon, Friedman…扒一扒非参数检验名称中的牛人

作者：flyfly、旭日阳光

责编：马腾飞 南京农业大学

审核：刘永鑫 中科院遗传发育所

猜你喜欢

10000+：菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发Nature Cell专刊 肠道指挥大脑

系列教程：微生物组入门 Biostar 微生物组  宏基因组

专业技能：学术图表 高分文章 生信宝典 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板: Shell  R Perl

生物科普:  肠道细菌 人体上的生命 生命大跃进  细胞暗战 人体奥秘

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读