我们将微生物群落组成研究分为两个主要部分：(1)分类多样性、OTUS和类群的假设检验；(2)类群间差异分析。第一个成分主要属于单变量群落分析。第二个成分可以进一步分为各种多变量技术，如聚类和排序，以及多变量差异分析的假设检验。
8.1 两组间的差异性比较
在我们的vdr−/−小鼠研究中，目的之一是测试两组(vdr−/−小鼠和野生型小鼠)在粪便和盲肠部位的差异。第6章利用粪便样品计算了香农多样性。这里，为了说明单变量群落分析，将使用各种测试统计数据比较计算出的香农多样性。
①Two-Sample Welch’s t-Test
T统计量是由William Sealy Gosset在1908年提出的。两样本t检验用于检验两总体均值相等。当测试统计数据服从正态分布时，它最常用。如果两组具有相同的方差，则t统计量可按如下方式计算
首先，加载和转置数据集。
> abund_table=read.csv("VdrGenusCounts.csv",row.names=1, check.names=FALSE)
> abund_table<-t(abund_table)
为了将数据集中的分组信息直接合并到比较中，我们需要进行数据管理。在数据集中，样本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章计算了香农多样性。为了方便起见，这里重复一遍。
> library(vegan)
> H<-diversity(abund_table, "shannon")
制作了香农多样性的数据框。
> df_H<-data.frame(sample=names(H),value=H,measure=rep("Shannon",length(H)))
然后，我们将数据框多样性和分组相结合，形成一个新的数据框。
> 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
接下来，我们将新数据框中的粪便数据子集提出
> 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
现在数据已经准备好进行统计分析了。在进行假设检验之前，让我们使用ggplot2包中的函数gglot()和plyr包中的ddply()函数来研究Shannon多样性值的分布。
> library(ggplot2)
我们使用facet_grid将绘图拆分成两个面板。
> p<-ggplot(Fecal_G, aes(x=value))+
+ geom_histogram(color="black", fill="black")+
+ facet_grid(Group * .)
plyr包用于计算每个组的平均香农多样性值。
> 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基因敲除导致更高的分集，因为该组的直方图相对于WT组向右移动(更高的分集值)。为了检验香农多样性无差异的零假设，使用Welch的t检验，结果是p值=0.01(t=3.6，df=5.9)。因此，我们拒绝无差异的零假设，而支持香农多样性在两组中不同的替代方案。
> 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秩和检验）
Wilcoxon秩和检验等价于Mann和Whitney(1947)提出的Mann-Whitney U检验。它是两样本t检验的一种非参数替代，它使用两个独立样本数据的秩来检验零假设：两个独立样本来自具有相同分布的总体(即，两个总体相同)。与t-检验不同的是，Wilcoxon秩和检验不需要正态分布的假设，并且几乎与t-检验一样有效。因此，它在微生物组研究中得到了广泛的应用。执行Wilcoxon秩和检验以求出检验统计量的值需要三个主要步骤：
步骤1. 给所有的观察值分配排名，最小值的排名为1。如果值是平分的，则分配平局中涉及的排名的平均值。
步骤2. 对两个样本中任何一个样本的秩进行求和。可以确定另一样本中的秩和，因为所有秩的和等于N(N+1)/2，其中N是观察的总数。如果两个测试总体具有相同的分布，则秩R的平均值为 ，标准偏差为 。当秩和R远离均值时，Wilcoxon秩和检验拒绝两个总体具有相同分布的假设。随着两个样本量的增加，秩和统计量变得近似正常。我们可以通过标准化秩和来形成统计量。
步骤3.使用以下公式计算z检验统计量的值：
R 编号为n1的样本的秩和；n1 找到秩和R的样本大小(如样本1)；n2 另一个样本大小(如样本2)。
以下代码用于进行Wilcoxon秩和检验。
> 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显示样本大小较小时分布是不对称的。Wilcoxon秩和检验可能更合适，但由Wilcoxon秩和检验给出的p=0.04与Welcht检验在0.05显著性水平上p=0.01得出的结论相同。
8.2 Comparisons of a Taxon of Interest Between Two Groups
①用Wilcoxon秩和检验比较相对丰度
当我们分析微生物组丰度数据时，从样品中OTUs的丰度或类群的丰度来推断生态系统的总丰度是不合适的。取而代之的是，我们可以利用样本中的相对丰度来推断其在生态系统中的分类单元的相对丰度。其背后的原因是它存在一个组成约束：样本和中的所有微生物相对丰度为1，这导致组成数据驻留在单纯形(Aitchison 1982,1986)而不是欧几里德空间。因此，通常将数据标准化到一个共同的尺度，以便于跨组的分类单元丰度的比较。方法是将分类单元计数除以读取总次数100，以将丰度转换为样本中读取的百分比，将数据缩放为“每100次读取的分类单元数量”。当我们选择一个特定的单个分类单元进行跨组测试时，重要的是要确保指定的分类单元是基于假设或理论的，以减少夸大假阳性率的机会(即，在不应该拒绝零假设的时候拒绝它)。小鼠体内的VDR显著影响β多样性，并持续影响单个细菌分类群，如副细菌分类群(Wang等人。2016)。在本节中，我们使用粪便样本说明Wilcoxon秩和检验来比较VDR小鼠数据集中的细菌类杆菌。
首先，检查每个样本的总丰度。
> 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
然后，通过将每个值除以样本总丰度来计算相对丰度。
> relative_abund_table<-decostand(abund_table, method = "total")
检查每个样品的总丰度，使上述计算正确。
> 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
请看一下转换后的数据。
> 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
我们感兴趣的细菌类杆菌在第8栏。让我们对这个细菌进行分类。
> (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
现在，将类杆菌和分组数据框组合在一起，并对粪便样本进行子集以供以后使用。
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()函数用于生成带有组的简单类杆菌箱线图(图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")
以下代码用于使用函数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检验表明，Vdr−/−小鼠和WT小鼠的Bacteriodes相对丰度有统计学意义。我们可以得出结论，VDR基因敲除富集了Bacteriodes。
②有无紫杉醇的卡方检验比较（Chi-Square Test）
卡方检验，也被写成X2 test，通常被用作皮尔逊卡方检验的缩写，被提出，并由卡尔皮尔逊在1900年(皮尔逊1900)首次研究其性质。X2test应用于分类数据集，以检验观察到的频率分布是否与声称或理论分布不同(t拟合优度)，并调查列联表中的行变量和列变量是否彼此独立(独立性检验)。拟合优度的测试统计量如下所示：
作为经验法则，它要求所有预期的细胞计数等于或超过5，以提供对卡方分布的充分近似(Wackerly等人。2002)，尽管Cochran(1952)指出，在某些情况下，该值可能低至1。在本节中，我们使用盲肠样本说明X2test来比较VDR小鼠数据集中的细菌Parabacteroides。为了说明X2检验，我们将Parabacteroides的丰度计数数据转换为二进制变量。如果样本中没有分类单元，则丰度表中的Parabacteroides分类单元的计数数据将转换为0，或者如果样本中存在分类单元，则转换为1。表8.1汇总了转换后的数据。
首先，查看丰富的数据来鉴定细菌副细菌，并对该细菌进行亚群划分。
> abund_table[1:16,1:27]
> (Parabacteroides <- abund_table[,27])
然后，将子集数据与分组数据框合并。
> Parabacteroides_G <-cbind(Parabacteroides, grouping)
> rownames(Parabacteroides_G)<-NULL
由于组合数据框既包括粪便样本又包括盲肠样本，因此让我们从该数据框中将盲肠数据子集提出。
> 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
重新编码用于卡方检验的二进制变量“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
以下代码用于进行卡方检验。
> 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个vdr−/−盲肠样本中有3个(60%)有Parabacteroides，而3个野生型样本中有3个(100%)有Parabacteroides。对两组发病率无差异的零假设进行卡方检验，P值为0.7(X-平方=0.18，df=1)，不能排除两组发病率无差异的零假设，得出两组发病率无差异的结论。请注意，由于样本量较小，输出中会出现一条警告消息。通常，如果列联表中的单元格值很小(例如<5)，卡方检验可能不正确，则应用Fisher精确检验。
> 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精确检验结果也不显著，p值为0.5，与卡方检验结果一致。然而，这项测试很难在无限的置信区间内得出结论。
8.3 Comparisons Among More than Two Groups Using ANOVA（ 两组以上组间比较的方差分析）
①One-Way ANOVA（单因素方差分析）
方差分析(ANOVA)是由罗纳德·费希尔(Ronald Fisher，1918)于1918年提出的，1925年费舍尔的著作“研究人员的统计方法”出版后，方差分析变得广为人知。方差分析将两样本t-test推广到两组以上。方差分析的零假设是：比较组的所有均值相等。使用方差分析的分析依赖于基础数据的正态性假设。然而，大多数微生物群落组成数据，特别是多变量数据并不是正态分布的，因此，本书仅使用方差分析比较α多样性测度的单变量分析。对于多变量群落组成数据，要么采用非参数形式的方差分析，要么采用其他合适的统计方法。检验统计量的形成是通过传统的平方和分割(方差分割)来实现的。样本方差的定义方程为
总偏差因素比较采用F检验。F值是通过将处理之间的方差除以处理内的方差来获得的。下面给出了单因素方差分析中的F检验统计量
为了说明方差分析在微生物群落组成研究中的作用，我们使用了我们的VDR−/−小鼠数据。本研究的一个假设是VDR状态和肠道位置对肠道细菌群落没有影响。我们使用方差分析来分析 Chao 1α多样性测定来解决这一假设。复制了第6章中使用的代码下面的Chao 1丰富度测定。结果在一个数据框中。合并用于方差分析检验的组信息。下面的代码构成了一个chao1丰富性的数据框，并将组信息添加到该数据框中。
> 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
新的四个级别的组是通过位置和组的交互作用生成的。
> 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指数(图8.4)。
> boxplot(value*Group4, data=df_CH_G, col=rainbow(4), main="Chao1 index")
下面的ggplot()生成高质量的框图以供发表使用(图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)
除了目视检查底层数据的正态性外，还可以测试方差的均匀性。Sokal和Rohlf(2011)描述了三个这样的检验：Bartlett的同质性检验、Hartley的Fmaxtest和log-anova检验，或SchefféBox检验。为了继续使用方差分析进行验证，我们必须首先检验方差的均匀性。软件R提供两个测试：Bartlett测试和Fligner-Killeen测试。为了说明方差的均匀性检验，我们使用了来自粪便和盲肠位置的Vdr−/−和WT小鼠数据的Chao 1丰富度度量。零假设(H0)是四组中的所有方差都是相同的。我们从Bartlett测验开始。为了方便处理Bartlett检验，我们使用dplyr包中select()函数来选择相关的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代码执行方差的同质性的Bartlett 检验：
> 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
该函数给出了统计检验的K平方值和p值。它表明，在5%的水平上可以拒绝零假设。或者，我们可以将Bartlett的K平方与卡方表格的值进行比较，在qchisq()函数中使用相同级别的α和自由度。如果X平方>Bartlett‘sK平方，我们接受零假设H0(方差齐性)，否则拒绝零假设。因为卡方小于Bartlett的K平方，所以我们拒绝零假设H0，并得出方差不相同的结论。我们现在使用Fligner-Killeen检验来检验同方差。如下所示的语法非常相似。
> 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
结论与Bartlett检验相似：方差不同。然而，为了说明的目的，我们继续对数据进行方差分析，而不考虑方差的齐性检验结果。以下R码符合该模型：
> fit = lm(formula = value*Group4,data=df_CH_G)
然后对方差分析模型进行了分析
> 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
或者只使用以下简洁的R代码：函数aov()嵌套在summary()函数中。
> 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代码输出截取结果。
> aov_fit <- aov(value*Group4,data=df_CH_G)
> summary(aov_fit, intercept=T)
该函数的输出是一个经典的方差分析表。
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
当p值>0.05时，我们接受零假设H0：这四个均值没有什么不同。您还可以将计算的F值与表格中的F值进行比较：
> qf(0.95, 12, 3)
[1] 8.745
因为列表的F值大于计算的F值，所以我们接受零点假设。方差分析的结果有点乱七八糟。您可以使用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多重比较
方差分析的结果给出了组间差异的总体检验(在这种情况下，4个组分别使用粪便、盲肠、Vdr−/−和WT组合)。我们的目的也是测试与Chao-1丰富度相关的每对差异。下面的步骤是为了说明R中的成对t检验和Tukey的即席多重比较的能力。让我们对所有四组进行未经调整的成对t检验。对于此测试，R中的默认设置是使用Holm方法作为后期调整p值，因此要获得未调整的p值，您需要指定p.adjust=“None”。R的默认值是假设方差齐次，因此不需要指定池。Sd=T。如果您的数据具有不等的方差，则需要使用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
如果我们不对p值进行任何调整，则粪便vdr−/−和盲肠vdr−/−之间存在统计学差异，盲肠wt和盲肠vdr−/−之间在统计学上略有差异。这些不同之处在上面的框图中可见一斑。正如我们注意到的，p.adjust()函数嵌套在pairwise.t.test()函数中。这是一个基本且非常有用的R函数。它可以用来控制家族I型误差。P.adjust()函数可以嵌套在其他函数中，也可以独立调用。在独立调用中，语法如下：p.adjust(p, method = p.adjust.methods, n = length(p))。 其中，p=p值的数值向量，method=校正方法，n=比较次数，必须至少是长度(P)。调整方式包括c(“Bonferroni”、“Holm”、“Hochberg”、“Hommel”、“BH”、“by”、“FDR”、“None”)。其中“bonferroni”是将p值乘以比较次数的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
上述所有四个调整都没有给出两两比较的显著差异。保守的Bonferroni和Benjamini&Yekutieli调整最大p值。用Benjamini&Hochberg方法比较也不显著，但调整后的p值较小。在这种情况下，Benjamini&Hochberg方法更为强大。Benjamini&Hochberg(BH)和Benjamini&Yekutieli(BY)方法都用于调整“错误发现率”。实际上，这并不是对家庭错误的真正控制。错误发现率方法得到相同的结果：所有的配对比较都没有显著差异。接下来，让我们展示如何使用TukeyHSD()函数对平均值进行Tukey多重比较，并获得它们的置信区间。调用此函数的方式类似于Summary()函数。它将原始方差分析计算中的变量作为其参数之一(图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))
此图表示所有可能的成对测试和p值，以及95%的置信区间。可以根据您的选择更改默认的95%置信级别。由于所有置信线均为0，因此对于此示例，使用Tukey多重比较进行调整后没有显著不同的项。
8.4 Comparisons Among More than Two Groups Using Kruskal-Wallis Test（ 使用Kruskal-Wallis检验对两组以上人群进行比较）
① Kruskal-Wallis Test
以William Kruskal和W.Allen Wallis命名的Kruskal-Wallis检验或单因素秩和方差分析是检验样本是否来自同一分布的非参数方法。Kruskal-Wallis检验的参数等价性是单向方差分析。它将Mann-Whitney U检验扩展到两组以上。Kruskal-Wallis检验的零假设是各组的平均等级相同。与类似的单因素方差分析不同，非参数Kruskal-Wallis检验不假定基础数据为正态分布。它已被广泛应用于微生物组的研究。例如，测序后的微生物组数据不是正态分布的，并且包含一些很强的异常值。因此，使用秩而不是实际值是合适的，以避免测试受到异常值的存在或非正态分布的影响。测试统计需要四个主要步骤：
步骤1.将所有组中的所有数据一起按升序排列在单个序列中，即忽略组成员身份从1到N进行排序。
步骤2.通过平均排名位置来分配任何并列的值。
步骤3.总结不同的等级，例如R1R2R3…。对于每个不同的组。
步骤4.应用以下公式计算测试统计量：
Kruskal-Wallis检验统计量近似为卡方分布，如果n值“大”，k−为1个自由度。当每个数值大于或等于5时，该近似通常被认为是足够的。
②Compare Diversities Among Groups（比较不同群体之间的差异）
Kruskal-Wallis检验或Kruskal-Wallis单因素方差分析用于比较数据不服从正态分布的多组。该检验类似于两个样本的Wilcoxon秩和检验。我们首先使用我们的VDR−/−鼠标数据来说明此测试。
> library(dplyr)
> Data <- mutate(df_CH_G, Group = factor(df_CH_G$Group4, levels=unique (df_CH_G$Group4)))
获取描述性统计信息
> 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))
直方图显示，在这种情况下，组之间的值分布是不同的。现在，我们使用kruskal.test()函数进行Kruskal-Wallis检验来比较中值的差异。
> #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
检验统计值为5.2，p值大于0.05，也低于卡方表法：
> qchisq(0.950, 3)
[1] 7.815
因此，我们接受零假设H0：在5%显着水平上，4组的中位数在统计上是相等的。一般情况下，如果Kruskal-Wallis检验有意义，则进一步进行事后分析，以找出哪些组的水平彼此不同。在这种情况下，Kruskal-Wallis测试并不重要。为了说明这一点，我们进行了两个事后检验：Nemenyi检验和Dunn检验。与方差分析类似，我们可以选择调整p值的方法来控制家族错误率或控制错误发现率。当您在R或RStudio中输入？p.adjust时，它会显示一个指向文档“调整多个比较的P值”的链接。您可以通过此链接查看调整方法的详细信息。
多重比较的Nemenyi检验：Nemenyi测试通过DescTools包中的Nemenyi Test()函数执行。我们首先加载DescTools包并调用函数Nemenyi Test()。这个调整p值的方法应该是“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检验表明，采用Tukey调整法，在粪便和盲肠VDR基因敲除样本中，不同地点和不同基因型的CHO1多样性的平均秩差异不显著。然而，当两组的观察次数不相等时，Nemenyi检验是不合适的，而Dunn检验是合适的(Zar 2010)。我们运行Dunn检验，如下所示。
多重比较的Dunn检验：最流行的后Kruskal-Wallis检验是Dunn检验。我们可以使用FSA包中的dunnTest()函数执行Dunn测试。下面，我们调用函数dunnTest()并使用Benjamini和Hochberg方法调整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检验表明，盲肠和粪便中VdR−/−样品的CHO1多样性差异有统计学意义。但是，用Benjamini-Hochberg方法调整多重比较p值后，样本中的位置和基因型之间没有统计学意义。
③ Find Significant Taxa Among Groups（在群组中找到重要的分类群）
在本节中，我们使用Kruskal-Wallis检验来说明如何在组之间寻找重要的分类群。假设我们想知道在来自粪便和盲肠的Vdr、−/−和WT小鼠样本中是否存在任何重要的分类群。我们使用Kruskal-Wallis测试对数据集中的248个分类(细菌)中的每一个进行测试。首先，对丰富的数据进行归一化，并将数据制成数据框。一种标准化方法是使用对数变换。
> data<-log((abund_table+1)/(rowSums(abund_table)+dim(abund_table) [2]))
> df<-as.data.frame(data)
另一种归一化方法是将当量计数转换为相对丰度。
> df <- as.data.frame(abund_table/rowSums(abund_table))
然后，使用kruskal.test()函数和迭代R函数执行248个测试(每个测试针对一个细菌)。Kruskal.test()函数有几个关键组件：代码为“for(i in 1：Dim(Df)[2])”的所有分类群(柱)的测试都是循环进行的。·对于每个循环，使用代码“kw_test<-kruskal.test(df[，i]，g=Group4)”运行Kruskal-Wallis测试。·结果存储在数据框中，每个样本一行，KW测试的每个p值一列。·报告cat函数为“cat(Paste(”Kruskal-Wallis test for“，Names(Df)[i]，”“，i，”/“，dim(Df)[2]，”；p-value=“，kw_test$p.value，”\n“，sep=”“)”的测试次数
> 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=""))
+ }
检查数据框表以确保函数工作正常：
> #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
④多重测试与E值、FWER和FDR
文献中存在几种不同类型的多重测试校正。其中，Bonferroni修正较为保守。修正的方法很简单，就是用阿尔法除以测试次数。我们已经说明了使用Bonferroni，Holm，Tukey方法和两种版本的方法调整成两个版本的“假发现率”的p值调整，在第二节中使用ANOVA进行配对比较和检验。8.3.2，Kruskal-Wallis检验在Sect.Kruskal-Wallis检验中的帖子。8.4.2.。在本节中，我们介绍一般的多重测试修正：E值、系列误差率(FWER)和FDR。
E-value：E值是您进行多次测试时偶然出现的错误阳性的预期数量。您可以简单地将p值与对其执行测试的分类群数相乘就可以得到它：E-value = p-value X the number of tests（测试次数）。请注意，在E值上，基本校正是使用原始α，即来自测试的p值，而不是标称p值。
> #E-value
> KW_table$E.value <-KW_table$p.value * dim(KW_table)[1]
> KW_table$E.value
由于E值只是将p值乘以测试次数，因此它可以大于1。如果数据框中有许多分类群用于测试，则此校正方法不容易找到有意义的分类群。重要分类群是E值远小于1的分类群。以下代码用于检查是否已将E值添加到结果数据框中：
> #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是至少出现一次误报(I类错误)的概率。换句话说，这是您没有拒绝零假设H0的概率：在进行多个测试时，不同组之间没有差异。该公式由下式给出
式中，T=测试次数。在R中，为避免直接使用上述公式计算FWER带来的舍入误差，最好采用右尾二项分布检验来计算FWER。R码如下所示。
> #FWER
> KW_table$FWER <- pbinom(q=0, p=KW_table$p.value,size=dim(KW_table)[1], lower.tail=FALSE)
检查数据框以查看是否将FWER添加到结果数据框：
> #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)将错误发现率定义如下：FDR=expected proportion of erroneous rejections among all rejections（错误拒绝在所有拒绝中的预期比例）。在本例中，FDR是在进行多次测试时被接受为阳性的分类群中假阳性的比例。 Benjamini-Hochberg校正包括以下步骤：首先，从最小到最大对p值进行排序，并进行排序(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
接下来，使用以下公式和代码计算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
然后，通过使用以下代码指定目标FDR，并识别q值等于或小于指定阿尔法的排名列表的最后一项：
> #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值，因此程序返回负无限。最后，显示结果框表和选择的分类群：
> #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"
因为在这种情况下没有小于或等于指定的alpha=0.05的Q值，所以上面显示的5个分类群不是基于FDR的。它们是p值最小的选定分类群。本贾米尼-霍希伯格校正比上面提出的其他多次测试校正不那么严格，因此具有更高的灵敏度。FDR广泛应用于微生物群系和其他研究领域和许多R函数。
8.5 Summary
在这一章中，我们介绍了各种研究领域中常用的和经典的方法。其中一些在微生物组研究中得到了广泛的应用。我们用分步实施的方式说明了这些微生物组数据的分析方法。R系统中的状态。数据集来自我们自己的出版物(金等人)。2015年；王等人。2016)。读者可以使用本章提供的R代码和解释来分析自己的微生物组数据。在本章中，我们重点研究了单变量群落微生物组数据的假设检验。在接下来的一章里。将着重对多变量群落微生物组数据进行假设检验。
网友评论