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

Posted on by

我在生信技能树分享了一个教程：不要怀疑，你的基因就是没办法富集到统计学显著的通路，强调了大家做生物信息学数据分析的同时，一定要加强统计学基础，比如把差异基因集（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也许会有不一样的结果哦！

未完，图表在生信菜鸟团推文：https://mp.weixin.qq.com/s/wNpjogN2edrtN-t2CbWvVA

​

Comments are closed.