第2章:把数据拽进R — OTU/ASV表、注释文件和元数据导入
做微生物组分析,数据源其实就三样东西:丰度表、注释文件、样本信息。
这三者就像拼图。少一块,你就没法往下走。
我记得刚入门时,拿到一个挺大的OTU表,稀里糊涂就传进R,结果行名对不上、样本ID乱码、注释信息全是乱的。折腾了两天才搞明白哪里出了问题。
所以这一章,我把常见的坑和正确姿势一块儿讲清楚。
2.1 读OTU/ASV表 — 别让行名列名卡住你
OTU表基本长这样:每一行是一个OTU(或ASV),每一列是一个样本。
表里是计数数据——也就是每个物种在对应样本里被检测到多少次。
最坑人的地方是:很多工具输出OTU表时会多一列分类注释,但你读进R之后,第一列反而成了行名。
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表格式对比:
| 工具 | 分隔符 | 带注释列 | 替换列名 |
|---|---|---|---|
| QIIME1 | 制表符 | 最后一列 | 否 |
| QIIME2 | 制表符 | 否 | 否 |
| mothur | 制表符 | 最后一列 | 否 |
| USEARCH | 制表符 | 否 | 否 |
| DADA2 | CSV/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)等层级做统计。
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会按字母顺序排——这会让你在画箱线图时,对照组和实验组的顺序莫名其妙地互换。
2.4 三表合一 — 让它们相互验证
数据都读进来之后,第一步就是做验证:
- 行名验证:OTU表的行名和注释表的行名要能一一对应。
- 列名验证:OTU表的列名(样本名)和元数据表的行名要能一一对应。
- 空值检查:确保没有缺失值或异常值。
我常用这样的脚本来做检查:
# 检查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")
现在你已经拿到了整齐的 list(otu, tax, meta) 结构。后面要做的,就是从这张清洁的底表出发,跑各种统计和画图了。