4. 表达矩阵预处理

拿到GEO的原始表达矩阵后，你会发现行名全是探针ID（比如“1007_s_at”、“AFFX-…”）。

这些ID工程师看不懂，生物学家也看不懂。说白了，它们只是芯片上的坐标编号。我们需要把它转成基因Symbol——也就是人人都认识的基因名字，比如TP53、EGFR。

这一步我习惯称为“探针映射”。

为什么要做映射？不做行不行？

你想想看，差异表达分析的最终目的是找基因。如果你用探针ID做分析，结果里出现一串“1007_s_at上调”，论文审稿人直接懵了。

我的经验是：探针映射是表达矩阵预处理中最容易出坑的一步。我之前在TCGA项目里就因为映射不全，硬生生丢了30%的基因……做了两天的分析全白费。

常见的坑：
  • 探针ID和基因Symbol不是一一对应的——一个探针可能映射到多个基因
  • 一个基因可能被多个探针覆盖——出现重复基因名
  • 部分探针没有对应注释——直接变成NA，如果不处理后面会报错

核心思路：两种映射方法

我一般推荐两种方法，看你手上有啥资源。

两种主流映射方法：
方法优点缺点
使用biomaRt（在线查询）数据新、覆盖全、代码简洁依赖网络、速度慢、可能断连
使用平台注释（GPL文件）离线可用、速度快、100%匹配该芯片平台文件格式混乱、需要清洗

方法一：用biomaRt做映射

biomaRt是Bioconductor的明星包。我个人最喜欢用，因为它不用下载任何本地注释文件。

# 加载包
library(biomaRt)

# 连接到Ensembl数据库
ensembl <- useMart("ensembl", dataset="hsapiens_gene_ensembl")

# 你的探针ID列表（假设已经读入）
probe_ids <- rownames(expr_matrix)

# 执行查询
mapping <- getBM(
  attributes = c("affy_hg_u133_plus_2", "hgnc_symbol"),
  filters = "affy_hg_u133_plus_2",
  values = probe_ids,
  mart = ensembl
)

这里要注意：attributes和filters的名字取决于你的芯片平台。比如你在处理GPL570平台，就要用affy_hg_u133_plus_2。如果是其他平台，得先去biomaRt里查对应名字。

小技巧：如果不确定用什么filter，可以运行listFilters(ensembl)查看所有可用的过滤条件。找到了就复制粘贴。

映射结果长这样：

   affy_hg_u133_plus_2 hgnc_symbol
1            1007_s_at        MIR4640
2             1053_at        RFC2
3              117_at        HSPA6
4              1494_f_at       MT-ND6
5             1552256_a_at        C17orf82
6             1552257_a_at        LINC02432

嗯，这里要注意一个问题：有些探针映射不到基因Symbol，会返回空值。如果你直接合并原矩阵，这些空值行需要处理。

方法二：用平台注释文件

这个方法虽然老派，但我在很多离线项目中不得不这样做。比如客户的内网服务器不可能连Ensembl。

# 假设平台文件是GPL570.annot，从GEO下载的
gpl <- read.table("GPL570.annot", 
                  header = TRUE, 
                  sep = "\t", 
                  quote = "",
                  stringsAsFactors = FALSE)

# 只提取探针ID和基因Symbol列
gpl_clean <- gpl[, c("ID", "Gene.Symbol")]

# 有些平台会把基因Symbol写成一堆别名，比如"TP53 // LINC01234"
# 我们只取第一个
gpl_clean$Gene.Symbol <- gsub(" // .*", "", gpl_clean$Gene.Symbol)

# 合并到表达矩阵
expr_matrix$Gene.Symbol <- gpl_clean$Gene.Symbol[match(rownames(expr_matrix), gpl_clean$ID)]
注意：平台文件的格式各不一样。有的叫“Symbol”，有的叫“Gene Symbol”，还有的直接写在“Description”里。我曾经因为写错列名查了一整天bug……建议先head(gpl)看一眼结构。

去除重复基因：如何处理“一基因多探针”

映射完成后，你会发现很多基因对应多个探针。比如TP53可能对应5个探针。差异分析时不能有重复行名，所以必须去重。

去重的核心问题是：保留哪个探针？

我的习惯是按最大表达量保留。说白了，哪个探针信号最强，就代表哪个探针最可靠。这在芯片数据中很常见——表达量高的探针通常更稳定。

# 假设expr_matrix已经合并了Gene.Symbol列
# 按基因Symbol分组，取每组中表达量最大值对应的行
library(dplyr)

expr_final <- expr_matrix %>%
  group_by(Gene.Symbol) %>%
  slice_max(order_by = rowMeans(across(where(is.numeric))), n = 1) %>%
  ungroup()

# 把Gene.Symbol设为行名
rownames(expr_final) <- expr_final$Gene.Symbol
expr_final <- expr_final[, -ncol(expr_final)]  # 去掉Symbol列

也可以按中位数或方差来选。我偏好最大值，因为直观。如果你做的是时间序列或需要更稳健的估计，可以试试中位数。

补充：另一种去重方式——直接取平均 
# 按基因Symbol取所有探针的平均值
expr_avg <- aggregate(. ~ Gene.Symbol, data = expr_matrix, FUN = mean)

取平均适合那种“我不想丢掉任何信息”的场景。但注意：如果某个探针是噪声，取平均反而会把好探针拖下水。

映射后质量检查

别急着跑下游分析。映射后一定要做三个检查：

  1. 检查映射率：计算多少探针成功映射到基因Symbol。如果低于80%，说明你得检查映射方法或平台文件。
  2. 检查重复基因占比：如果重复基因太多（比如超过20%），可能意味着芯片设计本身就偏冗余，可以考虑取平均而不是最大值。
  3. 检查NA值：映射失败的行在Gene.Symbol列会显示NA。我建议直接删掉，不要硬补成“Unknown”。
# 映射率检查
mapped_ratio <- sum(!is.na(expr_matrix$Gene.Symbol)) / nrow(expr_matrix)
cat("映射率：", round(mapped_ratio * 100, 2), "%")

# 重复基因检查
dup_ratio <- sum(duplicated(expr_matrix$Gene.Symbol)) / nrow(expr_matrix)
cat("重复基因比例：", round(dup_ratio * 100, 2), "%")

本章知识体系

我把探针映射和去重的整个流程画成了一张图，方便你理解。

探针ID → 基因Symbol 映射流程 原始表达矩阵 行名：探针ID 选择映射方法 ① biomaRt（在线） ② 平台注释（离线） ③ 自定义注释文件 映射后矩阵 含Gene.Symbol列 去除重复基因 按最大值保留 按中位数保留 最终表达矩阵 行名：基因Symbol 映射率检查 > 80% 重复基因去除后无NA

一点经验之谈

我做探针映射时，最容易被坑的是芯片平台不匹配。比如你明明下了GPL570的数据，结果biomaRt里却用了GPL96的过滤器——那映射率直接掉到20%。

所以我会建议：下载数据时就顺手记下平台GPL号。映射前用这个号跑一遍listFilters确认一下。

技巧：如果你对映射率不放心，可以在合并后跑一个简单检查：table(is.na(expr_matrix$Gene.Symbol))。如果NA太多，别硬跑下游，先回去修映射方法。

另一个我常被问的问题是：去重时能不能保留所有探针？我直接说不建议。差异分析工具如limma、DESeq2都要求唯一行名。如果不排重，后面会报错duplicate rownames not allowed。你想想看，与其报错再回来debug，不如一开始就解决好。

← 上一章 📖 返回目录 下一章 →