查看原文
其他

差异基因没办法富集到通路你就放弃了吗

生信技能树 生信菜鸟团 2022-06-07

我在生信技能树分享了一个教程:不要怀疑,你的基因就是没办法富集到统计学显著的通路,强调了大家做生物信息学数据分析的同时,一定要加强统计学基础,比如把差异基因集(500个左右的基因)富集到KEGG数据库通路,本质上就是对每个通路单独做一次超几何分布检验罢了!

如果你的差异基因集(500个左右的基因)本身确实没有一些生物学功能的偏好性,非常有可能是在统计学显著阈值条件限定下,就是拿不到富集的kegg通路。这个时候,大家可以修改代码,比如:

kk.down <- enrichKEGG(gene         =  gene_down,
                        organism     = 'hsa'
                        pvalueCutoff = 0.9,
                        qvalueCutoff =0.9)
# 需要自己差异分析筛选得到 gene_down 基因集 ,然后进行超几何分布检验

拿到通路后自己去调整阈值。

或者可以去MSigDB(Molecular Signatures Database)看看其它功能基因集是否可以富集到结果,MSigDB数据库中定义了已知的基因集合:http://software.broadinstitute.org/gsea/msigdb 包括H和C1-C7八个系列(Collection),每个系列分别是:

  • H: hallmark gene sets (癌症)特征基因集合,共50组,最常用;
  • C1: positional gene sets 位置基因集合,根据染色体位置,共326个,用的很少;
  • C2: curated gene sets:(专家)校验基因集合,基于通路、文献等:
  • C3: motif gene sets:模式基因集合,主要包括microRNA和转录因子靶基因两部分
  • C4: computational gene sets:计算基因集合,通过挖掘癌症相关芯片数据定义的基因集合;
  • C5: GO gene sets:Gene Ontology 基因本体论,包括BP(生物学过程biological process,细胞原件cellular component和分子功能molecular function三部分)
  • C6: oncogenic signatures:癌症特征基因集合,大部分来源于NCBI GEO  发表芯片数据
  • C7: immunologic signatures: 免疫相关基因集合。

但是这个分析的前提是拿到差异基因集(500个左右的基因)

有些时候,大家挑选到的差异基因集数量过少或者过多,又不知道如何调整阈值,得到的结果会稀奇古怪,不可理喻。也就是说,这个分析严重依赖于差异基因集的挑选,因人而异,这样不好。

所以我们会推荐大家走GSEA分析,它并不会依赖于差异分析本身,我在生信技能树多次讲解GSEA分析:

在R里面的实现方式也是非常的简单:

# 需要自己对全部的基因进行一个排序,排序后的变量是geneList
library(clusterProfiler)
    kk_gse <- gseKEGG(geneList     = geneList,
                      organism     = 'hsa',
                      nPerm        = 1000,
                      minGSSize    = 10,
                      pvalueCutoff = 0.9,
                      verbose      = FALSE)
    tmp=kk_gse@result
    gsea_kk=DOSE::setReadable(kk_gse, OrgDb='org.Hs.eg.db',keyType='ENTREZID')

也可以对MSigDB数据库的全部基因集进行批量gsea分析,然后细致的去挑选结果。

# 选择gmt文件(MigDB中的全部基因集)
# 自己在 MigDB 官网下载 gmt 文件哦!
  d='~/biosoft/MSigDB/symbols/'
  gmts <- list.files(d,pattern = 'all')
  gmts
  #GSEA分析
  library(GSEABase) # BiocManager::install('GSEABase')
  ## 下面使用lapply循环读取每个gmt文件,并且进行GSEA分析
  ## 如果存在之前分析后保存的结果文件,就不需要重复进行GSEA分析。
  f='gsea_results.Rdata'
  if(!file.exists(f)){
    gsea_results <- lapply(gmts, function(gmtfile){
      # gmtfile=gmts[2]
      geneset <- read.gmt(file.path(d,gmtfile)) 
      print(paste0('Now process the ',gmtfile))
      egmt <- GSEA(geneList, TERM2GENE=geneset, verbose=FALSE)
      head(egmt)
      # gseaplot(egmt, geneSetID = rownames(egmt[1,]))
      
      return(egmt)
    })
    # 上面的代码耗时,所以保存结果到本地文件
    save(gsea_results,file = f)
  }
  load(file = f)
  #提取gsea结果,熟悉这个对象
  gsea_results_list<- lapply(gsea_results, function(x){
    cat(paste(dim(x@result)),'\n')
    x@result
  })
  dev.off()
  gsea_results_df <- do.call(rbind, gsea_results_list)
   
  write.csv(gsea_results_df,file = 'gsea_results_df.csv')
  
  gseaplot(gsea_results[[2]],'KEGG_CELL_CYCLE'
  gseaplot(gsea_results[[2]],'FARMER_BREAST_CANCER_CLUSTER_6'
  gseaplot(gsea_results[[5]],'GO_CONDENSED_CHROMOSOME_OUTER_KINETOCHORE'
 

如果你差异基因没办法富集到通路,使用GSEA也许会有不一样的结果哦!

比如下面的这个通路,如果你使用上下调基因集去注释,通常呢你看不到它被富集到,因为top上下调基因都不是这个通路里面的,但是gsea分析就可以看到,它其实是显著下调的!


您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存