第四章：探索性数据分析——基础统计、箱线图、直方图、密度图

拿到数据第一件事，别急着调参数。

我做数据分析这几年，最深刻的体会就是：先看数据长什么样。你连数据的脾气都没摸清楚，上来就咔咔跑模型，多半要翻车。

这一章，咱们就聊聊怎么用R语言给丰度数据做“体检”。说白了就是四板斧：基础统计量、箱线图、直方图、密度图。别看都是老掉牙的东西，想用透真不容易。

4.1 用基础统计量摸清数据底细

我第一次接触宏基因组丰度数据时，上来就画热图。结果图是挺漂亮，可有一半的数值都是零，根本没法看。

从那以后，我养成了一个习惯：先跑一轮summary()。

# 假设你的丰度数据长这样
abundance_data <- read.csv("otu_table.csv", row.names = 1)

# 先看整体
summary(abundance_data[, 1:5])  # 挑前几个样本看看

# 常用的统计量
library(psych)
describe(abundance_data)

# 重点关注几个指标
# 1. 零值占比
zero_ratio <- apply(abundance_data, 1, function(x) sum(x == 0) / length(x))
# 2. 极值检查
max_vals <- apply(abundance_data, 1, max, na.rm = TRUE)

到这里，你就能回答几个关键问题了：

• 数据是稀疏还是稠密？（零值比例高不高）
• 有没有离谱的异常值？（最大值比第三四分位数大几个数量级）
• 数据是不是偏态分布？（均值和中位数差多少）

4.2 箱线图：一眼看穿离群样本

箱线图这玩意儿，说白了就是把summary()的结果画出来。但视觉冲击力完全不一样。

我记得有一次分析肠道菌群数据，summary一切正常，结果画了箱线图才发现：有个样本的中位丰度是其他样本的3倍。后来查记录，是标记错了。

你想想看，要是没这一步，后面所有结论都白搭。

# 基础箱线图
boxplot(abundance_data, 
        main = "样本丰度分布概览",
        xlab = "样本编号", 
        ylab = "OTU丰度",
        col = "#3498db",
        outline = TRUE)  # 保留异常值

# 如果样本太多，转置一下看OTU分布
boxplot(t(abundance_data[, 1:20]), 
        main = "前20个OTU的丰度分布",
        las = 2)  # 标签旋转

# 配合ggplot2更灵活
library(ggplot2)
library(reshape2)

abundance_melt <- melt(abundance_data)
ggplot(abundance_melt, aes(x = variable, y = value)) +
  geom_boxplot(fill = "#2ecc71", outlier.color = "#e74c3c") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  labs(x = "样本", y = "丰度")

看箱线图重点关注什么？

• 箱子高度——样本间的变异程度
• 须的长度——极端值多不多
• 中点位置——中位数在箱子里的位置，决定了偏态方向
• 异常点——那些超出1.5倍四分位距的点，逐个检查

4.3 直方图：判断分布形态

直方图能告诉你数据的“形状”。是正态、是偏态、还是双峰？这个信息直接决定你后续用什么统计方法。

比如，你要做差异分析。数据是正态分布，可以用t-test；不是的话，得用Wilcoxon。选错了，结果不可靠。

# 单个样本的直方图
hist(abundance_data[, 1], 
     breaks = 30,  # 调整组数
     main = "样本1的丰度分布",
     xlab = "丰度",
     col = "#9b59b6",
     border = "white")

# 多图对比，用par布局
par(mfrow = c(2, 3))
for(i in 1:6) {
  hist(abundance_data[, i], 
       breaks = 20,
       main = paste("样本", i),
       xlab = "丰度",
       col = "#1abc9c")
}

# 或者用ggplot2分面
ggplot(abundance_melt, aes(x = value)) +
  geom_histogram(bins = 30, fill = "#f39c12", color = "white") +
  facet_wrap(~ variable, scales = "free") +
  labs(x = "丰度", y = "频数")

看直方图时，我会问自己三个问题：

1. 数据是不是严重右偏？（大多数丰度数据都是，因为很多OTU丰度很低）
2. 有没有明显的双峰？（可能暗示不同的菌群类型）
3. 零值是不是单独成了一个峰？（需要决定要不要做零值处理）

4.4 密度图：平滑版的直方图

密度图可以理解为直方图的平滑版本。它用核密度估计（KDE）画出一条光滑曲线，能更清晰地展示数据的分布趋势。

# 基础密度图
plot(density(abundance_data[, 1]), 
     main = "样本1的丰度密度曲线",
     xlab = "丰度",
     col = "#e74c3c",
     lwd = 2)

# 多条曲线叠加
plot(density(abundance_data[, 1]), 
     col = "#3498db", 
     lwd = 2,
     main = "多个样本的密度对比",
     xlim = c(0, 100))
lines(density(abundance_data[, 2]), col = "#e74c3c", lwd = 2)
lines(density(abundance_data[, 3]), col = "#2ecc71", lwd = 2)
legend("topright", 
       legend = paste("样本", 1:3),
       col = c("#3498db", "#e74c3c", "#2ecc71"),
       lwd = 2)

# 对数转换后的密度图
library(vegan)
# 假设做了对数转换
log_data <- decostand(abundance_data, method = "log")
plot(density(log_data[, 1]), 
     main = "log转换后的丰度分布",
     col = "#9b59b6",
     lwd = 2)

密度图最实用的地方，是你能直观对比不同组样本的分布。比如：

• 处理组和对照组的分布有没有明显错位？
• 对数转换后，数据是不是更接近正态？
• 不同样本的曲线形状是否一致？（不一致可能意味着批次效应）

4.5 图形结合：完整的数据体检流程

到这里，基础统计、箱线图、直方图、密度图，四样工具我们都过了一遍。但实际的EDA不能孤立地用。

我的标准流程是这样的:先跑summary看数字，发现异常后用箱线图定位，判断分布用直方图，细节呈现用密度图。四者配合起来，就能在几分钟内掌握数据的全貌。

4.6 实战：从原始数据到分析决策

拿一个真实的土壤微生物数据来说。

原始数据有5000个OTU、30个样本。走完上述四步，我发现了几个问题：

1. summary结果显示，有40%的OTU在所有样本中都是零——这部分得过滤掉。
2. 箱线图显示，样本S03的箱子远高于其他样本——查发现是DNA提取量不同，需作标准化。
3. 直方图显示，原始丰度是典型的幂律分布——后续差异分析必须用非参数方法。
4. 密度图显示，对数转换后，数据接近正态——可以放心用线性模型。

# 完整的EDA流程（组合拳）
eda_report <- function(data) {
  cat("===== 1. 基础统计 =====\n")
  print(apply(data, 2, summary))
  
  cat("\n===== 2. 零值比例 =====\n")
  zero_pct <- apply(data, 2, function(x) sum(x == 0) / length(x) * 100)
  print(sort(zero_pct, decreasing = TRUE)[1:10])
  
  cat("\n===== 3. 箱线图 =====\n")
  boxplot(data, outline = FALSE, col = "#3498db")
  
  cat("\n===== 4. 直方图（汇总） =====\n")
  hist(rowMeans(data), breaks = 30, col = "#2ecc71")
  
  cat("\n===== 5. 密度图（对数尺度） =====\n")
  log_data <- log1p(data)
  plot(density(rowMeans(log_data)), col = "#e74c3c", lwd = 2)
}

# 跑一遍
eda_report(abundance_data)

这个函数我每次拿到新数据都会先跑一次。里面的参数都是调过的，箱子颜色、组数、尺度，都是多年经验积累出来的。