【原】继14种差异方法比对后用什么方法进行微生物组差异分析（ZicoSeq）

写在前面

这个方法在2022年的时候发表了Microbiome，比较了目前比较多的差异分析方法。并且提出了自己的方法。包括在他之前还有一个NC比较了14种差异分析方法。这些内容都表明了差异分析在目前组学分析中的重要程度。可以给他们加上一个标签：差异分析与生物标志物挑选。下面是官方教程，请大家查看。建议微生物组差异分析使用组成数据，无需任何标准化方法。

• 官方教程：

  https://cran./web/packages/GUniFrac/vignettes/ZicoSeq.html

介绍

ZicoSeq是一种线性模型和基于排列的方法，用于处理零膨胀组成测序数据（如微生物组测序数据）的差异丰度分析。它包含以下组件：

• Winsorization以减少离群值的影响；

• 基于Beta混合先验的后验抽样，以处理采样变异性和零膨胀；

• 基于参考的多阶段归一化，以处理组成效应；

• 处理不同的特征-协变量关系；

• 基于排列的假发现率控制/全家系错误率控制用于多重检验校正，考虑特征数据中的相关结构。

目前，ZicoSeq支持以下内容：

计数数据或比例数据。对于计数和比例数据，使用基于参考的比率方法来考虑组成效应。当提供计数矩阵时，它提供一个选项来抽取潜在比例的后验样本，以考虑测序过程中的采样变异性。测试结果将在这些后验样本上进行汇总。

其他数据类型。作为一种通用方法，ZicoSeq可用于其他高维数据集的差异分析，例如转录组学、表观遗传学、代谢组学和蛋白质组学数据。

安装R包

# install.packages("GUniFrac")

导入R包，并运行示例数据

这个示例数据集包含喉部微生物组的OTU计数数据，可在GUniFrac软件包中获得。它包含60个受试者，其中32个是非吸烟者，28个是吸烟者（Charlson，Emily S.等人，2010）。在这个示例中，我们将在调整性别的情况下识别与吸烟有关的OTU，因为性别是一个潜在的混杂因素（我们在吸烟者中看到更多的男性）。

library(GUniFrac)

data(throat.otu.tab)
data(throat.meta)
comm <- t(throat.otu.tab)
meta.dat <- throat.meta
meta.dat

ZicoSeq.obj <- ZicoSeq(meta.dat = meta.dat, feature.dat = comm, 
                    grp.name = 'SmokingStatus', adj.name = 'Sex', feature.dat.type = "count",
                    # Filter to remove rare taxa
                    prev.filter = 0.2, mean.abund.filter = 0,  
                    max.abund.filter = 0.002, min.prop = 0, 
                    # Winsorization to replace outliers
                    is.winsor = TRUE, outlier.pct = 0.03, winsor.end = 'top',
                    # Posterior sampling 
                    is.post.sample = TRUE, post.sample.no = 25, 
                    # Use the square-root transformation
                    link.func = list(function (x) x^0.5), stats.combine.func = max,
                    # Permutation-based multiple testing correction
                    perm.no = 99,  strata = NULL, 
                    # Reference-based multiple stage normalization
                    ref.pct = 0.5, stage.no = 6, excl.pct = 0.2,
                    # Family-wise error rate control
                    is.fwer = TRUE, verbose = TRUE, return.feature.dat = T)

结果解释

• p.raw: 基于置换的原始p值，如果perm.no很小，则其精度不高。

• p.adj.fdr: 基于置换的假阳性发现率调整后的p值。

• p.adj.fwer: 基于置换的簇内错误率调整后的p值（West-Young），仅当ZicoSeq中的is.fwer = T时才有此选项。

  建议将perm.no至少设置为999以获得准确的簇内错误率调整后的p值。

使用ZicoSeq注意

微生物组数据可能包含具有极大数的离群值。这些离群值可能降低统计功效，或更严重地增加第一类错误率。在ZicoSeq中，我们使用Winsorization将那些极大计数替换为某个分位数（例如，97%）。我们建议将is.winsor = TRUE设置为真实数据分析。所使用的特定分位数取决于预期异常值的比例。

后验抽样

通过设置is.post.sample = TRUE可以启用后验抽样。它只能用于数据类型为“计数”的情况。对于微生物组数据，大多数分类单元都很罕见，在有限的测序深度下，导致数据中存在大量的零。由于零的概率强烈依赖于测序深度，来自不同样本的零是不可比的。在ZicoSeq中，我们提供了一种通过跨样本汇集信息使用经验贝叶斯方法推断/插补潜在真实比例的方法。估计的比例取决于测序深度和样本间比例的估计（先验）分布。为了更好地建模先验分布，我们使用了Beta混合分布而非单组分Beta分布。当测序深度与感兴趣的协变量（例如，病例和对照在测序深度上有所不同）相关联时，后验抽样可以有效地降低第一类错误率，而且比稀释法更有力。当测序深度不是混杂因素时，is.post.sample会略微提高功效，并且可以禁用以分析大型数据集。

Omnibus测试

分类单元丰度与感兴趣的协变量之间的关系可能因分类单元而异，单一的转换函数可能不足以捕获各种关系。ZicoSeq允许指定多个可能的转换函数（link.func）并执行全路线测试。虽然对于微生物组来说，由于其可解释性，对数变换通常被广泛使用，但我们发现幂次变换可能更加有效。默认情况下，采用平方根变换。

R^2 (方差解释百分比)

将提供一个 R^2 值矩阵（特征数×转换函数数），R^2 是一种效应量度，可用于在调整其他协变量的情况下评估微生物丰度与感兴趣的协变量之间的关联强度。当多个转换函数被使用时，即多重转换测试时，将有多个 R^2 值，每个 R^2 对应于一个特定的转换函数。

ZicoSeq 输出可视化

ZicoSeq.plot 函数生成一张火山图，其中 y 轴为 log10(经过调整的) p 值，x 轴为带符号的 R^2，符号表示关联方向，基于回归系数的符号确定（对于多类别变量，不适用符号）。当使用多个转换函数时，将使用最大的 R^2 值。通过特定的显著性截断点的微生物分类名称将打印在图形上。当数据类型为计数和比例时，将显示平均丰度和流行度；当数据类型为其他时，将显示特征的平均值和标准差。在使用 plot 函数时，用户需要设置 return.feature.dat = T。

ZicoSeq.plot(ZicoSeq.obj, meta.dat, pvalue.type = 'p.adj.fdr', cutoff = 0.1, text.size = 10,             out.dir = NULL, width = 10, height = 6)

根际互作生物学研究室 简介