免疫浸润利器——ssGSEA纯代码实操（附可视化操作）

我们之前讲过计算免疫微环境的两种主流方法为CIBERSORT与ssGSEA，那么今天我们来学习一下ssGSEA为何方神圣吧。

1.ssGSEA简介

single sample Gene Set Enrichment Analysis(ssGSEA)?

虽然不认识他，但是他的兄弟GSEA总还是熟悉的吧。与他兄弟不同的是，它主要针对单样本无法做GSEA而提出的一种实现方法。那么他又是如何估计每个样本中的免疫细胞含量的呢？此处也用到了一种算法——随机漫步（此处略过，毕竟咱们也不是研究算法的嘛，这些工作就交给科研工作者就好）。这种算法精简版就是可以计算出每个基因集在这个样本中的富集分数。

大家之前应该看过CIBERSORT的代码操作了，那么可能各位小伙伴也会去选择CIBERSORT。CIBERSORT固然有他的优势：R包简便，也有网页在线操作优势显而易见。但是各位看官莫要忘记，“LM22”文件中只有22种免疫细胞，要是再想多一些那就只能sorry了。

ssGSEA其中最明显的特征就是可以高度定制化。根据输入的gmt文件即可对样本数据进行随机游走，之后对每个样本计算出其免疫细胞的富集得分。

话不多说，马上进入实操环节。##输入文件准备首先我们需要准备输入文件。其中文件主要为3个：

文件	说明
ssGSEA.R	本次操作所演示的代码，同学可在文末资料处领取
Immunecell.gmt	文件中对应了每个免疫细胞中的基因
symbol.txt	需要运算的样本文件（因文件过大，可在TCGA官网自行下载）

“symbol.txt”（仅展示部分）:数据为TCGA-LUSC的转录本FPKM文件

“immune cell.txt”（仅展示部分）:免疫细胞基因集

2.ssGSEA分析

大家可以点开Rstudio,按照如下方式导入资料脚本” ssGSEA.R”

接下来大家将光标放在第一行后，点击“Run”即可运行每行代码，此处不建议刷住多行代码直接点击Run,对于初学者应当感悟每行代码所代表的意义，而且万一出错会导致二次检查，费时费力

rm(list=ls())setwd("E:\rxys\03")          #设置工作目录#引用分析所用包。此处需要提前下载欧library(GSVA)library(limma)library(GSEABase)

#读取输入文件，并对输入文件处理symbol=read.table("symbol.txt",sep="  ",header=T,check.names=F)#读取输入文件，此处为TCGA-LUSC的FPKM标准化数据symbol=as.matrix(symbol)rownames(symbol)=symbol[,1]exp_data=symbol[,2:ncol(symbol)]dimnames=list(rownames(exp_data),colnames(exp_data))exp_data1=matrix(as.numeric(as.matrix(exp_data)),nrow=nrow(exp_data),dimnames=dimnames)exp_data1=avereps(exp_data1)exp_data1=exp_data1[rowMeans(exp_data1)>0,]geneSet=getGmt("immune cell.gmt", geneIdType=SymbolIdentifier())#读取我们所提供的免疫细胞基因文件ssGSEA_Score=gsva(exp_data1, geneSet, method='ssgsea', kcdf='Gaussian', abs.ranking=TRUE)#ssGSEA计算normalize=function(x){return((x-min(x))/(max(x)-min(x)))}#定义ssGSEA_Score矫正函数norm_ssGSEA_Score=normalize(ssGSEA_Score)#对ssGSEA_Score进行矫正norm_ssGSEA_Score=rbind(id=colnames(norm_ssGSEA_Score),norm_ssGSEA_Score)write.table(norm_ssGSEA_Score,file="norm_ssGSEA_Score.txt",sep="  ",quote=F,col.names=F)#此处的输出文件即为ssGSEA富集得分文件

此时会在文件夹中出现结果文件“norm_ssGSEA_Score.txt”，我们将其打开看一下内容。此时在文件中会对应展示每个样品中免疫细胞的富集分数，为了更好的展示，我们使用excel打开该文件（如下图：仅显示部分截图）

3.可视化操作

我们在论文中往往会对数据进行可视化操作，接下来就让我们对数据进行一下绘图操作。在此之前要对于结果文件进行处理，因为我们在ssGSEA中大多数情况下只对于肿瘤样本进行计算，因此我们要先筛掉那些正常样本

screen=sapply(strsplit(colnames(norm_ssGSEA_Score),"\-"),"[",4)screen=sapply(strsplit(screen,""),"[",1)screen=gsub("2","1",screen)data=norm_ssGSEA_Score[,screen==0]#对肿瘤样品分组，TCGA的命名规则可以复习一波

随后我们对样本进行聚类分组，并进行可视化热图绘制

h_cluster = hclust(dist(t(norm_ssGSEA_Score)))#计算样品距离并聚类cluster_result=cutree(h_cluster,3)#将聚类结果分成低，中，高免疫三组library(pheatmap)norm_ssGSEA_Score1=read.table("norm_ssGSEA_Score.txt",sep="  ",header=T,row.names=1,check.names=F)#读取结果文件V1<-names(cluster_result)V2<-as.integer(cluster_result)category<-data.frame(V1,V2)category$V1<-as.character(category$V1)category[,2]=paste0("Cluster",category[,2])category=category[order(category[,2]),]norm_ssGSEA_Score1=norm_ssGSEA_Score1[,as.vector(category[,1])]cluster=as.data.frame(category[,2])row.names(cluster)=category[,1]colnames(cluster)="Cluster"pdf("heatmap.pdf",height=5,width=9)color.key<-c("#3300CC","#3399FF","white","#FF3333","#CC0000")pheatmap(norm_ssGSEA_Score1, annotation=cluster,        color =colorRampPalette(color.key)(50),        cluster_cols =F,        fontsize=8,        fontsize_row=8,        scale="row",        show_colnames=T,        fontsize_col=3,        main = "Example heatmap")dev.off()

从图中可以看到低、中、高三组有较为明显的区分，接下来我们就可以利用这些分组完成一些自己课题的后续研究了。

4.实际应用

想必到了这里，大家应该对ssGSEA有了一个更为深刻的了解，趁热打铁，我们直接来到应用部分。我们今天介绍的是2019年发表在Oncology reports上的一篇文章“Tumor‑infiltrating M2 macrophages driven by specific genomic alterations are associated with prognosis in bladder cancer”

这篇文章采用了ssGSEA与CIBERSORT两种方法进行免疫浸润的相关研究。最终通过膀胱癌的组织学分级、病理分期与M2巨噬细胞浸润水平的关联性，得出膀胱癌特异性基因组改变可能是M2巨噬细胞浸润的重要驱动力在这里我们放两张论文中的图片：

是不是与我们的可视化结果较为相似，不过这个是利用CIBERSORT的结果做出来的。不过我们可以看一下下面这幅图

在这幅图中作者采用ssGSEA的方法对样品中的富集分数的计算（也就是我们的结果文件啦）再配合临床分期进行可视化操作。从而得出与低度膀胱癌组织相比，B细胞，巨噬细胞，嗜中性粒细胞等细胞在高度膀胱癌中富集更为显著。文献的全文我们已经打包在资料中，感兴趣的小伙伴可以复现一下这篇文献。

结语

无论是在线网站分析或者R语言代码分析，我们都不要忘记引用算法文献呀！！！

[1] GSVA: The Gene Set Variation Analysis package for microarray and RNA-seq data

[2] Xue Y, Tong L, LiuAnwei Liu F, et al. Tumor‑infiltrating M2 macrophages driven by specific genomic alterations are associated with prognosis in bladder cancer. Oncol Rep. 2019;42(2):581-594. doi:10.3892/or.2019.7196

[3] Jia Q, Wu W, Wang Y, et al. Local mutational diversity drives intratumoral immune heterogeneity in non-small cell lung cancer. Nat Commun. 2018;9(1):5361. Published 2018 Dec 18. doi:10.1038/s41467-018-07767-w

最后给小伙伴们留一个小作业：大家看到的这幅图可以依据我们的结果文件并利用pheatmap包完成其中热图的部分绘制，大家可以动手尝试一下了，看一看谁做的比较像吧:

回复"日行一膳03"得范例文献和数据~

欢迎大家关注解螺旋生信频道-挑圈联靠公号~

本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。