淺析單因素方差分析中的多重比較

本腳本側(cè)重于單因素方差分析中多重比較方法的運(yùn)用;
就不展示數(shù)據(jù)正態(tài)性及齊次性的運(yùn)算了(默認(rèn)都符合，一般理化數(shù)據(jù)是都符合的);
有的人喜歡用Tukey檢驗(yàn)，但會(huì)遇到一些不符合預(yù)期的問(wèn)題;
讓我們抽絲剝繭的來(lái)理清這些個(gè)問(wèn)題,尤其注重閱讀下面的討論說(shuō)明(說(shuō)不定你就遇到過(guò)這樣的問(wèn)題)；
這里用的數(shù)據(jù)涉及到多個(gè)α多樣性，多個(gè)的處理(若你是做基因你可以理解成多個(gè)采樣地的多個(gè)基因)同時(shí)進(jìn)行多重比較。

代碼和測(cè)試數(shù)據(jù)下載：http://210.75.224.110/github/Note/R/multcomp.zip

library(ggplot2)
library(ggprism)
dat <- read.table('./alpha.txt',row.names = 1,header = T,stringsAsFactors = F)#讀入α多樣性數(shù)據(jù)
head(dat, n = 3)
design <- read.table('./metadata.tsv',row.names = 1,header = T,stringsAsFactors = F)#讀入試驗(yàn)設(shè)計(jì)文件
head(design, n = 3)
dat <- merge(dat,design,by='row.names')#按照行名合并文件
head(dat, n = 3)
library(reshape2)
dat <- melt(dat,id.vars = -c(2:7),variable.name = 'alpha')#寬數(shù)據(jù)變長(zhǎng)數(shù)據(jù)
head(dat, n = 3)
dat$alpha <- as.factor(dat$alpha)#將α列轉(zhuǎn)化成因子
names(dat)[4] <- 'v'#給value重新賦列名
head(dat, n = 3)

函數(shù)和參數(shù)簡(jiǎn)介

函數(shù)參數(shù)設(shè)置：

• data就是上面整好的數(shù)據(jù),

• group是你的分組信息列，比如α多樣性的種類(lèi)(或不同的基因),

• compare是每個(gè)α多樣性要比較的不同處理(或每個(gè)gene要比較的不同處理),

• value 值就是要比較的α多樣性/gene拷貝數(shù)的數(shù)值。

整體思想如下(例如本數(shù)據(jù))：
首先給輸入數(shù)據(jù)dat，根據(jù)alpha列分成不同的小子集，每個(gè)小子集比較不同Group下v值的差異情況，最后匯總輸出。

# 1 -----------------------------------------------------------------------
ONE_Tukey_HSD1 <- function(data,group,compare,value){

library(multcomp)#Tukey檢驗(yàn)需要用到這個(gè)包來(lái)標(biāo)顯著性字母標(biāo)記

a <- data.frame(stringsAsFactors = F)#做一個(gè)空的數(shù)據(jù)框
type <- unique(data[,group])#統(tǒng)計(jì)需要運(yùn)行多重比較的次數(shù)
for (i in type)#進(jìn)行type次多重比較
{
g1=compare
sub_dat <- data[data[,group]==i,]#根據(jù)指定的i去取相應(yīng)的數(shù)據(jù)集出來(lái)
#fit <- aov(sub_dat[,value] ~ sub_dat[,compare] )
names(sub_dat)[names(sub_dat)==compare] <- 'g1' ## 重命名方便后面使用
names(sub_dat)[names(sub_dat)==value] <- 'value' ## 重命名方便后面使用
sub_dat$g1 <- factor(sub_dat$g1)#將列轉(zhuǎn)化成因子以進(jìn)行多重比較

fit <- aov(value ~ g1,data = sub_dat )#方差分析
#Tukey_HSD = TukeyHSD(fit, ordered = TRUE, conf.level = 0.95)
options(warn = -1)
tuk <- cld(glht(fit, alternative = 'two.sided', linfct = mcp(g1 = 'Tukey')), decreasing = TRUE)#Tukey檢驗(yàn)多重比較
Tukey.labels <- data.frame(Letters=tuk$mcletters$Letters, stringsAsFactors = FALSE)#獲取多重比較字母標(biāo)注
Tukey.labels$compare = rownames(Tukey.labels)## 提取字母分組行名為group組名
Tukey.labels$type <- i

mean_sd <- merge(aggregate(sub_dat[['value']],by=list(sub_dat[,'g1']),FUN=sd),#獲取數(shù)據(jù)標(biāo)準(zhǔn)差
aggregate(sub_dat[['value']],by=list(sub_dat[,'g1']),FUN=mean),by="Group.1"#獲取數(shù)據(jù)均值
)
names(mean_sd) <- c('compare','std','mean')#列名重命名

a <- rbind(a,merge(mean_sd,Tukey.labels,by='compare'))#合并數(shù)據(jù)
}

names(a) <- c(compare,'std','mean','Letters',group)#列名重命名
return(a)
}

# 2 -----------------------------------------------------------------------

ONE_Tukey_HSD2 <- function(data,group,compare,value){
library(multcompView)

a <- data.frame(stringsAsFactors = F)
type <- unique(data[,group])
for (i in type)
{
g1=compare
sub_dat <- data[data[,group]==i,]
#fit <- aov(sub_dat[,value] ~ sub_dat[,compare] )
## 重命名方便后面使用
names(sub_dat)[names(sub_dat)==compare] <- 'g1'
names(sub_dat)[names(sub_dat)==value] <- 'value'
sub_dat$g1 <- factor(sub_dat$g1)

fit <- aov(value ~ g1,data = sub_dat )
Tukey_HSD = TukeyHSD(fit, ordered = TRUE, conf.level = 0.95)
options(warn = -1)
tuk <- multcompLetters2(value ~ g1, Tukey_HSD$g1[,"p adj"], sub_dat)


#tuk <- cld(glht(fit, alternative = 'two.sided', linfct = mcp(g1 = 'Tukey')), decreasing = TRUE)
Tukey.labels <- data.frame(tuk['Letters'], stringsAsFactors = FALSE)
## 提取字母分組行名為group組名
Tukey.labels$compare = rownames(Tukey.labels)
Tukey.labels$type <- i

mean_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,Tukey.labels,by='compare'))
}

names(a) <- c(compare,'std','mean','Letters',group)
return(a)
}

ONE_Tukey_HSD1函數(shù)

這個(gè)函數(shù)核心
是cld(glht(fit, alternative = 'two.sided’, linfct = mcp(g1 = 'Tukey’)), decreasing = TRUE),
不會(huì)出現(xiàn)c>b>a情況(因?yàn)閐ecreasing = TRUE,當(dāng)然有的人喜歡這樣標(biāo),)和亂標(biāo)字母(比如對(duì)于mean最低的點(diǎn)
并不一定標(biāo)記成c(a>b>c時(shí))或并不一定標(biāo)記成a(c>b>a時(shí)),其只能保證有差異的數(shù)據(jù)一定是不同字母),
但是多重比較出現(xiàn)“ac”標(biāo)注，沒(méi)法解決。

ONE_Tukey_HSD2函數(shù)

而ONE_Tukey_HSD2核心是這個(gè) multcompLetters2(value ~ g1, Tukey_HSD$g1[,”p adj”], sub_dat),
multcompLetters2這個(gè)函數(shù)隸屬于multcompView包，與 multcompLetters不同的是
multcompLetters2可以接受formula，而multcompLetters只接受一個(gè)兩兩比較的p值的數(shù)據(jù)框,
且可能多重比較時(shí)出現(xiàn)“ac”標(biāo)注,以及出現(xiàn)c>b>a情況和亂標(biāo)字母(比如對(duì)于mean最低的點(diǎn)
并不一定標(biāo)記成c(a>b>c時(shí))或并不一定標(biāo)記成a(c>b>a時(shí)),其只能保證有差異的數(shù)據(jù)一定是不同字母)。

當(dāng)然多重比較好多方法，不要局限于一種方法,
例如下面的第三種可以用library(agricolae)包中的LSD檢驗(yàn)(用的“BH”校正),

當(dāng)然也可以用library(agricolae)包中的
【Duncan法】(新復(fù)極差法)(SSR);
【SNK法】(Student-Newman-Keuls);
【Scheffe檢驗(yàn)】;

這三種多重比較方法同LSD檢驗(yàn)的用法一樣都可以避免出現(xiàn)上面提到的三種情況即：

1、 a、b、c的順序不會(huì)出現(xiàn)c>b>a;

2、不會(huì)出現(xiàn)亂標(biāo)字母(比如對(duì)于mean最低的點(diǎn)并不一定標(biāo)記成c(a>b>c時(shí))或
并不一定標(biāo)記成a(c>b>a時(shí))，其只能保證有差異的數(shù)據(jù)一定是不同字母);

3、多重比較時(shí)出現(xiàn)“ac”標(biāo)注。

ONE_LSD函數(shù)

# 3 -----------------------------------------------------------------------
ONE_LSD <- function(data,group,compare,value){
library(agricolae)

a <- data.frame(stringsAsFactors = F)
type <- unique(data[,group])
for (i in type)
{
# sub_dat <- subset(data,group == i)
sub_dat <- data[data[,group]==i,]
# fit <- aov(value ~ compare,sub_dat)
fit <- aov(sub_dat[,value] ~ sub_dat[,compare] )
out <- LSD.test(fit,'sub_dat[, compare]',p.adj='BH')#進(jìn)行了p值校正
#out$groups就可獲取多重比較字母列表
out$groups$type <- i
out$groups$compare <- rownames(out$groups)

a <- rbind(a,merge(out$means[,1:2], out$groups,by='sub_dat[, value]'))
}
names(a) <- c('mean','std','lsd',group,compare)
return(a)
}

alpha多樣性在不同處理下的差別

運(yùn)行，這里拿alpha多樣性測(cè)試，看不同alpha多樣性在不同處理下的差別。

#1
#df1 <- ONE_Tukey_HSD1(data=dat,group='alpha',compare='Group',value='v')
df1 <- ONE_Tukey_HSD1(dat,'alpha','Group','v')

在此可以查看各個(gè)α多樣性下不同處理間的多重比較字母標(biāo)注結(jié)果，這也是本腳本的亮點(diǎn)之一

數(shù)據(jù)量很大的情況下，可以直接查看差異情況，不用一個(gè)個(gè)的做出圖再點(diǎn)開(kāi)看，很是方便。

df1

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))

本圖一張即可包含所有數(shù)據(jù)情況，方便查看

p1

#2
#df2 <- ONE_Tukey_HSD2(data=dat,group='alpha',compare='Group',value='v')
df2 <- ONE_Tukey_HSD2(dat,'alpha','Group','v')
df2

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

#3
#df3 <- ONE_LSD(data=dat,group='alpha',compare='Group',value='v')
df3 <- ONE_LSD(dat,'alpha','Group','v')
df3

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=lsd))+
facet_wrap(.~alpha,scales = "free_y")+ labs(x='Group',y='AlphaDiv')+
ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p3

Output figure width and height
Letter紙圖片尺寸為單欄89 mm，雙欄183 mm，頁(yè)面最寬為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")

參考資料

EasyAmplicon/script/alpha_boxplot.R

差異分析、顯著性標(biāo)記及統(tǒng)計(jì)作圖的自動(dòng)實(shí)現(xiàn)R代碼示例

multcompView: Visualizations of Paired Comparisons