4、基因列表准备：差异表达基因筛选方法、基因ID转换（Entrez/Symbol/Ensembl）

好，我们接着往下走。做GO富集分析，最核心的输入是什么？没错，就是差异表达基因列表。但这个问题没你想得那么简单——我见过太多人随便拿个p-value < 0.05的列表就开始跑富集，结果填进去一堆垃圾基因，出来的富集结果自然不可信。

这一章我们主要解决两件事：怎么筛选出可靠的差异表达基因，以及怎么把不同数据库的基因ID统一起来。后者其实是更烦人的事，因为公共数据源往往给你的是Symbol，而富集工具要的是Entrez ID，或者反过来。嗯，咱们一步步拆。

4.1 差异表达基因筛选的常见方法

差异表达分析（Differential Expression Analysis）的工具很多，但筛选逻辑大同小异。我个人习惯用两个阈值：|log2FC| 和 adjusted p-value（或FDR）。为什么不用原始p-value？因为多重检验校正必不可少。

实际操作中，我一般会再加一个表达量过滤：比如在至少一半样本中CPM（counts per million） > 1。为什么？因为低表达的基因统计检验不稳定，p-value再小也容易假阳性。

4.2 基因ID转换——躲不掉的坑

你拿到的差异基因列表，可能是Entrez ID（数字），可能是Gene Symbol（如TP53），也可能是Ensembl ID（ENSG开头）。而大多数GO富集工具（比如clusterProfiler）默认用Entrez ID做富集。所以，转换这一步，我必须做。

举个例子，你从GEO下载了一个表达数据集，里面注释是Symbol。你想用R包“clusterProfiler”做GO分析，但它的enrichGO函数要求输入Entrez ID。怎么办？

我推荐的转换途径：

• R包：org.Hs.eg.db（人类）、org.Mm.eg.db（小鼠）等。最靠谱。
• 在线工具：DAVID、bioDBnet。适合一次性的小列表。
• Ensembl Biomart：官方接口，支持批量查询。

注意一个非常重要的事：ID映射不是一一对应的。一个Symbol可能对应多个Entrez ID（比如基因家族成员），一个Entrez ID也可能没有对应Symbol。我做过一次批量转换，发现10%的Symbol无法映射到Entrez ID——如果不处理这些缺失，结果会偏。

4.3 代码实现：从差异表到富集输入

这部分我直接给一段R示例脚本。你运行的时候记得改一下文件路径和物种对应的数据库包。

# 安装必要的包（如果没装的话）
# install.packages("BiocManager")
# BiocManager::install("clusterProfiler")
# BiocManager::install("org.Hs.eg.db")

library(clusterProfiler)
library(org.Hs.eg.db)

# 假设你已经有一个差异分析结果 data.frame，列包括：
#   GeneSymbol, logFC, pvalue, padj
# 读取你的差异基因列表
deg <- read.table("deg_results.txt", header=TRUE, stringsAsFactors=FALSE)

# 筛选差异基因：|logFC| >= 1 且 padj <= 0.05
deg_sig <- deg[abs(deg$logFC) >= 1 & deg$padj <= 0.05, ]

# 从Symbol转换为Entrez ID
# 使用 bitr 函数，注意多对一时会保留多行
gene_entrez <- bitr(deg_sig$GeneSymbol, 
                    fromType = "SYMBOL", 
                    toType = "ENTREZID", 
                    OrgDb = org.Hs.eg.db)

# 移除未成功映射的
gene_entrez <- gene_entrez[!is.na(gene_entrez$ENTREZID), ]

# 输出用于GO富集的Entrez ID向量
go_input_genes <- unique(gene_entrez$ENTREZID)

# 查看有多少基因成功转换
cat("原始差异基因数：", nrow(deg_sig), "\n")
cat("成功映射到Entrez ID的基因数：", length(go_input_genes), "\n")

# 保存备用
write.table(go_input_genes, "go_input_genes.txt", 
            row.names=FALSE, col.names=FALSE, quote=FALSE)

4.4 Ensembl ID 的转换

有些公共数据直接给Ensembl ID，比如“ENSG00000141510”（人的TP53）。转换方法和前面类似，只是fromType改为"ENSEMBL"。

# 假设你有一个Ensembl ID向量
ens_ids <- c("ENSG00000141510", "ENSG00000146648")

# 转换成Entrez ID
gene_entrez_from_ens <- bitr(ens_ids,
                              fromType = "ENSEMBL",
                              toType = "ENTREZID",
                              OrgDb = org.Hs.eg.db)

# 也可以同时转成Symbol
gene_multimap <- bitr(ens_ids,
                       fromType = "ENSEMBL",
                       toType = c("ENTREZID", "SYMBOL"),
                       OrgDb = org.Hs.eg.db)

4.5 本章知识结构图

下面这张SVG帮你理清基因列表准备的整体流程：从原始数据到最终富集输入。

4.6 实战要点总结

• 先表达量过滤，再差异筛选：去掉低表达基因能减少很多假阳性。
• ID转换后必须检查丢失率：如果丢失超过20%，分析下原因——可能是物种包不对，也可能是基因名写错了。
• 保留原始ID和转换后ID的对应表：富集结果出来之后，你肯定想回查具体是哪个基因富集了。没有映射表会很痛苦。
• 考虑上下调分开导入：有些富集分析工具（如clusterProfiler的enrichGO）可以指定up或down基因集，分别做富集能得到更细致的生物学解释。