R语言phyloseq计算优势分类群平均相对丰度方法是否正确

问题背景 #

我想确认自己计算任意分类群平均相对丰度的方法是否正确：基于phyloseq对象的GlobalPatterns数据集，计算每个科（或任意分类阶元）的相对丰度（百分比），进而识别高丰度优势分类群。
编写的实现代码如下：

data("GlobalPatterns")

T <- GlobalPatterns %>% 
    tax_glom(., "Family") %>% 
    transform_sample_counts(function(x)100* x / sum(x)) %>% psmelt() %>% 
    arrange(OTU) %>% rename(OTUsID = OTU) %>% 
    select(OTUsID, Family, Sample, Abundance) %>%
    spread(Sample, Abundance)

T$Mean <- rowMeans(T[, c(3:ncol(T))])

FAM <- T[, c("Family", "Mean" ) ]

#order data frame  
FAM <- FAM[order(dplyr::desc(FAM$Mean)),]
rownames(FAM) <- NULL

head(FAM)
          Family          Mean
1    Bacteroidaceae     7.490944
2   Ruminococcaceae     6.038956
3   Lachnospiraceae     5.758200
4 Flavobacteriaceae     5.016402
5  Desulfobulbaceae     3.341026
6            ACK-M1     3.242808

按上述代码计算结果，GlobalPatterns数据集共包含26个样本、19216个OTU，其中Bacteroidaceae是丰度最高的科，在26个样本中的平均占比为7.49%。

核心疑问：使用T$Mean <- rowMeans(T[, c(3:ncol(T))])计算任意给定分类群跨样本平均丰度的操作是否正确？

结论 #

• 你的计算逻辑整体是正确的，符合算术平均相对丰度的计算要求：
  • 前置处理步骤符合phyloseq分析的标准流程：tax_glom(., "Family")完成科水平的分类单元聚合，transform_sample_counts(function(x)100* x / sum(x))将每个样本的丰度标准化为百分比形式的相对丰度（单样本内所有分类单元丰度和为100%）。
  • 经spread(Sample, Abundance)转换后，每一行对应一个聚合后的科，第3列到最后一列依次存储该科在每个样本中的相对丰度值，此时用rowMeans计算行均值，得到的就是该科在所有样本中的算术平均相对丰度，和你的分析目标匹配。
• 有两个可优化的细节：
  • 建议在rowMeans中加入na.rm = TRUE参数，即写为T$Mean <- rowMeans(T[, c(3:ncol(T))], na.rm = TRUE)，避免中途数据处理产生的NA值导致均值计算异常。
  • 该计算方式是给所有样本赋予相等权重的算术平均，如果你的分析需要考虑样本测序深度差异、或者分组样本量不均衡需要做加权计算，需要根据研究设计调整均值计算的权重逻辑，不能直接用等权平均。
  • 原始数据集19216个OTU经科水平聚合后行数会大幅减少，属于正常现象，是同科OTU合并导致的。

内容的提问来源于stack exchange，提问作者abraham