使用R在宏基因组数据树上聚合多变量的技术问询

Hey there! 我看你已经基于phyloseq的GlobalPatterns数据集搭好了分类树结构，现在需要在每个分类层级（界、门、纲…）上做数值聚合对吧？这在宏基因组分析里太常见了，用data.tree的工具就能轻松搞定，我给你几个实用的方案：

首先我们需要把OTU对应的样本丰度数据挂载到树的叶子节点上，这样后续才能往上聚合：

library(phyloseq)
library(dplyr)
library(data.tree)

# 你的初始代码
data(GlobalPatterns)
mydata <- merge(GlobalPatterns@tax_table, GlobalPatterns@otu_table, by = "row.names")
mydata <- mydata %>% select(-one_of("Row.names"))
mydata$pathString <- apply(mydata[,1:7], 1, paste, collapse="/")
tree <- as.Node(mydata)

# 提取样本列名称（也就是丰度列）
sample_cols <- colnames(mydata)[8:ncol(mydata)]

# 将每个OTU的丰度赋值给对应的叶子节点
for (i in seq(nrow(mydata))) {
  target_node <- tree$Get(mydata$pathString[i])
  for (col in sample_cols) {
    target_node[[col]] <- mydata[i, col]
  }
}

Aggregate是data.tree专门用来层级聚合的函数，配合**后序遍历（post-order）**可以实现从OTU（叶子节点）往上逐层求和：

# 对每个样本列，在所有分类层级上聚合求和
for (col in sample_cols) {
  tree$Do(function(node) {
    # 叶子节点（OTU）直接保留原始值，非叶子节点聚合子节点的数值
    node[[paste0(col, "_sum")]] <- if (node$isLeaf) {
      node[[col]]
    } else {
      Aggregate(node, attribute = col, aggFun = sum)
    }
  }, traversal = "post-order")
}

为什么用post-order？因为要先计算所有子节点的聚合值，再算父节点的，这样结果才准确

如果你的需求不是简单求和，而是计算均值、或者过滤掉低丰度OTU再聚合，可以自定义聚合函数：

# 示例：计算每个分类层级的平均丰度
tree$Do(function(node) {
  if (!node$isLeaf) {
    node$mean_abundance <- Aggregate(node, attribute = sample_cols[1], aggFun = mean)
  }
}, traversal = "post-order")

# 示例：只聚合丰度大于10的子节点
custom_agg <- function(x) sum(x[x > 10])
tree$Do(function(node) {
  if (!node$isLeaf) {
    node[[paste0(sample_cols[1], "_filtered_sum")]] <- Aggregate(node, attribute = sample_cols[1], aggFun = custom_agg)
  }
}, traversal = "post-order")

如果样本数量多，可以用purrr批量处理，避免写循环：

library(purrr)

map(sample_cols, function(col) {
  tree$Do(function(node) {
    node[[paste0(col, "_sum")]] <- if (node$isLeaf) {
      node[[col]]
    } else {
      Aggregate(node, col, sum)
    }
  }, traversal = "post-order")
})

聚合完成后，你可以轻松提取某个层级的结果，比如提取门（Phylum）级别的聚合丰度：

# 找到所有门级别的节点（假设Kingdom是第1层，Phylum是第2层）
phylum_nodes <- tree$Traverse(function(node) node$level == 2)

# 转成数据框方便后续分析
phylum_abundance_df <- map_dfr(phylum_nodes, function(node) {
  tibble(
    phylum_name = node$name,
    across(all_of(paste0(sample_cols, "_sum")), ~.x)
  )
})

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