---

背景

GSEA、GSVA、ssGSEA和GO、KEGG富集分析之间的区别？

从输入的表达矩阵入手
1. GSEA的输入矩阵：差异表达分析后得到的矩阵！
一列基因名/ENTREZ + logFC值（记得按照logFC值从大到小排列）
2. GSVA的输入矩阵：列是样本、行是基因、单元格内是表达量【没有差异分析！】
3. ssGSEA 名字长得像GSEA，其实是GSVA的好兄弟，因为（single sample GSEA）人如其名，都单样本了还怎么做差异分析呢？没有差异分析怎么做GSEA呢？所以ssGSEA就是单样本的GSVA
4. GO/KEGG的输入矩阵：它们连read counts都不要！只用提供基因名！

分析方法
1. GSEA：GSEA同时考虑了基因在整个表达谱中Rank of FoldChange & 同一基因集中的基因在表达谱中的相互之间的距离。 通俗来讲，GSEA基于如下假设： 一个基因集中的基因在表达谱中所处的Rank越极端（高/低FoldChange）& 基因之间的距离越短（Rank相近）= 该基因集越显著。
2. GSVA：更好理解！如果某个基因存在于某个通路，那就给它“一分”，不在就扣它“一分”，这样就能计算得到Enrichment Score（ES）。 这样，某个通路在某个样本中就会有一个最终的得分。所以看GSVA分析完之后的表达矩阵，变成了：列是样本，行是通路，单元格是Enrichment Score（ES）。

3. ssGSEA：只有一个样本，其他计算方法=GSVA。

GSEA分析需要输入gmt文件，这个文件可以来自[GSEA/MsigDB](GSEA | MSigDB (gsea-msigdb.org);
GSVA 和ssGSEA 需要提供geneset文件，他是一种算法，根据您提供的geneset进行相应的分析。
CIBERSORT是一个程序包，使用LM22仅仅对免疫浸润细胞进行分析。

ssGSEA分析

1 文件准备 （文件样式）

1. 表达矩阵
2.cellMarker 原始表格
3. cellMarker 整理成list样式

代码

#加载包
library(tidyverse)
library(data.table)
library(GSVA)

#1.2 准备细胞marker (geneset)
cellMarker <- data.table::fread("cellMarker.csv",data.table = F)
colnames(cellMarker)[2] <- "celltype"
#将cellMarker文件列名的第2个修改为celltype
type <- split(cellMarker,cellMarker$celltype)
#将cellMarker文件以celltype为分组拆分成list数据格式
#处理data.tables列表通常比使用group by参数按组对单个data.table进行操作要慢得多
cellMarker <- lapply(type, function(x){
  dd = x$Metagene
  unique(dd)
})
#将list中每个celltype中的基因进行合并
save(cellMarker,file = "cellMarker_ssGSEA.Rdata")#保存中间文件
load("immune_infiltration//cellMarker_ssGSEA.Rdata")

##1.3 表达量矩阵的准备
###行是基因，列是样本
expr <- data.table::fread("LIHC_fpkm_mRNA_01A.txt",data.table = F)   #读取表达文件
rownames(expr) <- expr[,1]   #将第一列作为行名
expr <- expr[,-1]   #去除第一列
expr <- as.matrix(expr)   #将expr转换为矩阵格式

2. 使用ssGSEA量化免疫浸润

gsva_data <- gsva(expr,cellMarker, method = "ssgsea")

a <- gsva_data %>% t() %>% as.data.frame()
identical(rownames(a),rownames(group))
a$group <- group$group
a <- a %>% rownames_to_column("sample")
write.table(a,"ssGSEA.txt",sep = "\t",row.names = T,col.names = NA,quote = F)

ssGSEA输出结果

3. 可视化

library(ggsci)
library(tidyr)
library(ggpubr)
b <- gather(a,key=ssGSEA,value = Expression,-c(group,sample)) # 提前准备好group信息
ggboxplot(b, x = "ssGSEA", y = "Expression",
         fill = "group", palette = "lancet")+
 stat_compare_means(aes(group = group),
                    method = "wilcox.test",
                    label = "p.signif",
                    symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
                                     symbols = c("***", "**", "*", "ns")))+
 theme(text = element_text(size=10),
       axis.text.x = element_text(angle=45, hjust=1)) 

image.png

CIBERSORT分析

#加载R包
library(e1071)
library(parallel)
library(preprocessCore)
source("CIBERSORT.R")   #加载CIBERSORT程序包
sig_matrix <- "LM22.txt"   
mixture_file = 'LIHC_fpkm_mRNA_01A.txt'   #肿瘤患者表达谱
res_cibersort <- CIBERSORT(sig_matrix, mixture_file, perm=100, QN=TRUE)
save(res_cibersort,file = "res_cibersort.Rdata")   #保存中间文件
load("res_cibersort.Rdata")
res_cibersort <- res_cibersort[,1:22]   
ciber.res <- res_cibersort[,colSums(res_cibersort) > 0]   #去除丰度全为0的细胞

可视化

显示同一个样本，不同免疫细胞的浸润比例，总和是1.

mycol <- ggplot2::alpha(rainbow(ncol(ciber.res)), 0.7) #创建彩虹色板（带70%透明度）
par(bty="o", mgp = c(2.5,0.3,0), mar = c(2.1,4.1,2.1,10.1),tcl=-.25,las = 1,xpd = F)
barplot(as.matrix(t(ciber.res)),
        border = NA, # 柱子无边框写
        names.arg = rep("",nrow(ciber.res)), # 无横坐标样本名
        yaxt = "n", # 先不绘制y轴
        ylab = "Relative percentage", # 修改y轴名称
        col = mycol) # 采用彩虹色板
axis(side = 2, at = c(0,0.2,0.4,0.6,0.8,1), # 补齐y轴添加百分号
 labels = c("0%","20%","40%","60%","80%","100%"))
legend(par("usr")[2]-20, # 这里-20要根据实际出图的图例位置情况调整
       par("usr")[4], 
       legend = colnames(ciber.res), 
       xpd = T,
       fill = mycol,
       cex = 0.6, 
       border = NA, 
       y.intersp = 1,
       x.intersp = 0.2,
       bty = "n")

结果输出：
每一列代表一个样本。可以清晰的看出，不同免疫细胞在肿瘤中的占比。

CIBERSORT免疫浸润

可视化

分组显示--箱型图

a <- read.table("CIBERSORT-Results.txt", sep = "\t",row.names = 1,check.names = F,header = T)
a <- a[,1:22]
identical(rownames(a),rownames(group))
b <- group
class(b$group)
a$group <- b$group
a <- a %>% rownames_to_column("sample")
library(ggsci)
library(tidyr)
library(ggpubr)
#install.packages("ggsci")
#install.packages("tidyr")
#install.packages("ggpubr")

b <- gather(a,key=CIBERSORT,value = Proportion,-c(group,sample))

ggboxplot(b, x = "CIBERSORT", y = "Proportion",
          fill = "group", palette = "lancet")+
  stat_compare_means(aes(group = group),
                     method = "wilcox.test",
                     label = "p.signif",
                     symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
                                      symbols = c("***", "**", "*", "ns")))+
  theme(text = element_text(size=10),
        axis.text.x = element_text(angle=45, hjust=1))

结果输出：

image.png

两种方法略有不同，但都是说明同一个问题，可以比较着看。