单细胞富集分析系列：
- 单细胞之富集分析-1：单细胞GSEA分析流程
- 单细胞之富集分析-2：批量GSEA和GSVA分析
- 单细胞之富集分析-3：GO和KEGG富集分析及绘图
- 单细胞之富集分析-4：分组水平GSVA
- 单细胞之富集分析-5：一张图展示多组富集分析结果
- 单细胞之富集分析-6：PROGENy
1. 数据准备
使用我们在hdWGCNA：单细胞WGCNA分析方法中使用的数据集。这个数据集包含了11个小鼠的样本。
seurat_obj <- readRDS('Zhou_2020.rds')
table(seurat_obj$Sample)
# C1 C11 C12 C2 C3 C4 C5 C6 C7 C8 C9
# 1788 3345 2333 1899 3002 711 6841 3320 2934 7618 2880
##但是该数据集的11个样本都是同一个组，因此需要重新设置一下分组（仅做演示）
seurat_obj$Group <- seurat_obj$Sample
seurat_obj$Group <- recode(seurat_obj$Group,
C1='HC',
C2='HC',
C3='HC',
C4='HC',
C5='HC',
C6='HC',
C7='Treat',
C8='Treat',
C9='Treat',
C11='Treat',
C12='Treat')
table(seurat_obj$Group)
# HC Treat
# 17561 19110
2. Pseudobulk数据整合
2.1 数据整合
这一步的目的主要是为了得到Pseudobulk的矩阵
- 方法1：提取矩阵出来做加和
bs = split(colnames(seurat_obj),seurat_obj$Sample)
ct = do.call(
cbind,lapply(names(bs), function(x){
# x=names(bs)[[1]]
kp =colnames(seurat_obj) %in% bs[[x]]
rowSums( as.matrix( seurat_obj@assays$RNA@counts[, kp] ))
})
)
colnames(ct) <- names(bs)
head(ct)
# C1 C11 C12 C2 C3 C4 C5 C6 C7 C8 C9
# MIR1302-2HG 3 13 15 8 12 4 2 6 3 42 6
# FAM138A 0 0 0 0 0 0 0 0 0 0 0
# OR4F5 0 0 1 0 1 0 1 13 0 1 7
# AL627309.1 262 243 387 120 322 46 515 859 296 616 648
# AL627309.3 0 0 0 0 0 0 0 0 0 0 0
# AL627309.2 0 0 0 0 0 0 0 0 0 0 0
- 方法2：使用
AggregateExpression函数
Idents(seurat_obj) <- 'Sample'
ct <- round(AggregateExpression(seurat_obj)[[1]])
head(ct)
# C1 C11 C12 C2 C3 C4 C5 C6 C7 C8 C9
# MIR1302-2HG 1 5 6 6 5 2 1 2 10 32 3
# FAM138A 0 0 0 0 0 0 0 0 0 0 0
# OR4F5 0 0 0 0 0 0 1 6 0 0 3
# AL627309.1 159 268 233 166 213 45 489 459 399 663 458
# AL627309.3 0 0 0 0 0 0 0 0 0 0 0
# AL627309.2 0 0 0 0 0 0 0 0 0 0 0
⚠️
呃...两个方法好像还是有差别的。
两种方法的差异可以参考：seurat-AverageExpression()源码解析
2.2 查看分组
phe = unique(seurat_obj@meta.data[,c('Sample','Group')])
phe
# Sample Group
# AGTCTTTGTTGATTCG-11 C1 HC
# CAGCTGGAGCCAGGAT-12 C11 Treat
# CACAAACAGTCATCCA-13 C12 Treat
# ACTGCTCGTACATGTC-14 C2 HC
# CAACTAGGTTGTGGAG-15 C3 HC
# CATGGCGGTTCTGTTT-16 C4 HC
# CCACCTAGTCTGCGGT-17 C5 HC
# AAACCTGGTGTGACGA-18 C6 HC
# GAGGTGAAGAAACGAG-19 C7 Treat
# TCTTTCCTCTCAACTT-20 C8 Treat
# CTTCTCTAGAGCTTCT-21 C9 Treat
2.3 测查数据并质控
group_list = phe[match(names(bs),phe$Sample),'Group']
table(group_list)
# group_list
# HC Treat
# 6 5
exprSet = ct
dim(exprSet)
# [1] 36601 11
exprSet=exprSet[apply(exprSet,1, function(x) sum(x>1) > 1),]
dim(exprSet)
# [1] 25684 11
table(group_list)
save(exprSet,group_list,file = 'input_for_deg.Rdata')
3. 差异分析
后续的分析可以参考3大差异分析r包：DESeq2、edgeR和limma。
load(file = 'input_for_deg.Rdata')
library(edgeR) #以edgeR为例
dge <- DGEList(counts=exprSet,group=group_list) #输入表达矩阵和分组信息数据
dge$samples$lib.size <- colSums(dge$counts)
dge <- calcNormFactors(dge)
design <- model.matrix(~0+group_list) #写不写0+是一样的
rownames(design)<-colnames(dge)
colnames(design)<-levels(group_list)
dge <- estimateGLMCommonDisp(dge, design)
dge <- estimateGLMTrendedDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)
fit <- glmFit(dge, design)
fit2 <- glmLRT(fit, contrast=c(-1,1)) #这里的contrast和DESeq2有差别，这里只需要输入c(-1,1)就好，-1对应着normal，是对照组，1对应着tumor，是实验组。
DEG=topTags(fit2, n=nrow(exprSet))
DEG=as.data.frame(DEG) #转化为数据框
View(DEG)
这样就得到了差异分析结果，可以继续做通路富集了。
参考：
