R语言进行单细胞转录组GSVA分析

R语言进行单细胞转录组GSVA分析

​ GSVA（基因集变异分析）是一种无监督的方法，用于评估样本中基因集的表达变化，并在单细胞转录组数据分析中特别有用。通过比较不同细胞亚群中基因集的活性差异，GSVA可以揭示细胞功能状态、识别疾病相关的细胞亚群，并预测对药物治疗的反应性。这些分析有助于理解组织中的细胞异质性和功能特异性，对疾病机理和治疗策略的研究至关重要。

library(Seurat)
library(SeuratData)
library(GSVA)
library(msigdbr)
library(clusterProfiler)
library(limma)

1. 载入数据集

# 之前经过seurat分析的单细胞数据集
pbmc <- readRDS('pbmc_Macrophages.rda')

# 查看分组信息
table(pbmc@meta.data$seurat_clusters) 
# 输出:   0   1   2   3   4 
#       819 365 154 125  47 

# 绘制UMAP降维图
DimPlot(pbmc, reduction = "umap", label = F, pt.size = 0.5) 

2. 查找每个细胞亚群的特征基因

查找每个细胞亚群的特征基因(Markers)在生物学研究和临床诊断中扮演着关键角色。这些基因的表达差异不仅助于精确定义和区分细胞的不同亚群，还可以揭示每种亚群的具体生物学功能和行为。例如，免疫细胞中的特征基因可以指示其在抗感染反应、炎症或免疫调节中的作用。

# 寻找每个数据集中的差异基因
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25, max.cells = 200)

# 筛选每种细胞首3位的Marker基因！！！
top3 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 3, wt = avg_log2FC)

# 绘制点阵图
DotPlot(pbmc, features = top3$gene[!duplicated(top3$gene)]) +  # 绘制点阵图，展示top3列表
  coord_flip() +  # 翻转坐标轴，使基因名在y轴，聚类信息在x轴
  theme_bw() +  # 使用黑白主题，简洁风格
  theme(panel.grid = element_blank(),  # 移除图面的网格线
        axis.text.x = element_text(hjust = 0.5, vjust = 0.5)) +  # 调整x轴文字对齐方式
  labs(x = NULL, y = NULL) +  # 移除x、y轴的标签
  guides(size = guide_legend(order = 3)) +  # 调整图例的显示顺序
  scale_color_gradientn(values = seq(0, 1, 0.2),  # 定义颜色渐变的值，从0到1
                        colours = c('#330066', '#336699', '#66CC66', '#FFCC33'))  

3. GSVA分析

GSVA（基因集变异分析）是一种用于生物信息学的统计方法，它可以评估一个预定义的基因集在样本间表达水平的变化。这个方法主要用于分析转录组数据（如RNA测序数据），帮助研究者了解在不同条件下，特定生物过程或信号通路的活动变化。

提取巨噬细胞平均值矩阵

# 设置巨噬细胞的Seurat聚类为标识符
Idents(pbmc) <- "seurat_clusters"

# 提取平均值
expr <- AverageExpression(pbmc, assays = "RNA",  slot = "data")[[1]] 

#选取非零基因
expr <- expr[rowSums(expr)>0,] 

# 转换成矩阵格式，以便后续分析
expr <- as.matrix(expr)

准备GSVA参考基因集

# 从Molecular Signatures Database获取小鼠生物过程(BP)基因集
genesets <- msigdbr(species = "Mus musculus",category = "C5",subcategory="BP") 

# 选择需要的列并转换为数据框
genesets <- subset(genesets, select = c("gs_name","gene_symbol"))%>% as.data.frame()

# 按基因集名称对基因符号进行分组
genesets <- split(genesets$gene_symbol,genesets$gs_name)

对表达矩阵进行GSVA分析

gsva_res <- gsva(expr,genesets, method="gsva")

提取GSVA结果

gsva.df <- data.frame(Genesets=rownames(gsva_res), gsva_res, check.names = F)

# 查看结果
print(gsva.df)

富集条目格式整理

library(stringr)

#去掉首列 
gsva.df <- gsva.df[,-1] 

# 将下划线替换为空格，并去除前缀
rownames(gsva.df) <- gsub("_"," ",substr(rownames(gsva.df), 6, nchar(rownames(gsva.df))))

# # 转换为句子格式（首字母大写）
rownames(gsva.df) <- str_to_sentence(str_to_lower(rownames(gsva.df))) 

# 设置细胞亚群名称
colnames(gsva.df) <- c("Cluster 0", "Cluster 1", "Cluster 2",
                       "Cluster 3", "Cluster 4")
# 查看结果
print(gsva.df)

保留每个单细胞亚群中 GSVA值最大的三个条目

# 使用apply函数获取每列中最大的三个值的行索引
rows_to_keep <- unique(unlist(lapply(gsva.df, function(col) {
  order(col, decreasing = TRUE)[1:3] # 保留前3位
})))

# 根据索引保留这些行
filtered_gsva.df <- gsva.df[rows_to_keep, ]

# 查看结果
print(filtered_gsva.df)

4. 使用pheatmap 包绘制热图

# 加载绘图库
library(pheatmap) 

# 使用pheatmap函数绘制热图
pheatmap(filtered_gsva.df,  # 输入的数据框，其中包含要展示的GSVA分析结果
         show_colnames = TRUE,  # 显示列名，即每个聚类的名称
         angle_col = "315",  # 列名的显示角度，设为315度使得列名从右上斜向左下显示
         cluster_rows = FALSE,  # 不进行行聚类，即不根据行数据对行进行聚类排序
         cluster_cols = FALSE,  # 不进行列聚类，即不根据列数据对列进行聚类排序
         border_color = "white",  # 设置单元格边框颜色为白色
         fontsize_col = 12,  # 设置列名的字体大小为12
         fontsize_row = 10)  # 设置行名的字体大小为10