淺析單因素方差分析中的多重比較本腳本側(cè)重于單因素方差分析中多重比較方法的運用;
就不展示數(shù)據(jù)正態(tài)性及齊次性的運算了(默認都符合，一般理化數(shù)據(jù)是都符合的);
有的人喜歡用Tukey檢驗，但會遇到一些不符合預期的問題;
讓我們抽絲剝繭的來理清這些個問題,尤其注重閱讀下面的討論說明(說不定你就遇到過這樣的問題)；
這里用的數(shù)據(jù)涉及到多個α多樣性，多個的處理(若你是做基因你可以理解成多個采樣地的多個基因)同時進行多重比較。 代碼和測試數(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)#讀入試驗設計文件
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ù)變長數(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ù)簡介函數(shù)參數(shù)設置： data就是上面整好的數(shù)據(jù), group是你的分組信息列，比如α多樣性的種類(或不同的基因), compare是每個α多樣性要比較的不同處理(或每個gene要比較的不同處理), value 值就是要比較的α多樣性/gene拷貝數(shù)的數(shù)值。
整體思想如下(例如本數(shù)據(jù))：
首先給輸入數(shù)據(jù)dat，根據(jù)alpha列分成不同的小子集，每個小子集比較不同Group下v值的差異情況，最后匯總輸出。 # 1 -----------------------------------------------------------------------
ONE_Tukey_HSD1 <- function(data,group,compare,value){
library(multcomp)#Tukey檢驗需要用到這個包來標顯著性字母標記
a <- data.frame(stringsAsFactors = F)#做一個空的數(shù)據(jù)框
type <- unique(data[,group])#統(tǒng)計需要運行多重比較的次數(shù)
for (i in type)#進行type次多重比較
{
g1=compare
sub_dat <- data[data[,group]==i,]#根據(jù)指定的i去取相應的數(shù)據(jù)集出來
#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)化成因子以進行多重比較
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檢驗多重比較
Tukey.labels <- data.frame(Letters=tuk$mcletters$Letters, stringsAsFactors = FALSE)#獲取多重比較字母標注
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ù)標準差
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ù)這個函數(shù)核心
是cld(glht(fit, alternative = 'two.sided’, linfct = mcp(g1 = 'Tukey’)), decreasing = TRUE),
不會出現(xiàn)c>b>a情況(因為decreasing = TRUE,當然有的人喜歡這樣標,)和亂標字母(比如對于mean最低的點
并不一定標記成c(a>b>c時)或并不一定標記成a(c>b>a時),其只能保證有差異的數(shù)據(jù)一定是不同字母),
但是多重比較出現(xiàn)“ac”標注，沒法解決。 ONE_Tukey_HSD2函數(shù)而ONE_Tukey_HSD2核心是這個 multcompLetters2(value ~ g1, Tukey_HSD$g1[,”p adj”], sub_dat),
multcompLetters2這個函數(shù)隸屬于multcompView包，與 multcompLetters不同的是
multcompLetters2可以接受formula，而multcompLetters只接受一個兩兩比較的p值的數(shù)據(jù)框,
且可能多重比較時出現(xiàn)“ac”標注,以及出現(xiàn)c>b>a情況和亂標字母(比如對于mean最低的點
并不一定標記成c(a>b>c時)或并不一定標記成a(c>b>a時),其只能保證有差異的數(shù)據(jù)一定是不同字母)。 當然多重比較好多方法，不要局限于一種方法,
例如下面的第三種可以用library(agricolae)包中的LSD檢驗(用的“BH”校正), 當然也可以用library(agricolae)包中的
【Duncan法】(新復極差法)(SSR);
【SNK法】(Student-Newman-Keuls);
【Scheffe檢驗】; 這三種多重比較方法同LSD檢驗的用法一樣都可以避免出現(xiàn)上面提到的三種情況即： 1、 a、b、c的順序不會出現(xiàn)c>b>a; 2、不會出現(xiàn)亂標字母(比如對于mean最低的點并不一定標記成c(a>b>c時)或
并不一定標記成a(c>b>a時)，其只能保證有差異的數(shù)據(jù)一定是不同字母); 3、多重比較時出現(xiàn)“ac”標注。 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')#進行了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多樣性在不同處理下的差別運行，這里拿alpha多樣性測試，看不同alpha多樣性在不同處理下的差別。 #1
#df1 <- ONE_Tukey_HSD1(data=dat,group='alpha',compare='Group',value='v')
df1 <- ONE_Tukey_HSD1(dat,'alpha','Group','v')在此可以查看各個α多樣性下不同處理間的多重比較字母標注結(jié)果，這也是本腳本的亮點之一數(shù)據(jù)量很大的情況下，可以直接查看差異情況，不用一個個的做出圖再點開看，很是方便。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，頁面最寬為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 差異分析、顯著性標記及統(tǒng)計作圖的自動實現(xiàn)R代碼示例 multcompView: Visualizations of Paired Comparisons
|
