单细胞之富集分析-7：Pseudobulk differenti

2022-10-13 本文已影响0人 Hayley笔记

单细胞富集分析系列：

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)

这样就得到了差异分析结果，可以继续做通路富集了。

参考：