GSVA或者GSEA各种算法都是可以自定义基因集的

Posted on by

表达矩阵的标准分析通常是不够的，定位到成百上千个有统计学显著变化的差异表达基因后，同样是可以有成百上千个生物学功能注释（最出名的是GO功能和KEGG通路），普通的超几何分布检验已经不能满足大家多元化的分析了，所以就有了大家耳熟能详的GSEA分析，以及绝大部分人比较陌生的GSVA分析。

GSVA分析的文章发表于2013年，GSVA: gene set variation analysis for microarray and RNA-Seq data 同样是broad 研究生出品，其在2005年PNAS发表的gsea已经高达1.4万的引用了，不过这个GSVA才不到300。去年我就介绍过一波它的分析流程，在：使用GSVA方法计算某基因集在各个样本的表现 非常简单的代码，所以各个培训机构，公司人员都开始学习和二次创作进而分享。

成百上千个生物学功能注释（最出名的是GO功能和KEGG通路）通常是从数据库里面下载，就是基因集，也可以叫pathway，如下：

pathway的具体含义

pathway在我这里是基因集的别名，其中msigdb有着丰富的基因集，MSigDB（Molecular Signatures Database）数据库中定义了已知的基因集合：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: 免疫相关基因集合。

但是是可以自定义基因集

之所以大家不知道可以自定义基因集，其实是因为，大家做数据分析的时候，习惯了软件包作者打包或者说封装好的函数，如下：

library(clusterProfiler)
data(gcSample)
enrichKK <- enrichKEGG(gcSample[[1]])
enrichKK
class(enrichKK)

简单的几句话就对你感兴趣的基因集，进行了KEGG数据库的普通的超几何分布检验，是不是很神奇，你压根就没有去下载KEGG数据库，也不知道里面有多少条通路，每个通路具体多少个基因，不同通路之间基因的重合情况。

实际上我给大家写的GEO数据库挖掘的标准代码教程里面是有GSEA分析的：

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

### 对 MigDB中的全部基因集 做GSEA分析。
# http://www.bio-info-trainee.com/2105.html
# http://www.bio-info-trainee.com/2102.html 
{
 load(file = 'anno_DEG.Rdata')
 geneList=DEG$logFC
 names(geneList)=DEG$symbol
 geneList=sort(geneList,decreasing = T)
 #选择gmt文件（MigDB中的全部基因集）
 d='../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
 })
 gsea_results_df <- do.call(rbind, gsea_results_list)
 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') 

}

看起来是稍微有点复杂了，其实如果你但凡是听过我的R语言课程，就很容易理解。核心就是2句话，读gmt格式的文件，通常是在，MSigDB（Molecular Signatures Database）数据库里面下载的文件，如下

image-20191214200909764

然后做GSEA或者GSVA分析罢了，如果是GSEA，就：

 library(GSEABase) # BiocManager::install('GSEABase')
 geneset <- read.gmt(file.path(d,gmtfile)) 
egmt <- GSEA(geneList, TERM2GENE=geneset, verbose=FALSE)
head(egmt)

如果是GSVA，就：

 library(GSVA) # BiocManager::install('GSVA')
 geneset <- getGmt(file.path(d,gmtfile)) 
 es.max <- gsva(X, geneset, 
 mx.diff=FALSE, verbose=FALSE, 
 parallel.sz=1)

可以看到，它们读取gmt格式文件的函数不一样，不过问题不大，只需要是标准的gmt格式文件即可，那么什么是标准的gmt格式文件呢，简单谷歌即可：

image-20191214201146450

其实就是，每个基因集是一行，然后第一列是基因集的ID，第二列是名字，后面的列就是基因集所含有的基因。

而不同的基因集，在不同的行，可以有不同数量的基因啦。

所以你只需要自己制作这样的gmt文件，就可以啦，使用上面我们提到的函数进行读取。

如果你自己手动制作gmt文件困难肿么办

几年前我写过一个函数，专门把基因集，写出成为gmt文件，你可以试试看，能不能达到我几年前的水平哈！

​

Comments are closed.