在R中将微生物组分类数据转换为Newick格式构建进化树

问题背景 #

我刚学生物信息学分析，现有约33000条序列的微生物组数据，行代表序列，列依次为domain、phylum、class、order、family、genus、species分类层级。需要将这些数据转换为Newick格式以构建进化树。此前尝试的方法要么无法适配大规模数据，要么依赖的包运行异常，还有的方案用Python实现难以理解。

样本数据 #

# 样本分类数据
tax_data <- structure(c("A", "C", "B", "A", "C", "B", "C", "B", "B", "C", "C", "B", "C", "C", "A", "C", "A", "B", "B", "B", "E", "G", "F", "H", "F", "I", "D", "H", "H", "E", "I", "D", "D", "H", "G", "H", "G", "I", "H", "E", "L", "K", "Q", "N", "O", "Q", "J", "K", "O", "Q", "M", "K", "M", "Q", "P", "N", "J", "Q", "M", "O", "R", "Y", "Z", "W", "V", "Y", "S", "V", "R", "T", "X", "R", "R", "W", "R", "X", "Y", "Z", "Z", "Z", "b", "e", "c", "c", "b", "j", "h", "d", "h", "a", "g", "d", "f", "a", "f", "h", "d", "g", "i", "j", "q", "p", "t", "n", "u", "p", "q", "u", "t", "s", "q", "k", "k", "k", "l", "s", "o", "v", "m", "m", "n H", "h A", "x J", "h T", "a H", "t X", "h O", "m J", "g B", "f Z", "w X", "l S", "s F", "w R", "v U", "z X", "e B", "c O", "n R", "b J"), dim = c(20L, 7L), dimnames = list(c("seq_1", "seq_2", "seq_3", "seq_4", "seq_5", "seq_6", "seq_7", "seq_8", "seq_9", "seq_10", "seq_11", "seq_12", "seq_13", "seq_14", "seq_15", "seq_16", "seq_17", "seq_18", "seq_19", "seq_20"), c("domain", "phylum", "class", "order", "family", "genus", "species")))

解决方案代码 #

以下代码基于dplyr和ape包实现，可高效处理大规模（33000条）数据，自动从分类层级生成Newick格式：

# 安装并加载所需包
if (!require("dplyr")) install.packages("dplyr")
if (!require("ape")) install.packages("ape")
library(dplyr)
library(ape)

# 将矩阵转换为数据框，并保留序列名
tax_df <- as.data.frame(tax_data) %>%
  mutate(seq_id = rownames(.))

# 1. 去重：同一物种的序列只保留一条分类路径（减少冗余）
unique_tax <- tax_df %>%
  distinct(domain, phylum, class, order, family, genus, species, .keep_all = TRUE)

# 2. 构建分类路径字符串（用分隔符连接各层级）
unique_tax <- unique_tax %>%
  mutate(tax_path = paste(domain, phylum, class, order, family, genus, species, sep = ";"))

# 3. 自定义函数：从分类路径生成Newick字符串
tax_path_to_newick <- function(paths) {
  # 拆分路径为层级列表
  path_list <- strsplit(paths, ";")
  
  # 构建嵌套的分类群结构
  tree_list <- list()
  for (path in path_list) {
    current_node <- tree_list
    for (level in path) {
      if (!level %in% names(current_node)) {
        current_node[[level]] <- list()
      }
      current_node <- current_node[[level]]
    }
  }
  
  # 递归生成Newick字符串
  generate_newick <- function(node) {
    if (length(node) == 0) {
      return("")
    }
    children <- lapply(names(node), function(child) {
      sub_newick <- generate_newick(node[[child]])
      if (sub_newick == "") {
        return(child)
      } else {
        return(paste0("(", sub_newick, ")", child))
      }
    })
    return(paste(children, collapse = ","))
  }
  
  newick_str <- paste0("(", generate_newick(tree_list), ");")
  return(newick_str)
}

# 4. 生成Newick字符串
newick_result <- tax_path_to_newick(unique_tax$tax_path)

# 5. 保存Newick到文件（可选）
writeLines(newick_result, "microbiome_tree.newick")

# 验证：读取并查看树结构
tree <- read.tree(text = newick_result)
plot(tree, cex = 0.6)

代码说明 #

• 去重步骤：避免同一物种的重复序列生成冗余分支，大幅提升处理33000条数据的效率。
• 分类路径构建：将层级转换为统一格式，方便后续嵌套结构生成。
• 递归生成Newick：自动识别分类层级的嵌套关系，生成符合标准的Newick格式，无需手动干预。
• 可扩展性：若后续增加分类层级（如subspecies），只需修改tax_path的拼接逻辑即可。

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