2. 数据导入与清洗：从原始数据到整洁数据框

大家好。这一章，我们正式开始动手处理数据。

宏基因组数据从测序平台出来之后，通常是各种格式的文件：OTU表、物种分类表、样本信息表。说白了，第一件事就是把它们读进R，变成我们能操控的数据框。

我见过不少初学者，读个文件都要折腾半天。其实没那么复杂，但细节确实不少。今天我们就来把这些事儿捋一遍。

2.1 数据导入：三大金刚函数

R里的数据导入函数很多，但常用的就这几个：read.table()、read.csv()，还有来自readr包的read_tsv()。它们各有脾性，用不对就容易出幺蛾子。

函数	默认分隔符	字符串处理	速度
read.table()	空白（空格、制表符）	默认将字符串转为因子	较慢
read.csv()	逗号	默认将字符串转为因子	较慢
read_tsv()	制表符	保留为字符，不自动转因子	快

我个人习惯，处理宏基因组数据时，首选read_tsv()。为什么？因为物种丰度表通常是制表符分隔的TSV文件，而且它不会自作主张把样本ID列变成因子。因子这个东西，有时候真的很坑。

典型用法

假设你有一个OTU表，叫做otu_table.tsv。它长这样：

OTU_ID   Sample1   Sample2   Sample3
OTU_1    123       456       789
OTU_2    234       567       890
OTU_3    345       678       901

用read_tsv()读它：

library(readr)
otu <- read_tsv("otu_table.tsv")

如果你想用read.table()，记得指定参数：

otu <- read.table("otu_table.tsv",
                  header = TRUE,
                  sep = "\t",
                  stringsAsFactors = FALSE)

嗯，这里要注意。很多人在read.table()里忘记写stringsAsFactors = FALSE。结果列名全变成了因子，后续操作处处碰壁。我在项目中遇到过好几次这种状况，排查起来挺心累的。

2.2 数据框操作：dplyr三板斧

数据读进来了，接下来就是清洗和整理。说实话，dplyr这个包真的太好用了。它的设计思路特别清晰：

• filter()：按行筛选，保留你想要的样本或物种
• select()：按列选，只要特定的变量
• mutate()：新增或修改列

你想想看，宏基因组数据的清洗过程，无非就是这三步来回组合。

2.2.1 filter() —— 给数据做减法

比如说，我们只想保留丰度总和大于1000的OTU：

library(dplyr)
otu_filtered <- otu %>%
  filter(rowSums(select(., -OTU_ID)) > 1000)

这里rowSums()计算了每个OTU在所有样本中的总和。select(., -OTU_ID)的意思是排除OTU_ID这一列，只计算数值列的和。

我曾经过度筛选，一下子把数据框筛成了空表。后来我学乖了，筛选前先用count()或summary()看看分布。

2.2.2 select() —— 只留关键变量

宏基因组数据经常有很多元数据列。有时候你只需要样本名和丰度值：

otu_selected <- otu %>%
  select(OTU_ID, Sample1, Sample2, Sample3)

如果你想剔除某列，用减号：

otu_no_meta <- otu %>%
  select(-Taxonomy, -Kingdom)

说白了，select()就是给数据框“瘦身”的利器。

2.2.3 mutate() —— 创造新变量

我们经常需要把原始计数转换成相对丰度。这时候mutate()就派上用场了：

otu_rel <- otu %>%
  mutate(across(starts_with("Sample"), ~ . / sum(.) * 100))

这个代码的意思是：对以"Sample"开头的所有列，每个值除以该列的总和，再乘以100，得到相对丰度百分比。

2.3 缺失值处理：别让NA毁了你的图

宏基因组数据里，缺失值（NA）是常客。测序深度不够、某些物种没检测到，都会产生NA。如果不处理，画图的时候就会出现奇怪的空白条。

常见的处理方式有三种：

• 删除：如果缺失值占比很小，直接删掉相关行或列
• 填充：用0、均值、中位数等填充
• 插补：用统计模型估计缺失值（宏基因组中较少用）

我个人的建议是，对于物种丰度数据，缺失通常意味着“未检测到”，所以用0填充是最合理的。我们来写一段代码：

otu_clean <- otu %>%
  mutate(across(everything(), ~ ifelse(is.na(.), 0, .)))

简单粗暴，但有效。注意，如果元数据里有NA，不要轻易填充0，那可能意味着样本信息丢失。

2.4 数据标准化：让不同样本可以比较

不同样本的测序深度不一样，直接比较原始计数是不公平的。标准化就是要解决这个问题。常见的标准化方法有：

• 相对丰度（百分比）：每个值除以该样本的总和，再乘以100
• CPM（Counts Per Million）：每个值除以该样本的总和，再乘以1,000,000
• TMM（Trimmed Mean of M-values）：来自edgeR包，考虑文库大小和组成

对于初学者，相对丰度是最直接的。我们刚刚用mutate()已经实现过了。来复习一下：

otu_rel_abund <- otu_clean %>%
  mutate(across(starts_with("Sample"), ~ . / sum(.) * 100))

标准化之后，就可以放心往下做了。画堆叠柱状图、热图，都不会被测序深度误导。

2.5 本章知识结构

为了方便大家理解，我用一张图整理了这个环节的思路：

2.6 本章小结

这一章讲的全是动手前的准备工作。数据导入用对函数，清洗用对方法，缺失值和标准化处理顺手，后面画图才能行云流水。

我自己的经验是，花在数据清洗上的时间，往往比画图本身还多。但这事急不得，一份整洁的数据框，是高质量可视化的地基。