📚 课程目录第4章数据预处理

第4章：数据导入——用phyloseq读取OTU/ASV表、分类注释和样本元数据

好，这一章我们来解决一个最实际的问题：怎么把手上那几张表格塞进phyloseq里。

我刚开始用R做宏基因组的时候，数据就是三个独立的CSV文件。一个物种丰度表，一个注释表，一个样本信息表。当时我还傻乎乎地手动拼数据框……嗯，后来发现phyloseq专门就是为了干这个活而设计的。你只需要把三个表准备好，几行代码就搞定了。

4.1 phyloseq对象长什么样？

在动手之前，我们先搞清楚目标：我们要造一个phyloseq对象。它本质上是一个容器，把三类数据打包在一起：

otu_table：物种丰度矩阵（行是OTU/ASV，列是样本）
tax_table：分类注释矩阵（行是OTU/ASV，列是界门纲目科属种）
sample_data：样本元数据（行是样本，列是实验条件）

说白了，这三张表必须能通过“行名”或“列名”互相连上。我见过有人直接把不同顺序的表塞进去，结果全部乱掉——这就是最常见的坑。

核心原则：OTU表列名 == 样本表行名；OTU表行名 == 注释表行名。名称必须一一对应，顺序无关紧要，phyloseq会自动对齐。

4.2 表格的格式要求

我们先看每张表需要长成什么样子。我直接拿我项目里整理好的结构展示给你：

表格类型	行（行名）	列（列名）	典型格式
OTU/ASV丰度表	OTU_ID（如 OTU_1）	样本名（如 S1, S2）	CSV 或 TSV，数值型
分类注释表	OTU_ID（与丰度表一致）	Kingdom/Phylum/Class/Order/Family/Genus/Species	CSV 或 TSV，字符型
样本元数据	样本名（与丰度表列名一致）	分组、时间、部位等	CSV 或 TSV，可含因子

你看看这三张表，是不是很像Excel里三个sheet？对，就是这么简单。

4.3 读取并构建phyloseq对象

假设你手头有三个文件：otu_table.csv、tax_table.csv、sample_metadata.csv。我给大家演示一下标准流程。

4.3.1 加载包

先加载phyloseq，它会顺带把依赖的ggplot2、vegan等一起拉进来。

library(phyloseq)
library(readr)      # 用readr读取CSV更稳定
library(tibble)     # 处理行名更方便

4.3.2 读取三张表

# 1. 读取OTU/ASV丰度表
otu_mat <- read.csv("otu_table.csv", row.names = 1)
# row.names = 1 表示第一列是OTU ID，会被自动设为行名

# 2. 读取分类注释表
tax_mat <- read.csv("tax_table.csv", row.names = 1)
# 注意：tax_mat 必须是字符矩阵，不能有因子

# 3. 读取样本元数据
samp_df <- read.csv("sample_metadata.csv", row.names = 1)
# row.names = 1 表示第一列是样本名

这里有个小细节。我遇到过好几次分类表里混进了因子，导致tax_table()报错。所以我建议你显式转成字符矩阵：

tax_mat <- as.matrix(tax_mat)
tax_mat <- apply(tax_mat, 2, as.character)  # 确保全是字符

4.3.3 转为phyloseq格式

OTU <- otu_table(otu_mat, taxa_are_rows = TRUE)
TAX <- tax_table(tax_mat)
SAM <- sample_data(samp_df)

taxa_are_rows = TRUE 表示OTU ID在行上。如果你的丰度表是“样本在行、OTU在列”——比如直接从某些流程导出的格式——那就设成 FALSE。我个人习惯永远是OTU在行，这样更符合生物信息学的直觉。

4.3.4 合并成phyloseq对象

ps <- phyloseq(OTU, TAX, SAM)
ps

如果一切正常，控制台会输出一个摘要信息：

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 856 taxa and 24 samples ]
tax_table()   Taxonomy Table:    [ 856 taxa by 7 taxonomic ranks ]
sample_data() Sample Data:       [ 24 samples by 6 sample variables ]

看到这个输出，恭喜你，数据导入成功了。

4.4 用简单方法来检验数据完整性

不要急着开始分析。先做几个快速检查：

sum(taxa_sums(ps) == 0)——看看有没有在所有样本中都是0的OTU
sum(sample_sums(ps) == 0)——看看有没有零reads的样本
ntaxa(ps) 和 nsamples(ps)——对比原始数据，确保数量一致

我曾经在导入后直接跑plot_bar，结果发现有个样本一条序列都没有……查了半天，原来是样本名末尾多了个空格。这种坑用上面的检查提前就能揪出来。

小技巧：用head(otu_table(ps))和head(tax_table(ps))看一眼前几行，比看整个对象更直观。

4.5 完整代码模板

下面是我每次新项目都会用的模板。你直接复制，改一下文件名就能用：

library(phyloseq)

# 第一步：读取数据
otu_mat <- read.csv("otu_table.csv", row.names = 1)
tax_mat <- read.csv("tax_table.csv", row.names = 1)
samp_df <- read.csv("sample_metadata.csv", row.names = 1)

# 第二步：转成正确格式
tax_mat <- as.matrix(tax_mat)
tax_mat <- apply(tax_mat, 2, as.character)

# 第三步：构建phyloseq组件
OTU <- otu_table(otu_mat, taxa_are_rows = TRUE)
TAX <- tax_table(tax_mat)
SAM <- sample_data(samp_df)

# 第四步：合并
ps <- phyloseq(OTU, TAX, SAM)

# 第五步：快速检查
ps
sum(taxa_sums(ps) == 0)
sum(sample_sums(ps) == 0)

如果你还在用旧版本的phyloseq（比如1.x），read.csv读东西可能会出现row.names不生效的问题。建议升级到3.x以上，或者直接用readr::read_csv后再手动设行名。

4.6 流程图：数据导入的核心逻辑

下面这张SVG图能让你一眼看清整个导入的流程和数据结构的关系：

你看，整个流程就是三个源头汇聚成一个对象。你做数据导入的时候，脑子里就保持这张图的结构。

4.7 常见问题与避坑

⚠ 警告：不要直接读BIOM格式时忘了指定参数

有些扩增子流程（如Qiime2）会输出.biom文件。用import_biom()可以读，但记得加上parseFunction = parse_taxonomy_greengenes或自建的注释解析函数，否则分类信息经常读成乱码。我曾经在这上面折腾了一个下午。

行名不对齐：用intersect检查一下共同名称。例如：intersect(rownames(otu_mat), rownames(tax_mat))。如果交集小于全集，那就是名称写错了。
样本名带特殊字符：减号、空格、点号都可能引发问题。我建议在前期数据清洗时就统一用下划线，比如“S01_Control”。
缺失分类信息：有些OTU可能注释到“Bacteria”就断了。这不是错误，但下游分析时你可能会丢失很多功能信息。提前用table(is.na(tax_table(ps)))看一下缺失情况。

4.8 小结一下

嗯，这一章内容其实就一个核心任务：把三张独立表格组装成一个phyloseq对象。你只需要三个步骤——读取、格式转换、合并——再加上一个快速检查。代码量很小，但这里面藏着很多经验性的细节：

始终检查行名/列名是否一致
分类表一定要转成字符矩阵
合并不报错不代表数据没问题，一定要用ps看一眼摘要

我刚开始做这个操作的时候，至少踩过五次“行名不匹配”的坑。后来我就养成了一个习惯：每次读表之后先用head扫一眼，然后再往下走。虽然多花了10秒，但能省掉后面复盘的一小时。你也试试看？