4. 数据导入

好,这一章我们来解决一个最实际的问题:怎么把手上那几张表格塞进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图能让你一眼看清整个导入的流程和数据结构的关系:

数据导入至phyloseq对象流程图 OTU / ASV 丰度表 行:OTU_ID,列:样本 分类注释表 行:OTU_ID,列:分类阶元 样本元数据 行:样本名,列:实验信息 合并:phyloseq(OTU, TAX, SAM) 自动对齐行名/列名 phyloseq 对象 包含 otu_table + tax_table + sample_data 快速检验:ntaxa、nsamples、零计数检查

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

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秒,但能省掉后面复盘的一小时。你也试试看?