算法原理不理解可以,但是请清楚一下概念

《道德经》“玄之又玄,众妙之门”

gsea和gsva算法大家应该是都很熟悉了,我也多次讲解:

不过里面有一个算法表格很难理解:

GSEA算法表格

所以最近又有学徒咨询这个细节了,她的疑问是,为什么ssGSEA得分,跟其基因集里面的每个基因的表达量,并不是完全相关的,我这次写了一个例子来讲解。

模拟表达矩阵后,获得基因排序表格

比如我们模拟一下 26个基因的4个样本的表达量矩阵,然后查看自己感兴趣的3个基因在4个样本的顺序。

dat=data.frame(s1=rnorm(26),
 s2=rnorm(26),
 s3=rnorm(26),
 s4=rnorm(26))
rownames(dat)=LETTERS
dat['A',]
pos=match(c('D','G','Z'),rownames(dat))
apply(dat, 2, function(x){
 order(x,decreasing = T)[pos]
})

这里面的rnorm函数是随机的,所以你复制粘贴运行我这个代码,结果会有点不一样。

然后对它进行ssGSEA算法计算,如下:

library(GSVA) 
library(GSEABase) 
gmtFile='test.gmt'
tmp=paste(c('test','test','D','G','Z'),sep = '\t')
write.table(t(as.data.frame(tmp)),
 file = gmtFile,sep = '\t',quote = F,row.names = F,col.names = F)
geneSet=getGmt(gmtFile, 
 geneIdType=SymbolIdentifier())
geneSet

ssgseaScore=gsva(as.matrix(dat), 
 geneSet, method='ssgsea', kcdf='Gaussian', abs.ranking=TRUE)
ssgseaScore

结果如下:

> ssgseaScore
 s1 s2 s3 s4
test 0.03809468 0.2781415 -0.7218585 -0.6190109
> apply(dat, 2, function(x){
+ sort(order(x,decreasing = T)[pos])
+ })
 s1 s2 s3 s4
[1,] 17 2 18 13
[2,] 23 23 20 17
[3,] 26 25 25 21

可以看到总共是26个基因的表达量,我们感兴趣的是3个基因,使用ssGSEA后:

  • s1样本的ssGSEA打分是0.04,基因都排在后面
  • s2样本的ssGSEA打分是0.27,因为有一个基因排的非常靠前!
  • s3样本的ssGSEA打分是 -0.72, 因为3个基因都排在最后面
  • s4样本的ssGSEA打分是 -0.62,同样的,3个基因都排在后面,但是略高于s3样本。

总之,如果要ssGSEA打分高,你的基因需要尽可能表达量排序靠前,反之,需要靠后使得ssGSEA打分降低,甚至为负值。

也就是说,仅仅是关心基因的表达量排序,而不是表达量本身。

真实数据

取一个项目的表达矩阵,然后提取11个基因组成的基因集的表达量热图,加上这11个基因的ssGSEA在这个表达量矩阵计算结果。

热图可视化如下:

n=t(scale(t( mat[treg,] ))) 
n[n>2]=2 
n[n< -2]= -2
n[1:4,1:4] 
ac=data.frame(gsea=as.numeric(ssgseaScore))
rownames(ac)=colnames(n)
pheatmap(n,annotation_col = ac,
 show_colnames =F,show_rownames = T)

image-20200508173043392

可以看到,ssGSEA得分,跟基因表达量,还是有关系的,毕竟表达量高基因排名也高啊!

你大概率上无需把最开始那个算法PPT里面的公式写出来

首先FC在单个样本里面是不存在的,但是可以把表达量进行zscore,代入公式,看看能不能计算得到差不多的ssGSEA值。

欢迎邮件交流你的代码,我的邮箱 jmzeng1314@163.com

文末友情宣传

强烈建议你推荐我们生信技能树给身边的博士后以及年轻生物学PI,帮助他们多一点数据认知,让科研更上一个台阶:

Comments are closed.