我們將微生物群落組成研究分為兩個(gè)主要部分：(1)分類多樣性表伦、OTUS和類群的假設(shè)檢驗(yàn)；(2)類群間差異分析。第一個(gè)成分主要屬于單變量群落分析。第二個(gè)成分可以進(jìn)一步分為各種多變量技術(shù)用踩，如聚類和排序，以及多變量差異分析的假設(shè)檢驗(yàn)忙迁。

8.1 兩組間的差異性比較

在我們的vdr?/?小鼠研究中脐彩，目的之一是測(cè)試兩組(vdr?/?小鼠和野生型小鼠)在糞便和盲腸部位的差異。第6章利用糞便樣品計(jì)算了香農(nóng)多樣性姊扔。這里惠奸，為了說明單變量群落分析，將使用各種測(cè)試統(tǒng)計(jì)數(shù)據(jù)比較計(jì)算出的香農(nóng)多樣性恰梢。

①Two-Sample Welch’s t-Test

T統(tǒng)計(jì)量是由William Sealy Gosset在1908年提出的佛南。兩樣本t檢驗(yàn)用于檢驗(yàn)兩總體均值相等。當(dāng)測(cè)試統(tǒng)計(jì)數(shù)據(jù)服從正態(tài)分布時(shí)嵌言，它最常用嗅回。如果兩組具有相同的方差，則t統(tǒng)計(jì)量可按如下方式計(jì)算

首先呀页，加載和轉(zhuǎn)置數(shù)據(jù)集妈拌。

> abund_table=read.csv("VdrGenusCounts.csv",row.names=1, check.names=FALSE)

> abund_table<-t(abund_table)

為了將數(shù)據(jù)集中的分組信息直接合并到比較中拥坛，我們需要進(jìn)行數(shù)據(jù)管理蓬蝶。在數(shù)據(jù)集中尘分，樣本ID和組信息是字符條紋。我們首先從那里提取它們丸氛。

> grouping<-data.frame(row.names=rownames(abund_table),t(as.data.frame (strsplit(rownames(abund_table),"_"))))

> grouping$Location <- with(grouping, ifelse(X3%in%"drySt-28F", "Fecal","Cecal"))

> grouping$Group <- with(grouping,ifelse(as.factor(X2)%in% c(11,12,13,14,15),c("Vdr-/-"), c("WT")))

> grouping <- grouping[,c(4,5)]

> grouping

Location? Group

5_15_drySt-28F? ? Fecal Vdr-/-

20_12_CeSt-28F? ? Cecal Vdr-/-

1_11_drySt-28F? ? Fecal Vdr-/-

2_12_drySt-28F? ? Fecal Vdr-/-

3_13_drySt-28F? ? Fecal Vdr-/-

4_14_drySt-28F? Fecal Vdr-/-

7_22_drySt-28F? ? Fecal? ? WT

8_23_drySt-28F? ? Fecal? ? WT

9_24_drySt-28F? ? Fecal? ? WT

19_11_CeSt-28F? ? Cecal Vdr-/-

21_13_CeSt-28F? ? Cecal Vdr-/-

22_14_CeSt-28F? ? Cecal Vdr-/-

23_15_CeSt-28F? ? Cecal Vdr-/-

25_22_CeSt-28F? ? Cecal? ? WT

26_23_CeSt-28F? ? Cecal? ? WT

27_24_CeSt-28F? ? Cecal? ? WT

第6章計(jì)算了香農(nóng)多樣性培愁。為了方便起見，這里重復(fù)一遍缓窜。

> library(vegan)

> H<-diversity(abund_table, "shannon")

制作了香農(nóng)多樣性的數(shù)據(jù)框定续。

> df_H<-data.frame(sample=names(H),value=H,measure=rep("Shannon",length(H)))

然后，我們將數(shù)據(jù)框多樣性和分組相結(jié)合禾锤，形成一個(gè)新的數(shù)據(jù)框私股。

> df_G <-cbind(df_H, grouping)

> rownames(df_G)<-NULL

> df_G

sample value measure Location? Group

1? 5_15_drySt-28F 2.461 Shannon? ? Fecal Vdr-/-

2? 20_12_CeSt-28F 2.340 Shannon? ? Cecal Vdr-/-

3? 1_11_drySt-28F 2.228 Shannon? ? Fecal Vdr-/-

4? 2_12_drySt-28F 2.734 Shannon? ? Fecal Vdr-/-

5? 3_13_drySt-28F 2.077 Shannon? ? Fecal Vdr-/-

6? 4_14_drySt-28F 2.467 Shannon? ? Fecal Vdr-/-

7? 7_22_drySt-28F 1.777 Shannon? ? Fecal? ? WT

8? 8_23_drySt-28F 2.000 Shannon? ? Fecal? ? WT

9? 9_24_drySt-28F 1.972 Shannon? ? Fecal? ? WT

10 19_11_CeSt-28F 1.345 Shannon? Cecal Vdr-/-

11 21_13_CeSt-28F 2.016 Shannon? ? Cecal Vdr-/-

12 22_14_CeSt-28F 1.955 Shannon? ? Cecal Vdr-/-

13 23_15_CeSt-28F 1.614 Shannon? ? Cecal Vdr-/-

14 25_22_CeSt-28F 1.959 Shannon? ? Cecal? ? WT

15 26_23_CeSt-28F 2.271 Shannon? ? Cecal? ? WT

16 27_24_CeSt-28F 2.002 Shannon? ? Cecal? ? WT

接下來，我們將新數(shù)據(jù)框中的糞便數(shù)據(jù)子集提出

> Fecal_G<- subset(df_G, Location=="Fecal")

> Fecal_G

sample value measure Location? Group

1 5_15_drySt-28F 2.461 Shannon? ? Fecal Vdr-/-

3 1_11_drySt-28F 2.228 Shannon? ? Fecal Vdr-/-

4 2_12_drySt-28F 2.734 Shannon? ? Fecal Vdr-/-

5 3_13_drySt-28F 2.077 Shannon? ? Fecal Vdr-/-

6 4_14_drySt-28F 2.467 Shannon? ? Fecal Vdr-/-

7 7_22_drySt-28F 1.777 Shannon? ? Fecal? ? WT

8 8_23_drySt-28F 2.000 Shannon? ? Fecal? ? WT

9 9_24_drySt-28F 1.972 Shannon? ? Fecal? ? WT

現(xiàn)在數(shù)據(jù)已經(jīng)準(zhǔn)備好進(jìn)行統(tǒng)計(jì)分析了恩掷。在進(jìn)行假設(shè)檢驗(yàn)之前倡鲸，讓我們使用ggplot2包中的函數(shù)gglot()和plyr包中的ddply()函數(shù)來研究Shannon多樣性值的分布。

> library(ggplot2)

我們使用facet_grid將繪圖拆分成兩個(gè)面板黄娘。

> p<-ggplot(Fecal_G, aes(x=value))+

+? geom_histogram(color="black", fill="black")+

+? facet_grid(Group * .)

plyr包用于計(jì)算每個(gè)組的平均香農(nóng)多樣性值峭状。

> library(plyr)

> mu <- ddply(Fecal_G, "Group", summarise, grp.mean=mean(value))

> head(mu)

Group grp.mean

1 Vdr-/-? 2.393

2? ? WT? ? 1.916

將平均線添加到繪圖中。

> p+geom_vline(data=mu, aes(xintercept=grp.mean, color="red"),

+? linetype="dashed")

該分布(圖8.1)顯示VDR基因敲除導(dǎo)致更高的分集逼争，因?yàn)樵摻M的直方圖相對(duì)于WT組向右移動(dòng)(更高的分集值)优床。為了檢驗(yàn)香農(nóng)多樣性無差異的零假設(shè)，使用Welch的t檢驗(yàn)誓焦，結(jié)果是p值=0.01(t=3.6胆敞，df=5.9)。因此杂伟，我們拒絕無差異的零假設(shè)竿秆，而支持香農(nóng)多樣性在兩組中不同的替代方案。

> fit_t <- t.test(value ~ Group, data=Fecal_G)

> fit_t

Welch Two Sample t-test

data:? value by Group

t = 3.6, df = 5.9, p-value = 0.01

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

0.1518 0.8026

sample estimates:

mean in group Vdr-/- mean in group WT

2.393? ? ? ? ? ? ? ? 1.916

② Wilcoxon Rank Sum Test（Wilcoxon秩和檢驗(yàn)）

Wilcoxon秩和檢驗(yàn)等價(jià)于Mann和Whitney(1947)提出的Mann-Whitney U檢驗(yàn)稿壁。它是兩樣本t檢驗(yàn)的一種非參數(shù)替代幽钢，它使用兩個(gè)獨(dú)立樣本數(shù)據(jù)的秩來檢驗(yàn)零假設(shè)：兩個(gè)獨(dú)立樣本來自具有相同分布的總體(即，兩個(gè)總體相同)傅是。與t-檢驗(yàn)不同的是匪燕，Wilcoxon秩和檢驗(yàn)不需要正態(tài)分布的假設(shè)，并且?guī)缀跖ct-檢驗(yàn)一樣有效喧笔。因此帽驯，它在微生物組研究中得到了廣泛的應(yīng)用。執(zhí)行Wilcoxon秩和檢驗(yàn)以求出檢驗(yàn)統(tǒng)計(jì)量的值需要三個(gè)主要步驟：

步驟1. 給所有的觀察值分配排名书闸，最小值的排名為1尼变。如果值是平分的，則分配平局中涉及的排名的平均值。

步驟2. 對(duì)兩個(gè)樣本中任何一個(gè)樣本的秩進(jìn)行求和嫌术“С海可以確定另一樣本中的秩和，因?yàn)樗兄鹊暮偷扔贜(N+1)/2度气，其中N是觀察的總數(shù)割按。如果兩個(gè)測(cè)試總體具有相同的分布，則秩R的平均值為? ? ? ? ? ? ? ? ? ? ?磷籍，標(biāo)準(zhǔn)偏差為? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?适荣。當(dāng)秩和R遠(yuǎn)離均值時(shí)，Wilcoxon秩和檢驗(yàn)拒絕兩個(gè)總體具有相同分布的假設(shè)院领。隨著兩個(gè)樣本量的增加弛矛，秩和統(tǒng)計(jì)量變得近似正常。我們可以通過標(biāo)準(zhǔn)化秩和來形成統(tǒng)計(jì)量比然。

步驟3.使用以下公式計(jì)算z檢驗(yàn)統(tǒng)計(jì)量的值：

R?編號(hào)為n1的樣本的秩和汪诉；n1? ?找到秩和R的樣本大小(如樣本1)；n2 另一個(gè)樣本大小(如樣本2)谈秫。

以下代碼用于進(jìn)行Wilcoxon秩和檢驗(yàn)扒寄。

> fit_w <- wilcox.test(value * Group, data=Fecal_G)

> fit_w

Wilcoxon rank sum test

data: value by Group

W = 15, p-value = 0.04

alternative hypothesis: true location shift is not equal to 0

圖8.1顯示樣本大小較小時(shí)分布是不對(duì)稱的。Wilcoxon秩和檢驗(yàn)可能更合適拟烫，但由Wilcoxon秩和檢驗(yàn)給出的p=0.04與Welcht檢驗(yàn)在0.05顯著性水平上p=0.01得出的結(jié)論相同该编。

8.2 Comparisons of a Taxon of Interest Between Two Groups

①用Wilcoxon秩和檢驗(yàn)比較相對(duì)豐度

當(dāng)我們分析微生物組豐度數(shù)據(jù)時(shí)，從樣品中OTUs的豐度或類群的豐度來推斷生態(tài)系統(tǒng)的總豐度是不合適的硕淑。取而代之的是课竣，我們可以利用樣本中的相對(duì)豐度來推斷其在生態(tài)系統(tǒng)中的分類單元的相對(duì)豐度。其背后的原因是它存在一個(gè)組成約束：樣本和中的所有微生物相對(duì)豐度為1置媳，這導(dǎo)致組成數(shù)據(jù)駐留在單純形(Aitchison 1982,1986)而不是歐幾里德空間于樟。因此，通常將數(shù)據(jù)標(biāo)準(zhǔn)化到一個(gè)共同的尺度拇囊，以便于跨組的分類單元豐度的比較迂曲。方法是將分類單元計(jì)數(shù)除以讀取總次數(shù)100，以將豐度轉(zhuǎn)換為樣本中讀取的百分比寥袭，將數(shù)據(jù)縮放為“每100次讀取的分類單元數(shù)量”路捧。當(dāng)我們選擇一個(gè)特定的單個(gè)分類單元進(jìn)行跨組測(cè)試時(shí)，重要的是要確保指定的分類單元是基于假設(shè)或理論的传黄，以減少夸大假陽性率的機(jī)會(huì)(即杰扫，在不應(yīng)該拒絕零假設(shè)的時(shí)候拒絕它)。小鼠體內(nèi)的VDR顯著影響β多樣性膘掰，并持續(xù)影響單個(gè)細(xì)菌分類群章姓，如副細(xì)菌分類群(Wang等人。2016)。在本節(jié)中凡伊，我們使用糞便樣本說明Wilcoxon秩和檢驗(yàn)來比較VDR小鼠數(shù)據(jù)集中的細(xì)菌類桿菌零渐。

首先，檢查每個(gè)樣本的總豐度窗声。

> apply(abund_table,1, sum)

5_15_drySt-28F 20_12_CeSt-28F 1_11_drySt-28F 2_12_drySt-28F 3_13_drySt-28F 4_14_drySt-28F 7_22_drySt-28F 8_23_drySt-28F 9_24_drySt-28F 19_11_CeSt-28F

1853? ? ? ? ? ? ? ? ? ? ? ? 3239? ? ? ? ? ? ? ? ? ? ? ? ? ?6211? ? ? ? ? ? ? ? ?5115? ? ? ? ? ? ? ?6016? ? ? ? ? ? ? ? ? ? ? ? ?2343? ? ? ? ? ? ? ? 2262? ? ? ? ? ? ? ? 7255? ? ? ? ? ? ? ? ? ? ?5502? ? ? ? ? ? ? ? ?5067

21_13_CeSt-28F 22_14_CeSt-28F 23_15_CeSt-28F 25_22_CeSt-28F 26_23_CeSt-28F? 27_24_CeSt-28F

2397? ? ? ? ? ? ? ? ? ? ? 3788? ? ? ? ? ? ? ? ? ? ? ? ? 9264? ? ? ? ? ? ? ? ? ?2072? ? ? ? ? ? ? ? ? ? 6903? ? ? ? ? ? ? ? ? ? ? ? ? ?6327

然后，通過將每個(gè)值除以樣本總豐度來計(jì)算相對(duì)豐度辜纲。

> relative_abund_table<-decostand(abund_table, method = "total")

檢查每個(gè)樣品的總豐度笨觅，使上述計(jì)算正確。

> apply(relative_abund_table, 1, sum)

5_15_drySt-28F 20_12_CeSt-28F 1_11_drySt-28F 2_12_drySt-28F 3_13_drySt-28F? 4_14_drySt-28F 7_22_drySt-28F 8_23_drySt-28F 9_24_drySt-28F 19_11_CeSt-28F

1? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1? ? ? ? ? ? ? ? ? ? ? ? ? ? 1? ? ? ? ? ? ? ? ? ? ? ? ? ?1? ? ? ? ? ? ? ?1? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1? ? ? ? ? ? ? ? ? ? ? ? ? 1? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1? ? ? ? ? ? ? ? ? ? 1? ? ? ? ? ? ? ? ? ? ? ? 1

21_13_CeSt-28F 22_14_CeSt-28F 23_15_CeSt-28F 25_22_CeSt-28F 26_23_CeSt-28F? ? 27_24_CeSt-28F

1? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1? ? ? ? ? ? ? ? ? 1? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1

請(qǐng)看一下轉(zhuǎn)換后的數(shù)據(jù)耕腾。

> relative_abund_table[1:16,1:8]

? ? ? ? ? ? ? ? ? ? ? ? Tannerella Lactococcus Lactobacillus Lactobacillus::Lactococcus

5_15_drySt-28F? 0.256881? ? 0.17593? ? ? 0.05073? ? ? ? ? ? ? ? ? 0.0005397

20_12_CeSt-28F? 0.020685? ? 0.22754? ? ? 0.18432? ? ? ? ? ? ? ? ? 0.0037048

1_11_drySt-28F? 0.088392? ? 0.36983? ? ? 0.06988? ? ? ? ? ? ? ? ? 0.0040251

2_12_drySt-28F? 0.113001? ? 0.10714? ? ? 0.14057? ? ? ? ? ? ? ? ? 0.0009775

3_13_drySt-28F? 0.165559? ? 0.39528? ? ? 0.05352? ? ? ? ? ? ? ? ? 0.0028258

4_14_drySt-28F? 0.172429? ? 0.20102? 0.08749? ? ? ? ? ? ? ? ? 0.0004268

7_22_drySt-28F? 0.141026? ? 0.38992? ? ? 0.28470? ? ? ? ? ? ? ? ? 0.0057471

8_23_drySt-28F? 0.072502? ? 0.27195? ? ? 0.32254? ? ? ? ? ? ? ? ? 0.0020675

9_24_drySt-28F? 0.077063? ? 0.41948? ? ? 0.18175? ? ? 0.0025445

19_11_CeSt-28F? 0.000000? ? 0.08328? ? ? 0.06513? ? ? ? ? ? ? ? ? 0.0013815

21_13_CeSt-28F? 0.002503? ? 0.07217? ? ? 0.26658? ? ? ? ? ? ? ? ? 0.0000000

22_14_CeSt-28F? 0.005280? ? 0.15312? ? ? 0.16711? ? ? ? ? ? ? ? ? 0.0007920

23_15_CeSt-28F? 0.003994? ? 0.52537? ? ? 0.19635? ? ? ? ? ? ? ? ? 0.0026986

25_22_CeSt-28F? 0.018340? ? 0.34122? ? ? 0.30164? ? ? ? ? ? ? ? ? 0.0043436

26_23_CeSt-28F? 0.011734? ? 0.20339? ? ? 0.19716? ? ? ? ? ? ? ? ? 0.0014486

27_24_CeSt-28F? 0.037142? ? 0.30235? ? ? 0.05769? ? ? ? ? ? ? ? ? 0.0020547

? ? ? ? ? ? ? ? ? ? ? ? ? ? ?Parasutterella Helicobacter Prevotella Bacteroides

5_15_drySt-28F? ? ? 0.0005397? ? 0.048030? 0.0652995? ? 0.147329

20_12_CeSt-28F? ? ? 0.0000000? ? 0.000000? 0.0021612? 0.010497

1_11_drySt-28F? ? ? 0.0001610? ? 0.000000? 0.0465303? ? 0.154242

2_12_drySt-28F? ? ? 0.0007820? ? 0.002542? 0.0193548? ? 0.073705

3_13_drySt-28F? ? ? 0.0003324? ? 0.003989? 0.0556848? ? 0.087434

4_14_drySt-28F? ? ? 0.0000000? ? 0.013658? 0.0610329? ? 0.085361

7_22_drySt-28F? ? ? 0.0000000? ? 0.001326? 0.0490716? ? 0.038019

8_23_drySt-28F? ? ? 0.0016540? ? 0.000000? 0.0122674? ? 0.058442

9_24_drySt-28F? ? ? 0.0001818? ? 0.000000? 0.0152672? ? 0.036714

19_11_CeSt-28F 0.0000000 0.000000 0.0000000 0.000000

21_13_CeSt-28F? ? ? 0.0000000? ? 0.000000? 0.0004172? ? 0.002086

22_14_CeSt-28F 0.0000000 0.000000 0.0007920 0.005280

23_15_CeSt-28F? ? ? 0.0002159? ? 0.000000? 0.0010794? ? 0.003346

25_22_CeSt-28F? ? ? 0.0000000? 0.000000? 0.0033784? ? 0.009170

26_23_CeSt-28F? ? ? 0.0002897? ? 0.000000? 0.0007243? ? 0.006664

27_24_CeSt-28F? ? ? 0.0000000? ? 0.000000? 0.0039513? ? 0.019282

我們感興趣的細(xì)菌類桿菌在第8欄见剩。讓我們對(duì)這個(gè)細(xì)菌進(jìn)行分類。

> (Bacteroides <-relative_abund_table[,8])

5_15_drySt-28F 20_12_CeSt-28F 1_11_drySt-28F 2_12_drySt-28F 3_13_drySt-28F??4_14_drySt-28F 7_22_drySt-28F 8_23_drySt-28F 9_24_drySt-28F 19_11_CeSt-28F

0.147329? ? ? ? ? ? ? ?0.010497? ? ? ? ? ? ? ?0.154242? ? ? ? ? ? ? ? ? 0.073705? ? ? 0.087434? ? ? ?0.085361? ? ? ? ? ? ? 0.038019? ? ? ? ? ? ? ? 0.058442? ? ? ? ? ? ? ? 0.036714? ? ? 0.000000

21_13_CeSt-28F 22_14_CeSt-28F 23_15_CeSt-28F 25_22_CeSt-28F 26_23_CeSt-28F?27_24_CeSt-28F

0.002086? ? ? ? ? ? ? ? ? 0.005280? ? ? ? ? ? ?0.003346? ? ? ? ? ? ? ? ?0.009170? ? ? ? ? ? ?0.006664? ? ? ? ? ? ? ? ?0.019282

現(xiàn)在扫俺，將類桿菌和分組數(shù)據(jù)框組合在一起苍苞，并對(duì)糞便樣本進(jìn)行子集以供以后使用。

Bacteroides Location Group

1? 0.14732866? ? Fecal Vdr-/-

3? 0.15424247? ? Fecal Vdr-/-

4? 0.07370479? ? Fecal Vdr-/-

5? 0.08743351? ? Fecal Vdr-/-

6? 0.08536065? ? Fecal Vdr-/-

7? 0.03801945? ? Fecal? ? WT

8? 0.05844245? ? Fecal? ? WT

9? 0.03671392? ? Fecal? ? WT

boxlot()函數(shù)用于生成帶有組的簡(jiǎn)單類桿菌箱線圖(圖8.2)狼纬。

Bacteroides Location Group

1? 0.14733? ? Fecal Vdr-/-

3? ? 0.15424? ? Fecal Vdr-/-

4? ? 0.07370? ? Fecal Vdr-/-

5? ? 0.08743? ? Fecal Vdr-/-

6? ? 0.08536? ? Fecal Vdr-/-

7? ? 0.03802? ? Fecal? ? WT

8? ? 0.05844? ? Fecal? ? WT

9? ? 0.03671? ? Fecal? ? WT

> boxplot(Bacteroides * Group,data=Fecal_Bacteroides_G, col=rainbow(2),main="Bacteroides in Vdr WT/KO mice")

以下代碼用于使用函數(shù)ggplot()生成箱線圖(圖8.3)羹呵。

> ggplot(Fecal_Bacteroides_G, aes(x=Group, y=Bacteroides,col=factor(Group)))+ geom_boxplot(notch=FALSE)

> ggplot(Fecal_Bacteroides_G, aes(x=Group, y=Bacteroides)) +

geom_boxplot(outlier.colour="red", outlier.shape=8, outlier.size=4) +

layer(stat_params = list(binwidth = 2))

箱線圖顯示野生型(WT，n=3)和VDR基因敲除小鼠(KO疗琉，n=5)的類群(Bacteriodes)冈欢。

> fit_w_b=wilcox.test(Bacteroides*Group,data=Fecal_Bacteroides_G)

> fit_w_b

Wilcoxon rank sum test

data: Bacteroides by Group

W = 15, p-value = 0.04

alternative hypothesis: true location shift is not equal to 0

上述Wolcoxon檢驗(yàn)表明，Vdr?/?小鼠和WT小鼠的Bacteriodes相對(duì)豐度有統(tǒng)計(jì)學(xué)意義盈简。我們可以得出結(jié)論凑耻，VDR基因敲除富集了Bacteriodes。

②有無紫杉醇的卡方檢驗(yàn)比較（Chi-Square Test）

卡方檢驗(yàn)柠贤，也被寫成X2 test香浩，通常被用作皮爾遜卡方檢驗(yàn)的縮寫，被提出臼勉，并由卡爾皮爾遜在1900年(皮爾遜1900)首次研究其性質(zhì)邻吭。X2test應(yīng)用于分類數(shù)據(jù)集，以檢驗(yàn)觀察到的頻率分布是否與聲稱或理論分布不同(t擬合優(yōu)度)宴霸，并調(diào)查列聯(lián)表中的行變量和列變量是否彼此獨(dú)立(獨(dú)立性檢驗(yàn))镜盯。擬合優(yōu)度的測(cè)試統(tǒng)計(jì)量如下所示：

作為經(jīng)驗(yàn)法則，它要求所有預(yù)期的細(xì)胞計(jì)數(shù)等于或超過5猖败，以提供對(duì)卡方分布的充分近似(Wackerly等人速缆。2002)，盡管Cochran(1952)指出恩闻，在某些情況下艺糜，該值可能低至1。在本節(jié)中，我們使用盲腸樣本說明X2test來比較VDR小鼠數(shù)據(jù)集中的細(xì)菌Parabacteroides破停。為了說明X2檢驗(yàn)翅楼，我們將Parabacteroides的豐度計(jì)數(shù)數(shù)據(jù)轉(zhuǎn)換為二進(jìn)制變量。如果樣本中沒有分類單元真慢，則豐度表中的Parabacteroides分類單元的計(jì)數(shù)數(shù)據(jù)將轉(zhuǎn)換為0毅臊，或者如果樣本中存在分類單元，則轉(zhuǎn)換為1黑界。表8.1匯總了轉(zhuǎn)換后的數(shù)據(jù)管嬉。

首先，查看豐富的數(shù)據(jù)來鑒定細(xì)菌副細(xì)菌朗鸠，并對(duì)該細(xì)菌進(jìn)行亞群劃分蚯撩。

> abund_table[1:16,1:27]

> (Parabacteroides <- abund_table[,27])

然后，將子集數(shù)據(jù)與分組數(shù)據(jù)框合并烛占。

> Parabacteroides_G <-cbind(Parabacteroides, grouping)

> rownames(Parabacteroides_G)<-NULL

由于組合數(shù)據(jù)框既包括糞便樣本又包括盲腸樣本胎挎，因此讓我們從該數(shù)據(jù)框中將盲腸數(shù)據(jù)子集提出。

> Cecal_Parabacteroides_G <- subset(Parabacteroides_G, Location=="Cecal")

> Cecal_Parabacteroides_G

Parabacteroides Location? Group

2? ? ? ? ? ? ? ? 0? ? Cecal Vdr-/-

10? ? ? ? ? ? ? 0? ? Cecal Vdr-/-

11? ? ? ? ? ? ? 1? ? Cecal Vdr-/-

12? ? ? ? ? ? ? 4? ? Cecal Vdr-/-

13? ? ? ? ? ? ? 15? ? Cecal Vdr-/-

14? ? ? ? ? ? ? 5? ? Cecal? ? WT

15? ? ? ? ? ? ? 4? Cecal? ? WT

16? ? ? ? ? ? ? 6? ? Cecal? ? WT

重新編碼用于卡方檢驗(yàn)的二進(jìn)制變量“Present”忆家。

> Cecal_Parabacteroides_G$Present <- ifelse((Cecal_Parabacteroides_G$Parabacteroides > 0), "Present","Absent")

> Cecal_Parabacteroides_G

Parabacteroides Location? Group Present

2? ? ? ? ? ? ? ? 0? ? Cecal Vdr-/-? Absent

10? ? ? ? ? ? ? 0? ? Cecal Vdr-/-? Absent

11? ? ? ? ? ? ? 1? ? Cecal Vdr-/- Present

12? ? ? ? ? ? ? 4? ? Cecal Vdr-/- Present

13? ? ? ? ? ? ? 15? ? Cecal Vdr-/- Present

14? ? ? ? ? ? ? 5? ? Cecal? ? WT Present

15? ? ? ? ? ? ? 4? ? Cecal? ? WT Present

16? ? ? ? ? ? ? 6? ? Cecal? ? WT Present

以下代碼用于進(jìn)行卡方檢驗(yàn)犹菇。

> library(MASS)

> tbl = table(Cecal_Parabacteroides_G$Group, Cecal_Parabacteroides_G$Present)

> tbl? ? ? ? ? ? ? ?

Absent Present

Vdr-/-? 2? ? ? 3

WT? ? ? ? ? 0? ? ? 3

> chisq.test(tbl)

Pearson's Chi-squared test with Yates' continuity correction

data:? tbl

X-squared = 0.18, df = 1, p-value = 0.7

Warning message:

In chisq.test(tbl) : Chi-squared approximation may be incorrect

表8.1顯示了這些比率的分布-5個(gè)vdr?/?盲腸樣本中有3個(gè)(60%)有Parabacteroides，而3個(gè)野生型樣本中有3個(gè)(100%)有Parabacteroides芽卿。對(duì)兩組發(fā)病率無差異的零假設(shè)進(jìn)行卡方檢驗(yàn)项栏，P值為0.7(X-平方=0.18，df=1)蹬竖，不能排除兩組發(fā)病率無差異的零假設(shè)沼沈，得出兩組發(fā)病率無差異的結(jié)論。請(qǐng)注意币厕，由于樣本量較小列另，輸出中會(huì)出現(xiàn)一條警告消息。通常旦装，如果列聯(lián)表中的單元格值很小(例如<5)页衙，卡方檢驗(yàn)可能不正確，則應(yīng)用Fisher精確檢驗(yàn)阴绢。

> fisher.test(tbl)

Fisher’s Exact Test for Count Data

data: tbl

p-value = 0.5

alternative hypothesis: true odds ratio is not equal to 1

95 percent confidence interval:

0.109? Inf

sample estimates:

odds ratio

Inf

Fisher‘s精確檢驗(yàn)結(jié)果也不顯著店乐，p值為0.5，與卡方檢驗(yàn)結(jié)果一致呻袭。然而眨八，這項(xiàng)測(cè)試很難在無限的置信區(qū)間內(nèi)得出結(jié)論。

8.3 Comparisons Among More than Two Groups Using ANOVA（ 兩組以上組間比較的方差分析）

①One-Way ANOVA（單因素方差分析）

方差分析(ANOVA)是由羅納德·費(fèi)希爾(Ronald Fisher左电，1918)于1918年提出的廉侧，1925年費(fèi)舍爾的著作“研究人員的統(tǒng)計(jì)方法”出版后页响，方差分析變得廣為人知。方差分析將兩樣本t-test推廣到兩組以上段誊。方差分析的零假設(shè)是：比較組的所有均值相等闰蚕。使用方差分析的分析依賴于基礎(chǔ)數(shù)據(jù)的正態(tài)性假設(shè)。然而连舍，大多數(shù)微生物群落組成數(shù)據(jù)没陡，特別是多變量數(shù)據(jù)并不是正態(tài)分布的，因此索赏，本書僅使用方差分析比較α多樣性測(cè)度的單變量分析盼玄。對(duì)于多變量群落組成數(shù)據(jù)，要么采用非參數(shù)形式的方差分析参滴，要么采用其他合適的統(tǒng)計(jì)方法强岸。檢驗(yàn)統(tǒng)計(jì)量的形成是通過傳統(tǒng)的平方和分割(方差分割)來實(shí)現(xiàn)的锻弓。樣本方差的定義方程為

總偏差因素比較采用F檢驗(yàn)砾赔。F值是通過將處理之間的方差除以處理內(nèi)的方差來獲得的。下面給出了單因素方差分析中的F檢驗(yàn)統(tǒng)計(jì)量

為了說明方差分析在微生物群落組成研究中的作用青灼，我們使用了我們的VDR?/?小鼠數(shù)據(jù)暴心。本研究的一個(gè)假設(shè)是VDR狀態(tài)和腸道位置對(duì)腸道細(xì)菌群落沒有影響。我們使用方差分析來分析 Chao 1α多樣性測(cè)定來解決這一假設(shè)杂拨。復(fù)制了第6章中使用的代碼下面的Chao 1豐富度測(cè)定专普。結(jié)果在一個(gè)數(shù)據(jù)框中。合并用于方差分析檢驗(yàn)的組信息弹沽。下面的代碼構(gòu)成了一個(gè)chao1豐富性的數(shù)據(jù)框檀夹，并將組信息添加到該數(shù)據(jù)框中。

> CH=estimateR(abund_table)[2,]

> df_CH <-data.frame(sample=names(CH),value=CH,measure=rep("Chao1",length(CH)))

> df_CH_G <-cbind(df_CH, grouping)

> rownames(df_G)<-NULL

> df_CH_G

sample? value measure Location? Group

5_15_drySt-28F 5_15_drySt-28F? 94.75? Chao1? ? Fecal Vdr-/-

20_12_CeSt-28F 20_12_CeSt-28F? 59.80? Chao1? ? Cecal Vdr-/-

1_11_drySt-28F 1_11_drySt-28F? 77.00? Chao1? ? Fecal Vdr-/-

2_12_drySt-28F 2_12_drySt-28F 103.27? Chao1? ? Fecal Vdr-/-

3_13_drySt-28F 3_13_drySt-28F? 85.67? Chao1? ? Fecal Vdr-/-

4_14_drySt-28F 4_14_drySt-28F? 55.14? Chao1? ? Fecal Vdr-/-

7_22_drySt-28F 7_22_drySt-28F? 62.75? Chao1? ? Fecal? WT

8_23_drySt-28F 8_23_drySt-28F? 67.67? Chao1? ? Fecal? ? WT

9_24_drySt-28F 9_24_drySt-28F? 80.50? Chao1? ? Fecal? ? WT

19_11_CeSt-28F 19_11_CeSt-28F? 52.17? Chao1? ? Cecal Vdr-/-

21_13_CeSt-28F 21_13_CeSt-28F? 55.00? Chao1? ? Cecal Vdr-/-

22_14_CeSt-28F 22_14_CeSt-28F? 59.00? Chao1? ? Cecal Vdr-/-

23_15_CeSt-28F 23_15_CeSt-28F? 60.88? Chao1? ? Cecal Vdr-/-

25_22_CeSt-28F 25_22_CeSt-28F? 51.00? Chao1? ? Cecal? ? WT

26_23_CeSt-28F 26_23_CeSt-28F 112.86? Chao1? ? Cecal? ? WT

27_24_CeSt-28F 27_24_CeSt-28F? 78.06? Chao1? ? Cecal? ? WT

新的四個(gè)級(jí)別的組是通過位置和組的交互作用生成的策橘。

> df_CH_G$Group4<- with(df_CH_G, interaction(Location,Group))

> df_CH_G

sample? value measure Location? Group? ? ? Group4

5_15_drySt-28F 5_15_drySt-28F? 94.75? Chao1? ? Fecal Vdr-/- Fecal.Vdr-/-

20_12_CeSt-28F 20_12_CeSt-28F? 59.80? Chao1? ? Cecal Vdr-/- Cecal.Vdr-/-

1_11_drySt-28F 1_11_drySt-28F? 77.00? Chao1? ? Fecal Vdr-/- Fecal.Vdr-/-

2_12_drySt-28F 2_12_drySt-28F 103.27? Chao1? ? Fecal Vdr-/- Fecal.Vdr-/-

3_13_drySt-28F 3_13_drySt-28F? 85.67? Chao1? ? Fecal Vdr-/- Fecal.Vdr-/-

4_14_drySt-28F 4_14_drySt-28F? 55.14? Chao1? ? Fecal Vdr-/- Fecal.Vdr-/-

7_22_drySt-28F 7_22_drySt-28F? 62.75? Chao1? ? Fecal? ? WT? ? Fecal.WT

8_23_drySt-28F 8_23_drySt-28F? 67.67? Chao1? ? Fecal? ? WT? ? Fecal.WT

9_24_drySt-28F 9_24_drySt-28F? 80.50? Chao1? ? Fecal? ? WT? ? Fecal.WT

19_11_CeSt-28F 19_11_CeSt-28F? 52.17? Chao1? ? Cecal Vdr-/- Cecal.Vdr-/-

21_13_CeSt-28F 21_13_CeSt-28F? 55.00? Chao1? ? Cecal Vdr-/- Cecal.Vdr-/-

22_14_CeSt-28F 22_14_CeSt-28F? 59.00? Chao1? ? Cecal Vdr-/- Cecal.Vdr-/-

23_15_CeSt-28F 23_15_CeSt-28F? 60.88? Chao1? ? Cecal Vdr-/- Cecal.Vdr-/-

25_22_CeSt-28F 25_22_CeSt-28F? 51.00? Chao1? ? Cecal? ? WT? ? Cecal.WT

26_23_CeSt-28F 26_23_CeSt-28F 112.86? Chao1? ? Cecal? ? WT? ? Cecal.WT

27_24_CeSt-28F 27_24_CeSt-28F? 78.06? Chao1? ? Cecal? ? WT? ? Cecal.WT

使用boxplot()探索Chao1指數(shù)(圖8.4)丽已。

> boxplot(value*Group4, data=df_CH_G, col=rainbow(4), main="Chao1 index")

下面的ggplot()生成高質(zhì)量的框圖以供發(fā)表使用(圖8.5)沛婴。

> p <- ggplot(df_CH_G, aes(x=Group4, y=value),col=rainbow(4), main="Chao1 index") + geom_boxplot()

> p + coord_flip()> ggplot(df_CH_G, aes(x=Group4, y=value,col=factor(Group4))) +

?+ geom_boxplot(notch=FALSE)

除了目視檢查底層數(shù)據(jù)的正態(tài)性外嘁灯，還可以測(cè)試方差的均勻性丑婿。Sokal和Rohlf(2011)描述了三個(gè)這樣的檢驗(yàn)：Bartlett的同質(zhì)性檢驗(yàn)藕夫、Hartley的Fmaxtest和log-anova檢驗(yàn)孽糖，或SchefféBox檢驗(yàn)。為了繼續(xù)使用方差分析進(jìn)行驗(yàn)證毅贮，我們必須首先檢驗(yàn)方差的均勻性办悟。軟件R提供兩個(gè)測(cè)試：Bartlett測(cè)試和Fligner-Killeen測(cè)試。為了說明方差的均勻性檢驗(yàn)滩褥，我們使用了來自糞便和盲腸位置的Vdr?/?和WT小鼠數(shù)據(jù)的Chao 1豐富度度量病蛉。零假設(shè)(H0)是四組中的所有方差都是相同的。我們從Bartlett測(cè)驗(yàn)開始瑰煎。為了方便處理Bartlett檢驗(yàn)铺然，我們使用dplyr包中select()函數(shù)來選擇相關(guān)的group and Chao 1值列。

> library(dplyr)

> df_CH_G4 <- select(df_CH_G, Group4,value)

> df_CH_G4

Group4? value

5_15_drySt-28F Fecal.Vdr-/-? 94.75

20_12_CeSt-28F Cecal.Vdr-/-? 59.80

1_11_drySt-28F Fecal.Vdr-/-? 77.00

2_12_drySt-28F Fecal.Vdr-/- 103.27

3_13_drySt-28F Fecal.Vdr-/-? 85.67

4_14_drySt-28F Fecal.Vdr-/-? 55.14

7_22_drySt-28F? ? Fecal.WT? 62.75

8_23_drySt-28F? ? Fecal.WT? 67.67

9_24_drySt-28F? ? Fecal.WT? 80.50

19_11_CeSt-28F Cecal.Vdr-/-? 52.17

21_13_CeSt-28F Cecal.Vdr-/-? 55.00

22_14_CeSt-28F Cecal.Vdr-/-? 59.00

23_15_CeSt-28F Cecal.Vdr-/-? 60.88

25_22_CeSt-28F? ? Cecal.WT? 51.00

26_23_CeSt-28F? ? Cecal.WT 112.86

27_24_CeSt-28F? ? Cecal.WT? 78.06

下面的R代碼執(zhí)行方差的同質(zhì)性的Bartlett 檢驗(yàn)：

> bartlett.test(df_CH_G4, Group4)

Bartlett test of homogeneity of variances

data: df_CH_G4

Bartlett’s K-squared = 62, df = 1, p-value = 3e-15

> qchisq(0.95, 1)

?[1] 3.841

該函數(shù)給出了統(tǒng)計(jì)檢驗(yàn)的K平方值和p值酒甸。它表明魄健，在5%的水平上可以拒絕零假設(shè)〔迩冢或者沽瘦，我們可以將Bartlett的K平方與卡方表格的值進(jìn)行比較，在qchisq()函數(shù)中使用相同級(jí)別的α和自由度农尖。如果X平方>Bartlett‘sK平方助隧，我們接受零假設(shè)H0(方差齊性)并村，否則拒絕零假設(shè)。因?yàn)榭ǚ叫∮贐artlett的K平方姐叁，所以我們拒絕零假設(shè)H0，并得出方差不相同的結(jié)論。我們現(xiàn)在使用Fligner-Killeen檢驗(yàn)來檢驗(yàn)同方差滔驾。如下所示的語法非常相似绕德。

> fligner.test(df_CH_G4, Group4)

Fligner-Killeen test of homogeneity of variances

data: df_CH_G4

Fligner-Killeen:med chi-squared = 21, df = 1, p-value = 6e-06

結(jié)論與Bartlett檢驗(yàn)相似：方差不同。然而臣咖，為了說明的目的，我們繼續(xù)對(duì)數(shù)據(jù)進(jìn)行方差分析刁赦，而不考慮方差的齊性檢驗(yàn)結(jié)果趴荸。以下R碼符合該模型：

> fit = lm(formula = value*Group4,data=df_CH_G)

然后對(duì)方差分析模型進(jìn)行了分析

> anova (fit)

Analysis of Variance Table

Response: value

Df Sum Sq Mean Sq F value Pr(>F)

Group4? ? 3? 1926? ? 642? ? 2.19? 0.14

Residuals 12? 3513? ? 293? ? ?

或者只使用以下簡(jiǎn)潔的R代碼：函數(shù)aov()嵌套在summary()函數(shù)中顿涣。

> summary(aov(value~Group4, data=df_CH_G))

Df Sum Sq Mean Sq F value Pr(>F)

Group4? ? ? 3? 1926? ? 642? ? 2.19? 0.14

Residuals? 12? 3513? ? 293

您可能還想使用下面的R代碼輸出截取結(jié)果。

> aov_fit <- aov(value*Group4,data=df_CH_G)

> summary(aov_fit, intercept=T)

該函數(shù)的輸出是一個(gè)經(jīng)典的方差分析表。

Df Sum Sq Mean Sq F value Pr(>F)

(Intercept)? 1? 83450? 83450? 285.08? 1e-09 ***

Group4? ? ? 3? 1926? ? 642? ? 2.19? 0.14? ?

Residuals? 12? 3513? ? 293? ? ? ? ? ? ? ? ?

---

Signif. codes:? 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

當(dāng)p值>0.05時(shí)揉阎，我們接受零假設(shè)H0：這四個(gè)均值沒有什么不同。您還可以將計(jì)算的F值與表格中的F值進(jìn)行比較：

> qf(0.95, 12, 3)

[1] 8.745

因?yàn)榱斜淼腇值大于計(jì)算的F值烙如，所以我們接受零點(diǎn)假設(shè)。方差分析的結(jié)果有點(diǎn)亂七八糟刀闷。您可以使用broom包來獲得整齊、信息更豐富的表格施蜜。

> library(broom)

> tidy(aov_fit)

term df sumsq meansq statistic p.value

1? ? Group4? 3? 1926? 641.9? ? 2.193? 0.1418

2 Residuals 12? 3513? 292.7? ? ? ? NA? ? ? NA

> augment(aov_fit)

? ? ? .rownames? ? ? ? ? ? ?value? ? ? Group4? ? ? ? ?.fitted? ?.se.fit? ? ?.resid? ? ? ? .hat? ?.sigma? ? ?.cooksd? ? ? ?.std.resid

1? ? 5_15_drySt-28F? ? 94.75? ? ?Fecal.Vdr-/-? ? 83.17? ? 7.651? 11.584? ? 0.2000? 17.44? ? 0.0358109? ? ?? 0.7570

2? ? 20_12_CeSt-28F? ?59.80? ? ?Cecal.Vdr-/-? ? 57.37? ?7.651? ?2.432? ? ?0.2000? ?17.85? ?0.0015781? ??? ?0.1589

3? ? 1_11_drySt-28F? ? ?77.00? ? ?Fecal.Vdr-/-? ? 83.17? ?7.651? ?-6.166? ? 0.2000? ?17.75? ?0.0101485? ? ?? -0.4030

4? ? 2_12_drySt-28F? ?103.27? ? ?Fecal.Vdr-/-? ? 83.17? ?7.651? ?20.106? ?0.2000? ?16.53? ?0.1078934? ? ? ?1.3139

5? ? 3_13_drySt-28F? ? 85.67? ? ? Fecal.Vdr-/-? ? 83.17? 7.651? ? 2.500? ? ?0.2000? ? 17.85? 0.0016683? ? ? ?0.1634

6? ? 4_14_drySt-28F? ? 55.14? ? ? Fecal.Vdr-/-? ? ?83.17? 7.651? -28.024? ?0.2000? ? 15.17? 0.2095941? ? ?-1.8313

7? ? 7_22_drySt-28F? ? 62.75? ? ? Fecal.WT? ? ? ? 70.31? ?9.878? ?-7.556? ?0.3333? ? ?17.65? ?0.0365658? ?- 0.5409

8? ? 8_23_drySt-28F? ? 67.67? ? ? Fecal.WT? ? ? ?70.31? ?9.878? ?-2.639? ? 0.3333? ? ?17.84? ?0.0044605? ? -0.1889

9? ? 9_24_drySt-28F? ? 80.50? ? ? Fecal.WT? ? ? 70.31? ? 9.878? ? 10.194? ?0.3333? ? 17.47? ?0.0665687? ? 0.7298

10? 19_11_CeSt-28F? ?52.17? ? ?Cecal.Vdr-/-? ? 57.37? ?7.651? ? -5.202? ?0.2000? ? ?17.78? ?0.0072213? ?-0.3399

11? 21_13_CeSt-28F? ?55.00? ? ?Cecal.Vdr-/-? ?57.37? ? 7.651? ? -2.368? ?0.2000? ? ? 17.85? ?0.0014970? ?-0.1548

12? 22_14_CeSt-28F? 59.00? ? ? Cecal.Vdr-/-? ? 57.37? ? 7.651? ?1.632? ?0.2000? ? ? 17.86? ?0.0007105? ? ?0.1066

13? 23_15_CeSt-28F? 60.88? ? ?Cecal.Vdr-/-? ? 57.37? ? 7.651? ? 3.507? 0.2000? ? ? ?17.83? ?0.0032819? ? ?0.2292

14? 25_22_CeSt-28F? 51.00? ? ? Cecal.WT? ? ? ?80.64? ?9.878? ?-29.639 0.3333? ? ? 14.13? ? 0.5626776? ??-2.1217

15? 26_23_CeSt-28F? 112.86? ? Cecal.WT? ? ? 80.64? ? ?9.878? ?32.218? 0.3333? ? ? 13.33? ? 0.6648948? ??2.3063

16? 27_24_CeSt-28F? 78.06? ? ?Cecal.WT? ? ? ?80.64? ? 9.878? ? -2.580? ?0.3333? ? ? 17.84? ? 0.0042631? ?-0.1847

> glance(aov_fit)

? ?r.squared? ?adj.r.squared? ?sigma? statistic? ?p.value? ?df? ?logLik? ? AIC? ? ?BIC? ? ?deviance? ??df.residual

1? ? 0.3541? ? ? ? 0.1926? ? ? ? ?17.11? ? 2.193? ? 0.1418? ? 4? ? -65.84? ?141.7? ?145.5? ? 3513? ? ? ? ? ?12

②兩兩比較和Tukey多重比較

方差分析的結(jié)果給出了組間差異的總體檢驗(yàn)(在這種情況下修械，4個(gè)組分別使用糞便、盲腸蹦渣、Vdr?/?和WT組合)。我們的目的也是測(cè)試與Chao-1豐富度相關(guān)的每對(duì)差異锄奢。下面的步驟是為了說明R中的成對(duì)t檢驗(yàn)和Tukey的即席多重比較的能力。讓我們對(duì)所有四組進(jìn)行未經(jīng)調(diào)整的成對(duì)t檢驗(yàn)堪滨。對(duì)于此測(cè)試遏乔，R中的默認(rèn)設(shè)置是使用Holm方法作為后期調(diào)整p值，因此要獲得未調(diào)整的p值，您需要指定p.adjust=“None”胞谭。R的默認(rèn)值是假設(shè)方差齊次丈屹，因此不需要指定池肤无。Sd=T竞漾。如果您的數(shù)據(jù)具有不等的方差坦仍，則需要使用pool.sd=F幔荒。

> #Pairwise tests of mean differences

> pairwise.t.test(df_CH_G$value, df_CH_G$Group4,p.adjust="none",pool.sd=T)

Pairwise comparisons using t tests with pooled SD

data:? df_CH_G$value and df_CH_G$Group4

Cecal.Vdr-/- Fecal.Vdr-/- Cecal.WT

Fecal.Vdr-/- 0.03? ? ? ? -? -

Cecal.WT? ? 0.09? ? ? 0.84? ? ? ? -

Fecal.WT? ? 0.32? ? ? ? 0.32? ? ? ? 0.47? ?

P value adjustment method: none

如果我們不對(duì)p值進(jìn)行任何調(diào)整，則糞便vdr?/?和盲腸vdr?/?之間存在統(tǒng)計(jì)學(xué)差異，盲腸wt和盲腸vdr?/?之間在統(tǒng)計(jì)學(xué)上略有差異硅确。這些不同之處在上面的框圖中可見一斑君编。正如我們注意到的，p.adjust()函數(shù)嵌套在pairwise.t.test()函數(shù)中兑燥。這是一個(gè)基本且非常有用的R函數(shù)。它可以用來控制家族I型誤差力崇。P.adjust()函數(shù)可以嵌套在其他函數(shù)中，也可以獨(dú)立調(diào)用茧吊。在獨(dú)立調(diào)用中，語法如下：p.adjust(p, method = p.adjust.methods, n = length(p))讶踪。?其中，p=p值的數(shù)值向量云石，method=校正方法，n=比較次數(shù)奖地，必須至少是長(zhǎng)度(P)。調(diào)整方式包括c(“Bonferroni”犬庇、“Holm”、“Hochberg”欢峰、“Hommel”、“BH”懊直、“by”、“FDR”、“None”)尝偎。其中“bonferroni”是將p值乘以比較次數(shù)的Bonferroni校正；“Holm”、“Hochberg”、“Hommel”寺晌、“BH”耘婚、“by”、“FDR”是指Holm(1979)、Hochberg(1988)兢榨、Hommel(1988)吵聪、Benjamini和Hochberg(1995)以及Benjamini和Yekutieli(2001)的別名，而“FDR”是指Holm(1979)、Hochberg(1988)局蚀、Hommel(1988)和Benjamini and Yekutieli(2001)的別名。它們是不那么保守的修正。

> #conservative Bonferroni adjustment

> pairwise.t.test(df_CH_G$value, df_CH_G$Group4, p.adjust="bonferroni", pool. sd = T)

Pairwise comparisons using t tests with pooled SD

data:? df_CH_G$value and df_CH_G$Group4

Cecal.Vdr-/- Fecal.Vdr-/- Cecal.WT

Fecal.Vdr-/- 0.2? ? ? ? ? -? -

Cecal.WT? ? 0.5? ? ? ? ? 1.0? ? ? ? ? -

Fecal.WT? ? 1.0? ? ? 1.0? ? ? ? ? 1.0? ?

P value adjustment method: bonferroni

> #Holm method

> pairwise.t.test(df_CH_G$value, df_CH_G$Group4, p.adjust="holm",pool.sd = T)

Pairwise comparisons using t tests with pooled SD

data:? df_CH_G$value and df_CH_G$Group4

? ? ? ? ? ? ? ? ? ? Cecal.Vdr-/- Fecal.Vdr-/- Cecal.WT

Fecal.Vdr-/- 0.2 - -

Cecal.WT? ? 0.4? ? ? ? ? 1.0? ? ? ? ? -

Fecal.WT? ? 1.0? ? ? ? ? 1.0? ? ? ? ? 1.0

P value adjustment method: holm

> #Benjamini & Hochberg(BH)

> pairwise.t.test(df_CH_G$value, df_CH_G$Group4, p.adjust="BH", pool.sd = T)

Pairwise comparisons using t tests with pooled SD

data:? df_CH_G$value and df_CH_G$Group4

Cecal.Vdr-/- Fecal.Vdr-/- Cecal.WT

Fecal.Vdr-/- 0.2? ? ? ? ? -? -

Cecal.WT? ? 0.3? ? ? ? ? 0.8? ? ? ? ? -

Fecal.WT? ? 0.5? ? ? ? ? 0.5? ? ? ? ? 0.6? ?

P value adjustment method: BH

> #Benjamini & Yekutieli

> pairwise.t.test(df_CH_G$value, df_CH_G$Group4, p.adjust="BY", pool.sd = T)

Pairwise comparisons using t tests with pooled SD

data:? df_CH_G$value and df_CH_G$Group4

Cecal.Vdr-/- Fecal.Vdr-/- Cecal.WT

Fecal.Vdr-/- 0.5? ? ? ? ? -? -

Cecal.WT? ? 0.6? ? ? ? ? 1.0? ? ? ? ? -

Fecal.WT? ? 1.0? ? ? ? ? 1.0? ? ? ? ? 1.0? ?

P value adjustment method: BY

上述所有四個(gè)調(diào)整都沒有給出兩兩比較的顯著差異。保守的Bonferroni和Benjamini&Yekutieli調(diào)整最大p值。用Benjamini&Hochberg方法比較也不顯著株茶，但調(diào)整后的p值較小。在這種情況下，Benjamini&Hochberg方法更為強(qiáng)大棍厂。Benjamini&Hochberg(BH)和Benjamini&Yekutieli(BY)方法都用于調(diào)整“錯(cuò)誤發(fā)現(xiàn)率”。實(shí)際上张漂，這并不是對(duì)家庭錯(cuò)誤的真正控制趴梢。錯(cuò)誤發(fā)現(xiàn)率方法得到相同的結(jié)果：所有的配對(duì)比較都沒有顯著差異憔狞。接下來簇抵，讓我們展示如何使用TukeyHSD()函數(shù)對(duì)平均值進(jìn)行Tukey多重比較，并獲得它們的置信區(qū)間焦履。調(diào)用此函數(shù)的方式類似于Summary()函數(shù)。它將原始方差分析計(jì)算中的變量作為其參數(shù)之一(圖8.6)屑宠。

> #Tukey multiple comparisons of means

> TukeyHSD(aov_fit, conf.level=.95)

Tukey multiple comparisons of means

95% family-wise confidence level

Fit: aov(formula = value ~ Group4, data = df_CH_G)

$Group4

diff? ? lwr? upr? p adj

Fecal.Vdr-/--Cecal.Vdr-/-? 25.798? -6.328 57.92 0.1334

Cecal.WT-Cecal.Vdr-/-? 23.270 -13.825 60.37 0.2935

Fecal.WT-Cecal.Vdr-/-? 12.937 -24.159 50.03 0.7328

Cecal.WT-Fecal.Vdr-/-? -2.528 -39.624 34.57 0.9969

Fecal.WT-Fecal.Vdr-/-? -12.861 -49.957 24.23 0.7362

Fecal.WT-Cecal.WT? ? ? ? -10.333 -51.807 31.14 0.8792

> plot(TukeyHSD(aov(df_CH_G$value~df_CH_G$Group4), conf.level=.95))

此圖表示所有可能的成對(duì)測(cè)試和p值，以及95%的置信區(qū)間仇让〉浞睿可以根據(jù)您的選擇更改默認(rèn)的95%置信級(jí)別。由于所有置信線均為0丧叽，因此對(duì)于此示例卫玖，使用Tukey多重比較進(jìn)行調(diào)整后沒有顯著不同的項(xiàng)。

8.4 Comparisons Among More than Two Groups Using Kruskal-Wallis Test（ 使用Kruskal-Wallis檢驗(yàn)對(duì)兩組以上人群進(jìn)行比較）

①?Kruskal-Wallis Test

以William Kruskal和W.Allen? Wallis命名的Kruskal-Wallis檢驗(yàn)或單因素秩和方差分析是檢驗(yàn)樣本是否來自同一分布的非參數(shù)方法溉躲。Kruskal-Wallis檢驗(yàn)的參數(shù)等價(jià)性是單向方差分析。它將Mann-Whitney U檢驗(yàn)擴(kuò)展到兩組以上滔以。Kruskal-Wallis檢驗(yàn)的零假設(shè)是各組的平均等級(jí)相同。與類似的單因素方差分析不同罚屋，非參數(shù)Kruskal-Wallis檢驗(yàn)不假定基礎(chǔ)數(shù)據(jù)為正態(tài)分布漆弄。它已被廣泛應(yīng)用于微生物組的研究渤愁。例如，測(cè)序后的微生物組數(shù)據(jù)不是正態(tài)分布的势篡，并且包含一些很強(qiáng)的異常值祝钢。因此，使用秩而不是實(shí)際值是合適的显晶，以避免測(cè)試受到異常值的存在或非正態(tài)分布的影響。測(cè)試統(tǒng)計(jì)需要四個(gè)主要步驟：

步驟1.將所有組中的所有數(shù)據(jù)一起按升序排列在單個(gè)序列中枉疼，即忽略組成員身份從1到N進(jìn)行排序侮措。

步驟2.通過平均排名位置來分配任何并列的值。

步驟3.總結(jié)不同的等級(jí)，例如R1R2R3…。對(duì)于每個(gè)不同的組呀狼。

步驟4.應(yīng)用以下公式計(jì)算測(cè)試統(tǒng)計(jì)量：

Kruskal-Wallis檢驗(yàn)統(tǒng)計(jì)量近似為卡方分布梆暖，如果n值“大”芒划，k?為1個(gè)自由度瘫析。當(dāng)每個(gè)數(shù)值大于或等于5時(shí)，該近似通常被認(rèn)為是足夠的呜达。

②Compare Diversities Among Groups（比較不同群體之間的差異）

Kruskal-Wallis檢驗(yàn)或Kruskal-Wallis單因素方差分析用于比較數(shù)據(jù)不服從正態(tài)分布的多組胎署。該檢驗(yàn)類似于兩個(gè)樣本的Wilcoxon秩和檢驗(yàn)占调。我們首先使用我們的VDR?/?鼠標(biāo)數(shù)據(jù)來說明此測(cè)試幔虏。

> library(dplyr)

> Data <- mutate(df_CH_G, Group = factor(df_CH_G$Group4, levels=unique (df_CH_G$Group4)))

獲取描述性統(tǒng)計(jì)信息

> library(FSA)

> Summarize(value ~ Group4, data = df_CH_G)

Group4 n? mean? ? sd? min? ? Q1 median? ? Q3? ? max

1 Cecal.Vdr-/- 5 57.37? 3.659 52.17 55.00? 59.00 59.80? 60.88

2 Fecal.Vdr-/- 5 83.17 18.494 55.14 77.00? 85.67 94.75 103.27

3? ? Cecal.WT 3 80.64 31.009 51.00 64.53? 78.06 95.46 112.86

4? ? Fecal.WT 3 70.31? 9.165 62.75 65.21? 67.67 74.08? 80.50

按組生成直方圖 (Fig. 8.7)。

> #Individual plots in panel of 2 columns and 2 rows

> library(lattice)

> histogram(* value|Group4, data=df_CH_G,layout=c(2,2))

直方圖顯示狼讨，在這種情況下，組之間的值分布是不同的∑饩海現(xiàn)在政供，我們使用kruskal.test()函數(shù)進(jìn)行Kruskal-Wallis檢驗(yàn)來比較中值的差異。

> #kruskal wallis test of Chao 1 richness

> kruskal.test(value * Group4, data = df_CH_G)

Kruskal-Wallis rank sum test

data: value by Group4

Kruskal-Wallis chi-squared = 5.2, df = 3, p-value = 0.2

檢驗(yàn)統(tǒng)計(jì)值為5.2朽基，p值大于0.05布隔，也低于卡方表法：

> qchisq(0.950, 3)

[1] 7.815

因此，我們接受零假設(shè)H0：在5%顯著水平上稼虎，4組的中位數(shù)在統(tǒng)計(jì)上是相等的衅檀。一般情況下，如果Kruskal-Wallis檢驗(yàn)有意義霎俩，則進(jìn)一步進(jìn)行事后分析哀军，以找出哪些組的水平彼此不同。在這種情況下打却，Kruskal-Wallis測(cè)試并不重要杉适。為了說明這一點(diǎn)，我們進(jìn)行了兩個(gè)事后檢驗(yàn)：Nemenyi檢驗(yàn)和Dunn檢驗(yàn)柳击。與方差分析類似猿推，我們可以選擇調(diào)整p值的方法來控制家族錯(cuò)誤率或控制錯(cuò)誤發(fā)現(xiàn)率。當(dāng)您在R或RStudio中輸入捌肴？p.adjust時(shí)蹬叭，它會(huì)顯示一個(gè)指向文檔“調(diào)整多個(gè)比較的P值”的鏈接。您可以通過此鏈接查看調(diào)整方法的詳細(xì)信息状知。

多重比較的Nemenyi檢驗(yàn)：Nemenyi測(cè)試通過DescTools包中的Nemenyi Test()函數(shù)執(zhí)行秽五。我們首先加載DescTools包并調(diào)用函數(shù)Nemenyi Test()。這個(gè)調(diào)整p值的方法應(yīng)該是“tukey”试幽、“ChiSq”之一筝蚕。這里我們選擇Tukey方法卦碾。

> library(DescTools)

> #Tukey method for adjusting p-values

> Test_N = NemenyiTest(x = df_CH_G$value,

+? ? ? ? ? ? ? ? ? g = df_CH_G$Group4,

+? ? ? ? ? ? ? ? ? dist="tukey")

> Test_N

Nemenyi's test of multiple comparisons for independent samples (tukey)?

mean.rank.diff? pval? ?

Fecal.Vdr-/--Cecal.Vdr-/-? 6.6000 0.1254? ?

Cecal.WT-Cecal.Vdr-/-? 4.7333 0.5237? ?

Fecal.WT-Cecal.Vdr-/-? 5.0667 0.4636? ?

Cecal.WT-Fecal.Vdr-/-? -1.8667 0.9501? ?

Fecal.WT-Fecal.Vdr-/-? -1.5333 0.9713? ?

Fecal.WT-Cecal.WT? ? ? ? ? ? ? ? 0.3333 0.9998? ?

---

Signif. codes:? 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Nemenyi‘s檢驗(yàn)表明，采用Tukey調(diào)整法起宽，在糞便和盲腸VDR基因敲除樣本中洲胖，不同地點(diǎn)和不同基因型的CHO1多樣性的平均秩差異不顯著。然而坯沪，當(dāng)兩組的觀察次數(shù)不相等時(shí)绿映，Nemenyi檢驗(yàn)是不合適的，而Dunn檢驗(yàn)是合適的(Zar 2010)腐晾。我們運(yùn)行Dunn檢驗(yàn)叉弦，如下所示。

多重比較的Dunn檢驗(yàn)：最流行的后Kruskal-Wallis檢驗(yàn)是Dunn檢驗(yàn)藻糖。我們可以使用FSA包中的dunnTest()函數(shù)執(zhí)行Dunn測(cè)試淹冰。下面，我們調(diào)用函數(shù)dunnTest()并使用Benjamini和Hochberg方法調(diào)整p值巨柒。

> library(FSA)

> # “bh” suggests Benjamini and Hochberg? method for adjusting p-values

> Test_N = dunnTest(df_CH_G$value ~ df_CH_G$Group4,data=df_CH_G, method="bh")

> Test_N

Dunn (1964) Kruskal-Wallis multiple comparison

p-values adjusted with the Benjamini-Hochberg method.

Comparison? ? ? ? Z P.unadj? P.adj

1? ? Cecal.Vdr-/- - Cecal.WT -1.36136 0.17340 0.3468

2 Cecal.Vdr-/- - Fecal.Vdr-/- -2.19190 0.02839 0.1703

3? ? Cecal.WT - Fecal.Vdr-/- -0.53688 0.59135 0.8870

4? ? Cecal.Vdr-/- - Fecal.WT -1.45723 0.14505 0.4352

5? ? ? ? Cecal.WT - Fecal.WT -0.08575 0.93167 0.9317

6? ? Fecal.Vdr-/- - Fecal.WT? 0.44100 0.65921 0.7911

Dunn檢驗(yàn)表明樱拴，盲腸和糞便中VdR?/?樣品的CHO1多樣性差異有統(tǒng)計(jì)學(xué)意義。但是洋满，用Benjamini-Hochberg方法調(diào)整多重比較p值后晶乔，樣本中的位置和基因型之間沒有統(tǒng)計(jì)學(xué)意義。

③ Find Significant Taxa Among Groups（在群組中找到重要的分類群）

在本節(jié)中牺勾，我們使用Kruskal-Wallis檢驗(yàn)來說明如何在組之間尋找重要的分類群正罢。假設(shè)我們想知道在來自糞便和盲腸的Vdr、?/?和WT小鼠樣本中是否存在任何重要的分類群驻民。我們使用Kruskal-Wallis測(cè)試對(duì)數(shù)據(jù)集中的248個(gè)分類(細(xì)菌)中的每一個(gè)進(jìn)行測(cè)試翻具。首先，對(duì)豐富的數(shù)據(jù)進(jìn)行歸一化回还，并將數(shù)據(jù)制成數(shù)據(jù)框呛占。一種標(biāo)準(zhǔn)化方法是使用對(duì)數(shù)變換。

> data<-log((abund_table+1)/(rowSums(abund_table)+dim(abund_table) [2]))

> df<-as.data.frame(data)

另一種歸一化方法是將當(dāng)量計(jì)數(shù)轉(zhuǎn)換為相對(duì)豐度懦趋。

> df <- as.data.frame(abund_table/rowSums(abund_table))

然后，使用kruskal.test()函數(shù)和迭代R函數(shù)執(zhí)行248個(gè)測(cè)試(每個(gè)測(cè)試針對(duì)一個(gè)細(xì)菌)疹味。Kruskal.test()函數(shù)有幾個(gè)關(guān)鍵組件：代碼為“for(i in 1：Dim(Df)[2])”的所有分類群(柱)的測(cè)試都是循環(huán)進(jìn)行的仅叫。·對(duì)于每個(gè)循環(huán)糙捺，使用代碼“kw_test<-kruskal.test(df[诫咱，i]，g=Group4)”運(yùn)行Kruskal-Wallis測(cè)試洪灯】茬裕·結(jié)果存儲(chǔ)在數(shù)據(jù)框中，每個(gè)樣本一行，KW測(cè)試的每個(gè)p值一列掏呼』悼欤·報(bào)告cat函數(shù)為“cat(Paste(”Kruskal-Wallis test for“，Names(Df)[i]憎夷，”“莽鸿，i，”/“拾给，dim(Df)[2]祥得，”；p-value=“蒋得，kw_test$p.value级及，”\n“，sep=”“)”的測(cè)試次數(shù)

> KW_table <- data.frame()

> for (i in 1:dim(df)[2]) {

+? #run KW test for each bacterium

+? KW_test <- kruskal.test(df[,i], g=df_CH_G$Group4)

+? # Store the result in the data frame

+? KW_table <- rbind(KW_table,

+? data.frame(id=names(df)[i],

+? p.value=KW_test$p.value

+ )? )

+? # Report number of bacteria tested

+ cat(paste("Kruskal-Wallis test for ",names(df)[i]," ", i, "/",

+? dim(df)[2], "; p-value=", KW_test$p.value,"\n", sep=""))

+ }

檢查數(shù)據(jù)框表以確保函數(shù)工作正常：

> #Check the data frame table

> head(kW_table)

id? p.value

1? ? ? ? ? ? ? ? Tannerella 0.005289

2? ? ? ? ? ? ? ? Lactococcus 0.407302

3? ? ? ? ? ? ? Lactobacillus 0.058626

4 Lactobacillus::Lactococcus 0.476355

5? ? ? ? ? ? Parasutterella 0.120519

6? ? ? ? ? ? ? Helicobacter 0.140053

④多重測(cè)試與E值额衙、FWER和FDR

文獻(xiàn)中存在幾種不同類型的多重測(cè)試校正饮焦。其中，Bonferroni修正較為保守入偷。修正的方法很簡(jiǎn)單追驴，就是用阿爾法除以測(cè)試次數(shù)。我們已經(jīng)說明了使用Bonferroni疏之，Holm殿雪，Tukey方法和兩種版本的方法調(diào)整成兩個(gè)版本的“假發(fā)現(xiàn)率”的p值調(diào)整，在第二節(jié)中使用ANOVA進(jìn)行配對(duì)比較和檢驗(yàn)锋爪。8.3.2丙曙，Kruskal-Wallis檢驗(yàn)在Sect.Kruskal-Wallis檢驗(yàn)中的帖子。8.4.2.其骄。在本節(jié)中亏镰，我們介紹一般的多重測(cè)試修正：E值、系列誤差率(FWER)和FDR拯爽。

E-value：E值是您進(jìn)行多次測(cè)試時(shí)偶然出現(xiàn)的錯(cuò)誤陽性的預(yù)期數(shù)量索抓。您可以簡(jiǎn)單地將p值與對(duì)其執(zhí)行測(cè)試的分類群數(shù)相乘就可以得到它：E-value = p-value X the number of tests（測(cè)試次數(shù)）。請(qǐng)注意毯炮，在E值上逼肯，基本校正是使用原始α，即來自測(cè)試的p值桃煎，而不是標(biāo)稱p值篮幢。

> #E-value

> KW_table$E.value <-KW_table$p.value * dim(KW_table)[1]

> KW_table$E.value

由于E值只是將p值乘以測(cè)試次數(shù)，因此它可以大于1为迈。如果數(shù)據(jù)框中有許多分類群用于測(cè)試三椿，則此校正方法不容易找到有意義的分類群缺菌。重要分類群是E值遠(yuǎn)小于1的分類群。以下代碼用于檢查是否已將E值添加到結(jié)果數(shù)據(jù)框中：

> #check E-value in result data frame

> head(kW_table)

id? p.value E.value

1? ? ? ? ? ? ? ? Tannerella 0.005289? 1.312

2? ? ? ? ? ? ? ? Lactococcus 0.407302 101.011

3? ? ? ? ? ? ? Lactobacillus 0.058626? 14.539

4 Lactobacillus::Lactococcus 0.476355 118.136

5? ? ? ? ? ? Parasutterella 0.120519? 29.889

6? ? ? ? ? ? ? Helicobacter 0.140053? 34.733

FWER：FWER是至少出現(xiàn)一次誤報(bào)(I類錯(cuò)誤)的概率搜锰。換句話說伴郁，這是您沒有拒絕零假設(shè)H0的概率：在進(jìn)行多個(gè)測(cè)試時(shí)，不同組之間沒有差異纽乱。該公式由下式給出

式中蛾绎，T=測(cè)試次數(shù)。在R中鸦列，為避免直接使用上述公式計(jì)算FWER帶來的舍入誤差租冠，最好采用右尾二項(xiàng)分布檢驗(yàn)來計(jì)算FWER。R碼如下所示薯嗤。

> #FWER

> KW_table$FWER <- pbinom(q=0, p=KW_table$p.value,size=dim(KW_table)[1], lower.tail=FALSE)

檢查數(shù)據(jù)框以查看是否將FWER添加到結(jié)果數(shù)據(jù)框：

> #check the dataframetable

> head(kW_table)

id? p.value E.value? FWER

1? ? ? ? ? ? ? ? Tannerella 0.005289? 1.312 0.7316

2? ? ? ? ? ? ? ? Lactococcus 0.407302 101.011 1.0000

3? ? ? ? ? ? ? Lactobacillus 0.058626? 14.539 1.0000

4 Lactobacillus::Lactococcus 0.476355 118.136 1.0000

5? ? ? ? ? ? Parasutterella 0.120519? 29.889 1.0000

6? ? ? ? ? ? ? Helicobacter 0.140053? 34.733 1.0000

FDR：最后但并非最不重要的是FDR顽爹。Benjamini和Hochberg(1995)將錯(cuò)誤發(fā)現(xiàn)率定義如下：FDR=expected proportion of erroneous rejections among all rejections（錯(cuò)誤拒絕在所有拒絕中的預(yù)期比例）。在本例中骆姐，F(xiàn)DR是在進(jìn)行多次測(cè)試時(shí)被接受為陽性的分類群中假陽性的比例镜粤。 Benjamini-Hochberg校正包括以下步驟：首先，從最小到最大對(duì)p值進(jìn)行排序玻褪，并進(jìn)行排序(1肉渴，2，3带射，…同规。，k窟社，…券勺。，T)灿里；代碼如下：

> #FDR

> #order p-values from smallest to largest

> kW_table <- kW_table[order(kW_table$p.value, decreasing=FALSE), ]

> head(kW_table)

id? p.value E.value? FWER

8? ? Bacteroides 0.003976? 0.986 0.6277

19? ? Alistipes 0.004637? 1.150 0.6842

7? ? Prevotella 0.005174? 1.283 0.7238

39 Butyricimonas 0.005174? 1.283 0.7238

1? ? Tannerella 0.005289? 1.312 0.7316

10? Odoribacter 0.008189? 2.031 0.8699

接下來关炼，使用以下公式和代碼計(jì)算q值：q-value = p - value *T/k;

> #calculate q-value

> kW_table$q.value.factor <- dim(kW_table)[1] / 1:dim(kW_table)[1]

> head(kW_table$q.value.factor)

[1] 248.00 124.00? 82.67? 62.00? 49.60? 41.33

> kW_table$q.value <- kW_table$p.value * kW_table$q.value.factor

> head(kW_table$q.value)

[1] 0.9860 0.5749 0.4277 0.3208 0.2623 0.3385

> #check to see if q-value added to the result data frame

> head(kW_table)

id? p.value E.value? FWER q.value.factor q.value

8? ? Bacteroides 0.003976? 0.986 0.6277? ? ? ? 248.00? 0.9860

19? ? Alistipes 0.004637? 1.150 0.6842? ? ? ? 124.00? 0.5749

7? ? Prevotella 0.005174? 1.283 0.7238? ? ? ? ? 82.67? 0.4277

39 Butyricimonas 0.005174? 1.283 0.7238? ? ? ? ? 62.00? 0.3208

1? ? Tannerella 0.005289? 1.312 0.7316? ? ? ? ? 49.60? 0.2623

10? Odoribacter 0.008189? 2.031 0.8699? ? ? ? ? 41.33? 0.3385

然后，通過使用以下代碼指定目標(biāo)FDR匣吊，并識(shí)別q值等于或小于指定阿爾法的排名列表的最后一項(xiàng)：

> #set up alpha value

> KW_alpha=0.05

> #identify the last item of the ranked list with a q-value =< alpha

> last.significant.item <- max(which(KW_table$q.value <= KW_alpha))

Warning message:

In max(which(KW_table$q.value <= KW_alpha)) :

no non-missing arguments to max; returning -Inf

> last.significant.item

[1] -Inf

在我們的示例中儒拂，沒有小于或等于指定的alpha的q值，因此程序返回負(fù)無限色鸳。最后侣灶，顯示結(jié)果框表和選擇的分類群：

> #display the chosen results

> selected <- 1:5

> #selected <- 1:last.significant.item

> print(kW_table[selected,])

? ? ? ? ? ?id? ? ? ? ? ? ? ?p.value E.value? FWER q.value.factor q.value

8? ? Bacteroides 0.003976? 0.986 0.6277? ? ? ? 248.00? 0.9860

19? ? Alistipes 0.004637? 1.150 0.6842? ? ? ? 124.00? 0.5749

7? Prevotella 0.005174? 1.283 0.7238? ? ? ? ? 82.67? 0.4277

39 Butyricimonas 0.005174? 1.283 0.7238? ? ? ? ? 62.00? 0.3208

1? ? Tannerella 0.005289? 1.312 0.7316? ? ? ? ? 49.60? 0.2623

> diff.taxa.factor <- kW_table$id[selected]

> diff.taxa <- as.vector(diff.taxa.factor)

> diff.taxa

[1] "Bacteroides"? "Alistipes"? ? "Prevotella"? ?"Butyricimonas" "Tannerella"

因?yàn)樵谶@種情況下沒有小于或等于指定的alpha=0.05的Q值，所以上面顯示的5個(gè)分類群不是基于FDR的缕碎。它們是p值最小的選定分類群。本賈米尼-霍希伯格校正比上面提出的其他多次測(cè)試校正不那么嚴(yán)格池户，因此具有更高的靈敏度咏雌。FDR廣泛應(yīng)用于微生物群系和其他研究領(lǐng)域和許多R函數(shù)凡怎。

8.5 Summary

在這一章中，我們介紹了各種研究領(lǐng)域中常用的和經(jīng)典的方法赊抖。其中一些在微生物組研究中得到了廣泛的應(yīng)用统倒。我們用分步實(shí)施的方式說明了這些微生物組數(shù)據(jù)的分析方法。R系統(tǒng)中的狀態(tài)氛雪。數(shù)據(jù)集來自我們自己的出版物(金等人)房匆。2015年；王等人报亩。2016)浴鸿。讀者可以使用本章提供的R代碼和解釋來分析自己的微生物組數(shù)據(jù)。在本章中弦追，我們重點(diǎn)研究了單變量群落微生物組數(shù)據(jù)的假設(shè)檢驗(yàn)岳链。在接下來的一章里。將著重對(duì)多變量群落微生物組數(shù)據(jù)進(jìn)行假設(shè)檢驗(yàn)劲件。