单细胞实战之pseudobulks分析，GSVA富集分析——入门到进阶(初级篇3）
原创单细胞实战之pseudobulks分析，GSVA富集分析——入门到进阶(初级篇3）
原创
该推文首发于公众号：单细胞天地
我们在上一讲内容中学习了“矫正”数据的三个工具~ 分别为细胞周期矫正，去除双胞体和RNA污染。接下来我们再来学习一下pseudobulks分析、GSVA富集分析。初级篇2最后得到的工程文件通过网盘分享：初级篇2 链接: https://pan.baidu.com/s/1ETVEkJNCJSe9NU7_aLEPVA 提取码: znmb
此外，可以向“生信技能树”公众号发送关键词‘单细胞’，直接获取Seurat V5版本的完整代码。
差异和富集分析
差异分析通过识别显著变化的基因，为我们提供了潜在的生物学信息，而通过将这些差异基因与生物学功能进行关联，富集分析进一步帮助我们理解基因在细胞过程中的作用。传统的富集分析方法通常依赖于基因集的统计学显著性，而基因集变异分析（GSVA）作为一种无监督的分析方法，可以基于样本的基因表达数据计算每个样本在某一基因集中的富集度，进而提供更加细致的功能注解。通过将差异分析与pseudobulks结合，我们能够有效地合并单细胞或小样本中的信息，获得更具代表性的群体特征，这对于高变异数据的分析尤其重要。GSVA在此过程中扮演着至关重要的角色，它可以准确地捕捉不同样本之间的功能差异，揭示潜在的功能性变化，帮助我们深入理解细胞功能的动态变化，推动生物学研究向更深层次发展。
1.pseudobulks分析
在Seurat官方流程中的Differential expression testing部分中展示了AggregateExpression函数的pseudobulks这差异分析。开发者提到经典的FindMarkers的方法是将每一个细胞看做独立重复，并忽略了源自同一个样品细胞之间的内在相关性，那么这种分析方式已经被多篇文献正式存在假阳性
那么接下来我们也尝试使用大模型来回答两个问题：
- pseudobulks分析是什么；
- 为什么需要进行 pseudobulks分析
大模型给出的回答如下：
什么是 Pseudobulk 分析？
在单细胞RNA测序（scRNA-seq）分析中，Pseudobulk分析 是一种将单个细胞的基因表达数据聚合成群体（或称为“伪样本”）的策略。这一方法将细胞分组（例如，同一细胞类型或同一实验组内的细胞）中的基因表达数据进行汇总，生成一个类似传统RNA测序（bulk RNA-seq）数据的“伪样本”。这些“伪样本”代表了群体内的基因表达模式，可以用于后续的差异表达分析和其他统计分析。
Pseudobulk分析的步骤：
- 细胞分组：将单细胞数据按照预定义的标准进行分组，如细胞类型、实验条件、样本组等。
- 表达数据汇总：对于每一组细胞，将同组内所有细胞的基因表达数据加总或平均，形成一个群体层次的表达谱（伪样本）。
- 差异表达分析：对伪样本进行标准的差异表达分析，以识别群体间的基因表达差异。
为什么需要进行 Pseudobulk 分析？
- 减少单细胞数据的噪声：单细胞RNA测序数据通常会受到技术噪声和批次效应的影响，这会使得基因表达量的测量具有较高的不确定性。通过将同一组细胞的表达数据聚合成伪样本，Pseudobulk分析可以有效降低这些噪声，提高分析的稳定性。
- 提高统计功效：单细胞数据通常由数千甚至数万个细胞组成，而每个细胞的基因表达水平可能有很大的变异，导致差异表达分析的统计功效较低。Pseudobulk分析通过聚合细胞数据，增强了信号与噪声的比率，从而增加了统计分析的功效。
- 克服单细胞RNA-seq的限制：由于单细胞RNA-seq的技术局限性，直接进行单细胞差异表达分析可能会受到较大的噪声干扰，导致结果不够可靠。Pseudobulk分析通过汇总信息，模拟了传统RNA-seq的数据特性，从而使得分析更加稳健。
- 适用于批量处理分析：在处理不同组或不同条件下的大量样本时，Pseudobulk分析可以减少计算资源的消耗，因为它不需要处理成千上万个单细胞数据，而是基于汇总后的数据进行分析。
- 便于后续分析：Pseudobulk数据可以直接与传统的bulk RNA-seq数据进行比较，或者用于与其他组学数据（如表观遗传学、蛋白质组学等）的联合分析。
总的来说，大模型确实提供了一定价值的回复。那么对于笔者来说，通常会使用这个工具查看不同生物学分组的差异分析结果。比如在今天的数据集中我们就尝试进行左右半结肠的Pseudobulk差异分析。
2.GSVA分析
至于GSVA分析，我认为不需要细讲了，简单来说就是一种用于评估基因集在不同样本中的活动或富集程度的方法。同时我们也可以多问一个问题，它跟GSEA分析方法有什么区别，让我们来问一下大模型。
大模型给出的回答如下：
主要区别：
1.分析对象不同
GSEA：评估整个基因集在不同分组的样本之间的富集情况。例如，通过比较肿瘤组与正常组中基因集的富集程度，判断该基因集是否在某一组中显著富集。
GSVA：评估每个单个样本中的基因集活动水平。计算每个基因集在每个样本中的活动评分，而不依赖于样本的分组。因此，即使没有预先定义的样本组，GSVA 也可以评估基因集在每个样本中的活跃度。
2.样本分组的需求
GSEA：需要预定义的样本分组，如疾病组与健康组、不同的临床亚型等。GSEA 计算基因集在这些分组之间的差异，基于每个基因在所有样本中的排序进行分析。
GSVA：不需要事先对样本进行分组，它为每个样本计算基因集的评分，从而评估该基因集在单一样本中的表现。适用于没有明确分组标签的数据，如单细胞 RNA-seq 数据。
3.评分方式不同
GSEA：使用排名统计方法。首先对基因进行排序，然后计算基因集在样本中的富集情况。GSEA 主要关注基因集在样本中的分布，尤其是基因表达水平的排序。
GSVA：使用加权排名分数。GSVA 通过评估基因集中的基因在样本中的表达情况，计算一个“活动评分”，反映每个样本中该基因集的活跃程度。
4.输出结果不同
GSEA：输出基因集在两组样本（或不同分组的样本）之间的显著性差异，通常包括：p 值（p-value）；NES（Normalized Enrichment Score，归一化富集分数）；FDR（False Discovery Rate，假发现率）；GSVA：输出每个基因集在每个样本中的活动评分，这些评分可用于比较不同样本中基因集的表现。
5.适用场景
GSEA：主要用于已分组的样本（如肿瘤 vs. 正常样本）之间的基因集富集分析。适用于样本组间差异分析，通常应用于疾病研究、临床分型等场景。
GSVA：适用于单细胞数据或没有明确分组的样本。GSVA 可在没有组别信息的情况下评估基因集在每个样本中的活动，适合大规模样本数据和细粒度差异分析。
感觉对于简单问题，大模型还是能够回答的有模有样。
分析步骤—pseudobulks及后续分析
1. 导入
再次提醒，我们这里用的数据只进行了双胞体+环境RNA污染过滤
rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
library(stringr)
library(patchwork)
library(ggsci)
library(ggpubr)
library(RColorBrewer)
library(qs)
sce.all <- qread("./4-Corrective-data/sce.all.qs")
colnames(sce.all@meta.data)
table(sce.all@meta.data$DF.classifications)
# Singlet
# 24080
range(sce.all@meta.data$contamination)
# [1] 6.286221e-05 1.997864e-01
dir.create("5-pseudobulks")
setwd("5-pseudobulks")1.1 pseudobulks+相关性分析
现在代码好简单，直接用AggregateExpression或者AverageExpression即可
具体内部参数调整可以直接看函数的帮助文档
# v4 是 AverageExpression
# 设置默认层为 counts（或 data，视需求而定）
# 提取想要的数据进行分析
# 相关性分析
av <-AggregateExpression(sce.all,
group.by = c("orig.ident","celltype"),
assays = "RNA",
#layer = "counts",
return.seurat = FALSE) # 返回总的计数
av=as.data.frame(av[[1]])
head(av)[1:3,1:3] # 可以看到是整数矩阵
# GSM5688706-WGC_B cells GSM5688706-WGC_endothelial cells GSM5688706-WGC_epithelial/cancer cells
# RP11-34P13.7 0 0 2
# FO538757.2 29 4 223
# AP006222.2 30 3 94
#write.csv(av,file = 'AverageExpression-0.8.csv')
# 相关性分析
# 找到sd最显著的排名前1000基因
cg=names(tail(sort(apply(log(av+1), 1, sd)),1000))
# 提前这些基因矩阵并进行log处理后计算不同样本对应的细胞的相关性值
df =cor(as.matrix(log(av[cg,]+1)))
colnames(df)
ac=as.data.frame(str_split(colnames(df),'_',simplify = T))
rownames(ac)=colnames(df)
ac$V1
ac$group = ifelse(grepl('GSM5688706-WGC|GSM5688707-JCA|GSM5688708-LS-CRC3',ac$V1),"left","right")
table(ac$group)
head(ac)
pheatmap::pheatmap(df ,
show_colnames = F,show_rownames = F,
annotation_col = ac,
filename = 'cor_celltype-vs-orig.ident.pdf' )
#dev.off()
save(av,file = 'av_for_pseudobulks.Rdata')检查一下相关性，基本上不同细胞亚群之间区分的挺好的，同时我们也可以看到免疫相关性细胞(B cells, mast cells, plasma, myeloid cells, proliferative cells, T/NK cells), 上皮/癌症细胞(epithelial/cancer cells), 其他细胞(endothelial cells, fibroblasts和VSMCs），分别存在整体的相关性。
1.2 pseudobulks+主成分分析
# 3.对于每个细胞亚群进行主成分分析##################
av <-AggregateExpression(sce.all ,
group.by = c("orig.ident","celltype"),
assays = "RNA")
av=as.data.frame(av[[1]])
df=log(av +1)
head(ac)
celltp = unique(ac$V2);celltp
library(patchwork)
library("FactoMineR")
library("factoextra")
library(ggstatsplot)
library(pheatmap)
pca_list <- lapply(celltp, function(x){
# x=celltp[1]
x
#～～～主成分分析图p2～～～
exp <- df[,rownames(ac[ac$V2==x,])]
cg=names(tail(sort(apply(exp, 1, sd)),1000))
exp=exp[cg,]
dat.pca <- PCA(as.data.frame(t(exp)) , graph = FALSE)
group_list=ac[ac$V2==x,'group']
this_title <- paste0(x,'_PCA')
p2 <- fviz_pca_ind(dat.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = group_list, # color by groups
palette = "Dark2",
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups")+
ggtitle(this_title)+
theme_ggstatsplot()+
theme(plot.title = element_text(size=12,hjust = 0.5))
p2
})
wrap_plots(pca_list, byrow = T, nrow = 2 )
ggsave('all_pca.pdf',width = 12,height = 6)
ggsave('all_pca.png',dpi = 300,width = 12,height = 6)主成分检查
1.3 pseudobulks + 差异分析
av <-AggregateExpression(sce.all,
group.by = c("orig.ident","location"),
assays = "RNA",
slot = "counts",
return.seurat = FALSE) # 返回总的计数
av=as.data.frame(av[[1]])
head(av)[1:3,1:3] # 可以看到是整数矩阵
# B cells_left B cells_right endothelial cells_left
# RP11-34P13.7 1 1 1
# FO538757.2 47 63 23
# AP006222.2 75 93 32
####### Deseq2
library(tibble)
library(DESeq2)
# 1. Get counts matrix
counts_res <- av
head(counts_res)
# GSM5688706-WGC_left GSM5688707-JCA_left GSM5688708-LS-CRC3_left GSM5688709-RS-CRC1_right GSM5688710-R-CRC3_right
# RP11-34P13.7 6 6 9 3 3
# FO538757.2 440 304 176 307 116
# AP006222.2 241 183 400 374 201
# RP4-669L17.10 5 0 4 15 0
# RP11-206L10.9 97 97 58 84 79
# LINC00115 122 64 110 84 81
# GSM5688711-R-CRC4_right
# RP11-34P13.7 2
# FO538757.2 280
# AP006222.2 131
# RP4-669L17.10 2
# RP11-206L10.9 57
# LINC00115 43
# 2. generate sample level metadata
colData <- data.frame(samples = colnames(counts_res))
colData <- colData %>%
mutate(condition = ifelse(grepl('left', samples), 'left', 'right')) %>%
column_to_rownames(var = 'samples')
# 3. perform DESeq2 --------
# Create DESeq2 object
dds <- DESeqDataSetFromMatrix(countData = counts_res,
colData = colData,
design = ~ condition) # condition 表示差异分析将基于colData的condition 变量进行
# filter
keep <- rowSums(counts(dds)) >=10
dds <- dds[keep,]
# run DESeq2
dds <- DESeq(dds)
# Check the coefficients for the comparison
resultsNames(dds)
# Generate results object
res <- results(dds, name = "condition_right_vs_left")
res
# log2 fold change (MLE): condition right vs left
# Wald test p-value: condition right vs left
# DataFrame with 19708 rows and 6 columns
# baseMean log2FoldChange lfcSE stat pvalue padj
# <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
# RP11-34P13.7 4.65497 -1.1013894 1.292952 -0.8518411 0.394302 NA
# FO538757.2 266.24931 -0.2441889 0.504740 -0.4837917 0.628534 0.999255
# AP006222.2 245.41566 -0.0125393 0.440271 -0.0284809 0.977279 0.999255
# RP4-669L17.10 3.71211 0.7664332 2.178159 0.3518720 0.724934 NA
# RP11-206L10.9 79.98878 0.0664181 0.517006 0.1284666 0.897780 0.999255
# ... ... ... ... ... ... ...
# NLGN4Y 7.54510 6.49146 3.23272 2.00805 0.0446379 0.822548
# TTTY14 6.19086 6.20600 3.36549 1.84401 0.0651814 NA
# RP11-424G14.1 6.19086 6.20600 3.36549 1.84401 0.0651814 NA
# KDM5D 82.60922 23.13708 3.90737 5.92140 NA NA
# AP001626.2 2.70850 5.01300 3.92520 1.27713 0.2015560 NA
resOrdered <- res[order(res$padj),]
head(resOrdered)
DEG =as.data.frame(resOrdered)
DEG_deseq2 = na.omit(DEG)
#添加上下调信息
DEG_deseq2 <- DEG_deseq2 %>%
mutate(Type = if_else(padj > 0.05, "stable",
if_else(abs(log2FoldChange) < 1, "stable",
if_else(log2FoldChange >= 1, "up", "down")))) %>%
arrange(desc(abs(log2FoldChange))) %>% rownames_to_column("Symbol")
# ggplot绘图
ggplot(DEG_deseq2, aes(log2FoldChange,-log10(padj))) +
geom_point(size = 3.5, alpha = 0.8,
aes(color = Type),show.legend = T) +
scale_color_manual(values = c("#00468B", "gray", "#E64B35")) +
ylim(0, 15) +
xlim(-10, 10) +
labs(x = "Log2(fold change)", y = "-log10(padj)") +
geom_hline(yintercept = -log10(0.05), linetype = 2, color = 'black',lwd=0.8) +
geom_vline(xintercept = c(-1, 1), linetype = 2, color = 'black',lwd=0.8)+theme_bw()+
theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())
ggsave("volcano.pdf",width = 9,height = 7)
分析步骤——GSVA分析
1.导入
rm(list = ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
library(stringr)
library(ggsci)
library(patchwork)
library(ggpubr)
library(RColorBrewer)
library(msigdbr)
library(qs)
library(GSVA)
sce.all <- qread("./4-Corrective-data/sce.all.qs")
colnames(sce.all@meta.data)
dir.create("6-GSVA")
setwd("6-GSVA")2.数据预处理
sce <- sce.all
genesets <- msigdbr(species = "Homo sapiens", category = "C2")
genesets <- subset(genesets, select = c("gs_name","gene_symbol")) %>% as.data.frame()
genesets <- split(genesets$gene_symbol, genesets$gs_name)
# 这里的分组会决定接下来的分析哦！
Idents(sce) <- sce$celltype
expr <- AverageExpression(sce, assays = "RNA", slot = "data")[[1]]
expr <- expr[rowSums(expr)>0,] #选取非零基因
expr <- as.matrix(expr)
head(expr)请注意这里AverageExpression用于计算每个基因在所有细胞中的平均表达值，
- assays = "RNA"：指定使用哪个assay。在Seurat对象中，一个对象可以包含多个assay（如 RNA、integrated、SCT 等）。这里指定的是 RNA assay，它通常包含基因表达数据。
- slot = "data"：指定从assay中提取哪个槽（slot）。对于 RNA assay，"data" 槽通常包含标准化后的基因表达值
3.GSVA分析
# gsva默认开启全部线程计算
gsvaPar <- gsvaParam(expr, genesets,maxDiff = TRUE)
gsvaPar
gsva.res <- gsva(gsvaPar)
dim(gsva.res)
head(gsva.res)[1:5,1:5]
write.csv(gsva.res,"gsva.res.csv")
gsva.df <- data.frame(Genesets=rownames(gsva.res), gsva.res, check.names = F)
gsva_d = gsva.res[sample(nrow(gsva.res),30),]
pheatmap::pheatmap(gsva_d, show_colnames = T,
scale = "row",angle_col = "45",
color = colorRampPalette(c("navy", "white", "firebrick3"))(50))
#数据来源是上面的
# 气泡图
library(reshape2)
gsva_long <- melt(gsva_d, id.vars = "Genesets")
# 创建气泡图
ggplot(gsva_long, aes(x = Var2, y = Var1, size = value, color = value)) +
geom_point(alpha = 0.7) + # 使用散点图层绘制气泡，alpha设置点的透明度
scale_size_continuous(range = c(1, 6)) + # 设置气泡大小的范围
theme_bw() +
scale_color_gradient(low = "#336699", high = "tomato") +
labs(x = "Gene Set", y = "Sample", size = "GSVA Score")+
ggtitle("GSVA analysis") +
theme(axis.text.x = element_text(angle = 45,vjust = 0.5,hjust = 0.5),
plot.title = element_text(hjust = 0.5))
ggsave('gsva_dot_C2.pdf',width = 10,height = 7)结果应该不需要重点去解释
参考资料：
- Seurat-AggregateExpression：https://satijalab.org/seurat/articles/de_vignette
- GSVA github：https://github.com/rcastelo/GSVA
- 生信技能树：https://mp.weixin.qq.com/s/jR2OdJQPfBAfxSLSYPyqUw
致谢：感谢曾老师以及生信技能树团队全体成员。更多精彩内容可关注公众号：生信技能树，单细胞天地，生信菜鸟团等公众号。
注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多相关内容可关注公众号：生信方舟
- END -
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
腾讯云开发者
