# 2、物种丰度数据导入:读取OTU/ASV表、分类学注释文件、样本元数据

第2章:把数据拽进R — OTU/ASV表、注释文件和元数据导入

做微生物组分析,数据源其实就三样东西:丰度表、注释文件、样本信息。

这三者就像拼图。少一块,你就没法往下走。

我记得刚入门时,拿到一个挺大的OTU表,稀里糊涂就传进R,结果行名对不上、样本ID乱码、注释信息全是乱的。折腾了两天才搞明白哪里出了问题。

所以这一章,我把常见的坑和正确姿势一块儿讲清楚。

2.1 读OTU/ASV表 — 别让行名列名卡住你

OTU表基本长这样:每一行是一个OTU(或ASV),每一列是一个样本。

表里是计数数据——也就是每个物种在对应样本里被检测到多少次。

最坑人的地方是:很多工具输出OTU表时会多一列分类注释,但你读进R之后,第一列反而成了行名。

个人习惯:我每次读OTU表都会指定 row.names = 1,确保第一列变成行名。

示例:

library(tidyverse)

# 读经典OTU表
otu_raw <- read.table(
  "otufile/otutab.txt",
  header = TRUE,
  row.names = 1,
  sep = "\t",
  check.names = FALSE
)

注意check.names = FALSE这个选项。

为什么加它?因为R默认会把列名里的「-」改成「.」,我不希望样本名被随意篡改。

我遇到过一次:Illumina给的样本编号里带横线,R自动把"Sample-1"改成了"Sample.1"。后面做元数据匹配时死活对不上,排查了大半天才发现是这里出了问题。

注意事项:列名(样本名)和行名(OTU/ASV编号)不要包含特殊字符,如空格、引号、括号。建议只用字母、数字、下划线。

再看看常用工具生成的OTU表格式对比:

工具分隔符带注释列替换列名
QIIME1制表符最后一列否
QIIME2制表符否否
mothur制表符最后一列否
USEARCH制表符否否
DADA2CSV/TSV否否

如果是QIIME1的输出,我一般会先做一步处理:把最后一列的注释拆出去。

# QIIME1 OTU表的处理
otu_no_tax <- otu_raw[, -ncol(otu_raw)]
# 如果最后一列名称包含 "taxonomy" 或 "Consensus"
otu_no_tax <- select(otu_raw, -matches("taxonomy|Consensus"))

2.2 读分类学注释文件 — 一条链上搞定多级注释

分类学注释文件通常包含两列:OTU/ASV编号和注释字符串。

注释字符串一般是分级的,用分号隔开。比如:

OTU001  d__Bacteria; p__Firmicutes; c__Clostridia; o__Clostridiales
OTU002  d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales

我会把这些注释拆开成单独的列,方便后续按门(Phylum)、纲(Class)、目(Order)等层级做统计。

关键点:R中可以用 separate() 把一列拆成多列,或者用正则表达式提取指定层级。我个人倾向 separate(),因为更直观。
library(tidyr)

# 读注释文件
tax_raw <- read.table(
  "otufile/taxonomy.txt",
  header = TRUE,
  row.names = 1,
  sep = "\t",
  stringsAsFactors = FALSE
)

# 假设注释在 tax 列,使用分号作为分隔符
tax_sep <- tax_raw %>%
  separate(tax, into = c("Domain","Phylum","Class","Order","Family","Genus","Species"),
           sep = ";", extra = "merge", fill = "right")

# 清理多余前缀
tax_clean <- tax_sep %>%
  mutate(across(everything(), ~ str_remove(.x, "^[a-z]__")))

注意extra = "merge"和fill = "right"这两个参数。

它们的组合能处理注释字符串长度不一致的问题——比如有的物种没有种水平的注释,就会自动留空。

避坑指南:我曾经接过一份注释文件,里面有很多 unclassified 和 uncultured 标记。如果你不想把这些丢弃,请保留它们,否则会导致丰度总和与原始不一致。

2.3 读样本元数据 — 绑定变量才有分析

元数据(metadata)就是你实验设计的表格。它告诉你:哪个样本是处理组、哪个是对照组、哪个是时间点1、哪个是时间点2。

没有元数据,OTU表只是一堆数字,没法做任何有意义的分组分析。

元数据示例(TSV格式):

sampleID  group   day   treatment
S1        Control  0     None
S2        Control  0     None
S3        Treatment 7   Antibiotic
S4        Treatment 7   Antibiotic

读法也很简单:

meta <- read.table(
  "metadata/metadata.txt",
  header = TRUE,
  row.names = 1,
  sep = "\t",
  stringsAsFactors = FALSE
)

# 注意:把分组变量转为因子
meta$group <- factor(meta$group, levels = c("Control", "Treatment"))
meta$day <- as.numeric(meta$day)

为什么要把group转成因子?

因为后续做统计分析和可视化时,R会根据因子水平确定顺序。如果不设置levels,R会按字母顺序排——这会让你在画箱线图时,对照组和实验组的顺序莫名其妙地互换。

注意:元数据表格中,样本ID(第一列)必须和OTU表的列名完全一致。包括大小写、横线、下划线。我建议在导入之前,先用Excel或文本编辑器统一一下命名规则。

2.4 三表合一 — 让它们相互验证

数据都读进来之后,第一步就是做验证:

  1. 行名验证:OTU表的行名和注释表的行名要能一一对应。
  2. 列名验证:OTU表的列名(样本名)和元数据表的行名要能一一对应。
  3. 空值检查:确保没有缺失值或异常值。

我常用这样的脚本来做检查:

# 检查OTU行名与注释行名的一致性
otu_in_tax <- sum(rownames(otu_no_tax) %in% rownames(tax_clean))
cat("OTU匹配注释数:", otu_in_tax, "/", nrow(otu_no_tax), "\n")

# 检查样本身份
samples_in_meta <- sum(colnames(otu_no_tax) %in% rownames(meta))
cat("样本匹配元数据数:", samples_in_meta, "/", ncol(otu_no_tax), "\n")

# 提取一致样本
common_samples <- intersect(colnames(otu_no_tax), rownames(meta))
otu_final <- otu_no_tax[, common_samples, drop = FALSE]
meta_final <- meta[common_samples, , drop = FALSE]

这里有个细节:不要在OTU表和注释表中保留完全匹配的行drop = FALSE,否则当结果只有1行或1列时,R会自动降维变成向量。

小技巧:做完数据验证后,可以把三者的维度信息汇总打印出来,方便后续查阅。

2.5 整个流程的整合 — 一个函数打通

如果你经常重复这个工作流,我个人建议把它封装成一个函数:

read_microbiome_data <- function(otu_file, tax_file, meta_file, sep = "\t") {
  
  # 读入三个表
  otu <- read.table(otu_file, header = TRUE, row.names = 1,
                    sep = sep, check.names = FALSE)
  tax <- read.table(tax_file, header = TRUE, row.names = 1,
                    sep = sep, stringsAsFactors = FALSE)
  meta <- read.table(meta_file, header = TRUE, row.names = 1,
                     sep = sep, stringsAsFactors = FALSE)
  
  # 如果OTU表最后一列是注释,就剔除
  if (grepl("taxonomy|Consensus", tail(names(otu), 1), ignore.case = TRUE)) {
    otu <- otu[, -ncol(otu)]
  }
  
  # 提取一致样本
  common <- intersect(colnames(otu), rownames(meta))
  otu <- otu[, common, drop = FALSE]
  meta <- meta[common, , drop = FALSE]
  
  # 提取一致OTU
  common_otu <- intersect(rownames(otu), rownames(tax))
  otu <- otu[common_otu, , drop = FALSE]
  tax <- tax[common_otu, , drop = FALSE]
  
  return(list(otu = otu, tax = tax, meta = meta))
}

# 使用示例
data <- read_microbiome_data("data/otutab.txt", "data/taxonomy.txt", "data/metadata.txt")
核心思维:数据导入不是简单的读文件,而是在做结构对齐和质量控制。每次分析前,务必确认三张表之间没有信息断裂。
第2章知识框架:文件读取与数据对齐 OTU/ASV 丰度表 (行=OTU, 列=样本) 分类学注释文件 (OTU→各级分类) 样本元数据 (实验分组与条件) 读入原始数据 分隔符、行名列名配置 注释拆分与清理 separate() 按分号拆分 元数据格式化 因子顺序、数值类型 三表验证与对齐 行名→行名,列名→行名 list(otu, tax, meta) 数据类型匹配:字符/因子/数值

现在你已经拿到了整齐的 list(otu, tax, meta) 结构。后面要做的,就是从这张清洁的底表出发,跑各种统计和画图了。

一句话总结:数据导入花20分钟认真对齐,能省后续两天排查bug的时间。我踩过的坑,你就不用再踩了。