ssGSEA是一种常用于免疫细胞浸润分析的方法。该方法通过将每个样本的基因表达数据与特定的基因集（免疫细胞基因集）进行比较，来估计该基因集在该样本中的相对富集程度。
2. 使用
反/去卷积
★ 单样本基因集富集分析（single sample gene set enrichment analysis, ssGSEA），是GSEA方法的扩展，计算每个样本和基因集配对的富集分数。首先，将所有基因按照其表达量从大到小进行排序，并计算在某个基因集内，基因表达量较高的基因的累积分布函数。这个累积分布函数被称为基因集富集得分（gene set enrichment score，GSE）。然后，对于每个样本，将该样本中的所有基因的表达量按照从大到小的顺序排列，计算每个位置上所对应的基因集富集得分。最后，将这些位置上的得分进行平均或加权平均，得到该样本在该基因集上的ssGSEA得分，用于估计该样本中该免疫细胞类型的相对丰度。
函数的输出包括每个样本的ssGSEA得分。通常来说，ssGSEA得分越高，表示该样本或基因集越富集相关基因。
★ 注意：ssGSEA得分可以取任何实数值，取值范围是从负无穷到正无穷。但在实际应用中，常?；岫缘梅纸泄橐换?，使其取值范围限定在0到1之间，方便进行比较和可视化。
2.1 准备基因表达矩阵及免疫基因集列表
# 采用CIBERSORT中处理的LUAD的FPKM数据作为基因表达矩阵
FPKM[1:5, 1:5]
## TCGA-35-5375-01A-01R-1628-07 TCGA-55-A4DF-01A-11R-A24H-07 TCGA-95-8039-01A-11R-2241-07 TCGA-MP-A4T4-01A-11R-A262-07 TCGA-62-A471-01A-12R-A24H-07
## A2M 40.9806 45.5385 226.9787 125.4271 23.3876
## A4GALT 6.3743 7.8681 2.7125 9.7885 20.6624
## AAAS 6.3150 8.2626 8.2346 8.3154 12.2582
## AACS 2.6022 3.5811 1.5165 2.0158 5.4609
## AADAC 28.5316 1.8770 10.4736 0.0450 4.0560
# 从TISIDB下载免疫基因集列表（http://cis.hku.hk/TISIDB/data/download/CellReports.txt）
library(tidyverse)
cellMarker <- read.csv("CellReports.txt", header = F, sep = "\t") # 用EXCEL打开删除NA列
cellMarker <- cellMarker %>% column_to_rownames("V1") %>% t()
a <- cellMarker
a <- a[1:nrow(a), ]
set <- colnames(a)
geneSet <- list()
# i <- "Activated CD8 T cell"
for (i in set) {
x <- as.character(a[,i])
x <- x[nchar(x)!=0]
x <- as.character(x)
geneSet[[i]] <-x
}
ssGSEA-1
2.2 ssGSEA分析
library(GSVA)
ssgsea <- gsva(FPKM, geneSet, method='ssgsea', kcdf='Gaussian', abs.ranking=TRUE)
## Estimating ssGSEA scores for 28 gene sets.
## |=======================================================================================================================| 100%
a <- ssgsea %>% t() %>% as.data.frame()
identical(rownames(a), rownames(group))
a$group <- LUAD_data$group # CIBERSORT中的数据
a <- a %>% rownames_to_column("sample")
write.table(a, "ssGSEA.txt", sep = "\t", row.names = T, col.names = NA, quote = F)
# Min-Max标准化是指对原始数据进行线性变换，将值映射到[0，1]之间
# 这里是将每个样本中不同的免疫细胞比例标准化到0-1之间
ssgsea.1 <- ssgsea
for (i in colnames(ssgsea)) {
#i <- colnames(ssgsea)[1]
ssgsea.1[,i] <- (ssgsea[,i] -min(ssgsea[,i]))/(max(ssgsea[,i] )-min(ssgsea[,i] ))
}
apply(ssgsea.1[,1:6], 2, range)
2.3 可视化
LUAD_ggsea <- gather(a,key = ssgsea, value = Expression, -c(group,sample))
ggplot(LUAD_ggsea, aes(x = ssgsea, y = Expression)) +
labs(y="Expression", x = NULL) +
geom_boxplot(aes(fill = group), position = position_dodge(0.5), width = 0.5, outlier.alpha = 0) +
scale_fill_manual(values = c("#096EA9", "#B33D27")) +
theme_bw() +
theme(plot.title = element_text(size = 12,color="black",hjust = 0.5),
axis.text.x = element_text(angle = 45, hjust = 1 ),
panel.grid = element_blank(),
legend.position = "top",
legend.text = element_text(size= 12),
legend.title= element_text(size= 12)) +
stat_compare_means(aes(group = group),
label = "p.signif",
method = "wilcox.test",
hide.ns = T)
ssGSEA-2
