正文

从癌症到罕见病基因富集分析实战 三步教你用GO和KEGG通路找出关键基因差异

/0 浏览量
0902

从癌症到罕见病基因富集分析实战 三步教你用GO和KEGG通路找出关键基因差异

说实话，我第一次接触基因富集分析的时候，整个人都是懵的。那些密密麻麻的基因名，什么GOTerm、KEGG pathway，看起来就像天书一样。但当你真正动手跑通一次之后，你会发现——这事儿其实没想象中那么可怕。

今天咱们就从头到尾，把一个完整的基因富集分析流程掰开揉碎讲清楚。不管你是研究癌症的，还是做罕见病基因的，这个套路都能用。

先搞明白：我们到底在干什么？

假设你现在手上有一份差异表达基因列表，比如从癌症样本和正常样本的RNA-seq数据里筛出来的500个显著变化的基因。你看着这500个名字，脑子里肯定有一万个问号：

“这些基因到底在干啥？” “它们跟什么生物学过程有关？” “能不能找到关键的信号通路？”

基因富集分析就是来回答这些问题的。它的基本思路很简单——不要只看单个基因，要看一群基因在干什么。

打个比方：你有一堆不同的人，你能通过观察他们从事的职业、活动的地点，来判断这是一个什么团队、在干什么事。基因也一样，如果一批差异基因都集中在”细胞凋亡”或者”PI3K-Akt信号通路”里，那说明这些过程在你的疾病模型里很重要。

第一步：准备你的差异基因数据

从哪里来？

差异表达基因通常来自转录组测序数据分析。常用的工具有DESeq2、edgeR、limma等。假设你已经用DESeq2跑出结果了，现在需要提取显著差异的基因。

这里用一个真实的例子。假设你的DESeq2分析结果长这样：

# 加载包
library(DESeq2)
library(dplyr)

# 读取count矩阵和样本信息
counts <- read.csv("count_matrix.csv", row.names = 1)
coldata <- read.csv("sample_info.csv")

# 创建DESeqDataSet
dds <- DESeqDataSetFromMatrix(
  countData = counts,
  colData = coldata,
  design = ~ condition
)

# 运行差异分析
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "tumor", "normal"))

# 筛选显著差异基因（FDR < 0.05，|log2FC| > 1）
sig_genes <- res %>%
  as.data.frame() %>%
  rownames_to_column(var = "gene_id") %>%
  filter(padj < 0.05, abs(log2FoldChange) > 1)

跑完这段代码，你会得到一个包含显著差异基因的表格。通常里面会有基因ID、log2FC、p值、padj等列。

基因ID的转换

这是很多人踩坑的地方。你的差异基因可能是Ensembl ID，比如”ENSG00000139618”，但做富集分析的时候，注释数据库通常用Gene Symbol（基因正式名称），比如”TP53”。

怎么转换？用biomaRt包：

library(biomaRt)

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

# 提取显著基因的Ensembl ID
ensembl_ids <- sig_genes$gene_id

# 转换ID
gene_info <- getBM(
  attributes = c("ensembl_gene_id", "external_gene_name", "gene_biotype"),
  filters = "ensembl_gene_id",
  values = ensembl_ids,
  mart = ensembl
)

# 看看转换结果
head(gene_info)

转换完之后，你需要检查一下转换效率。有些基因在数据库里找不到对应的Symbol，这是正常的。一般来说，转换率在80%以上就没问题。如果太低，可能需要检查你的物种或者数据库版本。

# 检查转换率
total_genes <- nrow(sig_genes)
converted_genes <- nrow(gene_info)
conversion_rate <- converted_genes / total_genes * 100

cat(sprintf("总差异基因数：%d\n转换成功数：%d\n转换率：%.1f%%\n", 
            total_genes, converted_genes, conversion_rate))

把转换好的基因Symbol保存下来，这就是你下一步分析的核心输入：

# 保存基因列表
write.table(gene_info$external_gene_name, 
            file = "significant_genes.txt", 
            quote = FALSE, 
            row.names = FALSE, 
            col.names = FALSE)

第二步：做GO功能富集分析

GO是什么？

基因本体（Gene Ontology，简称GO）是生物信息学里最常用的功能注释框架。它把基因功能分成三大类：

  • BP（Biological Process）：生物学过程，比如”细胞分裂”、”DNA修复”
  • CC（Cellular Component）：细胞组分，比如”细胞核”、”线粒体”
  • MF（Molecular Function）：分子功能，比如”ATP结合”、”激酶活性”

用clusterProfiler做GO分析

现在最主流的富集分析工具是R包clusterProfiler，它功能强大，文档也写得好。

library(clusterProfiler)
library(org.Hs.eg.db)  # 人类基因注释包

# 读取差异基因列表
gene_list <- read.table("significant_genes.txt", header = FALSE)
gene_symbol <- gene_list$V1

# 准备背景基因（所有检测到的基因，不只是显著的）
all_genes <- rownames(res)  # DESeq2结果里的所有基因
all_symbol <- mapIds(org.Hs.eg.db,
                     keys = all_genes,
                     column = "SYMBOL",
                     keytype = "ENSEMBL",
                     multiVals = "first")

# 构建带log2FC的基因向量（用于做ggsurv之类的图）
log2fc_vec <- setNames(res$log2FoldChange[all_genes], all_symbol)
log2fc_vec <- na.omit(log2fc_vec)

好，现在可以做GO富集分析了：

# GO富集分析
go_result <- enrichGO(
  gene = gene_symbol,
  OrgDb = org.Hs.eg.db,
  ont = "BP",          # 可以选BP/CC/MF，或者都做
  pAdjustMethod = "fdr",  # 用FDR校正
  pvalueCutoff = 0.05,
  qvalueCutoff = 0.05,
  readable = TRUE       # 把Ensembl ID转成基因名
)

# 看看结果
head(go_result)

输出结果长这样：

ID Description GeneRatio BgRatio pvalue p.adjust qvalue geneID Count
GO:0006915 apoptosis 45⁄500 120⁄20000 1.2e-10 3.5e-08 8.2e-08 TP53,BAX,… 45

每一行就是一个GO term。GeneRatio表示在你的差异基因里有多少个属于这个term，BgRatio表示在所有背景基因里有多少个属于这个term。如果GeneRatio远大于BgRatio，说明这个功能在你的基因集里显著富集。

可视化结果

# 条形图
barplot(go_result, showCategory = 20) + ggtitle("GO BP富集分析结果")

# 气泡图
dotplot(go_result, showCategory = 30) + ggtitle("GO富集气泡图")

气泡图里，点的大小代表GeneRatio（差异基因占比），颜色代表p.adjust值（越红越显著），X轴是GO term名称。一眼就能看出哪些通路最显著。

解读结果

假设你看到这些显著富集的term：

  • apoptosis（细胞凋亡）
  • cell cycle（细胞周期）
  • DNA repair（DNA修复）
  • regulation of cell proliferation（细胞增殖调控）

恭喜你，你找到方向了。这些过程跟癌症高度相关——癌细胞的核心特征就是不控制地增殖、逃避凋亡、积累DNA损伤。GO分析帮你从几百个基因里提炼出了关键的生物学故事。

第三步：做KEGG通路富集分析

KEGG是什么？

KEGG（Kyoto Encyclopedia of Genes and Genomes）是另一套流行的通路数据库。跟GO不一样，KEGG关注的是具体的信号通路和代谢通路，比如”p53信号通路”、”PI3K-Akt信号通路”、”MAPK信号通路”等。

KEGG富集分析

# KEGG富集分析
kegg_result <- enrichKEGG(
  gene = gene_symbol,
  organism = "hsa",      # 人类
  pAdjustMethod = "fdr",
  pvalueCutoff = 0.05,
  qvalueCutoff = 0.05,
  use_internal_data = FALSE  # 使用KEGG在线数据库
)

# 查看结果
head(kegg_result)

结果示例：

ID Description GeneRatio BgRatio pvalue p.adjust qvalue geneID Count
hsa04115 p53 signaling pathway 18⁄500 80⁄29664 2.3e-08 1.5e-06 3.2e-06 TP53,CDKN1A,… 18
hsa04151 PI3K-Akt signaling pathway 25⁄500 150⁄29664 1.1e-07 4.2e-06 8.5e-06 PIK3CA,AKT1,… 25

更高级的可视化：通路图

clusterProfiler还能把富集结果直接画到KEGG通路图上，这样你就能看到哪些基因在你的通路里：

# 选一个最显著的通路画图
p53_pathway <- kegg_result$ID[1]  # 最显著的通路

# 画通路图
plot_kegg_graph(p53_pathway, gene = gene_symbol, 
                organism = "hsa", 
                foldChange = log2fc_vec[gene_symbol])

基因在通路图上的颜色深浅代表表达变化的程度，红色是上调，蓝色是下调。这个图非常直观——你能看到p53通路里哪些节点被显著影响，哪些基因是协同变化的。

同时跑GO和KEGG：enricher函数

如果你觉得分别跑两次有点麻烦，或者你的基因列表是从其他分析来的（比如WGS筛选出的突变基因），可以用enricher函数做自定义富集：

# 准备基因集（每个term对应的基因列表）
gene_sets <- list(
  "apoptosis" = c("TP53", "BAX", "CASP3", "CASP9", "FAS", "TNFRSF10B"),
  "cell_cycle" = c("CDK4", "CDKN1A", "CCND1", "RB1", "MDM2"),
  "angiogenesis" = c("VEGFA", "ANGPT1", "TIE1", "KDR")
)

# 用enricher做自定义富集
enrich_result <- enricher(
  gene = gene_symbol,
  TERM2GENE = data.frame(
    term = rep(names(gene_sets), sapply(gene_sets, length)),
    gene = unlist(gene_sets)
  ),
  pAdjustMethod = "fdr",
  pvalueCutoff = 0.05
)

这个功能在研究罕见病时特别有用——你可以把自己领域已知的基因集放进去，看看你的差异基因里有没有 overlaps。

实战案例：从癌症到罕见病

讲了这么多，咱们来一个完整的实战。假设你在研究一种罕见病——脊髓性肌萎缩症（SMA），同时你也有癌症的数据可以做对比。

数据准备

# 罕见病SMA的差异化基因
sma_genes <- c("SMN1", "SMN2", "DCTN1", "KIF5A", "DYNC1H1", 
               "CNTNAP2", "NRXN1", "NLGN1", "SHANK3", "RBFOX1")

# 癌症（比如乳腺癌）的差异基因
cancer_genes <- c("TP53", "BRCA1", "BRCA2", "ERBB2", "PIK3CA", 
                  "PTEN", "EGFR", "MYC", "RB1", "CDH1",
                  "GATA3", "FOXA1", "ESR1", "MED1", "TCF3")

富集分析

# SMA的GO富集
sma_go <- enrichGO(
  gene = sma_genes,
  OrgDb = org.Hs.eg.db,
  ont = "BP",
  pAdjustMethod = "fdr",
  pvalueCutoff = 0.05,
  readable = TRUE
)

# 癌症的GO富集
cancer_go <- enrichGO(
  gene = cancer_genes,
  OrgDb = org.Hs.eg.db,
  ont = "BP",
  pAdjustMethod = "fdr",
  pvalueCutoff = 0.05,
  readable = TRUE
)

# 癌症的KEGG富集
cancer_kegg <- enrichKEGG(
  gene = cancer_genes,
  organism = "hsa",
  pAdjustMethod = "fdr",
  pvalueCutoff = 0.05
)

对比分析：找共同和独特的通路

# 提取两个GO结果的term
sma_terms <- unique(sma_go$Description)
cancer_terms <- unique(cancer_go$Description)

# 共同的通路
shared_terms <- intersect(sma_terms, cancer_terms)
cat("共同富集通路：", paste(shared_terms, collapse = ", "), "\n")

# SMA独有的通路
sma_unique <- setdiff(sma_terms, cancer_terms)
cat("SMA独有通路（前10个）：", paste(head(sma_unique, 10), collapse = ", "), "\n")

# 癌症独有的通路
cancer_unique <- setdiff(cancer_terms, sma_terms)
cat("癌症独有通路（前10个）：", paste(head(cancer_unique, 10), collapse = ", "), "\n")

结果解读

通过这样的对比，你可能会发现：

  • SMA和癌症都富集在”RNA splicing”（RNA剪接）相关的term里。这是有道理的——SMN1基因本身就是参与剪接体组装的，而癌症里也有很多剪接因子突变。
  • 癌症独有的通路主要是”cell cycle”、”p53 signaling”、”PI3K-Akt signaling”这些经典的肿瘤通路。
  • SMA独有的通路可能涉及”neuromuscular function”、”axon guidance”、”motor neuron development”等神经发育相关过程。

这个对比本身就很有意义——它帮你区分了疾病的共性机制和特异性机制。

避坑指南：这些错误我踩过

1. 背景基因集选错了

很多人直接用所有基因做背景，这是不对的。背景基因应该是你实验中能检测到的基因。比如在RNA-seq里，低表达的基因根本检测不到差异，不应该放进背景里。

# 正确的做法：用所有表达基因做背景
background_genes <- all_symbol  # 你测序能检测到的所有基因

2. 多重检验校正没做

富集分析同时检验几百上千个term，不做校正假阳性会爆炸。一定要用FDR（BH方法）校正，不要只看原始p值。

# clusterProfiler默认用fdr校正，但如果你手动算一定要记得
result$p.adjust <- p.adjust(result$pvalue, method = "fdr")

3. 基因ID混用

Ensembl ID、Gene Symbol、Entrez ID混着用是常见错误。同一份分析里统一用一种ID，需要转换就用biomaRt或clusterProfiler自带的readable参数。

4. 忽略基因数量太少的term

有些term只包含2-3个基因，就算全部差异也显得”富集”，这是假阳性。可以用minGSSize参数过滤：

go_result <- enrichGO(
  gene = gene_symbol,
  OrgDb = org.Hs.eg.db,
  ont = "BP",
  pAdjustMethod = "fdr",
  pvalueCutoff = 0.05,
  minGSSize = 10,      # 每个term至少10个基因
  maxGSSize = 500      # 每个term最多500个基因
)

最后：让结果说话

富集分析做完了，接下来就是解读。不要把结果扔在表格里就不管了，要真正去理解这些通路在你研究的问题里意味着什么。

一个实用的技巧是画venn图，看看不同条件之间的基因重叠：

library(VennDiagram)

# 比较SMA和癌症的差异基因重叠
venn.diagram(
  x = list(SMA = sma_genes, Cancer = cancer_genes),
  filename = "venn_diagram.png",
  height = 600,
  width = 600,
  col = c("blue", "red"),
  fill = c("blue", "red"),
  alpha = 0.5
)

如果发现有一些基因在SMA和癌症里都差异表达，那这些基因可能就是两个疾病的交叉靶点，值得深入研究。

基因富集分析不是终点，而是起点。它帮你从海量的数据里找到方向，真正的答案还需要你用实验去验证。但有了GO和KEGG的指引，你至少知道该往哪里走了。

希望这篇笔记能帮你少走弯路。有任何问题，随时交流。

-- 展开阅读全文 --