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