【无标题】

Live56789

已于 2025-05-20 11:47:32 修改

阅读量730

点赞数 26

文章标签：数据库

于 2025-05-19 23:21:59 首次发布

本文链接：https://blog.csdn.net/weixin_55795258/article/details/148062695

版权

20250519工作汇总

1.oxldl蛋白质组学数据汇总到AI

1.GO和KEGG分析差异基因富集到的通路里都没有D和H，用其他数据库重新作图。包括reactome，MSigDB ，wikipathways数据库。。明天可以专门运行代码看D和H在哪些信号通路中富集。

2.之前听了一个专门讲论文图表排布的讲座，他以nature的投稿要求讲的，感觉配色在AI里修改是认可的。目前我画图可能重心在图表呈现的内容本身，图的配色，字体等后期都可以在AI调整。

字体类型和大小去AI里调整，代码里调整字号等比较复杂（函数参数多，容易看的眼花缭乱）。

3.把可以用来展示的图的PDF矢量格式都放在一个文件夹内，方便后期取用和编辑。

4.？？？有一个问题是通路富集分析中D或H等基因不在呈现的通路中，展示的是其他基因，这些基因信息是否是富余的？没必要的？有注释通路内的基因的figure看到的比较少。

5.D H
在小鼠中的Entrenz ID
在人类中的Entrenz ID

genes <- c("D", "H")
entrez_ids <- mapIds(org.Hs.eg.db, 
                     keys = genes, 
                     column = "ENTREZID", 
                     keytype = "SYMBOL", 
                     multiVals = "first")

GSEA分析

GSEA 的核心思想是：
根据基因在整个表达谱中的排序（通常是 logFC 或 signal-to-noise ratio）评估某个基因集合是否富集在排序的顶部或底部。
enrichmentScore（ES）和 normalizedEnrichmentScore（NES）是两个核心指标，分别代表富集强度与富集强度的标准化程度。

ES:
ES > 0 → 基因集更集中于排序靠前 → 上调通路
ES < 0 → 基因集更集中于排序靠后 → 下调通路
|ES| 越大 → 富集趋势越强
🟡 但问题是：不同大小的 gene set（例如通路包含基因个数）会影响 ES，无法直接比较不同 gene set。

NES:
表示：标准化后的富集得分，可以用来比较不同 gene set 之间的富集强度。
它是 ES 经过 permutation 后的标准化结果，用于消除 gene set 大小对 ES 的影响。
解读：
NES > 0 → 上调通路
NES < 0 → 下调通路
|NES| 越大 → 趋势越强，更具生物学意义
NES 是 GSEA 推荐用于排序和显著性筛选的标准。

DEG_diff <- subset(DEG, p_value < 0.05 & abs(log2FC) > 0.263)
DEG_sorted <- DEG_diff[order(DEG_diff$log2FC, decreasing = TRUE), ]
DEG_fc=data.frame('Symbol'=DEG_sorted$X,'LogFC'=DEG_sorted$log2FC,stringsAsFactors=F)
DEG_fc[,'ENTREZID'] = mapIds(x = org.Hs.eg.db,
                             keys = DEG_fc$Symbol,
                             keytype = "SYMBOL",
                             column = "ENTREZID") #存在NA
DEG1_fc=na.omit(DEG_fc) #去除NA值

#####################LogFC和ENTREZID组成的数据框
gene_list=DEG1_fc$LogFC                       #提取logFC列
names(gene_list)=DEG1_fc$ENTREZID             #加上ENTREZID
gene_list = sort(gene_list, decreasing = T)  #降序排列

gsea_KEGG <- gseKEGG(gene_list,
                     organism = "hsa",
                     keyType = "kegg",
                     pvalueCutoff = 1)
as.data.frame(gsea_KEGG)
head(gsea_KEGG)

t_index=order(gsea_KEGG$NES,decreasing = T)   #这个位置更换选择的信息列
path=rownames(gsea_KEGG[t_index,])[1:5] #选择展示的pathway
write.csv(gsea_KEGG[t_index, ], file = "gsea_KEGG_result_ES.csv", row.names = TRUE)
write.csv(gsea_KEGG[t_index, ], file = "gsea_KEGG_result_NES.csv", row.names = TRUE)
write.csv(gsea_KEGG[t_index, ], file = "gsea_KEGG_result_pvalue.csv", row.names = TRUE)

install.packages("devtools")
devtools::install_github("kelsiereinaltt/GSEAplot",force = TRUE)
devtools::install_github("junjunlab/GseaVis")
library(GseaVis)
library(enrichplot)
library(ggplot2)

pastel_colors <- c(
  "#FBB4AE",  # light pink
  "#B3CDE3",  # light blue
  "#CCEBC5",  # light green
  "#DECBE4",  # light purple
  "#FED9A6"   # light orange
)
P<-gseaplot2(gsea_KEGG,
          path,
          color = pastel_colors,
          pvalue_table = T, #显示p值
          base_size = 10) #字体大小

cairo_pdf("GSEA_ES_top5.pdf", width = 10, height = 8, bg = "transparent")
print(P)
dev.off()

cairo_pdf("GSEA_NES_top5_2.pdf", width = 10, height = 8, bg = "transparent")
print(P)
dev.off()

cairo_pdf("GSEA_pvalue_top5.pdf", width = 10, height = 8, bg = "transparent")
print(P)
dev.off()

reactome富集分析

Reactome 是一个高度人工注释的通路数据库，涵盖多种物种（人类最全）。

if (!require("ReactomePA")) BiocManager::install("ReactomePA")
library(ReactomePA)
up_ids <- data$ENTREZID[data$group == "up"]
down_ids <- data$ENTREZID[data$group == "down"]
reactome_up <- enrichPathway(gene = up_ids,  # 或 down_entrez
                                 organism = "human", 
                                 pvalueCutoff = 1,   ##阈值放宽些，太严格了通路很少
                                 readable = TRUE)
reactome_down <- enrichPathway(gene = down_ids,  # 或 down_entrez
                             organism = "human", 
                             pvalueCutoff = 1,
                             readable = TRUE)
df_reactome_up <- as.data.frame(reactome_up)
df_reactome_down <- as.data.frame(reactome_down)
write.csv(df_reactome_up, "Reactome_Upregulated.csv", row.names = FALSE)
write.csv(df_reactome_down, "Reactome_Downregulated.csv", row.names = FALSE)
# 可视化
library(enrichplot)
P<-barplot(reactome_up, showCategory = 20, title = "Reactome Pathway Enrichment")
P<-dotplot(reactome_up, showCategory = 20, title = "Reactome Pathway Dotplot")

cairo_pdf("reactome_up_bar.pdf", width = 10, height = 8, bg = "transparent")
print(P)
dev.off()

P<-barplot(reactome_down, showCategory = 20, title = "Reactome Pathway Enrichment")
P<-dotplot(reactome_down, showCategory = 20, title = "Reactome Pathway Dotplot")

cairo_pdf("reactome_down_barplot.pdf", width = 10, height = 8, bg = "transparent")
print(P)
dev.off()

cairo_pdf("reactome_down_dotplot.pdf", width = 10, height = 8, bg = "transparent")
print(P)
dev.off()

在这里插入图片描述
第一张图是上调,第二章第一张图是下调的

MSigDB 数据库

MSigDB 是 GSEA 官方数据库，内容非常丰富，使用 msigdbr 加载通路集。

# 安装并加载
if (!require("msigdbr")) install.packages("msigdbr")
library(msigdbr)
library(clusterProfiler)

# 获取 MSigDB hallmark gene sets（也可以换成 C2/C5/C6）
msig_h <- msigdbr(species = "Homo sapiens", category = "H")  # category 也可设为 "C2" 等

# 提取成符合 enricher 要求的格式（TERM2GENE）
msigdb_term2gene <- msig_h[, c("gs_name", "entrez_gene")]
table(msigdb_term2gene$gs_name)
# up_ids 必须是字符型 ENTREZID
up_ids <- as.character(up_ids)

msig_result <- enricher(gene = up_ids,
                        TERM2GENE = msigdb_term2gene,
                        pvalueCutoff = 1)

# 查看结果
head(as.data.frame(msig_result))

dotplot(msig_result, showCategory = 15, title = "MSigDB Hallmark Pathways")
write.csv(as.data.frame(msig_result), "msigdb_enrichment.csv", row.names = FALSE)

category = "H"还可以替换为其他的值。"H"条件下，上调的差异基因在pvalueCutoff = 1都没有富集到一个通路，改为C2。



参数值
“H” Hallmark gene sets 凝练代表 50 个核心生物过程的“标志性通路”，适合概览变化全貌，最常用 ✅
“C1” Positional gene sets 染色体位置相关基因集（如 chr1q21），较少用
“C2” Curated gene sets 人工整理的通路，包括 KEGG / Reactome / PID / BioCarta / WikiPathways 等合集 ✅
“C3” Motif gene sets 含特定转录因子结合位点或 miRNA 作用位点的基因集
“C4” Computational gene sets (CGP) 来源于癌症细胞系表达谱的聚类分析结果（CGP），表示共表达模块
“C5” Gene Ontology gene sets GO 分类体系中的 BP（生物过程）/ MF（分子功能）/ CC（细胞组分） ✅
“C6” Oncogenic signatures 与特定致癌因素或突变状态相关的表达 signature
“C7” Immunologic signatures 来源于免疫细胞/刺激下的基因表达 signature，用于免疫研究 ✅
“C8” Cell type signature gene sets (v7+) 描述不同细胞类型特征表达的基因集（单细胞数据等）
“C9” Cancer hallmark gene sets (v7.5+) 进一步细化的肿瘤标志性特征通路，适合癌症研究
“C10” Age-related gene sets (v2023+) 新加入的与年龄或衰老相关的表达 signatures

下图是上调的差异基因C2，标题忘记改了，应该是Curated pathways.下调富集的信号通路就非常多了。。在这里插入图片描述

wikipathways数据库的富集分析也是同理

# 安装所需包
if (!require("clusterProfiler")) BiocManager::install("clusterProfiler")
if (!require("rWikiPathways")) BiocManager::install("rWikiPathways")

library(clusterProfiler)
library(rWikiPathways)

# 下载 WikiPathways -> 获取 pathway-gene 对应表
wikipathways_hs <- rWikiPathways::downloadPathwayArchive("Homo sapiens", format = "gmt")

# 读取 GMT 文件为 TERM2GENE
wikipathways_terms <- clusterProfiler::read.gmt(wikipathways_hs)

# 富集分析
wp_result <- enricher(gene = up_entrez,
                      TERM2GENE = wikipathways_terms,
                      pvalueCutoff = 0.05)

# 可视化
dotplot(wp_result, showCategory = 20, title = "WikiPathways Enrichment")
write.csv(as.data.frame(wp_result), "wikiPathways_enrichment.csv", row.names = FALSE)