基因富集分析（SPIA）

2025-02-22

生信

3.1k 词

引言
SPIA
实践

引言

一般来说，使用GO和KEGG数据库对差异基因进行基因富集分析，多是

ORA：基于对特定组内外（如通路）的差异基因、背景基因数量进行检验，看组内外的差异基因是否随机分布。比如ClusterProfiler包中的enrichGO()和enrichKEGG()
GSEA：差异基因按log2FoldChange降序排序后，顺序遍历，按差异基因是否在特定组（如通路）内进行加减分，并进行检验，最后NES可以看出组内的整体表达趋势。比如同样是ClusterProfiler包中的gseGO()和gseKEGG()

GO数据库条目之间至少可能还有Positively regulates和Negatively regulates两种关系，暂且不论。
但KEGG一个通路内的基因，极度简化，既有正向的也有负向的，而且还有更复杂的，有一些使用KEGG数据库的做法是把上下调的差异基因分别进行ORA，或者一起进行GSEA，仔细想想其实很没道理。以GSEA举例，如果一个通路内正向的基因上调，负向的基因下调，这样子整体表达趋势却可能是不变，会认为不显著。

（GSEA针对这个问题，有一种办法，可以取log2FC的绝对值进行单端富集）

SPIA

于是介绍一个也很古老的方法，SPIA，基于拓扑结构对KEGG通路进行评分。
简单来说，SPIA有2个指标

NDE：这个和上面的ORA差不多
PERT：通过激活/抑制关系，模拟基因表达变化对通路活性的影响

最后得到累计扰动tA，正负可以认为是激活或抑制，以及全局p值pG。

所以当然，因为NDE这部分的存在，SPIA的的结果没那么可信。
并且PERT部分的质量严重依赖于通路的质量，甚至一个通路的新旧数据能得出完全相反的结果。

~~这也就说明有很好的捏造能力~~

实践

R可以用SPIA包。

面对那么多通路，计算量有点大，选择先进行ORA，然后拿富集出的感兴趣的通路，比如Cellular Processes

BiocManager::install("SPIA")
# kegg ora
gene.df <- bitr(rownames(DEGs),fromType="SYMBOL",toType=c("ENTREZID"),OrgDb = org.Hs.eg.db,drop = F)
gene.df <- gene.df[!duplicated(gene.df$SYMBOL),]
DEGs$ENTREZID <- gene.df$ENTREZID
DEGs <- subset(DEGs,subset= !is.na(ENTREZID))
kegg_result <- enrichKEGG(unique(DEGs_filtered$ENTREZID))
# 筛选出Cellular Processes
kegg_result_cellular_processes <- kegg_result
kegg_result_cellular_processes@result <- subset(kegg_result@result,subset=category=="Cellular Processes")
# 加载相关通路信息
spia_kegg_dir <- system.file("extdata/keggxml/hsa",package="SPIA")
for (path_id in kegg_result_cellular_processes@result$ID) {
  download.file(url = paste0("http://rest.kegg.jp/get/", path_id, "/kgml"), 
                  destfile = file.path(spia_kegg_dir, paste0(path_id, ".xml")),
                  quiet = FALSE,
                  mode = "wb")
}
makeSPIAdata(kgml.path=spia_kegg_dir,organism="hsa",out.path="./")

或者筛选出所有Cellular Processes通路

pathway_ids <- gsub("path:", "", names(keggList("pathway", "hsa")))

cellular_processes <- sapply(pathway_ids, function(pid) {
  tryCatch({
    info <- keggGet(pid)[[1]]
    if (any(grepl("Cellular Processes", info))) {
      return(pid)
    }
    return(NA)
    }, error = function(e) {
      print(paste0("err on ", pid))
      return(NA)
  })
})

spia_kegg_dir <- system.file("extdata/keggxml/hsa",package="SPIA")
for (path_id in pathway_ids[!is.na(cellular_processes)]) {
  download.file(url = paste0("http://rest.kegg.jp/get/", path_id, "/kgml"), 
                  destfile = file.path(spia_kegg_dir, paste0(path_id, ".xml")),
                  quiet = FALSE,
                  mode = "wb")
}
makeSPIAdata(kgml.path=spia_kegg_dir,organism="hsa",out.path="./")

最后进行SPIA

# 数据准备
gene_list <- DEGs_filtered$avg_log2FC
names(gene_list) <- DEGs_filtered$ENTREZID
all_gene <- bitr(rownames(seu_obj),fromType="SYMBOL",toType=c("ENTREZID"),OrgDb = org.Hs.eg.db,drop = F)
# SPIA
spia_result <- spia(de = gene_list,all = all_gene$ENTREZID,organism = "hsa",data.dir = "./")
# 可视化
plotP(spia_result)
ggplot(spia_result, aes(x = tA, y = reorder(Name, -tA),fill = pG))+geom_bar(stat = "identity")

（有一些通路不会被SPIA计算）

分类

标签

基因富集分析（SPIA）

引言

SPIA

实践