非人物种的GSEA&GSVA分析

生信会客厅交流群里有几个朋友问我：MSigDBS数据库只有人的基因集gmt文件可供下载，小鼠的GSEA&GSVA分析如何做？解决问题的根本是做同源基因转换，代码能力强的朋友很容易实现，在线转换的网站也是有的。今天我给大家介绍一个R包——msigdbr，可以更轻松地帮大家得到部分非人物种的gmt文件或基因集列表。

演示数据

演示数据来自GSE109774中的一个小鼠心脏子数据集，整理好的seurat对象文件scRNAmm.RDS，我已上传到百度云。

链接：https://pan.baidu.com/s/1kWrL56O1EkAHPAY3RC4E2g

提取码：29pj

GSEA分析

准备gsea输入文件

##R语言准备gsea输入文件
library(Seurat)
library(tidyverse)
dir.create("GSEA")   
dir.create("GSEA/input")
dir.create("GSEA/output")
scRNA <- readRDS("scRNAmm.RDS")
#提取endothelial cell和fibroblast细胞做富集分析
tmp <- [e-mail]
tmp <- subset(tmp, subset = (tmp$cell_ontology_class=='endothelial cell'|tmp$cell_ontology_class=='fibroblast'))
sub.cells <- rownames(tmp)
scRNAsub <- subset(scRNA, cells=sub.cells)
scRNAsub$cell_ontology_class <- gsub('endothelial cell','endothelial_cell',scRNAsub$cell_ontology_class)
expr <- GetAssayData(scRNAsub, slot = 'data')
expr <- data.frame(NAME=rownames(expr), Description=rep('na', nrow(expr)), expr, stringsAsFactors=F)
write('#1.2', "GSEA/input/expr.gct", ncolumns=1)
write(c(nrow(expr),(ncol(expr)-2)), "GSEA/input/expr.gct", ncolumns=2, append=T, sep='\t')
write.table(expr, "GSEA/input/expr.gct", row.names=F, sep='\t', append=T, quote=F)
line.1 <- c((ncol(expr)-2), 2, 1)
tmp <- table(as.character([e-mail]$cell_ontology_class))
line.2 <- c("#", names(tmp))
line.3 <- c(rep(names(tmp)[1],tmp[1]), rep(names(tmp)[2],tmp[2]))
write(line.1, 'GSEA/input/group.cls', ncolumns=length(line.1), append=T, sep='\t')
write(line.2, 'GSEA/input/group.cls', ncolumns=length(line.2), append=T, sep='\t')
write(line.3, 'GSEA/input/group.cls', ncolumns=length(line.3), append=T, sep='\t')

准备小鼠gmt文件

##安装msigdbr包
install.packages("msigdbr")
#查看msigdbr包支持的物种
msigdbr_species()
   species_name             species_common_name      
 1 Bos taurus               cattle                   
 2 Caenorhabditis elegans   roundworm                
 3 Canis lupus familiaris   dog                      
 4 Danio rerio              zebrafish                
 5 Drosophila melanogaster  fruit fly                
 6 Gallus gallus            chicken                  
 7 Homo sapiens             human                    
 8 Mus musculus             house mouse              
 9 Rattus norvegicus        Norway rat               
10 Saccharomyces cerevisiae baker's or brewer's yeast
11 Sus scrofa               pig   

##准备小鼠gmt文件，以GO数据库的CC数据集为例
library(msigdbr)
#选择自己需要的基因集，category对应MSigDB的主分类，subcategory对应次级分类
genesets = msigdbr(species = "Mus musculus", category = "C5", subcategory = "CC")  
#整理成gmt文件
do_gmt <- select(genesets, c("gs_name","gs_url","gene_symbol")) %>% as.data.frame()
gs_url <- do_gmt[,c("gs_name","gs_url")] %>% unique.data.frame()
do_gmt <- split(do_gmt$gene_symbol, do_gmt$gs_name)
for(i in seq_along(do_gmt)){
    c1 <- names(do_gmt)[i]
    c2 <- gs_url[which(gs_url$gs_name==c1),2]
    c3 <- do_gmt[[i]]
    cb <- c(c1,c2,c3)
    write(cb, 'GSEA/input/genesets.gmt', ncolumns=length(cb), append=T, sep='\t')
}

推荐使用桌面版GSEA软件分析，使用方法参阅《单细胞转录组高级分析五：GSEA与GSVA分析》，部分结果如下：

GSVA分析

使用msigdbr包创建gsea分析使用的gmt文件还有点麻烦，但是用它创建gsva分析使用的基因集列表非常方便。

library(Seurat)
library(GSVA)
library(msigdbr)
library(tidyverse)
dir.create("GSVA")
##选择自己需要的基因集
genesets <- msigdbr(species = "Mus musculus", category = "C8") %>% 
select("gs_name","gene_symbol") %>% as.data.frame()
genesets <- split(genesets$gene_symbol, genesets$gs_name)

##提取表达矩阵
scRNA <- readRDS("scRNAmm.RDS")
expr <- as.matrix(scRNA@assays$RNA@data)
expr <- expr[rowSums(expr)>=1,]  #选取非零基因

##gsva默认开启全部线程计算，建议count用"Poisson"，data用"Gaussian"
es.matrix = gsva(expr, genesets, method="ssgsea", 
                 kcdf="Gaussian", parallel.sz=20) 
write.table(es.matrix, 'GSVA/gsva.xls', row.names=T, col.names=NA, sep='\t')

##差异分析
#挑选差异比较明显的genesets，过程略
write.table(s.gsva, 'GSVA/s_gsva.xls', row.names=T, col.names=NA, sep='\t')

##热图可视化
library(pheatmap)
s.gsva <- read.table('GSVA/s_gsva.xls', header = T, row.names = 1, 
                        sep = '\t', check.names = F)
pheatmap(s.gsva, show_rownames=1, show_colnames=0, fontsize_row=8,
        filename='GSVA/gsva_demo.png', width=15, height=10)

我特意挑选了MSigDB数据库的C8数据集——C8: cell type signature gene sets：Gene sets that contain curated cluster markers for cell types identified in single-cell sequencing studies of human tissue. 感觉它能辅助鉴定细胞类型，从我的结果来看不太靠谱。

交流探讨

如果您阅读此文有所疑惑，或有不同见解，亦或其他问题，欢迎添加我的微信探讨。我工作的公司实力雄厚，拥有10年测序服务经验，欢迎大家联系Kinesin洽谈测序及数据分析业务！