4. phyloseq包详解：从零搭起你的物种丰度分析骨架

说到 R 语言做宏基因组可视化，phyloseq 这个包你绕不开。它就像分析流程里的“集装箱”，把 OTU 丰度表、物种注释表、样本元数据全整合在一个对象里。我刚开始用 phyloseq 时，总觉得这步有点“黑盒”，后来踩过不少坑，才摸透了它到底是怎么把这些数据捏在一起的。

4.1 phyloseq对象创建：把三张表合并成一个“数据容器”

要建一个 phyloseq 对象，你手头至少得有这三样东西：

• OTU 丰度表（行是 OTU，列是样本）
• 物种注释表（行是 OTU，列是分类级别）
• 样本元数据表（行是样本，列是实验信息）

这三张表说白了，就是分开放着谁也用不上，必须塞进同一个结构里才方便分析。

4.1.1 数据准备：先把表格整成 phyloseq 能认的格式

phyloseq 不会凭空变出数据。你得先把原始的 OTU 表格读进来，再转成它要求的格式。我习惯用 readr::read_csv() 或 data.table::fread() 读原始数据，然后用 otu_table() 函数来构建正确格式：

library(phyloseq)
library(tidyverse)

# 读取原始 OTU 丰度表（行是 OTU，列是样本）
otu_raw <- read.csv("otutable_raw.csv", row.names = 1)

# 构建 OTU 丰度矩阵
otu_mat <- as.matrix(otu_raw)
OTU <- otu_table(otu_mat, taxa_are_rows = TRUE)

物种注释表也一样，把它转成 taxonomyTable 对象：

# 读取物种注释表（行是 OTU，列是 Kingdom/Phylum/Class 等）
tax_raw <- read.csv("taxonomy_raw.csv", row.names = 1)

# 构建物种注释矩阵
tax_mat <- as.matrix(tax_raw)
TAX <- tax_table(tax_mat)

至于样本元数据表，它必须是个 sample_data 对象，而且行名要与 OTU 表的列名完全一致——这就是个坑。我曾经因为排序不同，结果样本匹配不上，查了半天才发现是行名顺序的锅。

# 读取样本元数据（行是样本，列是分组/处理/时间等）
sam_raw <- read.csv("sample_metadata.csv", row.names = 1)

# 构建 sample_data 对象
SAM <- sample_data(sam_raw)

4.1.2 合并三张表：phyloseq 函数怎么做到“一键整合”？

三张表各自就绪后，用 phyloseq() 函数把它们包在一起，就跟打包行李一样：

# 合并 OTU + 物种注释 + 样本元数据
physeq <- phyloseq(OTU, TAX, SAM)

# 检查对象基本信息
physeq
# 输出示例：
# phyloseq-class experiment-level object
# otu_table()   OTU Table:         [ 200 taxa and 24 samples ]
# tax_table()   Taxonomy Table:    [ 200 taxa by 6 taxonomic ranks ]
# sample_data() Sample Data:       [ 24 samples by 8 sample variables ]

合并完成后，你就可以用 physeq 这个对象直接跑各种分析函数，比如 plot_bar()、estimate_richness() 等。

4.2 常用访问器函数：从 phyloseq 对象里“掏出”你想要的表

对象创建好之后，你不能总盯着 physeq 看。实际分析中，你得频繁从里面提取 OTU 表、物种注释、样本数据，或者修改它们。phyloseq 提供了一组配套的访问器函数，说白了就是“打开集装箱取零件”。

4.2.1 提取函数：怎么把数据从 phyloseq 里拿出来？

名字很好记，直接用和创建时一样的函数名就能提取：

# 提取 OTU 丰度表
otutab <- otu_table(physeq)
class(otutab)  # "otu_table" 对象
otutab_df <- as.data.frame(otutab)  # 转回普通 data.frame

# 提取物种注释表
taxtab <- tax_table(physeq)
taxtab_df <- as.data.frame(taxtab)

# 提取样本元数据
samtab <- sample_data(physeq)
samtab_df <- as.data.frame(samtab)

还有一个很少人用但我觉得很实用的函数 nsamples() 和 ntaxa()，能快速返回样本数和 OTU 数量：

nsamples(physeq)  # 样本数量
ntaxa(physeq)     # OTU 数量

4.2.2 修改函数：想在 phyloseq 对象里替换某张表？

你可能会想，我要是把 OTU 表经过去除低丰度 OTU 之后，怎么放回原对象？ 很多时候，你需要在不重建整个对象的情况下更新内部的数据。这时用相同的函数名，把 new 参数塞进去就行。

# 假设我们过滤了低丰度 OTU，得到了新的 OTU 表 new_otutab
# 把它替换回原 phyloseq 对象
physeq_new <- otu_table(physeq) <- new_otutab

# 或者用更明确的写法
physeq_new <- phyloseq(new_otutab, tax_table(physeq), sample_data(physeq))

我比较推荐第二种写法——显式地把三张表重新组合，不容易手滑改错。如果你只是临时修改某个表，可以用 sample_data(physeq)$ 语法直接操作：

# 给样本元数据增加一列：把 day 转为 factor
sample_data(physeq)$day_factor <- factor(sample_data(physeq)$day)

# 修改已有的列
sample_data(physeq)$group <- recode(sample_data(physeq)$group, 
                                     "Control" = "Ctrl", 
                                     "Treatment" = "Treat")

4.3 一个快速检查表：你该选哪种提取/修改方式？

场景	推荐操作	注意事项
只读 OTU 丰度表	otu_table(physeq)	结果是个 phyloseq 内部对象，不能直接运算
读完后做矩阵运算	as(otu_table(physeq), "matrix")	用 as() 强制转换，搞定再转回来
添加一列分组信息	sample_data(physeq)$newcol <- value	改完后要重新赋值：physeq <- physeq
替换整个物种注释表	重建：phyloseq(OTU, newTAX, SAM)	这样最安全，避免索引错乱
快速查看数据概况	直接用 physeq 打印	会显示各表的维度

4.4 嵌套操作与管道：更优雅的提取链

看过我代码的人都知道，我特别喜欢用 magrittr 的管道操作来连写提取和转换：

library(dplyr)

# 从 physeq 中提取 OTU 表，转为数据框，计算每行（每个OTU）的总丰度
otu_total <- otu_table(physeq) %>%
  as("matrix") %>%
  as.data.frame() %>%
  mutate(rowsum = rowSums(across(where(is.numeric))))

说白了，phyloseq 的强大之处不在于它有多少新奇的功能，而在于它把这三张表严格地约束在一起。你只要掌握了对象创建和访问器函数，后面那些可视化函数（比如 plot_bar、plot_heatmap）才能稳稳地建立在统一的数据结构上。