背景
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）。
- ssGSEA：只有一个样本，其他计算方法=GSVA。
GSEA分析需要输入gmt文件，这个文件可以来自[GSEA/MsigDB](GSEA | MSigDB (gsea-msigdb.org);
GSVA 和ssGSEA 需要提供geneset文件，他是一种算法，根据您提供的geneset进行相应的分析。
CIBERSORT是一个程序包，使用LM22仅仅对免疫浸润细胞进行分析。
ssGSEA分析
1 文件准备 （文件样式）
代码
#加载包
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)
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))
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")
结果输出：
每一列代表一个样本。可以清晰的看出，不同免疫细胞在肿瘤中的占比。
可视化
分组显示--箱型图
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))
结果输出：
两种方法略有不同，但都是说明同一个问题，可以比较着看。
