有粉丝提问：跟着我们《生信技能树》的教程：借鉴escape包的一些可视化GSVA或者ssGSEA结果矩阵的方法 做他自己的单细胞数据集的gsva发现内存不够，但他自己的个人电脑已经是32G的顶配了，一时间没办法搞到服务器。
如果这个时候你以为我要发一个服务器的宣传：《96核心384G内存的超级服务器(共享)使用权一年》，那你就错了，我们是认认真真教学！
其实个人电脑已经是32G的顶配了，不应该是存在gsva内存瓶颈，一般来说是自己的代码写的太烂了。
安装seurat-data包获取测试数据
代码其实一句话即可：
devtools::install_github('satijalab/seurat-data')
但是因为是在GitHub，所以中国大陆地区的小伙伴有时候会遇到：
> devtools::install_github('satijalab/seurat-data')
Error: Failed to install 'SeuratData' from GitHub:
Timeout was reached: [api.github.com] Operation timed out after 10002 milliseconds with 0 out of 0 bytes received
正常情况下应该是：
* installing *source* package ‘SeuratData’ ...
** using staged installation
** R
** exec
** inst
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
** building package indices
** testing if installed package can be loaded from temporary location
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (SeuratData)
如果你连GitHub都困难，你可以翻一下我们《生信技能树》很久以前的教程，其实把这个GitHub仓库自己备份到自己的gitee仓库就好了。
查看以及认识测试数据
library(SeuratData)
AvailableData()
# 假如你没有下载pbmc3k数据集 ，就使用下面的 代码
# InstallData("pbmc3k") # (89.4 MB)
# trying URL 'http://seurat.nygenome.org/src/contrib/pbmc3k.SeuratData_3.1.4.tar.gz'
# 上面的 InstallData 代码只需要运行一次即可哈
data("pbmc3k")
exp <- pbmc3k@assays$RNA@data
dim(exp)
exp[1:4,1:4]
table(is.na(pbmc3k$seurat_annotations ))
as.data.frame(table(pbmc3k$seurat_annotations))
这个测试数据pbmc3k，是已经做好了降维聚类分群和注释啦，这个数据集超级出名的，大家可以自行去了解它的背景知识哈。如下所示：
> exp[1:4,1:4]
4 x 4 sparse Matrix of class "dgCMatrix"
AAACATACAACCAC AAACATTGAGCTAC AAACATTGATCAGC AAACCGTGCTTCCG
AL627309.1 . . . .
AP006222.2 . . . .
RP11-206L10.2 . . . .
RP11-206L10.9 . . . .
> table(is.na(pbmc3k$seurat_annotations ))
FALSE TRUE
2638 62
> as.data.frame(table(pbmc3k$seurat_annotations))
Var1 Freq
1 Naive CD4 T 697
2 Memory CD4 T 483
3 CD14+ Mono 480
4 B 344
5 CD8 T 271
6 FCGR3A+ Mono 162
7 NK 155
8 DC 32
9 Platelet 14
如果是follow我这个教程，你完全不需要使用这个测试数据pbmc3k，只要是一个seurat的对象即可。
如果跑gsva面临内存问题
代码如下：
sce=pbmc3k
mat=sce@assays$RNA@counts
mat[1:4,1:4]
X=as.matrix(mat)
dim(X)
gmtfile ='../MSigDB/symbols/h.all.v7.2.symbols.gmt'
library(GSEABase)
library(GSVA)
geneset <- getGmt( gmtfile )
geneset
# method=c("gsva", "ssgsea", "zscore", "plage"),
# kcdf=c("Gaussian", "Poisson", "none"),
# mx.diff=FALSE: ES is calculated as the maximum distance of the random walk from 0.
# mx.diff=TRUE (default): ES is calculated as the magnitude difference between the largest positive and negative random walk deviations.
es.max <- gsva(X, geneset,
mx.diff=FALSE, verbose=FALSE,
parallel.sz= 8 )
其中 h.all.v7.2.symbols.gmt 文件需要自行下载哦，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: 免疫相关基因集合。
因为这个案例，仅仅是3千个单细胞，所以很难遇到内存不足的报错，其实真实情况下 往往是5万个单细胞以上，大家如果有自己的数据，就可以测试了。
其实解决内存不足问题超级简单
只需要把前面的单细胞的表达矩阵，根据细胞类型拆分开来，分开独立跑gsva，然后再合并，代码如下：
table(sce@meta.data$orig.ident)
es.max.list <- lapply(unique(sce@meta.data$orig.ident),function(i){
Y=X[,sce@meta.data$orig.ident == i ]
es.max <- gsva(Y, geneset,
mx.diff=FALSE, verbose=FALSE,
parallel.sz=8)
return(es.max)
})
es.max=do.call(cbind,es.max.list)
es.max[1:4,1:4]
当然了，这个时候表达矩阵的顺序其实是发生了变化哦，后续如果对这个gsva矩阵绘制热图的时候，需要调整回来哦。
es.max[1:4,1:4]
X[1:4,1:4]
identical(colnames(es.max),colnames(X))
pos=match(colnames(X),colnames(es.max))
es.max=es.max[,pos]
identical(colnames(es.max),colnames(X))
其实呢，也不一定要根据细胞类型来拆分表达量矩阵分开跑gsva，只需要拆分即可。毕竟啊， 按照细胞类型来拆分表达量矩阵，通常呢，并不是等价的。
as.data.frame(table(pbmc3k$seurat_annotations))
Var1 Freq
1 Naive CD4 T 697
2 Memory CD4 T 483
3 CD14+ Mono 480
4 B 344
5 CD8 T 271
6 FCGR3A+ Mono 162
7 NK 155
8 DC 32
9 Platelet 14
最后的可视化如下
ac=data.frame(celltype=pbmc3k$seurat_annotations)
rownames(ac)=rownames(sce@meta.data)
pheatmap::pheatmap(es.max,
annotation_col = ac,
show_colnames = F)
挺漂亮的：
如果是差异分析，基本上看我五年前的《数据挖掘》系列推文 就足够了；