第四章：表达矩阵质量控制

拿到RNA-seq的表达矩阵后，很多人急着开始差异分析、建模。我踩过这个坑——直接用原始count矩阵跑模型，结果模型训练时异常波动、验证集指标忽高忽低。后来才发现，是几个低质量样本污染了整体分布。

质量控制不是可有可无的步骤。它决定了后续所有分析的可靠性。这一章，我带你过一遍三个核心检查：样本相关性热图、PCA异常样本检测、library size分布。三个动作做完，你才能放心地把矩阵喂给模型。

4.1 样本相关性热图

样本相关性热图，说白了就是看所有样本之间是不是「抱团」。生物学重复之间应该高度相关（通常Pearson r > 0.9），不同处理组之间可以低一些。如果同一组的两个样本相关性只有0.6，那肯定有问题——可能是样本污染、测序深度差异过大，或者index hopping。

我习惯用 cor() 直接计算样本间相关性，再用 pheatmap 可视化。来个例子：

# 假设 count_matrix 是基因×样本的原始count矩阵
# 先做logCPM标准化，避免高表达基因主导
library(edgeR)
logcpm <- cpm(count_matrix, log=TRUE)

# 计算样本间Pearson相关
cor_mat <- cor(logcpm, method = "pearson")

# 热图展示
library(pheatmap)
pheatmap(cor_mat,
         cluster_rows = TRUE,
         cluster_cols = TRUE,
         display_numbers = TRUE,
         color = colorRampPalette(c("navy", "white", "firebrick3"))(50),
         main = "Sample Correlation Heatmap")

避坑指南： 我曾经在一个公开数据集上看到样本相关性全部>0.99，以为是质量好，结果发现是因为使用了没有去除低表达基因的原始count，相关性被大量零值拉高。后来我坚持过滤后（至少10个样本中CPM>1）再算相关性，才得到真实结果。

4.2 PCA异常样本检测

PCA把高维基因表达投影到前几个主成分。正常的样本应该按实验条件聚集。如果你发现某个样本孤零零在散点图边缘，很可能是个「 outlier 」。我通常在PCA图上标出样本名或编号，眼尖一点直接能看到。

# 继续用上面的 logcpm
pca_res <- prcomp(t(logcpm), center = TRUE, scale. = FALSE)

# 提取前两个PC
pca_df <- as.data.frame(pca_res$x[, 1:2])
pca_df$group <- sample_group  # 你需要事先定义的分组向量

library(ggplot2)
ggplot(pca_df, aes(PC1, PC2, color = group, label = rownames(pca_df))) +
  geom_point(size = 3) +
  geom_text(vjust = 1.5, size = 3) +
  stat_ellipse(level = 0.9) +
  theme_minimal() +
  labs(title = "PCA of samples colored by group")

判断标准： 至少要看两点：
1. 某个样本是否落在9成置信椭圆外。
2. 生物学重复是否没有被前两个主成分分开（如果分开但组间混在一起是合理的，但组内分离说明重复性差）。

还有一个实用的指标：查看每个样本对PC贡献的杠杆值（leverage）。具体可以算每个样本在PC1-PC5上的平方和，超过平均值的2倍值得警惕。

# 计算杠杆值
leverage <- rowSums(pca_res$x[, 1:5]^2) / sum(pca_res$sdev[1:5]^2)
threshold <- 2 * mean(leverage)
outliers <- which(leverage > threshold)
cat("Potential outliers:", names(outliers))

4.3 Library size分布检查

Library size（文库大小）就是每个样本的总read count。不同样本之间library size差异过大，会导致标准化后某些样本的基因表达被人为压低或抬高。我一般画一个箱线图+密度图，看一眼就清楚。

# 原始count矩阵列求和
lib_size <- colSums(count_matrix)
lib_df <- data.frame(sample = names(lib_size), lib_size = lib_size / 1e6)

# 箱线图展示分布
ggplot(lib_df, aes(x = "Library Size", y = lib_size)) +
  geom_boxplot(outlier.colour = "red") +
  geom_jitter(width = 0.1, alpha = 0.6) +
  labs(y = "Library Size (millions)", title = "Library Size Distribution") +
  theme_minimal()

# 或者按分组看
lib_df$group <- sample_group
ggplot(lib_df, aes(x = group, y = lib_size, fill = group)) +
  geom_boxplot() +
  geom_jitter(width = 0.15, size = 2) +
  labs(y = "Library Size (millions)") +
  theme_minimal()

一般来说，如果某个样本的library size与同组其他样本相差超过3倍，这个样本就需要打上问号。但也有项目故意从浅测序开始，这时候结合检测基因数（至少1 count的基因数量）一起看。我做过一个项目，有个样本library size只有中值的1/5，但基因数却只少了30%，说明测序偏向某些转录本，质量有问题。

4.4 综合判断与标记

三个指标要结合起来看，不能单一决定。比如某个样本PCA不异常，但library size极低，我会保留它但标记为「低深度」。在后续差异分析时对这类样本给予更高权重或剔除。我常用的方法是建立一个质量评分：

# 样本质量打分（简易版）
quality_score <- data.frame(
  sample = names(lib_size),
  lib_size_flag = ifelse(lib_size < median(lib_size)*0.5, 1, 0),
  pca_outlier = ifelse(sample %in% names(outliers), 1, 0),
  cor_outlier = ifelse(sample %in% low_cor_samples, 1, 0)
)
quality_score$total_flag <- rowSums(quality_score[,2:4])
quality_score$decision <- ifelse(quality_score$total_flag >= 2, "remove",
                                 ifelse(quality_score$total_flag == 1, "flag", "pass"))

经过这三道关卡，你留下的表达矩阵才是干净的。我敢保证，跳过QC直接跑模型的人，十有八九会在调参时被异常样本折腾得怀疑人生。质量控制不是浪费时间，是给模型训练铺一条平坦的路。