之前绘制过FindMarkers/FindAllmarkers差异分析后的单细胞差异基因火山图，除了FindMarkers/FindAllmarkers这种方法以外，pseudobulks是另一种单细胞差异基因分析的方法，这次就来学习和整理一下。
Pseudobulk 分析
概念：
● Pseudobulk分析将单细胞RNA测序数据中的细胞按特定的条件（如样本、群体、时间点等）聚合为“伪散装”样本，然后对这些聚合样本进行差异表达分析。
● 每个“伪散装”样本的表达量通常是将属于该组的细胞的表达数据求和或取平均值得到的。这种方法可以将单细胞数据转换为类bulk RNA-seq数据进行处理。
优点：
● 通过聚合细胞数据，减少了单细胞数据的稀疏性和高变异性，使分析更接近于传统的bulk RNA-seq差异分析方法。
● 可以直接应用成熟的bulk RNA-seq差异分析工具，如DESeq2或edgeR。
缺点：
● 聚合过程会丢失细胞间的异质性信息，只能反映群体水平的差异。
● 在细胞数较少的情况下，可能会导致伪散装样本的数据不稳定。
FindMarkers/FindAllmarkers 分析
概念：
● FindMarkers/FindAllmarkers是直接在单细胞RNA测序数据的细胞级别进行差异表达分析。它使用统计检验方法（如Wilcoxon秩和检验、MAST或t检验）来比较不同细胞群体之间的基因表达差异。
优点：
● 保留了单细胞数据的分辨率，能够捕捉到细胞群体内的异质性。
● 可以分析细胞亚群之间的差异，适合于细胞类型复杂的研究。
缺点：
● 由于单细胞数据的稀疏性和高噪声，差异表达分析的结果可能不如pseudobulk分析稳定。
● 分析结果依赖于选择的细胞群体，因此需要谨慎进行群体定义和数据预处理。
异同点总结
相同点：
● 两者都用于识别在不同条件或群体之间存在差异表达的基因。
● 都需要预处理和标准化单细胞RNA测序数据。
不同点：
● 数据处理方式：pseudobulk分析将细胞聚合为群体进行分析，而FindMarkers分析直接在单细胞层面上进行。
● 解析水平：pseudobulk分析关注的是群体间的差异，FindMarkers分析则关注细胞群体内及群体间的异质性。
● 工具和方法：pseudobulk分析可以使用传统的bulk RNA-seq分析工具，而FindMarkers/FindAllmarkers通常依赖于专门为单细胞数据设计的统计检验方法。
适用场景
● Pseudobulk分析：适用于样本数较多且希望降低单细胞数据噪声的研究，或希望利用传统bulk RNA-seq分析工具进行下游分析的场景。
● FindMarkers/FindAllmarkers分析：适用于细胞异质性较高的研究，或希望深入探索特定细胞亚群差异的场景。
分析流程
1、加载R包及读取数据
代码语言：javascript复制rm(list=ls())
library(qs)
library(ggplot2)
library(DESeq2)
library(Seurat)
library(tidyverse)
library(dplyr)
# 不知道R包有没有加载全，如果没有的话建议自行安装加载一下
scRNA = qread('./sce_celltype.qs')
dim(scRNA)
# [1] 20000 394912、数据处理
代码语言：javascript复制table(scRNA$tissue.type)
# CA NL
# 23082 16409
bs = split(colnames(scRNA),scRNA$sample.id)
ct = do.call(
cbind,lapply(names(bs), function(x){
# x=names(bs)[[1]]
kp =colnames(scRNA) %in% bs[[x]]
rowSums(as.matrix(scRNA@assays$RNA@layers$counts[, kp] ))
})
)
colnames(ct) <- names(bs)
head(ct)
rownames(ct) <- rownames(scRNA)split(colnames(scRNA), scRNA$sample.id)：将所有细胞的列名按 sample.id 进行分组。split 函数返回一个列表，其中每个元素包含属于同一 sample.id 的细胞列名。
bs 是一个列表，列表的每个元素代表一个组织类型（tissue.type），并包含所有属于该类型的细胞列名。
lapply(names(bs), function(x){...}): 对于每一个ID（即 names(bs) 的每一个元素），执行函数体内部的操作。
kp =colnames(scRNA) %in% bs[[x]]: 这一步通过布尔索引筛选出当前样本ID对应的所有细胞列。
rowSums(as.matrix(scRNA@assays$RNA@layers$counts[, kp])): 对选定的细胞列(不同组)中的基因表达矩阵进行行求和，得到每个基因在该样本中的总表达量。这里需要思考一下，我们使用的kp，这里的kp其实代表的是bs中的ID，所以按照这个数据而言，分别是对CA组和NL组的数据的基因表达矩阵进行行求和。
最终通过 cbind 函数将所有样本的基因表达总和结果列绑定（即按列组合），生成矩阵 ct，其中每一列对应一个样本，每一行对应一个基因。
不过此时需要注意的是，ct表格中没有行名，也就是没有基因名，因此我们需要把scRNA的行名加上去。也可以在差异分析前加上去
3、获取分组信息
代码语言：javascript复制# 获取分组信息
phe = scRNA@meta.data[,c('sample.id','tissue.type')]
phe = unique(scRNA@meta.data[,c('sample.id','tissue.type')])
phe[1:5,1:2]
# patient.id tissue.type
# AAACGGGCATGACGGA-1 P4 CA
# AAACCTGTCATCTGCC-2 P15 CA
# AAACCTGCATGTAGTC-3 P21 CA
# AAACCTGAGTCATGCT-4 P22 CA
# AAACCTGCAGCTCCGA-5 P26 CA
group_list = phe[match(names(bs),phe$sample.id),'tissue.type']
table(group_list) 第一行代码从 scRNA 对象的 meta.data 中提取两列数据：sample.id（样本ID）和 tissue.type（组织类型）。meta.data 是存储每个细胞对应的元数据信息的表格。提取后的结果 phe 是一个数据框，其中包含每个细胞的样本ID和对应的组织类型。
第二行代码使用 unique 函数对刚才提取的数据进行去重操作。unique 函数会移除数据框中重复的行，因此生成的 phe 数据框会包含每个样本ID唯一对应的一行记录，即每个样本ID对应的组织类型。这样处理后，phe 数据框的每一行代表一个样本，而不是一个细胞。
接下来的group_list代码是匹配样本ID并提取对应的组织类型：
names(bs): 这个部分提取的是之前创建的列表 bs 中的样本ID（样本的列名）。
match (names(bs), phe $ sample.id) : match函数用于在 phe$sample.id 中查找与 names(bs) 相匹配的样本ID，返回匹配到的位置索引。简单来说，它会告诉你每个 bs 列表中的样本ID在 phe 数据框中的位置。
phe[...]: 这里使用这些位置索引来从 phe 数据框中提取相应行的 tissue.type 列，最终得到的 group_list 是一个向量，包含了 bs 中样本ID对应的组织类型。
4、过滤数据
代码语言：javascript复制# 赋值并对每一行的
exprSet = ct
dim(exprSet)
exprSet=exprSet[apply(exprSet,1, function(x) sum(x>1) > 1),]
dim(exprSet)
table(group_list)
group_list <- factor(group_list,levels = c("NL","CA"))apply(exprSet, 1, function(x) sum(x > 1) > 1) 这段代码对 exprSet 矩阵的每一行（即每个基因）进行操作：
apply(exprSet, 1, function(x) ...)：apply 函数在矩阵的每一行（1 表示行操作）上应用给定的函数。sum(x > 1) > 1：对于每个基因（每行），计算在多少个样本（列）中该基因的表达量大于1，如果该数量大于1（即至少在两个样本中有表达量大于1），则保留该基因。
建议group_list设置成因子模式
5、接下来就是差异分析和绘图，除了DESeq2分析，还可以用limma，edgeR分析。
转录组差异分析方法可参考推文：https://mp.weixin.qq.com/s/dCgLjzOGKeSpZlpFMzIR9g
代码语言：javascript复制# DESeq2分析
colData <- data.frame(row.names=colnames(exprSet),group_list=group_list)
# exprSet <- apply(exprSet, 2, as.integer) 这一步不需要哦，可以想一想
# rownames(exprSet) <- rownames(scRNA) 如果这步上面没有做，可以再加上去
dds <- DESeqDataSetFromMatrix(countData = exprSet,
colData = colData,
design = ~ group_list)
dds <- DESeq(dds)
table(group_list)
res <- results(dds,
contrast=c("group_list",
levels(group_list)[2],
levels(group_list)[1]))
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)
参考资料：
1、生信技能树：https://mp.weixin.qq.com/s/jR2OdJQPfBAfxSLSYPyqUw
2、爱喝可乐的天天：https://mp.weixin.qq.com/s/lgQsfQ60uyQrQRk0gb2fLg
致谢：感谢曾老师以及生信技能树团队全体成员。
注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟
- END -
