蛋白质通路富集分析¶

方法一：ORA（过表示分析）= "数人头"
──────────────────────────────
输入：一组差异蛋白质列表（如100个上调蛋白）
原理：看这100个蛋白是不是"扎堆"出现在某个通路
统计：超几何检验（Fisher精确检验）
缺点：需要设阈值筛选差异蛋白，信息损失

方法二：GSEA（基因集富集分析）= "看趋势"
──────────────────────────────
输入：所有蛋白质的排序列表（按倍数变化排序）
原理：看某个通路的蛋白是集中在列表头部还是尾部
统计：富集得分ES（Enrichment Score）
优点：不需要设阈值，利用全部蛋白的信息

白话比喻：
ORA = 期末考试只看"不及格名单"在哪些班扎堆
GSEA = 看全年级成绩排名，某个班的学生是不是集中在前面/后面

# 安装clusterProfiler及相关包
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install(c(
    "clusterProfiler",  # 富集分析核心包
    "org.Hs.eg.db",     # 人类基因注释数据库
    "DOSE",             # 疾病富集分析
    "enrichplot",       # 富集结果可视化
    "pathview",         # KEGG通路可视化
    "ReactomePA"        # Reactome通路分析
))

library(clusterProfiler)  # 加载富集分析包
library(org.Hs.eg.db)     # 加载人类注释库
library(enrichplot)       # 加载可视化包
library(ggplot2)          # 加载画图包

#!/usr/bin/env Rscript
# ORA过表示分析（GO + KEGG）

# ========== 准备输入数据 ==========
# 读取limma差异分析结果
deg <- read.csv("limma_DEP_results.csv", row.names = 1)

# 筛选差异蛋白质
up_proteins <- rownames(deg)[deg$logFC > 1 & deg$adj.P.Val < 0.05]  # 上调蛋白
down_proteins <- rownames(deg)[deg$logFC < -1 & deg$adj.P.Val < 0.05]  # 下调蛋白
all_deg <- c(up_proteins, down_proteins)  # 所有差异蛋白

cat("上调蛋白:", length(up_proteins), "\n")
cat("下调蛋白:", length(down_proteins), "\n")

# 基因名转Entrez ID（GO/KEGG需要Entrez ID）
gene_ids <- bitr(
    all_deg,
    fromType = "SYMBOL",  # 输入类型：基因Symbol
    toType = "ENTREZID",  # 输出类型：Entrez ID
    OrgDb = org.Hs.eg.db  # 使用人类注释库
)
cat("成功映射:", nrow(gene_ids), "个基因\n")

# 背景基因集（所有检测到的蛋白质）
bg_genes <- bitr(
    rownames(deg),
    fromType = "SYMBOL",
    toType = "ENTREZID",
    OrgDb = org.Hs.eg.db
)

# ========== GO富集分析 ==========
# 分三个子本体：BP（生物学过程）、CC（细胞组分）、MF（分子功能）
go_bp <- enrichGO(
    gene = gene_ids$ENTREZID,  # 差异基因的Entrez ID
    universe = bg_genes$ENTREZID,  # 背景基因
    OrgDb = org.Hs.eg.db,  # 注释数据库
    ont = "BP",  # 子本体：生物学过程
    pAdjustMethod = "BH",  # p值校正方法
    pvalueCutoff = 0.05,  # p值阈值
    qvalueCutoff = 0.05,  # q值阈值
    readable = TRUE  # 把Entrez ID转回基因名
)
cat("GO-BP显著通路:", nrow(as.data.frame(go_bp)), "\n")

go_cc <- enrichGO(gene = gene_ids$ENTREZID, universe = bg_genes$ENTREZID,
                  OrgDb = org.Hs.eg.db, ont = "CC", readable = TRUE)  # 细胞组分
go_mf <- enrichGO(gene = gene_ids$ENTREZID, universe = bg_genes$ENTREZID,
                  OrgDb = org.Hs.eg.db, ont = "MF", readable = TRUE)  # 分子功能

# ========== KEGG富集分析 ==========
kegg <- enrichKEGG(
    gene = gene_ids$ENTREZID,  # 差异基因
    universe = bg_genes$ENTREZID,  # 背景基因
    organism = "hsa",  # 物种代码：hsa=人类
    pvalueCutoff = 0.05  # p值阈值
)
cat("KEGG显著通路:", nrow(as.data.frame(kegg)), "\n")

# ========== 可视化 ==========
# 1. 气泡图（GO-BP）
p1 <- dotplot(go_bp, showCategory = 15, title = "GO Biological Process") +
    theme(axis.text.y = element_text(size = 10))  # 调整字号
ggsave("go_bp_dotplot.png", p1, width = 10, height = 8, dpi = 300)

# 2. 柱状图（KEGG）
p2 <- barplot(kegg, showCategory = 15, title = "KEGG Pathway") +
    theme(axis.text.y = element_text(size = 10))
ggsave("kegg_barplot.png", p2, width = 10, height = 8, dpi = 300)

# 3. 网络图（GO term关联）
p3 <- emapplot(pairwise_termsim(go_bp), showCategory = 20)
ggsave("go_bp_emaplot.png", p3, width = 10, height = 10, dpi = 300)

# 4. 树状图
p4 <- treeplot(pairwise_termsim(go_bp), showCategory = 30)
ggsave("go_bp_treeplot.png", p4, width = 14, height = 10, dpi = 300)

cat("所有图片已保存\n")

#!/usr/bin/env Rscript
# GSEA基因集富集分析（不需要设阈值）

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

# ========== 准备排序列表 ==========
deg <- read.csv("limma_DEP_results.csv", row.names = 1)

# 转换基因ID
gene_map <- bitr(rownames(deg), fromType = "SYMBOL", toType = "ENTREZID",
                 OrgDb = org.Hs.eg.db)

# 合并
deg_mapped <- merge(deg, gene_map, by.x = "row.names", by.y = "SYMBOL")

# 创建排序向量：用log2FC排序
gene_list <- deg_mapped$logFC  # 排序指标：log2倍数变化
names(gene_list) <- deg_mapped$ENTREZID  # 名称：Entrez ID
gene_list <- sort(gene_list, decreasing = TRUE)  # 降序排列（上调在前）

cat("排序列表长度:", length(gene_list), "\n")
cat("范围:", range(gene_list), "\n")

# ========== GO-GSEA ==========
gsea_go <- gseGO(
    geneList = gene_list,  # 排序的基因列表
    OrgDb = org.Hs.eg.db,  # 注释数据库
    ont = "BP",  # 生物学过程
    minGSSize = 10,  # 最小基因集大小
    maxGSSize = 500,  # 最大基因集大小
    pvalueCutoff = 0.05,  # p值阈值
    eps = 0  # 允许精确到0的p值
)
cat("GO-GSEA显著结果:", nrow(as.data.frame(gsea_go)), "\n")

# ========== KEGG-GSEA ==========
gsea_kegg <- gseKEGG(
    geneList = gene_list,  # 排序的基因列表
    organism = "hsa",  # 人类
    minGSSize = 10,
    maxGSSize = 500,
    pvalueCutoff = 0.05,
    eps = 0
)
cat("KEGG-GSEA显著结果:", nrow(as.data.frame(gsea_kegg)), "\n")

# ========== GSEA可视化 ==========
# 1. GSEA山脊图（多通路对比）
p1 <- ridgeplot(gsea_go, showCategory = 15) +
    labs(title = "GSEA Ridge Plot - GO BP")
ggsave("gsea_ridgeplot.png", p1, width = 12, height = 8, dpi = 300)

# 2. 经典GSEA曲线（单通路）
if (nrow(as.data.frame(gsea_go)) > 0) {
    p2 <- gseaplot2(gsea_go, geneSetID = 1:3, pvalue_table = TRUE)  # 前3个通路
    ggsave("gsea_curve.png", p2, width = 10, height = 8, dpi = 300)
}

cat("GSEA可视化完成\n")

#!/usr/bin/env python3
"""Python版富集分析（使用gseapy包）"""

# pip install gseapy
import gseapy as gp  # GSEA Python实现
import pandas as pd  # 数据处理

# ========== 读取差异分析结果 ==========
deg = pd.read_csv("limma_DEP_results.csv", index_col=0)

# ========== ORA分析（Enrichr） ==========
# 获取差异蛋白列表
up_genes = deg[(deg["logFC"] > 1) & (deg["adj.P.Val"] < 0.05)].index.tolist()  # 上调基因列表

# 使用Enrichr在线API进行富集分析
enr = gp.enrichr(
    gene_list=up_genes,  # 基因列表
    gene_sets=["GO_Biological_Process_2023", "KEGG_2021_Human"],  # 基因集数据库
    organism="human",  # 物种
    outdir="enrichr_results"  # 输出目录
)

# 查看结果
print(enr.results.head(10))  # 显示前10个结果

# ========== GSEA分析（预排序） ==========
# 准备排序列表
rnk = deg[["logFC"]].copy()  # 提取log2FC列
rnk.columns = ["score"]  # 重命名
rnk = rnk.dropna().sort_values("score", ascending=False)  # 降序排列

# 运行预排序GSEA
gsea_res = gp.prerank(
    rnk=rnk,  # 排序列表
    gene_sets="GO_Biological_Process_2023",  # 基因集
    min_size=10,  # 最小基因集大小
    max_size=500,  # 最大基因集大小
    permutation_num=1000,  # 置换次数
    outdir="gsea_results",  # 输出目录
    seed=42  # 随机种子
)

# 查看显著结果
sig = gsea_res.res2d[gsea_res.res2d["FDR q-val"] < 0.05]  # 筛选显著结果
print(f"显著富集的通路: {len(sig)}")
print(sig.head(10))

========================================
蛋白质通路富集分析 速查表
========================================

【两种方法选择】
ORA（超几何检验）   → 有明确的差异蛋白列表时用
GSEA（排序分析）    → 不想设阈值，利用全部蛋白信息

【clusterProfiler核心函数】
enrichGO()          → GO过表示分析
enrichKEGG()        → KEGG过表示分析
gseGO()             → GO GSEA分析
gseKEGG()           → KEGG GSEA分析
enrichPathway()     → Reactome通路分析

【可视化函数】
dotplot()           → 气泡图（最常用）
barplot()           → 柱状图
cnetplot()          → 基因-通路网络图
emapplot()          → 通路关联网络图
treeplot()          → 树状聚类图
ridgeplot()         → 山脊图（GSEA）
gseaplot2()         → GSEA曲线图
pathview()          → KEGG通路着色图

【常用基因集数据库】
GO                  → Gene Ontology（BP/CC/MF）
KEGG                → 代谢和信号通路
Reactome            → 人工审核通路
WikiPathways        → 社区维护通路
MSigDB              → 分子标签数据库

【物种代码（KEGG）】
人类                → hsa
小鼠                → mmu
大鼠                → rno

【面试考点】
Q: ORA和GSEA有什么区别？
A: ORA需要差异基因列表（阈值选择影响结果），GSEA用全部基因排序分析

Q: p值和q值的区别？
A: p值是单次检验的显著性，q值是多重检验校正后的p值（FDR）

Q: 富集分析的背景基因集应该怎么选？
A: 应该用实验中检测到的所有蛋白质，而不是全基因组
========================================