 在進(jìn)行兩組數(shù)據(jù)間的差異分析時(shí)，我們通常會(huì)想到使用t檢驗(yàn)。但若數(shù)據(jù)不滿足執(zhí)行t檢驗(yàn)的參數(shù)假設(shè)（例如數(shù)據(jù)分布不符合正態(tài)性，變量在本質(zhì)上就嚴(yán)重偏倚或呈現(xiàn)有序關(guān)系），無法使用t檢驗(yàn)分析時(shí)，可以考慮使用非參數(shù)的方法來完成。 數(shù)據(jù)預(yù)處理及正態(tài)性假設(shè)檢驗(yàn) library(reshape2) #讀入文件，合并分組信息，數(shù)據(jù)重排 alpha <- read.table('alpha.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE) group <- read.table('group.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE) alpha <- melt(merge(alpha, group, by = 'sample'), id = c('sample', 'group')) #選擇要比較的分組（此處查看 group1 與 group2 在 shannon 指數(shù)上是否存在顯著差異） shannon_23 <- subset(alpha, variable == 'shannon' & group %in% c('2', '3')) shannon_23$group <- factor(shannon_23$group) head(shannon_23, 10) #Shapiro-Wilk 檢驗(yàn)數(shù)據(jù)是否符合正態(tài)分布（發(fā)現(xiàn)不符合正態(tài)分布） tapply(shannon_23$value, shannon_23$group, shapiro.test) 選取的數(shù)據(jù)框“shannon_23”內(nèi)容如下所示。第一列（sample），兩組數(shù)據(jù)中所含樣本名稱；第二列（group），兩組分組名稱，且分組列已轉(zhuǎn)化為因子類型；第三列（variable），alpha多樣性指數(shù)shannon指數(shù)；第四列（value），shannon指數(shù)的數(shù)值。  通過Shapiro-Wilk檢驗(yàn)得知數(shù)據(jù)分布不滿足正態(tài)性。這里p值小于0.05表明數(shù)據(jù)違背了正態(tài)性分布的零假設(shè)。  Wilcoxon檢驗(yàn) 不符合正態(tài)性前提的數(shù)據(jù)，無法應(yīng)用t檢驗(yàn)去比較差異。我們考慮使用非參數(shù)的方法作為替代，對(duì)于兩組數(shù)據(jù)的比較，可使用wilcoxon檢驗(yàn)。類似于t檢驗(yàn)，根據(jù)樣本間是否獨(dú)立，wilcoxon檢驗(yàn)分為wilcox秩和檢驗(yàn)以及wilcox符號(hào)秩和檢驗(yàn)。 wilcox秩和檢驗(yàn) 假設(shè)樣本間是相互獨(dú)立的，直接使用wilcox秩和檢驗(yàn)去處理。它是獨(dú)立樣本t檢驗(yàn)的一種非參數(shù)替代方法。 此處使用的wilcox.test()與t檢驗(yàn)t.test()的參數(shù)很相似。wilcox_test()中默認(rèn)兩組間相互獨(dú)立（默認(rèn)參數(shù)paired = FALSE），執(zhí)行獨(dú)立樣本的wilcox秩和檢驗(yàn)；同時(shí)默認(rèn)的備擇假設(shè)是雙側(cè)的（默認(rèn)參數(shù)alternative = 'two.sided'），即執(zhí)行雙側(cè)檢驗(yàn)，可分別使用“alternative = 'less'”或“alternative = 'greater'”執(zhí)行單側(cè)wilcox檢驗(yàn)。 ##wilcox 秩和檢驗(yàn)，我們執(zhí)行了一個(gè)雙側(cè)檢驗(yàn) wilcox_test <- wilcox.test(value~group, shannon_23, paired = FALSE, alternative = 'two.sided') wilcox_test wilcox_test$p.value 由于p值（約為0.003）小于0.05，即拒絕了原假設(shè)（原假設(shè)兩組間沒有差異），group2和group3的shannon指數(shù)間存在顯著不同。  wilcox符號(hào)秩和檢驗(yàn) 假設(shè)樣本間并非相互獨(dú)立的，可考慮wilcox符號(hào)秩和檢驗(yàn)，它是非獨(dú)立樣本t檢驗(yàn)的一種非參數(shù)替代方法。例如，非獨(dú)立組設(shè)計(jì)（dependent groups design）、重復(fù)測(cè)量設(shè)計(jì)（repeated measures design）等。盡管此時(shí)你選用獨(dú)立樣本的wilcox秩和檢驗(yàn)方法也是可行的，這種分析方法本身并沒問題（僅僅是在統(tǒng)計(jì)算法上存在一些不同，相較之下可能wilcox符號(hào)秩和檢驗(yàn)會(huì)更合適一些）。 此時(shí)在wilcox.test()中設(shè)定參數(shù)“paired = TRUE”即可執(zhí)行wilcox符號(hào)秩和檢驗(yàn)。 ##wilcox 符號(hào)秩和檢驗(yàn)，我們執(zhí)行了一個(gè)雙側(cè)檢驗(yàn) wilcox_test <- wilcox.test(value~group, shannon_23, paired = TRUE, alternative = 'two.sided') wilcox_test wilcox_test$p.value 根據(jù)p值（0.039，低于0.05）可知group2和group3的shannon指數(shù)間存在顯著不同。  可視化展示 考慮作圖將兩組差異進(jìn)行可視化展示。例如，一個(gè)簡(jiǎn)單的箱線圖示例。 #boxplot() 箱線圖
boxplot(value~group, data = shannon_23, col = c('blue', 'orange'), ylab = 'Shannon', xlab = 'Group', main = 'wilcox test: p-value = 0.00295')Wilcox秩和檢驗(yàn)的一個(gè)批處理示例 相較于參數(shù)分析的t檢驗(yàn)，wilcox秩和檢驗(yàn)不必事先驗(yàn)證數(shù)據(jù)分布的正態(tài)性，因此理論上來講，只要是兩組數(shù)據(jù)間的差異分析均可使用wilcox秩和檢驗(yàn)去完成，因此其適用范圍相較于t檢驗(yàn)也更廣泛。在數(shù)據(jù)量較大的情況下（可能會(huì)存在部分?jǐn)?shù)據(jù)滿足t檢驗(yàn)分析的條件，而另一部分?jǐn)?shù)據(jù)則不滿足，無法做到全部使用t檢驗(yàn)去實(shí)現(xiàn)），可以考慮使用循環(huán)逐一挑選分組后，直接執(zhí)行wilcox秩和檢驗(yàn)進(jìn)行各兩兩分組間的差異分析。盡管這種方法比較“粗暴”，但也不失為一種備選方案。 如下將展示一個(gè)批處理示例。 網(wǎng)盤附件中提供了另一示例數(shù)據(jù)集“gene.txt”。表格中每一行為一種基因，每一列為一個(gè)樣本，交叉區(qū)域?yàn)楦骰蛟诟鳂颖局械南鄬?duì)豐度。接下來，我們期望通過wilcox秩和檢驗(yàn)，找到在group1和group2組中，具有豐度差異的基因。 ##wilcox 檢驗(yàn)批處理示例 library(doBy) #使用其中的 summaryBy() 以方便按分組計(jì)算均值、中位數(shù) #讀取數(shù)據(jù) gene <- read.table('gene.txt', sep = '\t', row.names = 1, header = TRUE, stringsAsFactors = FALSE, check.names = FALSE) group <- read.table('group.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE) result <- NULL #使用循環(huán)，逐一對(duì)各基因進(jìn)行兩組間 wilcox 秩和檢驗(yàn) for (n in 1:nrow(gene)) { gene_n <- data.frame(t(gene[n,subset(group, group %in% c(1, 2))$sample])) gene_id <- names(gene_n)[1] names(gene_n)[1] <- 'gene' gene_n$sample <- rownames(gene_n) gene_n <- merge(gene_n, group, by = 'sample', all.x = TRUE) gene_n$group <- factor(gene_n$group) p_value <- wilcox.test(gene~group, gene_n)$p.value if (!is.na(p_value) & p_value < 0.05) { stat <- summaryBy(gene~group, gene_n, FUN = c(mean, median)) result <- rbind(result, c(gene_id, as.character(stat[1,1]), stat[1,2], stat[1,3], as.character(stat[2,1]), stat[2,2], stat[2,3], p_value)) } } #輸出統(tǒng)計(jì)結(jié)果 result <- data.frame(result) names(result) <- c('gene_id', 'group1', 'mean1', 'median1', 'group2', 'mean2', 'median2', 'p_value') result$p_adjust <- p.adjust(result$p_value, method = 'BH') #推薦加個(gè) p 值校正的過程 write.table(result, 'gene.wilcox.txt', sep = '\t', row.names = FALSE, quote = FALSE)  特別說明 既然參數(shù)檢驗(yàn)的前提條件有些苛刻，自己的數(shù)據(jù)不一定都滿足參數(shù)分析的條件，那么以后需要用到組間差異比較時(shí)，直接全部使用非參數(shù)的檢驗(yàn)就不可以了？ 雖然對(duì)全部數(shù)據(jù)直接使用非參數(shù)的檢驗(yàn)方式理論上沒啥問題，但還是有點(diǎn)粗暴了一些。兩種方法（此處比較了t檢驗(yàn)和wilcox秩和檢驗(yàn)）畢竟還是有差別的，非參數(shù)方法普遍沒有參數(shù)方法嚴(yán)格。對(duì)于符合參數(shù)檢驗(yàn)條件的數(shù)據(jù)來講，使用參數(shù)檢驗(yàn)還有可能會(huì)鑒別出非參數(shù)檢驗(yàn)鑒別不到的差異，此時(shí)需要特別關(guān)注。例如，某數(shù)據(jù)符合t檢驗(yàn)的條件，t檢驗(yàn)的p值顯著，但wilcox檢驗(yàn)p值不顯著，那么這時(shí)t檢驗(yàn)的結(jié)果會(huì)相對(duì)可靠一些。 |
|
|
