R语言计算肿瘤干细胞指数,容易上手的代码来了
在之前发表过的一篇文章中提到了研究肿瘤干细胞性的重要性以及相关的计算方法。2020年了,机械化自动化的时代了,敲敲代码,就可以实现对几十甚至上万的肿瘤样本“肿瘤干细胞指数“的刻画了,本篇总结主要想和大家交流如何通过ssgsea(single-samplegene set enrichment analysis)算法实现对肿瘤样本肿瘤干细胞指数的计算。
一、ssgsea(single-samplegene set enrichment analysis)算法的介绍该算法属于一种特殊的GSEA(Gene Set Enrichment Analysis),它主要针对单样本无法做GSEA而提出的一种实现方法,原理上与GSEA是类似,也是根据表达谱文件计算每个基因的rank值,再进行后续的统计分析,但ssgsea准备的不是gct格式的表达谱文件,最终得到的结果是每个样本在对应的背景基因集下自己有一个分值。如Fig 1所示从左往后的变化过程。
Fig 1. ssgsea算法思想
二、R运行ssgsea准备材料
(1)安装R包“GSVA”
source("http://bioconductor.org/biocLite.R")
biocLite("GSVA")
(2)肿瘤样本表达谱数据(行是基因,列是样本)
(3)肿瘤干细胞基因构成的基因集合(gene set):这点我已经给大家提前查过了,目前广泛被大家认可的肿瘤干细胞基因集合是来自“Cancer stemness, intratumoral heterogeneity, and immune responseacross cancer"整理的109个肿瘤干细胞基因集合。这个基因集是在该文章的补充材料中,具体下载网址如括号中展示的那样(https://www.pnas.org/content/pnas/suppl/2019/04/17/1818210116.DCSupplemental/pnas.1818210116.sd01.xlsx),如Fig 2所示,下述红色框框起来的那一列就是整合后的肿瘤干细胞基因。
Fig 2. 109个肿瘤干细胞基因
三、代码实操(小小尝试)
#加载相应的包和数据
rm(list=ls())##清空之前的所有变量
#加载相应的包
gene.exp=matrix(sample(1:1000,10000,replace=T),200,50)##表达谱数据的维度:200个基因,50个样本
dim(gene.exp)##可视化表达谱维度
rownames(gene.exp)=paste("g",1:200,sep="")##給行基因命名(g1,g2...)
colnames(gene.exp)=paste("s",1:50,sep="")##給列样本命名(s1,s2...)
head(gene.exp)
stem.gene.set=list(c("g1","g3","g5","g8","g9","g15"))##定义肿瘤干细胞基因
##运行ssgsea算法
gsva_score<- gsva(gene.exp, stem.gene.set, method="ssgsea",verbose=FALSE, parallel.sz=1)
class(gsva_score)
dim(gsva_score)
gsva_score
最终运行的结果如Fig 3所示:
Fig 3. ssgsea运行结果
以上,就是通过ssgsea算法得到样本肿瘤干细胞指数的代码。自己待计算样本的表达谱数据和背景基因集大家可以根据自己的需求进行替代。