R 语言分析 2-2：免疫渗透分析（ssGSEA）

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