第三章：数据导入与预处理——readr与data.table实战

数据导入这事，看起来简单对吧？实际上我见过太多分析翻车，根源都在数据读进来时就埋了雷。今天这章节，我们重点搞定两件事：用正确的工具把数据读进来，再在分析前把数据洗干净。

3.1 数据读取：两种主流方案

我个人习惯先看你数据量有多大。百万行以内，readr足够好用；超过这个量级，直接上data.table。

3.1.1 使用readr读取数据

readr是tidyverse家族的成员，最吸引我的一点——它不会自作主张把字符串转成因子。这个坑我踩过很多次，用base R的read.csv时，常常莫名其妙多了很多因子水平。

# 安装与加载
install.packages("readr")
library(readr)

# 读取OTU/ASV丰度表
otu_table <- read_delim("otutab.txt", delim = "\t")

# 读取物种注释表
taxonomy <- read_csv("taxonomy.csv")

# 检查数据类型
glimpse(otu_table)

为什么会推荐read_delim而非read.table？因为read_delim会自动推断数据类型，而且对日期格式的处理比base R智能很多。我在项目中处理过一列混合了“2023-01-15”和“2023/01/15”两种格式的数据，readr能直接识别成日期，base R则会全部读成字符。

特性	readr	base R	data.table
读取速度(1万行)	0.3秒	0.8秒	0.1秒
自动类型推断	是	否	有限
字符串转因子	否	默认是	否
内存占用	中等	高	低
返回对象类型	tibble	data.frame	data.table

3.1.2 使用data.table读取大文件

如果你面对的是几G的宏基因组丰度表，fread是你的救星。说实话，我第一次用fread读一个2.3G的文件，只用了11秒——当时base R已经跑了快两分钟还没读完。

# 安装与加载
install.packages("data.table")
library(data.table)

# fread读取 - 语法极其简洁
otu_dt <- fread("otutab_large.txt")

# 指定关键参数
otu_dt <- fread(
  "otutab_large.txt",
  sep = "\t",
  header = TRUE,
  nThread = 4,           # 并行读取，多核CPU好物
  na.strings = c("NA", "na", "", "N/A", "."),  # 统一缺失值表示
  skip = 0               # 跳过空行
)

3.2 缺失值处理：别让NA毁了你

我记得第一次做微生物组分析，一个样本的丰度表里全是NA，结果下游的alpha多样性直接报错。后来我养成了习惯，数据一读进来就先检查缺失值。

缺失值在物种丰度数据里常见的原因：

• 测序深度不够，某些物种没被检测到
• PCR扩增失败，导致某个样本整行缺失
• 数据合并时的对齐错误

# 检查缺失值的几种方式
# 方法1：看看整体情况
sum(is.na(otu_table))

# 方法2：按列统计
colSums(is.na(otu_table))

# 方法3：按行统计
rowSums(is.na(otu_table))

# 方法4：快速查看哪些列有缺失
library(dplyr)
otu_table %>%
  summarise(across(everything(), ~ sum(is.na(.)))) %>%
  glimpse()

处理策略，说白了就三种：删除、填补、标记。怎么选？

缺失比例	推荐策略	适用场景
< 1%	填补为0或均值	低丰度物种偶尔缺失
1% - 5%	按样本或按物种填补	技术重复中个别缺失
5% - 20%	考虑删除该行/列	测序批次效应导致
> 20%	删除特征或样本	质量差的样本或稀有物种

# 常用的缺失值处理代码

# 方案A：填补为0（丰度数据最常用）
otu_clean <- otu_table
otu_clean[is.na(otu_clean)] <- 0

# 方案B：删除含有缺失的行
otu_clean <- na.omit(otu_table)

# 方案C：删除缺失比例超过20%的列
threshold <- 0.2 * nrow(otu_table)
cols_to_keep <- colSums(is.na(otu_table)) <= threshold
otu_clean <- otu_table[, cols_to_keep]

# 方案D：按物种均值填补
library(dplyr)
otu_clean <- otu_table %>%
  mutate(across(where(is.numeric), ~ ifelse(is.na(.), mean(., na.rm = TRUE), .)))

3.3 过滤低丰度特征

这个问题我琢磨了很久。低丰度特征为什么需要过滤？说白了，那些只出现一两次的OTU/ASV，很可能是测序错误或者PCR产生的嵌合体。留着它们，不仅增加多重假设检验的负担，还会稀释真实信号。

常用的过滤标准：

• 总丰度阈值：在所有样本中的总reads数 < 10
• 出现频率阈值：在 < 20% 的样本中出现
• 平均丰度阈值：平均相对丰度 < 0.01%

# 过滤低丰度特征

# 假设OTU表：行是OTU，列是样本
# 先计算每个OTU的总丰度
otu_sum <- rowSums(otu_clean)

# 过滤：总丰度 < 10 的OTU去掉
otu_filtered <- otu_clean[otu_sum >= 10, ]

# 过滤：在少于20%样本中出现的OTU去掉
presence <- rowSums(otu_clean > 0)
min_samples <- 0.2 * ncol(otu_clean)
otu_filtered <- otu_clean[presence >= min_samples, ]

# 组合过滤：两者都满足才保留
keep <- (otu_sum >= 10) & (presence >= min_samples)
otu_filtered <- otu_clean[keep, ]

cat("过滤前OTU数量:", nrow(otu_clean), "\n")
cat("过滤后OTU数量:", nrow(otu_filtered), "\n")
cat("保留比例:", round(nrow(otu_filtered)/nrow(otu_clean)*100, 1), "%\n")

你想想看，如果读了上百个样本，每个样本有几千个OTU，不先过滤的话，后续的分析就像在垃圾堆里找宝藏——累得很。我曾经在一个项目里，原始OTU表有15,632个OTU，过滤后只剩3,421个。但这3,421个包含了98%的总reads数，也就是说我们去掉了大量无意义的噪音。

3.4 完整的数据预处理流程

把前面的步骤串起来，就是一套完整的预处理流水线：

# 完整的物种丰度数据预处理流程
library(data.table)
library(dplyr)

# 步骤1：读取数据
cat("[1/4] 读取数据...\n")
otu_raw <- fread("otutab.txt", sep = "\t", header = TRUE)

# 步骤2：检查数据质量
cat("[2/4] 检查缺失值...\n")
cat("缺失值总数:", sum(is.na(otu_raw)), "\n")

# 步骤3：处理缺失值（填补为0）
otu_raw[is.na(otu_raw)] <- 0

# 步骤4：过滤低丰度特征
cat("[3/4] 过滤低丰度特征...\n")
otu_sum <- rowSums(otu_raw[, -1])  # 假设第一列是OTU ID
presence <- rowSums(otu_raw[, -1] > 0)
threshold_sum <- 10
threshold_freq <- 0.2 * (ncol(otu_raw) - 1)

keep <- (otu_sum >= threshold_sum) & (presence >= threshold_freq)
otu_clean <- otu_raw[keep, ]

# 步骤5：输出统计信息
cat("[4/4] 完成！\n")
cat("原始OTU数:", nrow(otu_raw), "\n")
cat("清理后OTU数:", nrow(otu_clean), "\n")
cat("过滤比例:", round((1 - nrow(otu_clean)/nrow(otu_raw))*100, 1), "%\n")

# 保存清理后的数据
fwrite(otu_clean, "otutab_clean.txt", sep = "\t")