使用 ALDEx2 进行差异分析
ALDEx2 是进行微生物组差异分析较为常见的方法。该方法包含两个基本过程：
1.用原始输入数据生成每个分类单元的后验概率分布；然后将该分布进行中心对数变换。2.将变换后的值，用参数或非参数检验进行单变量统计检验，并返回 p 值和 Benjamini-Hochberg 校正后的 p 值。
安装 ALDEx2
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("ALDEx2")以 ALDEx2 自带的数据集为例，导入数据：
> library(ALDEx2)
>
> data(selex)
> #subset only the last 400 features for efficiency
> selex.sub <- selex[1:400,]
> selex.sub[1:4,1:4]
X1_ANS X1_BNS X1_CNS X1_DNS
A:D:A:D 347 271 396 317
A:D:A:E 436 361 461 241
A:E:A:D 476 288 378 215
A:E:A:E 513 307 424 381
> conds <- c(rep("NS", 7), rep("S", 7))
> conds
[1] "NS" "NS" "NS" "NS" "NS" "NS" "NS" "S" "S" "S" "S"
[12] "S" "S" "S" 我们既可以以模块化的方式分步进行 ALDEx2 分析，也可以直接一行命令完成，下面将分别介绍这两种方法。
分步运行 ALDEX
简单来说，这种方法的过程只是依次调用 aldex.clr()，aldex.ttest() 和 aldex.effect() 函数，然后将数据合并到一个对象中。
第一步：用 aldex.clr() 函数进行中心对数比转换
aldex.clr() 函数需要四个变量：
•OTU count 表•分组信息•蒙特卡罗实例数•是否输出函数运行过程（TRUE or FALSE）
作者建议若使用 t-test，则需 mc.samples >= 128 ；若需计算严格的效应值大小，建议 mc.samples >= 1000；若用 ANOVA，建议 mc.samples >= 16。
> #this operation is fast.
> x <- aldex.clr(selex.sub, conds, mc.samples=128, verbose=TRUE)第二步：用 aldex.ttest() 函数进行 Welch's t 检验和秩和检验
aldex.ttest() 函数需要三个输入变量：
•aldex.clr() 输出的 aldex 对象•分组信息•是否进行配对检验（TRUE or FALSE）
aldex.ttest() 函数将输出：
•we.ep (Welch's t 检验的 p 值)•we.eBH ( Welch's t 检验校正后的 p 值)•wi.ep (秩和检验的 p 值)•wi.eBH (秩和检验校正后的 p 值)。
> x.tt <- aldex.ttest(x, conds, paired.test=FALSE)
> head(x.tt)
we.ep we.eBH wi.ep wi.eBH
A:D:A:D 0.5506426178 0.748600650 0.3821386946 0.583273523
A:D:A:E 0.1980316100 0.454646926 0.0860513185 0.255437818
A:E:A:D 0.0001074824 0.002385507 0.0005827506 0.006854623
A:E:A:E 0.0001838241 0.003315360 0.0005827506 0.006854623
A:D:C:D 0.3176806833 0.565706881 0.4335846445 0.611279169
A:D:C:E 0.2929417329 0.561897763 0.4155876675 0.596003077除了 Welch's t 检验和秩和检验，我们还可以用 aldex.glm() 函数执行 glm 和Kruskal-Wallis 单因素方差分析。但运行起来比较慢。aldex.glm() 函数也会输出 kw.ep (Kruskal-Wallis 检验的 p 值)、kw.eBH (Kruskal-Wallis 检验校正后的 p 值)、glm.ep (glm 检验的 p 值)和 glm.eBH (glm 检验校正后的 p 值)。
> # based on the documentation this is slow and not evaluated
> x.kw <- aldex.kw(x)第三步：用 aldex.effect() 函数估计效应量大小。
aldex.effect() 函数会分别估计两种情况下的效应量大小以及条件值之间的范围。它需要四个输入变量：aldex.clr() 输出的 aldex 对象、分组信息、是否包括所有样本的 clr 中值，以及是否输出函数运行过程（TRUE or FALSE）。
> x.effect <- aldex.effect(x, verbose=FALSE)aldex.effect() 函数的输出包括:
•rab.all (特征中所有样本的 clr 中值)•rab.win.VdrFecal (VdrFecal 组样本的 clr 中值)•rab.win. VdrCecal (VdrCecal 组样本的 clr 中值)•dif.btw (VdrFecal 和 VdrCecal 组间的 clr 值 median difference)•dif.win (VdrFecal 和 VdrCecal 组内的 clr 值最大差异中值)•effect (效应量中值 diff.btw /max(diff.win)•overlap (效应量包含 0 的比例)
第四步：将所有数据合并到数据框中
> x.all <- data.frame(x.tt,x.effect)
> head(x.all)
we.ep we.eBH wi.ep wi.eBH rab.all rab.win.NS rab.win.S
A:D:A:D 0.5506426178 0.748600650 0.3821386946 0.583273523 1.6085267 1.5389689 2.3266274
A:D:A:E 0.1980316100 0.454646926 0.0860513185 0.255437818 1.9028527 1.7294433 4.0641478
A:E:A:D 0.0001074824 0.002385507 0.0005827506 0.006854623 3.9764334 1.6330604 10.7585071
A:E:A:E 0.0001838241 0.003315360 0.0005827506 0.006854623 3.9037136 1.8734491 13.7536087
A:D:C:D 0.3176806833 0.565706881 0.4335846445 0.611279169 0.5295745 0.6367628 -0.2396709
A:D:C:E 0.2929417329 0.561897763 0.4155876675 0.596003077 0.5029567 0.5999841 -0.1859935
diff.btw diff.win effect overlap
A:D:A:D 0.7986369 1.827900 0.3226198 0.3303572225
A:D:A:E 2.3164250 2.378344 0.7951613 0.1870826180
A:E:A:D 9.2302125 2.726573 3.3794301 0.0001566327
A:E:A:E 11.9669585 2.748247 4.0855989 0.0001566327
A:D:C:D -0.8688044 2.617255 -0.2949049 0.3608018441
A:D:C:E -0.6931246 2.419352 -0.2556930 0.3719376959
查看两组之间的差异：
> sig_by_both <- which(x.all$we.ep < 0.05 & x.all$wi.ep < 0.05)
> sig_by_both
[1] 3 4 7 8 11 12 15 16 21 22 23 24 40
[14] 59 60 63 64 71 72 75 76 101 103 164 181 182
[27] 183 184 244 262 263 264
> sig_by_both_fdr <- which(x.all$we.eBH < 0.05 & x.all$wi.eBH < 0.05)
> sig_by_both_fdr
[1] 3 4 7 8 11 12 15 16 21 22 23 24 59
[14] 60 63 64 71 72 75 103 164 183 184 244 263 264需要注意的是结果中除 p 值外，所有值都是经过 log2 转换后的值。此外，因为 aldex.clr() 函数使用蒙特卡罗方法对数据进行采样，所以所有结果中的值都是由 aldex.clr() 函数中的 mc.samples 变量给出的狄利克雷实例数的平均值。
一行命令进行 ALDEX 差异分析
目前，aldex 函数仅限于双样本检验和单因素方差分析。
> x.all <- aldex(selex.sub, conds, mc.samples=16, test="t", effect=TRUE,include.sample.summary=FALSE, denom="all", verbose=FALSE)
aldex.clr: generating Monte-Carlo instances and clr values
operating in serial mode
computing center with all features
aldex.ttest: doing t-test
aldex.effect: calculating effect sizes
Warning message:
In aldex.clr.function(reads, conds, mc.samples, denom, verbose, :
values are unreliable when estimated with so few MC smps
> head(x.all)
rab.all rab.win.NS rab.win.S diff.btw
A:D:A:D 1.6256967 1.5578277 2.2796353 0.6327612
A:D:A:E 1.9328492 1.7232726 3.9046551 2.2631860
A:E:A:D 4.2077070 1.6279158 10.7661169 9.2579901
A:E:A:E 4.0720951 1.8623388 13.7414099 11.8022926
A:D:C:D 0.5221376 0.5564199 0.2917089 -0.3592074
A:D:C:E 0.4919644 0.6149000 -0.1576435 -0.6227079
diff.win effect overlap we.ep
A:D:A:D 1.957942 0.3264417 0.350881391 0.5793775170
A:D:A:E 2.470586 0.8957338 0.160730987 0.1725091804
A:E:A:D 2.932941 2.9414700 0.001251685 0.0001001308
A:E:A:E 3.024319 3.3956943 0.001251685 0.0001684960
A:D:C:D 2.100723 -0.1510346 0.438598090 0.4208066570
A:D:C:E 2.541803 -0.2154117 0.339290476 0.3833382223
we.eBH wi.ep wi.eBH
A:D:A:D 0.758248318 0.4416885198 0.616740334
A:D:A:E 0.438015342 0.0725160256 0.215931319
A:E:A:D 0.002223689 0.0005827506 0.006975847
A:E:A:E 0.003049815 0.0005827506 0.006975847
A:D:C:D 0.685040496 0.6241258741 0.751821123
A:D:C:E 0.623946746 0.4764714452 0.642362543这里指定 test="t"，effect=TRUE。"t" 选项表示使用 Welch's t 检验和秩和检验。若test="glm"，则 effect 应指定为 FALSE。输出结果包括原始 p 值和校正后的值。
Difference Plot，Effect Size 和 Effect Plot
在下游分析中，我们将合并后的数据框进行可视化。
Bland-Altman Plot
Bland-Altman Plot，也叫做 difference plot 或 Tukey mean-difference。可分析两种不同方法之间的一致性。在 ALDEx2 可用 aldex.plot() 函数可视化比较 Median Log2 relative abundance 和 Median Log2 btw-Condition diff：
> aldex.plot(x.all, type="MA", test="welch", xlab="Log-ratio abundance",ylab="Difference")参数中，type = MA 指定要图的类型；test="welch" 表示使用 Welch's t 检验来计算显著性；cutoff=0.15 指定使用 0.15 作为 Benjamini-Hochberg FDR 的 cutoff，默认值为 0.1；all.cex=0.7 指定点的大小；called.cex=1.1 用来指定显著的点的大小，rare.col="grey" 为稀有分类单元指定颜色，默认为黑色；called.col="red" 指定用红点来表示显著的分类群。
Effect Size and Effect Size Plot
在 ALDEx2 中，效应量大小被定义为组间差异(diff.btw)和组内最大差异(diff.win或方差)的平均比率。
我们可用 aldex.plot() 函数绘制组间差异中值与组内差异中值，以可视化样本数据的差异丰度。
> par(mfrow=c(1,2))
> aldex.plot(vdr_all, type="MW", test="welch",cutoff=0.15, all.cex=0.7,
called.cex=1.1, rare.col="black", called.col="red")
> aldex.plot(vdr_all, type="MW",test="wilcox", cutoff=0.15, all.
cex=0.7, called.cex=1.1, rare.col="black", called.col="red")其中，type="MW" 指定图的类型为 MW：a difference between to a variance within。用 Welch's t 检验或秩和检验计算显著性。
左右两图分别使用 Welch's t 检验和秩和检验。红点代表差异的分类群。若图中没有红点，则表明未检测到显著性。灰点代表丰富的分类群，但不显著。灰色线条代表组内和组间值的等价线。黑点是比平均分类单元丰度低且不显著的分类单元。一般来说，这些分类群很难精确估计。
一般来说，p 值不如效应量大小显著。如果样本量足够大，则 0.5 或更大的效果量可被认为有一定生物学相关性。在 ALDEx2 中，建议效应大小 cutoff 设置为 1.5-2 ，overlap cutoff 设置为 0.01。
引用链接
[1] http://www.bioconductor.org/packages/release/bioc/vignettes/ALDEx2/inst/doc/ALDEx2_vignette.pdf
[2] Statistical Analysis of Microbiome Data with R